The Impact of Baryonic Physics on the Abundance, Clustering, and Concentration of Halos

Studying the connection between galaxies and their dark matter halos is one of the keys to understanding galaxy formation and evolution, as well as constraining cosmological models. In recent years, using hydrodynamic simulations has become a popular method for investigating this connection (e.g., Vogelsberger et al. 2014a). However, these simulations are not only computationally expensive but also inconsistent. We currently lack a consensus on the correct baryonic physics prescriptions to use; thus, different hydrodynamic simulations produce widely varying results.

By contrast, dark matter–only (DMO) simulations are much less computationally expensive, and although the only physics involved is gravity, they still allow us to predict the large-scale distribution of matter in the universe. However, a DMO simulation produces a halo mass function (HMF) that is, on average, shifted to higher masses than the HMF produced by a hydrodynamic simulation with the same cosmology and initial conditions (e.g., Cui et al. 2014). This is because in a hydrodynamic simulation, the presence of baryonic physics like stellar feedback, star formation, and feedback from active galactic nuclei (AGNs) has an impact on the masses of dark matter halos. This effect is nontrivial; the impact of baryonic physics varies with halo mass, as well as halo definition, environment, and redshift. And of course, the effect of baryonic physics depends on the details of the feedback prescriptions implemented, which varies from one hydrodynamic simulation to the next.

Several recent works have focused on comparing hydrodynamic simulations to models of the galaxy–halo connection (e.g., Hadzhiyska et al. 2020, 2021a, 2021b). Beltz-Mohrmann et al. (2020) examined the accuracy of halo occupation distribution (HOD) modeling compared to hydrodynamic simulations. They extracted HOD models from the galaxies in the Illustris and EAGLE simulations and subsequently applied these HOD models to the corresponding DMO simulations (i.e., Illustris-Dark and EAGLE-Dark). They found that these HOD models, when applied to the DMO halos, were unable to reproduce the galaxy clustering seen in the hydrodynamic simulations. This was in part because the HMFs in the DMO simulations were shifted to higher masses, leading to an overestimate in the overall number density of galaxies, as well as discrepancies in other clustering statistics (e.g., correlation function and void probability function). The authors found that applying a correction to the DMO halo masses led to a significant improvement in the galaxy clustering.

This issue with DMO HMFs could present a challenge for any work attempting to use halo modeling to constrain the galaxy–halo connection (e.g., Sinha et al. 2018), as well as any work attempting to constrain cosmology (e.g., Lange et al. 2019). If one's goal is merely to produce mock galaxy catalogs with the same clustering as seen in survey data, then one does not need to worry about this mass discrepancy; in this case, the parameterization of the galaxy–halo connection is unimportant, so long as the clustering mimics the real universe. However, if one's goal is to accurately constrain the galaxy–halo connection, then this mass discrepancy presents a problem. Moreover, if one is attempting to use DMO simulations plus halo modeling to constrain cosmology, then obtaining the correct HMF is essential, unless the HOD is able to fully absorb the effect without perturbing the resulting cosmological parameters.

Several previous works have investigated the effects of baryonic physics on the HMF. Cui et al. (2012) compared three simulations: a DMO simulation, a hydrodynamic simulation with nonradiative physics, and a hydrodynamic simulation with radiative processes. They found that the fractional difference between halo masses in the hydrodynamic and DMO simulations is almost constant for halos above log $({M}_{{{\rm{\Delta }}}_{c}}/{h}^{-1}\,{M}_{\odot })\gt 13.5$, but that for higher-overdensity halos, as well as smaller-mass halos, differences in halo mass appear that depend on halo mass, as well as baryonic physics.

Later, Cui et al. (2014) further examined the effects of baryonic physics on the HMF, focusing on the role of AGN feedback. They found that for both friends-of-friends (FoF) and spherical overdensity (SO) halos, AGN feedback suppresses the HMF to a level below that of the DMO simulations. They found that the ratio between the HMFs in the hydrodynamic and DMO simulations increases with overdensity but does not have any redshift or mass dependence. They found that halos in hydrodynamic simulations have shallower inner density profiles, which lends them to halo mass loss caused by "the sudden displacement of gas induced by thermal AGN feedback." The authors provided fitting functions to describe the halo mass variations between the full-physics and DMO simulations at different overdensities, which can recover the HMFs from hydrodynamic simulations for halo masses larger than 10¹³ h⁻¹ M_⊙.

Sawala et al. (2013) examined the effect of baryons on the abundance of structures and substructures and found that halo masses are reduced for halos smaller than 10¹² M_⊙, and the effect grows with decreasing mass. Later, van Daalen et al. (2014) looked at the impact of baryonic processes on the two-point correlation functions of galaxies, subhalos, and matter in large hydrodynamic simulations and found that the changes due to the inclusion of baryons are not limited to small scales. The authors found that the large-scale effects are due to the change in subhalo mass caused by feedback associated with galaxy formation. They concluded that predictions of galaxy–galaxy and galaxy–mass clustering from models based on collisionless simulations will have errors greater than 10% on scales <1 Mpc, unless the simulation results are modified to account for the effects of baryons on the distributions of mass and satellites.

Velliscig et al. (2014) used hydrodynamic simulations from the OverWhelmingly Large Simulations project to study how the physical processes related to galaxy formation (e.g., star formation, supernova and AGN feedback, etc.) impact the properties of halos. They found that the "gas expulsion and associated dark matter expansion induced by supernova-driven winds are important for halos with masses M₂₀₀ < 10¹³ M_⊙, lowering their masses by up to 20% relative a DM-only model." They also found that AGN feedback impacts halo masses up to cluster scales (M₂₀₀ ∼ 10¹⁵ M_⊙). Moreover, they found that baryonic physics alters the total mass profiles of halos out to several times the virial radius, which cannot be explained by a change in halo concentration. They concluded that the decrease in total halo mass leads to a decrease in the HMF of 20%. The authors provided analytic fitting formulae to correct halo masses and mass functions from DMO simulations.

Subsequently, Khandai et al. (2015) investigated the properties of halos in the MassiveBlack-II hydrodynamical simulation and found that baryons strongly affect the HMF when compared to DMO simulations, while Schaller et al. (2015) examined the effects of baryons on halos in the EAGLE simulation at low masses, and Chaves-Montero et al. (2016) used subhalo abundance matching to investigate the effect of baryons on the HOD and assembly bias in the EAGLE simulation.

Bocquet et al. (2016) used the Magneticum simulations to investigate the impact of baryons on the HMF and found that baryonic effects globally decrease the masses of galaxy clusters, which, at a given mass, results in a decrease of their number density. They found that this effect disappears at high redshift (z = 2) and for high halo masses (> 10¹⁴ h⁻¹ M_⊙). They concluded that when using a survey like eROSITA, ignoring the impact of baryonic physics on the HMF leads to an underestimate in Ω_m by about 0.01. The authors also provided HMF fitting formulae.

