The quasi-adiabatic relaxation of haloes in the IllustrisTNG and EAGLE cosmological simulations Free

Abstract

1 INTRODUCTION

A detailed understanding of the formation and evolution of galaxies and their interaction with their environment is a pressing open problem. In the Lambda-cold dark matter (ΛCDM) picture of the Universe, the primary environment of any galaxy is provided by a ‘halo’ of dark matter surrounding it (e.g. White & Rees 1978). Therefore, understanding the dynamical behaviour of these dark haloes during the evolution of the galaxies they host is a key ingredient needed for building a complete picture of galaxy evolution.

Dark haloes form through the gravitational collapse of overdensities that developed from tiny fluctuations in the initial distribution of matter (e.g. Press & Schechter 1974; for a review see Cooray & Sheth 2002). The properties of dark haloes identified in gravity-only cosmological N-body simulations have been extensively studied in the literature. For example, while these haloes are known to be triaxial (Frenk et al. 1988), their sphericalized mass profiles are found to have a universal form (Navarro, Frenk & White 1996b, 1997, hereafter, NFW; see also Navarro et al. 2010). However, as galaxies and clusters of galaxies form within these haloes from various non-gravitational baryonic interactions, their gravitational coupling to the dark matter can affect the spatial distribution and evolution of the latter. Understanding this response of a halo’s dark matter content to the baryons it hosts is then critical for understanding the coupled evolution of haloes and galaxies.

Early work treated this response using adiabatic invariants (Zel’dovich et al. 1980; Barnes & White 1984; Blumenthal et al. 1986; Ryden & Gunn 1987). Using simplifying assumptions such as spherical symmetry, no shell crossing, angular momentum conservation with circular orbits for dark matter particles, Blumenthal et al. (1986) derived a simple formula that quantifies the adiabatic relaxation of the dark matter mass profile in terms of the final baryonic distribution (we discuss this in detail later). Besides the change in their mass profiles, dark haloes can also become more spherical as a result of galaxy formation (Dubinski 1994). Currently, the most robust technique to understand the consequences of gas assembly and galaxy formation on dark matter structure is the use of high-resolution cosmological hydrodynamical (zoom) simulations, using ‘sub-grid’ recipes for modelling very small-scale astrophysics such as feedback from stellar/supernovae activity or the effects of active galactic nuclei (AGN) (see e.g. OWLS; Schaye et al. 2010, Illustris; Genel et al. 2014; FIRE; Hopkins et al. 2014; EAGLE; Schaye et al. 2015; Horizon-AGN; Kaviraj et al. 2017; SIMBA; Davé et al. 2019; IllustrisTNG; Nelson et al. 2019a). In such simulations, the response of the dark matter to the presence of baryons in a halo can be ascertained by comparing a halo in the full hydrodynamical simulation to a matched ‘partner’ halo in a collisionless gravity-only simulation performed using the same initial random fluctuations. Using this technique, it was found that a simple adiabatic contraction model like that of Blumenthal et al. (1986) is an inaccurate description of the response in a variety of simulations (see e.g. Gnedin et al. 2004; Gustafsson, Fairbairn & Sommer-Larsen 2006; Pedrosa, Tissera & Scannapieco 2010; Tissera et al. 2010; Duffy et al. 2010; Abadi et al. 2010; Teyssier et al. 2011; Dutton et al. 2016; Artale et al. 2019; Forouhar Moreno et al. 2022). Even in a simpler setting, where star formation and feedback effects are ignored, Abadi et al. (2010) found that the halo responds to the condensation of baryons to the centre by becoming more spherical and compact, while the change in its mass profile is found to be significantly less than the prediction of the idealized adiabatic relaxation model.

Several authors have investigated whether discarding some of the assumptions of the idealized model of Blumenthal et al. (1986) can help reconcile with the simulation results. Gnedin et al. (2004) considered non-spherical orbits for the dark matter particles and suggested a simple modification to the original model. This modified empirical formula shows wide variation in its parameters across haloes from different simulations and at different redshifts (Gustafsson et al. 2006; Duffy et al. 2010). On the other hand, Sellwood & McGaugh (2005) accounted for these random motions within the halo using invariant action integrals, following the method described in Young (1980). This model gives a reasonably good approximation of the response even in modern hydrodynamical simulations (Callingham et al. 2020); however, making predictions using this model requires access to orbital phase space information of the halo, which may not always be feasible.

Physically, one expects that the overall response of the halo is mediated by a combination of different astrophysical processes that occur in the galaxy. Feedback processes are known to reduce the contraction of the halo significantly; e.g. supernova-driven winds can completely transform the inner density profile of the dark matter halo (Navarro, Eke & Frenk 1996a). This may be the key in reconciling the observation of dark matter cores at the centre of various galaxies with the cuspy haloes found in gravity-only ΛCDM simulations (see Pontzen & Governato 2014, for a review). However, such feedback effects do not always produce dark matter cores from cusps, rather, this can depend on the amount of gas ejected, the mass-loss time-scale and the frequency of starburst events (see e.g. Ogiya & Mori 2011, 2014; Pontzen & Governato 2012; Benítez-Llambay et al. 2019). In massive haloes hosting galaxy groups or clusters, while the formation of powerful AGN in the central galaxy can strongly suppress star formation, it can still significantly reduce the adiabatic contraction of the halo (Teyssier et al. 2011). Moreover, the fluctuation in gravitational potential due to such feedback can expel the dark matter from the inner halo producing inner cores (Martizzi et al. 2012).

Different aspects of the halo response, such as the change in its mass profile, shape, phase space distribution, and substructure population, have been explored to date in a variety of hydrodynamical simulations (see e.g. Kazantzidis et al. 2004; Debattista et al. 2008; Di Cintio et al. 2014; Schaller et al. 2015; Zhu et al. 2017; Chua et al. 2017, 2019, 2021; Cataldi et al. 2021; Cautun et al. 2020; Freundlich et al. 2020; Riggs et al. 2022). Understanding the nature of this response (or ‘baryonic backreaction’) is critical in building accurate and robust models of halo shapes and sizes, for use in interpreting the results of upcoming large-volume surveys (Schneider & Teyssier 2015; Chisari et al. 2018; Aricò et al. 2021; see also Velliscig et al. 2014; Harnois-Déraps et al. 2015; Mead et al. 2015) as well as the detailed prediction of rotation curves and related statistics (Paranjape & Sheth 2021a, b). The goal of this work is to perform a systematic statistical study of this dark matter response in high-resolution hydrodynamical simulations incorporating realistic feedback and quantify it using simple analytical forms, including the sensitivity of this response to halo-centric distance and halo and galaxy properties. To this end, we use the publicly available suites of simulations from the IllustrisTNG and EAGLE projects.

The rest of this article is organized as follows. In Section 2, we briefly discuss the simulations used and the techniques employed in this work. In Section 3, we present our results. These are expected to be relevant for a variety of problems; for example, the change in the dark matter density profile of the halo caused by the galaxy formation affects the rotation curve of the galaxy. We discuss such applications in Section 4 and conclude in Section 5. Throughout, we will use log and |$\ln$| to denote the base-10 and natural logarithm, respectively.

