The effect of hydrodynamics alone on the subhalo population in a ΛCDM rich cluster sized dark matter halo

In the standard ΛCDM Cosmology, structure formation is hierarchical. Small dark matter haloes formed firstly, and then merge to form larger and larger systems. During this hierarchical clustering process, earlier accreted halo often survive as a subhalo to orbit its host. As the most massive subhaloes are expected to be the hosts of luminous satellite galaxies, properties of these subhaloes have been extensively investigated in the past two decades (e.g. Moore et al. 1998; Ghigna et al. 1998; Moore et al. 1999; Klypin et al. 1999; Ghigna et al. 2000; Springel et al. 2001a; Stoehr et al. 2002; Stoehr et al. 2003; De Lucia et al. 2004; Diemand et al. 2004; Gao et al. 2004; Kravtsov et al. 2004; Gao et al. 2012; Han et al. 2016; Han et al. 2018; van den Bosch & Ogiya 2018).

Limited by computational power, earlier works on subhaloes have been confined to use high resolution DMO simulations, whilst baryonic physics may play a sizeable impact on the subhalo population in dark matter haloes. Thanks to great advance in super computer power in recent years, properties of subhaloes have been widely studied with modern hydro-dynamical simulations with galaxy formation models (e.g. Duffy et al. 2010; Libeskind et al. 2010; Romano-Díaz et al. 2010; Garrison-Kimmel et al. 2017; Graus et al. 2018; Richings et al. 2020). Compared to DMO simulations, in general, the subhalo population is found to be less abundant in hydrodynamical simulations. However, the degree of the difference in the subhalo population between the two varies with the adopted different galaxy formation models, as structure of dark matter halo/subhalo is sensitive to the galaxy formation models. For example, Nagai & Kravtsov (2005), Macciò et al. (2006) and Weinberg et al. (2008) find that baryon has only a small impact on the subhalo population; Sawala et al. (2013), Schaller et al. (2015) and Sawala et al. (2017) show that the impact of baryonic physics depends on the subhalo mass, above 10¹² – 10¹³ M_⊙ the subhalo abundance ratio between the hydrodynamical and DMO runs is close to 1, but it drops below 1 for low mass subhaloes. Zhu et al. (2016) report a similar result in a set of simulations of Milky Way-size dark matter haloes. Chua et al. (2017) claim that baryonic physics changes the shape of the subhalo mass function: relative to the DMO simulations, subhalo is less abundant in the low-mass end and more abundant for relatively larger subhaloes in their hydrodynamical simulations.

In this paper, we compensate these studies by investigating the impact of hydrodynamics alone on the subhalo population in a ΛCDM rich cluster sized halo. Compared to earlier works on this topic (e.g., Lin et al. 2006), our NHD simulations have much higher mass and force resolutions, and have three different numerical resolutions. The latter allows us to carry out numerical convergence study. Our paper is organized as follows. In Section 2, we introduce the simulation sets used in this study. Section 3 compares the subhalo abundance, internal structure and the evolution of subhaloes in different sets of simulations.We draw our conclusion in Section 4.

The numerical simulations used in this work comprise two sets of ultra-high resolution re-simulation of a clustersized dark matter halo and its surrounding. This halo is selected from the Phoenix Project and is termed as the Ph-A halo. The Phoenix Project resimulated the formation and evolution of nine different cluster-sized dark matter haloes selected from the Millennium Simulation (Springel et al. 2005) at ultra-high resolution. We refer readers to Gao et al. (2012) for details of the Phoenix Project. The Ph-A halo has a virial mass M₂₀₀ ∼ 6.6 × 10¹⁴ h⁻¹ M_⊙ and a virial radius R₂₀₀ ∼ 1.4 × h⁻¹Mpc which is defined as the radius within which the enclosed mean density is 200 times the critical density of the Universe. In the Phoenix project the halo has been simulated at four different resolution levels, referred to as Ph-A-1, -2, -3 and -4 where 1 is the highest and 4 is the lowest level. In this study we re-run the Ph-A-2, -3 and -4 with the hydrodynamic code GADGET-3 (Springel et al. 2005) and denote them with a suffix "-g" to indicate the inclusion of hydrodynamics. The highest resolution simulation, Ph-A-2-g, has a mass resolution of 4.3 × 10⁶ h⁻¹ M_⊙ for dark matter and 8.1 × 10⁵ h⁻¹ M_⊙ for gas particles. All the simulations in this work adopt the cosmological parameters, Ω_m = 0.25, Ω_Λ = 0.75, h = 0.73, σ₈ = 0.9, n_s = 1, and Ω_b = 0.045. The detailed information of these simulations are summarized in Table 1.