Despali & Vegetti (2017) examined the impact of baryonic physics on the subhalo populations in EAGLE and Illustris and found that the presence of baryons reduces the number of subhalos, especially at the low-mass end. They found that the variations in the subhalo mass function depend on those in the HMF, which is shifted by the effect of stellar and AGN feedback.

Finally, Castro et al. (2021; and, earlier, Balaguera-Antolínez & Porciani 2013) addressed the impact of baryonic physics on clusters. They found that ignoring baryonic effects on the halo mass function and halo bias could significantly alter the cosmological parameter constraints, particularly in the upcoming generations of galaxy cluster surveys.

Most of these works concur that baryonic physics lead to a net reduction in the masses of halos and, consequently, the HMF. However, the magnitude and mass dependence of the effect differ between different works, likely because of the different hydrodynamic simulations used. This calls for a systematic study that compares the impact of baryons on the HMF across a variety of recent hydrodynamic simulations with different baryonic physics and feedback prescriptions and over a wide range of redshifts and halo definitions. Furthermore, no previous work has examined the environmental dependence of baryonic effects on the HMF or simultaneously examined the impact of baryonic physics on halo clustering and concentrations as a function of mass, all of which are important ingredients for halo models of galaxy clustering.

In this work, we investigate the impact of baryonic physics on the HMF in three different simulations: Illustris, IllustrisTNG, and EAGLE. We consider three different redshifts and five different halo definitions. Moreover, we investigate the environmental dependence of the effect of baryons on the HMF. For all of these cases, we provide formulae that can be used to correct the halo masses from DMO simulations to match the HMFs from their hydrodynamic counterparts. Finally, we investigate the effect of baryonic physics on halo clustering and the halo mass–concentration relation.

We discuss the simulation data in Section 2. In Section 3, we discuss the effects of baryons on halo populations, and in Section 4, we discuss our corrections. In Section 5, we discuss the environmental dependence of the HMF, and in Section 6, we discuss the halo mass–concentration relation. Finally, in Section 7, we summarize our results and conclusions.

For this analysis, we use three cosmological N-body simulations: Illustris (Genel et al. 2014; Vogelsberger et al. 2014a, 2014b; Nelson et al. 2015), IllustrisTNG (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Springel et al. 2018; Pillepich et al. 2018a, 2018b), and EAGLE (Springel 2005; Crain et al. 2015; Schaye et al. 2015; McAlpine et al. 2016; The EAGLE team 2017). Specifically, we use Illustris-1, TNG100-1, TNG300-1, and EAGLE RefL100N1504. Each of these hydrodynamic simulations has a corresponding DMO counterpart. A summary of the simulation parameters can be found in Table 1.

The Illustris and IllustrisTNG simulations were performed with the arepo code, while the EAGLE simulation was performed with gadget-3. While all three simulations model star formation, stellar evolution, gas cooling and heating, supernova feedback, black hole formation, and AGN feedback, gadget-3 and arepo contain different numerical hydrodynamical techniques, leading to large discrepancies in the galaxy populations they produce (Scannapieco et al. 2012). Additionally, though Illustris and TNG were performed with the same code, they vary in the strength of their feedback prescriptions. Illustris contains much stronger AGN feedback than either the TNG or EAGLE simulations (Weinberger et al. 2017), which leads to substantial differences in their results.

We choose to compare these three simulations because they are publicly available and all have DMO counterparts to their hydrodynamic simulations. Moreover, each simulation has adequate resolution and volume to allow us to examine its HMF between 10¹⁰ and 10¹⁴ h⁻¹ M_⊙. Finally, the three simulations have different baryonic physics prescriptions, which allows us to compare how variations in physics impact our results.

In this work, we utilize halo catalogs from each simulation at z = 0, 1, and 2. In all three simulations, halos were identified using a standard FoF algorithm (Davis et al. 1985) with a linking length of b = 0.2 times the mean interparticle separation. The FoF algorithm was run on the dark matter particles, and in the hydrodynamic simulations, baryonic particles were attached to the same FoF group as their nearest dark matter particle. Each FoF group has a mass, which we refer to in this work as M_FoF.

Additionally, for each FoF group, masses were calculated for several different SO definitions, which we will make use of in this paper. Specifically, these halo definitions are M_200b (the total mass of the group enclosed in a sphere whose mean density is 200 times the mean density of the universe at the time the halo is considered), M_200c and M_500c (the total mass of the group enclosed in a sphere whose mean density is 200 or 500 times the critical density of the universe at the time the halo is considered), and M_vir (the total mass of this group enclosed in a sphere whose mean density is Δ_c times the critical density of the universe at the time the halo is considered, where Δ_c derives from the solution of the collapse of a spherical top-hat perturbation from Bryan & Norman 1998).

For each halo definition and redshift, we only consider halos above 10¹⁰ h⁻¹ M_⊙. In Illustris, the number of halos above this threshold ranges from about 60,000 to 101,000, while in IllustrisTNG, the number ranges from about 70,000 to 120,000, and in EAGLE, the number ranges from about 53,000 to 85,000. The exact number of halos depends on the halo definition and redshift, but it is generally lowest for M_500c halos at z = 2 and highest for M_FoF halos at z = 1.

The most straightforward way to investigate the population of halos in a simulation is by looking at the HMF, which displays the abundance of halos as a function of mass. In Figure 1, we show the HMFs for M_200b halos at z = 0 for both the hydrodynamic and DMO versions of Illustris-1, TNG100-1, and EAGLE. Illustris is plotted in blue, TNG in green, and EAGLE in orange. It should be noted that the halo masses from the hydrodynamic simulations include both dark matter particles and baryonic particles. The hydrodynamic versions are plotted with dashed lines, while the DMO versions are plotted with solid lines. A corrected version of the DMO HMF is plotted with a dotted line and will be discussed further in the next section. In the residual panel, we plot the ratio of the hydrodynamic HMF to the DMO HMF (solid), as well as the ratio of the hydrodynamic HMF to the corrected DMO HMF (dotted) for each simulation.

Zoom In Zoom Out Reset image size

The residuals reveal sizable discrepancies between the DMO and hydrodynamic HMFs. In general, the hydrodynamic HMFs are shifted to lower masses in all three simulations, but this shift is mass-dependent and varies by simulation. In particular, the Illustris (blue) hydrodynamic HMF is consistently lower than the Illustris-Dark HMF above 10¹² h⁻¹ M_⊙, as well as below 10¹¹ h⁻¹ M_⊙. In TNG100 (green), the hydrodynamic HMF is lower than the DMO HMF between 10¹² and 10¹⁴ h⁻¹ M_⊙, as well as below 10¹² h⁻¹ M_⊙. In EAGLE (orange), the hydrodynamic HMF is lower than the DMO HMF at all halo masses below about 10¹³ h⁻¹ M_⊙.