2 SIMULATIONS AND TECHNIQUES

In this section, we describe the numerical simulations used in this work and the generation of matched halo catalogue as well as the methods employed in our study of halo response to galaxy formation.

2.1 Simulations

We use two suites of cosmological hydrodynamical simulations, namely the IllustrisTNG (Nelson et al. 2019a) and the EAGLE (The EAGLE team 2017) suites that were run with different (yet state-of-the-art) prescriptions for baryonic processes involved in galaxy formation. We isolate the effects of galaxy formation process on dark matter halo by comparing these hydrodynamical simulations with their corresponding gravity-only runs that evolve the same initial cosmological volumes but treating baryons as collisionless.

2.1.1 IllustrisTNG simulations

This suite of simulations by the TNG collaboration project were run with the arepo code (Weinberger, Springel & Pakmor 2020), simulating the hydrodynamics using moving mesh defined by Voronoi tessellation (Springel 2010). This follows an updated model of galaxy formation based on the results from the original Illustris simulation; including all major baryonic processes like cooling, star formation, stellar feedback, and AGN feedback (see Weinberger et al. 2017; Pillepich et al. 2018a, for details). In addition to these, this updated model of TNG includes a cosmic magnetic field. The TNG suite has three cosmological boxes TNG50, TNG100, and TNG300 with a periodic box size of 35 h⁻¹Mpc, 75 h⁻¹Mpc, and 200 h⁻¹Mpc, respectively, assuming a cosmology consistent with results from Planck Collaboration XIII (2016). Initial conditions of these simulations were generated with the Zel’dovich approximation (Zel’dovich 1970) at redshift z = 127 using the n-genic code (Springel 2015), for the realization selected by χ²-minimization for cumulative halo mass function. Results from this simulation have been compared against observations and found to have reasonably realistic galaxies (see Nelson et al. 2018; Pillepich et al. 2018b; Springel et al. 2018; Naiman et al. 2018; Marinacci et al. 2018; Pillepich et al. 2019; Nelson et al. 2019b). In order to statistically study the response of a wide range of haloes to galaxy formation, we use the highest resolution runs of all three cosmological boxes; while the smallest box TNG50 provides sufficient resolution to resolve low-mass haloes, we need the large box TNG300 to get sufficient number of cluster-scale haloes. Throughout this work, we use the redshift z = 0.01 data from IllustrisTNG for both hydrodynamical as well as the corresponding gravity-only runs.

2.1.2 EAGLE simulations

This suite of cosmological simulations was performed with a modified version of the gadget-3 code that evolves the gas using smoothed particle hydrodynamics (Springel 2005). Initial conditions for these simulations were generated using the ic_2lpt_gen code following Jenkins (2010). From this suite, we use the z = 0 data of L0025N0752 and L0100N1504 boxes that were run with EAGLE’s reference model of galaxy formation, along with their corresponding gravity-only runs. This reference model includes sub-grid prescriptions for various baryonic processes like cooling, star formation, and feedbacks (see Schaye et al. 2015; Crain et al. 2015, for more details). While this model differs from that of IllustrisTNG in many aspects (e.g. it does not include a cosmic magnetic field), this suite of simulations has also produced realistic galaxies (see e.g. Sawala et al. 2015; Furlong et al. 2015; Trayford et al. 2015)

2.1.3 Halo and galaxy properties

Along with the particle data, both these simulation projects provide a catalogue of Friends-Of-Friends (FOF) group haloes, found with a linking length of 0.2 times the interparticle spacing (see McAlpine et al. 2016; Nelson et al. 2019a, for specifics). For each halo, we take the comoving position of the minimum of the gravitational potential within that halo as its centre, and define its ‘virial’ radius R_vir as the halo-centric radius enclosing a total density of 200 times the critical density of the Universe: R_vir ≡ R_200c. Throughout, we consider the total mass M of each halo to be the mass enclosed inside its virial radius, i.e. M ≡ M_200c. In addition, we have a catalogue of gravitationally bound substructures identified by the subfind code (Springel et al. 2001) within those FOF haloes. A single FOF group can have more than one subhalo, the one containing the central particle is considered as the central subhalo associated with that FOF halo.

2.2 Halo matching

To study how a dark halo responded to the galaxy forming in it, we need to first reliably match the catalogue of haloes found in the hydrodynamic simulation (which includes galaxy formation physics) to those found in its gravity-only run. For various numerical reasons, a given halo in the hydro run may not have a true match in the halo catalogue of gravity-only run. So, we first try to obtain an exhaustive catalogue of matched haloes that will be used to build a statistical description of this halo response.

2.2.1 Matching procedure

We match the haloes using the particle data associated with the haloes; while the mass of each dark matter particle differs between the hydrodynamical and gravity-only runs, the number of dark matter particles within the same initial periodic box is the same for each of the five pairs of simulations that we consider in IllustrisTNG and EAGLE. So, a given particle in a hydro simulation has originated from the same region as the particle of the same ID in its corresponding gravity-only run. For any given pair of haloes, with one in the hydrodynamical simulation and the other in the corresponding gravity-only run, we define the matching fraction of each of those two haloes (with respect to the other) as the fraction of its dark matter particles that are also present in the other halo. Below, we describe how we use these matching fractions to decide if the given pair can be considered a valid matched pair.

Computing the matching fractions for every pair of haloes is computationally expensive with O(n²) for millions of haloes in the catalogue. We decrease the complexity to O(n) by supplying, for each halo in the hydro simulation, an ordered list of most probable match candidates in the gravity-only run. These candidate lists of gravity-only haloes are generated and ranked based on the spatial positions of the haloes and their masses using a KD-tree based neighbour finding algorithm, implemented using scipy.spatial.KDTree. For each halo in the hydro run, we test if the matching fraction of this halo with respect to any of the haloes in its match candidate list exceeds the value of 0.5. This ensures that at most one halo is selected as a match for each of the hydrodynamical haloes, so that our matched catalogue of halo pairs will be a subset of the source catalogue without repetitions. While we match as many haloes as possible, it is also important to ensure that false matches do not plague our study. Based on the results in Appendix  A, we therefore additionally require that in a valid matched pair, the gravity-only FOF halo must also have a matching fraction of more than 0.5 with respect to the hydrodynamical halo. Our FOF based matching can be compared with Lovell et al. (2018) where a similar matching algorithm has been followed for central subhaloes.

2.2.2 Halo pair catalogue

We generate a matched catalogue of haloes for each of the five simulations studied in this work, including FOF haloes resolved with more than 1000 particles.¹ The fraction of hydrodynamical haloes that fails to be part of the matched catalogue is shown in Fig. 1 in bins of halo mass for the IllustrisTNG simulations. In the mass range in which the halo samples are selected for this work (see Section 3.2), more than 96 per cent of the haloes in hydrodynamical simulation have been assigned a match in the gravity-only run. The small fraction of unmatched haloes primarily reside in dense environments where our algorithm presumably fails due to the inherent issues with 3D FOF algorithm in dealing with mergers. Similar results hold for the EAGLE simulations as well. For illustrative purpose, a visual representation of two randomly chosen halo pairs, one each from IllustrisTNG and EAGLE, is shown in Fig. 2.

