# The AGORA High-resolution Galaxy Simulations Comparison Project. V:
Satellite Galaxy Populations In A Cosmological Zoom-in Simulation of A Milky
Way-mass Halo

Minyong Jung<EMAIL_ADDRESS>Center for Theoretical Physics, Department of
Physics and Astronomy, Seoul National University, Seoul 08826, Korea  Santi
Roca-Fàbrega Code leaders<EMAIL_ADDRESS>Lund Observatory,
Department of Astronomy and Theoretical Physics, Lund University, Box 43,
SE-221 00 Lund, Sweden  Departamento de Física de la Tierra y Astrofísica,
Facultad de Ciencias Físicas, Plaza Ciencias, 1, 28040 Madrid, Spain  Ji-hoon
Kim Code leaders<EMAIL_ADDRESS>Center for Theoretical Physics, Department
of Physics and Astronomy, Seoul National University, Seoul 08826, Korea  Seoul
National University Astronomy Research Center, Seoul 08826, Korea  Anna Genina
Code leaders  Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1,
D-85748, Garching, Germany  Loic Hausammann  Code leaders  ITS High
Performance Computing, Eidgenössische Technische Hochschule Zürich, 8092
Zürich, Switzerland  Institute of Physics, Laboratoire d’Astrophysique, École
Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland  Hyeonyong
Kim Code leaders  Center for Theoretical Physics, Department of Physics and
Astronomy, Seoul National University, Seoul 08826, Korea  Department of
Aerospace Engineering, Seoul National University, Seoul 08826, Korea
Alessandro Lupi  Code leaders  DiSAT, Università degli Studi dell’Insubria,
via Valleggio 11, I-22100 Como, Italy  Dipartimento di Fisica “G. Occhialini”,
Università degli Studi di Milano-Bicocca, I-20126 Milano, Italy  Kentaro
Nagamine Code leaders  Department of Earth and Space Science, Graduate School
of Science, Osaka University, Toyonaka, Osaka, 560-0043, Japan  Kavli IPMU
(WPI), University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba, 277-8583, Japan
Department of Physics & Astronomy, University of Nevada Las Vegas, Las Vegas,
NV 89154, USA  Johnny W. Powell Code leaders  Department of Physics, Reed
College, Portland, OR 97202, USA  Yves Revaz  Code leaders  Institute of
Physics, Laboratoire d’Astrophysique, École Polytechnique Fédérale de
Lausanne, CH-1015 Lausanne, Switzerland  Ikkoh Shimizu  Code leaders  Shikoku
Gakuin University, 3-2-1 Bunkyocho, Zentsuji, Kagawa, 765-8505, Japan  Héctor
Velázquez  Code leaders  Instituto de Astronomía, Universidad Nacional
Autónoma de México, A.P. 70-264, 04510, Mexico, D.F., Mexico  Daniel Ceverino
Universidad Autónoma de Madrid, Ciudad Universitaria de Cantoblanco, E-28049
Madrid, Spain  CIAFF, Facultad de Ciencias, Universidad Autónoma de Madrid,
E-28049 Madrid, Spain  Joel R. Primack Department of Physics, University of
California at Santa Cruz, Santa Cruz, CA 95064, USA  Thomas R. Quinn
Department of Astronomy, University of Washington, Seattle, WA 98195, USA
Clayton Strawn Department of Physics, University of California at Santa Cruz,
Santa Cruz, CA 95064, USA  Tom Abel Kavli Institute for Particle Astrophysics
and Cosmology, Stanford University, Stanford, CA 94305, USA  Department of
Physics, Stanford University, Stanford, CA 94305, USA  SLAC National
Accelerator Laboratory, Menlo Park, CA 94025, USA  Avishai Dekel  Center for
Astrophysics and Planetary Science, Racah Institute of Physics, The Hebrew
University, Jerusalem 91904, Israel  Bili Dong  Department of Physics, Center
for Astrophysics and Space Sciences, University of California at San Diego, La
Jolla, CA 92093, USA  Boon Kiat Oh Department of Physics, University of
Connecticut, U-3046, Storrs, CT 06269, USA  Center for Theoretical Physics,
Department of Physics and Astronomy, Seoul National University, Seoul 08826,
Korea  Romain Teyssier  Department of Astrophysical Sciences, Princeton
University, Princeton, NJ 08544 USA  the AGORA Collaboration
http://www.AGORAsimulations.org The authors marked with * as code leaders
contributed to the article by leading the effort within each code group to
perform and analyze simulations.

(Received April XX, 2023; Revised April XX, 2023; Accepted May XX, 2023)

###### Abstract

We analyze and compare the satellite halo populations at $z\sim 2$ in the
high-resolution cosmological zoom-in simulations of a $10^{12}\,{\rm
M}_{\odot}$ target halo ($z=0$ mass) carried out on eight widely-used
astrophysical simulation codes (Art-I, Enzo, Ramses, Changa, Gadget-3, Gear,
Arepo-t, and Gizmo) for the AGORA High-resolution Galaxy Simulations
Comparison Project. We use slightly different redshift epochs near $z=2$ for
each code (hereafter “$z\sim 2$’) at which the eight simulations are in the
same stage in the target halo’s merger history. After identifying the matched
pairs of halos between the CosmoRun simulations and the DMO simulations, we
discover that each CosmoRun halo tends to be less massive than its DMO
counterpart. When we consider only the halos containing stellar particles at
$z\sim 2$, the number of satellite galaxies is significantly fewer than that
of dark matter halos in all participating AGORA simulations, and is comparable
to the number of present-day satellites near the Milky Way or M31. The so-
called “missing satellite problem’ is fully resolved across all participating
codes simply by implementing the common baryonic physics adopted in AGORA and
the stellar feedback prescription commonly used in each code, with sufficient
numerical resolution ($\lesssim 100$ proper pc at $z=2$). We also compare
other properties such as the stellar mass$-$halo mass relation and the
mass$-$metallicity relation. Our work highlights the value of comparison
studies such as AGORA, where outstanding problems in galaxy formation theory
are studied simultaneously on multiple numerical platforms.

cosmology: theory – galaxies: formation – galaxies: evolution – galaxies:
kinematics and dynamics – galaxies: structure – galaxies: ISM – methods:
numerical – hydrodynamics

##  1 Introduction

Studied extensively by cosmologists, the $\Lambda$-Cold Dark Matter
($\Lambda$CDM) model is considered the standard model of Big Bang cosmology,
encompassing dark energy and dark matter. However, there is a certain tension
between theory and observed galaxies, especially on a small scale (for
reviews, see e.g., Bullock & Boylan-Kolchin, 2017; Del Popolo & Le Delliou,
2017). For example, the observed number of dwarf galaxies around the Local
Group is significantly fewer than that of the dark matter halos found in
$N$-body simulations when compared based on their circular velocity. This so-
called “missing satellite problem” is one of the long-standing challenges of
the contemporary $\Lambda$CDM model (Kauffmann et al., 1993; Klypin et al.,
1999; Moore et al., 1999; Benson et al., 2002). Reproducing satellite galaxies
and small-scale substructures in a simulation within the $\Lambda$CDM
framework is a nontrivial task because it requires high numerical resolution
and sophisticated baryonic physics.

The mismatch between the theory and the observation on a small scale has
motivated a great deal of theoretical modeling such as the warm dark matter
(WDM; e.g., Bode et al., 2001), fuzzy dark matter (e.g., Hu et al., 2000), and
self-interacting dark matter (SIDM; e.g., Spergel & Steinhardt, 2000). By
suppressing the small-scale matter power spectrum in the early universe and/or
stimulating halo disruptions at later times, these alternative dark matter
models have shown to reduce the number of subhalos around the Milky Way
(MW)-mass halos (e.g., Dunstan et al., 2011; Nadler et al., 2021).

On the other hand, it is possible that baryonic processes could suppress the
formation of some dwarf galaxies or make them difficult to observe, which
could explain the missing satellite problem (D’Onghia et al., 2010; Brooks et
al., 2013; Brooks & Zolotov, 2014; Sawala et al., 2016a; Wetzel et al., 2016;
Applebaum et al., 2021). In such cases, dark matter halos may still exist, but
may not have formed visible dwarf galaxies due to the effects of baryonic
physics. This is a possible solution to the missing satellite problem within
the framework of the $\Lambda$CDM paradigm. In addition, many authors have
shown that low-mass halos could be easily disrupted by baryon-induced physics
such as cosmic reionization, tidal stripping, ram pressure stripping, and
stellar feedback (Zhu et al., 2016; Simpson et al., 2018). In the Latte
simulations with the FIRE star formation and feedback model, the dwarf galaxy
population near the MW/M31-mass halo was found to agree well with the observed
population in the Local Group (Wetzel et al., 2016). Meanwhile, the “Mint”
resolution DC Justice League suite of MW-like zoom-in simulations showed that
the number of satellite galaxies matches the observed population of the dwarf
galaxies around MW-sized galaxies down to the ultrafaint dwarf regime (UFD;
Applebaum et al., 2021). Some studies have also shown that $\Lambda$CDM
simulations can reproduce the radial distribution of MW satellites (Santos-
Santos et al., 2018; Garrison-Kimmel et al., 2019; Samuel et al., 2020).
Moreover, the number of observed faint galaxies has increased recently (for
reviews, see Simon, 2019), which partially mitigates the missing satellite
problem. These findings suggest that the missing satellite problem is very
close to being solved. In fact, some researchers such as, Kim et al. (2018)
and Sales et al. (2022), argue that the problem is resolved.

Ideally, we would then expect the satellite galaxy populations to be
consistent regardless of the simulation code utilized. Nevertheless, due to
differences in the inherent properties of the simulations such as the adopted
physics models and the implementations of the gravity solver, discrepancies
may arise between codes (O’Shea et al., 2005; Heitmann et al., 2008). For
example, Elahi et al. (2016) studied subhalos and galaxies in a galaxy cluster
produced by multiple simulation codes, and found that in dark matter-only
(DMO) simulations, the population and properties of subhalos show good
agreements across code platforms. Nevertheless, they also discovered that the
codes produce significantly different galaxy populations when baryonic physics
models were included. While they found both similarities and disparities in
the galaxy population, the comparison of dwarf galaxy populations ($M_{\rm
halo}<10^{10}\,{\rm M}_{\odot}$) was not feasible due to the limited
resolution of their simulations. Indeed, there is an urgent need for
controlled comparisons of the dwarf galaxy populations produced by different
simulation codes. Such comparisons will be essential to understand the
robustness of the satellite galaxy populations predicted in the simulations,
and how sensitive they are with respect to the specific numerical methods and
assumptions adopted in the simulations.