In all three simulations, the discrepancies between the hydrodynamic and DMO HMFs reach the 20% level for much of the halo mass range. We emphasize that not only is this a significant difference, it is also not trivial to correct. It is not equivalent to a difference in halo definition, which would be a constant offset and would not vary by simulation; this effect exhibits a trend with halo mass and depends on the feedback implemented in the hydrodynamic simulation.

Our HMF results for Illustris are consistent with the findings of Vogelsberger et al. (2014b), who found that the HMF in Illustris is most affected at halo masses below 10¹⁰ and above 10¹² h⁻¹ M_⊙, where baryonic feedback processes (e.g., reionization and supernova/AGN feedback) are strongest, leading to a reduction in halo mass compared to their DMO counterparts. They also found that halos around 10¹¹ h⁻¹ M_⊙, where star formation is most efficient, tend to be more massive than their DMO counterparts.

Our HMF results for TNG100 are consistent with the findings of Springel et al. (2018), who found that baryons in the TNG simulation have a larger impact on low-mass halos and a smaller impact on high-mass halos compared to Illustris. The TNG has weaker AGN feedback than the original Illustris simulation, which leads to there being less discrepancy between the DMO and hydrodynamic HMFs at the high-mass end. The TNG also has a different wind model than Illustris, which leads to stronger feedback effects on lower-mass halos, leading to more discrepancy at the lower-mass end of the HMF in TNG.

Our HMF results for the EAGLE simulation are consistent with Schaller et al. (2015) and Desmond et al. (2017), who examined the differences between the halo masses in the EAGLE DMO and hydrodynamic runs and found the halos to be less massive, on average, in the hydrodynamic run. Desmond et al. (2017) found that, at low halo masses, stellar feedback in EAGLE strips baryons from the halo, which reduces the growth rate of the halo. At higher halo masses, stellar feedback becomes less effective, but AGN feedback is still capable of expelling baryons in all but the most massive halos.

Another way to visualize the discrepancy between the halo mass distributions from hydrodynamic and DMO simulations is to plot the fractional difference between their masses. In Figures 2 and 3, we plot this fractional difference for the M_200b (top) and M_500c (bottom) halos in Illustris, TNG100, and EAGLE at z = 0 and 2, respectively. Illustris halos are plotted in blue, TNG100 halos in green, and EALGE halos in orange. The x-axis is the DMO halo mass in units of h⁻¹ M_⊙, and the y-axis is the fractional difference of the hydrodynamic and DMO halo mass. In the left column, we use the total hydrodynamic halo mass (i.e., dark matter and baryonic particles) to calculate this fractional difference. In the middle column, we use only the dark matter component of the hydrodynamic halo to calculate this fractional difference, and we normalize by 1 minus the universal baryonic mass fraction, 1 − Ω_b /Ω_m , or (Ω_m − Ω_b )/Ω_m . In the right column, we use only the baryonic component (i.e., gas and star particles) of the hydrodynamic halo to calculate this fractional difference, and we normalize by the universal baryonic mass fraction Ω_b /Ω_m .

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In each panel of Figures 2 and 3, hydrodynamic and DMO halos are paired based on their ranked masses, rather than spatial positions or particle IDs. Thus, the most massive DMO halo is paired with the most massive hydrodynamic halo, and so on. In other words, we essentially "abundance match" the halos in the hydrodynamic and DMO simulations (we do this separately for each redshift and halo definition, as well as each column in the figures; therefore, one cannot compare an individual point from one column or row to the next). We adopt this procedure because the fractional mass differences calculated in this way represent the correction one would need to apply to the halo masses from one of the DMO simulations in order to exactly recover the global HMF from the corresponding hydrodynamic simulation. While this abundance-matching technique does not produce the exact same results as position or particle matching would, the overall trends seen are the same, although the scatter is drastically reduced by abundance matching (Vogelsberger et al. 2014b; Schaller et al. 2015; Springel et al. 2018).

Examining the top left panel of Figure 2, we can see that the same trends are present here as in the HMFs shown in Figure 1. The results at z = 0 for M_200b halos in Illustris (consistent with Vogelsberger et al. 2014b and Springel et al. 2018) indicate that stellar feedback slightly reduces the masses of the lowest-mass halos (by up to ∼ 10%), while star formation efficiency slightly increases the masses of halos around 10¹¹ h⁻¹ M_⊙, and AGN feedback severely reduces the masses of high-mass halos (by up to ∼20%). The results for TNG100 (consistent with Springel et al. 2018) indicate that stellar feedback reduces the masses of low-mass halos (by up to ∼15%), while star formation efficiency peaks at slightly higher masses than in Illustris (but is not quite so efficient as to increase halo masses), while AGN feedback is less strong than in Illustris and does not affect the highest-mass halos but reduces the masses of intermediate-mass halos by up to ∼10%. The results for EAGLE (consistent with Schaller et al. 2015) indicate that stellar feedback severely reduces the masses of low-mass halos (by up to ∼20%), while AGN feedback is similar to that in TNG and reduces the masses of intermediate-mass halos (by ∼10%) but does not impact the highest-mass halos.

Looking at the second and third columns of Figure 2, we can see that feedback has a much more extreme effect on the baryons in a halo than it does on the dark matter particles. Thus, most of the mass difference is due to a loss or gain of baryons, and not dark matter. For example, in the EAGLE simulation, there is an extreme lack of baryons at lower halo masses, while in the Illustris simulation, there is a significant lack of baryons at higher halo masses, but for both simulations, the amount of dark matter in each halo remains relatively constant.

The differences between the three simulations in the third column of Figure 2 indicate that these three hydrodynamic simulations disagree in terms of the baryon content of their halos. Precise observational constraints on the baryon fraction as a function of halo mass could, in principle, allow us to differentiate between these three simulations. Gonzalez et al. (2013) measured the baryons contained in both the stellar and hot-gas components for galaxy clusters and groups with M₅₀₀ > 10¹⁴ M_⊙ at z = 0.1 and found that the weighted mean baryon fraction for halos with M₅₀₀ > 2 × 10¹⁴ M_⊙ is 7% below the universal value when using a Planck cosmology. This is consistent with the results from the TNG and EAGLE simulations at the high-mass end but is not in agreement with the results from the Illustris simulation. The more significant differences between the three simulations, however, occur at lower halo masses. Eckert et al. (2017) examined the baryonic content of halos in the ECO and RESOLVE galaxy surveys (Moffett et al. 2015). While these results do extend to lower halo masses, they only include cold baryonic content, so we cannot make a direct comparison to the baryonic content in the three hydrodynamic simulations. In the future, more precise observational constraints on the baryon content of low-mass halos could allow us to rule out certain hydrodynamic simulations.

Comparing the top and bottom rows of Figure 2, we can see that the discrepancy between hydrodynamic and DMO halo masses is more extreme for M_500c halos than for M_200b halos, indicating that feedback has a stronger impact on the inner regions of a halo. This is most evident when looking at the baryonic component of Illustris M_500c halos (third panel of the second row of Figure 2), where feedback and star formation have an extreme effect on the percentage of baryons in a halo.