Figure 1.

Fraction of haloes in the TNG hydrodynamical simulations that have not found a match is shown as a function of mass M_200c. coloured vertical bands represent the mass range relevant for this work in each of the three TNG simulation boxes.

Open in new tabDownload slide

Figure 2.

Visually inspecting two matched FOF halo pairs, one each from IllustrisTNG (top row) and EAGLE (bottom row) using 2D-projected dark matter density field around the centre of the hydrodynamical halo in a thick slice. The left-hand panel shows the halo in the hydrodynamical simulation and the right-hand panel shows the corresponding matched pair in the gravity-only simulation. The black circle shows the virial boundary of the gravity-only halo and blue circle shows that of the hydrodynamical halo. In both cases, the hydrodynamical halo is noticeably more spherical and compact than its gravity-only counterpart, with a spatially offset centre. See the text for a discussion.

Open in new tabDownload slide

2.3 Methods to study the halo response

The response of dark matter halo to galaxy formation has primarily two aspects, contraction or expansion of the halo towards the centre and a change in its triaxial shape. In this work, we study the former aspect of the halo response, by focusing on spherically averaged mass profiles. The illustrative haloes shown in Fig. 2 become more compact and spherical in the hydrodynamical simulation that includes galaxy formation. Also notice that there is an offset between the centre-of-potential locations of matched pairs of haloes. These offsets are likely correlated with the halo tidal environment and will be interesting to follow-up in future work.

2.3.1 Mass profiles

The overall expansion and/or contraction of dark matter in response to galaxy formation can be studied through the differences in spherically averaged mass profiles between matched haloes. For the dark matter, these radial profiles are obtained by adding up the mass of all dark matter particles contained within concentric spherical shells. In addition to these, we also need baryon mass profiles in modelling the dark matter response. While stellar mass profiles are computed in a similar fashion as dark matter, for the gas mass profiles, we use a Gaussian kernel to assign mass enclosed to each of the spherical shells.² The width of this Gaussian kernel was taken to match the SPH smoothing length for the EAGLE simulation, whereas for IllustrisTNG, we use the cube root of the Voronoi cell volume to define the kernel smoothing scale. We have tested that our results are robust to differences in the choice of this kernel.

2.3.2 Quasi-adiabatic relaxation model

2.3.3 The relaxation relation

Our focus in this work is to characterize the relaxation relation equation (4) as a function of halo and galaxy properties over a wide dynamic range; e.g, we would like to ask whether equation (5) is a good description of this relation. To study this, we must extract this relation for individual haloes in hydrodynamical simulations. For a given hydrodynamical halo in the matched catalogue, we can obtain this relaxation relation by considering its matched halo in the gravity-only run to represent its unrelaxed state. We find it convenient to work with r_f as a control variable. In this case, the values of r_i, M_i(r_i) and M_f(r_f) must be obtained from the matched halo pair, which can be done as follows.

Thus, for any shell defined by its relaxed radius r_f, we can obtain both the relaxation ratio r_f/r_i and the mass ratio M_i/M_f from its unrelaxed radius computed from mass profiles. Hence, we can obtain the relaxation relation by placing multiple concentric shells around the halo all the way to its virial radius R_vir. In Appendix  B, we have tested this algorithm on mock halo + galaxy systems generated with fixed known relaxation relations.

3 RESULTS

3.1 Relaxation of haloes in the IllustrisTNG

We find that the relaxation relation estimated as described in Section 2.3 varies widely across haloes in the matched catalogue. In Fig. 3, we show the relaxation relation for four different samples of haloes selected by their unrelaxed mass from the IllustrisTNG simulations. The first two samples are from TNG50, with masses M ∼ 10^11.5  and 10¹²h⁻¹M_⊙, respectively. Similarly the other two samples are from TNG100 and TNG300 with masses M ∼ 10^12.5 and 10^13.5 h⁻¹M_⊙, respectively. The relaxation relations of few individual randomly chosen haloes from each sample are shown by grey lines; we also show stacked relaxation relations for each of the sample (see below for measurement details). The quasi-adiabatic relaxation model equation (5) with q = 0.68 and q = 0.33 is shown by the dot-dashed and dashed purple lines, respectively, in each panel. The value q = 0.68 was proposed by Schneider & Teyssier (2015) as being a reasonable description of cluster-sized haloes, while Paranjape & Sheth (2021b) argued that q = 0.33 leads to a good description of the radial acceleration relation of Milky Way-sized spiral galaxies (see their Appendix A1). We will use these two models as reference points in the comparisons below. Since the samples shown are representative of the haloes in IllustrisTNG over a large mass range, it is clear that equation (5) with a constant q does not work for the majority of haloes in IllustrisTNG. Similar results hold for EAGLE haloes as well. This motivates a systematic study of the relaxation relation as a function of halo mass and other properties.

Figure 3.

Relaxation relation for 4 different sample haloes selected by mass from the IllustrisTNG simulations. The large coloured circles denote the stacked relaxation ratio and total mass ratio at 20 different shells, whose radii is indicated by the colour. Small coloured markers joined by grey lines show the relaxation relation of a few randomly chosen individual haloes in each of the samples. The black curves denote the radius-independent stack of relaxation relation for each sample (see the text). The quasi-adiabatic relaxation model equation (5) with q = 0.68 and q = 0.33 are shown by the dot–dashed and dashed purple lines, respectively, in each panel.

Open in new tabDownload slide

For each of the four samples selected by halo mass, we compute the relaxation ratio r_f/r_i and the enclosed total mass ratio M_i/M_f at 20 concentric spherical shells for all haloes in the sample. We take the largest shell at the relaxed virial radius r_f = R_vir, while the remaining 19 shells are taken at fixed values of r_f/R_vir for each of the halo. This allows us to stack the relaxation relation by simply taking the mean and standard deviation of the relaxation ratio and mass ratio at each of the 20 shells. While the physical size of the shell differs from halo to halo, we ensure that the smallest shell has a radius of at least 10 times the force smoothing length of the simulation. In Fig. 3, we show this stacked relaxation relation in large coloured markers, where the colour denotes the relaxed radius of the shell; and the error bar shown in red corresponds to the statistical error in the estimate of the mean value. By comparing with the small markers of the same colour, we can see that there is a significant scatter not only in the relaxation ratio but also in the mass ratio at fixed r_f/R_vir across haloes in each sample. To assess the level of systematic error introduced by our default choice of stacking technique, we also tested an alternate stacking definition, wherein we interpolate the relaxation relation of individual haloes to obtain the relaxation ratio at fixed values of mass ratios and stack them by ignoring the value of corresponding relaxed radii. However, this stacking method ignores radius information completely; we discuss the consequence of this later in Section 3.3.

3.2 Trend in relaxation relation with halo mass