The AGORA High-resolution Galaxy Simulations Comparison Project (Assembling
Galaxies of Resolved Anatomy) has aimed at collectively raising the predictive
power of numerical galaxy formation simulations, by comparing high-resolution
galaxy-scale calculations across multiple code platforms, using a DMO galaxy
formation simulation (Kim et al., 2014, hereafter Paper I), an idealized disk
galaxy formation simulation (Kim et al., 2016, hereafter Paper II), and a
fully cosmological zoom-in galaxy formation simulation (Roca-Fàbrega et al.,
2021; Roca-Fàbrega, 2023, hereafter Papers III and IV). In this paper, we
analyze the satellite halos around the target MW-like halo in the AGORA
“CosmoRun” simulation suite introduced and studied in Papers III and IV.
Specifically, we compare the eight hydrodynamic CosmoRuns and eight DMO
simulations, all performed with the state-of-the-art galaxy simulation codes
widely used in the numerical galaxy formation community, and study the
populations of their satellite halos and galaxies. We choose slightly
different redshift epochs near $z=2$ for each code in order to compare the
runs at the same dynamical stage in the target halo’s evolution history (see
Section 2.1 for details). We then compare the number of satellite halos in
CosmoRuns with its counterpart in the DMO simulations. We also explore the
consistency between the codes in other properties of satellite galaxies,
including the stellar mass$-$halo mass relation as well as the
mass$-$metallicity relation.

This paper is organized as follows. Section 2 describes the AGORA CosmoRun and
the DMO simulation, as well as the definition of a satellite halo. In Section
3, the satellite halo and galaxy populations in the CosmoRuns are presented in
comparison with those in the DMO runs. In Section 4, based on our results we
predict the satellite galaxy population at $z\sim 0$, and test inter-code
convergence in other satellite properties. Finally, we conclude the paper in
Section 5.

##  2 Methodology

###  2.1 The AGORA “CosmoRun” Simulation Suite

The CosmoRun described in Paper III is a suite of high-resolution cosmological
zoom-in simulations of a MW-mass halo ($10^{12}\,{\rm M}_{\odot}$ at $z=0$) on
multiple code platforms.111For publicly available datasets, visit
http://www.AGORAsimulations.org or http://flathub.flatironinstitute.org/agora.
The simulations analyzed herein started from a cosmological initial condition
at $z=100$ and reached $z\lesssim 2$. The adopted cosmological parameters are
$\Omega_{\Lambda}$ = 0.728, $\Omega_{\rm matter}=0.272$, $\Omega_{\rm
DM}=0.227$, $\sigma_{8}=0.807$, $n_{\rm s}=0.961$, and $h=0.702$. The code
groups participating in this particular comparison encompass both particle-
based and mesh-based codes: Art-I, Enzo and Ramses are mesh-based codes
whereas Changa, Gadget-3, Gear, Arepo-t and Gizmo are particle-based
codes.222We classify the SPH codes (Changa, Gadget-3, Gear) and the arbitrary
Lagrangian-Eulerian codes (Arepo and Gizmo) as particle-based codes. Galaxy
formation has been studied using both approaches, each with its own advantages
and disadvantages. After a series of calibration steps, all the codes in Paper
III reached an overall agreement in the stellar properties of the target halo,
and in its mass assembly history. The final CosmoRun suite includes common
baryonic physics modules in AGORA such as the Grackle radiative gas cooling
(Smith et al., 2017), cosmic ultraviolet background radiation (Haardt & Madau,
2012), and star formation, as well as the code-dependent physics including —
most notably — stellar feedback prescriptions. Both code-independent and code-
dependent physics implemented in each code are explained in great detail in
Paper III (some in Paper II). We update the two models from Paper III, Art-I
and Changa, to include weaker stellar feedback. In Art-I, we change the
condition for the minimum time step at high redshifts to achieve better
convergence in the halo growth history. We also incorporate a new model using
the Arepo code into our analysis. We refer to this as Arepo-t, which
represents the Arepo code with thermal feedback. The differences between the
old and new Art-I and Changa models, as well as the details of the Arepo-t
model, are illustrated in Paper IV.333In the analysis presented in Section
4.1, we used the older Art-I model, labeled as Art-I (old), which is described
in Paper III, because the new model has not reached $z\lesssim 1$. The results
with both models are mostly consistent.

The gravitational force softening length for the particle-based codes in the
highest-resolution region is 800 comoving pc until $z=9$ and 80 proper pc
afterward. Meanwhile, the finest cell size of the mesh-based codes is set to
163 comoving pc, or 12 additional refinement levels for a $128^{3}$ root
resolution in a $(60\,\,\,{\rm comoving}\,\,h^{-1}\,{\rm Mpc})^{3}$ box. A
cell is adaptively refined into 8 child cells on particle over-densities of 4.
For details on runtime parameters, we refer the readers to Paper III.

While all the AGORA CosmoRun simulations were calibrated to produce similar
stellar masses in the host halo by $z=4$ (see Section 5.4 and Figure 12 in
Paper III), we find that the host halo in some codes’ CosmoRun are at a
different stage in its dark matter accretion history from others’ at $z=2$.
This is likely due to the inter-code “timing discrepancy” (see Section 5.3 in
Paper I for more information). Because the halos in different codes are at
different evolutionary stages, the satellite halo abundances are also
different among the CosmoRun. To resolve this timing discrepancy, we have
created a merger tree for each code and selected an epoch near $z=2$
(hereafter called “$z\sim 2$”) for each code so that the target halo is in the
same stage in its merger history (for more information, see Paper IV). The
list of epochs for each code used for the present paper is in Table 1.
Snapshots of the CosmoRun simulations at $z\sim 2$ are shown in Figure 1.

Code | Redshift epoch
---|---
CosmoRun | Dark matter-only (DMO) run
Art-I | 1.85 | 2.18
Enzo | 2.29 | 2.15
Ramses | 2.21 | 2.12
Changa | 2.08 | 2.09
Gadget-2/3 | 2.13 | 2.05
Gear | 1.88 | 1.87
Arepo-t | 1.98 | 2.11
Gizmo | 2.02 | 2.11
Table 1: The redshift epoch selected for each code to be analyzed in this
paper. At these epochs, the eight CosmoRuns are in the same stage in the
target halos’ merger history. See Section 2.1 for details.

###  2.2 The Dark Matter-Only (DMO) Simulations

In order to investigate the role of baryonic physics adopted for AGORA in the
satellite halo population, we have also performed DMO simulations using the
same zoom-in initial condition generated with Music (Hahn & Abel, 2011) but
with no gas component. Accordingly, the mass of the dark matter particles in
the DMO runs is $\Omega_{\rm matter}/\Omega_{\rm DM}=1.20$ times heavier than
that in the CosmoRun. While Paper I found that the dark matter properties and
the satellite halo populations are nearly identical across all participating
codes in AGORA, there remained a systematic discrepancy in the satellite halo
populations in the low-mass end. Therefore, we have employed all eight codes
in AGORA to run DMO simulations to check their consistency.444In terms of the
gravity solver for collisionless components, Gadget-2 (latest version in 2011)
and Gadget-3 (first introduced in Springel et al., 2008) are identical for our
purpose, and will produce practically identical results in the DMO runs.
Snapshots of these DMO runs $z\sim 2$ are also included in Figure 1. The
runtime parameters governing the collisionless dynamics in the DMO runs are
set to be identical to those used for the CosmoRun.

Figure 1: The dark matter surface densities at $z\sim 2$ (the exact redshift
in each code in Table 1) for eight hydrodynamic “CosmoRun” simulations and
eight dark matter-only (DMO) simulations, projected through a 1.8 comoving Mpc
thick slab, with the target host halo’s virial radius $R_{\rm vir}$ drawn in a
white circle. See Section 2 for more information on these simulations.
Simulations performed by: Santi Roca-Fàbrega (Art-I, Ramses, and Art-I-dmo),
Ji-hoon Kim (Enzo), Johnny Powell and Héctor Velázquez (Changa and Changa-
dmo), Kentaro Nagamine and Ikkoh Shimizu (Gadget-3), Loic Hausammann and Yves
Revaz (Gear and Gear-dmo), Anna Genina (Arepo-t, and Arepo-dmo), Alessandro
Lupi and Bili Dong (Gizmo), Hyeonyong Kim (Enzo-dmo, Ramses-dmo, Gadget-2-dmo,
and Gizmo-dmo). Note that the mean dark matter surface densities in DMO runs
are $\Omega_{\rm matter}/\Omega_{\rm DM}=1.20$ times higher since it includes
the contribution from baryons. The high-resolution versions of this figure and
article are available at the Project website,
http://www.AGORAsimulations.org/. Figure 2: The dark matter surface densities
at $z\sim 2$ with the halos identified by the Rockstar halo finder drawn in
white circles whose radii indicate $0.5R_{\rm vir}$. Only the halos located
within 300 comoving kpc from the host halo’s center and those more massive
than $10^{7}h^{-1}\,{\rm M}_{\odot}$ in dark matter are drawn. Readers can
readily see that the DMO runs have more satellite halos than the CosmoRuns.
See Section 3.1 for more information.

###  2.3 Halo Finding

Halos in the CosmoRun and the DMO runs are identified with the Rockstar halo
finder (Behroozi et al., 2013) using only the highest-resolution dark matter
particles (i.e., not stellar or gas particles). We further narrow down the
identified halos to satellite halos using the following criteria: (i) it must
reside within 300 comoving kpc from the target host halo (100 proper kpc at
$z=2$; similar to the virial radius, $R_{\rm vir}$, of our host halo at
$z=0$), and (ii) it must be more massive than $10^{7}h^{-1}\,{\rm M}_{\odot}$
in dark matter (equivalent to 45 dark matter particles in the DMO
runs).555Note that there is no velocity criteria or requirement when
identifying satellites. Therefore, some halos may be counted as satellites
despite not being gravitationally bound to the host halo. We follow Bryan &
Norman (1998) definition of virial radius and mass.

