# Baryonic Imprints on DM Halos: The concentration–mass relation in the CAMELS
simulations

Mufan (Jon) Shao,1 Dhayaa Anbajagane,2, 3⋆ Chihway Chang2, 3

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 Author emails<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The physics of baryons in halos, and their subsequent influence on the total
matter phase space, has a rich phenomenology and must be well understood in
order to pursue a vast set of questions in both cosmology and astrophysics. We
use the Camels simulation suite to quantify the impact of four different
galaxy formation parameters/processes (as well as two cosmological parameters)
on the concentration–mass relation, $c_{\rm vir}-M_{\rm vir}$. We construct a
simulation-informed nonlinear model for concentration as a function of halo
mass, redshift, and 6 cosmological/astrophysical parameters. This is done for
two galaxy formation models, IllustrisTNG and Simba, using 1000 simulations of
each. We extract the imprints of galaxy formation across a wide range in mass
$M_{\rm vir}\in[10^{11},10^{14.5}]\,{\rm M}_{\odot}/h$ and in redshift
$z\in[0,6]$ finding many strong mass- and redshift-dependent features.
Comparisons between the IllustrisTNG and Simba results show the astrophysical
model choices cause significant differences in the mass and redshift
dependence of these baryon imprints. Finally, we use existing observational
measurements of $c_{\rm vir}-M_{\rm vir}$ to provide rough limits on the four
astrophysical parameters. Our nonlinear model is made publicly available and
can be used to include Camels-based baryon imprints in any halo model-based
analysis.

###### keywords:

galaxies: haloes – galaxies: statistics – dark matter

††pubyear: 2022††pagerange: Baryonic Imprints on DM Halos: The
concentration–mass relation in the CAMELS simulations–B

##  1 Introduction

The distribution of matter in our Universe contains information on both the
initial conditions of the Universe, i.e. the primordial density field, and the
subsequent cosmological and astrophysical evolution of the field to late times
(e.g., Allen et al., 2011; Weinberg et al., 2013; Kravtsov & Borgani, 2012).
Many efforts have been made in the past decades to measure this structure and
use it to understand processes of both cosmological and/or astrophysical
nature (e.g., Asgari et al., 2021; Pandey et al., 2021; Gatti et al., 2021;
Secco et al., 2022; Amon et al., 2022; Abbott et al., 2022; Tröster et al.,
2022; Schneider et al., 2022)

Efforts to understand the astrophysical impact in particular have grown
significantly in recent times, especially as it becomes clearer that our lack
of knowledge on these astrophysical processes will impact our ability to infer
precise and unbiased cosmological constraints from future measurements
(Semboloni et al., 2011; Chisari et al., 2018). The impact of baryons on
structure formation is strongest on “nonlinear” length scales corresponding to
the sizes of massive halos, and thus is intimately connected to galaxy
formation processes (e.g., Rudd et al., 2008; Schneider et al., 2019). The
influence of such non-gravitational processes on the baryon matter
distribution is indirectly felt by the dark matter (DM) distribution via the
gravitational coupling between the two components.

These astrophysical impacts on the matter distribution — either directly in
the case of gas and stellar distributions, and indirectly in the case of the
DM distribution — will be reflected in the density profiles of the halos.
These profiles generally follow an NFW form (Navarro et al., 1997), where for
a halo of a given mass the only free parameter is the halo concentration,
$c_{\rm vir}$. The changes in the matter distribution of a halo due to
astrophysics will result in changes in its $c_{\rm vir}$ value.

On nonlinear scales, many aspects of the observed structure can be understood
and modelled using these profiles of halos (Cooray & Sheth, 2002); in recent
times such halo-based approaches have been used to model the nonlinear galaxy
and shear fields (e.g., Zacharegkas et al., 2021), the outskirts of halos
(e.g., Shin et al., 2021; Anbajagane et al., 2022c), and the cross-correlation
between matter and thermodynamic fields (e.g., Vikram et al., 2017; Pandey et
al., 2019, 2021; Sánchez et al., 2022). Thus, the question of understanding
astrophysical impacts on nonlinear structure formation can be reformulated
into understanding such impacts on the halo concentration, or more generally
its impacts on the halo density profile (e.g., Schneider et al., 2019).

The impact of astrophysics on halo concentration has a long history of study,
first performed in smaller samples (e.g., Gnedin et al., 2004; Bryan et al.,
2013), and more recently in samples numbering millions of halos and spanning
multiple decades in halo mass (Ragagnin et al., 2019; Beltz-Mohrmann &
Berlind, 2021; Anbajagane et al., 2022a). A caveat however is that different
hydrodynamic simulations employ varying prescriptions for galaxy formation,
including the actual equations that are solved in the simulation. While this
variety in prescriptions has been useful in setting bounds on the impact of
astrophysics on halo properties, it does not enable a systematic study of
astrophysics in halos. There is indication that the form of the mass-dependent
impact of baryons on the halo matter distribution is similar across state-of-
the-art galaxy formation models but the quantitative behavior is still
different (Chua et al., 2021; Anbajagane et al., 2022b).

The Camels suite of hydrodynamic simulations (Villaescusa-Navarro et al.,
2021) was designed for quantitative studies of astrophysical impacts on
structure formation. The suite contains $\mathcal{O}(10^{3})$ simulations of a
relatively small volume of halos while also varying four astrophysical
simulation parameters — two corresponding to supernovae (SN) feedback, and two
corresponding to active galactic nuclei (AGN) feedback — in addition to the
two cosmological parameters most relevant for structure formation;
$\Omega_{\rm m}$, the fraction of matter energy density at $z=0$, and
$\sigma_{8}$, the root-mean-squared fluctuation of the $z=0$ density field
smoothed on the $8\,{\rm Mpc}\,h^{-1}$ length scale. These simulations are
also run with two different galaxy formation models, IllustrisTNG and Simba,
which then provides a way to estimate model-dependent differences in the
emergent astrophysical impacts.

A limitation of the Camels suite, however, is its small box size of
$L=25\,{\rm Mpc}/h$, which does not contain many halos of group and cluster
mass-scales ($M_{\rm vir}\sim 10^{13}-10^{14}\,{\rm M}_{\odot}/h$). The
astrophysical impacts on halos of these mass-scales are known to be distinct
from the impacts in less massive counterparts, so accurately capturing these
effects at high mass is of interest to this study. In addition, such massive
halos have a notable impact on the observed structure in our Universe given
they are large halos that form in the most overdense regions, even though
their number counts are far lower than their less massive counterparts.

In this work, we combine 1000 Camels simulations per galaxy formation model to
create an “uber” halo sample. This sample is then used to build a non-
parameteric model, constructed using a local linear regression technique, that
can predict the 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})$, which is a
function of halo mass, redshift, two SN feedback parameters, two AGN feedback
parameters, and two cosmological parameters. The choice of a local linear
regression allows us to capture non-linear behaviors with halo mass while
still using the familiar concept of power-law scaling relations. Such non-
linearity with mass is an expected feature in galaxy formation given the many
different mass scales of different processes, and this has been shown in
various works across wide mass ranges (e.g., Farahi et al., 2018; Chua et al.,
2019; Anbajagane et al., 2020; Chua et al., 2021; Anbajagane et al., 2022a;
Anbajagane et al., 2022b).

Our halo concentration model is then examined to understand the mass-dependent
impacts of astrophysical feedback on the $c_{\rm vir}-M_{\rm vir}$ relation;
specifically to identify certain features/trends in the relation and motivate
the physical processes that lead to such behavior. We do this over a wide
range in redshift, spanning from $z\sim 3$, near the peak of star formation,
all the way to the present epoch. While we generally will discuss results for
both the IllustrisTNG and Simba models, our more detailed discussions will
lean towards the IllustrisTNG results as we use the existing work of
Anbajagane, Evrard & Farahi (2022a, henceforth referred to as A22) on baryon
imprints in the fiducial IllustrisTNG simulations to identify the different
physical processes that may underlie the features found in our results here.
In fact, our work is an extension on an aspect of A22, where we now focus on
just the concentration–mass relation but use Camels to extract the
quantitative impact of galaxy formation on this relation.

To summarize, our goals in this work are to (i) construct and examine a model
of the concentration as a function of astrophysical processes, with the
dependence on redshift, halo mass, and cosmology folded in, (ii) to compare
and contrast results from different formation models provided by the Camels
suite (IllustrisTNG and Simba), and (iii) provide qualitative, observational
constraints on the astrophysical Camels simulation parameters using available
data on the $c_{\rm vir}-M_{\rm vir}$ relation.

This work is organized as follows. We present our simulation data and the
local linear regression method in §2. We showcase our model, including the
mass- and redshift-dependent impact of astrophysics, and the difference
between IllustrisTNG and Simba, in §3. We then discusses uses of our model,
including simulation “translations” and observational constraints, in §4. We
conclude and discuss future extensions of this work in §5.

##  2 Data & Methods

###  2.1 CAMELS Simulation Suite

Figure 1: The halo counts as a function of halo mass for the uber sample of
1000 Camels simulations at three redshifts $z\in\\{0,0.5,1\\}$. We only show
results for the IllustrisTNG formation model. The halo counts are computed in
bins of $0.1{\,\,\rm dex}$. The dashed lines show the theoretical expectation
using the halo mass function model of Seppi et al. (2021) and the theory
accounts for the 1000 different cosmologies combined in the uber sample. The
measured counts generally agree with the theory, validating that the uber
sample is still a cosmologically consistent sample.

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 with varying cosmological and astrophysical parameters, and
multiple galaxy formation models. The specific cosmological parameters are
$\Omega_{\rm m}$ and $\sigma_{8}$, which are the two parameters that have the
strongest influence on structure formation, and the astrophysical parameters
are $A_{\rm SN1},A_{\rm SN2},B_{\rm SN1},B_{\rm SN2}$ controlling supernovae
feedback, and $A_{\rm AGN1},A_{\rm AGN2},B_{\rm AGN1},B_{\rm AGN2}$ which
control active galactic nuclei (AGN) feedback. The specific meaning of each
parameter differs based on the chosen galaxy formation model — IllustrisTNG or
Simba — and these differences are detailed further below in §2.1.1 and §2.1.2
respectively. To highlight these differences explicitly, we will refer to
these parameters with different notations, as $A_{X}$ for IllustrisTNG and
$B_{X}$ for Simba respectively.

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 was designed specifically for training
machine learning algorithms but the dataset is generally well-suited for non-
parameteric studies of astrophysics and cosmology. There are different types
of runs within the full suite, of which three are of interest to us. The
Cosmic Variance (CV) runs are 27 simulations run with fiducial settings —
$\Omega_{\rm m}=0.3$, $\sigma_{8}=0.8$, and $A_{X}=B_{X}=1$, with $X\in\\{\rm
SN1,SN2,AGN1,AGN2\\}$ — and with varying initial conditions. 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. Finally, 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, with the points in parameter space
chosen via an optimal scheme called a Latin Hypercube, and this sampling
technique has become a frequent tool used in building emulators for cosmology
(Heitmann et al., 2009). Note that in addition to the IllustrisTNG and Simba
versions, Camels also has a dark matter only (DMO) version of every run (CV,
1P and LH). In this case, only the cosmological parameters matter and the
astrophysical parameter values have no impact.

We make use of the Rockstar halo catalogs (Behroozi et al., 2013a) from all
these simulations. The halo concentration in these catalogs is computed by
fitting the total matter density profiles of halos with an NFW form (Navarro
et al., 1997),

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

