# Baryonic Imprints on DM Halos: the concentration–mass relation and its
dependence on halo and galaxy properties

Mufan (Jon) Shao1,⋆ Dhayaa Anbajagane2,3,⋆<EMAIL_ADDRESS><EMAIL_ADDRESS>1 Department of Physics, University of Chicago, Chicago,
IL 60637, USA  2 Department of Astronomy and Astrophysics, University of
Chicago, Chicago, IL 60637, USA  3 Kavli Institute for Cosmological Physics,
University of Chicago, Chicago, IL 60637, USA

###### Abstract

The halo concentration–mass relation has ubiquitous use in modeling the matter
field for cosmological and astrophysical analyses, and including the imprints
from galaxy formation physics is tantamount to its robust usage. Many
analyses, however, probe the matter around halos selected by a given
halo/galaxy property — rather than by halo mass — and the imprints under each
selection choice can be different. We employ the Camels simulation suite to
quantify the astrophysics and cosmology dependence of the concentration–mass
relation, $c_{\rm vir}-M_{\rm vir}$, when selected on five properties: (i)
velocity dispersion, (ii) formation time, (iii) halo spin, (iv) stellar mass,
and (v) gas mass. We construct simulation-informed nonlinear models for all
properties as a function of halo mass, redshift, and 6
cosmological/astrophysical parameters, with a mass range $M_{\rm
vir}\in[10^{11},10^{14.5}]\,{\rm M}_{\odot}/h$. There are many mass-dependent
imprints in all halo properties, with clear connections across different
properties and non-linear couplings between the parameters. Finally, we
extract the $c_{\rm vir}-M_{\rm vir}$ relation for subsamples of halos that
have scattered above/below the mean property–$M_{\rm vir}$ relation for a
chosen property. Selections on gas mass or stellar mass have a significant
impact on the astrophysics/cosmology dependence of $c_{\rm vir}$, while those
on any of the other three properties have a significant (mild) impact on the
cosmology (astrophysics) dependence. We show that ignoring such selection
effects can lead to errors of $\approx 25\%$ in the baryon corrections of
$c_{\rm vir}$. Our nonlinear model for all properties is made publicly
available.

##  1 Introduction

The halo concentration, which parameterizes the density profile of halos, is a
ubiquitous tool for modeling many astrophysical and cosmological observables.
Under the halo model prescription (Cooray & Sheth 2002), the
concentration–mass relation can be used, in conjunction with a model of the
halo abundance and halo bias, to model the matter density field on quasi-
linear and non-linear scales. This has been verified through extensive
validation of this model against predictions from N-body simulations which
solve the full non-linear dynamics of the density field. Gaining the ability
to model, and thus analyze, such scales allows us to extract significant
information on cosmological/astrophysical processes (e.g., Pandey et al. 2021;
Gatti et al. 2021; Zacharegkas et al. 2021; Miyatake et al. 2023).

While N-body simulations precisely solve the gravitational dynamics of the
matter density field, they do not model the presence of baryons and the
associated galaxy formation processes. Once such processes are included, the
density field has qualitatively different behaviors on these quasi-linear and
non-linear scales (e.g., Chisari et al. 2018; Schneider et al. 2019; Giri &
Schneider 2021) which then impacts small-scale measurements of observables
connected to the density field (Gatti et al. 2020; Secco et al. 2022;
Anbajagane et al. 2023a). The evolution of gas and stars in the halo
significantly alter the halo density profile. This evolution, in turn, depends
on the rich physics of galaxy formation. The downstream effects of such
processes on the concentration–mass relation have been studied previously, and
found to have strong mass- and redshift-dependence (Duffy et al. 2010; Beltz-
Mohrmann & Berlind 2021; Anbajagane et al. 2022a; Shao et al. 2023). The work
of Shao et al. (2023) calibrated a part of these baryonic effects —
parameterized with four amplitudes corresponding to various processes
associated with supernovae (SNe) and supermassive black holes (SMBH) — using
the Camels simulation suite (Villaescusa-Navarro et al. 2021).

While these calibrated relations are already useful tools for theoretical
modelling — as alternative, baryonified relations to the fiducial N-body based
concentration–mass relations — a potentially missing piece in the robust
application of these relations to data is the impact of halo selection effects
on the model. While these selection effects are completely negligible for
modelling the density field auto-correlations111The density field is the sum
of matter contributions from all halos across the full range of masses, and
thus there is no sense of a “halo selection” when working with only the
density field., they become important when correlating the density field with
any halo or galaxy field. For example, the correlation between galaxy
positions and weak gravitational lensing (commonly denoted as galaxy-galaxy
lensing or “gg-lensing” is a powerful tool in extracting cosmology from
optical surveys). The galaxies observed in these surveys are not a fair sample
of all possible galaxies in the Universe, but some specific subsample
determined by survey configuration (e.g., the wavelengths observed by the
survey, the observing depth etc.), and requirements on the redshift
estimation, survey completeness, and other systematic effects (Porredon et al.
2021; Rodríguez-Monroy et al. 2022; Bautista et al. 2021). An even simpler
selection is a magnitude threshold (or a signal-to-noise threshold) — which is
always inherent in astronomical observations — as it is fundamentally
different from a halo mass threshold. Thus, it is apt to check whether these
selections result in a halo sample for which the baryon imprints look
different from those of the halo mass-selected sample; an example of a
concrete question is whether the concentration of stellar mass-selected halos
depends differently on galaxy formation processes when compared to the
concentration in a gas mass-selected sample.

This question can be directly probed by subsampling halo catalogs from
cosmological hydrodynamical simulations to match the selections of a given
survey, and then calibrate the baryon imprints on the concentration–mass
relation for this specific subsample. In practice, the robust implementation
of such selections is highly non-trivial and requires data products (such as
the observed magnitudes of galaxies in a given survey) that are not easily
available/created. Additionally, such an analysis would only be useful for the
specific survey being considered given it is tuned for that specific survey’s
selection properties. In this work, we instead take a simpler and more general
approach of splitting the halo sample according to various halo and
galaxy222In this work, we denote halo properties as those that summarize the
full matter distribution, i.e. including the dark matter, whereas galaxy
properties include only the baryonic components. properties and quantifying
the differences in baryon imprints of the concentration–mass relation caused
by this split. Such a choice minimizes the direct applicability of the
calibrated relations to survey data, but in favor of enabling explorations of
a broader range of selection effects and thus types of surveys (optical
selection, X-ray selection etc.). The properties we choose to split the halo
sample by are the: (i) total matter velocity dispersion, $\sigma_{\rm 3D}$,
which is connected to X-ray luminosity (ii) halo formation time, $a_{\rm
form}$, which changes the star formation rate, and thus the galaxy colors
(iii) the halo spin $\lambda_{\rm h}$, (iv) stellar mass, $M_{\rm star}$,
which is directly proportional to galaxy luminosity and (v) gas mass, $M_{\rm
gas}$, which is directly proportional to luminosity in X-ray and millimeter
wavelengths. These properties were chosen by looking at all available
properties in the Rockstar catalogs of the Camels simulations and picking
those connected to different halo/galaxy selections.

previously used the combination of 1000 Camels simulations to analyze the
$c_{\rm vir}-M_{\rm vir}$ relation of dark matter halos and the imprints of
baryonic processes on it. That work built non-parameteric models of the halo
concentration $c_{\rm vir}(M_{\rm vir},z,\Omega_{\rm m},\sigma_{8},A_{\rm
SN1},A_{\rm AGN1},A_{\rm SN2},A_{\rm AGN2})$ as a function of halo mass,
redshift, two cosmological parameters, two supernovae feedback parameters, and
two active galactic nuclei (AGN) feedback parameters. In this work we build
upon those results by: (i) constructing models for the five additional
properties described above — $\sigma_{\rm 3D}$, $\lambda_{\rm h}$, $a_{\rm
form}$, $M_{\rm star}$, and $M_{\rm gas}$ — using the same simulations and
methods as , (ii) comparing mass-dependent features across all properties (the
five above in conjunction with $c_{\rm vir}$, as analyzed in ) to infer the
origin of different baryon imprints, and; (iii) probing changes in the
baryonic imprints on $c_{\rm vir}$ due to subselecting halos by a given
property.

Our work is organized as follows: in §2 we describe the Camels simulation
data, including the properties we study, and the kernel-localized linear
regression (Kllr) method used to analyze them. In §3 we extract the full
astrophysics and cosmology dependence of the property–mass relations for the
five properties, and in §4 we probe the baryon imprints on halo concentration
for halos under different halo and galaxy selection effects, and quantify the
error from ignoring such selections. We conclude in Section 5.

##  2 Data and Methods

###  2.1 The Camels simulations

The Cosmology and Astrophysics with MachinE Learning Simulations (Camels,
Villaescusa-Navarro et al. 2021), is a simulation suite of $V=(25\,{\rm
Mpc}/h)^{3}$ boxes, each with varying cosmological and astrophysical
parameters. The cosmological parameters include the fractional matter energy
density in the present epoch ($\Omega_{\rm m}$), and the root-mean-square of
deviations in the linear density field smoothed on a scale of $8\,{\rm Mpc}/h$
($\sigma_{8}$). The astrophysical parameters are $A_{\rm SN1},A_{\rm SN2}$,
which control the supernovae feedback, and $A_{\rm AGN1},A_{\rm AGN2}$ which
control the active galactic nuclei (AGN) feedback.

The Camels simulations were run with $256^{3}$ dark matter particles of mass
$6.5\times 10^{7}(\Omega_{\rm m}-\Omega_{\rm b})/0.251\,{\rm M}_{\odot}/h$,
and then $256^{3}$ gas resolution elements with an initial mass of $1.27\times
10^{7}\,{\rm M}_{\odot}/h$. The suite is specialized for training machine
learning algorithms but the dataset can also be used for non-parameteric
studies of astrophysics and cosmology. The full suite contains different types
of runs, of which two are of interest to this work. The 1P runs are 61
datasets (10 simulations per each of the 6 parameters and 1 fiducial run) that
systematically vary a single parameter. All runs share the same initial
conditions seed. The Latin Hypercube (LH) simulations are the main driver of
this work. These are 1000 simulations that span a 6D parameter space of the
input parameters of the simulation. The points in parameter space are chosen
via an optimal scheme called a Latin Hypercube, which is a sampling technique
that has become a frequent tool for building emulators for cosmology (Heitmann
et al. 2009). The priors for the parameters are $0.1<\Omega_{\rm m}<0.5$,
$0.6<\sigma_{8}<1.0$, $0.25<A_{\rm SN1,AGN1}<4$, and $0.5<A_{\rm SN2,AGN2}<2$.