For our analysis in Sections 3.4 and 4, we assign a stellar particle to a halo
following the process in Samuel et al. (2020). We first identify all stellar
particles located within $0.8R_{\rm vir}$ from the halo, with their velocities
relative to the halo less than twice the halo’s maximum circular velocity. We
then calculate the radius that encompasses 90% of the stellar particles
($R_{90}$) and the stellar velocity dispersion ($\sigma_{\rm vel}$). To
further refine our selection, we narrow down the stellar particle list to
those satisfying two more conditions: (1) they are located within
$1.5\,R_{90}$ from the center of mass of the halo and stellar particles, and
(2) their velocities relative to the halo is less than $2\sigma_{\rm vel}$. We
then iterate the analysis, recalculating $R_{90}$ and $\sigma_{\rm vel}$ for
the selected member particles until they converge within 99% of the previous
values. We start from the most massive satellite halos to lower ones, making
sure not to reassign stellar particles that have already been allocated.
Finally we define satellite “galaxies” as those whose stellar masses are at
least six times the approximate mass resolution of stellar particles (i.e.,
$M_{\rm star}>6m_{\rm gas,\,IC}=2.38\times 10^{5}h^{-1}\,{\rm M}_{\odot}$; see
Section 3.1 of Paper III).

##  3 Results

###  3.1 Satellite Halo Populations At $z\sim 2$

Figure 2 shows the dark matter surface density plots in a (600 comoving kpc)2
box, with the target host halo and the satellite halos drawn in white circles
(whose radii indicate half the virial radii, $0.5R_{\rm vir}$). One can
already observe that the eight hydrodynamic CosmoRuns have produced similar
numbers of dark matter halos with similar $R_{\rm vir}$’s for the host halo.
Readers can also see that the DMO runs clearly have more satellite halos than
the CosmoRuns.

To quantitatively study the differences in the participating simulations, in
Figure 3 we plot the cumulative number of satellite halos at $z\sim 2$ in
their dark matter mass, $N_{{\rm halo}}(>M)$ (left panels), and in radial
distance from the host halo’s center, $N_{{\rm halo}}(<r)$ (right panels). It
is worth noting several points:

  * •

First, we find that all eight hydrodynamic CosmoRuns have fewer satellite
halos than the DMO runs do across all halo masses and radii. In the halo mass
function (left panels of Figure 3), the numbers of satellite halos in all
CosmoRuns are systematically fewer than those in the DMO runs by a factor of
$\sim\,2$ for $M_{\rm halo}$ (halo dark matter mass) $<10^{8.5}h^{-1}\,{\rm
M}_{\odot}$. To put it differently, the ratios of the number of the CosmoRun
satellite halos to that in the DMO run (the mean number of halos in the eight
DMO runs) in each mass bin, $N_{\rm halo}/\langle N_{\rm halo,\,DMO}\rangle$,
is $\sim\,$0.5 (bottom left panel).

  * •

Second, the ratio of the satellite halos’ radial distribution function in the
CosmoRun to that in the DMO run, $N_{\rm halo}/\langle N_{\rm
halo,\,DMO}\rangle$, tends to become small — often zero — in the bin closest
to the host halo’s center, $r<40\,\,\,{\rm comoving\,\,kpc}$ (bottom right
panel of Figure 3). This implies that the causes of the deficit — the effect
of baryonic physics which we will explore in depth in Section 3.2 — have a
stronger influence near the host halo’s center. This is consistent with the
findings of earlier studies (e.g., Brooks & Zolotov, 2014; Wetzel et al.,
2016; Sawala et al., 2017; Garrison-Kimmel et al., 2017; Kelley et al., 2019).

  * •

The satellite halo populations in the eight DMO runs are slightly different
but are in general agreement with one another in both mass and space (upper
panels of Figure 3; lines with a reduced stroke width and darker colors).
Among the DMO runs, no systematic difference exists between the mesh-based and
particle-based codes, a result somewhat different from the earlier studies
(e.g., O’Shea et al., 2005; Heitmann et al., 2008) or from our findings in
Paper I.666In Paper I, systematic difference between particle-based and mesh-
based codes at the low-mass end was observed. It was because mesh-based codes
tend to have coarser force resolution as they attempt to resolve minute
density fluctuations in the outskirts of the target halo at high redshift,
resulting in a difference in the abundance of low-mass satellite halos. Note
also that the analysis in Paper I was carried out with the HOP halo finder, a
different choice from the Rockstar halo finder for what is presented here,
which could produce different numbers of halos identified. Instead, Enzo-dmo
has a slightly higher number of halos compared to the particle-based codes,
and Ramses-dmo is in the middle of the pack of particle-based codes. The
numerical resolution in the highest-resolution region of the CosmoRun — and
correspondingly, our new DMO runs — is chosen to resolve the interstellar
medium (ISM) and the star-forming regions in them (i.e., $\lesssim 100$ proper
pc at all times between $z=100$ and 2). The high resolution in our simulation
suite might have been sufficient for Enzo-dmo and Ramses-dmo to alleviate any
discrepancy previously observed in the satellite halo population between the
particle-based and mesh-based DMO runs. Note that Art-I-dmo seems to have had
a slightly harder time fully resolving the outskirts of our target
halo.777Note that runtime parameters are not chosen to match the refinement
structure between the CosmoRun and the DMO run. In a typical mesh-based DMO
run, cells are adaptively refined only by dark matter mass, whereas in a
CosmoRun they are refined not only by dark matter mass but also by baryon mass
and others. As a result, Art-I and Art-I-dmo may have different refinement
structure, especially in the outskirts of the target halo. Even so, the
difference is not as severe as what was seen in the previous studies.

  * •

Two CosmoRuns, Art-I and Enzo, have smaller satellite halo populations than
the rest of the participating codes do, especially in the low-mass end
($M_{\rm halo}<10^{7.5}h^{-1}\,{\rm M}_{\odot}$) and in the outskirts of the
host halo ($r>200\,\,\,{\rm comoving\,\,kpc}$). And the number of halos in the
three mesh-based CosmoRuns tends to be lower in the range of halo masses
$M_{\rm halo}\sim 10^{8.5}h^{-1}\,{\rm M}_{\odot}$ and radial distances
$\left[100,\,150\right]$ comoving kpc. The inter-code difference among the
CosmoRuns, which their counterpart DMO runs do not exhibit, should be
attributed to how the same (or similar) baryonic physics are treated
differently in the two hydrodynamics approaches.

Figure 3: The cumulative number of satellite halos at $z\sim 2$ in their dark
matter mass, $N_{{\rm halo}}(>M)$ (left), and in radial distance from the host
halo’s center, $N_{{\rm halo}}(<r)$ (right). We distinguish mesh-based codes
(solid lines) and particle-based codes (dashed lines) with different line
styles, and CosmoRuns and DMO runs with different brightness. The bottom
panels display the ratio of the number of the CosmoRun’s satellite halos to
that in the DMO run (the mean value of the eight DMO runs) in each mass/radius
bin. All hydrodynamic CosmoRuns have fewer satellite halos than the DMO runs
do across all halo masses and radii. See Section 3.1 for more information.

One of the most notable findings among the above is that all hydrodynamic
CosmoRuns have produced fewer satellite halos than the DMO runs have by $z\sim
2$ across all halo masses and radii. We further study in Section 3.4 that the
so-called “missing satellite problem” (over-abundance of satellite halos in
simulations; Kauffmann et al., 1993; Klypin et al., 1999; Moore et al., 1999;
Benson et al., 2002) could be easily resolved in all participating codes
simply by implementing the baryonic physics adopted for AGORA in simulations
with sufficient numerical resolution ($\lesssim 100$ proper pc at $z=2$) by
examining the satellite galaxy populations around the target host halo. For
now, in Sections 3.2 and 3.3 we focus on the causes of differences seen in
Figure 3 — between the satellite halos in the CosmoRuns and in the DMO runs.

Figure 4: The evolution of the number of satellite halos across cosmic time,
for two different dark matter mass bins, $10^{7}h^{-1}\,{\rm M}_{\odot}<M_{\rm
halo}<10^{8}h^{-1}\,{\rm M}_{\odot}$ (left) and $M_{\rm
halo}>10^{8}h^{-1}\,{\rm M}_{\odot}$ (right). We only count halos that reside
within 300 comoving kpc from the target host halo. The top panels highlight
the CosmoRuns with the DMO runs shown at reduced opacity, while the bottom
panels emphasize the DMO runs with the CosmoRuns at reduced opacity. All
hydrodynamic CosmoRuns have fewer satellite halos than the DMO runs do in both
mass bins at nearly all redshifts. The mesh-based CosmoRuns tend to host
slightly fewer lower-mass satellite halos than the particle-based CosmoRuns do
at most redshifts (top left panel). See Section 3.2 for more information.

###  3.2 Evolution of Satellite Halo Populations

We now study the halo populations at multiple epochs from $z=15$ to $\sim 2$
to understand and discriminate various causes that have affected the satellite
halo populations. Figure 4 shows the evolution of the number of satellite
halos in two different mass bins, $10^{7}h^{-1}\,{\rm M}_{\odot}<M_{\rm
halo}<10^{8}h^{-1}\,{\rm M}_{\odot}$ (left panels) and $M_{\rm
halo}>10^{8}h^{-1}\,{\rm M}_{\odot}$ (right panels). Readers can notice that
the hydrodynamic CosmoRun simulations and the DMO runs show systematic
differences in both mass bins at nearly all redshifts. In the left panels of
Figure 4, one may spot the systematic disagreements between the DMO runs, the
group of particle-based CosmoRuns, and the group of mesh-based CosmoRuns,
already at $z=15$. From $z=8$ to 4, the numbers of halos in the CosmoRuns
barely grows in both mass bins, while those in the DMO runs increase steadily.
The number of halos in the CosmoRuns in the higher-mass bin (right panels;
$M_{\rm halo}>10^{8}h^{-1}\,{\rm M}_{\odot}$) remain approximately constant
from $z=3$ to $\sim 2$, whereas those in the DMO runs continue to increase. In
the meantime, just as in Figure 3, there exists disagreement between the
particle-based codes (colored dashed lines in the top panels of Figure 4) and
the mesh-based codes (colored solid lines), especially in the lower-mass bin
(top left panel; $10^{7}h^{-1}\,{\rm M}_{\odot}<M_{\rm
halo}<10^{8}h^{-1}\,{\rm M}_{\odot}$).

Now we investigate various causes for these differences in time:

Figure 5: The ratio of the dark matter mass of an individual halo in the
CosmoRun to that of its DMO counterpart ($M_{\rm halo}/M_{\rm DMO}$). Each
symbol represents a CosmoRun – DMO run pair at four different epochs, from
$z=12$ to $\sim 2$. A thick solid line is the median value in each mass bin
for each CosmoRun. The solid horizontal line denotes $\Omega_{\rm
DM}/\Omega_{\rm matter}=0.83$. In all CosmoRuns we analyzed, the median value
of $M_{\rm halo}/M_{\rm DMO}$ is less than 0.83. This means that the halos
originated from the same patch in the initial condition but under the
influence of baryonic physics did not grow as much in mass as their DMO
counterparts did. See Section 3.3 for more information.

  * •

As early as at $z=12$, the CosmoRuns tend to have fewer halos than the DMO
runs, even before the cosmic reionization begins or the extragalactic
ultraviolet background radiation is turned on in Grackle.888We set the Grackle
parameters UVbackground_redshift_on $=15$ and UVbackground_redshift_fullon
$=15$. It is especially true in the higher-mass bin (right panels; $M_{\rm
halo}>10^{8}h^{-1}\,{\rm M}_{\odot}$). At $z\sim 15$, smaller dark matter
halos have difficulties at keeping baryons because the density fluctuation of
gas is smoother than that of dark matter on small scales (Gnedin & Hui, 1998).
Thus, these smaller halos have a smaller enclosed baryon mass than the cosmic
average (O’Leary & McQuinn, 2012). This leads to lower halo masses in the
CosmoRuns at $z=15$, therefore, fewer halos in Figure 4. From $z=15$ to 8 the
difference between the CosmoRuns and the DMO runs persists in both mass bins.

  * •

Again as early as at $z=12$, one can observe a discrepancy in the CosmoRuns
between particle-based codes and mesh-based codes in the lower-mass bin (top
left panel; $10^{7}h^{-1}\,{\rm M}_{\odot}<M_{\rm halo}<10^{8}h^{-1}\,{\rm
M}_{\odot}$). Particle-based codes show a tendency to have more halos than
mesh-based codes do until $z\sim 6$ when Ramses begins to behave like
particle-based codes. It is well documented that the particle-based codes may
produce more satellite halos due to the so-called “gas – dark matter particle
coupling” in the early universe (Yoshida et al., 2003; O’Leary & McQuinn,
2012). When a gas particle is close to a nearby dark matter particle, the gas
particle could be captured in the dark matter particle’s potential well. This
gas particle now obtains an artificial velocity that follows that of the dark
matter particle, resulting in an increased power on small scales. Even a
minute gas – dark matter two-particle coupling could be a source of
numerically-driven fluctuation, particularly in the early universe that is
nearly homogeneous. While this artifact may be alleviated with adaptive
gravitational softening, the particle-based codes in the AGORA CosmoRun suite
adopted a fixed gravitational softening length (see Section 2), prone to
overproduction of satellite halos.999A new type of cosmological initial
conditions generated with a higher-order Lagrangian perturbation theory may
provide another solution to this problem (Michaux et al., 2021). It will
enable us to start our simulation at $z\simeq 15$, much later than $z=100$ as
in the CosmoRun, bypassing the gas – dark matter coupling problem at high
redshift entirely. On the other hand, some DMO runs in mesh-based codes may
have coarser force resolution at high redshift as they attempt to resolve
small density fluctuations in the outskirts of the target halo (O’Shea et al.,
2005; Heitmann et al., 2008), leading to smaller numbers of halos in e.g.,
ARt-I-dmo and Ramses-dmo.

  * •

From $z=8$ to 4, reionization plays an important role in suppressing the
growth of satellite halos in the CosmoRun (see Section 2 and Paper III), when
other local baryonic physics mechanisms are yet to become effective. The
extragalactic photoionizing background radiation heats and removes the gas
prior to infall, and efficiently inhibits the growth of halos at $z\lesssim 8$
(Sawala et al., 2015; Qin et al., 2017).101010It is worth to remind the
readers that, in DMO runs, the mass of the gas is included in the dark matter
component, effectively making $\Omega_{\rm matter}=0.272=\Omega_{\rm DM}$ (see
Section 2.2). Therefore, while the gas experiences hydrodynamic forces such as
reionizing radiation and may “evaporate” in the CosmoRun, the gas mass
contribute in whole to the growth of the halo in the DMO run. Reionization is
relatively more effective on the low-mass halos ($v_{\rm circ,\,max}<20\,{\rm
km\,s^{-1}}$; Sawala et al., 2015; Zhu et al., 2016).

  * •

At later times, other baryonic effects enhance the depletion of substructures
when compared to the DMO counterparts. Gas in low-mass halos is removed by
ram-pressure stripping before the infall, along with the extragalactic
radiation field. Tidal stripping in the steep gravitational potential of the
host halo becomes important now, and significantly affects the satellite halo
population, especially in the intermediate-mass range ($20\,{\rm
km\,s^{-1}}<v_{\rm circ,\,max}<35\,{\rm km\,s^{-1}}$; D’Onghia et al., 2010;
Brooks et al., 2013; Brooks & Zolotov, 2014; Sawala et al., 2016b; Zhu et al.,
2016; Sawala et al., 2017; Garrison-Kimmel et al., 2017; Kelley et al., 2019).
However, tidal disruption induced by stellar bulge and disk in cosmological
simulations could be overestimated due to insufficient resolution (see Webb &
Bovy, 2020; Green et al., 2022).

Stellar feedback such as supernovae also expels the gas and impedes the halos’
mass growths (Brooks et al., 2013; Munshi et al., 2013; Velliscig et al.,
2014; Schaller et al., 2015; Fitts et al., 2017). These late-time baryonic
processes can explain the widening gap between the CosmoRuns and the DMO runs
at $z\lesssim 4$ in both left and right panels of Figure 4.

In summary, we find that baryonic processes cause all hydrodynamic CosmoRun
simulations to have fewer satellite halos than the DMO runs at nearly all
redshifts. While baryonic physics left only indirect signatures in the halos’
growth histories of dark matter masses ($M_{\rm halo}$), in Section 3.4 we
will see its more direct impact on the halos’ stellar components.

###  3.3 How Baryonic Physics Affects Each Individual Halo

Until now we have mainly focused on the population of satellite halos, and how
it changes with the inclusion of baryonic physics. To investigate how each
individual halo is actually affected by the baryonic processes, we now match
and compare halos in hydrodynamic simulations (e.g., Enzo CosmoRun) to their
counterparts in the DMO simulation (e.g., Enzo-dmo run). The matching process
is adapted from Schaller et al. (2015) and Lovell et al. (2021). For every
satellite halo in e.g., Enzo CosmoRun, we first identify the 40 dark matter
particles that are closest to the halo’s center. Since the particle IDs are
shared by the Enzo and Enzo-dmo run, we can locate these 40 particles in the
Enzo-dmo run. Then we search for a halo containing 50% or more of these
counterpart particles. Finally, by carrying out the same procedure in reverse,
another link is obtained — i.e., first find the 40 most bound particles in the
Enzo-dmo run, and then locate these particles in the Enzo CosmoRun. A pair of
two halos that are bijectively mapped (bidirectionally connected) between the
two simulations are considered as a “matched” pair. Particle IDs are
identically assigned in the initial conditions of Changa, Gagdet-3, Gear,
Arepo-t, Gizmo and their DMO counterparts, so we can similarly find matched
pairs in between the two codes. But because particle IDs in Art-I, Ramses and
their DMO counterpart simulations are not identically assigned in their
initial conditions, halos in these two simulations need to be matched with a
different method based on the distribution of dark matter particles at $z=100$
(see Appendix A for details).

We conjecture that various baryonic processes have slowed down the growth of
halos in hydrodynamic simulations compared to their DMO counterparts. To test
this hypothesis, in Figure 5, we plot the ratio of the dark matter mass of an
individual halo in the CosmoRun to that of its matched DMO counterpart
($M_{\rm halo}/M_{\rm DMO}$). A few observations to note:

  * •

In all eight CosmoRuns matched to their respective DMO counterpart, the ratio
$M_{\rm halo}/M_{\rm DMO}$ is on average less than $\Omega_{\rm
DM}/\Omega_{\rm matter}=0.83$ (marked with solid horizontal lines in Figure
5), where $\Omega_{\rm matter}=0.272$ and $\Omega_{\rm DM}=0.227$ (see Section
2.1). If the halo in the CosmoRun had followed the identical mass growth
history as that of its DMO counterpart, the dark matter mass of the CosmoRun
halo, $M_{\rm halo}$, should have been $M_{\rm DMO}\times\Omega_{\rm
DM}/\Omega_{\rm matter}=0.83\,M_{\rm DMO}$. The fact that the ratio lies below
0.83 means that the halos have smaller masses and smaller virial radii in all
hydrodynamic simulations when compared with their DMO counterparts. Although
the halos originated from the same patch in the initial condition, the ones
under the influence of baryonic physics did not grow as much in mass as their
DMO counterparts did.

  * •

The baryonic effects are present at all redshifts, slightly more so at later
times (i.e., the $M_{\rm halo}/M_{\rm DMO}$ ratio is smaller at $z\sim 2$ than
at $z=12$). The baryonic effects are the combination of early- and late-time
processes, such as reionization inhibiting the growth of small satellite halos
with shallow gravitational wells, and the host halo’s tidal field stripping
the halos of existing gas, as discussed in Section 3.2.

  * •

Among the CosmoRuns, the mesh-based codes tend to have lower $M_{\rm
halo}/M_{\rm DMO}$ values than the particle-based codes do in general,
especially at low-mass end. The discrepancy between the two code groups is
considerable already at $z=12$, which is in line with the over-abundance of
satellite halos in particle-based codes discussed in Section 3.2. Furthermore,
consistent with the patterns observed in Figures 3 and 4, Art-I and Enzo have
smaller ratios compared to the other codes.

Figure 6: The fraction of halos whose counterparts in the DMO simulations are
identified at $z\sim 2$, averaged over all eight CosmoRuns. Error bars
indicate one standard deviation. The “match fraction” is lower for low-mass
halos because the matching criteria is more demanding for them. It implies
that the low-mass matched halos are likely biased towards the ones with
quiescent accretion histories. See Section 3.3 for more information.