where $\rho_{s}$ and $r_{s}$ are the characteristic density and the scale
radius, respectively. The concentration is then derived as,

$c_{\rm vir}=R_{\rm vir}/r_{s},$ (2)

where the 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.

Note that the concentration and virial mass here are computed using all
available particles rather than just the DM alone. Thus, this work studies the
NFW profiles of the total matter component and not just the DM only. Though,
given the DM constitutes 85% of the matter content, the total matter
distribution — on length scales of interest to us, namely excluding the halo
core where the stellar mass distribution is most prominent — would to good
approximation be given by just the DM distribution.111Note that Rockstar
obtains $c_{\rm vir}$ by fitting density profiles down to the force softening
scale of the simulation, which is $\epsilon\approx 2\,{\rm kpc}$ (comoving)
for the simulations discussed here. In practice, the smallest informative
length scale of the profile is 2-3 times the force softening as densities on
scales close to $\epsilon$ are downweighted in the Rockstar fitting procedure.
The total matter distribution is observationally more relevant, and is
conveniently the definition used in these Rockstar catalog.

The total matter distribution has been shown to follow an NFW profile for
$r>0.1R_{\rm vir}$ but deviates from it below that radius scale due to the
presence of the stellar component (Duffy et al., 2010, see their Figure 5).
The Rockstar fitting algorithm uses the measured profile down to a minimum
scale set by the force softening scale of the simulation, $2\rm kpc$. This is
within $0.1R_{\rm vir}$ for all halos in our sample. So the Rockstar
concentration estimates will be affected by the fact that the true density
profile does not follow an NFW form on small-scales. Thus, we treat these
estimates as NFW-based approximations of the true density profile.

We also perform some quality cuts on the raw catalogs. In particular, the
Rockstar catalog shows a notable outlier halo population — 2% (5%) at $z=0$
($z=1$) — that have $c_{\rm vir}\lesssim 1$. We treat these objects as
outliers and remove them from our analysis. In short, this is done by forming
a histogram of $c_{\rm vir}$, identifying where the low-$c_{\rm vir}$ peak
begins to rise and removing all halos below that limit. We detail this
procedure more in Appendix A.

The uber halo samples — obtained by combining halos from all 1000 LH
simulations — span 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 and varies by redshift. In all our analyses, we set
the maximum $M_{\rm vir}$ for each redshift as the 10th most massive halo mass
in the uber sample at that redshift. The distribution of masses for some
redshifts is displayed in Figure 1, and shows that the uber sample follows a
cosmological distribution of halos, as evidenced by the $\approx 20\%$
agreement between the measured counts and the theoretical predictions. The
measurements come from aggregating counts from all 1000 LH simulations. The
theory line is constructed by obtaining predictions for each of the 1000
individual simulations — taken from the model of Seppi et al. (2021) as
implemented in Colossus (Diemer, 2018) — and aggregating them the same way we
aggregate the counts in the LH simulations. In this way the theoretical
prediction accounts for the cosmology variation across the 1000 LH
simulations.

We now briefly detail the specific meaning of each astrophysical parameter in
the two galaxy formation models, IllustrisTNG and Simba. More detail on the
simulation models can be found in Villaescusa-Navarro et al. (2021, see their
Section 3).

####  2.1.1 IllustrisTNG model

IllustrisTNG, which we will interchangeably refer to as TNG, 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 same model and
parameters from the fiducial simulations runs — which are described in detail
in Weinberger et al. (2017); Pillepich et al. (2018a) — are adopted in the
Camels simulations for this galaxy formation model. The resolution level of
the Camels simulation roughly corresponds to the TNG100-2 run.

The total energy injection rate (power) per unit star-formation, $e_{w}$, and
the SN wind speed, $v_{w}$, are modulated by two normalization factors,
$A_{\rm SN1}$ and $A_{\rm SN2}$, respectively. Thus $A_{\rm SN1}$ controls the
energy output rate, and $A_{\rm SN2}$ controls the wind speed. The mass
loading factor, $\eta$, which sets how much mass is transferred away through
the SN 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 nearby gas is thermally heated, and a radio/kinetic
mode, where the gas is given a momentum kick. The kinetic mode is turned on
during the lower accretion rate-phase of the SMBH, and is the dominant
mechanism for halos at the Milky Way-scale and above (Weinberger et al., 2018,
A22). Additionally, the kinetic feedback is not a continuous output into the
halo; instead the SMBH “accumulates” energy and expels it once an 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 control the rate of the bursts, while the latter controls the
speed of the kinetic feedback jets.

####  2.1.2 SIMBA model

Simba is another galaxy formation model (Davé et al., 2019) that contains a
rich variety of astrophysical processes, such as SN, SMBH feedback etc. The
implementation of the scaling parameters, $B_{X}$, in Simba is roughly similar
to $A_{X}$ in IllustrisTNG but has some differences in the details of the
implementation.

While $B_{\rm SN2}$ continues to correspond to the SN wind velocity like
$A_{\rm SN2}$, $B_{\rm SN1}$ directly scales the mass-loading factor now (as
opposed to in TNG where both SN parameters controlled the mass-loading). There
are no parameters that directly scale the energy output of SN events. We
stress that the specific equations governing the winds and the feedback still
differ significantly between IllustrisTNG and Simba. Of particular note is
that the Simba model is informed more strongly by results from existing ultra
high-resolution galaxy-scale zoom-in simulations (Muratov et al., 2015).

As for the SMBH model, the $B_{\rm AGN1}$ parameter in Simba sets the scale of
the momenta kicks provided to the gas through the kinetic feedback mechanism,
which is different from the energy accumulation rate that $A_{\rm AGN1}$
controls in TNG. Then, $B_{\rm AGN2}$ sets the velocity scale of the jets that
result from the kinetic feedback, and this is similar to the role of $A_{\rm
AGN2}$ in TNG. Now, once $B_{\rm AGN2}$ sets the jet velocity scale, the value
of $B_{\rm AGN1}$ determines the mass in the AGN outflows that is required to
obtain the momentum scale set by $B_{\rm AGN1}$.

###  2.2 Observational $c_{\rm vir}-M_{\rm vir}$ relation data

We use observational data for the $c_{\rm vir}-M_{\rm vir}$ relation to place
constraints on the astrophysical parameters in Camels ; namely to ask whether
the constraints lie above/below the fiducial values. We take data mainly from
two existing works — Johnston et al. (2007); Mandelbaum et al. (2008) — both
of which measured halo concentration from the Sloan Digital Sky Survey (SDSS,
York et al., 2000). The values for $c_{\rm vir}$ were inferred from average
lensing profiles of ensembles of halos. Mandelbaum et al. (2008) impose scale
cuts to extract parts of the profile that can be more easily compared to
N-body simulations, while Johnston et al. (2007) used much smaller scales but
accounted for the baryonic components when extracting the halo concentration.
The estimates were made at various redshifts across works but can be brought
to a common redshift of $z=0$ using a simple scaling as was done in Child et
al. (2018). We also convert the chosen overdensity definitions of the
measurements, using the method of Hu & Kravtsov (2003), to match the virial
mass definition of Bryan & Norman (1998) which is the default output
definition in the Rockstar halo catalog. This conversion is done by using the
mean mass and concentration of a given datapoint to compute the offets to
apply in mass and concentration and thus convert the definition. Notably, we
do not incorporate any additional uncertainties on concentration and mass that
arise from performing this procedure. This is adequate for the precision level
of our analysis.

We also use three additional points — one from Brimioulle et al. (2013) at
galaxy-scales using the Canada-France-Hawaii telescope (CFHT) data, and two
points from Mandelbaum et al. (2006) at group scales, using an older version
of SDSS data than the ones used in Johnston et al. (2007); Mandelbaum et al.
(2008).

The two works of Johnston et al. (2007) and Mandelbaum et al. (2008) both use
SDSS data but have notable differences. They use different cluster catalogs
that vary by $1.5{\,\,\rm dex}$ in the mass range they cover. The latter also
uses additional galaxy samples beyond those used by the former. There are
furthermore many analysis choices that differ between the two in how they
extract profiles, and how they account for different effects in the
measurements. A more detailed description can be found in Mandelbaum et al.
(2008, Section III. C)

When using this data to constrain model parameters, we must know the
covariance of the data in order to compute the chi-squared, $\chi^{2}=\delta
x\mathcal{C}^{-1}\delta x$, where $\delta x$ is the residual between
measurement and prediction. In our case, we assume our datapoints are
uncorrelated, i.e. $\mathcal{C}$ is a diagonal matrix. This is certainly not
true given some datapoints across different works have some shared data
underlying them (notably, SDSS), but our limitation is the lack of any
information on this point-to-point covariance. However, as we will stress
later, our goal with using this observational data is not to make precision
estimates of the parameter constraints but rather use the data to make
qualitative statements, i.e. do the extracted parameter posteriors prefer
values greater/lesser than the fiducial parameter values. For such a purpose,
we find it suffices to approximate the covariance as a diagonal matrix,
especially alongside other approximations above such as the mass definition
conversion and the redshift re-scaling.

###  2.3 Kernel-Localized Linear Regression

Traditionally, scaling relations are fit with a simple power law, i.e. a log-
linear relation, or are extended to include a broken power-law model with two
power-law indices instead of one. This is a perfectly sound technique, but
does not capture all the non-linearity that can be present in the scaling
relations. Using simple power-laws to represent scaling relations is well
motivated in gravity-only simulations where there are very few unique scales
in the problem that can break self-similar scaling. With the inclusion of
galaxy formation, however, many scales now arise and lead to strong deviations
from the power law behavior in the matter phase space (A22). Given the surplus
of simulation data available now, one can explore models with more freedom
that can accurately capture such non-linearities without assuming fixed
functional forms.

A common model-independent approach has been to use machine learning (ML)
techniques, such as random forest and neural networks (e.g., Machado Poletti
Valle et al., 2020; Stiskalek et al., 2022). However, this has the issue of
being difficult to interpret and extract physical explanations, though there
are approaches one can take to make it easier (e.g., Ntampaka et al., 2019;
Machado Poletti Valle et al., 2020, A22). An alternative, is to go back to our
well-known, well-understood linear models but extend them 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 is weighted
according to a kernel, which for studying halo scaling relations is most often
implemented as a gaussian kernel in $\log_{10}M_{\rm vir}$, and a linear
regression is fit to the kernel-weighted data. As the kernel is moved over the
log-mass direction, the method maps out the non-linear scaling relation using
connected linear pieces. The technique has already been used extensively to
show nonlinearity in different halo property–mass scaling relations (Farahi et
al., 2018; Anbajagane et al., 2020; Anbajagane et al., 2022a; Anbajagane et
al., 2022b). In this work, we use Kllr to study a multivariate scaling
relation that is still modelled as being locally linear in halo mass; a
version of Kllr that was first introduced and utilized in A22.

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

$\displaystyle\log_{10}c_{\rm vir}=$ $\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,$ (3)