The Camels simulation suite contains multiple galaxy formation models. In this
work, we use the simulations from the galaxy formation model of IllustrisTNG,
which we will interchangeably refer to as TNG, where the latter is a state-of-
the-art high-resolution suite of simulations (Nelson et al. 2018; Pillepich et
al. 2018b; Springel et al. 2018; Marinacci et al. 2018; Naiman et al. 2018;
Nelson et al. 2019; Pillepich et al. 2019). The model choices made in the
fiducial TNG simulations — which are described in detail in Weinberger et al.
(2017); Pillepich et al. (2018a) — are adopted in the Camels simulations. The
resolution level of the Camels simulation roughly corresponds to the TNG100-2
run.

There are four astrophysical parameters varied in Camels. The total energy
injection rate (power) per unit star-formation, $e_{w}$, and the SN wind
speed, $v_{w}$, are controlled by two normalization factors, $A_{\rm SN1}$ and
$A_{\rm SN2}$, respectively. Thus $A_{\rm SN1}$ corresponds to the energy
output rate, and $A_{\rm SN2}$ corresponds to the wind speed. The mass loading
factor, $\eta$, which corresponds to the amount of mass transferred away
through the SNe winds/feedback, is set by the combination $v_{w}^{-2}e_{w}$
and so is proportional to the combination $\eta\propto A_{\rm SN1}A_{\rm
SN2}^{-2}$.

The supermassive black hole (SMBH) model has different modes of feedback — a
quasar/thermal mode, where the gas in the vicinity of the SMBH is thermally
heated, and a radio/kinetic mode, where the gas is given a momentum kick
instead. The kinetic mode is turned on during the phase of the SMBH where it
has lower accretion rates, and is the dominant mechanism for halos at the
Milky Way-scale and above (Weinberger et al. 2018, A22). Note that the kinetic
feedback is not continuously output into the halo. Rather, in the TNG model,
the SMBH “accumulates” energy and expels it once a chosen energy scale is
reached, which results in a “bursty” feedback process. The parameter $A_{\rm
AGN1}$ controls the rate at which energy is accumulated by the SMBH and
$A_{\rm AGN2}$ controls the minimum energy needed to initiate such a burst.
Both AGN parameters control the rate of the bursts, while the latter also
controls the speed of the kinetic feedback jets.

The uber halo sample used in this work — which is obtained by combining halos
from all 1000 LH simulations — spans a mass range of $10^{11}\,{\rm
M}_{\odot}/h<M_{\rm vir}<10^{14.5}\,{\rm M}_{\odot}/h$; the minimum mass is
set by requiring at least $\approx 1000$ particles in a halo, and the maximum
mass is set by the number of available halos. We set the maximum $M_{\rm vir}$
as the 10th most massive halo mass in the uber sample. The uber sample, as
well as the mass range we use, follows the choices from .

In this work, we focus on six halo properties available in the Rockstar halo
catalogs as provided by Camels. The halo virial radius, $R_{\rm vir}$, is
defined as the radius of the sphere containing an enclosed density that is
$\rho=\Delta_{\rm vir}(z)\rho_{c}(z)$, where $\rho_{c}$ is the critical
density at a given epoch, and $\Delta_{\rm vir}(z)$ is the density contrast
described in Bryan & Norman (1998). The virial mass, $M_{\rm vir}$, is defined
as the mass enclosed within this radius, $M_{\rm vir}=4/3\pi R_{\rm
vir}^{3}\times\Delta_{\rm vir}\rho_{c}(z)$ and is the mass definition used
throughout this work. The other halo properties we use in this work are

  * •

Concentration ($c_{\rm vir}$), defined as $c_{\rm vir}=R_{\rm vir}/r_{s}$,
where $r_{s}$ is the scale radius defined in the NFW profile (Navarro et al.
1997),

$\rho(r)=\frac{\rho_{s}}{\bigg{(}\frac{r}{r_{s}}\bigg{)}\bigg{(}1+\frac{r}{r_{s}}\bigg{)}^{2}}$
(1)

and is computed in Rockstar through fits to the total matter density profile
of the halo.

  * •

Velocity dispersion ($\sigma_{\rm 3D}$), the standard deviation of particle
velocities with the halo, with the different cartesian components added in
quadrature to obtain the isotropic quantity. The isotropic velocity dispersion
definition follows the definition of Yahil & Vidal (1977) and is larger than
that of Rockstar by a simple, constant factor of $\sqrt{3}$. We implement this
latter factor in our work. When presenting the $\sigma_{\rm 3D}-M_{\rm vir}$
relation, we normalize the dispersion by $V_{\rm vir}=\sqrt{GM_{\rm
vir}/R_{\rm vir}}$ for visualization purposes.

  * •

Spin ($\lambda_{\rm h}$), defined as the modulus of the halo angular momentum
vector, normalized by $\sqrt{2}M_{\rm vir}V_{\rm vir}R_{\rm vir}$ as defined
in Bullock et al. (2001, see their Equation 5)

  * •

Formation time ($a_{\rm form}$), the scale factor at which the halo first
achieved half of its present day mass, $M_{\rm vir}(a=a_{\rm form})=0.5M_{\rm
vir}(a=1)$.

  * •

Stellar mass ($M_{\rm star}$), the total mass of stars within $R_{\rm vir}$.

  * •

Gas mass ($M_{\rm star}$), the total mass of gas within $R_{\rm vir}$. The
quantity uses all gas and has no temperature cuts.

###  2.2 Kernel-localized linear regression

Many halo property scaling relations are adequately fit with a simple power
law, i.e. a log-linear relation, while some are fit with extended parametric
models such as a broken power-law model with two power-law indices instead of
one. These parameterizations in general do not capture all the non-linearity
in a scaling relation. The use of simple power-laws to represent a scaling
relation is motivated by the observed behaviors in gravity-only simulations,
where there are few unique scales that can break a self-similar scaling; a
self-similar scaling would manifest in the form of a constant, scale-
independent power law index. The inclusion of galaxy formation introduces many
unique mass/time-scales (as we will show further below) and thereby cause
scaling relations for most halo properties to deviate from the power law
behavior (A22). With the surplus of simulation data available to us, we can
now explore data-driven models that accurately capture non-linearities without
assuming fixed functional forms.

One approach to building model-independent representations of scaling
relations is machine learning (ML) techniques, such as random forest, neural
networks, and more (e.g., Machado Poletti Valle et al. 2020; Stiskalek et al.
2022). Such methods are also more difficult to interpret, and thereby more
difficult to extract physical explanations from. Though there are approaches
that can alleviate the interpretability issue (e.g., Ntampaka et al. 2019;
Machado Poletti Valle et al. 2020, ). An alternative then is to continue with
well-known, well-understood linear models but extend the models to capture
some or most of the non-linearity.

Kernel-localized linear regression (Kllr, Farahi, Anbajagane & Evrard 2022b)
is an extension that allows such flexibility. The input data (i.e. halos) are
weighted according to a kernel, which is chosen to be a Gaussian kernel in
$\log_{10}M_{\rm vir}$, and a linear regression is fit to the kernel-weighted
data. As the kernel is shifted in the log-mass direction, the method measures
the non-linear scaling relation using multiple, connected linear pieces. Kllr
has had extensive use in extracting nonlinearity in different halo
property–mass scaling relations (Farahi et al. 2018; Anbajagane et al. 2020,
2022b, , ). In this work, we continue on the efforts of  and use Kllr to study
a multivariate scaling relation that is still modelled as being locally linear
in halo mass.

The locally linear, multivariate scaling relation of concentration is written
in Kllr as,

$\displaystyle\log_{10}Y=$ $\displaystyle\,\,\pi(M_{\rm
vir},z)+\alpha_{M}(M_{\rm vir},z)\log_{10}M_{\rm vir}$
$\displaystyle+\sum_{X}\alpha_{X}(M_{\rm vir},z)\log_{10}X,$ (2)

where $Y\in\\{c_{\rm vir},\sigma_{\rm 3D},a_{\rm form},\lambda_{\rm h},M_{\rm
star},M_{\rm gas}\\}$, $\pi$ is the intercept, $\alpha_{M}$ is the slope with
halo mass, $\alpha_{X}$ is the slope with other parameters, which are
$X\in\\{\Omega_{\rm m},\sigma_{8},A_{\rm SN1},A_{\rm AGN1},A_{\rm SN2},A_{\rm
AGN2}\\}$. The slopes and intercepts are all functions of halo mass and
redshift. The property, $Y$, is thus modelled as a function of mass, redshift,
and 6 cosmological/astrophysical parameters. The dependence on the
cosmological/astrophysical parameters is also a function of mass. When
analyzing an individual, single simulation (e.g., the individual 1P
simulations), we drop all terms with $\alpha_{X}$ (since none of the
parameters $X$ vary in these scenarios) and only keep the slope with mass. In
Section 4 below, we will discuss extensions to the model described above,
where the slopes and intercepts are functions of other halo properties as
well.

Figure 1: The Kllr $\sigma_{\rm 3D}-M_{\rm vir}$ scaling relation. We
normalize $\sigma_{\rm 3D}$ by $V_{\rm vir}$, so the variation across mass can
be more easily seen. The colors show results for different simulation
parameters. The solid lines are the prediction from our model, made using the
1000 LH simulations. There is a rich variety of trends in how the halo
property responds to changes in both the astrophysical and cosmological
parameters. The dashed lines are the Kllr result from the 1P simulations, and
we show our model captures the trends in the 1P runs.

We compute a Kllr model for each of the 34 available snapshots in Camels,
using the uber sample at each redshift between $0<z<6$. In practice, this
means we compute Equation (2.2) separately for the halo sample from each
snapshot. Our results in this work will focus on the model for just $z=0$ but
models for the full range of redshifts are provided publicly. We compute the
scaling relations for every $0.1{\,\,\rm dex}$ in mass, with the minimum mass
set at $10^{11}\,{\rm M}_{\odot}/h$ and a maximum mass set by the 10th most
massive halo of the uber sample. The upper limit choice excludes regions where
the model solutions are imprecise given the low number of halos. The Kllr
technique uses a Gaussian kernel in $\log_{10}M_{\rm vir}$, with a width of
$0.3{\,\,\rm dex}$. Uncertainties on all the Kllr parameters of equation (2.2)
are obtained by estimating regression coefficients for 100 bootstrap
realizations of the data. Thus, we have 100 versions of equation (2.2) that
then quantify the uncertainty in the predictions.

Other techniques, such as Gaussian processes, can also be used to capture the
non-linearity of a relation in a non-parametric way with minima assumptions.
For example, Kllr and Gaussian processes have been previously used to extract
the same non-linear relations in halo quantities (e.g., Farahi et al. 2021,
2022a). In this work, we continue using Kllr given its prior, extensive use in
characterizing halo scaling relations.