It should be noted that only a small fraction of halos are matched with their
DMO counterparts in the low-mass end ($M_{\rm DMO}\lesssim 10^{8}h^{-1}\,{\rm
M}_{\odot}$). Figure 6 illustrates the fraction of halos whose counterparts in
the DMO simulations are identified at $z\sim 2$, averaged over all eight
CosmoRuns. Since our halo matching process is more demanding for halos with
fewer member particles, the “match fraction” is lower for low-mass halos. This
suggests that the low-mass matched halos are likely biased towards the ones
with quiescent accretion histories and without major mergers or disruptions in
the past, potentially resulting in an overestimated $M_{\rm halo}/M_{\rm DMO}$
ratio. Readers may also find it interesting that massive satellite halos
($M_{\rm DMO}\gtrsim 10^{9}h^{-1}\,{\rm M}_{\odot}$) have disappeared between
$z=7$ and 4 (second and third panel from the right in Figure 5). According to
the accretion history of the host halo, multiple mergers occur between $z=7$
and 4, which explains the disappearance of the massive satellites by $z=4$.

To summarize, we have shown that each individual halo tends to have a slower
mass accretion history until $z\sim 2$ in the CosmoRun than in its counterpart
DMO run. The discrepancy can be explained by early- and late-time baryonic
physics that slows down the growth of satellite halos.

Figure 7: The cumulative number of satellite galaxies at $z\sim 2$ in their
stellar mass, $N_{{\rm galaxy}}(>M)$ (left), and in their 3-dimensional
stellar velocity dispersion, $N_{{\rm galaxy}}(>\sigma_{\rm vel})$ (right). We
define galaxies as satellite halos that contain stellar particles. The number
of satellite galaxies is significantly fewer than that of satellite halos in
all eight CosmoRuns. Readers may compare $N_{\rm galaxy}$ with $N_{{\rm
halo}}$ (in Figure 3), or with the gray dotted line in the right panel above
denoting the average number of satellite halos in all CosmoRuns (plotted with
their dark matter velocity dispersion). The thick black solid line and dashed
line in both panels indicate the known present-day satellites around the Milky
Way (MW) and M31, respectively, which of course are lower limits to the true
numbers. See Section 3.4 for more information.

###  3.4 Satellite Galaxy Populations At $z\sim 2$

In Section 3.1 we have demonstrated that all hydrodynamic CosmoRun simulations
produce fewer satellite halos around our host halo than the DMO runs do across
all satellite masses and radii. To verify that this finding naturally leads to
the baryonic solution to the “missing satellite problem” (see Section 3.1)
that is independent of the numerical platform utilized, in this section, we
examine the satellite galaxy populations around the target host halo. Here we
define “galaxies” as satellite halos that contain stellar particles (see
Section 2.3 for more information).

In Figure 7 we plot the cumulative number of satellite galaxies at $z\sim 2$
in their stellar mass, $N_{{\rm galaxy}}(>M)$ (left), and in their
3-dimensional stellar velocity dispersion, $N_{{\rm galaxy}}(>\sigma_{\rm
vel})$ (right). In both panels, we restrict the satellite galaxies to those
with $M_{\rm star}>6m_{\rm gas,\,IC}=2.38\times 10^{5}h^{-1}\,{\rm M}_{\odot}$
(see Section 2.3). Several notable points are as follows:

  * •

By comparing with the mass function of satellite halos, $N_{{\rm halo}}$, in
Figure 3, or with the gray dotted line in the right panel of Figure 7 denoting
the average number of satellite halos in all CosmoRuns, one can readily see
that the number of satellite galaxies is significantly fewer than that of
satellite halos in all participating CosmoRuns. While baryonic physics left
indirect signatures in the halos’ growth histories of dark matter masses in
Sections 3.1 and 3.2, we can now see its more direct impact on the halos’
stellar components. All the baryonic processes discussed in Section 3.2 — such
as cosmic reionization, tidal stripping, ram pressure stripping, and stellar
feedback — act efficiently to halt or impede the stellar mass growth inside
the halo (Bullock et al., 2001; Revaz & Jablonka, 2018). For example, gas in a
low-mass halo with a shallow potential well is removed by ram pressure
stripping before its infall to the host halo, and by stellar feedback as
supernovae explode. Further, these processes can interact; for example,
supernova feedback can expel gas from halos which is then more easily removed
by ram pressure stripping. As a result, gas is depleted in most low-mass
satellite halos in the CosmoRuns, which end up with few stellar particles.

  * •

The thick black solid line and dashed line in both panels of Figure 7 indicate
the present-day satellites around the MW and M31, respectively (McConnachie,
2012).111111The latest compilation in 2021 can be found in https://www.cadc-
ccda.hia-iha.nrc-cnrc.gc.ca/en/community/nearby/. For the stellar velocity
dispersion of M33, however, we followed Quirk et al. (2022). To estimate the
stelar mass of galaxies, we assume a mass-to-light ratio of 1. Recent
observations find more satellite galaxies in the Local Group, yielding 50$-$60
satellites around the MW. However, most of these newly discovered galaxies are
UFDs, which are fainter than $M_{\rm V}=-7.7$ (or $L=10^{5}\,{\rm L}_{\odot}$;
Simon, 2019). Therefore, these newly discovered UFDs are beyond the scope of
the present study owing to the limited numerical resolution.,121212The
3-dimensional stellar velocity dispersions, $\sigma_{\rm vel}$, of the MW and
M31 satellite galaxies are estimated by the line-of-sight stellar velocity
dispersions multiplied by $\sqrt{3}$. The number of satellite galaxies around
M31 shown in two panels of Figure 7 differ slightly due to the lack of stellar
velocity information for two satellites. That is, among the 19 observed
satellites around M31, 16 have available stellar velocity information.
Although readers should be cautioned that we are comparing two datasets at
different epochs, one can observe that the satellite galaxy populations in the
CosmoRuns at $z\sim 2$ are largely consistent with those of the MW and M31 at
$z=0$ in their stellar masses and velocity dispersions. For more on how we
attempt to compare the satellite galaxy populations at $z\sim 0$, see Section
4.1 and Figure 8.

  * •

The agreement amongst the satellite galaxy populations of the eight CosmoRuns
is better in the right panel ($N_{{\rm galaxy}}(>\sigma_{\rm vel})$ in stellar
velocity dispersion) than in the left panel ($N_{{\rm galaxy}}(>M)$ in stellar
mass). It is because the velocity dispersion serves as a better and more
useful proxy for the dynamical mass of a system, as it reflects the
gravitational impact of the underlying dark matter halo (not just the stellar
component of a galaxy). While the difference in the satellite galaxy
population amongst the CosmoRuns is more pronounced when considering their
stellar mass, there remains a good overall agreement. In Section 4.2, we
further investigate the relationship between the dark matter mass and stellar
mass of the satellite halos.

To sum up, we have found that the number of satellite galaxies is
significantly fewer than that of dark matter halos in all CosmoRun
simulations, and is comparable to the number of present-day satellites near
the MW or M31. The so-called “missing satellite problem” is resolved in all
participating codes simply by implementing the baryonic physics adopted for
AGORA in simulations with sufficient numerical resolution ($\lesssim 100$
proper pc at $z=2$). We argue that various baryonic processes make the
CosmoRuns have far fewer satellite galaxies than the satellite dark matter
halos in the DMO runs. Future studies tracing the star formation history and
the trajectory of each halo will tell us which baryonic mechanism acts most
prominently and when.

Figure 8:  The cumulative number of satellite galaxies at $z=0.3$ in their
stellar mass, $N_{{\rm galaxy}}(>M)$ (left), and in their 3-dimensional
stellar velocity dispersion, $N_{{\rm galaxy}}(>\sigma_{\rm vel})$ (right),
for the Art-I (old), Enzo, Gadget-3 and Gear CosmoRuns. The grey area
represents the minimum and maximum values of the cumulative number of
satellite galaxies at $z\sim 2$ from those codes. The numbers of satellite
galaxies at $z=0.3$ are only a factor of $\lesssim 2$ larger than those at
$z\sim 2$, or almost identical in abundance in some codes. As the latest Art-I
CosmoRun has not reached $z=0.3$, we use the older model for the Art-I code,
denoted as Art-I (old), which is represented in Paper III. See Section 4.1 for
more information.

##  4 Discussion

###  4.1 Predicting Satellite Galaxy Populations At Lower Redshifts

In Section 3, we choose to study the satellite halo populations at $z\sim 2$
because not all the AGORA CosmoRuns have yet reached $z=0$. This approach, of
course, has limitations, as the majority of satellite halos at $z\sim 2$ are
likely to undergo mergers or be disrupted by $z=0$. Here we describe what form
of satellite galaxy populations we expect to see at lower redshifts.

In the DMO simulations, the number of satellite halos tends to increase over
time, which may lead one to conclude that a higher satellite halo abundance is
expected at $z\sim 0$ than at $z\sim 2$. However, as illustrated in Section
3.2, baryonic physics may disrupt halos, considerably reducing the number of
satellite halos in the hydrodynamic simulations at lower redshift. And
particularly because the AGORA CosmoRun adopted an initial condition of a halo
with a quiescent merger history after $z=2$, the number of newly accreted
satellite halos may be small after $z=2$. Therefore, we expect that the
satellite galaxy population in the CosmoRuns would not change dramatically
from $z\sim 2$ to $z\sim 0$.

In order to verify this, in Figure 8 we plot the satellite halo population at
$z=0.3$ as a function of stellar mass and stellar velocity dispersion for the
five codes that reached the epoch already.131313Art-I (old) in Figure 8
represents the older Art-I model used in Paper III. The satellite galaxy
population at $z\sim 2$ is almost identical between the two models. However,
the older model, Art-I (old), exhibits a more severe inter-code timing
discrepancy. The numbers of satellite galaxies at $z=0.3$ for these codes are
either only a factor of $\lesssim 2$ larger than those at $z\sim 2$, or almost
identical (as in the case of Enzo and Arepo-t). However the inter-code
differences have greatly increased. For instance, Enzo and Arepo-t have no
satellite galaxies with $M_{\rm star}>10^{7}h^{-1}\,{\rm M}_{\odot}$ at that
epoch, while Gadget-3 and Gear each have six satellite galaxies exceeding that
stellar mass. Compared by 3-dimensional stellar velocity dispersion, Enzo has
no satellite galaxies with $\sigma_{\rm vel}>40$ km/s, while Gear has five
satellite galaxies exceeding that velocity dispersion. Despite these
differences, the results still align well with the observed satellite galaxy
populations for the Milky Way (MW) and M31, except that Enzo shows a slightly
lower population in both stellar mass and velocity dispersion, and Arepo-t
exhibits a reduction in stellar mass. We conclude that our findings in Section
3 for $z\sim 2$ will likely also hold at $z\sim 0$.