x where $\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},Y_{\rm SN1},Y_{\rm AGN1},Y_{\rm SN2},Y_{\rm AGN2}\\}$, where
$Y\in\\{A,B\\}$ depending on whether the model is for IllustrisTNG or Simba.
The slopes and intercepts are all functions of halo mass and redshift. The
concentration, $c_{\rm vir}$, is thus modelled as a function of mass,
redshift, and 6 cosmological/astrophysical parameters. The dependence on the
cosmological/astrophysical parameters is notably also a function of mass. When
we analyze the DMO simulations, we drop all astrophysical terms in equation
(2.3) and keep only the cosmology ones. When analyzing the CV simulations or
any individual, single simulation, we drop all terms with $\alpha_{X}$ (since
none of the parameters $X$ vary in these scenarios) and only keep the slope
with mass.

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$. We compute the mean
$c_{\rm vir}$ at 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 (<10). The
Kllr technique uses a Gaussian kernel in $\log_{10}M_{\rm vir}$, with a width
of $0.3{\,\,\rm dex}$.222In general, wider kernels smooth out more features,
and narrower kernels lead to noisy estimates of the scaling relation. Our
choice of $0.3{\,\,\rm dex}$ is guided by the semi-optimized choices of the
previous works that used Kllr. We obtain uncertainties on all the Kllr
parameters of equation (2.3) by obtaining the regression coefficients for 100
bootstrap realizations of the data. Thus, we have 100 versions of equation
(2.3) that can then be used to quantify the uncertainty in our Kllr model
predictions.

Note that Kllr is by no means a unique method for capturing non-linearity in
multivariate scaling relations via a method that is not an ML technique
(random forest, neural networks etc.). For example, one could utilize Gaussian
processes to perform the regression instead; similar to Kllr, the Gaussian
processes assumes no apriori information about the form of the scaling
relation, and will capture non-linearities present in the data. Indeed, both
Kllr and gaussian processes have already been used interchangeably to extract
non-linear relations in halo quantities (e.g., Farahi et al., 2021, 2022a).

##  3 Results I: Baryon-dependent $c_{\rm vir}-M_{\rm vir}$ relation

We present and examine the Kllr-based model built using the various Camels
simulations. In particular, we introduce the different mean relations (§3.1),
detail the dependence of $c_{\rm vir}-M_{\rm vir}$ on the astrophysical and
cosmology parameters and the physical processes they originate from (§3.2),
differences between TNG and SIMBA in the baryon imprints found in each (§3.3),
the change with redshift of the astrophysical and cosmoogical dependence in
$c_{\rm vir}$ (§3.4), and finally the difference in behaviors in the
hydrodynamic simulations with galaxy formation and their dark matter only
counterparts (§3.5).

All results in this section come from our Kllr model results. Uncertainties
are obtained from 100 bootstrap resamplings and represent the 68% confidence
interval.

###  3.1 Fiducial $c_{\rm vir}-M_{\rm vir}$ relation in IllustrisTNG and Simba

Figure 2: The $z=0$ Kllr scaling relation of $\rm c_{\rm vir}-M_{\rm vir}$ for
TNG, SIMBA, and DMO runs. The bands show the 68% confidence interval obtained
from 100 bootstrap realizations for the data and for the model. The solid
lines are predictions for the fiducial simulation parameters using our model
that was constructued using the 1000 LH sims. The dashed lines are the Kllr
scaling relations estimated on a combined halo sample from the 27 CV
simulations. TNG and SIMBA agree well at higher halo masses, with both
deviating significantly from DMO, and have bigger differences at low masses.
In all cases, the prediction from the model agrees with the CV results.

Since the efficiency of dissipative processes peaks near MW masses, the
concentration in the FP runs is boosted at that scale. However, competing
feedback processes, from SN at lower masses and SMBHs at higher masses,
disperse baryons out of the core and into the outer halo, thus decreasing the
DM concentration at these mass scales. Above $\sim 10^{14}\,{\rm M}_{\odot}$,
the shift in mean concentration is very small.

As mentioned before, the impact of baryons on halo concentration has a long
history of study in simulations. The first discovery was that the presence of
gas — a component that could radiate away energy via cooling mechanisms — led
to adiabatic contraction, and this change in the gas distribution, through
gravitational back-reaction on the DM component, caused halos to be more
concentrated (Blumenthal et al., 1986; Gnedin et al., 2004) and others have
found the effect amplitude can depend on the history of how baryons assmebled
in the halo (Abadi et al., 2010; Tissera et al., 2010; Pedrosa et al., 2010;
Artale et al., 2019; Forouhar Moreno et al., 2021). The subsequent
implementation of stellar feedback, and then AGN feedback, in the simulations
led to the discovery that the halo is made more diffuse through these
explosive processes and this has been shown before for cluster-scale haloes
(Duffy et al., 2010; Cui et al., 2016), MW-scale haloes (Pedrosa et al., 2010;
Artale et al., 2019; Forouhar Moreno et al., 2021) and for dwarf galaxy-scale
haloes (Wetzel et al., 2016; Wheeler et al., 2019). More recently, A22 showed
the full mass-dependence of these effects in TNG across 6 decades in halo
mass, from dwarf galaxies to galaxy clusters, and that these different effects
interact to create a non-trivial “wiggle” feature in the concentration–mass
relation. Anbajagane et al. (2022b) showed — using a different tracer of the
halo potential — that such mass-dependence from baryon imprints can be
generally found in many of the latest hydrodynamic simulations, such as
IllustrisTNG, Magneticum, Bahamas and The300 (Nelson et al., 2018; Pillepich
et al., 2018b; Springel et al., 2018; Marinacci et al., 2018; Naiman et al.,
2018; Hirschmann et al., 2014; McCarthy et al., 2017; Cui et al., 2018).
However, all of the existing works summarizing large samples of halos have
made comparisons between three to four specific galaxy formation models. In
this work, we have extended this to utilize $\mathcal{O}(10^{3})$ variations
of a given galaxy formation model, and we explicitly compare the impact of
galaxy formation in a mass-dependent way.

We start in Figure 2, showing the comparison of several fiducial $\rm c_{\rm
vir}-M_{\rm vir}$ relations at $z=0$; we show TNG (blue), SIMBA (green), and
DMO (gray) versions. The solid lines are the predictions from our Kllr-based
models — which were developed using the 1000 LH simulations — given fiducial
values for the six input parameters of the simulations. The dashed lines are
the Kllr scaling relation of a combined halo sample from all 27 CV
simulations. We see that there are notable differences between the TNG, SIMBA,
and DMO version. At group and cluster mass-scales, TNG and SIMBA both have a
lower concentration than the DMO result, but this difference is subdued when
looking toward the most massive objects, and this behavior is consistent with
results from A22 (see their Figure 3). It is also promising to see that both
SIMBA and TNG agree at the highest masses. This implies the fiducial model of
the feedback prescriptions in TNG and SIMBA leads to broadly similar effects
on $c_{\rm vir}$. At smaller masses, $M_{\rm vir}<10^{12}\,{\rm M}_{\odot}/h$,
the scaling relations have notable differences.

In all three cases, our Kllr model has captured the overall dependence of
$c_{\rm vir}$ on simulation parameters in each version, as shown by the
agreement between the solid and dashed lines. This is not a trivial statement;
even though both dashed lines and solid lines ultimately derive from using
Kllr, they are coming from completely different sets of simulations. The solid
lines are a model built using 1000 LH simulations with varying cosmological
and astrophysical parameters, whereas the dashed lines are 27 CV sims at a
single fiducial cosmology. Thus the agreement in Figure 2 is one validation of
our methodology.

###  3.2 Astrophysics and cosmology dependence in IllustrisTNG

Figure 3: The $z=0$ Kllr slopes of concentration $c_{\rm vir}$ with halo mass
$M_{\rm vir}$ and the 6 astrophysical and cosmological parameters. The bands
show the 68% confidence interval. The left (right) panel shows results from
the TNG (SIMBA) sample. In all cases, the slopes have strongly nonlinear
behaviors with halo mass and are mostly non-monotonic. The cosmological
parameters $\Omega_{\rm m}$ and $\sigma_{\rm 8}$ have the strongest influence
on halo concentrations, and are followed by the astrophysical parameters
$A_{\rm SN1}$ and $A_{\rm SN2}$. The AGN parameters generally have a non-zero
but low influence on $c_{\rm vir}$.

Figure 3 shows the key component of our model — the slopes of $c_{\rm vir}$
with $M_{\rm vir}$ and the 6 simulation parameters, for both TNG (left) and
SIMBA (right). We first provide some intuition for how to interpret the
slopes. Looking at the expression for our Kllr-based model in equation (2.3),
we see there is an explicit mass-dependence for the slopes. In the limit that
there is no such dependence, i.e. the slopes are independent of mass, they
will be flat, horizontal lines in Figure 3, and indicate that a simple power-
law is adequate to represent the relationship between the two variables
connected by the given slope parameter. A slope that changes linearly in
Figure 3 shows preference for a linear mass term in the power-law index, e.g.,
$X^{\alpha_{X}+\beta_{X}M_{\rm vir}}$, where $X$ is any of the 7 independent
variables in the regression (halo mass or the 6 simulation parameters).

Going back to Figure 3, we find strong mass-dependent features in both the TNG
and SIMBA models, an expectation that was motivated by the existing works on
baryon imprints that were described earlier. We see that the slope of $c_{\rm
vir}$ and $M_{\rm vir}$ has a clear, non-monotonic relation with a minimum
roughly around the Milky Way (MW) mass-scale, though the exact location
appears to vary by a factor of 2-3 in mass. This non-monotonic behavior was
expected given the form of the scaling relations in Figure 2. The slopes for
both TNG and SIMBA also have the same form of approaching zero at low and high
mass-limits and having a distinct minima at MW mass-scales. It has been
identified that this form arises in TNG due to gas cooling and thermal
feedback competing to dominate the halo energetics (A22), with some indication
from just group and cluster-scales that other hydrodynamic simulations show
similar features (Anbajagane et al., 2022b), and we now explicitly confirm the
same qualitative behavior in SIMBA across three and a half decades in halo
mass.

Moving on to the astrophysics dependence, we focus on the TNG model first. We
find that the SN parameters have the strongest astrophysical impact on the
halo concentration, as evidenced by their large slope amplitudes compared to
the AGN slopes, and also have the most mass-dependent effect given their
slopes have clear non-monotonic features (including zero-crossings) with halo
mass. The SN wind speed parameter, $A_{\rm SN2}$, is the most impactful,
particularly at the MW mass-scale. Following A22 (see their Figure 7), our
expectation is for SN feedback to matter more at masses smaller than this MW-
scale, as the MW mass-scale is where gas cooling and thermal (quasar) feedback
become comparable; stellar feedback is subdominant.

Thus, it may seem unexpected that the maximum impact of SN winds is at this
mass scale where SN feedback is subdominant and gas cooling is dominant.
However, our assumption underlying this statement is that the impact of SN is
solely in how their feedback blows out the matter component and makes the halo
diffuse. SN, and their winds in particular, are also a necessity for
dispersing heavier elements into the halo gas component and increasing the
metallicity through the halo. Increasing metallicity enhances gas cooling,
which in turn increases the concentration through adiabatic contraction. For
example, Strickland & Heckman (2009, see their Equation 15) show that the gas
cooling rate is proportional to the mass loading factor. In the case of TNG,
increasing $A_{\rm SN2}$ (with all other parameters fixed) causes the wind
speed to increase, and the mass loading factor to drop. It is therefore
possible for the SN model to have a noticeable baryonic imprint in halos over
a wider mass-range due to its additional impact on gas cooling. The full set
of processes that cause this specific mass-dependence is unclear and is beyond
the scope of this work.