We can also see that the trends present in each simulation at z = 0 are quite different from those at z = 2, which are shown in Figure 3. For example, we can see clearly from the M_500c Illustris halos that at z = 2, stellar and AGN feedback has not kicked in yet, and star formation efficiency is quite strong, resulting in a strong overabundance of gas and star particles. Because of this, essentially all M_500c Illustris halos at z = 2 have a higher mass in the hydrodynamic simulation than in the DMO simulation (by ∼10%). However, by z = 0, stellar feedback has kicked in for halos at or below 10¹⁰ h⁻¹ M_⊙, while AGN feedback has kicked in for halos above 5 ×10¹¹ h⁻¹ M_⊙, resulting in these halos having lower masses in the hydrodynamic simulation than their DMO counterparts. Meanwhile, in TNG100 and EAGLE, it appears that stellar feedback has already kicked in by z = 2, but AGN feedback has not, leading to less of a reduction in halo mass at the high-mass end but a similar result at low masses.

Overall, Figures 2 and 3 emphasize the fact that the effect of baryonic physics on the HMF is to shift the HMF to lower masses. This shift can be as large as 25%. However, it is clear that this effect varies dramatically from one simulation to the next. This presents a problem for anyone using DMO halo catalogs to constrain the galaxy–halo connection or cosmology (Beltz-Mohrmann et al. 2020). The solution to this problem is not as simple as adding a parameter to the model, because it is a problem with the dark halo population itself. Additionally, the solution is not as simple as applying a constant offset to all DMO halo masses, because the discrepancy depends on the mass regime. Finally, once again, the solution depends on which hydrodynamic simulation is regarded as having the correct baryonic physics.

One possible solution is to apply a correction to each of the halos in a DMO simulation, so that the HMF better mimics the mass function from a hydrodynamic simulation. This correction should serve to adjust each DMO halo so that it has the mass of its corresponding hydrodynamic halo. We can use the fractional difference in halo mass shown in Figures 2 and 3 to identify a functional form for this correction.

To do this, we fit a polynomial to the fractional difference in halo masses between our hydrodynamic and DMO simulations (i.e., first column of Figures 2 and 3). The polynomial fits to each halo mass relationship were found using NumPy's polyfit function. After examining the effectiveness of several polynomials to correct the HMF, we found that a seventh-order polynomial was the lowest order that could accurately capture the halo mass trend and allow us to correct the DMO HMF. We have additionally looked at the Bayesian information criterion (BIC) for polynomials ranging from third to 12th order and found that the mean BIC (averaged over each halo definition and redshift combination) continues to decrease until we reach seventh order but does not continue to decrease significantly for higher orders. Thus, we decided that a seventh-order polynomial was the lowest order we could use for our fits and still accurately correct the HMF.

The corrected DMO halo masses (in units of 10¹⁰ h⁻¹ M_⊙) are given by

$$\begin{eqnarray}&&{M}_{{\rm{h}},\mathrm{corrected}}=(y+1)\times {M}_{{\rm{h}},\mathrm{DMO}},\end{eqnarray} \tag{ 1 }$$

where M_h,DMO is the unlogged original halo mass in units of 10¹⁰ h⁻¹ M_⊙, and

$$\begin{eqnarray}&&y={{ax}}^{7}+{{bx}}^{6}+{{cx}}^{5}+{{dx}}^{4}+{{ex}}^{3}+{{fx}}^{2}+{gx}+h,\end{eqnarray} \tag{ 2 }$$

where $x={\mathrm{log}}_{10}({M}_{{\rm{h}},\mathrm{DMO}})$ and a–h are the polynomial coefficients for a given simulation and halo definition. These fits (for M_200b and M_500c halos at redshifts zero and 2) are plotted with solid lines in the first column of Figures 2 and 3. The fits at redshift zero for all halo definitions are listed in Table 2. Additionally, fits for these same five halo definitions at redshifts 1 and 2 are available in machine-readable format.

Once again, this correction is based on "abundance matching" the halos between hydrodynamic and DMO simulations. Thus, applying this correction will assign to the most massive DMO halo the mass of the most massive hydrodynamic halo, and so on. After applying our corrections to each DMO simulation, we can examine our new corrected DMO HMFs. In Figure 1, we plot these corrected HMFs for M_200b halos at z = 0 with dotted lines. Looking at the residual panel, we can see that the corrections lead to significant improvement and almost perfectly reproduce the HMFs from the hydrodynamic simulations (deviations are less than 5%), which was precisely the goal of our abundance-matching method.

Unfortunately, the Illustris, TNG100, and EAGLE simulations do not contain any halos with masses above about 4 × 10¹⁴ h⁻¹ M_⊙ at z = 0. This upper limit decreases for z = 1 and 2 halos and is also somewhat dependent on halo definition. We do not make any assumptions about the masses of our halos beyond our data. Thus, our mass corrections should not be applied to any DMO halos above the limits given in Table 2. Rather, any DMO halo above our limits should be left unaltered. (This is very important; because our mass corrections are seventh-order polynomials, extending them beyond our mass limits will lead to very large changes in halo mass.)

It is noteworthy that in TNG100 and EAGLE, the mass correction is already almost zero at the high-mass end. However, in Illustris, for z = 0 halos, this is not the case. This means that applying the Illustris halo mass correction will lead to a slight discontinuity in the HMF. This discontinuity could be alleviated by, for example, extrapolating the fit at the high-mass end; however, because doing this would not be based on any data, we do not provide any extrapolations of our fits here. Additionally, our mass corrections should only be applied to halos with masses above 10¹⁰ h⁻¹ M_⊙. We do not present corrections for halos below 10¹⁰ h⁻¹ M_⊙ in this work due to the mass resolutions of the simulations that we use. (Once again, it is very important not to apply the corrections to any DMO halos below 10¹⁰ h⁻¹ M_⊙ due to the nature of the seventh-order polynomial fits.)

We recommend that any future work using halos from a DMO simulation do the following: for a given halo definition and redshift, apply at least one of our halo mass corrections to correct the HMF. Once the mass corrections are applied, the remaining analysis (e.g., HOD, conditional luminosity function (CLF), etc.) can be performed on the corrected halo catalogs. In this way, one can investigate the robustness of their results to changes in the HMF. Ideally, we recommend applying all three corrections (i.e., each correction based on Illustris, TNG100, and EAGLE), since the variation among them represents the theoretical uncertainty in baryonic physics.

We have created a python module for implementing our halo mass corrections, which is publicly available at https://github.com/gbeltzmo/halo_mass_correction. This module takes in an array of halo masses, a halo definition (M_200b , M_FoF, M_vir, M_200c , or M_500c ), a redshift (0, 1, or 2), and a simulation (Illustris, TNG, or EAGLE) and returns the corresponding corrected halo masses. If a given halo mass is outside the accepted mass range, the code will issue a warning and return the original (uncorrected) halo mass.