##  3 Dependence on astrophysics & cosmology

We detail the baryon imprints on each halo property in the following sections.
The halo concentration is omitted, as the results have already been presented
in . All results in this section come from our Kllr model fits, and
uncertainties are 68% confidence intervals obtained from 100 bootstrap
resamplings.

###  3.1 Particle velocity dispersion, $\sigma_{\rm 3D}$

Figure 1 shows that the slope of $\sigma_{\rm 3D}$ with $M_{\rm vir}$ is a
strongly non-monotic function of mass, consistent with  and with Anbajagane et
al. (2022b) who showed this behavior for a variety of simulation models. The
velocity dispersion increases with $\Omega_{\rm m}$ and $\sigma_{8}$ across
all masses, and the influence of $\sigma_{8}$ (as determined by the amplitude
of the slope between $\sigma_{\rm 3D}-\sigma_{8}$, shown in Figure 6) is
greatest, compared to that of $\sigma_{\rm 3D}$ with the other five simulation
parameters, for most of the mass range before dropping by a factor of 3 at
$M_{\rm vir}\sim 10^{14}\,{\rm M}_{\odot}$ (see Figure 6).

Singh et al. (2020, henceforth ) used 15 simulations from the Magneticum suite
(Hirschmann et al. 2014) that span a variety of cosmological parameters and
found that $\sigma_{\rm 3D}$ is independent of $\sigma_{8}$; in particular,
they determine the mass-independent slope between $\sigma_{\rm 3D}$ and
$\sigma_{8}$ to be $\alpha_{\rm\sigma_{8}}<0.01$ at 95% confidence (see their
Table 5). Our results indicate the mass-dependent slope drops towards
$\alpha_{\rm\sigma_{8}}\approx 0.05$ at $M_{\rm vir}=10^{14}\,{\rm
M}_{\odot}/h$ and suggests the slope decreases further towards higher mass.
The mass sample in  extends to $M_{\rm vir}=10^{15.5}\,{\rm M}_{\odot}$ (about
$1{\,\,\rm dex}$ higher than our sample). Therefore, we expect their estimate
of the mass-independent slope to be lower due to the presence of these massive
halos in their sample.  also show that the the slope between $\sigma_{\rm 3D}$
and the baryon fraction $f_{b}=\Omega_{\rm b}/\Omega_{\rm m}$ is
$\alpha_{f_{\rm b}}=-0.05\pm 0.01$. The Camels suite varies $\Omega_{\rm m}$
while fixing $\Omega_{\rm b}=0.049$, which means changes in $\Omega_{\rm m}$
can be directly translated to those in $f_{\rm b}$. The slopes of $\sigma_{\rm
3D}$ with these two properties are inter-related through a simple negative
sign; $\alpha_{\rm\Omega_{\rm m}}=-\alpha_{\rm f_{\rm b}}$. We find that at
high mass, our slope estimates are higher than those of , with our estimates
nearing $-\alpha_{\rm\Omega_{m}}=-0.1$, but are consistent with their work
when considering the behavior of halos with lower masses.

Figure 2: Same as Figure 1 but for the $a_{\rm form}-M_{\rm vir}$ scaling
relation.

Turning to the SNe feedback parameters, $A_{\rm SN1}$ and $A_{\rm SN2}$, we
see $\sigma_{\rm 3D}$ scales inversely with their amplitude. This is expected
as SNe feedback causes the halo matter distribution to become more disperse
(Bryan et al. 2013, , ), which in turn will lower $\sigma_{\rm 3D}$ due to the
strong correlation between the $\sigma_{\rm 3D}$ and $c_{\rm vir}$ ().
Weakening SNe feedback causes the maximum and minimum value of $\sigma_{\rm
3D}/V_{\rm vir}$ to shift to higher masses, and vice-verse when strengthening
the feedback. Lau et al. (2010) compare the velocity dispersion (of dark
matter) in hydrodynamic simulations that are either non-radiative or contain
gas cooling plus star formation. They find that the presence of cooling and
star formation increase $\sigma_{\rm 3D}$ (see their Figure 1) and this is due
to the adiabatic contraction of the halo potential enabled by the radiative
processes (Blumenthal et al. 1986; Gnedin et al. 2004). Increasing the
amplitude of SNe feedback hinders this cooling process, and will reduce
$\sigma_{\rm 3D}$, as observed in Figure 1.

Finally, the AGN parameters have a weak impact on this relation, as expected
from their weak impact on the halo concentration (). The $A_{\rm AGN1}$
parameter has no significant impact across the whole mass range, while $A_{\rm
AGN2}$ influences the relation above Milky way (MW) masses. At these larger
mass-scales, the thermal feedback of the AGN becomes less dominant, and the
kinetic feedback takes over ( A22, see their Figure 7). $A_{\rm AGN2}$
controls the energy of the kinetic bursts, and so its strongest impact is
expected to be in this mass range. Unlike with the SNe parameters, this AGN
parameter increases $\sigma_{\rm 3D}$. However, this is expected as $A_{\rm
AGN2}$ is directly proportional to the energy (and thus momentum) input into
the gas. Stronger bursts result in more significant momentum increases, which
will then increase the velocity dispersion within the halo. See Weinberger et
al. (2017) for details on the AGN implementation in IllustrisTNG. Note that
these results do not imply that the AGN as a whole has only a mild impact on
the halo internal dynamics. Instead, it implies the specific AGN processes
being varied in Camels — as parameterized by $A_{\rm AGN1}$ and $A_{\rm AGN2}$
— have a mild impact.

###  3.2 Formation time, $a_{\rm form}$

Figure 2 shows the scaling relation for $a_{\rm form}$ with $M_{\rm vir}$,
which is mostly log-linear/monotonic with the exception of variations induced
by changing the SNe feedback amplitudes. The shape of the fiducial relation
(best seen in the top right panel) is log-linear with a deviation at MW masses
also seen in ( , see their Figure 4). Larger values of $a_{\rm form}$ indicate
halos formed later in time, and the quantity is limited to the range $0\leq
a_{\rm form}\leq 1$.

First, we discuss the cosmology dependence, where $a_{\rm form}$ decreases
with $\sigma_{8}$ and $\Omega_{\rm m}$. The slope between $a_{\rm form}$ and
the former is nearly a factor of $6$ larger than the slopes with all other
simulation parameters (see Figure 6). It is nearly constant below MW masses
and increases with mass for halos above MW masses. The slope with $\Omega_{\rm
m}$ increases over the whole mass-range, from $-0.1$ at small masses to $0$ at
$10^{14}\,{\rm M}_{\odot}/h$, and increasing to positive values at higher
masses. At the highest masses, the positive slope is found to be non-zero at
$2\sigma$. If this is a robust feature, it may arise from a nonlinear coupling
between $\Omega_{\rm m}$ and other simulation parameters; such nonlinear
coupling has been found in other analyses of Camels, though in the context of
the astrophysical parameters (Gebhardt et al. 2023) and we discuss other such
couplings below. In general, we expect $a_{\rm form}$ to decrease with
$\sigma_{8}$ and $\Omega_{\rm m}$ as any increases in these parameters’
amplitudes (which change the amount of matter and amount of clustering of this
matter) make it easier for structures to form, which then causes their
formation to happen earlier and in turn cause $a_{\rm form}$ to decrease (Gao
et al. 2005; Li et al. 2008).

Figure 3: Same as Figure 1 but for the $\lambda_{\rm h}-M_{\rm vir}$ scaling
relation.

The SNe feedback parameters have a rich variety of impacts on $a_{\rm form}$.
At lower mass, $M_{\rm vir}<10^{11.5}\,{\rm M}_{\odot}/h$, they either have no
influence or they reduce $a_{\rm form}$. This is expected as in small halos
with shallow gravitational potentials, SNe feedback can blow out the gas. The
halo will then require more time to grow and achieve its present-day mass. At
higher masses, $10^{11.5}\,{\rm M}_{\odot}/h<M_{\rm vir}<10^{13}\,{\rm
M}_{\odot}/h$ this behavior qualitatively changes as increasing SNe feedback
causes halos to form earlier. This is an example of the non-linear coupling
between the many astrophysical processes that occur within a halo, and
showcases how varying a single parameter can change the impact of multiple
processes in the halo rather than the impact of just the single process
associated with that parameter.

In this case, the non-linear coupling is between the SNe feedback and the AGN
feedback, as discussed in Gebhardt et al. (2023). For strong SNe feedback, the
gas in these more massive halos is dispersed into the further outskirts of the
halo. The SNe are generally not strong enough to eject the gas out of the
halo, so the total mass within the halo remains conserved, while the
distribution of matter within it changes (e.g., see results on $c_{\rm vir}$
from ). When AGN feedback is a prominent process in the halo, which begins at
$M_{\rm vir}\approx 10^{11.5}\,{\rm M}_{\odot}/h$ and is maximized at $M_{\rm
vir}\approx 10^{12.5}\,{\rm M}_{\odot}/h$ (see , their Figure 7), it can eject
gas far outside the halo boundary and thereby alter the halo mass
distribution. However, the impact of the AGN jets depends crucially on the
amount of gas present in the vicinity of the AGN. If the gas is more dispersed
before the AGN turns on, then the impact of the AGN on the halo is weaker.
Since SNe affect this dispersion of gas, they also control the impact of the
AGN feedback.

Thus, increased SNe feedback enhances the dispersion of the gas and reduces
the impact of ejection from AGN jets. Halos then retain more of their gas and
require less time to achieve their present-day mass. They can then form later
in cosmic history. This causal picture is corroborated by the $M_{\rm gas}$
behavior detailed in Figure 5 below. Note that the conventional behavior of
increased SNe feedback reducing $a_{\rm form}$ is recovered at galaxy
group/cluster mass-scales, $M_{\rm vir}>10^{13.5}\,{\rm M}_{\odot}/h$. Note
that $A_{\rm SN1}$ has a further cross-over point at $M_{\rm vir}\approx
10^{14}\,{\rm M}_{\odot}/h$, which is within a factor of 2 of the mass-scale
of a similar crossover point found in the $\Omega_{\rm m}$ dependence. This
suggests the $\Omega_{\rm m}$ and $A_{\rm SN1}$ may be linked via a non-linear
coupling, similar to that of the SNe and AGN feedback.