As can be already noted in Fig. 3, the relaxation relation shows very different behaviour at different mass scales. In this section, we focus on the stacked relation (using our default stacking definition) and study how it varies as a function of unrelaxed halo mass. For this, we consider nine mass bins starting from log (M/h⁻¹ M_⊙) = 10 to 14 in steps of 0.5 dex. We list the colour labels used for these mass bins in Fig. 4; this colour-coding will be used in all subsequent plots. None of the five simulations considered, simultaneously provides a sufficiently large sample of cluster-scale haloes and well-resolved low-mass haloes. In the IllustrisTNG suite, we use the smallest box TNG50 to study haloes with mass 10¹⁰ h⁻¹M_⊙ < M < 10¹² h⁻¹M_⊙, whereas we use TNG100 and TNG300 to study haloes with mass 10¹¹ h⁻¹ M_⊙ < M < 10^12.5h⁻¹M_⊙ and 10¹²h⁻¹M_⊙ < M < 10¹⁴h⁻¹M_⊙, respectively. At those mass bins where multiple IllustrisTNG boxes provide halo samples, the smaller box provides a smaller sample but with better resolution. For computational ease, we limit the size of each sample to be ≤500 haloes, as we find that the statistics are well converged with this number.

Figure 4.

Representative colours we use to denote each of the halo mass bins.

Open in new tabDownload slide

By repeating the procedure described in Section 3.1, we obtain both the fixed-radius (default) stack and radius-independent (alternate) stack of the relaxation relation for each of the halo samples taken from IllustrisTNG at these nine mass bins (see left-hand panel of Fig. 5). For reference, note that the case of no relaxation would correspond to a horizontal line at unity in this plot. Relaxation is strongest for Milky Way-scale haloes, as indicated by the small values of the relaxation ratio for M ∼ 10¹² h⁻¹ M_⊙; we discuss the physical implications of this result later. Note that the simple quasi-adiabatic relaxation model equation (5) with q = 0.68 used in Schneider & Teyssier (2015) fails to explain the relaxation relation for any of the halo masses considered; however this model with q = 0.33 is reasonably close to the relaxation relation at M ∼ 10¹³ h⁻¹ M_⊙. And while the quadratic model proposed by Abadi et al. (2010) matches with the relaxation relation of 10^12.5h⁻¹M_⊙ haloes in IllustrisTNG, this is possibly a coincidence given that this model was built using zoom simulation of haloes in the mass range 10^11.5–10¹² h⁻¹ M_⊙, which shows a very different relaxation relation in IllustrisTNG.³

Figure 5.

The stacked relation between relaxation ratio and mass ratio as a function of halo mass in IllustrisTNG (left-hand panel) and EAGLE (right-hand panel) simulations. Here, the points and solid lines represent two different stacking methods as in Fig. 3. The colour-coding follows Fig. 4.

Open in new tabDownload slide

We also take six samples of haloes from the EAGLE simulation, in mass bins log (M/h⁻¹ M_⊙) = 10.5, 11, 11.5 from the small, high-resolution L25 box and in mass bins 12, 12.5, 13 from the main L100 box. Here too, the q = 0.68 model fails for all masses, but q = 0.33 model works reasonably for M ∼ 10¹³ h⁻¹M_⊙ haloes (see the right-hand panel of Fig. 5). We find that, despite having a different galaxy formation model, the relaxation relation for haloes found in the primary EAGLE run L100 is consistent with the results from IllustrisTNG. IllustrisTNG samples reach lower values of the relaxation ratio and mass ratio than EAGLE because of the better resolution available. For M₂₀₀ = 10¹² h⁻¹ M_⊙, the mean relaxation relation shown in Fig. 5, does not seem to be very different between IllustrisTNG and EAGLE L100, at least not anymore than the difference between different boxes of the IllustrisTNG. However, the haloes from EAGLE L25 simulation shows a unique behaviour where the relaxation ratio increases with decrease in mass ratio in the innermost regions. We expect that this might be due to the fact that the EAGLE reference model required recalibration at this higher resolution. In a future work, we will explore how the dark matter response depends on such variations in the baryonic prescription.

3.3 Parametrized model of quasi-adiabatic relaxation

In this section, we model the relaxation relations discussed above, with a focus on conveniently quantifying this response across a wide range of halo masses.

3.3.1 Expectations from simulation measurements

In both IllustrisTNG and EAGLE, for all masses other than 10¹³h⁻¹ M_⊙, the simple quasi-adiabatic relaxation model equation (5) fails to explain the measured relation with any value of q. As seen in Fig. 5, an important aspect of this mismatch is caused by the model’s requirement that shells which hold their baryonic mass fixed (i.e. for which M_i/M_f = 1) must necessarily also hold their radius fixed (r_f/r_i = 1), and vice-versa. The measurements, however, show substantial offsets in the relaxation ratio from unity for shells with M_i/M_f = 1, and also substantial offsets in the mass ratio from unity for shells with r_f/r_i = 1, across nearly the entire range of halo mass. One way of understanding this effect physically is due to feedback-related baryonic outflows: a particular shell which maintains its radius after relaxation (r_f/r_i = 1), could still lose its baryonic mass due to outflows, resulting in M_i/M_f > 1 (Forouhar Moreno et al. 2022). Alternatively, the interplay between cooling-related condensation (which increases baryonic mass in a given shell) and feedback-related outflows (which decrease baryonic mass) could result in a situation where the baryonic mass after relaxation retains its initial value despite an overall relaxation, e.g. due to approximate angular momentum conservation, leading to r_f/r_i < 1 while M_i/M_f = 1. These trends are visible in Fig. 5 for haloes with M < 10¹³h⁻¹ M_⊙. Fig. 3 shows that the former trend (M_i/M_f > 1 when r_f/r_i = 1) occurs in the halo outskirts (r_f ∼ R_vir) and the latter (r_f/r_i < 1 when M_i/M_f = 1) in the inner halo (⁠|$r_f\lesssim 0.3\, R_{\rm vir}$|⁠), for M < 10¹³ h⁻¹M_⊙. On the other hand, more massive haloes show little to no relaxation in both inner halo (where there is net baryonic inflow, M_i/M_f < 1) and outer halo (where there is net baryonic outflow, M_i/M_f > 1).

However, the above holds true only when the dark matter particles are in circular orbits. When galactic processes lead to changes in the baryonic mass profile, even the dark matter particles in circular orbit can go into elliptical orbits (see e.g. Sellwood & McGaugh 2005). For example, when there is a sudden expulsion of gas due to feedback events, the total mass enclosed decreases and the particles start moving radially outward. During this period the mass ratio M_i/M_f can become greater than one and still have no relaxation (i.e. r_f/r_i = 1) as discussed in the start of this section.