Figure 9: The stellar mass$-$halo mass relation of the satellite (left) and
field (right) galaxies at $z\sim 2$. $y$-axis indicates the mean value of
stellar masses in each mass bin. The thick grey dotted, dashed and dot-dashed
lines are for dwarf galaxies in other zoom-in simulations at $z=0$ (FIRE-2,
Auriga and DC Justice League, respectively; Hopkins et al., 2018; Grand et
al., 2021; Munshi et al., 2021), while the thin black and grey solid lines
without markers are semi-empirical models for $2<z<2.5$ with extrapolation to
low-mass galaxies (Girelli et al., 2020; Legrand et al., 2019). The
relationship with a 68% confidence interval, as constrained by Nadler et al.
(2020) using Milky Way satellites, is represented in the blue shaded region.
All CosmoRuns produce similar relations with no systematic discrepancy between
mesh-based codes and particle-based codes. See Section 4.2 for more
information. Figure 10: The mass$-$metallicity relation of the satellite
(left) and field (right) galaxies at $z\sim 2$ using stellar (top) and gas
(bottom) metallicity. $y$-axis indicates the mean value of stellar
metallicities in each mass bin. A systematic difference exists between mesh-
based codes (solid lines) and particle-based codes (dashed lines). For
comparison, the thick grey (red) dotted line represents the stellar mass$-$
stellar (gas) metallicity relation at $z=2$ that best fits the field galaxies
in the FIRE-2 simulation (Ma et al., 2016), while the thick grey (brown) dot-
dashed line represents the same relation at $z=2$ in the TNG50 simulation
(Pillepich et al., 2019; Nelson et al., 2019). We also include the observed
mass-metallicity relation of dwarf galaxies in the Local Group, represented by
grey crosses with error bars (top panels; Sanati et al., 2023), as well as the
mass-metallicity relation observed by JWST in a $z\sim 2$ galaxy cluster
field, shown as a thick light blue line (bottom panels; Li et al., 2022). See
Section 4.2 for more information.

###  4.2 Testing Inter-code Convergence In Satellite Properties: The Stellar
Mass$-$Halo Mass Relation And The Mass$-$Metallicity Relation

We now explore the inter-platform convergence in satellite properties amongst
the CosmoRuns by studying the two relations that probe the baryonic physics:
the stellar mass$-$halo mass relation and the mass$-$metallicity relation.
This will allows us to verify the realism of the AGORA baryonic physics in the
CosmoRun.

First, in the left panel of Figure 9, we show the stellar mass$-$halo mass
relation at $z\sim 2$ of the satellite galaxies identified in Sections 2.3 and
3.4. For the completeness of our analysis, in the right panel of Figure 9, we
also draw the same plot using the field galaxies found in our simulations.
This is possible thanks to the sufficiently large zoom-in region around the
host halo that contains 20$-$30 field galaxies at $z\sim 2$.141414In contrast
to the satellite halos defined in Section 2.3, we define field halos using the
following criteria: (i) a field halo must reside beyond 300 comoving kpc of
our target host halo (or 100 proper kpc at $z=2$; a value similar to the
virial radius of our host halo at $z=0$), (ii) it must be more massive than
$10^{7}h^{-1}\,{\rm M}_{\odot}$ in dark matter, and (iii) it must not have a
parent halo in the Rockstar halo catalog (i.e., satellites of other halos are
excluded). And after assigning stellar particles to these halos using the
method described in Section 2.3, we plot only the field galaxies whose stellar
masses are heavier than $6m_{\rm gas,\,IC}=2.38\times 10^{5}h^{-1}\,{\rm
M}_{\odot}$, just as in Section 3.4. One may notice that, on average, the dark
matter halos of field galaxies are about 2.5 times more massive than the dark
matter halos of satellite galaxies for a given luminosity. The satellites’
halo masses do not grow after their infall to the host, or rather, decrease
due to tidal stripping. In the meantime, their stellar masses continue to grow
(Gunn & Gott, 1972; Behroozi et al., 2019). Thus, satellite galaxies tend to
have more stellar masses at a given halo mass than field galaxies do.

Different simulation codes display varied behaviors, but there is also
evidence of remarkable convergence. To begin, inter-code differences in the
mass-metallicity relation are evident. Art-I, Ramses, and Gear show a
relatively large $M_{\rm star}/M_{\rm halo}$ value, while the ratios for Enzo
and Gizmo are slightly smaller than those of other codes. This trend reflects
what was already discovered in the satellite galaxy populations (left panel of
Figure 7). Changa also shows a large $M_{\rm star}/M_{\rm halo}$ for the
satellite galaxies, which could arise from the insufficient number of
satellite galaxies. The field galaxies in Changa exhibit a relation consistent
with other codes (right panel of Figure 9).

However, despite these initial differences, the overall picture reveals
convergence. The differences in stellar mass$-$halo mass relations are within
1 dex for the field galaxies across all CosmoRuns with no visible systematic
discrepancy between mesh-based and particle-based codes. The common baryonic
physics adopted in AGORA and the stellar feedback prescription typically used
in each code group (calibrated to produce a similar stellar mass at $z=4$) are
responsible for this convergence, particularly because the simulations are
performed with sufficient resolution ($\lesssim 100$ proper pc at $z=2$).

We then compare our result with previous studies. The thick grey dotted,
dashed, and dot-dashed lines represent the relation for dwarf galaxies at
$z=0$ in the FIRE-2, Auriga, and the DC Justice League simulation (Hopkins et
al., 2018; Grand et al., 2021; Munshi et al., 2021, respectively).151515These
lines represent both satellite and field galaxies in both panels. In contrast
to the previous studies listed here that employ the total halo mass, the halo
mass $M_{\rm halo}$ in the present paper specifically refers to the mass of
dark matter in the halo. However, as the majority of mass in satellite
galaxies is dark matter, this slight difference in the mass definition does
not significantly affect Figure 9.161616Hopkins et al. (2023) show that the
latest FIRE model, FIRE-3, predicts a stellar mass up to a factor of ten
higher compared to the FIRE-2 model for dwarf galaxies with $M_{\rm
peak}\approx 10^{9}\,{\rm M}_{\odot}$. For the dwarf galaxies with stellar
masses $\gtrsim 10^{6}-10^{7}\,{\rm M}_{\odot}$, in contrast, there is little
difference in galaxy stellar masses (Hopkins et al., 2023). The blue shaded
region represents the relation inferred from Milky Way satellites (Nadler et
al., 2020). The thin black and grey solid lines are for the semi-empirical
models at $2<z<2.5$ with extrapolation to dwarf-sized galaxies (Girelli et
al., 2020; Legrand et al., 2019, respectively). The stellar masses in the
AGORA CosmoRuns are on average $\sim$ 0.5 dex higher than the empirical
predictions, but are largely consistent with the previous simulation studies
at $z=0$. The inter-code scatters in the low-mass halos ($M_{\rm halo}\lesssim
10^{9}h^{-1}\,{\rm M}_{\odot}$) is due to the complex interplay between
baryonic physics and different merger history of the halos, which cannot
easily be reproduced by abundance matching in the empirical models (Revaz &
Jablonka, 2018).171717In Figure 9 some codes do not reach the highest dark
matter mass bin, indicating an absence of satellite halos in that range.
Similarly, some codes do not reach the lowest dark matter mass bin because
there is no halo with a sufficient number of stellar particles in that mass
range. The galaxy$-$halo connection seen in Figure 9 indicates not only the
robustness and reproducibility of the participating simulations, but also the
realism of the AGORA CosmoRun baryonic physics.

Second, in Figure 10 we present the mass$-$metallicity relation at $z\sim 2$
for the satellite and field galaxies. Stellar masses and metallicities are
used to draw the plots in the top panels, while stellar masses and gas
metallicities are used in the bottom panels. For the bottom panels, we assign
a gas parcel (cell or particle) to a halo if it is within $0.15R_{\rm vir}$
from the halo’s center. Only galaxies whose gas mass is greater than three
times the approximate mass resolution of stellar particles, $3m_{\rm
gas,\,IC}=1.19\times 10^{5}h^{-1}\,{\rm M}_{\odot}$ (see Section 3.1 of Paper
III), are included in the bottom panels of Figure 10.

We note that there exists a small, but systematic difference between the mesh-
based and particle-based codes. Some mesh-based codes, Art-I and Ramses, tend
to show higher stellar metallicities than the particle-based codes do for the
satellite galaxies (left panels in Figure 10), and a similar trend exists for
the field galaxies, too (right panels). However, the differences between codes
are mitigated the high stellar mass end, $M_{\rm star}\gtrsim 10^{8}\,{\rm
M}_{\odot}$, where the mean values of the relation converge within $\sim 0.5$
dex for the field galaxies (right panels). Since our careful calibration of
stellar feedback for the CosmoRun has yielded similar star formation histories
across the participating codes (see Papers III and IV for detailed
discussion), the difference in metallicities is most likely due to the
difference in the metal transportation scheme each simulation has adopted. We
have already reported a systematic discrepancy in the metal distribution
between the two hydrodynamics approaches in the isolated galaxy comparison
(Paper II). And Shin et al. (2021) quantified the difference in metal
distribution caused by different metal diffusion schemes and different
numerical resolutions (especially in galactic halos). We will further
investigate the circumgalactic and intergalactic media of the CosmoRuns,
providing clues to the origin of the discrepancy in their metal content.

Comparing the results with previous studies, the best fit to the FIRE-2
simulations sits right in the middle of our eight simulations (thick dotted
lines; Ma et al., 2016), while that to the TNG50 simulation sits $\sim$ 1.0
dex higher in metallicity than our CosmoRuns do (thick dot-dashed lines;
Pillepich et al., 2019; Nelson et al., 2019). Additionally, we include the
observed mass$-$metallicity relation of present-day dwarf galaxies in the
Local Group (grey crosses with error bars; Sanati et al., 2023), and the
median value of the relation for 29 galaxies in a galaxy cluster field at
$z\sim 2$ (thick light blue line; JWST; Li et al., 2022).181818We adopt
$\log{(Z_{\rm gas}/Z_{\odot})}=12+\log{\rm(O/H)}-9.0$ (Ma et al., 2016). Both
of these observations show $\sim 0.5$ dex higher metallicity than the
satellite galaxies in the CosmoRuns do. Considering that the metallicity of
dwarf galaxies tends to increase from $z=2$ to 0 (by $\sim 0.4$ dex in the
FIRE-2 simulations), this difference between the CosmoRuns and Local Group
dwarfs may be less pronounced at $z=0$. The mass$-$metallicity relation seen
in Figure 10 is an important test of the realism of the feedback prescriptions
used in the CosmoRuns. While all CosmoRuns at $z\sim 2$ reproduce the stellar
masses of satellite galaxies similar to that of the MW and M31, the
differences in their metallicities indicate that the baryon physics models
implemented in the CosmoRuns have limitations.