The AGN feedback has very little influence on the $a_{\rm form}-M_{\rm vir}$
relation. The only notable dependence is on $A_{\rm AGN2}$ at the high mass
end, measured to be non-zero at $3\sigma$, and this implies that increasing
$A_{\rm AGN2}$ leads to halos needing less time to achieve their present-day
mass. Note that $A_{\rm AGN2}$ does not affect the total energy injected into
the halo, and instead affects the energy in each burst. Thus, the parameter
impacts the “burstiness” or stochasticity of the AGN jet feedback. In Figure 5
below we show that increasing $A_{\rm AGN2}$ reduces the gas mass in the halo
at the present time. This reduction will then change the corresponding
estimate of $a_{\rm form}$ to be higher. As a pedagogical example, consider a
halo of $M_{\rm vir}=10^{14.5}\,{\rm M}_{\odot}/h$ at $z=0$ with a given mass
accretion history $M_{\rm vir}(a)$. The formation time is defined as the scale
factor where $M_{\rm vir}(a_{\rm form})=0.5M_{\rm vir}(a=1)$. Thus, a process
that is most impactful at $a\sim 1$ will affect the present-day mass and
thereby affect the $a_{\rm form}$ estimate. If the process in question is a
large-scale ejection event, such as one generated by AGN feedback, then
$M_{\rm vir}(a=1)$ is reduced while the overall history $M_{\rm vir}(a)$ is
left unaffected, especially for the earlier epochs. In this case, $a_{\rm
form}$ is reduced as the mass to be achieved ($M_{\rm vir}$ at $a=1$) has been
reduced, and so the halo can achieve it at earlier times.

###  3.3 Spin, $\lambda_{\rm h}$

Figure 3 shows the halo spin–mass relation, which has a monotonic, but
nonlinear, behavior across the full mass range. Increases in $\sigma_{8}$
increase the $\lambda_{\rm h}$ for all masses considered, and the slope
$\alpha_{\rm\sigma_{8}}$ is the largest amongst the six simulation parameters
(see slopes in Figure 6). The influence of $\Omega_{\rm m}$ has a crossover
point at $M_{\rm vir}\approx 10^{12}\,{\rm M}_{\odot}/h$; the impact is
negative below this mass and positive above it. Given $\Omega_{\rm b}$ is
fixed in Camels, the variation in $\Omega_{\rm m}$ is a variation in $1/f_{\rm
b}$. Thus, variation in this parameter will also amplify/suppress any of the
baryonic effects in the simulation. This could explain the existence of a
crossover point, as well the mass scale of $10^{12}\,{\rm M}_{\odot}/h$, which
coincidences with the peak of the stellar mass fraction (see Figure 4) and
with the onset of AGN feedback in TNG (). Macciò et al. (2008) find that spin
in dark matter-only simulations has a weak dependence on cosmology and mass,
and Bryan et al. (2013) also find a weak dependence on cosmology in full
hydrodynamic simulations. Both tests of cosmology explore a narrower range of
cosmology parameters than the range probed in our work. The lack of any mass
dependence quoted in these two works is inconsistent with our findings, which
show a factor of two differences in $\lambda_{\rm h}$ across three decades in
$M_{\rm vir}$. Knebe & Power (2008), on the other hand, find the spin to
decrease with halo mass, consistent with our findings. The scatter in the
$\lambda_{\rm h}-M_{\rm vir}$ relation is $\approx 60\%$ which is similar in
magnitude to the change observed across the three decades in mass.

The SNe feedback parameters scale inversely with the spin at low masses: more
feedback results in less coherent rotation, and reduces the spin. As shown in
, we expect SNe feedback to be most impactful in halos below the MW mass. For
$A_{\rm SN1}$, the impact reverses sign at higher masses. This could be
related to the non-linear coupling of SNe and AGN discussed in Section 2,
where higher SNe feedback reduces the AGN feedback. Note also the similarity
in the variation of spin with $A_{\rm SN1}$ and with $\Omega_{\rm m}$ (better
thought of as $1/f_{\rm b}$ for this discussion): both have a crossover point,
though at slightly different masses. In $A_{\rm SN2}$, the inverse scaling —
where $\lambda_{\rm h}$ decreases with increases to $A_{\rm SN2}$ — is found
throughout the full mass range, with no crossover point. As we have noted
before, the impact of $A_{\rm SN2}$ on the halo property is through its
coupling to AGN feedback.

The AGN feedback parameters themselves have only a weak impact on the spin.
There is a slight variation at high masses, and we measure the slopes with
both parameters, $\alpha_{\rm AGN1/2}$, to be non-zero at $2\sigma$. At
$M_{\rm vir}\approx 10^{12}\,{\rm M}_{\odot}/h$, the AGN jet feedback begins
to play a role ( , see their Figure 7), and will inject kinetic energy into
the gas component.

Bryan et al. (2013, see their Figure 3) find that the inclusion of baryons
does not lead to a significant change in the $\lambda_{\rm h}-M_{\rm vir}$
relation. However, their sample has a median mass of $\approx 10^{12}\,{\rm
M}_{\odot}/h$, which does not probe the $M_{\rm vir}=10^{11}\,{\rm
M}_{\odot}/h$ halos whose spins have a stronger dependence on baryons as we
find here, and their halo sample is $\mathcal{O}(10^{3})$ times smaller than
the one we use here, meaning their uncertainties on the estimated
slopes/behaviors is larger than the trends we observe.

###  3.4 Stellar mass, $M_{\rm star}$

Figure 4: Same as Figure 1 but for the $M_{\rm star}-M_{\rm vir}$ scaling
relation. We normalize $M_{\rm star}$ by $M_{\rm vir}$, so the variation
across mass can be more easily seen. Thus, we plot the stellar mass fraction
rather than the stellar mass.

Moving to the baryonic quantities, Figure 4 shows the stellar mass
fraction–halo mass relation. The qualitative shape of the relation is similar
under variations of all six simulation parameters: there is a steep increase
at lower masses, a peak around MW mass scales, and a shallow dropoff towards
higher masses. The exact mass scales of each feature vary with changes to some
of the simulation parameters. The slope of $M_{\rm star}$ (i.e. the stellar
mass, not the stellar fraction) with $M_{\rm vir}$ approaches $2$ at the
lowest masses and then tends toward 1 from below starting around $M_{\rm
vir}=10^{13}\,{\rm M}_{\odot}/h$ (Figure 6). The latter, asymptotic behavior
is consistent with previous works that used the Kllr to find this trend in
multiple hydrodynamical simulations (Farahi et al. 2018; Anbajagane et al.
2020).

The impact of $\Omega_{\rm m}$ on the stellar fraction is approximately scale-
independent with a slope $\alpha_{\rm\Omega_{m}}\sim-1$, though it tends
towards $\alpha_{\rm\Omega_{m}}\sim-1.5$ at the lowest masses probed here (see
Figure 6). Recall that in Camels, $\alpha_{\rm\Omega_{m}}=-\alpha_{\rm f_{b}}$
as discussed above and so $\alpha_{\rm f_{b}}\sim 1$. Thus, the variation with
$\Omega_{\rm m}$ mentioned above can be removed through a simple rescaling of
$M_{\rm star}\rightarrow M_{\rm star}f_{\rm b}$, and this a commonly used
technique when comparing baryonic quantities across different simulations
(e.g., Lim et al. 2021; Lee et al. 2022). The dependence on $\sigma_{8}$ is
$\alpha_{\sigma_{8}}\approx 1$ for the full mass range.  finds that $M_{\rm
star}$ scales as $M_{\rm star}\propto f_{b}^{2.5}$ and $\sigma_{8}^{1.7}$, and
both scaling exponents are nearly twice as large as those found in our work at
similar masses ($\alpha_{f_{b}}\approx 1$ and $\alpha_{\sigma_{8}}\approx 1$).
The two estimates for $\alpha_{\sigma_{8}}$ are within $2\sigma$ of each other
and thus still in the regime of statistical consistency. However, as mentioned
before, the “cosmology dependence” studied in baryonic simulations includes
non-linear coupling between many galaxy formation processes. Thus, differences
in the simulated galaxy processes will cause differences in the manifested
cosmology dependence. Anbajagane et al. (2020, see their Figure 7) found that
the galaxy formation implementation in the Magneticum simulations, which are
used in , resulted in correlations between stellar properties that were
sometimes qualitatively different from those of other simulations, including
the IllustrisTNG model discussed in this work.

The SNe feedback parameters both scale inversely with $M_{\rm star}$ and
differ in the mass range where they are relevant. The scaling with $A_{\rm
SN1}$ follows the same trend as that with $\Omega_{\rm m}$, and is nearly
constant above MW masses with $\alpha_{\rm\Omega_{\rm m}}\approx-0.6$ and
drops to a slope of $-1.5$ below that mass-scale. The scaling with $A_{\rm
SN2}$ is also roughly constant above that mass-scale, and rapidly drops to a
slope, $\alpha_{A_{\rm SN2}}=0$, below it. Increases in the SNe feedback
amplitudes cause the stellar material to be more dispersed into the halo, and
lead to weaker subsequent star formation. This negative feedback loop of
higher $M_{\rm star}$ leading to more SNe events is detailed in Muratov et al.
(2015), where a suite of high-resolution, zoom-in simulations is used to show
that gas dynamics and star formation combine into a regulatory mechanism that
preserves the overall scaling of the $M_{\rm star}-M_{\rm vir}$ relation.

The slopes of AGN feedback parameters with $M_{\rm star}$ are essentially zero
for the full mass range. We reiterate that this only implies an insensitivity
of $M_{\rm star}$ to the specific processes being varied here; baryonic matter
is still significantly impacted by the presence of AGN feedback.

###  3.5 Gas mass, $M_{\rm gas}$

Figure 5: Same as Figure 1 but for the $M_{\rm gas}-M_{\rm vir}$ scaling
relation. We normalize $M_{\rm gas}$ by $M_{\rm vir}$, so the variation across
mass can be more easily seen. Thus, we are actually plotting the gas mass
fraction rather than the gas mass.

Finally, in Figure 5 we show the gas mass fraction–halo mass scaling relation.
The slope with $M_{\rm vir}$ changes significantly over the mass range but
with a median of around $\alpha=1$ (see Figure 6). At the high mass end, the
slopes is $\alpha\sim 1.5$ and tends towards $\alpha\sim 1$, consistent with
the results of Farahi et al. (2018). The gas fraction dips to a minimum at MW
masses.