Dark matter haloes of the simulations are identified by standard friends-of-friends (FOF) algorithm (Davis et al. 1985) with a linking length of 0.2 times mean interparticle separation. Based on dark matter halo catalogue, we further identify locally overdense and selfbound subhaloes with the SUBFIND algorithm (Springel et al. 2001b). The merger trees of the dark matter haloes and subhaloes are constructed using the D-Trees algorithm (Jiang et al. 2014), which identifies the descendant of an object at next output time by tracing the most bound particles of the object.

In Figure 1, we present a visual impression of the Ph-A-2 and Ph-A-2-g FOF haloes at z = 0. The upper panel compares the density distribution of total matter in the Ph-A-2 and Ph-A-2-g simulations. Note that only particles belonging to the FOF groups are used in making the projections. While the total matter distribution of the group is overall similar, some notable differences can still be found, massive subhaloes are more pronounced in the DMO simulation, and relatively central region is round in the Ph-A-2-g run. The lower panel shows separately the dark matter and the gas distribution in the Ph-A-2-g. Compared to the dark matter, the gas distribution is less clumpy, in agreement with many existing results (e.g., Fang et al. 2009).

Zoom In Zoom Out Reset image size

3.1. The Subhalo Mass and V_max Function

The left panel of Figure 2 shows the cumulative subhalo mass function for both DMO (solid lines) and NHD runs (dotted lines). Different colours are used to distinguish different resolutions as shown in the label. The cumulative subhalo mass functions are multiplied by subhalo mass in order to remove the dominant mass dependence. Apparently, massive subhaloes are more abundant in the DMO than in the NHD run, in agreement with many existing studies (e.g., Dolag et al. 2009, Schaller et al. 2015, Zhu et al. 2016, Sawala et al. 2017), while the difference becomes smaller with decreasing subhalo mass. Below 2 × 10⁹ h⁻¹ M_⊙, interestingly, the NHD runs contain more abundant subhaloes than their DMO counterparts. The other noticeable fact is that the slope of the cumulative subhalo mass function is −1.1, in the NHD runs, steeper than that of the DMO simulations (−1.0). Note that the cumulative subhalo mass functions between different resolution simulations converge down to subhaloes containing about 150 dark matter particles for both DMO and NHD simulations (vertical dotted line in the left panel of Fig. 2).

Zoom In Zoom Out Reset image size

The suppression in the subhalo mass function in hydrodynamic simulations of rich cluster sized dark matter haloes in the subhalo mass range 10¹⁰ h⁻¹ M_⊙ < M_sub < 10¹² h⁻¹ M_⊙ is also found by Dolag et al. (2009). The difference may be understood as the subhaloes in the NHD simulations suffer from additional ram pressure stripping in relative to the DMO runs (Tormen et al. 2004; Puchwein et al. 2005).