##  5 Conclusion

We have studied the satellite halo populations near $z=2$ in the high-
resolution cosmological zoom-in simulations carried out on eight widely-used
astrophysical simulation codes (Art-I, Enzo, Ramses, Changa, Gadget-3, Gear,
Arepo-t, and Gizmo) for the AGORA High-resolution Galaxy Simulations
Comparison Project. We use different redshift epochs near $z=2$ for each code
(“$z\sim 2$”) at which the eight CosmoRuns are in the same evolutionary stage
in the target halo’s merger history, in order to alleviate the timing
discrepancy. Our key results are as follows:

  * •

All hydrodynamic CosmoRuns have fewer satellite halos than the DMO runs do at
$z\sim 2$ across all halo masses. The numbers of satellite halos in all
CosmoRuns are fewer than those in the DMO runs by a factor of $\sim\,2$ for
$M_{\rm halo}<10^{8.5}h^{-1}\,{\rm M}_{\odot}$ (Section 3.1).

  * •

The difference between CosmoRuns and DMO runs exists as early as at $z=12$.
The discrepancies in the early universe can be explained by the “gas – dark
matter particle coupling” in the particle-based codes and/or by the coarse
force resolution in the mesh-based codes in the outskirts of the target halo.
Other late-time baryonic effects such as reionization, tidal stripping, ram
pressure stripping, and stellar feedback enhance the depletion of
substructures when compared to the DMO counterparts (Sections 3.2 and 3.3).

  * •

When we consider only the halos containing stellar particles at $z\sim 2$, the
number of satellite galaxies is significantly fewer than that of dark matter
halos in all participating AGORA simulations. The populations of satellite
galaxies in all eight CosmoRuns are indeed comparable to that of present-day
satellites near the MW or M31 in their stellar masses and in their
3-dimensional stellar velocity dispersions. This finding is in line with
previous studies (Section 3.4; see also Brooks & Zolotov, 2014; Sawala et al.,
2016a; Wetzel et al., 2016; Applebaum et al., 2021).

  * •

Using the five CosmoRuns that reached $z=0.3$, we also show that the number of
satellite galaxies at $z=0.3$ are expected to be only a factor of $\lesssim 2$
larger than that at $z\sim 2$. Thus, our conclusion that the number of
satellite galaxies is significantly fewer than that of satellite halos will
likely also hold at $z\sim 0$ (Section 4.1).

  * •

We also find small, but systematic differences in other galaxy properties such
as the stellar mass$-$halo mass relation and the mass$-$metallicity relation.
Art-I, Ramses, and Gear show a relatively large $M_{\rm star}/M_{\rm halo}$
value, while the ratios for Enzo and Gizmo are slightly smaller than those of
the other codes. Similarly, Art-I and Ramses exhibit a relatively large
mass$-$metallicity relation (Section 4.2). We observe the differences in the
metallicities between CosmoRuns and the observations, which indicate that the
baryon physics models implemented in the CosmoRuns have limitations.

Overall, it is notable that the so-called “missing satellite problem” is fully
and easily resolved across all participating codes simply by implementing the
common baryonic physics adopted in AGORA and the stellar feedback prescription
commonly used in each code group, with sufficient numerical resolution
($\lesssim 100$ proper pc at $z=2$). We have demonstrated that the baryonic
solution to the decade-old problem in the $\Lambda$CDM model is effective in
all eight AGORA participating codes at $z\sim 2$. Because the results of our
numerical experiment are reproduced by one another through the AGORA
framework, the solution is independent of the numerical platform adopted —
excluding the possibility that it is an artifact of any one particular
numerical implementation. Note that the stellar feedback prescriptions in the
CosmoRun suite were calibrated to produce similar stellar masses in the host
halo by $z=4$ (see Section 5.4 in Paper III) which remains true to $z\sim 2$
(see Paper IV), but they were never specifically aimed or designed to suppress
the satellite galaxy population.

We thank all of our colleagues participating in the AGORA Project for their
collaborative spirit which has enabled the Collaboration to remain strong as a
platform to foster and launch multiple science-oriented comparison efforts. We
particularly thank Oscar Agertz, Kirk Barrow, Oliver Hahn, Desika Narayanan,
Eun-jin Shin, Britton Smith, Ben Tufeld and Matthew Turk for their insightful
comments during the work presented in this paper, and Yongseok Jo and Seungjae
Lee for their helpful feedback on the earlier version of this manuscript. We
are also grateful to Volker Springel for making the Gadget-2 code public, and
for providing the original versions of Gadget-3 to be used in the AGORA
Project. This research used resources of the National Energy Research
Scientific Computing Center, a DOE Office of Science User Facility supported
by the Office of Science of the U.S. Department of Energy under Contract No.
DE-AC02-05CH11231. Ji-hoon Kim acknowledges support by Samsung Science and
Technology Foundation under Project Number SSTF-BA1802-04. His work was also
supported by the National Research Foundation of Korea (NRF) grant funded by
the Korea government (MSIT) (No. 2022M3K3A1093827 and 2023R1A2C1003244). His
work was also supported by the National Institute of Supercomputing and
Network/Korea Institute of Science and Technology Information with
supercomputing resources including technical support, grants
KSC-2020-CRE-0219, KSC-2021-CRE-0442 and KSC-2022-CRE-0355. The publicly
available Enzo and yt codes used in this work are the products of
collaborative efforts by many independent scientists from numerous
institutions around the world. Their commitment to open science has helped
make this work possible.

##  Appendix A Halo Matching Process Between Art-I/Ramses and Art-I-
dmo/Ramses-dmo

To investigate how each individual halo is affected by the baryonic processes,
in Section 3.3 we match and compare halos in hydrodynamic simulations (e.g.,
Enzo CosmoRun) to their counterparts in the DMO simulation (e.g., Enzo-dmo
run). But because particle IDs in Ramses and Ramses-dmo (or Art-I and Art-I-
dmo) are not identically assigned in their initial conditions, we need to
employ a method different from what is described in Section 3.3 to match the
halos between these two simulations. Here we introduce an alternative approach
to find a pair of matched halos between two simulations that only share the
initial condition, but not their particle IDs. The idea is that we find a pair
of halos originating from the same dark matter patch at a nearly homogenous
early universe. First, we choose 40 particles closest to a target halo’s
center in e.g., Ramses CosmoRun at $z\sim 2$. We trace each dark matter
particle in a halo back in time, and find its position at $z=100$ (the initial
condition). Now for each of the 40 particles in the Ramses run, a particle in
the Ramses-dmo run is randomly assigned. For each pair of particles we can
compute the distances between them at $z=100$. Among the 40 edges, we sum up
the 20 smallest distances. A group of particles with the smallest distance sum
is chosen in the Ramses-dmo initial condition, and we trace them forward in
time to find a “matched” halo at $z\sim 2$. Finally, by carrying out the same
procedure in reverse, another link is obtained — i.e., first find the 40 most
bound particles in the Ramses-dmo run, and then locate their counterpart
particles in the Ramses CosmoRun. A pair of two halos that are bijectively
mapped (bidirectionally connected) in between the two simulations are
considered as a “matched” pair.