Similar to the case with $M_{\rm star}$, we find the slope between $M_{\rm
gas}$ and $\Omega_{\rm m}$ to be $\alpha_{\rm\Omega_{m}}\approx-1$, which can
be removed through a simple rescaling of $M_{\rm gas}\rightarrow M_{\rm
gas}f_{\rm b}$. The gas mass also scales inversely with $\sigma_{8}$, and
shows a very strong mass trend, where the slope is $-0.5$ at both high and low
masses and then drops to $-1.2$ near MW masses. The gas dynamics in and around
the halo depend on the large-scale inflows of cold gas from the surrounding
environment (e.g., Baxter et al. 2021; Aung et al. 2020; Anbajagane et al.
2022c, 2023b), and at fixed halo mass the rate of inflow can be parameterized
by $\sigma_{8}$. Interestingly, increasing $\sigma_{8}$ decreases the gas mass
fraction whereas it increases the stellar mass fraction.  find $M_{\rm
gas}\propto f_{b}^{0.8}$, consistent with our results. They also indicate
$M_{\rm gas}$ is only weakly dependent on $\sigma_{8}$ when considering the
uncertainty on their measurement, with $M_{\rm gas}\propto\sigma_{8}^{-0.14}$.
However, their halo sample is limited to $M_{\rm vir}>2\times 10^{14}\,{\rm
M}_{\odot}$. Our result for this mass range is approximately $\alpha=-0.2$
(See Figure 6), which is consistent with their scaling to within $1\sigma$.

The SNe parameters have strongly non-linear relations with mass. Both
parameters decrease $M_{\rm gas}$ at low halo masses and increase it towards
high halo masses. The former behavior is related to the ability of SNe
explosions to eject gas out of the shallow halo potential of low mass objects.
The latter behavior is related to the non-linear coupling of SNe and AGN
discussed above for various properties, and also in Gebhardt et al. (2023).
Stronger AGN energetics naturally lead to a lower baryon fraction; see Lim et
al. (2021, see Figure 2) for a comparison of baryon fractions in various
simulations. SNe events also increase the available gas via the conversion of
stars into disrupted material sent back out into the halo.

Finally, the two AGN parameters continue to have the weakest impact of the six
simulation parameters. Though, their slope with $M_{\rm gas}$ is larger than
their slopes with any other halo properties, as expected from the tight
connection between AGN physics and gas dynamics. The impact of the $A_{\rm
AGN2}$ is strongest at $M_{\rm vir}=10^{13}\,{\rm M}_{\odot}/h$ while that of
$A_{\rm AGN1}$ is close to zero. Mitchell et al. (2020, see their Figure 1)
show that the mass outflow rate into the interstellar medium has a maximum at
$M_{\rm vir}\approx 10^{12}\,{\rm M}_{\odot}$, whereas the rate of outlow out
of the halo and into the large-scale structure is nearly mass-independent.
Their work also shows that above (below) the mass-scale quoted above, the mass
outflow is dominated by the AGN (SNe). This is consistent with the results and
the picture presented above, where the AGN impact on the gas increases above
MW masses.

Figure 6: The Kllr slopes of different halo and galaxy properties with halo
mass $M_{\rm vir}$ and with the 6 astrophysical and cosmological parameters in
the TNG model. Each panel corresponds to one halo property as indicated by the
label on top, each color represents the slope with respect to one parameter,
and the bands show the 68% confidence interval. The slopes show strong non-
monotonic features for every property.

Uniqueness of the MW mass-scale. For halos of masses $M_{\rm
vir}=10^{11.5}-10^{12.5}\,{\rm M}_{\odot}/h$ we find key features in all
properties discussed in the subsections above. At this mass scale, the SNe
dependence of $c_{\rm vir}$, $\sigma_{\rm 3D}$, $\lambda_{\rm h}$ and $M_{\rm
star}$ is maximized (see Figure 6). This mass coincides with crossover points
in the sign of the scaling of $\lambda_{\rm h}$ with $\Omega_{\rm m}$, and of
$M_{\rm gas}$ and $a_{\rm form}$ with SNe parameters. This mass is where the
stellar mass fraction peaks, which has been observationally confirmed through
many works (Behroozi et al. 2010; Moster et al. 2010; Yang et al. 2012; Wang
et al. 2013; Behroozi et al. 2013; Reddick et al. 2013; Moster et al. 2013;
Birrer et al. 2014; Lu et al. 2015; Rodríguez-Puebla et al. 2017; Shankar et
al. 2017; Kravtsov et al. 2018; Behroozi et al. 2019). In the IllustrisTNG
model, it is the mass scale where gas cooling and star formation are most
efficient and coincides with the onset of AGN feedback processes ().

##  4 Baryon imprints on $c_{\rm vir}$: Dependence on other halo/galaxy
properties

The halo properties discussed above all have significant correlations with one
another, and we see similarities in baryonic imprints across different
properties. Thus, the selection of halos on one of the five properties above
will affect not just the scaling relation of other halo properties, but also
the baryonic imprints on said properties. While the former is an effect in the
mean $X-M_{\rm vir}$ relation of the selected halo sample, the latter is an
effect in the change of the $X-M_{\rm vir}$ relation with changes in
parameters that control galaxy formation processes. In this case, we can
rewrite Equation (2.2) as

$\displaystyle\log_{10}Y=$ $\displaystyle\,\,\pi(M_{\rm
vir},W,z)+\alpha_{M}(M_{\rm vir},W,z)\log_{10}M_{\rm vir}$
$\displaystyle+\sum_{X}\alpha_{X}(M_{\rm vir},W,z)\log_{10}X,$ (3)

where the intercept $\pi$ and the slopes $\alpha$ are now all functions of the
halo properties $W\in\\{c_{\rm vir},\sigma_{\rm 3D},a_{\rm form},\lambda_{\rm
h},M_{\rm star},M_{\rm gas}\\}$ as well and where the set $W$ does not contain
the property $Y$. We have thus far assumed the slopes $\alpha_{X}$, which
define the dependence of property $Y$ on cosmology and galaxy formation
parameters, do not depend on the halo properties other than $M_{\rm vir}$.
However, these slopes are generically expected to depend on other properties.
For example, the baryon fraction determines the specifics of different
processes within the halo (e.g., Pandey et al. 2023; Grandis et al. 2023).

The Kllr model defined by Equation (4) is computationally expensive to build,
as the slopes must be computed for a grid of dimension N + 1, where N is the
number of halo properties it is a function of (the extra dimension is
redshift). In Equation (2.2) this was just a 2D grid spanned by redshift and
mass. In this work, we do not calibrate the full N + 1 dimensional grid and
instead take the simpler approach of computing the slopes in 2D grids but for
pairs of different halo samples,

$\displaystyle\log_{10}Y^{W\pm}=$ $\displaystyle\,\,\pi^{W\pm}(M_{\rm
vir},z)+\alpha_{M}^{W\pm}(M_{\rm vir},z)\log_{10}M_{\rm vir}$
$\displaystyle+\sum_{X}\alpha^{W\pm}_{X}(M_{\rm vir},z)\log_{10}X,$ (4)

where the superscript ${W\pm}$ denotes a model corresponding to (and built
with) halos where the value of a single property $W$ is higher or lower than
some chosen value.

In this section, we extract the relation $c_{\rm vir}(M_{\rm
vir},z,\Omega_{\rm m},\sigma_{8},A_{\rm SN1},A_{\rm AGN1},A_{\rm SN2},A_{\rm
AGN2})$ for halos selected on the five properties, split into selections on
halo properties (which include DM and baryons) in §4.1 and on galaxy
properties (which include only baryons) in §4.2. We compute $\alpha^{W\pm}$
and discuss the changes in the estimated slopes between pairs of halo
subsamples. The subsample pair is constructed by taking the uber sample of all
halos from the 1000 LH simulation. Then, for each halo and chosen selection
property, $W$, we predict the mean property value given the halo mass,
redshift, and six input simulation parameters. The models used to predict the
mean are the same ones shown above in Section 3. Each halo is then placed into
one of two subsamples (denoted by $W\pm$) depending on whether the value of
$W$ for that halo is above/below the predicted mean. We then rebuild the
$c_{\rm vir}-M_{\rm vir}$ model separately for each subset, which involves
estimating the slopes $\alpha^{W\pm}$ in Equation (4).

Note that such a split will cause a difference in the mean $c_{\rm vir}-M_{\rm
vir}$ relation, i.e. the intercept $\pi^{W+}(M_{\rm vir},z)$ will differ from
$\pi^{W-}(M_{\rm vir},z)$, in a way defined by the correlation between the
property $Y$ and the split variable $W$. The impact of this for correlated
variables is detailed in the formalism of Evrard et al. (2014). While we will
briefly discuss the impact of selection effects on $\pi^{W\pm}(M_{\rm vir})$
in Figure 8, our discussions focus on $\alpha^{W\pm}(M_{\rm vir})$ as this
quantifies changes in the baryonic imprints due to different halo/galaxy
selections.

###  4.1 Halo properties

Figure 7: The slopes of $c_{\rm vir}$ with the six parameters, for halo
samples split by a given property (colored lines), where the split is on halos
that lie above/below the mean relation of a given property. The fiducial
slope, estimated using the full, unsplit halo sample, is shown in the black
dotted line. The scaling with astrophysics parameters is mostly unchanged
through this selection, while that with cosmology parameters is modified
significantly.

Figure 7 presents the slopes of the halo concentration with respect to the
four astrophysical and two cosmological parameters. The slopes are shown for
halos split on different properties, as defined in Equation (4). The changes
in mass dependence are shown in Figure 8. In Figure 7, the solid (dashed)
curves correspond to the subset of halos with a specific property value
greater (lesser) than the mean. The average of the two slopes is not
guaranteed to match the slope of the fiducial sample. The mean relations can
follow such a behavior but the slopes, which are the derivative of the mean
relation with respect to mass and thus an estimate of the shape of the
relation, need not follow such behavior.

Focusing first on the cosmology dependence, selecting halos on $a_{\rm form}$
and $\lambda_{\rm h}$ has a significant impact on $\alpha_{\Omega_{\rm m}}$
while $a_{\rm form}$, $\lambda_{\rm h}$, and $\sigma_{\rm 3D}$ all affect
$\alpha_{\sigma_{8}}$. The difference in $\alpha_{\Omega_{\rm m}}$ is nearly
mass-independent. The results for $\alpha_{\sigma_{8}}$, however, show a
strong mass-dependence when split on $\lambda_{\rm h}$ and $a_{\rm form}$,
while the split on $\sigma_{\rm 3D}$ still leaves a nearly mass-independent
effect. All three of $a_{\rm form}$, $\lambda_{\rm h}$, and $\sigma_{\rm 3D}$
are highly correlated with the mass accretion rate (Lau et al. 2020;
Anbajagane et al. 2022a) and in turn on the local environment of the halo.
Halos in more dense regions can form earlier (higher $a_{\rm form}$) and have
more concentrated potentials (Wechsler et al. 2002), leading to higher
velocity dispersions. The spin depends on the level of coherent velocities in
the halo, and high mass accretion present in dense environments will suppress
such coherent rotational velocities (as the incoming matter comes from
multiple different directions). Thus, high $\lambda_{\rm h}$ objects will be
found in less dense environments. Thus, these selections are highly correlated
with environment-dependent selections. The evolution of halos in more dense
environments is more strongly affected by the density of matter in the
Universe, or by the amplitude of fluctuations in this matter field. Thus,
environment-driven physics can explain the impact of the three property
selections on the slopes $\alpha_{\sigma_{8}}$ and $\alpha_{\Omega_{\rm m}}$.