While this extended linear model can describe the relaxation relation at few other halo masses (see for example 10^12.5h⁻¹ M_⊙ haloes in both left and right-hand panels of Fig. 5), even this model fails at many halo masses. Moreover, we have checked that there is no simple polynomial model favoured by standard information criteria such as AICC (Liddle 2007) to describe the relaxation relation at all masses. Rather, we find that, if we simply elevate q₀ and q₁ in equation (7) to functions of r_f/R_vir, then this model can be applied at all the mass scales that we consider. For this, we need the relaxation relation at fixed relaxed radius, which we obtain as follows. We measure the relaxation ratio and mass ratio of shells having fixed r_f/R_vir for all haloes in a selected sample, then stack them in bins of mass ratio at each spherical shell separately. We find that equation (7) is consistent with the relaxation relation of all halo masses considered, where we infer the values of q₀ and q₁ at each r_f using standard least-squares fitting (the reduced χ² values are always close to unity, for all masses and radial shells).

In Fig. 6, we show the measured relaxation relation for six different sample of haloes at shells of selected radii, compared with the best-fitting model equation (7) for each case; the model clearly describes these measurements extremely well. We already noted from Fig. 5, that the haloes in the small volume EAGLE simulation with the reference model shows a different relaxation behaviour. This is also apparent in Fig. 6, between the 10¹¹h⁻¹M_⊙ haloes from that L25 simulation (last row, left-hand panel) and the haloes of the same mass from IllustrisTNG (first row, left-hand panel); however they both follow the linear relaxation relation at fixed radii. An interesting feature to note is the dramatic change in slope q₁ for the most massive haloes (middle row, right-hand panel), with q₁ changing sign as one moves outwards through the halo. This also qualitatively explains the non-monotonicity and multivalued nature of the default stacks in the lower right-hand panel of Fig. 3.

Figure 6.

Relaxation relation stacked separately at five different radii indicated by colour in six samples of haloes selected by mass from IllustrisTNG and EAGLE simulations. We also show linear polynomial fit to this relation, following equation (7) with the best-fitting values for the parameters q₀ and q₁ at each of the selected radii for each of the six halo samples.

Open in new tabDownload slide

3.3.2 Modelling the radial dependence of the relaxation relation

As described above, we obtained the best fit parameters q₀ and q₁ of the relaxation relation at each relaxed radius r_f, this is shown in Fig. 7 as a function of r_f/R_vir.

Figure 7.

Linear quasi-adiabatic relaxation model parameters as a function of the radius of relaxed halo at different halo masses in IllustrisTNG (upper panel) and EAGLE (lower panel). The colour-coding follows Fig. 4. See the text for details.

Open in new tabDownload slide

For cluster-scale haloes, this simple monotonic dependence of q₁(r_f) is replaced with oscillatory behaviour. In fact, some of the peaks in q₁ correspond to q₁ ≈ 1; combined with |q₀| ≪ 1, this indicates that these peaks are shells which nearly perfectly conserve angular momentum. (e.g. this happens at |$r_f/R_{\rm vir}\simeq 0.5\, (0.2)$| for |$M=10^{13.5}\, (10^{14})\,{h^{-1}\,\mathrm{ M}_{\odot }}$|⁠.) This is in strong contrast to Fig. 5 where the slope of the relaxation relation represented by a globally defined q₁ (e.g. equation 7 without explicit r_f dependence in the parameters) is close to zero for these haloes, indicating maximum deviation from the adiabatic relaxation. This complicated behaviour in the cluster-scale haloes could be due to the presence of substructures. The q₀ parameter now also shows interesting behaviour; it is only slightly negative in the outer halo but becomes close to zero in the inner halo for these haloes (q₀ was relatively constant with more negative value for less massive haloes).

Figure 8.

Fitting values for the three parameters namely q₀, q₁₀, and q₁₁ in radially dependent quasi-adiabatic relaxation model described by equation 12 as a function of halo mass in the three IllustrisTNG simulations and the two EAGLE simulations as listed in Section 2.

Open in new tabDownload slide

3.4 Effect of properties beyond halo mass

In this section, we study our parametrized model of the response of dark matter in the halo as a function of halo and galaxy properties beyond halo mass. This is motivated by the fact that there is a significant scatter in the relaxation relation (see Fig. 3), even within halo samples selected by mass. Such a study would also have implications for halo assembly bias and similar environmental correlations predicted in the ΛCDM framework (see e.g. the discussion in Paranjape & Sheth 2021a). We focus on the four halo properties defined in Section 2.1.3; while these halo properties show overall trend with halo mass, they take a wide range of values even at fixed mass scale (see Fig. 9).

Figure 9.

The 10th, 50th, and 90th percentile of the four different halo properties in each of the sample selected by mass from three different cosmological boxes of the IllustrisTNG. This includes the concentration (c) of the unrelaxed halo, and the stellar fraction (f*), specific star formation rate (SSFR), and gas fraction (f_g) of the hydrodynamical halo, all of which are defined in Section 2.1.3.

Open in new tabDownload slide

For all subsamples selected by a secondary halo/galaxy property at fixed halo mass, we use the 3-parameter model discussed above in the mass range 10¹⁰ to 10^12.5 h⁻¹M_⊙, while for cluster-scale haloes, where our 3-parameter model fails, we directly compare the linear quasi-adiabatic relaxation model parameters q₀ and q₁ as a function of the scaled halo-centric radius r_f/R_vir. The results for low-mass (massive) haloes are shown in Fig. 10 (Fig. D1).

Figure 10.

Radially dependent quasi-adiabatic relaxation model parameters q₀, q₁₀, and q₁₁ estimated as a function of halo properties in IllustrisTNG simulations. In the top row panels, only halo mass dependence is shown; whereas in the next three rows, we show the dependence on halo concentration, stellar mass fraction, specific star formation rate, and gas fraction in terms of percentiles, respectively.

Open in new tabDownload slide

For reference, the upper panels of Fig. 10 show the best-fitting values of q₀, q₁₀, and q₁₁ as a function of halo mass alone for haloes with M ≤ 10¹³h⁻¹M_⊙, which repeat the corresponding curves in Fig. 8. The upper panels of Fig. D1 similarly repeat the results for q₀(r_f) and q₁(r_f) for the mass bins M = 10¹³, 10^13.5, 10¹⁴ h⁻¹M_⊙ from Fig. 7. By displaying both, the full radial dependence as well as the 3-parameter description for the mass bin 10¹³ h⁻¹ M_⊙, we can assess the reliability of the latter around the mass scale where it begins to fail. We repeat this for subsamples split by secondary halo/galaxy properties below.

3.4.1 Dependence on unrelaxed halo concentration

Then, we select haloes from each of the mass bins in three separated c_s percentile bins |$(10\pm 10, ~50\pm 10 ~\&~ 90\pm 10)$| and compute the relaxation relation as described in previous sections for each of those samples. We find that q₀ shows strong dependence on the concentration, with more concentrated haloes having higher value of q₀ at most halo masses with M ≤ 10¹³h⁻¹M_⊙ (see the second row in Fig. 10 and the second row, first column in Fig. D1). This can be understood in terms of the formation time of the halo, since concentration is correlated with the formation time. More concentrated haloes, that have formed earlier might have had enough time for the dark matter to respond to the baryonic feedback, and hence, there is less offset. On the other hand, q₁₀ or q₁₁ displays a more complex dependence at low mass, with no clear monotonic trend. Meanwhile, for cluster-scale haloes, q₁ shows a very different behaviour as a function of r_f at different halo concentrations (see the second row of Fig. D1). We leave a fuller exploration of these trends, particularly their dependence on substructure properties, to future work.