## References

  * Applebaum et al. (2021) Applebaum, E., Brooks, A. M., Christensen, C. R., et al. 2021, ApJ, 906, 96, doi: 10.3847/1538-4357/abcafa
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143, doi: 10.1093/mnras/stz1182
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, ApJ, 762, 109, doi: 10.1088/0004-637X/762/2/109
  * Benson et al. (2002) Benson, A. J., Lacey, C. G., Baugh, C. M., Cole, S., & Frenk, C. S. 2002, MNRAS, 333, 156, doi: 10.1046/j.1365-8711.2002.05387.x
  * Bode et al. (2001) Bode, P., Ostriker, J. P., & Turok, N. 2001, ApJ, 556, 93, doi: 10.1086/321541
  * Brooks et al. (2013) Brooks, A. M., Kuhlen, M., Zolotov, A., & Hooper, D. 2013, ApJ, 765, 22, doi: 10.1088/0004-637X/765/1/22
  * Brooks & Zolotov (2014) Brooks, A. M., & Zolotov, A. 2014, ApJ, 786, 87, doi: 10.1088/0004-637X/786/2/87
  * Bryan & Norman (1998) Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80, doi: 10.1086/305262
  * Bullock & Boylan-Kolchin (2017) Bullock, J. S., & Boylan-Kolchin, M. 2017, ARA&A, 55, 343, doi: 10.1146/annurev-astro-091916-055313
  * Bullock et al. (2001) Bullock, J. S., Kravtsov, A. V., & Weinberg, D. H. 2001, ApJ, 548, 33, doi: 10.1086/318681
  * Del Popolo & Le Delliou (2017) Del Popolo, A., & Le Delliou, M. 2017, Galaxies, 5, 17, doi: 10.3390/galaxies5010017
  * D’Onghia et al. (2010) D’Onghia, E., Springel, V., Hernquist, L., & Keres, D. 2010, ApJ, 709, 1138, doi: 10.1088/0004-637X/709/2/1138
  * Dunstan et al. (2011) Dunstan, R. M., Abazajian, K. N., Polisensky, E., & Ricotti, M. 2011, arXiv e-prints, arXiv:1109.6291, doi: 10.48550/arXiv.1109.6291
  * Elahi et al. (2016) Elahi, P. J., Knebe, A., Pearce, F. R., et al. 2016, MNRAS, 458, 1096, doi: 10.1093/mnras/stw338
  * Fitts et al. (2017) Fitts, A., Boylan-Kolchin, M., Elbert, O. D., et al. 2017, MNRAS, 471, 3547, doi: 10.1093/mnras/stx1757
  * Garrison-Kimmel et al. (2017) Garrison-Kimmel, S., Wetzel, A., Bullock, J. S., et al. 2017, MNRAS, 471, 1709, doi: 10.1093/mnras/stx1710
  * Garrison-Kimmel et al. (2019) Garrison-Kimmel, S., Hopkins, P. F., Wetzel, A., et al. 2019, MNRAS, 487, 1380, doi: 10.1093/mnras/stz1317
  * Girelli et al. (2020) Girelli, G., Pozzetti, L., Bolzonella, M., et al. 2020, A&A, 634, A135, doi: 10.1051/0004-6361/201936329
  * Gnedin & Hui (1998) Gnedin, N. Y., & Hui, L. 1998, MNRAS, 296, 44, doi: 10.1046/j.1365-8711.1998.01249.x
  * Grand et al. (2021) Grand, R. J. J., Marinacci, F., Pakmor, R., et al. 2021, MNRAS, 507, 4953, doi: 10.1093/mnras/stab2492
  * Green et al. (2022) Green, S. B., van den Bosch, F. C., & Jiang, F. 2022, MNRAS, 509, 2624, doi: 10.1093/mnras/stab3130
  * Gunn & Gott (1972) Gunn, J. E., & Gott, J. Richard, I. 1972, ApJ, 176, 1, doi: 10.1086/151605
  * Haardt & Madau (2012) Haardt, F., & Madau, P. 2012, ApJ, 746, 125, doi: 10.1088/0004-637X/746/2/125
  * Hahn & Abel (2011) Hahn, O., & Abel, T. 2011, MNRAS, 415, 2101, doi: 10.1111/j.1365-2966.2011.18820.x
  * Heitmann et al. (2008) Heitmann, K., Lukić, Z., Fasel, P., et al. 2008, Computational Science and Discovery, 1, 015003, doi: 10.1088/1749-4699/1/1/015003
  * Hopkins et al. (2018) Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800, doi: 10.1093/mnras/sty1690
  * Hopkins et al. (2023) Hopkins, P. F., Wetzel, A., Wheeler, C., et al. 2023, MNRAS, 519, 3154, doi: 10.1093/mnras/stac3489
  * Hu et al. (2000) Hu, W., Barkana, R., & Gruzinov, A. 2000, Phys. Rev. Lett., 85, 1158, doi: 10.1103/PhysRevLett.85.1158
  * Kauffmann et al. (1993) Kauffmann, G., White, S. D. M., & Guiderdoni, B. 1993, MNRAS, 264, 201, doi: 10.1093/mnras/264.1.201
  * Kelley et al. (2019) Kelley, T., Bullock, J. S., Garrison-Kimmel, S., et al. 2019, MNRAS, 487, 4409, doi: 10.1093/mnras/stz1553
  * Kim et al. (2014) Kim, J.-h., Abel, T., Agertz, O., et al. 2014, ApJS, 210, 14, doi: 10.1088/0067-0049/210/1/14
  * Kim et al. (2016) Kim, J.-h., Agertz, O., Teyssier, R., et al. 2016, ApJ, 833, 202, doi: 10.3847/1538-4357/833/2/202
  * Kim et al. (2018) Kim, S. Y., Peter, A. H. G., & Hargis, J. R. 2018, Phys. Rev. Lett., 121, 211302, doi: 10.1103/PhysRevLett.121.211302
  * Klypin et al. (1999) Klypin, A., Kravtsov, A. V., Valenzuela, O., & Prada, F. 1999, ApJ, 522, 82, doi: 10.1086/307643
  * Legrand et al. (2019) Legrand, L., McCracken, H. J., Davidzon, I., et al. 2019, MNRAS, 486, 5468, doi: 10.1093/mnras/stz1198
  * Li et al. (2022) Li, M., Cai, Z., Bian, F., et al. 2022, arXiv e-prints, arXiv:2211.01382.  https://arxiv.org/abs/2211.01382
  * Lovell et al. (2021) Lovell, C. C., Wilkins, S. M., Thomas, P. A., et al. 2021, MNRAS, doi: 10.1093/mnras/stab3221
  * Ma et al. (2016) Ma, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2016, MNRAS, 456, 2140, doi: 10.1093/mnras/stv2659
  * McConnachie (2012) McConnachie, A. W. 2012, AJ, 144, 4, doi: 10.1088/0004-6256/144/1/4
  * Michaux et al. (2021) Michaux, M., Hahn, O., Rampf, C., & Angulo, R. E. 2021, MNRAS, 500, 663, doi: 10.1093/mnras/staa3149
  * Moore et al. (1999) Moore, B., Ghigna, S., Governato, F., et al. 1999, ApJ, 524, L19, doi: 10.1086/312287
  * Munshi et al. (2021) Munshi, F., Brooks, A. M., Applebaum, E., et al. 2021, ApJ, 923, 35, doi: 10.3847/1538-4357/ac0db6
  * Munshi et al. (2013) Munshi, F., Governato, F., Brooks, A. M., et al. 2013, ApJ, 766, 56, doi: 10.1088/0004-637X/766/1/56
  * Nadler et al. (2021) Nadler, E. O., Banerjee, A., Adhikari, S., Mao, Y.-Y., & Wechsler, R. H. 2021, ApJ, 920, L11, doi: 10.3847/2041-8213/ac29c1
  * Nadler et al. (2020) Nadler, E. O., Wechsler, R. H., Bechtol, K., et al. 2020, ApJ, 893, 48, doi: 10.3847/1538-4357/ab846a
  * Nelson et al. (2019) Nelson, D., Pillepich, A., Springel, V., et al. 2019, MNRAS, 490, 3234, doi: 10.1093/mnras/stz2306
  * O’Leary & McQuinn (2012) O’Leary, R. M., & McQuinn, M. 2012, ApJ, 760, 4, doi: 10.1088/0004-637X/760/1/4
  * O’Shea et al. (2005) O’Shea, B. W., Nagamine, K., Springel, V., Hernquist, L., & Norman, M. L. 2005, ApJS, 160, 1, doi: 10.1086/432645
  * Pillepich et al. (2019) Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196, doi: 10.1093/mnras/stz2338
  * Qin et al. (2017) Qin, Y., Duffy, A. R., Mutch, S. J., et al. 2017, MNRAS, 467, 1678, doi: 10.1093/mnras/stx083
  * Quirk et al. (2022) Quirk, A. C. N., Guhathakurta, P., Gilbert, K. M., et al. 2022, AJ, 163, 166, doi: 10.3847/1538-3881/ac5324
  * Revaz & Jablonka (2018) Revaz, Y., & Jablonka, P. 2018, A&A, 616, A96, doi: 10.1051/0004-6361/201832669
  * Roca-Fàbrega (2023) Roca-Fàbrega, S. 2023, in prep.
  * Roca-Fàbrega et al. (2021) Roca-Fàbrega, S., Kim, J.-H., Hausammann, L., et al. 2021, ApJ, 917, 64, doi: 10.3847/1538-4357/ac088a
  * Sales et al. (2022) Sales, L. V., Wetzel, A., & Fattahi, A. 2022, Nature Astronomy, doi: 10.1038/s41550-022-01689-w
  * Samuel et al. (2020) Samuel, J., Wetzel, A., Tollerud, E., et al. 2020, MNRAS, 491, 1471, doi: 10.1093/mnras/stz3054
  * Sanati et al. (2023) Sanati, M., Jeanquartier, F., Revaz, Y., & Jablonka, P. 2023, A&A, 669, A94, doi: 10.1051/0004-6361/202244309
  * Santos-Santos et al. (2018) Santos-Santos, I. M., Di Cintio, A., Brook, C. B., et al. 2018, MNRAS, 473, 4392, doi: 10.1093/mnras/stx2660
  * Sawala et al. (2017) Sawala, T., Pihajoki, P., Johansson, P. H., et al. 2017, MNRAS, 467, 4383, doi: 10.1093/mnras/stx360
  * Sawala et al. (2015) Sawala, T., Frenk, C. S., Fattahi, A., et al. 2015, MNRAS, 448, 2941, doi: 10.1093/mnras/stu2753
  * Sawala et al. (2016a) —. 2016a, MNRAS, 457, 1931, doi: 10.1093/mnras/stw145
  * Sawala et al. (2016b) —. 2016b, MNRAS, 456, 85, doi: 10.1093/mnras/stv2597
  * Schaller et al. (2015) Schaller, M., Frenk, C. S., Bower, R. G., et al. 2015, MNRAS, 451, 1247, doi: 10.1093/mnras/stv1067
  * Shin et al. (2021) Shin, E.-J., Kim, J.-H., & Oh, B. K. 2021, ApJ, 917, 12, doi: 10.3847/1538-4357/abffd0
  * Simon (2019) Simon, J. D. 2019, ARA&A, 57, 375, doi: 10.1146/annurev-astro-091918-104453
  * Simpson et al. (2018) Simpson, C. M., Grand, R. J. J., Gómez, F. A., et al. 2018, MNRAS, 478, 548, doi: 10.1093/mnras/sty774
  * Smith et al. (2017) Smith, B. D., Bryan, G. L., Glover, S. C. O., et al. 2017, MNRAS, 466, 2217, doi: 10.1093/mnras/stw3291
  * Spergel & Steinhardt (2000) Spergel, D. N., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 84, 3760, doi: 10.1103/PhysRevLett.84.3760
  * Springel et al. (2008) Springel, V., Wang, J., Vogelsberger, M., et al. 2008, MNRAS, 391, 1685, doi: 10.1111/j.1365-2966.2008.14066.x
  * Velliscig et al. (2014) Velliscig, M., van Daalen, M. P., Schaye, J., et al. 2014, MNRAS, 442, 2641, doi: 10.1093/mnras/stu1044
  * Webb & Bovy (2020) Webb, J. J., & Bovy, J. 2020, MNRAS, 499, 116, doi: 10.1093/mnras/staa2852
  * Wetzel et al. (2016) Wetzel, A. R., Hopkins, P. F., Kim, J.-h., et al. 2016, ApJ, 827, L23, doi: 10.3847/2041-8205/827/2/L23
  * Yoshida et al. (2003) Yoshida, N., Sugiyama, N., & Hernquist, L. 2003, MNRAS, 344, 481, doi: 10.1046/j.1365-8711.2003.06829.x
  * Zhu et al. (2016) Zhu, Q., Marinacci, F., Maji, M., et al. 2016, MNRAS, 458, 1559, doi: 10.1093/mnras/stw374