Moving onto the SNe feedback, we see that most selections cause very little
change as the slopes of the subsamples overlap with that of the fiducial
relation (black dotted line). All halo samples show the same phenomenological
features, such as local minima and maxima. The main differences occur from
selecting halos on $\sigma_{\rm 3D}$. We see that in high $\sigma_{\rm 3D}$
systems, the concentration is less dependent on $A_{\rm SN1}$ for
$10^{12}\,{\rm M}_{\odot}/h<M_{\rm vir}<10^{13}\,{\rm M}_{\odot}/h$, and in
low $\sigma_{\rm 3D}$ systems it is less dependent on $A_{\rm SN2}$ for
$M_{\rm vir}<10^{12.5}\,{\rm M}_{\odot}/h$. The SNe wind speed modeled in
IllustrisTNG relies on the local (dark matter) velocity dispersion around the
stars (Pillepich et al. 2018a, see their Equation 1). This implementation is
based on previous semi-analytic models that found such a scaling was needed to
reproduce the observed stellar mass functions and restframe B-band and K-band
luminosity functions across redshift (Henriques et al. 2013). Given this
choice, it is natural for the scaling with $A_{\rm SN2}$, where this parameter
controls the wind speed amplitude (Villaescusa-Navarro et al. 2021, see their
Equation 4), to depend $\sigma_{\rm 3D}$ selections.

Finally, focusing on AGN feedback, the halo selections cause no statistically
significant difference in the scaling of $c_{\rm vir}$ with the AGN-related
parameters. The strongest feature (which is still only a $2\sigma$ difference)
is that halos of high $\sigma_{\rm 3D}$ or low $\lambda_{\rm h}$ have $c_{\rm
vir}$ increase with $A_{\rm AGN}$.

In Figure 8, we also compare the mean $c_{\rm vir}-M_{\rm vir}$ scaling
relations of the halo samples split on $a_{\rm form}$, $\sigma_{\rm 3D}$ or
$\lambda_{\rm h}$. These splits result in clear differences in the $c_{\rm
vir}-M_{\rm vir}$ relations that are nearly constant across the whole mass
range. The shape of the relation is similar across all three subsamples,
meaning the slope/shape of $c_{\rm vir}-M_{\rm vir}$ does not change much with
selections on these three properties. Similar to our discussion above, we find
$c_{\rm vir}$ is positively correlated with $\sigma_{\rm 3D}$ and negatively
correlated with $a_{\rm form}$ and $\lambda_{\rm h}$, consistent with previous
work (e.g., Wechsler et al. 2002; Lau et al. 2020; Hearin et al. 2021;
Anbajagane et al. 2022a).

Figure 8: The $c_{\rm vir}-M_{\rm vir}$ relation for two halo subsamples,
defined as halos with a property value above/below the mean value at that
mass. All splits show clear differences in at least some part of the halo mass
range. The split on $M_{\rm gas}$ causes a crossover in $c_{\rm vir}$ towards
higher halo masses.

###  4.2 Galaxy properties

Finally, we consider the impact of baryonic property selections on the slopes
$\alpha_{X}$. Under this selection, all slopes experience significant changes
as seen by the differences between the dashed and solid lines in each plot; as
expected given the direct connection between galaxy formation processes and
baryonic properties.

Figure 9: Similar to Figure 7 but for selections/splits on $M_{\rm star}$ and
$M_{\rm gas}$.

There are significant differences in the scaling of $c_{\rm vir}$ with
cosmology parameters when comparing the two subsamples. Starting with
$\Omega_{\rm m}$, we find that the $M_{\rm gas}$ selection is most impactful,
particularly for $M_{\rm vir}<10^{13.5}\,{\rm M}_{\odot}/h$. The $M_{\rm
star}$ selection mimics the $M_{\rm gas}$ one in the range $10^{12}\,{\rm
M}_{\odot}/h<M_{\rm vir}<10^{13.5}\,{\rm M}_{\odot}/h$, suggesting that halo
baryon selection is the key factor here (since we find halos with high $M_{\rm
star}$ or high $M_{\rm gas}$ to cause similar changes in the slopes). The
change in $\alpha_{\sigma_{8}}$ is fairly similar across most of the mass
range. Below $M_{\rm vir}<10^{12.5}\,{\rm M}_{\odot}/h$, we see that halos
with high $M_{\rm star}$ are qualitatively similar to those with low $M_{\rm
gas}$, or vice-versa, as they both have a stronger dependence on $\sigma_{8}$.
This sign flip in the behavior is consistent with Figure 4 and 5, where we
show that increasing $\sigma_{8}$ leads to increases in $M_{\rm star}$ but
decreases in $M_{\rm gas}$.

Moving onto the changes in the scaling between $c_{\rm vir}$ and SNe
parameters, we find some variations due to selections on stellar/gas mass and
they are all statistically significant. For $A_{\rm SN1}$, though, the changes
are small in amplitude (the slope varies by $\Delta\alpha<0.05$). The $A_{\rm
SN2}$ dependence is much more significant, where halos with low stellar/gas
mass have a lower dependence on this parameter. Since SNe winds depend on the
amount of stars, and the impact of baryons on the concentration is through
redistribution of the gas, reducing $M_{\rm star}$ or $M_{\rm gas}$ reduces
the impact of this parameter on $c_{\rm vir}$.

The AGN parameters once again have the smallest variations. The gas selection
impacts the $A_{\rm AGN1}$ and $A_{\rm AGN2}$ dependence at low masses,
$M_{\rm vir}<10^{12}\,{\rm M}_{\odot}/h$. We do not discuss this further given
the low amplitudes of these parameters.

In Figure 8, we also compare the $c_{\rm vir}-M_{\rm vir}$ scaling relations
of the halo samples split on $M_{\rm star}$ and $M_{\rm gas}$. These splits
result in clear differences in the $c_{\rm vir}-M_{\rm vir}$ relation, and
differences that have unique a mass-dependence. Splitting on $M_{\rm star}$
induces a larger change in the $c_{\rm vir}$ normalization for $M_{\rm
vir}<10^{13}\,{\rm M}_{\odot}/h$, and almost no significant change above that
mass. Splitting on $M_{\rm gas}$ causes smaller changes at $M_{\rm
vir}<10^{12}\,{\rm M}_{\odot}/h$ and large changes above this mass scale. This
latter split also causes a crossover at $M_{\rm vir}\approx 10^{13.5}\,{\rm
M}_{\odot}/h$. Above (below) this mass, halos with more gas have a higher
(lower) concentration. In the TNG model, the internal energetics of halos
above (below) this mass-scale are dominated by gas cooling (AGN feedback); see
Figure 7 in . When selecting halos with more gas at fixed halo mass, we are
selecting halos whose gas fraction is higher, which means the gas can affect
the internal dynamics more. At high (low) masses the dominance of gas cooling
(AGN feedback) means more of the matter can be concentrated in the halo core
(dispersed out of the halo core), which leads to an increase (decrease) in
$c_{\rm vir}$. The same behavior does not occur for the $M_{\rm star}$ split
as stellar matter is both a much smaller fraction of the total matter, and
does not directly couple to feedback.

Estimate of potential modelling error. As expected, selections on gas mass or
stellar mass result in different slopes, $\alpha_{X}$, between $c_{\rm vir}$
and the galaxy formation parameters $X$. An accurate treatment of baryon
imprints on the density field around halos, where the halos are selected by an
observable connected to emission from the gas or stellar component, would
require a more general model than a mass-dependent one alone, i.e. the slopes,
$\alpha_{X}$, cannot be represented as just functions of halo mass. We make a
quantitative estimate of this by computing the differences between two
predictions for various parameters. This can be written as

$\displaystyle\Delta\log_{10}c_{\rm vir}(X,W)=$
$\displaystyle\,\,(\alpha_{M}^{W+}-\alpha_{M}^{W-})\log_{10}M_{\rm vir}$
$\displaystyle+\sum_{X}(\alpha_{X}^{W+}-\alpha_{X}^{W-})\log_{10}X$ (5)

where $X$ is a set of values for the six simulation parameters, and $W$ is the
property we perform the split with. The individual slopes, $\alpha^{W\pm}$,
are shown in Figure 7 and 9. By computing $\Delta\log_{10}c_{\rm vir}$ for
various choices of input $X$, we can quantify the possible modeling errors
caused by the selection, $W^{\pm}$. These can be thought of as errors arising
from an improperly specified model e.g., the theoretical error from a model
built on halos with high $W$ being applied to halos with low $W$. Equation
(4.2) only includes the slopes and not the normalization, $\pi^{W\pm}$, as
only the former captures the response of the $c_{\rm vir}-M_{\rm vir}$
relation to changes in galaxy formation and cosmology. We have also already
shown in Figure 8 the selection-based differences in the normalization, $\pi$.
Thus, the error defined in Equation (4.2) is for models that specify how the
$c_{\rm vir}-M_{\rm vir}$ relation responds to changes in the six simulation
parameters.

Figure 10 presents the modelling errors in the baryon imprints due to halo
selection effects. We show $\Delta\ln c_{\rm vir}=\ln(10)\Delta\log_{10}c_{\rm
vir}$, as it can be interpreted as a fractional difference in $c_{\rm vir}$.
These differences depend on the choice of input parameters and so we first
sample $10^{4}$ points in the 6D parameter space, which is bound by
$0.1<\Omega_{\rm m}<0.5$, $0.6<\sigma_{8}<1.0$, $0.25<A_{\rm SN1},A_{\rm
AGN1}<4$, and $0.5<A_{\rm SN2},A_{\rm AGN2}<2$. This bound is simply the range
of parameters sampled in the LH suite of the Camels simulations (see Section
2.1). The $10^{4}$ points are sampled from a uniform distribution in the
linear (not logarithmic) parameters. For each point, we compute $\Delta\ln
c_{\rm vir}$, the fractional modelling error in $c_{\rm vir}$, as denoted in
Equation 4.2. Figure 10 shows, as colored bands, the 68% bounds of this error
across the $10^{4}$ points in parameter space. The dashed lines are $\Delta\ln
c_{\rm vir}$ for the fiducial parameter values of $\Omega_{\rm m}=0.3$,
$\sigma_{8}=0.8$ and $A_{X}=1$.