The other SN parameter, $A_{\rm SN1}$ also has a decent impact on $c_{\rm
vir}$, albeit at lower masses. This is the more “traditional” result of SN
blowing material out the halo potential, and is reasonable behavior in CAMELS
as $A_{\rm SN1}$ controls the energy output of SN and thus directly impacts
the strength of these explosive events. While there appears to be a mild
positive slope/correlation between $A_{\rm SN1}$ and $c_{\rm vir}$ at the
highest masses, the statistical significance of this slope deviating from 0 is
$\lesssim 2\sigma$, and so we ignore this feature.

Next, both AGN parameters – $A_{\rm AGN1}$ and $A_{\rm AGN2}$ — have a weak
impact on $c_{\rm vir}$. In particular, $A_{\rm AGN1}$, which controls the
rate at which the kinetic feedback mode accumulates energy before expelling it
as a jet, has no visible impact on this relation. This is once again not
particularly surprising as results in A22 (see Figure 7) show that in TNG the
kinetic (radio) feedback is at best comparable with other galaxy formation
energetics. Thus, we expect weak to no dependence here. Our slope in this case
is consistent with 0 at $\lesssim 1\sigma$.

The second parameter, $A_{\rm AGN2}$, changes the “burstiness” of the kinetic
feedback and does have an impact on $c_{\rm vir}$. While this may seem
counter-intuitive, since both AGN parameters impact the same kinetic feedback
mode and $A_{\rm AGN1}$ clearly has no impact on $c_{\rm vir}$, there are some
key differences. Recall our previous discussion in Section 2.1.1 that the
kinetic feedback in TNG works by accumulating energy until a set energy scale
is reached and then by expelling it as kinetic energy in the form of momentum
kicks. This is what creates the “burstiness” in the kinetic jets. $A_{\rm
AGN1}$ changes how quickly/slowly the SMBH achieves this energy threshold to
then create a jet. Whereas $A_{\rm AGN2}$ changes the exact energy scale of
the jet, which also sets the minimum kinetic energy injected into the halo.
Although kinetic feedback is subdominant in the fiducial TNG model, varying
$A_{\rm AGN2}$ changes the strength of individual feedback events, which can
alter this behavior. As a toy example, 100 temporally spaced out AGN events of
$10^{50}\rm ergs$ can have a subdominant/incremental impact compared to a
single event of $10^{52}\rm ergs$. The parameter $A_{\rm AGN2}$ sets where we
are on this scale of “many, weaker events” and “fewer, stronger events”. We
also stress that the dependence of $c_{\rm vir}$ on $A_{\rm AGN2}$ is still
mild for the mass ranges explored in this work, though it is possible it is
stronger at the most massive halo scales $M_{\rm vir}>10^{15}\,{\rm
M}_{\odot}/h$ that we were unable to probe here.

Finally, we come to the slopes of the different cosmology parameters
$\Omega_{\rm m}$ and $\sigma_{8}$, which have the strongest impact, i.e. the
highest slopes, on $c_{\rm vir}$. The slopes do have some mass-dependence,
especially for $\Omega_{\rm m}$ whose slope shows a “wiggle”-type feature.
When running a full hydrodynamic simulation, we expect the cosmology
dependence is “convolved” with the suite of baryon evolution effects, and
given galaxy formation physics has many different mass/energy/length scales,
we will expect the cosmology dependence to also become strongly scale-
dependent, with potentially non-monotonic features. This is exactly what we
see. However, to zeroth-order the slopes’ change with mass is consistent with
that seen in DMO simulations, and we will discuss this further in Section 3.5.

Looking at the slopes between $c_{\rm vir}$ and the cosmology parameters, one
may wonder why the astrophysical dependence is deemed important given the
amplitude of slopes between $c_{\rm vir}$ and astrophysical parameters can be
a factor of 2-3 lower than that between $c_{\rm vir}$ and the cosmological
parameters. However, we know the value of cosmological parameters very well.
Measurements of $\sigma_{8}$ and $\Omega_{\rm m}$ have resulted in 1-2%
constraints (e.g., Planck Collaboration et al., 2016; Asgari et al., 2021;
Abbott et al., 2022). This causes a shift of 1% (0.5%) in $c_{\rm vir}$ due to
any shifts in $\sigma_{8}$ ($\Omega_{\rm m}$) resulting from the uncertainty.
On the other hand, we do not have a good handle for priors on the
astrophysical parameters. Going by the range of values chosen in Camels, they
can vary by $150-250\%$, and so the concentration can change by up to $100\%$.
This is just a back-of-the-envelope estimate meant to motivate why the
astrophysical impacts are still important (and in fact, likely dominant) here,
even though the slopes with cosmological parameters have larger amplitudes.

###  3.3 Comparison of astrophysics and cosmology dependence in IllustrisTNG
and Simba

As was described before, the TNG and SIMBA galaxy formation models are
different, and the Camels simulation parameters modify different equations
within each model. Thus, the connection between $c_{\rm vir}$ and the
simulation parameters is altered across models. We detail the differences
here.

First, the slopes with $M_{\rm vir}$ have different detailed behaviors. This
difference is a combination of all differences in the physics of each galaxy
formation model. We do not focus on the details of these differences in this
work as it is difficult to disentagle all the input effects to make causal
statements. The general point we make, however, is that the mass-dependence of
the slope with $M_{\rm vir}$ in both simulations is qualitatively the same —
the slope is closer to zero at both high and low masses, but has a strongly
negative value at Milky Way mass scales.

Now focusing on the 6 simulation parameters, we first find that the cosmology
dependence is quantitatively — and to some extent qualitatively — different
between the TNG and SIMBA models. This is expected given our previous
discussion about how these slopes with $\Omega_{\rm m}$ and $\sigma_{8}$
capture the effect of cosmology convolved with all the physics of a galaxy
formation model. As the galaxy formation model changes, the cosmological
dependence will be affected as well. In general, both the TNG and SIMBA models
show the impact of cosmology is more strongly felt at lower masses,
particularly in SIMBA, and we show below that this is a feature that exists in
the DMO results as well (Section 3.5). The “wiggle” feature in the
$\Omega_{\rm m}$ slope of TNG is not seen in SIMBA.

Moving onto the four astrophysical parameters, we find more differences here.
Notably, the slope with $B_{\rm SN1}$ takes a positive value in SIMBA, even at
small masses, whereas for $A_{\rm SN1}$ in TNG this slope was negative. This
can once again be attributed to cooling — $B_{\rm SN1}$ in SIMBA directly
controls the wind mass-loading (whereas in TNG the combination of $A_{\rm
SN1}$ and $A_{\rm SN2}$ controlled this), so as $B_{\rm SN1}$ is increased,
the cooling rate increases (e.g., Strickland & Heckman, 2009, see their
Equation 15). This in turn causes a higher $c_{\rm vir}$, leading to the
positive slope between $c_{\rm vir}$ and $B_{\rm SN1}$. Next, $B_{\rm SN2}$
alters the wind speed in SIMBA, just like $A_{\rm SN2}$ in TNG, and so the
behavior is qualitatively similar in the mass-dependence as $\\{A,B\\}_{\rm
SN2}$ in both TNG and SIMBA are consistent with zero at high mass and become
more negative towards low mass. Differences in the actual amplitude of the
slopes will arise from more technical differences in the implemented model and
we do not focus on those here.

Moving to the AGN parameters, we find once again that the parameters do not
impact the $c_{\rm vir}-M_{\rm vir}$ relation much, as evidenced by slope
amplitudes being close to zero. Of note is the fact that in SIMBA, the $B_{\rm
AGN1}$ parameter changes the kinetic energy scale of the jets, whereas $A_{\rm
AGN1}$ in TNG changes the energy accumulation rate. The $B_{\rm AGN2}$ in
SIMBA controls the speed of kinetic feedback jets, which is directly analogous
to the $A_{\rm AGN2}$ parameter in TNG. Thus, we see that $B_{\rm AGN2}$ is
the more impactful of the two SIMBA AGN parameters, is mostly impactful at
high halo mass, and works to reduce the concentration.

###  3.4 Redshift dependence in IllustrisTNG

Figure 4: The Kllr slopes between halo concentration $c_{\rm vir}$ and the six
astrophysical and cosmological parameters (denoted by text in each panel) at
redshifts $0<z<3$. The slopes are interpolated across redshift to get the
values at the exact redshift of interest. The bands show the 68% confidence
intervals. There is a clear and strong mass-dependence in all parameters, and
a significant redshift evolution in most.

We have thus far described the behavior of the relations between halo
concentration and the simulation parameters at $z=0$. However, these
relationships will not be constant over redshift, given that galaxy formation
processes vary strongly across redshift (see Figure 1 in Weinberger et al.,
2018); for example, star formation efficiency is known to peak at $z=2$ and is
not constant over time (e.g., Behroozi et al., 2013b, Figure 3).

In Figure 4, we show the slopes between $c_{\rm vir}$ and the six simulation
parameters in IllustrisTNG as a function of redshift for the range $0<z<3$.
The results for Simba are extracted in Appendix B. In practice, to get these
results, we use the 34 Kllr models built at different redshifts to create an
interpolator333We use the interpn function from the scipy package (Virtanen et
al., 2020), and interpolate over a regular, rectangular grid of the mass and
redshift input variables. over $\log_{10}M_{\rm vir}$ and $\log_{10}(1+z)$.
This interpolator is used to get predictions at the exact mass and redshifts
of interest. We use the 100 realizations of the Kllr model to create multiple
interpolators, which then give us multiple predictions and quantify
uncertainties in these predictions.

In general, the shape of the slopes’ mass-dependence seen at $z=0$ is
preserved across a wide range in redshift, though the exact values of the
slopes differ significantly. We now examine the specific behavior in each
slope.

In $\boldsymbol{A_{\rm SN1}}$, the slopes become more negative with time
across all masses. This is reasonable as star formation becomes much more
efficient around $z\sim 2$, so the energy scale of SN feedback is more
relevant/impactful at those redshifts. We see a characteristic dip at $M_{\rm
vir}\approx 10^{11.5}\,{\rm M}_{\odot}/h$ that is mostly constant over
redshift. The slopes also have a plateauing phase around $10^{13}\,{\rm
M}_{\odot}/h$ which is seen at all redshifts, but appears to happen at higher
masses at higher redshifts. Multiple redshifts also show the characteristic
rise after the plateauing that gets stronger towards galaxy cluster mass-
scales.

In $\boldsymbol{A_{\rm SN2}}$, there appear to be two different mass regimes.
At low masses, the slopes show the same shape for the mass-dependence, but get
more shallow with redshift. At higher masses, the mass-dependence is nearly
constant with redshift, and appears to tend towards positive values at the
highest masses. These two regimes are split at around the same mass of
$10^{12.8}\,{\rm M}_{\odot}/h$ for all redshifts. For the low-mass dependence,
the exact minima of the slopes shifts in mass as well; at higher redshift, it
occurs at higher mass.