One question worth investigating is whether applying our mass corrections to a box with a very different resolution produces accurate results. To investigate this, we applied the TNG correction (based on the TNG100-1 simulation) to the TNG300-1-Dark box (which is about a factor of 8 lower in resolution than the TNG100-1-Dark simulation) to see if we could reproduce the HMF from the TNG300-1 simulation. Looking only at the 10¹⁰–10¹⁴ h⁻¹ M_⊙ mass range, this correction almost perfectly reproduces the HMF from TNG300-1, with all deviations less than 5% (compared to 17% deviations without the correction). This indicates that as long as the mass corrections are only applied within the appropriate mass regime, the corrections are accurate even when used with simulations of different volumes and resolutions.

A. O. Szewciw et al. (2021, in preparation) applied our mass corrections to M_vir halos (z = 0) from a large DMO simulation (Las Damas; McBride et al. 2009) and examined how their HOD parameter constraints on Sloan Digital Sky Survey galaxies varied with the different halo mass corrections. They found that for both of their luminosity samples (${M}_{r}^{-19}$ and ${M}_{r}^{-21}$), the different halo mass corrections lead to changes in all of their HOD parameters. The biggest changes are seen in $\mathrm{log}{M}_{\min }$ and $\mathrm{log}{M}_{1}$, and in most cases, the Illustris-based mass correction leads to the biggest change, although for the ${M}_{r}^{-19}$ sample, it is the EAGLE-based correction that leads to the biggest difference in $\mathrm{log}{M}_{\min }$. While it is to be expected that the halo mass corrections lead to changes in the mass parameters of the HOD, these changes are simulation-dependent and are not trivial to predict. The halo mass corrections ultimately do not lead to better χ² values for the best-fitting models, nor do they lead to the models being ruled out.

Our abundance-matching technique does not explicitly take the halo environment into account when correcting the halo masses. This means that while our mass corrections successfully reproduce the global HMF from our hydrodynamic simulations, they will not correct the conditional HMF if baryonic effects on halo mass are environment-dependent.

Fitting halo mass corrections in which halos are matched between DMO and hydrodynamic based on position or particle IDs would inherently take the halo environment into account. However, there are several issues with this technique. First, matching halos based on position or particle IDs does not guarantee that every halo in the hydrodynamic simulation has a match in the DMO simulation. Second, this method of matching introduces a significant amount of scatter into the hydrodynamic-to-DMO halo mass relationship. Therefore, if one were to use this relationship to correct halos from a large DMO simulation, one would have to either ignore the scatter (in which case, the result is essentially the same as the abundance-matching correction) or take the scatter into account by binning halos by mass and then drawing their corrected mass from a distribution within that bin. While this accounts for the scatter, it still does not account for the halo environment. Additionally, this is not the cleanest method for correcting DMO halo masses; when applied to any of the three simulations used in this work, this method does not successfully reproduce the correct global or conditional HMF from the hydrodynamic simulations.

One alternate possibility is to use our original abundance-matching technique but to separate the halos by environment. We can do this by measuring the large-scale environment around each of our DMO and hydrodynamic halos. We can then split our halos into those with high-density environments and those with low-density environments and subsequently abundance match between DMO and hydrodynamic, matching only halos within similar environments. For example, the most massive DMO halo in a high-density environment would be matched with the most massive hydrodynamic halo in a high-density environment, while the most massive DMO halo in a low-density environment would be matched with the most massive hydrodynamic halo in a low-density environment, and so on. This procedure will yield mass corrections that are guaranteed to recover the correct conditional HMF.

To do this for each of our simulations, we first have to measure the large-scale environment around each of our halos. To measure the halo environment, we calculate the total mass of halos within 5 Mpc spheres centered on each halo of interest (excluding the mass of the halo of interest itself). In other words, we do not sum up all particles but rather sum up the masses of all halos in the halo catalog whose centers fall within the sphere (excluding the main halo). We do not impose any lower mass limit on the halos included in this sum. We measure environments for all DMO halos above 10¹⁰ h⁻¹ M_⊙ and all hydrodynamic halos that are matched to them using the abundance-matching method described above.

We can thus define an environment factor δ for each halo, such that $\delta =({\rho }_{\mathrm{sphere}}/{\rho }_{\mathrm{box}})-1$, where ${\rho }_{\mathrm{sphere}}$ is the mass of halos in a 5 Mpc sphere around the halo divided by the volume of a 5 Mpc sphere, and ρ_box is the sum of all halo masses in the box divided by the volume of the box. We measure the halo environment in this way for all three of our simulations, hydrodynamic and DMO alike. We repeat this measurement for each of our different halo definitions and redshifts. We then split our halos into high- and low-density environments based on the median environment (δ_med) for that simulation. This is done separately for each halo definition and redshift and for the hydrodynamic and DMO simulations. Thus, each simulation/halo definition/redshift combination has a slightly different δ_med. Subsequently, DMO halos in high-density environments are abundance-matched with hydrodynamic halos also in high-density environments, and likewise for halos in low-density environments.

Shown in Figure 4 are the results of this environment-dependent abundance-matching technique for M_200b halos at z = 0 in Illustris-1 (blue), TNG100-1 (green), and EAGLE (orange). Halos in high-density environments are plotted in darker colors, and halos in low-density environments are plotted in lighter colors. Each of these relationships is once again fit with a seventh-order polynomial, which is plotted here with a solid line. The fits for M_200b and M_vir halos at z = 0, along with the δ_med and mass limit, are given in Table 3.

Zoom In Zoom Out Reset image size

We can see from Figure 4 that for each simulation, the relationship between hydrodynamic and DMO halos closely resembles that seen in Figure 2 for M_200b halos, but we do detect a difference between the high- and low-density environments. In Illustris and TNG100, the differences between environments appear for all DMO halos above 5 × 10¹⁰ h⁻¹ M_⊙, where hydrodynamic halos in high-density environments are slightly more massive relative to DMO halos than those in low-density environments. In EAGLE, below 2 × 10¹⁰ h⁻¹ M_⊙, the high- and low-density trends are the same. Between 2 × 10¹⁰ and 5 ×10¹¹ h⁻¹ M_⊙, the hydrodynamic halos in high-density environments are slightly more massive relative to DMO halos than those in low-density environments; at higher masses, this relationship is reversed. Additionally, for all three simulations, the highest-mass hydrodynamic and DMO halos are exclusively found in high-density environments, as expected.

While our original abundance-matching halo mass correction reproduces the global HMF from a given hydrodynamic simulation, our environment-dependent halo mass correction reproduces both the global and the conditional HMF. To examine whether it is important to implement this more complicated halo mass correction, we can see how well our different mass corrections reproduce the halo clustering found in hydrodynamic simulations by measuring the halo correlation function, ξ(r).