The cumulative V_max functions weighted by subhalo ${V}_{\max }^{4}$ are presented in the right panel of Figure 2. In agreement with many existing results (e.g., Zhu et al. 2016), V_max functions in the DMO runs have higher amplitude than their NHD counterparts for relatively large V_max subhaloes, while below V_max ∼ 70 km s⁻¹, the relation is reversed that the abundance of subhaloes in our NHD simulation exceeds the DMO runs. Also the slope of the V_max function is much steeper in the NHD than DMO simulations, which is −4.4 for the NHD runs and −3.4 for the DMO runs, respectively. The scale dependence of differences in the subhalo mass/V_max functions may be caused by different physics playing at different scales, e.g. subhalo mass/V_max dependence of interplay between tidal force and ram pressure as will be discussed later.

3.2. The Structure of Subhaloes

It is interesting to see how the hydrodynamics affect the internal structure of subhaloes. The left panel of Figure 3 shows V_max − M_sub relation of subhaloes in both DMO and NHD simulations with the level-2 and level-3 resolutions at z = 0, median values of V_max in each mass bin are shown. The results for the level-2 and level-3 runs are shown with solid and dashed lines, respectively; and the results for the DMO and NHD runs are distinguished with different colors as shown in the label. Interestingly, for subhaloes more massive than about 2 × 10¹⁰ h⁻¹ M_⊙, V_max − M_sub relation is nearly identical between two sets of simulations, while below this mass, subhaloes tend to have higher V_max in the NHD simulations. The difference becomes larger with decreasing subhalo mass. At M_sub = 10⁹ h⁻¹ M_⊙, the median value of V_max of subhaloes in the NHD simulations is about 30 percent larger than that of the DMO runs. These results suggest that low-mass subhaloes in the NHD runs are more concentrated than their DMO counterparts. The good agreements between different resolutions suggest that our results are robust to changes in numerical resolution.

Zoom In Zoom Out Reset image size

To more clearly see the difference in the internal structure of subhaloes between the DMO and NHD simulation, we select two 20 random low- and highmass subhalo samples from both Ph-A-2 and Ph-A-2-g simulations. The mass range of the low-mass sample is [5 × 10⁹,10¹⁰]h⁻¹ M_⊙, and it is [5 × 10¹⁰,10¹¹]h⁻¹ M_⊙ for the high-mass sample. Then we stack rotational curves of all subhaloes in each sample and display them in Figure 4, the vertical dotted lines indicate 2.8 times softening length at which the force calculation in the Gadget code becomes Newtonian. While the outer parts of rotation curves are quite similar, differences in the inner part are noticeable between the Ph-A-2 and Ph-A-2-g. In particular for the low mass subhalo samples, inner density is substantially higher for the subhaloes in the Ph-A-2-g than Ph-A-2 simulation. For the high mass subhaloes sample, even its V_max is slightly lower in the Ph-A-2-g than Ph-A-2, the inner density is still slightly higher in the Ph-A-2-g. In the Appendix, we provide one-to-one comparisons of rotational curves of four randomly selected matched individual subhaloes in both DMO and NHD runs, results are consistent with the statistical results presented here.

Zoom In Zoom Out Reset image size

Is the difference in the V_max–M_sub relation of subhalo between the DMO and NHD simulations due to interactions of subhaloes with their host halo? We trace the progenitors of the present-day subhaloes just before accretion in both DMO and NHR simulations, and plot their V_max as a function of the progenitor mass in the right panel of Figure 3. Apparently, the relation for the progenitors is quite similar to that of subhaloes, namely low-mass progenitors are statistically more concentrated in the NHD than in DMO simulations. As the accretion time distribution of the progenitors is quite broad, the low-mass progenitors are more concentrated in our NHD than in DMO simulations over a quite broad redshift range, and thus the different subhalo V_max–M_sub relation between the DMO and NHD simulations simply reflects the different structural properties between two kinds of simulations.