In $\boldsymbol{A_{\rm AGN1}}$ and $\boldsymbol{A_{\rm AGN2}}$, the slopes are
mostly consistent with zero across all mass ranges and redshifts. The
strongest impact of $A_{\rm AGN1}$ is around $M_{\rm vir}\approx
10^{12.5}\,{\rm M}_{\odot}/h$, where the slope deviates the most from zero;
but this is only at $z=0$, as this slope returns back to 0 at $z=3$. For
$A_{\rm AGN2}$, the slopes above $M_{\rm vir}>10^{13}\,{\rm M}_{\odot}/h$ move
away from zero and become negative. At the highest masses, the level of
impact, shown by the amplitude of the slope, is comparable to the maximum
impact $A_{\rm SN1}$ has. The lack of redshift dependence is consistent with
how, at the highest halo mass, kinetic AGN feedback in TNG is quite constant
across redshift (Weinberger et al., 2018, Figure 1).

In $\boldsymbol{\Omega_{\rm m}}$, we see the strongest redshift dependence.
Once again, the shape of the mass-dependence is mostly constant, but the
amplitude varies significantly, where at higher redshift there is a lower
impact on $c_{\rm vir}$. There is a peak in the slopes at $M_{\rm
vir}=10^{11}\,{\rm M}_{\odot}/h$ that is prominent at low redshifts and
diminishes at higher redshifts. The slope also seems to tend towards
$\alpha\sim 0.25$ at the highest masses, and this is true for all redshifts
where we can probe the mass range $M_{\rm vir}>10^{13.5}\,{\rm M}_{\odot}/h$.
There is also a local minima at $M_{\rm vir}=10^{11.5}\,h^{-1}\,{\rm
M}_{\odot}$, which is interestingly the same mass-scale where we found a
minimum in the $A_{\rm SN1}$ slopes. For this work, we do not pursue any
causal link connecting the two features.

Finally, in $\boldsymbol{\sigma_{8}}$, we have differences in the shape of the
mass-dependence across redshift. At low masses, there is a clear plateau in
the slopes at all redshifts, with the plateau value decreasing with redshift.
At higher masses, the slope decreases more steeply with mass as we increase
redshift, though the increase in steepness reaches an asymptotic value after
$z=1$, especially for $M_{\rm vir}>10^{12.5}\,{\rm M}_{\odot}/h$.

In general, Figure 4 shows a rich phenomenology for how different galaxy
formation processes impact $c_{\rm vir}$ as a function of redshift. It also
provides an example of the scale-dependence of galaxy formation discussed
previously. We have thus far seen this scale-dependence as a non-monotonic
mass-dependence of the slopes, but we now explicitly see this in the redshift-
dependence as well. Notably, there is no simple redshift-based rescaling of
$c_{\rm vir}$ or the slopes $\alpha_{X}$ that can give us the right redshift-
dependence for all six parameters at once. Each parameter must be treated
differently, given it is addressing different galaxy formation processes which
are all mass- and redshift-dependent in their own unique way.

While we have discussed only the TNG here, this statement on the mass- and
redshift-dependence of the impact of galaxy formation on halos is also true in
SIMBA. We show and briefly discuss those results in Appendix B.

###  3.5 Comparison of redshift evolution in IllustrisTNG and DMO

Figure 5: The Kllr slopes of halo concentration $c_{\rm vir}$ with halo mass
$M_{\rm vir}$, $\Omega_{\rm m}$ and $\sigma_{8}$ at $0<z<1$. The format is the
same as Figure 4. We compare the TNG slopes (solid) with the slopes from the
DMO simulations (dashed) and predictions from the model of Ishiyama et al.
(2021) using the Colossus package (Diemer, 2018). The mass-dependence of the
slopes in the DMO model and Colossis have similar trends with mass, with the
slopes decreasing towards larger mass-scales. The deviations between the TNG
and DMO runs, however, are qualitatively different for $M_{\rm vir}$ and
$\Omega_{\rm m}$, and less so for $\sigma_{8}$.

Finally, we discuss the differences between the slopes in TNG and those in the
DMO simulations of Camels. Figure 5 shows the slopes for TNG and DMO, in
addition to theoretical predictions using the model of Ishiyama et al. (2021)
as implemented in Colossus (Diemer, 2018).444We compute the slopes for the
Ishiyama et al. (2021) model by taking numerical derivatives, $\Delta\ln
c_{\rm vir}/\Delta\ln p=(\ln c_{\rm vir}(p^{+},...)-\ln c_{\rm
vir}(p^{-},...))/\Delta\ln p$, where $p$ is one of the parameters of interest.
For slopes of $c_{\rm vir}$ with $M_{\rm vir}$ and $\Omega_{\rm m}$, the TNG
results shows more complexity than either the DMO or Ishiyama model. The
slopes with $M_{\rm vir}$ have a minimum near $10^{12}\,{\rm M}_{\odot}/h$ and
slopes with $\Omega_{\rm m}$ have a maxima near $\rm 10^{13.5}\,{\rm
M}_{\odot}/h$. Both features are not presented in the DMO or Ishiyama models,
which always show monotonic trends for the dependence of the slope on halo
mass. For the slopes of $\sigma_{8}$, the TNG model has qualitatively similar
trends as the DMO and Ishiyama model in that it uniformly decreases with mass.
In each of the three slopes (panels), the redshift dependence is similar among
the three models shown. At higher redshift, the slopes of $M_{\rm vir}$
increase and those of $\Omega_{\rm m}$ and $\sigma_{\rm 8}$ decrease.

Our findings confirm that the slopes with mass and cosmological parameters is
to zeroth order given by the DMO behavior. However, there are notable
deviations — with the deviations’ amplitudes varying across $M_{\rm vir}$,
$\Omega_{\rm m}$ and $\sigma_{8}$ — for the slopes of $c_{\rm vir}$ with all
three quantities and these are caused by galaxy formation processes.

##  4 Results II: Model applications & Parameter Constraints

Having now discussed the physics extracted from examining our Kllr-based
model, we now move towards applications of the model. We first start by
describing some potential applications in §4.1. Then we show two specific
applications — a simulation “translation” to match TNG and SIMBA in §4.2 and
an observational constraint in §4.3.

###  4.1 Model features and applications

Figure 6: The $z=0$ TNG Kllr scaling relations for $c_{\rm vir}$ with various
simulation parameters (shown in the legends) and the color tones show these
relations as we vary the parameters. The solid lines are the prediction from
our model, made using the 1000 LH simulations, and the dashed lines are the
Kllr result from the 1P simulations. We find that the model generally picks up
the same trends as those found in the simulation results. Note that we do not
expect precise agreement between the model and the 1P simulations here given
the 1P simulations all share the exact same initial conditions and so cosmic
variance (which is particularly impactful here given small boxsize) can shift
a single simulation’s result high/low. The relative difference between
different colors of the same linestyle is similar across the two linsetyles.

By extracting the slopes of the halo concentration with mass and the six
simulation parameters, we have built a non-linear (but locally linear) model
for the $c_{\rm vir}-M_{\rm vir}$ relation as a function of astrophysical and
cosmological parameters. While Kllr is a non-parameteric non-linear model
(i.e. it apriori knows nothing about the scaling relations) its connection to
linear regression, through the common language of intercepts and slopes,
allows for an easy and familiar physical interpretation of the model results.
In Figure 6 we show that this model correctly captures the dependence of the
nonlinear $c_{\rm vir}-M_{\rm vir}$ relation on cosmological and astrophysical
parameters, as determined by the TNG galaxy formation model. In all panels,
the mass-dependent evolution of $c_{\rm vir}$ as a parameter is varied is
similar across our model (solid lines) and the 1P simulations (dashed lines).
There are some offsets, but we are primarily focused on the relative trends
between different colors of a given linestyle rather than offsets between
linestyles. This is because the 1P simulations are a single run with the same
initial conditions rather than the “true mean” obtained by averaging over many
realizations.

As discussed before, our model consists of the intercepts and slopes needed
for equation (2.3) and this is done for 34 redshifts and for both galaxy
formation models. These outputs can be used to construct an interpolator that
efficiently make predictions at the exact mass and redshift needed for any
analysis; we have performed such an exercise for in the above sections and
will further make use of this tool in the sections below.

In many aspects of cosmology and astrophysics, we use the $c_{\rm vir}-M_{\rm
vir}$ relation to forward model different observables. Some well-known
examples on the structure formation side are the small-scale matter power
spectra (e.g., Zacharegkas et al., 2021), and cross-correlations between
various fields (e.g., Pandey et al., 2019, 2021; Gatti et al., 2021). It is
also used extensively in analytic and numerical studies of satellite halo
dynamics (e.g., van den Bosch et al., 2018; Shen et al., 2022). All of these
works use relations calibrated from DMO simulation suites (e.g., Bhattacharya
et al., 2013; Diemer & Kravtsov, 2015; Diemer & Joyce, 2019; Ishiyama et al.,
2021). The relations calibrated in our work, however, extend these models in a
complementary way to include baryon imprints from galaxy formation processes.
In doing so, we provide a simple, yet powerful, method for incorporating
astrophysical effects into any halo model-based analysis. Of course, one could
directly measure the relevant observable, such as the matter power spectrum or
the gas-matter cross power spectrum, in Camels and build a calibration
framework for that observable specifically. However, calibrating the halo
model (i.e. the halo profiles) is inherently more intuitive given baryonic
effects occur within halos, and it is also more flexible in that it allows any
halo model-based prediction for a matter density field observable to
incorporate these baryonic effects.

In general, our model captures the simulation results to the $10\%$ accuracy
level. This is estimated by taking each LH sim and computing the residuals
between the sim’s Kllr fit — which is taken as the ground truth — and our
model’s prediction for that LH sim. We then take the residuals from the
ensemble of 1000 available LH simulations and compute the 68% interval for
these residuals as a function of mass, which leads us to the estimate of 10%
accuracy. This number is best representative of the accuracy at low masses
($M_{\rm vir}<10^{13}\,{\rm M}_{\odot}/h$). At higher masses, we cannot
reliably compute the accuracy as the Kllr fit used a ground truth has a cosmic
variance component which dominates the model accuracy and thus would dominate
the residuals. We have also tested the accuracy improvement gained by
including higher, quadratic terms, such as $A_{\rm X}^{2}$ and $A_{\rm
X}A_{\rm Y}$, and found improvements of only 1-2%. Thus, we optimize for
simplicity in the model and do not include these higher-order terms in the
local linear regression described in equation (2.3).

One subtle limitation of our model is our mass range; we have an upper bound
of $M_{\rm vir}\approx 10^{14.5}\,{\rm M}_{\odot}/h$ at $z=0$, and this
decreases with redshift. Some observables, such as the cross power spectra
between gas pressure and total matter fields, have a significant contribution
from halos above this mass range (see Figure 2 of Pandey et al., 2021, for an
example). Though, the exact mass range that is relevant depends on the
redshift/length-scale being studied. Other observables such as the lensing
power spectra or total matter power spectra have dominant contributions from
lower mass halos that are well-represented in our sample and have only
diminished contributions from halos at the upper edge of our sample’s mass
range. The lower bound of $M_{\rm vir}=10^{11}\,{\rm M}_{\odot}/h$ in this
work, set by resolution constraints, is adequate for most astrophysical and
all large-scale structure modeling. Our redshift range, with sufficient
statistics up to $z=3$ and a strict cut at $z=6$, is adequate for almost all
studies of late-time structure formation.

###  4.2 Simulation translation