3.4.2 Dependence on baryonic halo properties

We now shift our focus to hydrodynamical halo properties. In this regard, we study the response as a function of the three halo properties defined in Section 2.1.3; namely the specific star formation rate (SSFR) at current redshift z = 0 and the total stellar mass fraction (f_*), and gas fraction (f_g) at redshift z = 0 which respectively represent the integrated star formation activity and gas content of the central subhalo. At each halo mass bin, we take three subsamples selected by bins of percentiles in f_*, SSFR, and f_g, in a similar fashion as with concentration significance.

From the third column of Fig. 10, we note that f_* does not affect the relaxation response significantly; in particular, the q₀ parameter is relatively least dependent on f_*, compared to other halo properties. This is consistent with the fact that the q₀ parameter clearly converges between three TNG boxes (see the upper panel of Fig. 10), despite large differences in f_* with resolution (see the upper right-hand panel of Fig. 9). On the other hand, we can see a clear trend in q₀ parameter with SSFR (see the first column in the third row of Fig. 10); at a given halo mass, the q₀ value is closer to zero when the star formation activity is lower. To recall, q₀ ≃ 0 would mean no offset in the relaxation relation, and in that case for shells having no relaxation, the mass ratio is unity indicating that the enclosed baryonic mass also remains the same. This result is consistent with our argument in Section 3.3 that the q₀ < 0 is caused by the recent baryonic outflows due to feedback which is lower in these low mass haloes when SSFR is low. From the first panel in the last row of Fig. 10, we can see a similar trend in q₀ with gas fraction f_g of the halo; this is likely due to the fact that the FOF haloes with more gas have relatively higher active star formation with larger recent baryonic outflows. However, even the haloes with similar SSFR and different gas fraction may show different relaxation behaviour. In a future work, we will study such effects using hydrodynamic simulations with different baryonic prescriptions, that produces haloes with same f_g but very different SSFR and vice versa. On the other hand, for the cluster-scale haloes shown in Fig. D1, the q₀ parameter does not vary significantly with any of the halo property that we considered.

Turning to q₁, while we see no clear dependence of its constituent parameters q₁₀ and q₁₁ on the hydrodynamical halo properties of low-mass haloes in Fig. 10, for cluster-scale haloes, we do see a strong albeit complex dependence of q₁ on SSFR (see the third row of Fig. D1). Like the case of the concentration significance, it will be interesting in future work to understand the physical mechanisms driving some of the stronger correlations of the halo response with properties such as SSFR and f_g seen above.

4 APPLICATIONS

In this section, we briefly discuss a few (potential) applications of our analysis.

4.1 Baryonification schemes

The results above show that the response of a halo’s dark matter content to the galaxy and gas evolving in it depends not only on the integrated properties of the halo and galaxy (such as mass, concentration, etc.) but also on halo-centric distance, even at fixed mass ratio. This is in stark contrast to analytical approximations employed in the literature which typically use simplified relations between the relaxation ratio and mass ratio, ignoring the radial dependence. These analytical approximations are now commonly employed in baryonification schemes to predict the total matter power spectrum for a given cosmological model using only the results of gravity-only N-body simulations (Schneider & Teyssier 2015; Chisari et al. 2018; Aricò et al. 2021). Our results above can directly impact such predictions by modifying the small-scale (deep 1-halo regime) behaviour of the power spectrum.

For example, to model the effect of baryons in low- and intermediate-mass haloes (≲ 10¹³h⁻¹M_⊙), we advocate the use of our fitting function equation (12) for the relaxation relation, with parameters set to q₀ ≃ −0.05, q₁₀ ≃ 1.1, and q₁₁ ≃ 0.5,⁴ which gives a good description of the results of both IllustrisTNG and EAGLE haloes (see Fig. 8). For larger (cluster-sized) haloes, the response is still accurately described by the relation equation (7), but with more complex behaviour for the parameters q₁(r_f) and q₀(r_f), which presently needs to be accounted for numerically (see e.g. Fig. D1). We discuss this further in Section 5.

4.2 Mass profiles

The primary utility of an analytical model (or fitting function) such as equation (7) for the relaxation relation is to be able to predict the relaxed dark matter profile |$M_f^d(r_f)$| of a halo which has responded to its baryonic content. The procedure for obtaining this profile is straightforward (see e.g. Appendix A of Paranjape et al. 2021b): The relaxation relation is solved iteratively using the unrelaxed mass profile and the baryonic mass profile as inputs, until a converged answer for |$M_f^d(r_f)$| is achieved.⁵ In our case, the procedure to obtain the relaxed dark matter profile using equation (7) works identically. The additional radial dependence of q₁ and q₀ is not an issue, since the radius r_f itself is used as the control variable in solving for the enclosed mass.

As an example, we compare the relaxed profiles predicted by this procedure – using the unrelaxed and baryonic mass profiles and the fits to the relaxation relation equation (7) as inputs – with the dark matter profiles actually measured in the hydrodynamical simulations for the same haloes. For simplicity, in this analysis, we ignore the dependence of the dark matter response on halo properties other than the mass; we use q₀ and q₁ as a function of r_f/R_vir as shown in the Fig. 7 for each halo mass bin. In the upper panel of Fig. 11, we show this estimated mass profile along with the actual mass profiles found in the IllustrisTNG simulation. For comparison, we also show the results of replacing the relaxation relation with simpler approximations from the literature (while still using the unrelaxed and baryonic mass profiles from the simulation as inputs in the iterative procedure). Our model produces significantly better estimates of the relaxed dark matter profile, especially in the inner halo, we obtain an order of magnitude better accuracy in comparison to such simple models (see lower panel of Fig. 11). In Appendix  C, we show that even the simple three parameter model gives a reasonably good prediction of the relaxed mass profile, while also being easier to incorporate into the existing procedures that use an adiabatic relaxation model.

Figure 11.

(Top row:) For the haloes in IllustrisTNG simulations, the mean radial mass profiles are shown in bins of halo mass for the baryonic component (dash–dotted curves) and dark matter component in hydrodynamic (dashed curves) and gravity-only (dotted curves) runs. The relaxed dark matter mass profile predicted as described in Section 4.2 is shown by solid curves. The colour-coding follows Fig. 4; for clarity, we use two panels to show the averaged mass profiles for the nine mass bins. (Bottom row:) The ratio of the relaxed dark matter mass profile predicted by our model to that from the hydrodynamic simulation is shown by solid curves. For comparison, the corresponding ratio for quasi-adiabatic relaxation model with q = 0.33 is shown by dashed curves and the ratio of dark matter mass profile between gravity-only simulation to the full hydrodynamic simulation is shown by dotted curves, representing the case of no relaxation.

Open in new tabDownload slide

4.3 Rotation curves