Indeed, previous studies suggest that, compared to the DMO simulation, dark matter haloes are slightly more concentrated in hydrodynamic simulation (Rasia et al. 2004 and Lin et al. 2006). For instance, Lin et al. (2006) found that about 3 – 8 percent difference in halo concentration parameter c between the NHD and DMO simulations for haloes with virial mass from 10¹³ to 7 × 10¹⁴ h⁻¹ M_⊙. In Figure 5, we examine the concentrationmass relation of isolated dark matter haloes with our own simulations and extend to much lower halo mass. Here the isolated haloes are selected within a cluster-centric radius [1.5, 3]R₂₀₀, and concentrated parameters of dark matter haloes are represented with V_max/V₂₀₀. Results for two different resolutions and two redshifts are shown. Black and red dots show the results for Ph-A-2 and Ph-A-2-g respectively, the different lines show median values of V_max/V₂₀₀ in each halo mass bin. Clearly, concentration parameter is larger for low-mass haloes with masses less massive than 2 × 10¹⁰ h⁻¹ M_⊙ in the NHD than DMO simulations. Above this mass, the concentration parameter of dark matter haloes are similar between two sets of simulations, slightly inconsistent with previous studies such as Lin et al. (2006). This inconsistency is mainly due to the different derivation of the concentration parameter. Lin et al. (2006) derived the concentration parameter of a dark matter halo by fitting its density profile. When we fit c rather than using the proxy V_max/V₂₀₀, our result is consistent with Lin et al. (2006).

Zoom In Zoom Out Reset image size

3.3. The Evolution of Subhaloes

In this subsection, we compare the evolution of subhaloes once they were accreted into the Ph-A halo with our matched highest resolution simulations–the Ph-A-2 and Ph-A-2-g. Following Xie & Gao (2015), in Figure 6 we present the survival number and retained mass fraction we present the survival number and retained mass fraction of the progenitors accreted at a fixed epoch z = 2 as a function of redshift. Here we consider two progenitor samples according to their halo mass, M_infall > 10¹⁰ h⁻¹ M_⊙ and 10⁹ h⁻¹ M_⊙ < M_infall < 10¹⁰ h⁻¹ M_⊙. Apparently, the survival number fraction is larger in the Ph-A-2-g in relative to the Ph-A-2 simulation. For relatively more massive progenitors, more than 90 percent of them survive to the present day in the Ph-A-2-g, while only 75 percent of them survive in the Ph-A-2. For low-mass progenitors, the survival fraction is quite similar to their high-mass counterparts in the Ph-A-2-g run, while the fraction decreases to 60 percent in the Ph-A-2 simulation. This is mainly a consequence of higher concentration of low-mass dark matter haloes in the NHD than DMO simulations. Note that the survival number fraction of subhaloes in our Ph-A-2-g simulation is consistent with that of hydrodynamic simulation with cooling and star formation of Bahé et al. (2019).

Zoom In Zoom Out Reset image size

The right panel of Figure 6 shows the evolution of the retained mass fraction of the same progenitor samples as a function of redshift. Here the retained mass fraction is defined as the ratio of the total retained mass to the total progenitor mass at each redshift. For relatively massive progenitors, their mass-loss is more rapid in the Ph-A-2-g than Ph-A-2 simulation, at intermediate redshift, say z = 1, the retained mass fraction is about three times higher in the Ph-A-2 than Ph-A-2-g simulation. This is expected because subhaloes suffer from additional rampressure stripping in the hydrodynamic simulation. For the low-mass progenitors, the relation is reverse, the mass-loss of subhaloes is less efficient in the Ph-A-2-g than Ph-A- 2. As low-mass haloes have higher concentration in NHD simulation, which counteracts the effect of ram-pressure and results in less efficient mass stripping. These subhalo mass dependence survival/retained fractions should largely account for the different shapes in the subhalo mass/V_max functions between the DMO and NHD simulations.

We perform a set of matched hydrodynamic simulations of a rich cluster-sized dark matter halo from the Phoenix Project, and compare in detail the subhalo population in the halo simulated with and without non-radiative hydrodynamics, and with three different numerical resolutions. Our results can be summarized as follows.