Figure 7: Results from doing a “translation” between simulations using the
Kllr model. Left: The fiducial $c_{\rm vir}-M_{\rm vir}$ relation for each
model A (denoted as “Data”) and the accompanying fit using the alternative
model B (denoated as “Fit”). Results are obtained from running an MCMC. The
uncertainties shown in the “Data” come from the boostrap realizations of the
interpolated Kllr predictions, whereas for the fit they come from the MCMC
posteriors. Right: The posterior distributions for the 4 astrophysical
parameters of Camels in each galaxy formation model. Note that the meaning of
the parameters is different between TNG and SIMBA. These plots show that our
model complexity retains enough flexibility to have the results of one galaxy
formation model match that of another for a large part of the mass range.

An interesting application of our model is to take the predictions from galaxy
formation model A, obtained using the fiducial parameter values, and see how
the astrophysical parameters in a model B must deviate from their fiducial
setting so the prediction of model B matches that of model A. This is in
effect probing an effective “translation” between model B and model A.

In Figure 7 we show results from this procedure. In specific, we run a
Bayesian analysis, using Markov chain Monte Carlo (MCMC), to determine the
astrophysical parameters that make TNG and SIMBA match each other. We fix the
cosmology parameters to $\Omega_{\rm m}=0.3$ and $\sigma_{8}=0.8$, which are
the fiducial values. The astrophysical parameters are given a prior of
$0.25<\\{A,B\\}_{\rm SN1},\\{A,B\\}_{\rm AGN1}<4$ and $0.5<\\{A,B\\}_{\rm
SN2},\\{A,B\\}_{\rm AGN2}<2$, which come from the range of values used in the
LH sims in Camels.

The MCMC is run using the emcee package (Foreman-Mackey et al., 2013). For the
covariance matrix of the “data” (which is the prediction of model A), we use a
modified matrix. In general, one would get the covariance matrix using the 100
bootstraps, and that is perfectly valid. Note that the variance on the mean
$c_{\rm vir}$ increases significantly as a function of $M_{\rm vir}$ given the
halo counts drops exponentially with mass. This means the MCMC is incentivized
to fit the $c_{\rm vir}-M_{\rm vir}$ relation at the smallest mass-scales,
where even moderate improvements can greatly increase the chi-squared.

Our goal, on the other hand, is to create a toy scenario where we first assume
the mean prediction of model A is the “true mean” (i.e. we know the right mean
$c_{\rm vir}$ at every $M_{\rm vir}$ with no measurement uncertainty) and then
ask how model B can replicate that “true mean”. In this context, it is
undesirable to have the MCMC fit just the lowest masses. Thus, we “normalize”
the covariance matrix so the concentration at all masses are prioritized
equally in the fit. We first start by computing the correlation matrix
$r_{ij}$ and then convert it to a covariance matrix assuming a constant
scatter of $1\%$, in the appropriate $\rm dex$ units. By forcing the variance
of the measurements to be mass-independent, we make the MCMC procedure fit all
masses equally.555This can also be viewed as regularizing the optimization
function, which is simple least squares regression, to include some mass-
dependence where measurements at lower masses are downweighted to offset their
high precision relative to those at higher masses. Note that with this
procedure we still keep the correlation structure of the covariance matrix,
i.e. we do not ignore any off-diagonal terms. We have also checked that
varying the value of the $1\%$ constant scatter does not affect our final
results

We perform this analysis in two settings, each time swapping which of TNG and
SIMBA serve as model A and model B. Both results are shown in Figure 7.
Qualitatively we find that for TNG to match SIMBA, we must reduce the SN
energy scale ($A_{\rm SN1}$) and reduce the AGN jet speed and burstiness
($A_{\rm AGN2}$). The other parameters are within $<2\sigma$ of their fiducial
values. To match SIMBA with TNG, we need to reduce the mass loading factor
($B_{\rm SN1}$) and increase the supernova wind speed ($B_{\rm SN2}$). The AGN
parameters are also within $<2\sigma$ of their fiducial values

In general the preference of SN parameters in both cases is expected. A look
at the dashed lines on left panels of Figure 7, which are the fiducial
relations from TNG and SIMBA, shows that TNG has a distinct overturn at low
masses which is not present in SIMBA. In TNG this feature is caused by the
interplay between cooling increasing concentration at MW-scales and SN
feedback reducing it at lower masses (A22). A similar feature is not present
at the same mass-scale in SIMBA, likely due to a weaker SN feedback or more
dominant cooling, which will keep increasing the concentration without any
overturn. Notably, the mass loading factor in SIMBA is roughly an order of
magnitude higher at the smallest masses when compared to the corresponding
values in TNG (Davé et al., 2019, see their Figure 1), and more mass loading
increases the metal enrichment of the halo gas, leading to more cooling and
efficient adiabatic contraction. So to match the predictions of the two
simulations, one must either modify TNG to reduce SN power (reduce $A_{\rm
SN1}$), or change SIMBA to increase SN power (increase $B_{\rm SN2}$) and/or
reduce gas cooling (decrease $B_{\rm SN1}$).

For TNG to match SIMBA we see that $A_{\rm AGN2}$ must take values below the
fiducial setting of $A_{\rm AGN2}=1$. Given our current results, we are unable
to form a more informative picture beyond realizing that TNG prefers to have
weaker but more frequent kinetic jet feedback events. It is also worth noting
that simulations run with versions of the SIMBA model modified for galaxy
cluster simulations (Cui et al., 2022) have been found to heat the gas
considerably more compared to the models of TNG and other simulations (Lee et
al., 2022). This suggests the AGN feedback in SIMBA has a stronger impact on
its environment than in TNG. This is also supported by studies of the baryon
radial distribution in TNG and in the original, unmodified SIMBA model
(Ayromlou, Nelson & Pillepich, 2022).

In general, a key takeaway is we are able to make one model match the other in
the mass range of $M_{\rm vir}>10^{11.5}\,{\rm M}_{\odot}/h$ or $M_{\rm
vir}>10^{12}\,{\rm M}_{\odot}/h$, depending on which model is the “data”.
Right at and below these mass ranges the difference in galaxy formation
physics is far too large to allow us to “translate” between the models.
Conversely, this means if we can obtain precise observational measurements of
the concentration at these lower masses ($M_{\rm vir}<10^{12}\,{\rm
M}_{\odot}/h$), we will be able to better determine which of the TNG and SIMBA
models more closely match the real Universe. Note that while the TNG and SIMBA
fiducial relations do have mass ranges where they are generally in good
agreement, this happens above $M_{\rm vir}>10^{12.5}\,{\rm M}_{\odot}/h$ but
still with some visible deviations (Figure 2). By varying the astrophysical
parameters input to the model, we get the agreement to improve over all mass-
scales and also to stretch down in mass by an order of magnitude.

In addition, while we are varying four highly-relevant astrophysical
parameters, there are additional processes with more nuanced effects that we
do not vary given Camels does not vary the equations. So we apriori do not
expect to be able to precisely convert results between simulations given the
modest degrees of freedom in our model. For example, A22 found that the
quasar/thermal feedback mode was the most relevant process at $M_{\rm
vir}\approx 10^{12}\,{\rm M}_{\odot}/h$ scales, and also that the SMBH mass-
scale setting the transition from thermal to kinetic feedback played a role in
shifting the nonlinear mass-dependence of $c_{\rm vir}$; both are parameters
that are not varied in the TNG version of Camels. Finally, we note that we
have only performed this “translation” exercise for $z=0$.

Note that using TNG as the model in the fitting procedure results in tighter
posteriors for the parameters (i.e. blue posteriors are tighter than green in
Figure 7), and this is simply because the slopes between $c_{\rm vir}$ and the
parameters in TNG are generally of a larger amplitude, which naturally results
in tighter constraints on said parameters.

###  4.3 Observational constraints

While we have used our Kllr-based model so far to get parameter constraints
“translating” between simulation predictions, we can also use observational
measurements on the $c_{\rm vir}-M_{\rm vir}$ relation to obtain constraints
on the simulation parameters for a given galaxy formation model. The data we
use for this analysis is described in Section 2.2 and is composed of results
from multiple works (Mandelbaum et al., 2006; Johnston et al., 2007;
Mandelbaum et al., 2008; Brimioulle et al., 2013).

We analyze this data using both the TNG model and the SIMBA model, and results
are shown in Figure 8. We find, unsurprisingly, that we are still prior
dominated for multiple parameters, where our prior is simply the range of
values considered in the LH simulations and thus the widest range of values
our model would be valid over. A large part of why we are prior dominated is
due to the sample size of available data, which is primarily from older
surveys. This prior-domination can also be a sign of a non-optimal model; the
4 parameters varied in Camels were not optimized to find processes that most
impact $c_{\rm vir}$. While most of them do have a strong impact on $c_{\rm
vir}$, we showed in Figure 3 that one of the AGN parameters had weak impact on
$c_{\rm vir}$, and this lends itself to the resulting posteriors being prior-
dominated.

Even though many of our constraints are prior-dominated, we are still able to
discuss qualitative behaviors for most parameters. Similar to our discussion
on the simulation “translation” analysis, we re-emphasize here that our
conservative claims are on whether the constraints prefer the feedback
parameters to be larger/smaller than their fiducial value. A precise estimate
of the value they must take is both not achievable given our data quality, and
will also require significant, additional robustness checks that are beyond
the scope of this work. Many of the approximations used for the observational
data in this work (see Section 2.2 for details) are inaccurate and thus will
need to be changed to enable precise estimates.

Figure 8: Top: The posteriors on the Camels astrophysical parameters from
fitting the two galaxy formation models to observational data. Bottom: The
$c_{\rm vir}-M_{\rm vir}$ relations for TNG and SIMBA models as determined by
the MCMC fits, and the data points used in the analysis. The two models
generally agree above $10^{12}\,{\rm M}_{\odot}/h$ and start disagreeing more
at lower masses where we have no data.

Figure 8 shows our observational constraints. We find $\\{A,B\\}_{\rm SN1}$
prefers lower values than the fiducial, though we do not actually see a
constraint here as the results are cutoff by the lower prior boundary. In TNG
this reduces the energy scale of the SN feedback, and in SIMBA it reduces the
mass-loading factor, which at fixed wind velocity reduces the SN energy scale
as well. We do find interesting constraints on $\\{A,B\\}_{\rm SN2}$, which
has a stronger impact on $c_{\rm vir}$ and at a slightly higher halo mass than
$\\{A,B\\}_{\rm SN1}$; TNG (SIMBA) appears to prefer values of $A_{\rm SN2}$
($B_{\rm SN2}$) that are the same (higher) than the fiducial. This is
consistent with our previous statement in Section 4.2 that the stellar
feedback in TNG, which causes the downturn in $c_{\rm vir}-M_{\rm vir}$ at
lower masses, is stronger than that in SIMBA. For the first AGN parameter, we
see that TNG (SIMBA) appears to prefer values of $A_{\rm AGN1}$ ($B_{\rm
AGN1}$) that are the lower (higher) than the fiducial value, while for the
second one (which controls jet speed) both simulations prefer higher values
than the fiducial. We also note that the $\chi^{2}$ of the results show an
equal preference for TNG and SIMBA, which is expected given that the largest
differences between the two simulation predictions occur at the lowest halo
masses, where we do not have any constraining data,

Finally, one may worry to what extent other effects on the $c_{\rm vir}-M_{\rm
vir}$ relation (either physical or observational systematics-related) can
appear as a feedback-related effect. In general, this is a very critical point
to understand, especially if one were pursue precise, robust constraints on
the parameters. Our model here can be thought of as a “fiducial” $c_{\rm
vir}-M_{\rm vir}$ relation from a hydrodynamic simulation, with additional
“basis” functions — in the form of the mass-dependent slopes — to capture the
impact of different parameters. Given these bases have a specific mass-
dependence that is taken from simulations, a non-astrophysical effect will
need to have a similar mass-dependence in order to mix their signal with an
astrophysical effect.