Since our model can predict relaxed mass profiles using unrelaxed and baryonic mass profiles as inputs, it can also predict rotation curves of galaxies using the same inputs, along with some assumptions regarding the geometry of various mass components. The interpretation of observed rotation curves and related statistics such as the radial acceleration relation, using data from spatially resolved spectroscopy of nearby galaxies, forms a key aspect of discussions in the literature regarding the nature of gravity at galactic scales (e.g. Lelli, McGaugh & Schombert 2016; Lelli et al. 2017).

In the ΛCDM context, such studies typically model the relaxed dark matter profile using a generalized NFW profile, with or without a core, but unconnected to the baryonic mass (e.g. Li et al. 2020). Previous work has suggested that the use of a parametrized model of dark matter response, rather than the relaxed profile itself, should lead to more robust results (Paranjape & Sheth 2021b; Paranjape et al. 2021a). For example, it is known that the use of smooth NFW-like profiles does not produce formally good fits in cases where the observed rotation curve shows oscillatory behaviour. Rather, these oscillations in the rotation curve correlate with similar oscillations seen in the measured baryonic mass profiles (see e.g. figs. 4 and 6 of Li et al. 2020).⁶ It is then reasonable to speculate that a model which smoothly parametrizes the physics of the dark matter response, rather than the profile of dark matter itself, might account for such correlations naturally. More generally, such a model is more physically motivated than one which directly parametrizes the dark matter profile itself.

In future work, we plan to confront observed rotation curves for low-mass systems with the 3-parameter relaxation model presented above. Our specific results for the values of these parameters in IllustrisTNG and EAGLE can then provide useful priors for the statistical comparison with data.

5 CONCLUSION

In this work, we have explored in detail the response of the dark matter content of a halo to the galaxy and gas it hosts. Understanding and accurately modelling this response is important for a number of applications including baryonification schemes for small-scale power spectrum emulation, rotation curve modelling, and constraining the nature of dark matter using inner halo mass profiles, etc.

Using haloes and galaxies identified in the IllustrisTNG and EAGLE simulations and matched to their gravity-only counterparts, our analysis demonstrates that the simplified analytical schemes used thus far to model the dark matter response (e.g. Blumenthal et al. 1986; Abadi et al. 2010; Schneider & Teyssier 2015) are inadequate in describing its detailed behaviour across a variety of halo and galaxy types. Specifically, we showed that the dark matter response, or relaxation relation (see equation 4), which connects the relaxation ratio r_f/r_i to the mass ratio M_i/M_f between unrelaxed (gravity-only) and relaxed (hydrodynamical) haloes, explicitly depends on halo-centric distance r_f in the relaxed halo, in addition to being sensitive to a number of halo and galaxy properties including halo mass, halo concentration, stellar and gas mass fraction, and specific star formation rate. These effects, especially the dependence on halo-centric distance, have been typically neglected by existing quasi-adiabatic relaxation models.

We presented a simple physically motivated extension (equation 7) of the existing models which accurately captures the dark matter response over 4 orders of magnitude in halo mass (10¹⁰ ≲ M/(h⁻¹M_⊙) ≲ 10¹⁴) and ∼2 orders of magnitude in relative halo-centric distance (0.02 ≲ r_f/R_vir ≤ 1). Apart from an explicit radial dependence of the relaxation relation (e.g. equation 12 for low-mass haloes), a second novelty of our model is the inclusion of a parameter q₀ which characterizes feedback-induced offsets seen in the relaxation relation measured in IllustrisTNG and EAGLE haloes in which, e.g. shells that do not show an overall change in radius (r_f/r_i ≃ 1) nevertheless have M_i/M_f > 1 (indicating loss of baryonic material). The existing quasi-adiabatic relaxation models do not allow for the existence of such shells, which are however captured well by our new null-offset parameter q₀ (see Section 3.3.1 for a detailed discussion). We argued that our results could have a significant impact on the applications listed above.

Our analysis also raises some interesting questions, which we briefly discuss before concluding. We noted in Section 3.3 that, unlike low-mass haloes whose relaxation relation is well described by equation (12), the radial dependence of the relaxation parameters q₁ and q₀ in equation (7) for haloes with M ≳ 10¹³ h⁻¹M_⊙ shows non-trivial features and oscillations that are not easily captured by simple fitting functions (see Fig. D1). These features, which typically occur in the halo outskirts, are likely due to the presence of substructure or recent mergers, which would generically lead to a disturbed dynamical state of the halo. In this work, we have not attempted to model these features; it will be interesting in the future to systematically study the dependence of these features on substructure fraction, merger history, locations of shocks, etc.

At the other extreme, in the inner halo of low-mass systems, it is very interesting to ask whether the simple quasi-adiabatic relaxation prescription we have calibrated in this work can naturally lead to cored inner dark matter profiles. Previous attempts at coupling the relaxed dark matter profile to baryonic physics using simple prescriptions have focused on introducing a baryonic dependence of the parameters describing the dark matter profile itself (e.g. Di Cintio et al. 2014). Our approach, on the other hand, parametrizes the physics of the dark matter response, and it will be interesting to see whether this leads to more robust results for cored inner profiles. For such an exercise, it will also be important to understand the dependence of our calibrated parameters on technical choices defining the sub-grid physics models used in the simulations, which can significantly impact the formation of cores (Benítez-Llambay et al. 2019).

Finally, building a more in-depth understanding of our results will need a physical, preferably analytical, model. One possibility is to use the self-similar approximation (Fillmore & Goldreich 1984; Bertschinger 1985; Lau et al. 2015; Shi 2016) to model the combined dynamical evolution of dark matter, gas, and stars in a halo. We will report the results of such studies in future work.

ACKNOWLEDGEMENTS