The bands in the top panel show the range of errors in a 6D parameter space
comprising all six simulation parameters, the middle and bottom panels show
the same in 2D spaces varying either only cosmology parameters or only SNe
parameters, respectively. In all panels, the dashed lines show the error when
assuming the fiducial simulation parameters, and the bands show the range of
possible errors if we vary certain sets of parameters. The modeling errors are
found to be $\approx 25\%$ in amplitude under selections on $M_{\rm gas}$ or
$a_{\rm form}$. This is generally lower than, but of the same order as, the
scatter in $c_{\rm vir}$ (which is $\approx 40\%$) and it is also comparable
to differences between the $c_{\rm vir}-M_{\rm vir}$ relation in dark-matter
only simulations and full hydrodynamic simulations from IllustrisTNG (). Thus,
selection-based errors in the baryon imprints model are not guaranteed to be
subdominant, though their exact importance will depend on the dataset being
modeled. The above suggests that halo selection may be a relevant factor to
consider when building simulation-calibrated models for the impact of galaxy
formation on the density field.

Figure 10: The fractional modeling error, $\Delta\ln c_{\rm vir}$, due to
ignoring selection effects on a given halo property. If the value is positive
(negative), then the fiducial baryon imprints model, which accounts for no
selection effects, would underpredict (overpredict) the $c_{\rm vir}$ of the
property-selected halo sample. The bands show the 68% distribution of errors
across a wide prior in the simulation parameters. The dashed lines are the
error for the fiducial simulation parameter values. The top panel shows
results from varying all six parameters, middle panel from varying cosmology
alone, and bottom panel from varying SNe parameters alone. Selection on
$M_{\rm gas}$ causes differences of the order 25% at high and low masses.
Selection on $a_{\rm form}$ or $M_{\rm star}$ (which impacts galaxy colors and
luminosity) has an effect of a similar magnitude but at galaxy group scales.

##  5 Conclusion

The concentration–mass relation has ubiquitous use in modeling the matter
fields in our Universe on quasi-linear and non-linear scales. The impact of
baryons on the concentration, and the matter distribution, is an active area
of study. An equally pertinent question is the impact of baryons around
certain types of halos, such as those selected by a specific kind of
observational technique — X-ray, optical, Sunyaev Zeldovich decrement etc. In
this work, we provide a first attempt at quantifying the impact of halo/galaxy
selection on the scaling of $c_{\rm vir}$ with galaxy formation parameters. We
use 1000 simulations from the Camels suite to build a model of halo
concentration ($c_{\rm vir}$), velocity dispersion ($\sigma_{\rm 3D}$), spin
($\lambda_{\rm h}$), formation time ($a_{\rm form}$), stellar mass ($M_{\rm
star}$), and gas mass ($M_{\rm gas}$). We first analyze the dependence of each
individual property on galaxy formation parameters, and then also estimate the
change in the slopes $\alpha_{X}$ — which quantify how $c_{\rm vir}$ scales
with these parameters — after selecting halos on one property. Our results are
as follows:

  * •

The relationship between the halo properties and the simulation parameters is
captured well through a simple local multi-linear regression. The locality of
the regression is necessary as the trends are always nonlinear in halo mass
and often non-monotonic as well (Figures 1, 2, 3, etc.).

  * •

All halo properties have the shape of the mass-dependence altered via changes
to the simulation parameters. In general, changes to the SNe feedback
parameters $A_{\rm SN1/2}$ and the cosmology parameters $\Omega_{\rm m}$ and
$\sigma_{8}$ have a significant impact on scaling relations. The AGN feedback
parameters have a negligible impact.

  * •

We find $M_{\rm vir}\approx 10^{12}\,{\rm M}_{\odot}/h$ to be a key mass
scale, as it is associated with (i) the strongest dependence of $c_{\rm vir}$,
$\sigma_{\rm 3D}$, $M_{\rm star}$ on $A_{\rm SN2}$, and of $\sigma_{\rm 3D}$
on $\sigma_{8}$, (ii) a sign-flip in the dependence of $\lambda_{\rm h}$ on
$\Omega_{\rm m}$, (iii) the peak of the stellar mass fraction, (iv) a sign
flip in the dependence of $a_{\rm form}$ and $M_{\rm gas}$ on $A_{\rm SN1}$
and $A_{\rm SN2}$.

  * •

The $a_{\rm form}-M_{\rm vir}$ and $M_{\rm gas}-M_{\rm vir}$ scaling relation
show a highly non-monotonic dependence on $A_{\rm SN1}$ and $A_{\rm SN2}$, and
that increasing SNe feedback increases the baryon fraction (Figure 5). We
connect this to the nonlinear coupling of SNe and AGN, where the SNe takes gas
from the vicinity of the AGN and sends it further out in the halo (though not
out of the halo), which then prevents the AGN from ejecting material out of
the halo (Figure 2, 5).

  * •

The baryon imprints in $c_{\rm vir}$ — captured by the scaling of $c_{\rm
vir}$ with galaxy formation parameters — depend only loosely on any selections
in the halo properties $a_{\rm form}$, $\lambda_{\rm h}$, and $\sigma_{\rm
3D}$ (Figure 7). Instead, these selections have a large impact on the
cosmology-dependence of the $c_{\rm vir}-M_{\rm vir}$ relation.

  * •

These imprints are more sensitive to selections on $M_{\rm gas}$ or $M_{\rm
star}$ (Figure 9). In general, halo samples with higher $M_{\rm star}$ show a
weaker dependence of $c_{\rm vir}$ with $\Omega_{\rm m}$ and $\sigma_{8}$ and
stronger dependence with $A_{\rm SN2}$. Samples of higher $M_{\rm gas}$ follow
the same trend except for $\sigma_{8}$, where the higher $M_{\rm gas}$ sample
has a weaker dependence on $\sigma_{8}$.

  * •

The mean $c_{\rm vir}-M_{\rm vir}$ relation of the split subsamples generally
have the same shape but with slightly higher/lower amplitudes (Figure 8).
Selections on higher $\sigma_{\rm 3D}$ and $M_{\rm star}$ results in higher
$c_{\rm vir}$, since both properties are connected to the core of the halo
potential. Selecting high $a_{\rm form}$ or $\lambda_{\rm h}$ results in lower
$c_{\rm vir}$, while selecting high $M_{\rm gas}$ leads to lower $c_{\rm vir}$
at MW masses and high $c_{\rm vir}$ at cluster masses.

  * •

If we ignore the impact of selection effects, the subsequent error in modeling
of baryon imprints in the $c_{\rm vir}-M_{\rm vir}$ relation is of the order
$25\%$ for group-scale and cluster-scale halos (Figure 10).

Our results show that the baryon imprints on $c_{\rm vir}$ — defined in this
work through the slopes, $\alpha_{X}$, which parameterize the scaling of
$c_{\rm vir}$ with the galaxy formation parameters — have a clear dependence
on the selection function of the sample. The most impactful selections include
those on $M_{\rm gas}$ and $M_{\rm star}$, which are also the most
observationally relevant as they correspond to selections on X-ray luminosity
and Sunyaev Zeldovich decrement, or on the galaxy brightness/magnitude,
respectively. The importance of this dependence between baryon imprints and
halo selection will vary from application to application, but not accounting
for this selection-based effect in the baryon modeling can lead to biases.
However, the exact bias depends on the range of masses and physical scales
being probed in the analysis of interest.

We also build on previous works to show evidence of the non-linear coupling
between different feedback mechanisms. Thus, SNe feedback may be naively
thought to have no impact on the ejection of material out of the halo, whereas
in practice it has a significant impact. Such couplings can be useful in the
building and implementation of any semi-analytic models of galaxy formation,
where AGN feedback is the dominant modelling component. We have only explored
six simulation parameters (two cosmological parameters and four galaxy
formation parameters). Exploring the scaling with other processes/parameters
will uncover more non-linear couplings between different galaxy formation
processes, and shed light on the range of processes that can alter the impact
of AGN feedback in the halo’s internal dynamics. Suites such as the
28-parameter variation in Camels will enable such analyses (Ni et al. 2023).

The improvements in computational resources have enabled the production of
large simulation suites to study various physical phenomena, of both
astrophysical and cosmological nature (e.g., Villaescusa-Navarro et al. 2020,
2021; Kacprzak et al. 2023; Anbajagane et al. 2023c). This plethora of data
has moved forward our ability to construct data-driven, non-parametric models
that capture the complexities of the scaling relations. Summarizing an aspect
of these complexities, such as the non-linear mass dependence, can also help
build phenomenological models that can be utilized in analyses of
observational data. Our ability to build data-driven models of halo properties
will grow exponentially with the expected influx of both simulated and
observational data, and so studies of the robust construction and application
of such models are paramount to the efficient use of simulations in
constraining astrophysics and cosmology.

## Acknowledgements

We thank Andrew Hearin for discussions on this work, and Chihway Chang for
suggestions on an earlier version of this analysis and draft. We also thank
the Camels collaboration for making their data publicly available and easily
accessible to the community. DA is supported by NSF grant No. 2108168.

All analysis in this work was enabled greatly by the following software:
Pandas (McKinney 2011), NumPy (Van der Walt et al. 2011), SciPy (Virtanen et
al. 2020), and Matplotlib (Hunter 2007). We have also used the Astrophysics
Data Service (ADS) and arXiv preprint repository extensively during this
project and the writing of the paper.

## Data Availability

The tabulated Kllr parameters and the interpolator framework used in this work
are made available at https://github.com/mufanshao/CAMELS_ScalingRelations.
The public Camels data, including the Rockstar catalogs we use in this work,
can be found at https://camels.readthedocs.io.