In Figure 5, we plot the halo correlation functions for M_200b halos at z = 0 in TNG300-1. We use the TNG300-1 box for this analysis because the smaller hydrodynamic boxes contain too few halos to accurately examine the clustering of halos. We calculate the halo correlation function in two different mass bins: halos greater than 10¹¹ h⁻¹ M_⊙ and halos greater than 10¹² h⁻¹ M_⊙. For the lower-mass sample, we measure ξ in 13 bins of separation r between 0.49 and 15 h⁻¹ Mpc, and for the higher-mass sample, we measure ξ in 10 bins of separation r between 1.07 and 15 h⁻¹ Mpc. We compute ξ with the blazing fast code corrfunc (Sinha & Garrison 2017, 2020). We measure ξ on halos from the following versions of TNG300: (i) the hydrodynamic simulation, (ii) the DMO simulation, (iii) the DMO simulation with the global abundance-matching halo mass correction, (iv) and the DMO simulation with the environment-dependent halo mass correction. In this figure, the corrections are done on a halo-by-halo basis; in other words, we do not use our polynomial fits to make these corrections, but rather we directly assign DMO halos the exact masses of the corresponding abundance-matched hydrodynamic halos. In this way, we can assess the ability of a given mass correction to reproduce the bias seen in the full-physics version of TNG300 without introducing any error due to our fits. In the bottom panel of Figure 5, we plot the ratio of each subsequent ξ (DMO, DMO with abundance-matching correction, or DMO with environment-dependent correction) to the hydrodynamic ξ.

Zoom In Zoom Out Reset image size

For both mass samples, the DMO halos do not exhibit the same correlation function as the hydrodynamic halos; the residual panel shows that the difference in ξ for the low- (high-) mass sample is about 8% (7%) at the smallest scales, and for both mass samples, it is about 4% for all scales above 2 h⁻¹ Mpc. After applying the global abundance-matching halo mass correction, the difference in ξ reduces to about 2% for both mass samples for all scales above 2 h⁻¹ Mpc (the difference is slightly larger for the high-mass sample). This difference reaches about 6% on the smallest scales for both mass samples. We emphasize that this lingering discrepancy in ξ after applying the abundance-matching mass correction is not due to the accuracy of our seventh-order polynomial fit, because in this figure, we are not using our fits but rather matching individual halos. In other words, at this point, the global HMF has been perfectly corrected, so the remaining differences in ξ must be due to other factors.

After applying the environment-dependent halo mass correction, the discrepancies in ξ shrink once again. For both mass samples, the discrepancy between the environment-dependent corrected DMO halo correlation function and the hydrodynamic halo correlation function is less than 1% for all scales above 2 h⁻¹ Mpc. At smaller scales, the environment-dependent correction is still an improvement over the original correction or no correction, although the discrepancy still reaches about 5% or 6% at the smallest scales we consider. This is due to the fact that halo exclusion occurs at smaller scales in the hydrodynamic simulation than in the DMO simulation. Because our halo mass corrections have no impact on the sizes or positions of halos, they cannot improve the halo correlation function until it reaches a scale where both the DMO and hydrodynamic simulations are not lacking halo pairs due to halo exclusion.

These results indicate that the impact of baryonic physics on the HMF is dependent on the environment of the halo, which in turn affects the ability of a DMO simulation to reproduce the halo clustering observed in hydrodynamic simulations. Our original abundance-matching halo mass corrections are able to reproduce the global HMF seen in hydrodynamic simulations and the halo clustering to within a few percent. However, if higher accuracy is needed, it is important to account for environment-dependent baryonic effects. Our environment-dependent halo mass corrections reproduce both the global and the conditional HMF from hydrodynamic simulations, and they reproduce the halo clustering to within less than 1%.

In Table 3, we provide our environment-dependent halo mass corrections for M_200b and M_vir halos in Illustris, TNG, and EAGLE at z = 0. Additionally, corrections for the remaining halo definitions and redshifts are available in machine-readable format, as well as in our python module (https://github.com/gbeltzmo/halo_mass_correction).

Figure 6 shows the results of applying these environment-dependent halo mass corrections to Illustris (left), TNG (middle), and EAGLE (right). Each panel shows the discrepancy between the hydrodynamic correlation function and that in DMO (dark), DMO with the original halo mass correction (medium), and DMO with the environment-dependent halo mass correction (light). In this figure, the mass corrections used are the fits given in Tables 2 and 3 for M_200b halos at z = 0. The top panels show the correlation function discrepancies for halos above 10¹¹ M_⊙, and the bottom panels show the same for halos above 10¹² M_⊙. (We do not show halos below 10¹¹ M_⊙ because in all three simulations, the environment-dependent corrections do not deviate from the original corrections below 10¹¹ M_⊙. Additionally, the clustering of low-mass halos is more complicated, and its investigation is beyond the scope of this paper.)

Zoom In Zoom Out Reset image size

In each simulation, the DMO halos exhibit some clustering discrepancy compared to the hydrodynamic halos. For the lower-mass sample (top panels), DMO underestimates the clustering by about 3%, on average, in Illustris and TNG and about 5%, on average, in EAGLE. In Illustris, the original mass correction actually makes the clustering slightly worse, while the environment-dependent correction is an improvement on small scales and a slight overcorrection on large scales. In TNG, the original mass correction provides little to no improvement over DMO, but the environment-dependent correction improves the clustering almost completely. In EAGLE, the original and environment-dependent corrections are in very close agreement, and both improve the clustering almost entirely compared to DMO.

For the higher-mass sample (bottom panels), the error bars are large due to a lack of high-mass halos. However, in each case, the DMO simulations underestimate the clustering of halos compared to the hydrodynamic simulations. This discrepancy is smallest in Illustris, and as a result, both the original and the environment-dependent mass corrections yield a slight overestimation of clustering. In TNG, the original mass correction provides little to no improvement, while the environment-dependent correction once again improves the clustering almost completely. In EAGLE, both the original and the environment-dependent corrections provide a slight improvement compared to DMO.

Applying these corrections to a large DMO simulation is slightly more complicated, because it requires calculating the large-scale environment for each halo of interest (only those above 10¹⁰ h⁻¹ M_⊙) and then separating these halos into "high-" and "low"-density environments based on the median environment δ_med and applying the corresponding mass correction. (In our module, an environment argument must be passed, which can be "all," "high," or "low," wherein the code assumes it is applying the correction to all halos, only halos in high-density environments, or only halos in low-density environments, respectively).

We have provided the δ_med that we found for each simulation and halo definition in Table 3. If applying our environment-dependent corrections, we suggest using the δ_med that we provide for a given correction, rather than finding a new δ_med that is specific to your sample. This is because our corrections are based on splitting halos into high- and low-density environments using a δ_med for all halos above 10¹⁰ h⁻¹ M_⊙. However, δ_med does increase with halo mass. Thus, for a halo catalog that only includes halos with masses above 10¹² h⁻¹ M_⊙, for example, the pertinent δ_med would be higher. This means that halos that we identified as being in high-density environments might now be identified as being in low-density environments. In this case, it would not be appropriate to apply our environment-dependent halo mass corrections to this sample. This is admittedly a limitation of our environment-dependent halo mass corrections. The effect of using the "wrong" value of δ_med is very slight, but we still recommend using the value of δ_med that we provide to ensure that the mass corrections are applied to the appropriate halos.

We once again provide upper mass limits for our halo mass corrections and emphasize the importance of not applying the corrections to any halos above these mass limits or below 10¹⁰ h⁻¹ M_⊙.

Halo modeling often uses the internal structure of halos to determine the placement of satellite galaxies within the halo (e.g., Zehavi et al. 2011; Zentner et al. 2014). This often involves assuming that the spatial distribution of satellite galaxies within halos follows a Navarro–Frenk–White (NFW) profile (Navarro et al. 1996), which includes a parameter for the concentration of the halo. Up to this point, we have investigated and quantified the impact of baryonic physics on the HMF, as well as the halo correlation function. We now examine the impact of baryonic physics on halo concentration.

The concentration c of a halo is defined as the ratio of the virial radius R_vir of the halo to the scale radius R_s (Navarro et al.1997). The relationship between halo mass and concentration has been previously studied in simulations, and it has been found that halo concentration has a weak power-law dependence on halo mass, with a slope of approximately −0.13 (Bullock et al. 2001). This relationship has a great deal of scatter and is predicated on the assumption that the density within a dark matter halo follows an NFW profile.

Ragagnin et al. (2019) investigated the dependence of halo concentration on mass and redshift in the Magneticum hydrodynamic simulations, and later, Ragagnin et al. (2021) examined the cosmology dependence of halo masses and concentrations in hydrodynamic simulations. Bose & Loeb (2021) examined the mass and concentration of host halos in IllustrisTNG as they relate to velocity dispersion. Additionally, Wang et al. (2020) studied the density profiles of early-type galaxies (ETGs) in IllustrisTNG at z = 0 and found that the profiles are steeper in the hydrodynamic simulation than their counterparts in the DMO simulation. They also found that the density profiles of the ETG dark matter halos are well described by steeper-than-NFW profiles.

In order to investigate the effect of baryonic physics on the halo mass–concentration relation in Illustris, TNG100, and EAGLE, we must first calculate the concentration for each halo (above 10¹⁰ h⁻¹ M_⊙) in these simulations. While the concentration is not directly provided in any of these simulations' halo catalogs, we can estimate it from the various halo mass definitions that are provided: M_200b , M_200c , and M_500c . The method works as follows. Assuming that each halo obeys an NFW profile, we can integrate the profile and use the virial mass and radius of the halo to determine the scale density ρ₀ of the halo, substituting R_vir/c for R_s :

$$\begin{eqnarray*}&&{M}_{\mathrm{vir}}={\int }_{0}^{{R}_{\mathrm{vir}}}4\pi {r}^{2}\rho (r){dr},\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\rho }_{0}=\tfrac{{M}_{\mathrm{vir}}{c}^{3}}{4\pi {R}_{\mathrm{vir}}^{3}[{ln}(1+c)-\tfrac{c}{1+c}]}.\end{eqnarray*}$$