We thank Nishant Singh, Nishikanta Khandai and Kandaswamy Subramanian for useful discussions in the early phases of this work. We thank the anonymous referee for useful comments that improved the clarity of the presentation. We gratefully acknowledge the use of high performance computing facilities at IUCAA (http://hpc.iucaa.in). This work made extensive use of the open source computing packages numpy (Van Der Walt, Colbert & Varoquaux 2011),⁷scipy (Virtanen et al. 2020),⁸matplotlib (Hunter 2007),⁹pandas (The pandas development team 2022),¹⁰schwimmbad (Price-Whelan & Foreman-Mackey 2017),¹¹h5py,¹² Colossus (Diemer 2018b),¹³jupyter notebook¹⁴, and Code-OSS.¹⁵

DATA AVAILABILITY

The IllustrisTNG simulations are publicly available at https://www.tng-project.org/. The EAGLE simulations are publicly available at https://icc.dur.ac.uk/Eagle/.

Footnotes

REFERENCES

APPENDIX A: CHOICE OF MATCHING ALGORITHM

In this work, we studied the response of dark matter halo to galaxy by comparing the haloes in hydrodynamical simulations to their counterparts in gravity-only runs. We described the matching procedure in Section 2.2, here, we discuss some of the specific choices in that procedure. For this purpose, let us consider the FOF group haloes in TNG300 simulation with log M(h⁻¹ M_⊙) > 10.5. There are 543588 FOF groups satisfying this criterion in the hydrodynamical run of TNG300, with each of them having more than 500 particles within their R_vir in the highest resolution run. Following the matching procedure described in Section 2.3 (requiring only that the matching fraction of the hydrodynamical halo with respect to the gravity-only halo is greater than 0.5), we get a matching halo in the gravity-only run for 541 594 of them, leaving out only 1994 haloes unmatched, that is a negligible 0.4 per cent spread across the mass range (Fig. A1). However, if we follow the same procedure using matching fraction between the central subhaloes instead of FOF group themselves, then we get an order of magnitude more unmatched haloes as can be seen in the Fig. A1.

Figure A1.

Fraction of haloes in the TNG300-1 hydrodynamical simulation that have not found a match is shown as a function of mass M_200c.

Open in new tabDownload slide

While obtaining matches for as many haloes as possible, we also have to ensure the quality of match. To check how well the haloes in each of the pairs in our catalogue are matching, we look at the particle matching fraction of each halo in the pair with respect to the other as shown in Fig. A2. By definition the gravity-only halo in every pair has at least 50 per cent of the dark matter particles of the hydrodynamical halo in that pair. But note that in FOF group matching based catalogue, for a significant number of pairs the hydrodynamical halo in the pair has less than 10 per cent of the gravity-only halo’s particles. By visually inspecting some of those halo pairs, we find that, they represent haloes with ongoing merger events. This explains the significant loss in matching efficiency when we used central subhalo matching.

Figure A2.

Histogram of matching accuracy of haloes in the FOF group matched catalogue (left) and central subhalo matched catalogue (right). The x-axis is the fraction of dark matter particles in the hydrodynamical halo that is also in the gravity-only halo.

Open in new tabDownload slide

In the Fig. A3, matching loss as a function of halo mass is shown after removing all those pairs in which the particle matching fraction of the gravity-only halo with respect to hydrodynamical halo is less than 50 per cent. Since this symmetrical matching condition produces consistent matched catalogue of haloes, we apply this additional matching condition but stick to matching whole FOF groups. With this procedure, our final matched catalogue for TNG300 contains 524 841 pairs, leaving 18 747 out of 543 588 hydrodynamical haloes unmatched. The additional unmatched haloes primarily reside in dense environments, and we cannot find match for those haloes because of the inherent issues with 3D FOF algorithm in dealing with mergers.

Figure A3.

Fraction of haloes in the TNG300-1 hydrodynamical simulation that have not found a match is shown in left-hand panel as a function of mass M_200c after removing asymmetrical matches (see the text in appendix  A).

Open in new tabDownload slide

APPENDIX B: MOCK PARTICLES IN A GALAXY-HALO SYSTEM

Here, we validate our method used in extracting the relaxation relation (see Section 2.3.3) using mock galaxy-halo systems. In this, we use Hernquist profile for the initial unrelaxed dark matter halo and similar analytical mass profiles for the baryon components (see appendix of Paranjape & Sheth 2021b). We then get the relaxed dark matter profile using quasi-adiabatic relaxation model with different values of q. Mock particle data is then generated with these mass profiles for different components of hydrodynamical halo and the corresponding gravity-only halo.

For those mock halo pairs, we compute the mass profiles and obtain the relaxation ratio as a function of the dark matter shell radius as discussed in Section 2.3. We then repeat this for mock halo pairs generated with different particle resolutions. We find that the computed relaxation relation shown by coloured markers matches with expected relaxation relation shown in solid lines; however, the choice of radial bins is limited by the particle resolution (see Fig. B1).

Figure B1.

Mock halo with Hernquist profile for unrelaxed dark matter and with the analytical baryon mass profile from Paranjape & Sheth (2021b) for gas and star components. The relaxed dark matter profile is modelled by quasi-adiabatic relaxation with different choices of q parameter. Haloes are sampled by 10⁴ particles in left-hand panel and 10⁵ particles in the right-hand panel. (a) Mass enclosed within spherical shell as a function of the shell radius for different components. (b) Relaxation ratio as a function of the ratio of total mass enclosed. The linestyle indicates the different values for the model parameter q used in generating the mock data, while the colourbar shows the final radius. While the solid lines represent the expected relaxation relation for each case, the coloured markers show the computed relaxation relation from the mock particle data parametrized by the relaxed radius in units of the virial radius.

Open in new tabDownload slide

APPENDIX C: USING THE THREE PARAMETER MODEL

Here, we show the relaxed mass profile predicted by our 3-parameter model equation (12) for the haloes with mass, M ≤ 10¹³h⁻¹M_⊙. We follow a similar procedure as described in Section 4.2; once again, we consider the response as a function of only the halo mass and ignore the dependence on other halo properties. For a given halo, we obtain the values of the three parameters namely q₀, q₁₀, and q₁₁ by simply interpolating the fitted parameters shown in Fig. 8 as a function of mass. We find that accounting for radial dependence through a simple equation (12), gives a mass profile that is within |$10{{\ \rm per\ cent}}$| of the simulation even upto |$2 {{\ \rm per\ cent}}$| of virial radii (see upper left-hand panel of Fig. C1) and within |$3 {{\ \rm per\ cent}}$| for the low-mass haloes. The upper right-hand panel of Fig. C1, shows the corresponding profiles assuming a mass independent fixed values for the parameters q₀ ≃ −0.05, q₁₀ ≃ 1.1, and q₁₁ ≃ 0.5. In addition to these, we also compare with mass profiles predicted by the standard adiabatic relaxation Blumenthal et al. 1986 and few other relaxation models from Gnedin et al. (2004), Paranjape & Sheth (2021b), and Cautun et al. (2020) for the IllustrisTNG haloes.

Figure C1.

Ratio of the relaxed dark matter mass profile predicted by various models to the dark matter profile found in the hydrodynamical simulation IllustrisTNG is shown by solid curves. Here, the mass profiles are stacked across multiple haloes selected by their mass and the colour coding follows Fig. 4. Results from our three parameter model equation (12) is shown in the upper panel, while the left-hand panel accounts for mass dependence in the parameters, the right-hand panel assumes fixed values namely q₀ ≃ −0.05, q₁₀ ≃ 1.1, and q₁₁ ≃ 0.5. In the rest of the panels, the corresponding ratio is shown for the mass profiles predicted by few existing models as mentioned in the appendix  C. The result from this work, as shown in bottom panel of Fig. 11 is shown in dashed curves for reference.

Open in new tabDownload slide

APPENDIX D: DEPENDENCE ON GALAXY PROPERTIES IN THE CLUSTER SCALE

The effect of different halo and galaxy properties on the response of the dark matter is presented here; the relaxation parameters q₀(r) and q₁(r) are shown in Fig. D1 as our three-parameter description fails for these cluster scale haloes (see Section 3.4 for details).

Figure D1.

Similar to upper panel of Fig. 7 with cluster-scale haloes further split by other properties. In top row, only halo mass dependence is shown, whereas, in the next three rows, we further show the dependence on halo concentration, stellar mass fraction, specific star formation rate and gas fraction in terms of percentiles, respectively.

Open in new tabDownload slide