## References

  * Anbajagane et al. (2020) Anbajagane D., Evrard A. E., Farahi A., Barnes D. J., Dolag K., McCarthy I. G., Nelson D., Pillepich A., 2020, MNRAS, 495, 686
  * Anbajagane et al. (2022a) Anbajagane D., Evrard A. E., Farahi A., 2022a, MNRAS, 509, 3441
  * Anbajagane et al. (2022b) Anbajagane D., et al., 2022b, MNRAS, 510, 2980
  * Anbajagane et al. (2022c) Anbajagane D., et al., 2022c, MNRAS, 514, 1645
  * Anbajagane et al. (2023a) Anbajagane D., et al., 2023a, arXiv e-prints, p. arXiv:2308.03863
  * Anbajagane et al. (2023b) Anbajagane D., et al., 2023b, arXiv e-prints, p. arXiv:2310.00059
  * Anbajagane et al. (2023c) Anbajagane D., Chang C., Lee H., Gatti M., 2023c, arXiv e-prints, p. arXiv:2310.02349
  * Aung et al. (2020) Aung H., Nagai D., Lau E. T., 2020, arXiv e-prints, p. arXiv:2012.00977
  * Bautista et al. (2021) Bautista J. E., et al., 2021, MNRAS, 500, 736
  * Baxter et al. (2021) Baxter E. J., Adhikari S., Vega-Ferrero J., Cui W., Chang C., Jain B., Knebe A., 2021, arXiv e-prints, p. arXiv:2101.04179
  * Behroozi et al. (2010) Behroozi P. S., Conroy C., Wechsler R. H., 2010, ApJ, 717, 379
  * Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57
  * Behroozi et al. (2019) Behroozi P., Wechsler R. H., Hearin A. P., Conroy C., 2019, MNRAS, 488, 3143
  * Beltz-Mohrmann & Berlind (2021) Beltz-Mohrmann G. D., Berlind A. A., 2021, arXiv e-prints, p. arXiv:2103.05076
  * Birrer et al. (2014) Birrer S., Lilly S., Amara A., Paranjape A., Refregier A., 2014, ApJ, 793, 12
  * Blumenthal et al. (1986) Blumenthal G. R., Faber S. M., Flores R., Primack J. R., 1986, ApJ, 301, 27
  * Bryan & Norman (1998) Bryan G. L., Norman M. L., 1998, ApJ, 495, 80
  * Bryan et al. (2013) Bryan S. E., Kay S. T., Duffy A. R., Schaye J., Dalla Vecchia C., Booth C. M., 2013, MNRAS, 429, 3316
  * Bullock et al. (2001) Bullock J. S., Dekel A., Kolatt T. S., Kravtsov A. V., Klypin A. A., Porciani C., Primack J. R., 2001, ApJ, 555, 240
  * Chisari et al. (2018) Chisari N. E., et al., 2018, MNRAS, 480, 3962
  * Cooray & Sheth (2002) Cooray A., Sheth R., 2002, Phys. Rep., 372, 1
  * Duffy et al. (2010) Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C., Battye R. A., Booth C. M., 2010, MNRAS, 405, 2161
  * Evrard et al. (2014) Evrard A. E., Arnault P., Huterer D., Farahi A., 2014, MNRAS, 441, 3562
  * Farahi et al. (2018) Farahi A., Evrard A. E., McCarthy I., Barnes D. J., Kay S. T., 2018, MNRAS, 478, 2618
  * Farahi et al. (2021) Farahi A., Nagai D., Chen Y., 2021, AJ, 161, 30
  * Farahi et al. (2022a) Farahi A., Nagai D., Anbajagane D., 2022a, arXiv e-prints, p. arXiv:2204.13578
  * Farahi et al. (2022b) Farahi A., Anbajagane D., Evrard A. E., 2022b, ApJ, 931, 166
  * Gao et al. (2005) Gao L., Springel V., White S. D. M., 2005, MNRAS, 363, L66
  * Gatti et al. (2020) Gatti M., et al., 2020, MNRAS, 498, 4060
  * Gatti et al. (2021) Gatti M., et al., 2021, arXiv e-prints, p. arXiv:2108.01600
  * Gebhardt et al. (2023) Gebhardt M., et al., 2023, arXiv e-prints, p. arXiv:2307.11832
  * Giri & Schneider (2021) Giri S. K., Schneider A., 2021, J. Cosmology Astropart. Phys, 2021, 046
  * Gnedin et al. (2004) Gnedin O. Y., Kravtsov A. V., Klypin A. A., Nagai D., 2004, ApJ, 616, 16
  * Grandis et al. (2023) Grandis S., Arico’ G., Schneider A., Linke L., 2023, arXiv e-prints, p. arXiv:2309.02920
  * Hearin et al. (2021) Hearin A. P., Chaves-Montero J., Becker M. R., Alarcon A., 2021, arXiv e-prints, p. arXiv:2105.05859
  * Heitmann et al. (2009) Heitmann K., Higdon D., White M., Habib S., Williams B. J., Lawrence E., Wagner C., 2009, ApJ, 705, 156
  * Henriques et al. (2013) Henriques B. M. B., White S. D. M., Thomas P. A., Angulo R. E., Guo Q., Lemson G., Springel V., 2013, MNRAS, 431, 3373
  * Hirschmann et al. (2014) Hirschmann M., Dolag K., Saro A., Bachmann L., Borgani S., Burkert A., 2014, MNRAS, 442, 2304
  * Hunter (2007) Hunter J. D., 2007, Computing in Science and Engineering, 9, 90
  * Kacprzak et al. (2023) Kacprzak T., Fluri J., Schneider A., Refregier A., Stadel J., 2023, J. Cosmology Astropart. Phys, 2023, 050
  * Knebe & Power (2008) Knebe A., Power C., 2008, ApJ, 678, 621
  * Kravtsov et al. (2018) Kravtsov A. V., Vikhlinin A. A., Meshcheryakov A. V., 2018, Astronomy Letters, 44, 8
  * Lau et al. (2010) Lau E. T., Nagai D., Kravtsov A. V., 2010, ApJ, 708, 1419
  * Lau et al. (2020) Lau E. T., Hearin A. P., Nagai D., Cappelluti N., 2020, arXiv e-prints, p. arXiv:2006.09420
  * Lee et al. (2022) Lee E., et al., 2022, MNRAS, 517, 5303
  * Li et al. (2008) Li Y., Mo H. J., Gao L., 2008, MNRAS, 389, 1419
  * Lim et al. (2021) Lim S. H., Barnes D., Vogelsberger M., Mo H. J., Nelson D., Pillepich A., Dolag K., Marinacci F., 2021, MNRAS, 504, 5131
  * Lu et al. (2015) Lu Z., Mo H. J., Lu Y., Katz N., Weinberg M. D., van den Bosch F. C., Yang X., 2015, MNRAS, 450, 1604
  * Macciò et al. (2008) Macciò A. V., Dutton A. A., van den Bosch F. C., 2008, MNRAS, 391, 1940
  * Machado Poletti Valle et al. (2020) Machado Poletti Valle L. F., Avestruz C., Barnes D. J., Farahi A., Lau E. T., Nagai D., 2020, arXiv e-prints, p. arXiv:2011.12987
  * Marinacci et al. (2018) Marinacci F., et al., 2018, MNRAS, 480, 5113
  * McKinney (2011) McKinney W., 2011, Python for High Performance and Scientific Computing, 14
  * Mitchell et al. (2020) Mitchell P. D., Schaye J., Bower R. G., Crain R. A., 2020, MNRAS, 494, 3971
  * Miyatake et al. (2023) Miyatake H., et al., 2023, arXiv e-prints, p. arXiv:2304.00704
  * Moster et al. (2010) Moster B. P., Somerville R. S., Maulbetsch C., van den Bosch F. C., Macciò A. V., Naab T., Oser L., 2010, ApJ, 710, 903
  * Moster et al. (2013) Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121
  * Muratov et al. (2015) Muratov A. L., Kereš D., Faucher-Giguère C.-A., Hopkins P. F., Quataert E., Murray N., 2015, MNRAS, 454, 2691
  * Naiman et al. (2018) Naiman J. P., et al., 2018, MNRAS, 477, 1206
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Nelson et al. (2018) Nelson D., et al., 2018, MNRAS, 475, 624
  * Nelson et al. (2019) Nelson D., et al., 2019, MNRAS, 490, 3234
  * Ni et al. (2023) Ni Y., et al., 2023, arXiv e-prints, p. arXiv:2304.02096
  * Ntampaka et al. (2019) Ntampaka M., et al., 2019, ApJ, 876, 82
  * Pandey et al. (2021) Pandey S., et al., 2021, arXiv e-prints, p. arXiv:2108.01601
  * Pandey et al. (2023) Pandey S., et al., 2023, MNRAS, 525, 1779
  * Pillepich et al. (2018a) Pillepich A., et al., 2018a, MNRAS, 473, 4077
  * Pillepich et al. (2018b) Pillepich A., et al., 2018b, MNRAS, 475, 648
  * Pillepich et al. (2019) Pillepich A., et al., 2019, MNRAS, 490, 3196
  * Porredon et al. (2021) Porredon A., et al., 2021, Phys. Rev. D, 103, 043503
  * Reddick et al. (2013) Reddick R. M., Wechsler R. H., Tinker J. L., Behroozi P. S., 2013, ApJ, 771, 30
  * Rodríguez-Monroy et al. (2022) Rodríguez-Monroy M., et al., 2022, MNRAS, 511, 2665
  * Rodríguez-Puebla et al. (2017) Rodríguez-Puebla A., Primack J. R., Avila-Reese V., Faber S. M., 2017, MNRAS, 470, 651
  * Schneider et al. (2019) Schneider A., Teyssier R., Stadel J., Chisari N. E., Le Brun A. M. C., Amara A., Refregier A., 2019, J. Cosmology Astropart. Phys, 2019, 020
  * Secco et al. (2022) Secco L. F., et al., 2022, Phys. Rev. D, 105, 023515
  * Shankar et al. (2017) Shankar F., et al., 2017, ApJ, 840, 34
  * Shao et al. (2023) Shao M. J., Anbajagane D., Chang C., 2023, MNRAS, 523, 3258
  * Singh et al. (2020) Singh P., Saro A., Costanzi M., Dolag K., 2020, MNRAS, 494, 3728
  * Springel et al. (2018) Springel V., et al., 2018, MNRAS, 475, 676
  * Stiskalek et al. (2022) Stiskalek R., Bartlett D. J., Desmond H., Anbajagane D., 2022, MNRAS, 514, 4026
  * Van der Walt et al. (2011) Van der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science and Engineering, 13, 22
  * Villaescusa-Navarro et al. (2020) Villaescusa-Navarro F., et al., 2020, ApJS, 250, 2
  * Villaescusa-Navarro et al. (2021) Villaescusa-Navarro F., et al., 2021, ApJ, 915, 71
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Wang et al. (2013) Wang L., et al., 2013, MNRAS, 431, 648
  * Wechsler et al. (2002) Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., Dekel A., 2002, ApJ, 568, 52
  * Weinberger et al. (2017) Weinberger R., et al., 2017, MNRAS, 465, 3291
  * Weinberger et al. (2018) Weinberger R., et al., 2018, MNRAS, 479, 4056
  * Yahil & Vidal (1977) Yahil A., Vidal N. V., 1977, ApJ, 214, 347
  * Yang et al. (2012) Yang X., Mo H. J., van den Bosch F. C., Zhang Y., Han J., 2012, ApJ, 752, 41
  * Zacharegkas et al. (2021) Zacharegkas G., et al., 2021, arXiv e-prints, p. arXiv:2106.08438