We can then substitute this to solve for any other halo mass M_h and radius R_h (e.g., M_200b and R_200b ):

$$\begin{eqnarray*}&&{M}_{h}=\tfrac{{M}_{\mathrm{vir}}}{{ln}(1+C)-\tfrac{c}{1+c}}[{ln}(1+{{cR}}_{h}/{R}_{\mathrm{vir}})-\tfrac{{{cR}}_{h}}{{{cR}}_{h}+{R}_{\mathrm{vir}}}].\end{eqnarray*}$$

We can rearrange this equation to find the ratio of M_vir to M_h as a function of R_vir, R_h , and c:

$$\begin{eqnarray*}&&\tfrac{{M}_{\mathrm{vir}}}{{M}_{h}}=\tfrac{{ln}(1+c)-\tfrac{c}{1\,+\,c}}{{ln}(1+{{cR}}_{h}/{R}_{\mathrm{vir}})-\tfrac{{{cR}}_{h}}{{{cR}}_{h}+{R}_{\mathrm{vir}}}}.\end{eqnarray*}$$

For each of our simulations, we have three known values of M_h —M_200b , M_200c , and M_500c —and three corresponding values of R_h . To find the value of concentration for each halo, we loop over possible values of c between zero and 2000 and find the one that minimizes the sum of the squared fractional difference between the left- and right-hand sides of the previous equation over the three halo definitions. In other words,

$$\begin{eqnarray*}&&A=\tfrac{{M}_{\mathrm{vir}}}{{M}_{h}}\,\mathrm{and}\,B=\tfrac{{ln}(1+c)-\tfrac{c}{1\,+\,c}}{{ln}(1+{{cR}}_{h}/{R}_{\mathrm{vir}})-\tfrac{{{cR}}_{h}}{{{cR}}_{h}+{R}_{\mathrm{vir}}}},\end{eqnarray*}$$

and we find the value of c for each halo that minimizes

$$\begin{eqnarray*}&&{\left(\tfrac{{B}_{200b}-{A}_{200b}}{{A}_{200b}}\right)}^{2}+{\left(\tfrac{{B}_{200c}-{A}_{200c}}{{A}_{200c}}\right)}^{2}+{\left(\tfrac{{B}_{500c}-{A}_{500c}}{{A}_{500c}}\right)}^{2}.\end{eqnarray*}$$

We have tested this method of determining concentration on a DMO halo catalog for which we know the concentration values for each halo (calculated by rockstar; Behroozi et al. 2013). For this particular DMO halo catalog, using the given values of concentration led to a halo mass–concentration relation with a slope of −0.104 ± 0.004, while using our method of determining concentration led to a slope of −0.098 ± 0.006. Therefore, we can conclude that our method leads to the correct overall halo mass–concentration relation.

After confirming its accuracy, we applied this method to each of our hydrodynamic and DMO simulations to determine the concentration of each halo. In Figure 7, we plot log(c) as a function of log(M_vir) (in units of h⁻¹ M_⊙) for the Illustris (left), TNG100 (middle), and EAGLE (right) halos. The results for the DMO simulations are plotted as gray points, while the results for the hydrodynamic simulations are plotted in blue, green, and orange. We then bin the points by mass and plot the mean and standard deviation of these bins (in black for DMO or blue/green/orange for the hydrodynamic simulations). We subsequently fit a line to these means, which we plot with a dashed line. (Note that the dashed line is not connecting the points but rather is a fit to the points.) The slope and y-intercept for each simulation are given in Figure 7, along with the standard error in the slope.

Zoom In Zoom Out Reset image size

Based on these results, we can see that all three DMO simulations have halo mass–concentration relations that are consistent with each other and with what we expect from previous studies. The Illustris DMO simulation has a halo mass–concentration relation with a slope of −0.125 ± 0.004, while the TNG DMO simulation has a halo mass–concentration relation with a slope of −0.122 ± 0.005, and the EAGLE DMO simulation has a halo mass–concentration relation with a slope of −0.113 ± 0.010 (only slightly shallower than expected).