(1) The subhalo mass function of the cluster-sized dark matter halo has different shapes between the DMO and NHD simulations. The cumulative subhalo mass function can be approximated with N(M) ∝ M⁻¹ for the DMO simulation, while it is N(M) ∝ M^−1.1 for the NHD run. Subhaloes are more abundant in the DMO runs for subhaloes more massive than 2 × 10⁹ h⁻¹ M_⊙, below this mass, there are more subhaloes in the NHD. Similarly, the subhalo V_max function is also steeper (n ∼ −3.4 versus −4.4) in the NHD hydrodynamic than in the DMO simulation. Above V_max ∼ 70 km s⁻¹, subhaloes are more abundant in the DMO simulation, while below this value there are more subhaloes in the NHD simulation.

(2) For subhaloes less massive than 10¹⁰ h⁻¹ M_⊙, subhaloes have larger concentration in the NHD than DMO simulation, above this mass subhaloes have similar concentration in two sets of simulations. We show that this is mainly because that progenitors of the present-day subhaloes are more concentrated before accretion in the NHD than in the DMO simulation. At z = 0 and z = 1, the low-mass field dark matter haloes have larger concentration parameters in the NHD than DMO simulation.

(3) The relatively larger concentration parameters of the low mass field haloes also result in larger survival number and larger retained mass fraction of the low mass progenitors in the NHD than DMO simulation.

Based on these results, a simple picture to interpret the differences in the subhalo population between the DMO and NHD simulations can be summarized as follows. Difference in structural properties of dark matter haloes varies with halo mass between the DMO and NHD simulations, with the difference being large with the decreasing halo mass. As a result, when these dark matter halo merge into their primary halo, their survival ability is larger in the NHD than DMO simulation, leading to the more abundant lowmass subhaloes in the NHD simulation.

Note that our main results are based on a set of simulations of a single cluster-sized dark matter halo, statistical results are necessary and will be presented in a future work.

We are grateful to useful discussions with Dr. Shihong Liao. This work is supported by the National Key R&D Program of China (No. 2017YFB0203300), and the Key Program of the National Natural Science Foundation of China (NFSC, Grant No. 11733010).

As we will make one-to-one comparison of subhaloes in the DMO and NHD simulations here, we follow Schaller et al. (2015) to match subhaloes according to their particle Ids. In practise, we select subhaloeswithmore than 200 dark matter particles in each level NHD simulation. For each of them, we identify 50 most bound dark matter particles. If a subhalo in its DMO counterpart contains 25 of them, we consider the subhalo a potential candidate for a matched one in the DMO run. Now we further require the subhalo in the NHD simulation also containing half of 50 most bound darkmatter of the candidate. If both conditions are satisfied, we consider them a matched subhalo pair. For our level 2 simulations, more than 53.2% subhaloes with M_sub > 10⁹ h⁻¹ M_⊙ can be matched at z = 0, and 81.87% for subhaloes above 10¹⁰ h⁻¹ M_⊙, in consistent with Schaller et al. (2015).

In Figure A.1, we select four random pairs of matched subhaloes in [5 × 10⁹,10¹⁰]h⁻¹ M_⊙ and [5 × 10¹⁰,10¹¹]h⁻¹ M_⊙ subhalo samples, and plot their rotation curves. The overall results are consistent with the statistical ones presented in Figure 4, namely subhaloes are more centrally concentrated in the NHD than DMO runs. It is also quite interesting that, even for subhaloes with even larger V_max in the DMO run, central density is still higher in the NHD run. This should account for the fact that the survival number fraction is higher in the NHD than DMO runs even they have similar concentration.

Zoom In Zoom Out Reset image size

In Figure A.2, we compare V_max of all our matched subhalo samples. Clearly, subhaloes are more concentrated in the NHD simulations, in particularly the low mass subhaloes, consistentwith the results presented in Figure 3.

Zoom In Zoom Out Reset image size