Now, the caveat here is that feedback processes can still have a relatively
mass-localized impact — for example, Weinberger et al. (2018) and A22 show the
strongest impact for a given process (SN feedback, thermal AGN feedback etc.)
is localized within $\sim 1\rm{\,\,\rm dex}$ in halo mass — and so it is
possible that any systematic that affects just halos at the lowest or highest
masses (where the measurements come from different samples, and thus have
different systematics) could more easily contaminate the astrophysical
constraints. Both better data quality and a dedicated analysis of any effects
would be required to disentagle all these factors. The purpose of this section
is to demonstrate and place some observational constraints and then
qualitatively interpret their meaning (i.e. are the parameter values
higher/lower than the chosen fiducial values in TNG and SIMBA), so a detailed
exploration is beyond the scope of our work.

A simple forecast – where we take the mean $c_{\rm vir}-M_{\rm vir}$
prediction from our TNG model with fiducial parameter values as the simulated
datavector, add a 5% measurement uncertainty, and run the MCMC — shows that we
can get constraints of around $\sigma(A_{X})\in\\{0.25,1.0,0.1,0.3\\}$ for the
four astrophysical parameters. The median uncertainty of the existing
measurements in Figure 8 is $\sim 20\%$.

##  5 Conclusion

The halo concentration–mass relation is a key quantity for understanding the
formation and evolution of halos, and is actively being used in modelling many
astrophysical and cosmological observables. Existing works have provided semi-
analytic models for this relation as a function of cosmology using high-
resolution DMO simulations. However, the quantitative impact of astrophysics
on this relation is relatively unexplored, though many qualitative analyses
already exist.

In this work, we use the Camels simulation suite to study the mass- and
redshift-dependent response of halo concentration from galaxy formation
processes, in particular SN and AGN feedback. We do so by building a locally
linear but globally non-linear model of $c_{\rm vir}$ as a function of
redshift, mass, and cosmology/astrophysics parameters for two different galaxy
formation models. Our choice allows for easy interpretability of our results,
and the identification of physical processes leading to specific features,
while capturing most of the nonlinearity in the relations. By using an uber
halo sample constructed from 1000 simulations, we achieve higher statistical
power and thereby analyze higher mass-scales that are normally unachievable in
the small $L=25\,h^{-1}\,{\rm Mpc}$ boxes of the individual Camels
simulations.

Our main results are summarized as follows:

  * •

Figure 3 shows the dependence of $c_{\rm vir}$ on various quantities. The SN
feedback in TNG is most impactful at the $10^{11}-10^{12}\,{\rm M}_{\odot}/h$
mass range, given the shallower halo potentials at lower masses, and works to
make the halo more diffuse. There continues to be some minor impact at larger
mass scales, where gas cooling is known to be comparable/larger than stellar
feedback. This is likely caused by changes to the mass loading factor —
controlled by both SN parameters — which impacts the gas metal enrichment and
thus metal cooling, which alters how efficiently the halo potential can
contract.

  * •

AGN feedback in TNG becomes prominent at $M_{\rm vir}>10^{13}\,{\rm
M}_{\odot}/h$, which is the mass at which the energy scale of kinetic feedback
(radio jets) starts dominating the halo (A22).

  * •

Figure 3 also shows that the response of $c_{\rm vir}$ to SN and AGN feedback
is quantitatively different in SIMBA compared to the response in TNG. This is
expected as the equations being modified by the Camels astrophysics parameters
are different. Notably, the SN mechanism is significantly weaker in SIMBA
(reflected in how the slopes are closer to 0).

  * •

The dependence of $c_{\rm vir}$ on cosmology in the two simulations is
generally consistent but differs in detail, especially at low masses. This
highlights how the cosmology dependence, when conditioned on a given galaxy
formation model, can have notable variations. The dependence on $\sigma_{8}$
has the same mass-dependence in all three of TNG, DMO, and semi-analytical
predictions (Figure 5). The dependence on $\Omega_{\rm m}$ and $M_{\rm vir}$
is also similar but has notable non-monotonic deviation.

  * •

The redshift dependence of the qualitative statements above is generally weak
(Figure 4). Some interesting features are that the slope between $c_{\rm
vir}-A_{\rm SN1}$ shows a zero-crossing at low redshift, which then disappears
at higher redshift. The slopes for $c_{\rm vir}-A_{\rm SN2}$ also has a
stationary point around $M_{\rm vir}=5\times 10^{12}\,{\rm M}_{\odot}/h$, with
strong redshift dependence below this mass-scale and weak dependence above it.

  * •

Using our model to “translate” between simulations we can check how the Camels
astrophysics parameters must be changed in order to get galaxy formation model
A to match the fiducial result of model B. We find either that TNG must have
lower SN feedback energy output, or that SIMBA must have lower mass loading
and/or higher SN wind speeds. There is also a mild preference in TNG to have
less bursty AGN kinetic jets.

  * •

Finally, we do a simplistic test of our model on observational data for the
$c_{\rm vir}-M_{\rm vir}$ relation, and find that in most cases the
astrophysical parameters prefer values away from the fiducial setting of
$A_{X}=B_{X}=1$ (Figure 8). The parameters from the TNG- and SIMBA-model
analysis imply a preference in both for having less energetic SN feedback and
stronger AGN kinetic jet speeds.

Through performing a qualitative and quantitative study of $c_{\rm vir}$ in
simulations, we extract the specific impact of galaxy formation on the matter
distribution, identifying the mass- and redshift-dependent impact of specific
processes. From a halo formation standpoint, this allows us to closely examine
the exact impacts on the halo density distribution due to these processes.
From an application standpoint, it allows us to understand the signature of
specific processes on traditional, cosmological observables that can be
forward-modelled using the halo model approach (eg. nonlinear power spectra).

Focusing on cosmological applications in particular, a highly relevant
extension of this work would be to analyze a suite of simulations of a larger
box size (while accommodating coarser mass resolution). Our main limitation in
this work is the exact upper mass limit we can use in our scaling relations,
and the cosmic variance component at these higher masses. While the current
setup still allows for meaningful studies of fields like the total matter or
lensing fields, which are most sensitive to halo masses well within our range,
other studies of structure may require model predictions for more massive
halos.

An ideal scenario of this type of work would also be to recover the DMO result
by “turning off” the galaxy formation processes in a model for $c_{\rm
vir}-M_{\rm vir}$. However, this is possible only if we can turn off all
baryonic processes, which would require extending Camels to include on the
order of $\mathcal{O}(100)$ parameters, rather than the 4 being currently
considered. This is computationally infeasible at the moment. However,
increasing the number of parameters by reasonable levels will continue to give
us greater insight into the physics of galaxy formation and its impact on the
larger structures.

The matter distribution within halos is an interesting problem. It has a rich
phenomenology involving processes from gravitational collapse to gas cooling
and feedback to metal enrichment to baryon backreactions and so forth.
Understanding these processes is on its own interesting from a galaxy
formation and galaxy evolution perspective, but it is then also a necessary
component in understanding how the observed structure in our Universe behaves
on nonlinear scales. Existing work has already extracted the qualitative
features of baryon imprints. This work takes the first step in performing a
more quantitative analysis of the same question, by leveraging a large number
of simulations and answering how those baryon imprints vary quantitatively as
specific galaxy formation processes are altered. As more observational
constraints are obtained and more galaxy formation models are refined (through
quantitative predictive models such as the ones explored here), we will be
able to narrow the possible galaxy formation processes that take place in our
Universe and have a more definitive picture of how baryons impact the matter
distributions.

## Acknowledgements

We thank Andrew Hearin for many helpful comments on an earlier version of this
draft. We also thank the Camels collaboration for making their data publicly
available and easily accessible to the community. We thank the referee for
their helpful comments.

DA is supported by the National Science Foundation Graduate Research
Fellowship under Grant No. DGE 1746045. CC is supported by the Henry Luce
Foundation and DOE grant DE-SC0021949.

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/CAMELSconcentration. The
public Camels data, including the Rockstar catalogs we use in this work, can
be found at https://camels.readthedocs.io.