The Illustris hydrodynamic simulation has a halo mass–concentration relation with a slope of −0.239 ± 0.011, which is steeper than we would expect and not consistent with that of the DMO simulation. The TNG100 hydrodynamic simulation has a halo mass–concentration relation with a slope of −0.118 ± 0.012, and the EAGLE hydrodynamic simulation has a slope of −0.102 ± 0.014, both of which are consistent with their corresponding DMO simulations. In each of these cases, the scatter among both the DMO and hydrodynamic halos is quite large, but there does not appear to be a systematic offset between the DMO and hydrodynamic distributions for any simulation. Based on these results, we can conclude that in TNG and EAGLE, baryonic physics does not significantly impact the halo mass–concentration relation, while in Illustris, baryonic physics results in a slightly steeper halo mass–concentration relation compared to that in Illustris-Dark.

We would also like to know what the halo mass–concentration relation would look like in each of these simulations if we corrected the DMO halo masses but did not alter their concentrations. Another way of saying this is if someone were to apply our mass correction to a large DMO box but not change the halo concentration parameter, how would that impact the halo mass–concentration relation? To investigate this, we apply our abundance-matching halo mass correction to our DMO halos for each simulation and then plot the original DMO halo concentrations as a function of these corrected masses. In this case, we once again do a direct halo-by-halo correction rather than applying our seventh-order polynomial fit, in order to perfectly reproduce the mass function from the hydrodynamic simulation. Thus, the concentration (which we found via the method outlined above using the original halo masses and radii) of the most massive DMO halo is plotted against the mass of the most massive hydrodynamic halo (i.e., the "corrected" DMO halo mass), and so on. We then find the mean concentration in bins of mass and fit a new line to this relationship. This new fit is plotted as a dotted line in each panel of Figure 7, and the slope and y-intercept are given in the legend of each panel.

In each case, the "corrected" halo mass–concentration relation has a slope that is still consistent with what we would expect for halos obeying an NFW profile. In all three cases, applying the correction achieves a slope that is slightly closer to the hydrodynamic results. This is most significant in Illustris, where the slope becomes steeper, although it still does not agree fully with the hydrodynamic results. In TNG and EAGLE, the DMO, hydrodynamic, and corrected results are all consistent with one another.

Based on this, for any work relying on halo masses and concentrations from a DMO simulation, we believe that, depending on how one wants to employ the halo mass–concentration relation, one can either only correct the halo masses (and leave the concentrations untouched, with little impact on the results) or correct the halo concentrations as well using the slopes from the full-physics versions of Illustris, TNG, or EAGLE (shown in Figure 7). Either way, the impact of baryonic physics on the concentration–mass relation is likely to be small.

The implementation of baryonic physics in hydrodynamic simulations results in HMFs that are generally shifted to lower masses than those produced by DMO simulations. This is because stellar and AGN feedback removes baryonic particles (as well as some dark matter particles) from halos over time. This effect varies with halo mass; stellar feedback has more of an impact on lower-mass halos, while AGN feedback has more of an impact on higher-mass halos. Additionally, particularly efficient star formation can serve to increase the masses of some hydrodynamic halos (compared to their DMO counterparts).

In this work, we have quantified the relationship between the masses of halos in hydrodynamic simulations and those in corresponding DMO simulations. The impact of baryonic physics on the HMF depends on redshift, as well as halo definition. Additionally, because different hydrodynamic simulations contain different baryonic physics prescriptions, the halo mass discrepancy between hydrodynamic and DMO halos varies widely from one simulation to the next. Furthermore, the impact of baryonic physics on halo mass depends somewhat on the large-scale environment of the halo; in Illustris and TNG, halos in low-density environments are more impacted by baryonic physics, leading to greater discrepancies in their masses compared to their DMO counterparts. In EAGLE, halos in high-density environments exhibit a greater mass discrepancy between hydrodynamic and DMO simulations. This indicates that in hydrodynamic simulations in general, the strength of feedback (and its ability to impact halo mass) has a dependence on the density of the halo's environment.

We have found that this relationship as a function of DMO halo mass is well fit with a seventh-order polynomial. We provide these fits for Illustris, IllustrisTNG100, and EAGLE halos for M_200b , M_FoF, M_vir, M_200c , and M_500c halos at z = 0, 1, and 2. These fits are based on matching halos by mass (i.e., abundance matching) across hydrodynamic and DMO simulations. In other words, these are the corrections one would need to apply to the halos from Illustris-Dark, for example, to reproduce the HMF in the full-physics version of Illustris. We also provide these same fits after taking the halo environment into account, which can be used to reproduce the conditional as well as the global mass function. Furthermore, we have shown that these corrections for halo mass also reproduce the large-scale clustering of halos, though the environment-dependent corrections are required to achieve an accuracy better than 2%. Finally, we have shown that baryonic effects do not impact the halo concentration–mass relation substantially. Our halo mass corrections are publicly available as a python module at https://github.com/gbeltzmo/halo_mass_correction.

Any work relying on halo catalogs from DMO simulations (e.g., HOD modeling, conditional luminosity function modeling, stellar-to-halo mass relation modeling, etc.) could potentially be impacted by inaccuracies in the HMF. In particular, as these types of analyses start to be used more to constrain cosmological parameters, it is imperative that any conclusions are robust to changes in the HMF on the order of what we find in this work. For example, a DMO simulation created with a given set of cosmological parameters can produce an HMF that, after being adjusted with one of our mass corrections, resembles an HMF from a different cosmology. Thus, without understanding the uncertainty in the HMF due to baryonic physics, it may be challenging to distinguish between these cosmological models.

We recommend that any future work utilizing a halo catalog from a DMO simulation repeat its analysis after applying at least one of our halo mass corrections (and ideally more than one). This will provide a rough estimate of the systematic uncertainty in one's results due to baryonic effects on the HMF. We make no assumptions about which of the three hydrodynamic simulations used in this work produces the "correct" HMF but rather provide our halo mass corrections as a method for determining the robustness of one's modeling results to changes in the HMF.

We acknowledge the Virgo Consortium for making their simulation data available. The EAGLE simulations were performed using the DiRAC-2 facility at Durham, managed by the ICC, and the PRACE facility Curie based in France at TGCC, CEA, Bruyères-le-Châtel. Some of the computational facilities used in this project were provided by the Vanderbilt Advanced Computing Center for Research and Education (ACCRE). This project has been supported by the National Science Foundation (NSF) through award AST-1909631. This research has made use of NASA's Astrophysics Data System, as well as python (https://www.python.org/), the ipython package (Pérez & Granger 2007), scipy (Jones et al. 2001), numpy (Van Der Walt et al. 2011), and matplotlib, a Python library for publication-quality graphics (Hunter 2007). These acknowledgments were compiled using the Astronomy Acknowledgment Generator (http://astrofrog.github.io/acknowledgment-generator/). We especially thank Jolanta Zjupa for her help and our anonymous referee for constructive comments that helped us improve this paper.