## References

  * Abadi et al. (2010) Abadi M. G., Navarro J. F., Fardal M., Babul A., Steinmetz M., 2010, MNRAS, 407, 435
  * Abbott et al. (2022) Abbott T. M. C., et al., 2022, Phys. Rev. D, 105, 023520
  * Allen et al. (2011) Allen S. W., Evrard A. E., Mantz A. B., 2011, ARA&A, 49, 409
  * Amon et al. (2022) Amon A., et al., 2022, Phys. Rev. D, 105, 023514
  * 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
  * Artale et al. (2019) Artale M. C., Pedrosa S. E., Tissera P. B., Cataldi P., Di Cintio A., 2019, A&A, 622, A197
  * Asgari et al. (2021) Asgari M., et al., 2021, A&A, 645, A104
  * Ayromlou et al. (2022) Ayromlou M., Nelson D., Pillepich A., 2022, arXiv e-prints, p. arXiv:2211.07659
  * Behroozi et al. (2013a) Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013a, ApJ, 762, 109
  * Behroozi et al. (2013b) Behroozi P. S., Wechsler R. H., Conroy C., 2013b, ApJ, 770, 57
  * Beltz-Mohrmann & Berlind (2021) Beltz-Mohrmann G. D., Berlind A. A., 2021, arXiv e-prints, p. arXiv:2103.05076
  * Bhattacharya et al. (2013) Bhattacharya S., Habib S., Heitmann K., Vikhlinin A., 2013, ApJ, 766, 32
  * Blumenthal et al. (1986) Blumenthal G. R., Faber S. M., Flores R., Primack J. R., 1986, ApJ, 301, 27
  * Brimioulle et al. (2013) Brimioulle F., Seitz S., Lerchster M., Bender R., Snigula J., 2013, MNRAS, 432, 1046
  * 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
  * Child et al. (2018) Child H. L., Habib S., Heitmann K., Frontiere N., Finkel H., Pope A., Morozov V., 2018, ApJ, 859, 55
  * Chisari et al. (2018) Chisari N. E., et al., 2018, MNRAS, 480, 3962
  * Chua et al. (2019) Chua K. T. E., Pillepich A., Vogelsberger M., Hernquist L., 2019, MNRAS, 484, 476
  * Chua et al. (2021) Chua K. T. E., Vogelsberger M., Pillepich A., Hernquist L., 2021, arXiv e-prints, p. arXiv:2109.00012
  * Cooray & Sheth (2002) Cooray A., Sheth R., 2002, Phys. Rep., 372, 1
  * Cui et al. (2016) Cui W., et al., 2016, MNRAS, 458, 4052
  * Cui et al. (2018) Cui W., et al., 2018, MNRAS, 480, 2898
  * Cui et al. (2022) Cui W., et al., 2022, MNRAS, 514, 977
  * Davé et al. (2019) Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827
  * Diemer (2018) Diemer B., 2018, ApJS, 239, 35
  * Diemer & Joyce (2019) Diemer B., Joyce M., 2019, ApJ, 871, 168
  * Diemer & Kravtsov (2015) Diemer B., Kravtsov A. V., 2015, ApJ, 799, 108
  * 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
  * 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
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Forouhar Moreno et al. (2021) Forouhar Moreno V. J., Benítez-Llambay A., Cole S., Frenk C., 2021, arXiv e-prints, p. arXiv:2107.14245
  * Gatti et al. (2021) Gatti M., et al., 2021, arXiv e-prints, p. arXiv:2108.01600
  * Gnedin et al. (2004) Gnedin O. Y., Kravtsov A. V., Klypin A. A., Nagai D., 2004, ApJ, 616, 16
  * Heitmann et al. (2009) Heitmann K., Higdon D., White M., Habib S., Williams B. J., Lawrence E., Wagner C., 2009, ApJ, 705, 156
  * Hirschmann et al. (2014) Hirschmann M., Dolag K., Saro A., Bachmann L., Borgani S., Burkert A., 2014, MNRAS, 442, 2304
  * Hu & Kravtsov (2003) Hu W., Kravtsov A. V., 2003, ApJ, 584, 702
  * Hunter (2007) Hunter J. D., 2007, Computing in Science and Engineering, 9, 90
  * Ishiyama et al. (2021) Ishiyama T., et al., 2021, MNRAS, 506, 4210
  * Johnston et al. (2007) Johnston D. E., et al., 2007, arXiv e-prints, p. arXiv:0709.1159
  * Kravtsov & Borgani (2012) Kravtsov A. V., Borgani S., 2012, ARA&A, 50, 353
  * Lee et al. (2022) Lee E., et al., 2022, MNRAS, 517, 5303
  * 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
  * Mandelbaum et al. (2006) Mandelbaum R., Seljak U., Cool R. J., Blanton M., Hirata C. M., Brinkmann J., 2006, MNRAS, 372, 758
  * Mandelbaum et al. (2008) Mandelbaum R., Seljak U., Hirata C. M., 2008, J. Cosmology Astropart. Phys., 2008, 006
  * Marinacci et al. (2018) Marinacci F., et al., 2018, MNRAS, 480, 5113
  * McCarthy et al. (2017) McCarthy I. G., Schaye J., Bird S., Le Brun A. M. C., 2017, MNRAS, 465, 2936
  * McKinney (2011) McKinney W., 2011, Python for High Performance and Scientific Computing, 14
  * 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
  * Ntampaka et al. (2019) Ntampaka M., et al., 2019, ApJ, 876, 82
  * Pandey et al. (2019) Pandey S., et al., 2019, Phys. Rev. D, 100, 063519
  * Pandey et al. (2021) Pandey S., et al., 2021, arXiv e-prints, p. arXiv:2108.01601
  * Pedrosa et al. (2010) Pedrosa S., Tissera P. B., Scannapieco C., 2010, MNRAS, 402, 776
  * 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
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Ragagnin et al. (2019) Ragagnin A., Dolag K., Moscardini L., Biviano A., D’Onofrio M., 2019, MNRAS, 486, 4001
  * Rudd et al. (2008) Rudd D. H., Zentner A. R., Kravtsov A. V., 2008, ApJ, 672, 19
  * Sánchez et al. (2022) Sánchez J., et al., 2022, arXiv e-prints, p. arXiv:2210.08633
  * 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
  * Schneider et al. (2022) Schneider A., Giri S. K., Amodeo S., Refregier A., 2022, MNRAS, 514, 3802
  * Secco et al. (2022) Secco L. F., et al., 2022, Phys. Rev. D, 105, 023515
  * Semboloni et al. (2011) Semboloni E., Hoekstra H., Schaye J., van Daalen M. P., McCarthy I. G., 2011, MNRAS, 417, 2020
  * Seppi et al. (2021) Seppi R., et al., 2021, A&A, 652, A155
  * Shen et al. (2022) Shen X., Xiao H., Hopkins P. F., Zurek K. M., 2022, arXiv e-prints, p. arXiv:2207.11276
  * Shin et al. (2021) Shin T., et al., 2021, arXiv e-prints, p. arXiv:2105.05914
  * 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
  * Strickland & Heckman (2009) Strickland D. K., Heckman T. M., 2009, ApJ, 697, 2030
  * Tissera et al. (2010) Tissera P. B., White S. D. M., Pedrosa S., Scannapieco C., 2010, MNRAS, 406, 922
  * Tröster et al. (2022) Tröster T., et al., 2022, A&A, 660, A27
  * Van der Walt et al. (2011) Van der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science and Engineering, 13, 22
  * Vikram et al. (2017) Vikram V., Lidz A., Jain B., 2017, MNRAS, 467, 2315
  * 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
  * Weinberg et al. (2013) Weinberg D. H., Mortonson M. J., Eisenstein D. J., Hirata C., Riess A. G., Rozo E., 2013, Phys. Rep., 530, 87
  * Weinberger et al. (2017) Weinberger R., et al., 2017, MNRAS, 465, 3291
  * Weinberger et al. (2018) Weinberger R., et al., 2018, MNRAS, 479, 4056
  * Wetzel et al. (2016) Wetzel A. R., Hopkins P. F., Kim J.-h., Faucher-Giguère C.-A., Kereš D., Quataert E., 2016, ApJ, 827, L23
  * Wheeler et al. (2019) Wheeler C., et al., 2019, MNRAS, 490, 4447
  * York et al. (2000) York D. G., et al., 2000, AJ, 120, 1579
  * Zacharegkas et al. (2021) Zacharegkas G., et al., 2021, arXiv e-prints, p. arXiv:2106.08438
  * van den Bosch et al. (2018) van den Bosch F. C., Ogiya G., Hahn O., Burkert A., 2018, MNRAS, 474, 3043

##  Appendix A Outlier rejection

The concentration, $c_{\rm vir}$, for halos $M_{\rm vir}>10^{11}\,{\rm
M}_{\odot}/h$ at a given redshift has a notable fraction of outliers with
$c_{\rm vir}\lesssim 1$. For these objects the scale radius is equivalent to
$R_{\rm vir}$, and this is not a physical population of halos. We do not
explore the exact reasons underlying the existence of this population, as it
is beyond the scope of our work. However, we must remove such outliers before
computing our scaling relations as they will bias the mean concentration to
lower values. Thus, we require an adequate rejection algorithm.

To reject outliers within a given halo sample, we first create a histogram of
the halo concentration in 200 bins. A gaussian filter with a smoothing scale
equal to the bin width then smooths the histogram to remove any noise
fluctuations. We interpolate the function via a cubic spline. The interpolated
distributions for $z\in\\{0,1\\}$ for TNG are shown in Figure 9, with a clear
peak at low $c_{\rm vir}$. We use the inflection point at low $c_{\rm vir}$ as
our cutoff; all halos with $c_{\rm vir}$ below this point are thrown out of
the analysis. This inflection point is obtained by getting the first root of
the derivative of the interpolated histogram. For $z=0$ ($z=1$), we find that
$\approx 2\%$ ($\approx 6\%$) of halos are rejected using this procedure.

Figure 9: The distribution of $c_{\rm vir}$ at $z=0$ and $z=1$, for halos with
$M_{\rm vir}>10^{11}\,{\rm M}_{\odot}/h$. We see the distributions are
generally well-behaved but have a sharp increase at low values of $c_{\rm
vir}\lesssim 1$. The vertical dashed lines indicate the inflection point of
the distribution, which is where the histogram rises to the left to form the
low $c_{\rm vir}$ peak. All halos with $c_{\rm vir}$ below this threshold are
rejected.

##  Appendix B Redshift dependence in SIMBA

Figure 10: The Kllr slopes between halo concentration $c_{\rm vir}$ and the
six astrophysical and cosmological parameters. Same format as Figure 4. The
bands show the 68% confidence interval from the model uncertainty. Similar to
the case in IllustrisTNG, there is a clear and strong mass-dependence in all
parameters, and a significant redshift evolution in most. Figure 11: The Kllr
slopes of halo concentration $c_{\rm vir}$ with halo mass $M_{\rm vir}$,
$\Omega_{\rm m}$ and $\sigma_{8}$ at $0<z<1$. The format is the same as Figure
5. We compare the SIMBA slopes (solid) with the slopes from the DMO
simulations (dashed) and predictions from the model of Ishiyama et al. (2021)
using the Colossus package (Diemer, 2018).

The redshift dependence discussed in this work has so far been about the
IllustrisTNG model (Figure 4); a choice made for brevity sake. In this
appendix, we extract the redshift dependence in Simba and discuss it. Once
again, we find that the slopes depend on mass in a non-monotonic way across
all redshifts.

In $\boldsymbol{B_{\rm SN1}}$, the slopes are generally low in amplitude
across all masses. There is a local maxima around $10^{12}\,{\rm M}_{\odot}/h$
and the mass-scale of this maxima increases with redshift. The fact that the
slopes are overall positive is unsurprising — $B_{\rm SN1}$ controls the wind
mass-loading factor which, as we have discussed previously, impacts the gas
cooling efficiency through metal enrichment. Increases in this efficiency
increase $c_{\rm vir}$, leading to the generally positive slope, especially at
lower redshifts where cooling has had more time to impact the matter
distribution.

In $\boldsymbol{B_{\rm SN2}}$, we find behavior similar to $A_{\rm SN2}$;
there is a minima in the slopes at $10^{11.5}\,{\rm M}_{\odot}/h$ that
increases with redshift. This is also unsurprising as both $A_{\rm SN1}$ and
$B_{\rm SN2}$ control the SN wind speeds of their respective galaxy formation
models, so even though the wind speeds are implemented differently in either
model (i.e. different equations) both parameters change the same underlying
quantity.

In $\boldsymbol{B_{\rm AGN1}}$ and $\boldsymbol{B_{\rm AGN2}}$, the slopes
deviate more from 0 compared to the TNG counterparts. However, the amplitude
of the slope is still quite weak. $B_{\rm AGN2}$ controls the jet speed, which
is qualitatively similar to the “burstiness” controlled by $A_{\rm AGN2}$ and
so we see similar behavior of slopes being zero at low mass but becoming non-
zero at high mass. In fact the AGN slopes in Simba have a higher amplitude
than their TNG counterpart.

Finally, $\boldsymbol{\Omega_{\rm m}}$ and $\boldsymbol{\sigma_{8}}$ have
similar trends in Simba as they did in IllustrisTNG. The dependence on
$\sigma_{8}$ is quite similar, showing lower variation with redshift at low
masses and stronger variations at high mass. The slopes in Simba seem to keep
rising towards the lower masses, especially at low redshift, whereas in TNG
they asymptoted to a flat, mass-independent trend. The redshift dependence of
the $\Omega_{\rm m}$ slopes differs from that of TNG. At high masses there
appears to be a constant multiplicative factor accounting for the redshift
dependence. At low masses, the mass-dependence of the slope steepes more
strongly at low redshift, and flattens at high redshift. Even though these are
cosmological parameters that are detached from galaxy formation, the actual
impact of cosmology on structure formation is some complicated convolution of
cosmological evolution with galaxy formation. Given Simba and IllustrisTNG
have different galaxy formation models, the dependence of $c_{\rm vir}$ on
cosmology will vary as well.

Lastly, for completeness, we also show in Figure 11 the redshift-dependent
comparison of slopes from SIMBA and DMO runs, similar to Figure 5. The
behavior in these plots follows from the discussions above. In general, there
continue to be notable differences between the DMO and hydrodynamic runs, as
is expected.