# Modeling the Kinematics of Central and Satellite Galaxies Using Normalizing
Flows

K. J. Kwon Department of Physics, University of California, Santa Barbara,
Broida Hall, Santa Barbara, CA 93106, USA  ChangHoon Hahn Department of
Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ
08544, USA

###### Abstract

Galaxy clustering contains information on cosmology, galaxy evolution, and the
relationship between galaxies and their dark matter hosts. On small scales,
the detailed kinematics of galaxies within their host halos determines the
galaxy clustering. In this paper, we investigate the dependence of the central
and satellite galaxy kinematics on $\bm{\theta}$, the intrinsic host halo
properties (mass, spin, concentration), cosmology ($\Omega_{\textrm{m}}$,
$\sigma_{8}$), and baryonic feedback from active galactic nuclei and
supernovae ($A_{\textrm{AGN1}}$, $A_{\textrm{AGN2}}$, $A_{\textrm{SN1}}$,
$A_{\textrm{SN2}}$). We utilize 2,000 hydrodynamic simulations in CAMELS run
using IllustrisTNG and SIMBA galaxy formation models. Focusing on central and
satellite galaxies with $M>10^{9}M_{\ast}$, we apply neural density estimation
(NDE) with normalizing flows to estimate their $p(\Delta r\,|\,\bm{\theta})$
and $p(\Delta v\,|\,\bm{\theta})$, where $\Delta r$ and $\Delta v$ are the
magnitudes of the halo-centric spatial and velocity offsets. With NDE, we
accurately capture the dependence of galaxy kinematics on each component of
$\bm{\theta}$. For central galaxies, we identify significant spatial and
velocity biases dependent on halo mass, concentration, and spin. For satellite
distributions, we find significant deviations from an NFW profile and evidence
that they consist of distinct orbiting and infalling populations. However, we
find no significant dependence on $\bm{\theta}$ besides a weak dependence on
host halo spin. For both central and satellite galaxies, there is no
significant dependence on cosmological parameters and baryonic feedback. These
results provide key insights for improving the current halo occupation
distribution (HOD) models. This work is the first in a series that will re-
examine and develop HOD frameworks for improved modeling of galaxy clustering
at smaller scales.

Cosmology (343), Galaxies (573), Large-scale structure of the universe (902),
Galaxy dark matter halos (1880)

##  1 Introduction

Galaxy clustering is one of the most important probes for large-scale
structure (e.g., Percival et al., 2001; Peacock et al., 2001; Tegmark et al.,
2004; Cole et al., 2005; Zehavi et al., 2011; White et al., 2011; Tinker et
al., 2012; Reid et al., 2014; Alam et al., 2017; Ivanov et al., 2020;
Kobayashi et al., 2022). It encodes key cosmological information about e.g.,
the growth and expansion history of the Universe (Ivanov & Philcox 2023; see
Huterer 2023 for a recent review). It also offers a unique opportunity to
study the connection between galaxies and their host halos as well as the
complex physics of galaxy formation and evolution (e.g., Abazajian et al.
2005; Masjedi et al. 2006; Zheng et al. 2007; Zheng & Weinberg 2007; White et
al. 2007; van den Bosch et al. 2013; Guo et al. 2015a; Hahn & Villaescusa-
Navarro 2021a; Yuan et al. 2022a). Numerous galaxy redshift surveys have thus
been conducted to analyze galaxy clustering: e.g., Sloan Digital Sky Survey
(SDSS; York et al., 2000), the Baryon Oscillation Spectroscopic Survey (BOSS;
Dawson et al., 2013a), and the Extended Baryon Oscillation Spectroscopic
Survey (eBOSS; Dawson et al., 2016). Many more upcoming surveys are planned
for the near future, e.g., the Dark Energy Spectroscopic Instrument (DESI;
DESI Collaboration et al., 2022), the Subaru Prime Focus Spectrograph (PFS;
Takada et al., 2014), the ESA Euclid satellite mission (Laureijs et al.,
2011), and the Nancy Grace Roman Space Telescope (Roman; Spergel et al., 2015;
Wang et al., 2022), to probe galaxies over even larger cosmic volumes.

Simulations, used to analyze galaxy clustering, are most commonly modeled
using the halo occupation framework, where galaxies occupy dark matter halos
predicted by cosmological $N$-body simulations (Benson et al. 2000; Peacock &
Smith 2000a; Scoccimarro et al. 2001; Yang et al. 2003; Reid et al. 2012; Zhai
et al. 2019; see Wechsler & Tinker 2018 for a recent review). The Halo
Occupation Distribution (HOD) model, one of the main halo occupation models,
provides an empirical prescription of populating halos with galaxies, without
computationally modeling the physics of galaxy formation (e.g., Seljak, 2000;
Peacock & Smith, 2000b; Berlind & Weinberg, 2002; Zheng et al., 2005, 2007;
Watson et al., 2012b; Hearin & Watson, 2013). It consists of two procedures:
determining the number of galaxies within a halo (occupation number) and
assigning them positions and velocities relative to the halo. In earlier HOD
models, the occupation number of central and satellite galaxies was determined
solely by the host halo mass. More recent models, based on observational
evidence, now include additional dependence on other halo properties or their
environment (e.g., Croton et al., 2007; Zentner et al., 2014; Hearin et al.,
2016; Zehavi et al., 2018).

For the positions and velocities of galaxies, the simplest model places
central galaxies at the center of the host halos with zero relative velocity.
The positions of satellite galaxies are sampled to trace the Navarro-Frenk-
White (NFW; Navarro et al., 1996) profiles of host halos and their velocities
are sampled using the velocity dispersion profiles of host halos, derived from
solving the Jean’s equation using the assumed halo density profile. State-of-
the-art HOD models include additional flexibility through galaxy concentration
bias and velocity bias (e.g., Guo et al., 2015a, b; Hearin et al., 2016).
Concentration bias allows the concentration of satellite galaxies to deviate
from the dark matter halo (e.g., Nagai & Kravtsov, 2005; Croton et al., 2007)
and velocity bias allows for the systematic discrepancy between the velocities
of galaxies and their host halos (e.g., Berlind et al., 2003; Reid et al.,
2014; Guo et al., 2015a, b).

Despite the continual developments, there are still major open questions in
HOD modeling. For instance, there is no consensus as to how closely satellites
trace the host density profiles (e.g., Nagai & Kravtsov 2005; Tal et al. 2012;
Watson et al. 2012a; Budzynski et al. 2012; Guo et al. 2015a; McDonough &
Brainerd 2022). In addition, given the evidence of galaxy assembly bias (e.g.,
Zentner et al., 2014; Hearin et al., 2016; Zehavi et al., 2018; Xu et al.,
2021; Hadzhiyska et al., 2021) in the occupation number, it remains to be
investigated whether halo properties and environment also impact the
kinematics. Furthermore, studies have consistently demonstrated that the
presence of baryons can leave an imprint on matter or halo clustering in small
scales (e.g., Masjedi et al. 2006; White et al. 2007; van Daalen et al. 2011;
Chan et al. 2015; Hellwing et al. 2016; Chisari et al. 2019; Schneider et al.
2019; Springel et al. 2018; Foreman et al. 2020; van Daalen et al. 2020). In
particular, van Daalen et al. (2011, 2020) and Foreman et al. (2020) used
hydrodynamic simulations to show that baryonic processes such as feedback from
active galactic nuclei (AGN) or supernovae (SN) can alter the matter power
spectrum and bispectrum from what is predicted by dark-matter-only (DMO)
simulations by $>1$% at $k\gtrsim 1\ h\,{\rm Mpc}^{-1}$. Beyond the matter
distribution, it remains to be seen whether these baryonic processes also
impact the kinematics of galaxies within the halo.

Whether and to what extent the kinematics of galaxies depend on secondary halo
properties, cosmology, or baryonic processes remains to be investigated in
detail. Quantifying these dependencies will enable us to develop HOD models
that more accurately model galaxy clustering on small scales. This is
particularly timely as recent works have now firmly established the advantages
of analyzing galaxy clustering down to small scales (e.g., Verde et al., 2002;
Sefusatti & Scoccimarro, 2005; Hahn & Villaescusa-Navarro, 2021b; Massara et
al., 2021; Gualdi et al., 2021; Eickenberg et al., 2022; Wang et al., 2022;
Yuan et al., 2022b; Massara et al., 2023; Hou et al., 2023; Porth et al.,
2023). Simulation-based analyses with emulators or simulation-based inference
using the HOD model have now been successfully applied to observations (Zhai
et al., 2019; Hahn et al., 2022; Valogiannis & Dvorkin, 2022; Storey-Fisher et
al., 2022; Zhai et al., 2023; Paillas et al., 2023; Lange et al., 2023; Thiele
et al., 2023; Hahn et al., 2023a, b; Lemos et al., 2023b; Valogiannis et al.,
2023). Improved HOD models and a better understanding of the galaxy-halo
connection will enable us to more robustly exploit this additional information
on small scales.

In this paper, we investigate the dependence of central and satellite galaxy
kinematics on internal halo properties (mass, concentration, spin), cosmology
($\Omega_{\rm m}$, $\sigma_{8}$), and baryonic feedback processes by AGN and
SNe. We use galaxies from 2,000 hydrodynamic simulations of the Cosmology and
Astrophysics with Machine Learning Simulations (CAMELS; Villaescusa-Navarro et
al. 2021, 2022) suite: 1,000 simulations using the IllustrisTNG (Weinberger et
al., 2017; Pillepich et al., 2018; Nelson et al., 2019) and 1,000 from the
SIMBA (Davé et al., 2019) galaxy formation models. CAMELS simulations are
constructed using a wide range of cosmology and baryonic processes. This
provides an ideal dataset for examining the kinematic properties of galaxies
and their dependence on cosmological and baryonic feedback parameters. To
quantify the dependencies between the galaxy kinematics and the parameters, we
use neural density estimation (NDE) with normalizing flows, which can
accurately estimate complex high-dimensional conditional probability
distributions. With NDE, we quantify the impact of internal halo properties,
cosmological parameters, and baryonic feedback on galaxy kinematics.
Furthermore, we examine deviations in galaxy kinematics from the standard HOD
models.

The structure of this paper is as follows. In Section 2, we describe CAMELS
and halo catalogs that we use to extract information on galaxy kinematics and
halo properties. Section 3 summarizes how we use NDEs to estimate central and
satellite galaxy kinematics from CAMELS. We present our results in Section 4.
Based on our findings, we discuss the implications for informing and improving
HOD models in Section 5. We summarize our analysis in Section 6.

##  2 Data

In this section, we briefly describe CAMELS and the halo catalogs constructed
from CAMELS that we use to extract the galaxy kinematics and the intrinsic
properties of the host halos.

###  2.1 CAMELS

CAMELS is a suite of cosmological simulations designed for machine-learning-
based cosmological analyses. It consists of 4,033 simulations, each with a
comoving volume of (25 $h^{-1}{\rm Mpc}$)3. 2,184 of these are the
(magneto-)hydrodynamic simulations that implement the subgrid physics models
of IllustrisTNG or SIMBA. Each hydrodynamic simulation evolves $256^{3}$ dark
matter particles of mass $6.49\times 10^{7}(\Omega_{\rm m}-\Omega_{\rm
b})/0.251\ h^{-1}M_{\odot}$ with $256^{3}$ fluid elements of mass $1.27\times
10^{7}\ h^{-1}M_{\odot}$ from $z=127$ to $z=0$ for given values of the
cosmological and hydrodynamic parameters.

For our study, we use the 2,000 hydrodynamic simulations at $z=0$ constructed
using the IllustrisTNG and SIMBA models (1,000 each). These simulations vary
two cosmological parameters, $\Omega_{\rm m}$ and $\sigma_{8}$, uniformly
across the range $[0.1,0.5]$ and $[0.6,1.0]$, respectively. The other
cosmological parameters are fixed to $\Omega_{\rm b}=0.049$, $h=0.6711$,
$n_{s}=0.9624$, $M_{\nu}=0.0$ eV, $w=-1$, and $\Omega_{K}=0$. The simulations
also vary four hydrodynamic parameters, $A_{\rm hydro}=\\{A_{\rm SN1},A_{\rm
SN2},A_{\rm AGN1},A_{\rm AGN2}\\}$, that control the energy scale and the
injection rate of SN and AGN feedback. $A_{\rm SN1}$ and $A_{\rm SN2}$ provide
normalization factors for the galactic wind flux and speed, respectively.
$A_{\rm AGN1}$ and $A_{\rm AGN2}$ control the energy output normalization and
specific energy of AGN feedback. They are all log-uniformly sampled across the
range $A_{\rm SN1}$, $A_{\rm AGN1}\in[0.25,4.0]$ and $A_{\rm SN2}$, $A_{\rm
AGN2}\in[0.5,2.0]$. We note that although they share the $A_{\rm hydro}$
parameterization, each subgrid physics model implements $A_{\rm hydro}$
differently (see Villaescusa-Navarro et al., 2021, for more details). In
total, six parameters—$\Omega_{\rm m}$, $\sigma_{8}$, $A_{\rm SN1},A_{\rm
SN2},A_{\rm AGN1},A_{\rm AGN2}$—of our CAMELS simulations are varied and
arranged in a Latin Hypercube (LH).

###  2.2 Rockstar Halos

For each snapshot of the simulations, CAMELS provides multiple halo catalogs
constructed using different halo-finding algorithms. In this work, we use halo
catalogs generated by the Robust Overdensity Calculation using K-Space
Topologically Adaptive Refinement (Rockstar; Behroozi et al., 2013), a 6D
phase space temporal halo-finding method. Rockstar implements a variant of the
3D friends-of-friends (FoF) algorithm to find overdense regions. Within each
FoF group, it identifies the hierarchy of FoF subgroups by adaptively
establishing a 6D phase-space metric at each depth. This process terminates
when the number of particles in a subgroup falls below a minimum threshold.
Once all 6D FoF subgroups have been identified, a seed halo is assigned to
each subgroup at the deepest level of the hierarchy. Subsequently, particle
membership is determined based on the phase-space distance to seed halos to
determine the halo substructures. Rockstar then calculates subhalo and halo
properties using the particle membership. In this work, we only include host
halos and subhalos from IllustrisTNG and SIMBA suites that host galaxies with
stellar mass, $M_{*}$, greater than $10^{9}M_{\odot}$. This is to focus on
intermediate- to high-mass galaxies that are typically probed by cosmological
galaxy surveys (e.g., Eisenstein et al., 2001; Dawson et al., 2013b; Hahn et
al., 2023c; Zhou et al., 2023). We also remove a small fraction ($\ll 1\%$) of
central galaxies that are identified outside of their host halos. These are
mostly galaxies that have recently experienced major mergers, for which
determining the host halos is ambiguous.

From the galaxy/subhalo properties, we derive the kinematics of central and
satellite galaxies using their halo-centric radial position and velocity. For
central galaxies, we use $X_{\rm off}$ and $V_{\rm off}$ of host halos
provided by Rockstar. $X_{\rm off}$ is the comoving distance between the
density peak, the position of the particle with the minimum gravitational
potential, and the center of mass of a halo, calculated using all particles
within the virial radius, $R_{\rm vir}$111Rockstar determines the halo center
using $n$ innermost particles such that it minimize the Poisson error of the
center.. $V_{\rm off}$ is the magnitude of the velocity offset between the
bulk (all particles within $R_{\rm vir}$) and the core (particles within
0.1$R_{\rm vir}$) of the halo. For satellite galaxies, we derive $\Delta r$ by
taking the magnitude of the difference between the comoving position vectors
of its subhalo and its host halo, $|\bm{r}_{\rm sat}-\bm{r}_{\rm host}|$.
Similarly, $\Delta v$ is taken to be the magnitude of the difference between
the bulk velocities of the satellite and host halos, $|\bm{v}_{\rm
sat}-\bm{v}_{\rm host}|$.

In addition to the kinematic properties, we also use three internal properties
of host halos, which we collectively refer to as $\bm{\theta}_{\rm h}$. We
include halo virial mass ($M_{\rm vir}$), spin ($\lambda$), and concentration
($c$). $\lambda$ is given by $J\sqrt{|E|}/GM_{\rm vir}^{2.5}$, where $J$ is
the magnitude of the halo angular momentum and $|E|$ is the magnitude of the
total energy of the halo (Peebles, 1969). We derive $c$, which has been shown
to be a proxy of halo assembly history (Castorina & Sheth, 2013; Mao et al.,
2015), as the ratio of the fitted NFW scale radius and $R_{\rm vir}$. In
Figure 1, we present the distribution of $\bm{\theta}_{\rm h}$ for the central
(blue) and satellite galaxies (gray) in SIMBA. The contours are smoothed using
a Gaussian kernel with 1.5$\sigma$. We only show the $\bm{\theta}_{\rm h}$
distribution of SIMBA galaxies for clarity, as IllustrisTNG galaxies have a
similar distribution.

In total, our dataset includes 359,494 central galaxies and 22,244 satellites
from SIMBA, and 276,527 central galaxies and 34,695 satellites from
IllustrisTNG. We randomly select 10% of the galaxies and reserve them for
testing; we use the remaining galaxies for training222Each of 2,000
simulations is not specifically set aside for training or testing; i.e.,
different galaxies from the same simulation can be used in training,
validation, and testing. For each galaxy $i$, we compile $(\Delta r_{i},\Delta
v_{i},\bm{\theta}_{{\rm h},i},\Omega_{{\rm m},i},\sigma_{8,i},A_{{\rm
hydro},i})$.

Figure 1: The distribution of host halo properties, $M_{\rm vir},\lambda,c$,
for central galaxies (blue) and satellite galaxies (gray) in the SIMBA
hydrodynamic suite. The contours are marked for 68th, 95th, and 99th
percentiles. In total, we have 359,494 central galaxies and 22,244 satellite
galaxies with $M_{*}>10^{9}M_{\odot}$ in SIMBA. IllustrisTNG has a similar
distribution, with 276,527 central and 34,695 satellite galaxies. The black
dashed lines mark the fiducial values of halo properties. We examine the
distribution of galaxy kinematics at the points marked in orange (‘x’) to
identify the dependence on each halo property.

##  3 Neural Density Estimation with Normalizing Flows

Our objective is to identify how the host halo properties, cosmological
parameters, and baryonic processes affect $\Delta r$ and $\Delta v$ of
galaxies within a halo. We also seek to quantify any deviations in the $\Delta
r$ and $\Delta v$ distributions from the standard HOD prescription. To do
this, we estimate $p(\Delta r\,|\,\bm{\theta})$ and $p(\Delta
v\,|\,\bm{\theta})$, where $\bm{\theta}=\\{\bm{\theta}_{\rm h},\Omega_{\rm
m},\sigma_{8},A_{\rm hydro}\\}$ for CAMELS using neural density estimation
(NDE), which has been used extensively in the literature to accurately
estimate complex high-dimensional distributions (e.g., Alsing et al., 2019;
Hassan et al., 2022; Ruhe et al., 2022; Hahn & Melchior, 2022; Dai & Seljak,
2022; Lovell et al., 2023).

With NDE, we train neural networks ($q$) with free parameters $\bm{\phi}$ such
that $q_{\bm{\phi}}(\Delta a\,|\,\bm{\theta})\approx p(\Delta
a\,|\,\bm{\theta})$, where $\Delta a$ refers to either $\Delta r$ and $\Delta
v$. Once trained, $q_{\bm{\phi}}(\Delta a\,|\,\bm{\theta})$ can be evaluated
and sampled at different values of $\bm{\theta}$. This means that we can
isolate the effect of a single parameter by comparing $q_{\bm{\phi}}(\Delta
a\,|\,\bm{\theta})$ to $q_{\bm{\phi}}(\Delta a\,|\,\bm{\theta}^{\prime})$,
where we only vary a single parameter in $\bm{\theta}$ to obtain
$\bm{\theta}^{\prime}$. In this work, we use NDE based on normalizing flows
(Tabak & Vanden-Eijnden, 2010; Tabak & Turner, 2013), which estimates a target
distribution by mapping it to a base distribution that is simple and easy to
evaluate (e.g., a multivariate Gaussian distribution). The mapping between the
target and base distributions is set to a series of invertible and
differentiable transformations. This enables us to calculate and sample the
target distribution using the base distribution with a change of variables.

Among various normalizing flow architectures, e.g., Masked Autoregressive Flow
(Papamakarios et al., 2017), Neural Spline Flow (Durkan et al., 2019), we use
Masked Autoencoder for Distribution Estimation (MADE; Germain et al., 2015;
Durkan et al., 2020), based on experimentation. MADE decomposes
$p(\bm{x}\,|\,\bm{\theta})$ into a product of nested conditional
probabilities,
$\prod_{i=1}^{\text{dim}(\bm{x})}p(x_{i}\,|\,\bm{x}_{<i},\bm{\theta})$, where
$\bm{x}_{<i}=(x_{1},x_{2},\cdots,x_{i-1})$. Given the ordering of input
$\bm{x}$, an autoregressive model parameterizes each
$p(x_{i}\,|\,\bm{x}_{<i},\bm{\theta})$ with $\bm{x}_{<i}$ and $\bm{\theta}$ as
the only input. MADE preserves the autoregressive property by using binary
masks to drop out connections between $p(x_{i}\,|\,\bm{x}_{<i},\bm{\theta})$
and $\bm{x}_{\geq i}$. We use the MADE implementation from the
sbi333https://github.com/sbi-dev/sbi python package (Greenberg et al., 2019;
Tejero-Cantero et al., 2020).

To train our flow, we first divide our dataset $\\{(\Delta
a_{i},\bm{\theta}_{i})\\}$ from Section 2 into a training and validation set
with a 90/10 split. Then, we maximize the total log-likelihood
$\sum_{i=1}^{n}\log q_{\bm{\phi}}(\Delta a_{i}\,|\,\bm{\theta}_{i})$, for
galaxies in our training set. This is equivalent to minimizing the Kullback-
Leibler divergence between $q_{\bm{\phi}}(\Delta a\,|\,\bm{\theta})$ and
$p(\Delta a\,|\,\bm{\theta})$. For the maximization, we use the Adam optimizer
(Kingma & Ba, 2017) with a learning rate of $5\times 10^{-4}$. We impose early
stopping to prevent overfitting by stopping training when the log-likelihood
evaluated on the validation set fails to increase after 20 consecutive epochs.
We determine the number of hidden layers and units of $q_{\bm{\phi}}$ through
experimentation to achieve the lowest validation losses. In total, we train 8
separate flows to estimate $p(\Delta r\,|\,\bm{\theta})$ and $p(\Delta
v\,|\,\bm{\theta})$444In practice, we estimate $p(\log\Delta
r\,|\,\bm{\theta})$ and $p(\log\Delta v\,|\,\bm{\theta})$ to reduce the
dynamic range of the kinematics data. However, we find that the choice between
training NDEs on linear or logarithmic scales does not have a significant
impact on their accuracy., of central and satellite galaxies, and for
IllustrisTNG and SIMBA. In Appendix A, we validate the accuracy of our trained
NDEs using the 10% of the CAMELS galaxies reserved for testing.

##  4 Results

In this section, we present the dependence of the kinematics of central and
satellite galaxies on their host halo properties, cosmology, and baryonic
processes using $q_{\bm{\phi}}(\Delta r\,|\,\bm{\theta})$ and
$q_{\bm{\phi}}(\Delta v\,|\,\bm{\theta})$. To focus on the effect of
individual parameters, we first select a fiducial $\bm{\theta}_{\rm fid}$:

$\bm{\theta}_{\rm fid}=(\log M_{\rm vir}=13.5,\log
c=0.77,\log\lambda=-1.34,\Omega_{\rm m}=0.3,\sigma_{8}=0.8,A_{\rm hydro}=1),$
(1)

where $M_{\rm vir}$ is in units of $h^{-1}M_{\odot}$. We then compare
$q_{\bm{\phi}}(\Delta r\,|\,\bm{\theta}_{\rm fid})$ to $q_{\bm{\phi}}(\Delta
r\,|\,\bm{\theta}^{\prime})$, where $\bm{\theta}^{\prime}$ has the same
parameter values as $\bm{\theta}_{\rm fid}$, except for a single parameter.

We set the fiducial $M_{\rm vir}$ based on the typical host halo mass of
luminous red galaxies (LRG; e.g., Zheng et al., 2009). We choose the fiducial
concentration near the center of the training data (Figure 1) and the fiducial
$\lambda$ based on the typical spin of halos in numerical simulations
(Steinmetz & Bartelmann, 1995). The fiducial $\Omega_{\rm m}$ and $\sigma_{8}$
values are at the center of the CAMELS prior range, which is consistent with
Planck Collaboration et al. (2020) constraints. The fiducial $A_{\rm hydro}$
are set to unity, which corresponds to the original feedback settings of
IllustrisTNG and SIMBA. We note that although we do not include $R_{\rm vir}$
in $\bm{\theta}_{\rm h}$, it varies based on $M_{\rm vir}$ and $\Omega_{\rm
m}$555$M_{\rm vir}=4\pi\Delta(x)\,\rho_{\rm crit}/3R_{\rm vir}^{3}$, where
$x=1-\Omega_{\rm m}$ and $\Delta(x)=18\pi^{2}+82x-39x^{2}$ (Bryan & Norman,
1998).. In this work, we present results at $\bm{\theta}_{\rm fid}$ specified
in Equation 1. However, our results do not vary significantly for different
$\bm{\theta}_{\rm fid}$.

Figure 2: Dependence of the central galaxy kinematics, $q_{\bm{\phi}}(\Delta
r_{\rm cen}/R_{\rm vir}\,|\,\bm{\theta})$ and $q_{\bm{\phi}}(\Delta v_{\rm
cen}/V_{\rm vir}\,|\,\bm{\theta})$, on $\bm{\theta}_{\rm h}$ ($M_{\rm vir}$,
$c$, and $\lambda$, from top to bottom) in IllustrisTNG (red) and SIMBA
(blue). The $y$-axis denotes the probability density on a linear scale. The
parameters that do not appear in the legends are kept at the fiducial values.
The dashed lines mark the medians of the distributions. We find a significant
dependence of $q_{\bm{\phi}}(\Delta r_{\rm cen}/R_{\rm vir}\,|\,\bm{\theta})$
not only on $M_{\rm vir}$ but also on $c$ and $\lambda$ for both IllustrisTNG
and SIMBA: $q_{\bm{\phi}}(\Delta r_{\rm cen}/R_{\rm vir}\,|\,\bm{\theta})$ is
skewed towards lower $\Delta r_{\rm cen}/R_{\rm vir}$ for higher $M_{\rm vir}$
and $c$ while the trend is reversed for $\lambda$. The trends in both
simulation suites are consistent with each other; however,
$q_{\bm{\phi}}(\Delta r_{\rm cen}/R_{\rm vir}\,|\,\bm{\theta})$ is broader in
SIMBA than in IllustrisTNG.

###  4.1 Central Galaxies

_Host halo properties_ : In Figure 2, we present the dependence of central
galaxy kinematics, $q_{\bm{\phi}}(\Delta r_{\rm cen}/R_{\rm
vir}\,|\,\bm{\theta})$ (a) and $q_{\bm{\phi}}(\Delta v_{\rm cen}/V_{\rm
vir}\,|\,\bm{\theta})$ (b), on $\log M_{\rm vir}$ (top), $\log c$ (center),
and $\log\lambda$ (bottom) for IllustrisTNG (red) and SIMBA (blue). For
clarity we present $\Delta r$ and $\Delta v$ in units of $R_{\rm vir}$ and
$V_{\rm vir}$. The value of the varied parameter increases from light to dark.
The dashed vertical lines indicate the medians of the distributions to
highlight any trends. Overall, we find that the dependence of the central
kinematics on each halo property is qualitatively consistent between
IllustrisTNG and SIMBA. Furthermore, we find that central galaxies are
radially displaced from the center of the halos and have non-zero relative
velocities. This validates the use of velocity bias of central galaxies in
current state-of-the-art HOD models.

Figure 2 also confirms a number of expected dependencies of $\Delta r_{\rm
cen}$ and $\Delta v_{\rm cen}$ with halo properties. For instance, the medians
of $q_{\bm{\phi}}(\Delta r_{\rm cen}/R_{\rm vir}\,|\,\bm{\theta})$ and
$q_{\bm{\phi}}(\Delta v_{\rm cen}/V_{\rm vir}\,|\,\bm{\theta})$ decrease as
$M_{\rm vir}$ or $c$ increase. The gravitational potential is deeper for
higher $M_{\rm vir}$ and steeper towards the center for higher $c$. In
contrast, the medians of the distributions increase with increasing $\lambda$.
At fixed fiducial $M_{\rm vir}$ and $R_{\rm vir}$, increasing $\lambda$ is
equivalent to increasing $J$, i.e., increased rotational support against the
gravitational infall, which would increase $\Delta r_{\rm cen}$ and $\Delta
v_{\rm cen}$. Furthermore, $q_{\bm{\phi}}(\Delta r_{\rm cen}/R_{\rm
vir}\,|\,\bm{\theta})$ is slightly broader in SIMBA than IllustrisTNG, most
noticeably in the $\log c$ and $\log\lambda$ panels.

_Cosmological parameters_ : In Figure 3, we present the dependence of central
galaxy kinematics, $q_{\bm{\phi}}(\Delta r_{\rm cen}/R_{\rm
vir}\,|\,\bm{\theta})$ (a) and $q_{\bm{\phi}}(\Delta v_{\rm cen}/V_{\rm
vir}\,|\,\bm{\theta})$ (b), on $\Omega_{\rm m}$ (top) and $\sigma_{8}$
(bottom) for IllustrisTNG (red) and SIMBA (blue). Overall, the cosmological
parameters have a negligible effect on the central galaxy kinematics. For both
IllustrisTNG and SIMBA, there is a marginal dependence of $\Omega_{\rm m}$ on
$q_{\bm{\phi}}(\Delta r_{\rm cen}/R_{\rm vir}\,|\,\bm{\theta})$, where central
galaxies are slightly farther away from the halo centers at higher
$\Omega_{\rm m}$. This trend remains even when we vary halo mass with
$\Omega_{\rm m}$, based on the fixed number density of the halo mass function
(Appendix B). Nevertheless, the effect of $\Omega_{\rm m}$ on
$q_{\bm{\phi}}(\Delta r_{\rm cen}/R_{\rm vir}\,|\,\bm{\theta})$ and
$q_{\bm{\phi}}(\Delta v_{\rm cen}/V_{\rm vir}\,|\,\bm{\theta})$ is small.

Figure 3: Dependence of the central galaxy kinematics, $q_{\bm{\phi}}(\Delta
r_{\rm cen}/R_{\rm vir}\,|\,\bm{\theta})$ and $q_{\bm{\phi}}(\Delta v_{\rm
cen}/V_{\rm vir}\,|\,\bm{\theta})$, on $\Omega_{\rm m}$ (top) and $\sigma_{8}$
(bottom) in IllustrisTNG (red) and SIMBA (blue). The parameters that do not
appear in the legends are kept at the fiducial values. The dashed lines mark
the medians of the distributions. In both IllustrisTNG and SIMBA,
$q_{\bm{\phi}}(\Delta r_{\rm cen}/R_{\rm vir}\,|\,\bm{\theta})$ is slightly
skewed to higher $\Delta r_{\rm cen}/R_{\rm vir}$ at higher $\Omega_{\rm m}$.
Overall, we find no significant dependence of galaxy kinematics on the
cosmological parameters.

_Baryonic feedback_ : In Figure 4, we present the dependence of central galaxy
kinematics, $q_{\bm{\phi}}(\Delta r_{\rm cen}/R_{\rm vir}\,|\,\bm{\theta})$
(a) and $q_{\bm{\phi}}(\Delta v_{\rm cen}/V_{\rm vir}\,|\,\bm{\theta})$ (b),
on $A_{\rm hydro}$ for IllustrisTNG (left) and SIMBA (right). In each panel,
we include the fiducial distributions (gray-shaded) for comparison. For the
other distributions, we vary a single feedback parameter: $A_{\rm SN1}=0.5$
(blue), $A_{\rm AGN1}=0.5$ (orange), $A_{\rm SN2}=1.5$ (green), and $A_{\rm
AGN2}=1.5$ (red). We set the parameters to either 0.5 or 1.5 to ensure that we
sufficiently vary baryonic feedback. For reference, in SIMBA, setting the
parameters to 0.5 or 1.5 modifies the stellar mass function (SMF) by $\gtrsim
25\%$ at $10^{11}M_{\odot}$. Despite the significant impact they have on the
galaxy population, the baryonic feedback parameters overall do not
significantly impact the kinematics of central galaxies.

Figure 4: Dependence of the central galaxy kinematics, $q_{\bm{\phi}}(\Delta
r_{\rm cen}/R_{\rm vir}\,|\,\bm{\theta})$ and $q_{\bm{\phi}}(\Delta v_{\rm
cen}/V_{\rm vir}\,|\,\bm{\theta})$, on the baryonic feedback parameters. The
distributions vary a single feedback parameter as indicated in the legend
(blue, orange, green, and red), while the rest are kept at the fiducial
values. We include the distributions at $\bm{\theta}_{\rm fid}$ for reference
(gray). The dashed lines mark the medians of the distributions. The baryonic
feedback parameters do not significantly impact central galaxy kinematics.
Figure 5: Dependence of the satellite galaxy kinematics, $q_{\bm{\phi}}(\Delta
r_{\rm sat}/R_{\rm vir}\,|\,\bm{\theta})$ and $q_{\bm{\phi}}(\Delta v_{\rm
sat}/V_{\rm vir}\,|\,\bm{\theta})$, on $\bm{\theta}_{\rm h}$ ($M_{\rm vir}$,
$c$, and $\lambda$, from top to bottom) in IllustrisTNG (red) and SIMBA
(blue). The parameters that do not appear in the legends are kept at the
fiducial values. For reference, we include NFW profiles using the three values
of $c$ (dashed; $\log c$ panel). $q_{\bm{\phi}}(\Delta r_{\rm sat}/R_{\rm
vir}\,|\,\bm{\theta})$ differs between IllustrisTNG and SIMBA. They also
differ from the NFW profiles (Section 5.2). In addition, $q_{\bm{\phi}}(\Delta
v_{\rm sat}/V_{\rm vir}\,|\,\bm{\theta})$ exhibits bimodality, which is due to
infalling and orbiting populations of satellite galaxies (Appendix C).
Overall, the effect of $\bm{\theta}_{\rm h}$ on the satellite kinematics is
not significant, except for the effect of $\lambda$ on $\Delta v_{\rm
sat}/V_{\rm vir}$, which changes the infalling fraction.

###  4.2 Satellite Galaxy Kinematics

_Host halo properties_ : In Figure 5, we present the dependence of satellite
galaxy kinematics, $q_{\bm{\phi}}(\Delta r_{\rm sat}/R_{\rm
vir}\,|\,\bm{\theta})$ (a) and $q_{\bm{\phi}}(\Delta v_{\rm sat}/V_{\rm
vir}\,|\,\bm{\theta})$ (b), on $\log M_{\rm vir}$ (top), $\log c$ (center),
and $\log\lambda$ (bottom) for IllustrisTNG (red) and SIMBA (blue). We note
that for $q_{\bm{\phi}}(\Delta r_{\rm sat}/R_{\rm vir}\,|\,\bm{\theta})$, we
impose an upper limit of $\Delta r_{\rm sat}/R_{\rm vir}=1$, since Rockstar
only identifies satellites within $R_{\rm vir}$. This only discards a small
fraction of satellites and does not affect our overall results. We find
significant discrepancies in $q_{\bm{\phi}}(\Delta r_{\rm sat}/R_{\rm
vir}\,|\,\bm{\theta})$ between IllustrisTNG and SIMBA. Compared to
IllustrisTNG, the satellite distribution in SIMBA is skewed toward $R_{\rm
vir}$. Despite the discrepancy, we find no strong dependence of
$q_{\bm{\phi}}(\Delta r_{\rm sat}/R_{\rm vir}\,|\,\bm{\theta})$ on any
$\bm{\theta}_{\rm h}$ for either subgrid model. The lack of dependence on
$M_{\rm vir}$ is consistent with Hadzhiyska et al. (2023), where they examined
the radial distribution of satellite subhalos in MilleniumTNG (MTNG; Pakmor et
al., 2023) as a function of host halo mass.

Beyond the lack of $\bm{\theta}_{\rm h}$ dependence, we find that
$q_{\bm{\phi}}(\Delta r_{\rm sat}/R_{\rm vir}\,|\,\bm{\theta})$ deviates
significantly from the NFW profile (dashed). This can be driven by numerical
effects that unphysically disrupt subhalos and make them unidentifiable by
halo-finding algorithms at the inner radii of host halos (e.g., van Kampen,
1995, 2000; Moore et al., 1996; van Kampen et al., 1999; Klypin et al., 1999;
van den Bosch, 2017; van den Bosch et al., 2018; Diemer et al., 2023). We
discuss this in more detail in Section 5.2.

For $q_{\bm{\phi}}(\Delta v_{\rm sat}/V_{\rm vir}\,|\,\bm{\theta})$, we first
note that the distribution is significantly bimodal. This bimodality comes
from separate orbiting and infalling populations of satellites and the latter
approximately accounts for $10\%$ of our satellites in IllustrisTNG and SIMBA
(see Appendix C). The infalling populations are distributed around $\Delta
v_{\rm sat}/V_{\rm vir}\approx 1.5$. In addition, the bimodality becomes
slightly more pronounced with higher $M_{\rm vir}$ in SIMBA, while the effect
of $M_{\rm vir}$ is negligible in IllustrisTNG. There is also a slight effect
of $\lambda$ on $q_{\bm{\phi}}(\Delta v_{\rm sat}\,|\,\bm{\theta})$, where the
bimodality becomes stronger with higher $\lambda$ in both IllustrisTNG and
SIMBA.

Figure 6: Dependence of the satellite galaxy kinematics, $q_{\bm{\phi}}(\Delta
r_{\rm sat}/R_{\rm vir}\,|\,\bm{\theta})$ and $q_{\bm{\phi}}(\Delta v_{\rm
sat}/V_{\rm vir}\,|\,\bm{\theta})$, on $\Omega_{\rm m}$ (top) and $\sigma_{8}$
(bottom) in IllustrisTNG (red) and SIMBA (blue). The parameters that do not
appear in the legends are kept at the fiducial values. The effects of the
cosmological parameters are small in both hydrodynamic suites.

_Cosmological parameters_ : In Figure 6, we present the dependence of
satellite galaxy kinematics, $q_{\bm{\phi}}(\Delta r_{\rm sat}/R_{\rm
vir}\,|\,\bm{\theta})$ (a) and $q_{\bm{\phi}}(\Delta v_{\rm sat}/V_{\rm
vir}\,|\,\bm{\theta})$ (b), on $\Omega_{\rm m}$ (top) and $\sigma_{8}$
(bottom) for IllustrisTNG (red) and SIMBA (blue). For $q_{\bm{\phi}}(\Delta
r_{\rm sat}/R_{\rm vir}\,|\,\bm{\theta})$, we find no significant dependence
on $\Omega_{\rm m}$ and $\sigma_{8}$. For $q_{\bm{\phi}}(\Delta v_{\rm
sat}/V_{\rm vir}\,|\,\bm{\theta})$, the distribution slightly shifts to lower
$\Delta v_{\rm sat}$ with higher $\Omega_{\rm m}$ in both IllustrisTNG and
SIMBA. Additionally, exclusively in IllustrisTNG, the velocity distribution
shifts to lower $\Delta v_{\rm sat}$ with higher $\sigma_{8}$. In Appendix B,
we vary $M_{\rm vir}$ with $\Omega_{\rm m}$ using the fixed number density in
the halo mass function. We find a similar weak dependence. Overall, varying
$\Omega_{\rm m}$ and $\sigma_{8}$ does not significantly affect the satellite
kinematics as well.

Figure 7: Dependence of the satellite galaxy kinematics, $q_{\bm{\phi}}(\Delta
r_{\rm sat}/R_{\rm vir}\,|\,\bm{\theta})$ and $q_{\bm{\phi}}(\Delta v_{\rm
sat}/V_{\rm vir}\,|\,\bm{\theta})$, on the baryonic feedback parameters. The
distributions vary a single feedback parameter as indicated in the legend
(blue, orange, green, and red), while the rest are kept at the fiducial
values. We include the distributions at $\bm{\theta}_{\rm fid}$ for reference
(gray). The effect of baryonic feedback is negligible for satellite galaxy
kinematics.

_Baryonic feedback_ : In Figure 7, we present the dependence of satellite
galaxy kinematics, $q_{\bm{\phi}}(\Delta r_{\rm sat}/R_{\rm
vir}\,|\,\bm{\theta})$ (a) and $q_{\bm{\phi}}(\Delta v_{\rm sat}/V_{\rm
vir}\,|\,\bm{\theta})$ (b), on $A_{\rm hydro}$ for IllustrisTNG (left) and
SIMBA (right). We find no significant dependence of $q_{\bm{\phi}}(\Delta
r_{\rm sat}/R_{\rm vir}\,|\,\bm{\theta})$ or $q_{\bm{\phi}}(\Delta v_{\rm
sat}/V_{\rm vir}\,|\,\bm{\theta})$ on the change of the baryonic feedback
parameters. This is in spite of the fact that we vary $A_{\rm hydro}$
sufficiently to modify the SMF (Section 4.1). While studies find that baryonic
processes can affect small-scale matter clustering, we find that they alone do
not significantly impact the kinematics of satellite galaxies.

##  5 Discussion

###  5.1 Galaxy Kinematics

Galaxy kinematics estimated by our NDEs have several implications for HOD
models. First, central galaxies are not aligned with the centers of their host
halos with zero halo-centric velocity. Instead, on average, they are located
at $\Delta r_{\rm cen}\sim 0.03R_{\rm vir}$ with $\Delta v_{\rm cen}\sim
0.1V_{\rm vir}$666This spatial bias is not due to gravitational softening, as
the halo density profile at $r=0.03R_{\rm vir}$ is weakly affected by the
choice of the softening length (Power et al., 2003; Zhang et al., 2019).. Our
spatial and velocity biases of central galaxies are consistent with
observational evidence. van den Bosch et al. (2005) found that brightest halo
galaxies in the Two-Degree Field Galaxy Redshift Survey (2dFGRS; Cole et al.,
2005) and SDSS have non-zero specific kinetic energy in the coordinate frame
defined with respect to the satellites and the velocity bias is larger in more
massive halos. Guo et al. (2015b) reached similar conclusions using SDSS-III
BOSS central galaxies. Additionally, Zitrin et al. (2012) showed that bright
cluster galaxies in 10,000 SDSS clusters have a spatial bias from the halo
center. In addition to spatial and velocity biases, we also find that the
kinematics of central galaxies depend significantly on $c$ and $\lambda$.
Incorporating these biases and dependencies into HOD models will enable more
accurate modeling of central galaxy kinematics.

Beyond these effects, we find little evidence that there are additional
dependencies that need to be included in HOD models for central galaxies. In
particular, we find that AGN and SN feedback do not significantly impact the
central galaxy kinematics. This is consistent with previous works. For
instance, Hellwing et al. (2016) matched halos in a hydrodynamic simulation
and its dark matter only counterpart in the Evolution and Assembly of GaLaxies
and their Environments (EAGLE; Schaye et al., 2015) simulation. They found
peculiar velocity offsets between galaxies are consistent with zero.

For satellite galaxies, we find that the radial distribution differs
significantly between IllustrisTNG and SIMBA. The distributions of both models
also deviate from the NFW profile, with distributions skewed towards higher
$\Delta r_{\rm sat}/R_{\rm vir}$ than NFW. Similar studies by Yuan et al.
(2022c) and Hadzhiyska et al. (2023) find LRG satellites in IllustrisTNG and
MTNG, selected using sliding color-magnitude cuts and stellar mass cuts
respectively, have consistent shapes of radial distributions as our
$q_{\bm{\phi}}(\Delta r_{\rm sat}/R_{\rm vir}\,|\,\bm{\theta})$. In this work,
we confirm that the deviation from NFW is not unique to IllustrisTNG but is
also evident in SIMBA. However, as mentioned above, the deviation from NFW may
be driven by known numerical limitation of $N$-body simulations, where
satellite subhalos near the center of host halos can experience unphysical
mass stripping and disruption. We discuss this limitation further in the
following section.

In addition, we find signatures of accreted, infalling satellite galaxies in
$q_{\bm{\phi}}(\Delta v_{\rm sat}/V_{\rm vir}\,|\,\bm{\theta})$ (see also
Appendix C). The velocity distribution of the satellite galaxies exhibits
bimodality due to the coexistence of infalling and orbiting populations. We
estimate in Appendix C that the infalling population accounts for $\sim 10\%$
of satellite galaxies with $M_{\ast}>10^{9}M_{\odot}$. The distinct
populations of infalling and orbiting satellites pose the question of whether
the current HOD models should explicitly include an infalling population.
Since their kinematics are distinct from orbiting satellites, the infalling
population can significantly impact small-scale galaxy clustering, especially
in redshift space. Furthermore, including infalling satellites may enable us
to probe galaxy formation and galaxy-halo connection (e.g., Mihos & Hernquist,
1994; Gill et al., 2005; More et al., 2011; Rocha et al., 2012; Lange et al.,
2019; Wright et al., 2022) or test general relativity (e.g., Zu & Weinberg,
2013; Zu et al., 2014). If one were to include infalling satellites in HOD
modeling, our NDEs would provide a direct way of including them by using the
dynamics of simulations.

Aside from the shape of the satellite galaxy kinematics distributions, we do
not find any significant parameter dependence on $\bm{\theta}_{\rm h}$,
$\Omega_{\rm m}$, $\sigma_{8}$, or $A_{\rm hydro}$. We find a dependence of
$q_{\bm{\phi}}(\Delta v_{\rm sat}/V_{\rm vir}\,|\,\bm{\theta})$ on $\lambda$;
however, the effect is small and the impact is primarily on the relative
amplitude between the orbiting and infalling fraction. Interestingly, baryonic
feedback does not significantly impact satellite positions or velocities, even
when we vary $A_{\rm hydro}$ enough to induce $\gtrsim 25\%$ changes to the
SMF at $10^{11}M_{\odot}$.

Overall, our results indicate that HOD models can assign the kinematics of
satellite galaxies without considering a detailed dependence on
$\bm{\theta}_{\rm h}$, $\Omega_{\rm m}$, $\sigma_{8}$, and AGN and SN
feedback. However, for both central and satellite galaxies, we find
differences in $q_{\bm{\phi}}(\Delta r/R_{\rm vir}\,|\,\bm{\theta})$ between
IllustrisTNG and SIMBA, which indicates that hydrodynamic simulations have yet
to converge in their modeling of galaxy kinematics. This underscores the
necessity for HOD models to incorporate more flexible parameterizations of
galaxy kinematics that can describe the range of predictions in simulations.

###  5.2 Satellite Incompleteness

The radial distribution of our satellites significantly deviates from NFW in
both IllustrisTNG and SIMBA. The lack of satellites near the host halo centers
compared to the NFW profile is partly anticipated by known numerical
limitations in $N$-body simulations. The disruption of subhalos by dynamical
friction and tidal stripping can be unrealistically enhanced by gravitational
softening (Moore et al., 1996; Klypin et al., 1999; Peñarrubia et al., 2010;
van den Bosch, 2017; van den Bosch & Ogiya, 2018; van den Bosch et al., 2018).
Additional numerical limitations in the two-body relaxation, particle-halo
heating, and particle-halo tidal heating can also quickly disintegrate
subhalos (van Kampen, 1995, 2000; Moore et al., 1996; Klypin et al., 1999).
Even when subhalos are not fully disrupted, numerical effects may influence
the recovered subhalo properties (Mansfield & Avestruz, 2021) and halo-finding
algorithms may lose subhalos experiencing tidal deformation (Diemer et al.,
2023). These causes of satellite incompleteness are relevant even to subhalos
with a large number of particles (e.g., van den Bosch, 2017; van den Bosch &
Ogiya, 2018; Diemer et al., 2023).

To assess the degree of subhalo incompleteness, we compare the galaxy
kinematics derived from CAMELS with the UniverseMachine (UM; Behroozi et al.,
2019). UM is a semi-empirical model that employs analytic prescription to
account for both the numerical and physical disruption of satellite subhalos.
It is built on halo merger trees of Bolshoi-Planck $N$-body simulations
(Klypin et al., 2016; Rodríguez-Puebla et al., 2016) of comoving box size
equal to 250 $h^{-1}$Mpc. Once satellite subhalos become unidentifiable by
Rockstar, UM mitigates satellite incompleteness by analytically integrating
softened force-law to track “orphan” satellites until they are considered
disrupted or merged (see Appendix B of Behroozi et al., 2019).

In Figure 8, we compare the halo-centric radial and velocity distributions of
satellite galaxies in UM (gray shade), IllustrisTNG (red), and SIMBA (blue).
From UM, we take satellite galaxies at $z=0$ with $M_{\ast}>10^{9}M_{\odot}$
and host halo mass around our fiducial value, ${13.45}\leq\log M_{\rm
vir}\leq{13.55}$. For IllustrisTNG and SIMBA, we use our NDEs with
$\bm{\theta}_{\rm fid}$, which has consistent cosmology with the UM. Compared
to our $q_{\bm{\phi}}(\Delta r_{\rm sat}/R_{\rm vir}\,|\,\bm{\theta})$ and
$q_{\bm{\phi}}(\Delta v_{\rm sat}/V_{\rm vir}\,|\,\bm{\theta})$ from CAMELS,
the satellite $\Delta r$ and $\Delta v$ distributions are skewed toward the
center of the halo in UM. The contrast between the satellite distributions of
UM and the CAMELS simulations suggests that numerical disruption may indeed
affect satellites in hydrodynamic simulations. Moreover, this indicates that
we cannot take hydrodynamic simulations as the ground truth and that HOD
models should incorporate additional flexibility to account for our current
lack of understanding regarding the kinematics of satellites near the halo
center.

Figure 8: Comparison of the $\Delta r_{\rm sat}/R_{\rm vir}$ and $\Delta
v_{\rm sat}/V_{\rm vir}$ of satellite galaxies in IllustrisTNG (red), and
SIMBA (blue), and the UniverseMachine (UM, gray shade). The IllustrisTNG and
SIMBA distributions are derived from our NDEs with $\bm{\theta}_{\rm fid}$
while the UM distribution is derived from $M_{*}>10^{9}M_{\odot}$ UM satellite
galaxies at $z=0$ with ${13.45}\leq\log M_{\rm vir}\leq{13.55}$. UM uses
analytic prescriptions to mitigate satellite incompleteness. The comparison
suggests that IllustrisTNG and SIMBA may experience satellite incompleteness,
which significantly reduces the fraction of satellite galaxies close to the
host halo center.

###  5.3 Limitations and Outlook

The primary goal of this study is to investigate how the positions and
velocities of galaxies depend on the internal halo properties, cosmological
parameters, and feedback parameters. Currently, our NDEs focus solely on the
magnitudes of halo-centric position and velocity offsets, $\Delta r$ and
$\Delta v$, without considering correlations between them. We also do not
include information on anisotropies in satellite galaxy kinematics. Some
studies using observations from 2dFGRS and SDSS (e.g., Sales & Lambas, 2004;
Bailin et al., 2008) as well as simulations (e.g., Zentner et al., 2005;
Hadzhiyska et al., 2023) find evidence for anisotropies among satellite
galaxies. While these effects may affect how galaxies should be placed within
a halo in practice, they do not impact our results on the parameter
dependence. We also note that, in principle, our NDE-based approach can
incorporate correlations between $\Delta r$ and $\Delta v$ as well as
satellite anisotropies.

Additionally, while this work focuses solely on the internal properties of
host halos, many studies have demonstrated the important role of environmental
properties and assembly history of halos in the galaxy-halo connection (e.g.,
Sheth & Tormen, 2004; Shen et al., 2006; Croton et al., 2007; Hahn et al.,
2007, 2009; Borzyszkowski et al., 2017; Paranjape et al., 2018; Xu et al.,
2021; Delgado et al., 2022; Hadzhiyska et al., 2022). The shear field defined
by the tidal tensor of ambient mass distribution, for instance, may play a
pivotal role in shaping the satellite kinematics by affecting the orbits and
survival of satellites within halos (e.g., Hahn et al., 2007, 2009; Paranjape
et al., 2018). In future works, we plan to incorporate environment information
at both the dark matter particle level and halo catalog level using, e.g.,
$k^{\rm th}$ nearest neighbor (e.g., Banerjee & Abel, 2021) and graph neural
networks (e.g., Cranmer et al., 2019; Jespersen et al., 2022; Villanueva-
Domingo et al., 2022; Wu & Kragh Jespersen, 2023).

This work is the first in a series of re-evaluating the HOD using CAMELS to
develop a more accurate framework for modeling galaxy clustering at small
scales. In subsequent work, we will incorporate additional parameters, such as
the environmental properties of halos and angular information of galaxy
kinematics, into our NDE-based approach. This will provide a more complete
understanding of the galaxy-halo connection and whether the HOD model needs to
be expanded for additional parameter dependency.

Once the importance of additional parameters have been explored, we will
derive new analytic HOD models using symbolic regression (e.g., Cranmer,
2023). We will construct these models to also incorporate current modeling
uncertainties in the satellite kinematics near the halo center by using a
functional form that can describe satellites in hydrodynamic simulations and
with corrections for satellite incompleteness (Section 5.2). This will provide
a more flexible and comprehensive description of galaxy kinematics within
halos than current HOD models, effectively capturing any additional parameter
dependencies. Finally, we will demonstrate that our new HOD models can more
accurately describe higher-order galaxy clustering down to small scales and
test whether they can meet the precision level required to analyze upcoming
galaxy surveys.

In the upcoming work, we also plan to extend our work to additional
simulations. We extensively validated our NDEs and demonstrated their
robustness for our analysis (see Appendix A). However, incorporating more
simulations will further enhance their accuracy, especially by increasing the
number of satellites and massive host halos. We will also expand the
comparison of galaxy kinematics to other hydrodynamic simulations in CAMELS
that employ different subgrid prescriptions for modeling baryonic processes,
such as Astrid (Bird et al., 2022), EAGLE, Ramses (Teyssier, R., 2002), and
Magneticum (Hirschmann et al., 2014). Future CAMELS suites will also increase
the number of simulation parameters to 28, including $\Omega_{\rm b}$, $h$,
and additional stellar and galactic wind parameters. By using additional
simulations with more cosmological and astrophysical parameters, we will
further validate our findings and develop an even more robust and flexible HOD
model.

##  6 summary

A better understanding of galaxy kinematics will enable us to model galaxy
clustering more accurately and, thus, extract more constraining power from
galaxy surveys for cosmology, galaxy evolution, and galaxy-halo connection
studies. In this paper, we investigate the dependence of the central and
satellite galaxy kinematics on $\bm{\theta}$, the intrinsic host halo
properties (mass, spin, concentration), cosmology ($\Omega_{\rm m}$,
$\sigma_{8}$), and baryonic feedback from AGN and SNe ($A_{\rm SN1},A_{\rm
SN2},A_{\rm AGN1},A_{\rm AGN2}$).

We use simulated galaxies from 2,000 CAMELS hydrodynamic simulations
constructed using the IllustrisTNG and SIMBA galaxy formation models. We focus
on central and satellite galaxies with $M_{*}>10^{9}M_{\odot}$, identified
using the Rockstar halo finder. We apply NDE based on normalizing flows to
estimate the $p(\Delta r\,|\,\bm{\theta})$ and $p(\Delta v\,|\,\bm{\theta})$
of central and satellite galaxies. This enables us to accurately capture the
dependence of galaxy kinematics on $\bm{\theta}$ and quantify the impact of
each parameter on $p(\Delta r\,|\,\bm{\theta})$ and $p(\Delta
v\,|\,\bm{\theta})$.

We summarize our main findings below:

  1. (i)

Central galaxies have significant spatial bias and velocity bias with respect
to their host halos, consistent with previous observational evidence. We find
significant dependence in these biases on host halo mass, concentration, and
spin. In contrast, we find no significant dependence of central galaxy
kinematics on $\Omega_{\rm m}$, $\sigma_{8}$, and AGN and SN feedback.

  2. (ii)

Satellite galaxies in CAMELS have radial distributions that strongly deviate
from an NFW profile, skewing away from the halo center. We find a deviation
from NFW in both IllustrisTNG and SIMBA satellites. This discrepancy, which is
most significant near the halo center, suggests that hydrodynamic simulations
may experience satellite incompleteness. Additionally, the significant
differences in the radial distributions of satellites in IllustrisTNG and
SIMBA indicate that HOD models should incorporate additional flexibility when
modeling satellite positions.

  3. (iii)

The velocity distribution of satellite galaxies is bimodal due to orbiting and
infalling satellites (Appendix C). The infalling population accounts for $\sim
10\%$ of $M_{\ast}>10^{9}M_{\odot}$ satellites.

  4. (iv)

The velocity distribution of satellites has a slight dependence on host halo
spin which affects the bimodality. Besides this weak dependence, the
kinematics of satellite galaxies do not depend strongly on other host halo
properties, cosmology, or feedback parameters.

Our results provide timely insights into improving HOD models and our
understanding of the galaxy-halo connection. Galaxy clustering analyses are
currently pushing to small scales to exploit additional constraining power. At
the same time, upcoming galaxy surveys like DESI, PFS, and Euclid, with their
unprecedented statistical power, will have significantly higher accuracy
requirements for modeling galaxy clustering. Our efforts seek to develop a new
framework for HOD that will meet the upcoming requirements in modeling galaxy
clustering on small scales.

This paper is the first of a series that will re-examine the HOD using CAMELS.
In subsequent works, we will incorporate additional parameters, such as the
environmental properties of halos and angular information of galaxy kinematics
to reveal any additional parameter dependencies. We will then derive new
analytic HOD models that account for current uncertainties in the satellite
kinematics near the halo center due to satellite incompleteness in
simulations. Finally, we will demonstrate that our new HOD model can more
accurately describe the higher-order galaxy clustering down to small scales at
the precision level required by upcoming galaxy surveys. Ultimately, these
works will develop a new HOD framework that can more accurately, flexibly, and
robustly model galaxy clustering down to small scales.

##  Appendix A NDE Validation

Figure 9:  Comparison between the $\Delta r/R_{\rm vir}$ (top) and $\Delta
v/V_{\rm vir}$ (bottom) distributions drawn from our NDEs (solid line) and
from the test sample (shaded). We compare the distributions for central (a)
and satellite (b) galaxies in IllustrisTNG (red) and SIMBA (blue). The NDE
distributions are constructed using $\Delta r^{\prime}/R_{\rm vir}$ and
$\Delta v^{\prime}/V_{\rm vir}$ drawn from the NDEs at the $\bm{\theta}$
values of every test galaxy. We find excellent agreement between the NDE and
test sample distributions, which demonstrates the overall accuracy of our
NDEs.

In this section, we validate the accuracy of the NDEs, $q_{\bm{\phi}}(\Delta
r_{\rm cen}/R_{\rm vir}\,|\,\bm{\theta})$, $q_{\bm{\phi}}(\Delta r_{\rm
sat}/R_{\rm vir}\,|\,\bm{\theta})$, $q_{\bm{\phi}}(\Delta v_{\rm cen}/V_{\rm
vir}\,|\,\bm{\theta})$, and $q_{\bm{\phi}}(\Delta v_{\rm sat}/V_{\rm
vir}\,|\,\bm{\theta})$, using the randomly selected 10% of CAMELS galaxies
that we reserve for testing (Section 2). This test sample consists of 27,652
central galaxies and 3,469 satellites for IllustrisTNG and 35,949 central
galaxies and 2,224 satellites for SIMBA. For each test galaxy $i$, we compile
its parameters: $\bm{\theta}_{i}=(M_{{\rm vir},i}$, $c_{i}$, $\lambda_{i}$,
$\Omega_{{\rm m},i}$, $\sigma_{8,i}$, $A_{{\rm hydro},i})$. We then use the
NDEs to draw $\Delta r_{{\rm cen},i}^{\prime}$, $\Delta r_{{\rm
sat},i}^{\prime}$, $\Delta v_{{\rm cen},i}^{\prime}$, and $\Delta v_{{\rm
sat},i}^{\prime}$ for the given $\bm{\theta}_{i}$. Afterwards, we compare the
distributions of the kinematic quantities drawn from our NDEs to $p(\Delta
r_{{\rm cen},i}/R_{\rm vir})$, $p(\Delta r_{{\rm sat},i}/R_{\rm vir})$,
$p(\Delta v_{{\rm cen},i}/V_{\rm vir})$, and $p(\Delta v_{{\rm sat},i}/V_{\rm
vir})$ of the test galaxies.

In Figure 9, we present the comparison between the $p(\Delta r^{\prime}/R_{\rm
vir})$ and $p(\Delta r/R_{\rm vir})$ (top) and $p(\Delta v^{\prime}/V_{\rm
vir})$ and $p(\Delta v/V_{\rm vir})$ (bottom) for IllustrisTNG (red) and SIMBA
(blue), and for central (a) and satellite (b) galaxies. $p(\Delta
r^{\prime}/R_{\rm vir})$ and $p(\Delta v^{\prime}/V_{\rm vir})$ correspond to
the distributions drawn from our NDEs (solid line). We draw 15 samples from
our NDEs for each test galaxy to reduce the noise in the distribution.
Meanwhile, $p(\Delta r/R_{\rm vir})$ and $p(\Delta v/V_{\rm vir})$ correspond
to the distribution of all test galaxies (shaded). For both central and
satellite galaxies, and for both subgrid models, the distributions from our
NDEs are in excellent agreement with the test sample distributions. This
demonstrates the overall accuracy of our NDEs.

We further validate the NDEs, this time focusing on the parameter dependence.
Instead of comparing $p(\Delta r/R_{\rm vir})$ and $p(\Delta v/V_{\rm vir})$
over the full test sample, we split the sample into thirds along each
parameter of $\bm{\theta}$. Then, for each of the test subsamples, we repeat
the validation above. In Figure 10, we present the validation results for
central galaxies divided using $\log M_{\rm vir}$ for IllustrisTNG (red) and
SIMBA (blue). As in Figure 9, we mark the distributions of our NDE draws in
solid and the distributions from the test sample in the shaded regions. From
top to bottom, we show the comparison in low, intermediate, and high $M_{\rm
vir}$ ranges. In all of the $M_{\rm vir}$ ranges, for both IllustrisTNG and
SIMBA, the distributions are in excellent agreement. We similarly find good
agreements for all other parameters in $\bm{\theta}$ and for satellites.

Figure 10:  Comparison between the $\Delta r_{\rm cen}/R_{\rm vir}$ (a) and
$\Delta v_{\rm cen}/V_{\rm vir}$ (b) distributions drawn from our NDEs (solid
line) and from the test sample (shaded) divided into low, intermediate, and
high $M_{\rm vir}$ bins. We compare the distributions for IllustrisTNG (red)
and SIMBA (blue). The NDE distributions are constructed in the same fashion as
in Figure 9 using the test galaxies in each $M_{\rm vir}$ bin. By splitting
the test sample by $M_{\rm vir}$, our comparison examines the accuracy of our
NDEs, while accounting for parameter dependencies. We find good agreement
between the NDE and test sample distributions in every panel. We similarly
find good agreements for other components of $\bm{\theta}$ and for satellites.

Finally, we validate our NDEs using the Test of Accuracy with Random Points
(TARP; Lemos et al., 2023a) coverage test. TARP assesses the accuracy of NDEs
by estimating the expected coverage probabilities (ECPs) for given credibility
levels. For the test, we randomly select 10,000 central galaxies and use all
satellites from the IllustrisTNG and SIMBA test datasets. For each selected
test galaxy $i$ with $(\Delta a_{i},\bm{\theta}_{i})$, we draw 1,000
kinematics estimates, $\\{\Delta a_{i}^{\prime}\\}$, using our NDEs. We sample
a random point $\Delta a_{i,r}$ uniformly from the range of $\Delta a$ spanned
by the selected test galaxies. For a given credibility level, we estimate the
ECP from the fraction of $\\{\Delta a_{i}^{\prime}\\}$ with $|\Delta
a^{\prime}_{i}-\Delta a_{i,r}|<|\Delta a_{i}-\Delta a_{i,r}|$. For an optimal
NDE, ECP is equal to the credibility level.

In Figure 11, we compare the ECPs to the credibility levels for
$q_{\bm{\phi}}(\Delta a_{\rm cen}\,|\,\bm{\theta})$ (solid lines) and
$q_{\bm{\phi}}(\Delta a_{\rm sat}\,|\,\bm{\theta})$ (dashed lines) in
IllustrisTNG (a) and SIMBA (b), derived from TARP. The excellent agreement
between the ECPs and credibility levels demonstrates that all of our NDEs
provide accurate estimates of kinematics distributions. We therefore conclude
that our NDEs estimate $p(\Delta a|\bm{\theta})$ in CAMELS with sufficient
accuracy.

Figure 11: Comparison between credibility levels and corresponding expected
coverage probabilities for $q_{\bm{\phi}}(\Delta a_{\rm cen}\,|\,\bm{\theta})$
(solid lines) and $q_{\bm{\phi}}(\Delta a_{\rm sat}\,|\,\bm{\theta})$ (dashed
lines) in IllustrisTNG (a) and SIMBA (b). We compute the expected coverage
using randomly selected 10,000 central galaxies and all satellites in test
datasets of both hydrodynamic suites. The good agreement between the expected
coverage probability and the credibility level demonstrates that our NDEs
provide accurate estimates of kinematics distributions.

##  Appendix B $\Omega_{\rm m}$ dependence

In Section 4, we examine the dependence of central and satellite kinematics on
$\Omega_{\rm m}$. However, to estimate the dependence on $\Omega_{\rm m}$
alone, we must account for the fact that the definition of halo mass, $M_{\rm
vir}$, changes with $\Omega_{\rm m}$. In this section, we examine whether the
kinematics dependence, or lack thereof, on $\Omega_{\rm m}$ is affected if we
vary $M_{\rm vir}$ according to $\Omega_{\rm m}$.

To disentangle the dependence on $\Omega_{\rm m}$ from $M_{\rm vir}$, for a
given $\Omega_{\rm m}$ we set $M_{\rm vir}$ to the halo mass that matches the
fiducial number density of $M_{\rm vir}=10^{13.5}M_{\odot}$ halos at the
fiducial $\Omega_{\rm m}=0.3$ cosmology. We do this by measuring the halo mass
functions (HMF) for IllustrisTNG and SIMBA at $\Omega_{\rm m}=0.2$, 0.3, and
0.4, then matching them to the fiducial halo number density. For $\Omega_{\rm
m}=0.2$ and 0.4 we use $M_{\rm vir}\approx 10^{13.25}M_{\odot}$ and
$10^{13.72}M_{\odot}$ for IllustrisTNG, and $M_{\rm vir}\approx
10^{13.14}M_{\odot}$ and $10^{13.67}M_{\odot}$ for SIMBA.

In Figure 12, we present the central galaxy kinematics (a) and satellite
kinematics (b), $q_{\bm{\phi}}(\Delta r/R_{\rm vir}\,|\,\bm{\theta})$ (top)
and $q_{\bm{\phi}}(\Delta v/V_{\rm vir}\,|\,\bm{\theta})$ (bottom) for
different $\Omega_{\rm m}$ while varying $M_{\rm vir}$ to match the fiducial
number density. We include the distributions for IllustrisTNG (red) and SIMBA
(blue). Even after fixing number density, i.e., disentangling the potential
dependence on $M_{\rm vir}$, there is only a slight dependence on $\Omega_{\rm
m}$ for both central and satellite galaxy kinematics. This is consistent with
our findings in Section 4 and confirms that $\Omega_{\rm m}$ does not
significantly impact the kinematics of both central and satellite galaxies.

Figure 12: Dependence of the central galaxy kinematics (a) and satellite
kinematics (b), $q_{\bm{\phi}}(\Delta r/R_{\rm vir}\,|\,\bm{\theta})$ (top)
and $q_{\bm{\phi}}(\Delta v/V_{\rm vir}\,|\,\bm{\theta})$ (bottom), on
$\Omega_{\rm m}$ in IllustrisTNG (red) and SIMBA (blue), while varying $M_{\rm
vir}$ to a fixed fiducial halo number density. The parameters that do not
appear in the legends are kept at the fiducial values. The dashed lines in (a)
mark the medians of the distributions. Overall, varying $\Omega_{\rm m}$ does
not impact our findings in Section 4 and we recover the same weak $\Omega_{\rm
m}$ dependence for both central and satellite galaxy kinematics.  Figure 13:
(a) Satellite radial velocity, $v_{\rm rad}/V_{\rm vir}$ distribution as a
function of satellite radial position, $\Delta r_{\rm sat}/R_{\rm vir}$ for
all satellites in IllustrisTNG (top) and SIMBA (bottom). Negative $v_{\rm
rad}$ is toward the halo center. We classify infalling galaxies from orbiting
using an approximate criterion (dashed line) derived from Aung et al. (2023).
(b) $p(\Delta r_{\rm sat}/R_{\rm vir})$ of the satellite galaxies. We present
the total satellite population in black and the infalling population in dark
green. Infalling satellites contributes significantly near $\Delta r_{\rm
sat}=R_{\rm vir}$. (c) $p(\Delta v_{\rm sat}/V_{\rm vir})$ of the satellite
galaxies. The bimodality in $p(\Delta v_{\rm sat}/V_{\rm vir})$ is primarily
driven by infalling galaxies distributed around $\Delta v_{\rm sat}\approx
1.5V_{\rm vir}$.

##  Appendix C Infalling Satellites

We consistently find a bimodality in $q_{\bm{\phi}}(\Delta v_{\rm sat}/V_{\rm
vir}\,|\,\bm{\theta})$ for both IllustrisTNG and SIMBA (Section 4.2). One of
the modes is centered at $\Delta{v}_{\rm sat}/V_{\rm vir}\approx 1$ while the
other is centered at $\approx 1.5$. In this section, we examine whether these
modes come from separate orbiting and infalling populations.

We begin by selecting infalling galaxies using an approximate cut derived from
the position of infalling galaxies on $v_{\rm rad}$–$\Delta r_{\rm sat}$ plane
in Aung et al. (2023) (see Figure 2 therein): $v_{\rm rad}/V_{\rm vir}>\Delta
r_{\rm sat}/R_{\rm vir}-2$. We use this approximate selection as CAMELS
snapshots have low time resolution so we cannot accurately trace the
pericentric passages of our satellite galaxies. In Figure 13a, we show the
distribution of $\Delta r_{\rm sat}$ and $v_{\rm rad}$ for all IllustrisTNG
(red) and SIMBA (blue) satellite galaxies with $M_{\ast}>10^{9}M_{\odot}$. We
mark our approximate infall cut (black dashed) and the infalling population in
dark green.

In Figure 13b and c, we show the satellite radial distribution and velocity
offset distribution of the same IllustrisTNG (red) and SIMBA (blue) satellite
galaxies. We denote the total satellite population in black and the infalling
satellites in dark green. For the satellite radial distribution, we find that
infalling satellites contribute significantly near $\Delta r_{\rm sat}=R_{\rm
vir}$ and account for some of the differences between the IllustrisTNG and
SIMBA distributions. Furthermore, for the velocity offset distribution, we
find the second mode at $\Delta v_{\rm sat}/V_{\rm vir}\approx 1.5$ is
primarily due to the infalling population.

$\sim 10\%$ of the total satellite population consists of infalling
satellites. Despite their significant contribution, infalling satellites are
currently not included in standard HOD models. As we show in this section, the
kinematics of infalling galaxies are distinct from the orbiting satellites.
This suggests that excluding them in HOD models may significantly limit how
accurately we can model small-scale galaxy clustering, especially in redshift
space.

## acknowledgement

It is a pleasure to thank Andrew Hearin, Francisco Villaescusa-Navarro, Romain
Teyssier, and Tjitske Starkenburg for discussions on the comprehensive aspects
of hydrodynamic simulations as well as constructive feedback on the paper. The
authors also benefited from detailed discussions with Chirag Modi, Christopher
Lovell, John Wu, Peter Behroozi, Shy Genel, and Zixuan Peng. KJK acknowledges
support from the National Science Foundation under Grant No. 2012086. CH was
supported by the AI Accelerator program of the Schmidt Futures Foundation.

## References

  * Abazajian et al. (2005) Abazajian, K., Zheng, Z., Zehavi, I., et al. 2005, ApJ, 625, 613, doi: 10.1086/429685
  * Alam et al. (2017) Alam, S., Ata, M., Bailey, S., et al. 2017, MNRAS, 470, 2617, doi: 10.1093/mnras/stx721
  * Alsing et al. (2019) Alsing, J., Charnock, T., Feeney, S., & Wandelt, B. 2019, MNRAS, 488, 4440, doi: 10.1093/mnras/stz1960
  * Aung et al. (2023) Aung, H., Nagai, D., Rozo, E., Wolfe, B., & Adhikari, S. 2023, MNRAS, 521, 3981, doi: 10.1093/mnras/stad601
  * Bailin et al. (2008) Bailin, J., Power, C., Norberg, P., Zaritsky, D., & Gibson, B. K. 2008, MNRAS, 390, 1133, doi: 10.1111/j.1365-2966.2008.13828.x
  * Banerjee & Abel (2021) Banerjee, A., & Abel, T. 2021, MNRAS, 500, 5479, doi: 10.1093/mnras/staa3604
  * 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. (2000) Benson, A. J., Cole, S., Frenk, C. S., Baugh, C. M., & Lacey, C. G. 2000, MNRAS, 311, 793, doi: 10.1046/j.1365-8711.2000.03101.x
  * Berlind & Weinberg (2002) Berlind, A. A., & Weinberg, D. H. 2002, ApJ, 575, 587, doi: 10.1086/341469
  * Berlind et al. (2003) Berlind, A. A., Weinberg, D. H., Benson, A. J., et al. 2003, ApJ, 593, 1, doi: 10.1086/376517
  * Bird et al. (2022) Bird, S., Ni, Y., Di Matteo, T., et al. 2022, MNRAS, 512, 3703, doi: 10.1093/mnras/stac648
  * Borzyszkowski et al. (2017) Borzyszkowski, M., Porciani, C., Romano-Díaz, E., & Garaldi, E. 2017, MNRAS, 469, 594, doi: 10.1093/mnras/stx873
  * Bryan & Norman (1998) Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80, doi: 10.1086/305262
  * Budzynski et al. (2012) Budzynski, J. M., Koposov, S. E., McCarthy, I. G., McGee, S. L., & Belokurov, V. 2012, MNRAS, 423, 104, doi: 10.1111/j.1365-2966.2012.20663.x
  * Castorina & Sheth (2013) Castorina, E., & Sheth, R. K. 2013, MNRAS, 433, 1529, doi: 10.1093/mnras/stt824
  * Chan et al. (2015) Chan, T. K., Kereš, D., Oñorbe, J., et al. 2015, MNRAS, 454, 2981, doi: 10.1093/mnras/stv2165
  * Chisari et al. (2019) Chisari, N. E., Mead, A. J., Joudaki, S., et al. 2019, The Open Journal of Astrophysics, 2, 4, doi: 10.21105/astro.1905.06082
  * Cole et al. (2005) Cole, S., Percival, W. J., Peacock, J. A., et al. 2005, MNRAS, 362, 505, doi: 10.1111/j.1365-2966.2005.09318.x
  * Cranmer (2023) Cranmer, M. 2023, arXiv e-prints, arXiv:2305.01582, doi: 10.48550/arXiv.2305.01582
  * Cranmer et al. (2019) Cranmer, M. D., Xu, R., Battaglia, P., & Ho, S. 2019, arXiv e-prints, arXiv:1909.05862, doi: 10.48550/arXiv.1909.05862
  * Croton et al. (2007) Croton, D. J., Gao, L., & White, S. D. M. 2007, MNRAS, 374, 1303, doi: 10.1111/j.1365-2966.2006.11230.x
  * Dai & Seljak (2022) Dai, B., & Seljak, U. 2022, MNRAS, 516, 2363, doi: 10.1093/mnras/stac2010
  * Davé et al. (2019) Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, Mon. Not. Roy. Astron. Soc., 486, 2827, doi: 10.1093/mnras/stz937
  * Dawson et al. (2013a) Dawson, K. S., Schlegel, D. J., Ahn, C. P., et al. 2013a, AJ, 145, 10, doi: 10.1088/0004-6256/145/1/10
  * Dawson et al. (2013b) —. 2013b, AJ, 145, 10, doi: 10.1088/0004-6256/145/1/10
  * Dawson et al. (2016) Dawson, K. S., Kneib, J.-P., Percival, W. J., et al. 2016, AJ, 151, 44, doi: 10.3847/0004-6256/151/2/44
  * Delgado et al. (2022) Delgado, A. M., Wadekar, D., Hadzhiyska, B., et al. 2022, MNRAS, 515, 2733, doi: 10.1093/mnras/stac1951
  * DESI Collaboration et al. (2022) DESI Collaboration, Abareshi, B., Aguilar, J., et al. 2022, AJ, 164, 207, doi: 10.3847/1538-3881/ac882b
  * Diemer et al. (2023) Diemer, B., Behroozi, P., & Mansfield, P. 2023, arXiv e-prints, arXiv:2305.00993, doi: 10.48550/arXiv.2305.00993
  * Durkan et al. (2019) Durkan, C., Bekasov, A., Murray, I., & Papamakarios, G. 2019, arXiv e-prints, arXiv:1906.04032, doi: 10.48550/arXiv.1906.04032
  * Durkan et al. (2020) Durkan, C., Bekasov, A., Murray, I., & Papamakarios, G. 2020, nflows: normalizing flows in PyTorch, v0.14, Zenodo, doi: 10.5281/zenodo.4296287
  * Eickenberg et al. (2022) Eickenberg, M., Allys, E., Moradinezhad Dizgah, A., et al. 2022, arXiv e-prints, arXiv:2204.07646, doi: 10.48550/arXiv.2204.07646
  * Eisenstein et al. (2001) Eisenstein, D. J., Annis, J., Gunn, J. E., et al. 2001, AJ, 122, 2267, doi: 10.1086/323717
  * Foreman et al. (2020) Foreman, S., Coulton, W., Villaescusa-Navarro, F., & Barreira, A. 2020, MNRAS, 498, 2887, doi: 10.1093/mnras/staa2523
  * Germain et al. (2015) Germain, M., Gregor, K., Murray, I., & Larochelle, H. 2015, arXiv e-prints, arXiv:1502.03509, doi: 10.48550/arXiv.1502.03509
  * Gill et al. (2005) Gill, S. P. D., Knebe, A., & Gibson, B. K. 2005, MNRAS, 356, 1327, doi: 10.1111/j.1365-2966.2004.08562.x
  * Greenberg et al. (2019) Greenberg, D. S., Nonnenmacher, M., & Macke, J. H. 2019, arXiv e-prints, arXiv:1905.07488, doi: 10.48550/arXiv.1905.07488
  * Gualdi et al. (2021) Gualdi, D., Gil-Marín, H., & Verde, L. 2021, J. Cosmology Astropart. Phys, 2021, 008, doi: 10.1088/1475-7516/2021/07/008
  * Guo et al. (2015a) Guo, H., Zheng, Z., Zehavi, I., et al. 2015a, MNRAS, 453, 4368, doi: 10.1093/mnras/stv1966
  * Guo et al. (2015b) —. 2015b, MNRAS, 446, 578, doi: 10.1093/mnras/stu2120
  * Hadzhiyska et al. (2021) Hadzhiyska, B., Liu, S., Somerville, R. S., et al. 2021, MNRAS, 508, 698, doi: 10.1093/mnras/stab2564
  * Hadzhiyska et al. (2022) Hadzhiyska, B., Eisenstein, D., Hernquist, L., et al. 2022, arXiv e-prints, arXiv:2210.10072, doi: 10.48550/arXiv.2210.10072
  * Hadzhiyska et al. (2023) Hadzhiyska, B., Hernquist, L., Eisenstein, D., et al. 2023, MNRAS, 524, 2524, doi: 10.1093/mnras/stad279
  * Hahn & Melchior (2022) Hahn, C., & Melchior, P. 2022, ApJ, 938, 11, doi: 10.3847/1538-4357/ac7b84
  * Hahn & Villaescusa-Navarro (2021a) Hahn, C., & Villaescusa-Navarro, F. 2021a, J. Cosmology Astropart. Phys, 2021, 029, doi: 10.1088/1475-7516/2021/04/029
  * Hahn & Villaescusa-Navarro (2021b) —. 2021b, J. Cosmology Astropart. Phys, 2021, 029, doi: 10.1088/1475-7516/2021/04/029
  * Hahn et al. (2022) Hahn, C., Eickenberg, M., Ho, S., et al. 2022, arXiv e-prints, arXiv:2211.00723, doi: 10.48550/arXiv.2211.00723
  * Hahn et al. (2023a) Hahn, C., Lemos, P., Parker, L., et al. 2023a, arXiv e-prints, arXiv:2310.15246.  https://arxiv.org/abs/2310.15246
  * Hahn et al. (2023b) Hahn, C., Eickenberg, M., Ho, S., et al. 2023b, arXiv e-prints, arXiv:2310.15243.  https://arxiv.org/abs/2310.15243
  * Hahn et al. (2023c) Hahn, C., Wilson, M. J., Ruiz-Macias, O., et al. 2023c, AJ, 165, 253, doi: 10.3847/1538-3881/accff8
  * Hahn et al. (2007) Hahn, O., Carollo, C. M., Porciani, C., & Dekel, A. 2007, MNRAS, 381, 41, doi: 10.1111/j.1365-2966.2007.12249.x
  * Hahn et al. (2009) Hahn, O., Porciani, C., Dekel, A., & Carollo, C. M. 2009, MNRAS, 398, 1742, doi: 10.1111/j.1365-2966.2009.15271.x
  * Hassan et al. (2022) Hassan, S., Villaescusa-Navarro, F., Wandelt, B., et al. 2022, ApJ, 937, 83, doi: 10.3847/1538-4357/ac8b09
  * Hearin & Watson (2013) Hearin, A. P., & Watson, D. F. 2013, MNRAS, 435, 1313, doi: 10.1093/mnras/stt1374
  * Hearin et al. (2016) Hearin, A. P., Zentner, A. R., van den Bosch, F. C., Campbell, D., & Tollerud, E. 2016, MNRAS, 460, 2552, doi: 10.1093/mnras/stw840
  * Hellwing et al. (2016) Hellwing, W. A., Schaller, M., Frenk, C. S., et al. 2016, MNRAS, 461, L11, doi: 10.1093/mnrasl/slw081
  * Hirschmann et al. (2014) Hirschmann, M., Dolag, K., Saro, A., et al. 2014, MNRAS, 442, 2304, doi: 10.1093/mnras/stu1023
  * Hou et al. (2023) Hou, J., Moradinezhad Dizgah, A., Hahn, C., & Massara, E. 2023, J. Cosmology Astropart. Phys, 2023, 045, doi: 10.1088/1475-7516/2023/03/045
  * Huterer (2023) Huterer, D. 2023, A&A Rev., 31, 2, doi: 10.1007/s00159-023-00147-4
  * Ivanov & Philcox (2023) Ivanov, M. M., & Philcox, O. H. E. 2023, arXiv e-prints, arXiv:2305.07977, doi: 10.48550/arXiv.2305.07977
  * Ivanov et al. (2020) Ivanov, M. M., Simonović, M., & Zaldarriaga, M. 2020, J. Cosmology Astropart. Phys, 2020, 042, doi: 10.1088/1475-7516/2020/05/042
  * Jespersen et al. (2022) Jespersen, C. K., Cranmer, M., Melchior, P., et al. 2022, ApJ, 941, 7, doi: 10.3847/1538-4357/ac9b18
  * Kingma & Ba (2017) Kingma, D. P., & Ba, J. 2017, Adam: A Method for Stochastic Optimization.  https://arxiv.org/abs/1412.6980
  * Klypin et al. (1999) Klypin, A., Gottlöber, S., Kravtsov, A. V., & Khokhlov, A. M. 1999, ApJ, 516, 530, doi: 10.1086/307122
  * Klypin et al. (2016) Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340, doi: 10.1093/mnras/stw248
  * Kobayashi et al. (2022) Kobayashi, Y., Nishimichi, T., Takada, M., & Miyatake, H. 2022, Phys. Rev. D, 105, 083517, doi: 10.1103/PhysRevD.105.083517
  * Lange et al. (2023) Lange, J. U., Hearin, A. P., Leauthaud, A., et al. 2023, MNRAS, 520, 5373, doi: 10.1093/mnras/stad473
  * Lange et al. (2019) Lange, J. U., van den Bosch, F. C., Zentner, A. R., Wang, K., & Villarreal, A. S. 2019, MNRAS, 482, 4824, doi: 10.1093/mnras/sty2950
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints, arXiv:1110.3193
  * Lemos et al. (2023a) Lemos, P., Coogan, A., Hezaveh, Y., & Perreault-Levasseur, L. 2023a, 40th International Conference on Machine Learning, 202, 19256, doi: 10.48550/arXiv.2302.03026
  * Lemos et al. (2023b) Lemos, P., Parker, L., Hahn, C., et al. 2023b, arXiv e-prints, arXiv:2310.15256.  https://arxiv.org/abs/2310.15256
  * Lovell et al. (2023) Lovell, C. C., Hassan, S., Villaescusa-Navarro, F., et al. 2023, in Machine Learning for Astrophysics, 21, doi: 10.48550/arXiv.2307.06967
  * Mansfield & Avestruz (2021) Mansfield, P., & Avestruz, C. 2021, MNRAS, 500, 3309, doi: 10.1093/mnras/staa3388
  * Mao et al. (2015) Mao, Y.-Y., Williamson, M., & Wechsler, R. H. 2015, ApJ, 810, 21, doi: 10.1088/0004-637X/810/1/21
  * Masjedi et al. (2006) Masjedi, M., Hogg, D. W., Cool, R. J., et al. 2006, ApJ, 644, 54, doi: 10.1086/503536
  * Massara et al. (2021) Massara, E., Villaescusa-Navarro, F., Ho, S., Dalal, N., & Spergel, D. N. 2021, Phys. Rev. Lett., 126, 011301, doi: 10.1103/PhysRevLett.126.011301
  * Massara et al. (2023) Massara, E., Villaescusa-Navarro, F., Hahn, C., et al. 2023, ApJ, 951, 70, doi: 10.3847/1538-4357/acd44d
  * McDonough & Brainerd (2022) McDonough, B., & Brainerd, T. G. 2022, ApJ, 933, 161, doi: 10.3847/1538-4357/ac752d
  * Mihos & Hernquist (1994) Mihos, J. C., & Hernquist, L. 1994, ApJ, 425, L13, doi: 10.1086/187299
  * Moore et al. (1996) Moore, B., Katz, N., & Lake, G. 1996, ApJ, 457, 455, doi: 10.1086/176745
  * More et al. (2011) More, S., van den Bosch, F. C., Cacciato, M., et al. 2011, MNRAS, 410, 210, doi: 10.1111/j.1365-2966.2010.17436.x
  * Nagai & Kravtsov (2005) Nagai, D., & Kravtsov, A. V. 2005, ApJ, 618, 557, doi: 10.1086/426016
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563, doi: 10.1086/177173
  * Nelson et al. (2019) Nelson, D., Springel, V., Pillepich, A., et al. 2019, Computational Astrophysics and Cosmology, 6, 2, doi: 10.1186/s40668-019-0028-x
  * Paillas et al. (2023) Paillas, E., Cuesta-Lazaro, C., Percival, W. J., et al. 2023, arXiv e-prints, arXiv:2309.16541, doi: 10.48550/arXiv.2309.16541
  * Pakmor et al. (2023) Pakmor, R., Springel, V., Coles, J. P., et al. 2023, MNRAS, 524, 2539, doi: 10.1093/mnras/stac3620
  * Papamakarios et al. (2017) Papamakarios, G., Pavlakou, T., & Murray, I. 2017, arXiv e-prints, arXiv:1705.07057, doi: 10.48550/arXiv.1705.07057
  * Paranjape et al. (2018) Paranjape, A., Hahn, O., & Sheth, R. K. 2018, MNRAS, 476, 3631, doi: 10.1093/mnras/sty496
  * Peñarrubia et al. (2010) Peñarrubia, J., Benson, A. J., Walker, M. G., et al. 2010, MNRAS, 406, 1290, doi: 10.1111/j.1365-2966.2010.16762.x
  * Peacock & Smith (2000a) Peacock, J. A., & Smith, R. E. 2000a, MNRAS, 318, 1144, doi: 10.1046/j.1365-8711.2000.03779.x
  * Peacock & Smith (2000b) —. 2000b, MNRAS, 318, 1144, doi: 10.1046/j.1365-8711.2000.03779.x
  * Peacock et al. (2001) Peacock, J. A., Cole, S., Norberg, P., et al. 2001, Nature, 410, 169, doi: 10.1038/35065528
  * Peebles (1969) Peebles, P. J. E. 1969, ApJ, 155, 393, doi: 10.1086/149876
  * Percival et al. (2001) Percival, W. J., Baugh, C. M., Bland-Hawthorn, J., et al. 2001, MNRAS, 327, 1297, doi: 10.1046/j.1365-8711.2001.04827.x
  * Pillepich et al. (2018) Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077, doi: 10.1093/mnras/stx2656
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Porth et al. (2023) Porth, L., Bernstein, G. M., Smith, R. E., & Lee, A. J. 2023, MNRAS, 518, 3344, doi: 10.1093/mnras/stac3225
  * Power et al. (2003) Power, C., Navarro, J. F., Jenkins, A., et al. 2003, MNRAS, 338, 14, doi: 10.1046/j.1365-8711.2003.05925.x
  * Reid et al. (2014) Reid, B. A., Seo, H.-J., Leauthaud, A., Tinker, J. L., & White, M. 2014, MNRAS, 444, 476, doi: 10.1093/mnras/stu1391
  * Reid et al. (2012) Reid, B. A., Samushia, L., White, M., et al. 2012, MNRAS, 426, 2719, doi: 10.1111/j.1365-2966.2012.21779.x
  * Rocha et al. (2012) Rocha, M., Peter, A. H. G., & Bullock, J. 2012, MNRAS, 425, 231, doi: 10.1111/j.1365-2966.2012.21432.x
  * Rodríguez-Puebla et al. (2016) Rodríguez-Puebla, A., Behroozi, P., Primack, J., et al. 2016, MNRAS, 462, 893, doi: 10.1093/mnras/stw1705
  * Ruhe et al. (2022) Ruhe, D., Wong, K., Cranmer, M., & Forré, P. 2022, arXiv e-prints, arXiv:2211.09008, doi: 10.48550/arXiv.2211.09008
  * Sales & Lambas (2004) Sales, L., & Lambas, D. G. 2004, MNRAS, 348, 1236, doi: 10.1111/j.1365-2966.2004.07443.x
  * Schaye et al. (2015) Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521, doi: 10.1093/mnras/stu2058
  * Schneider et al. (2019) Schneider, A., Teyssier, R., Stadel, J., et al. 2019, J. Cosmology Astropart. Phys, 2019, 020, doi: 10.1088/1475-7516/2019/03/020
  * Scoccimarro et al. (2001) Scoccimarro, R., Sheth, R. K., Hui, L., & Jain, B. 2001, ApJ, 546, 20, doi: 10.1086/318261
  * Sefusatti & Scoccimarro (2005) Sefusatti, E., & Scoccimarro, R. 2005, Phys. Rev. D, 71, 063001, doi: 10.1103/PhysRevD.71.063001
  * Seljak (2000) Seljak, U. 2000, MNRAS, 318, 203, doi: 10.1046/j.1365-8711.2000.03715.x
  * Shen et al. (2006) Shen, J., Abel, T., Mo, H. J., & Sheth, R. K. 2006, ApJ, 645, 783, doi: 10.1086/504513
  * Sheth & Tormen (2004) Sheth, R. K., & Tormen, G. 2004, MNRAS, 350, 1385, doi: 10.1111/j.1365-2966.2004.07733.x
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, Wide-Field InfrarRed Survey Telescope-Astrophysics Focused Telescope Assets WFIRST-AFTA 2015 Report
  * Springel et al. (2018) Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676, doi: 10.1093/mnras/stx3304
  * Steinmetz & Bartelmann (1995) Steinmetz, M., & Bartelmann, M. 1995, MNRAS, 272, 570, doi: 10.1093/mnras/272.3.570
  * Storey-Fisher et al. (2022) Storey-Fisher, K., Tinker, J., Zhai, Z., et al. 2022, arXiv e-prints, arXiv:2210.03203, doi: 10.48550/arXiv.2210.03203
  * Tabak & Turner (2013) Tabak, E. G., & Turner, C. V. 2013, Communications on Pure and Applied Mathematics, 66, 145, doi: https://doi.org/10.1002/cpa.21423
  * Tabak & Vanden-Eijnden (2010) Tabak, E. G., & Vanden-Eijnden, E. 2010, Communications in Mathematical Sciences, 8, 217
  * Takada et al. (2014) Takada, M., Ellis, R. S., Chiba, M., et al. 2014, Publications of the Astronomical Society of Japan, 66, R1, doi: 10.1093/pasj/pst019
  * Tal et al. (2012) Tal, T., Wake, D. A., & van Dokkum, P. G. 2012, ApJ, 751, L5, doi: 10.1088/2041-8205/751/1/L5
  * Tegmark et al. (2004) Tegmark, M., Strauss, M. A., Blanton, M. R., et al. 2004, Phys. Rev. D, 69, 103501, doi: 10.1103/PhysRevD.69.103501
  * Tejero-Cantero et al. (2020) Tejero-Cantero, A., Boelts, J., Deistler, M., et al. 2020, Journal of Open Source Software, 5, 2505, doi: 10.21105/joss.02505
  * Teyssier, R. (2002) Teyssier, R. 2002, A&A, 385, 337, doi: 10.1051/0004-6361:20011817
  * Thiele et al. (2023) Thiele, L., Massara, E., Pisani, A., et al. 2023, arXiv e-prints, arXiv:2307.07555, doi: 10.48550/arXiv.2307.07555
  * Tinker et al. (2012) Tinker, J. L., Sheldon, E. S., Wechsler, R. H., et al. 2012, ApJ, 745, 16, doi: 10.1088/0004-637X/745/1/16
  * Valogiannis & Dvorkin (2022) Valogiannis, G., & Dvorkin, C. 2022, Phys. Rev. D, 105, 103534, doi: 10.1103/PhysRevD.105.103534
  * Valogiannis et al. (2023) Valogiannis, G., Yuan, S., & Dvorkin, C. 2023, arXiv e-prints, arXiv:2310.16116, doi: 10.48550/arXiv.2310.16116
  * van Daalen et al. (2020) van Daalen, M. P., McCarthy, I. G., & Schaye, J. 2020, MNRAS, 491, 2424, doi: 10.1093/mnras/stz3199
  * van Daalen et al. (2011) van Daalen, M. P., Schaye, J., Booth, C. M., & Dalla Vecchia, C. 2011, MNRAS, 415, 3649, doi: 10.1111/j.1365-2966.2011.18981.x
  * van den Bosch (2017) van den Bosch, F. C. 2017, MNRAS, 468, 885, doi: 10.1093/mnras/stx520
  * van den Bosch et al. (2013) van den Bosch, F. C., More, S., Cacciato, M., Mo, H., & Yang, X. 2013, MNRAS, 430, 725, doi: 10.1093/mnras/sts006
  * van den Bosch & Ogiya (2018) van den Bosch, F. C., & Ogiya, G. 2018, MNRAS, 475, 4066, doi: 10.1093/mnras/sty084
  * van den Bosch et al. (2018) van den Bosch, F. C., Ogiya, G., Hahn, O., & Burkert, A. 2018, MNRAS, 474, 3043, doi: 10.1093/mnras/stx2956
  * van den Bosch et al. (2005) van den Bosch, F. C., Weinmann, S. M., Yang, X., et al. 2005, MNRAS, 361, 1203, doi: 10.1111/j.1365-2966.2005.09260.x
  * van Kampen (1995) van Kampen, E. 1995, MNRAS, 273, 295, doi: 10.1093/mnras/273.2.295
  * van Kampen (2000) —. 2000, arXiv e-prints, astro, doi: 10.48550/arXiv.astro-ph/0002027
  * van Kampen et al. (1999) van Kampen, E., Jimenez, R., & Peacock, J. A. 1999, MNRAS, 310, 43, doi: 10.1046/j.1365-8711.1999.02955.x
  * Verde et al. (2002) Verde, L., Heavens, A. F., Percival, W. J., et al. 2002, MNRAS, 335, 432, doi: 10.1046/j.1365-8711.2002.05620.x
  * Villaescusa-Navarro et al. (2021) Villaescusa-Navarro, F., Anglés-Alcázar, D., Genel, S., et al. 2021, ApJ, 915, 71, doi: 10.3847/1538-4357/abf7ba
  * Villaescusa-Navarro et al. (2022) Villaescusa-Navarro, F., Genel, S., Anglés-Alcázar, D., et al. 2022, arXiv e-prints, arXiv:2201.01300.  https://arxiv.org/abs/2201.01300
  * Villanueva-Domingo et al. (2022) Villanueva-Domingo, P., Villaescusa-Navarro, F., Anglés-Alcázar, D., et al. 2022, ApJ, 935, 30, doi: 10.3847/1538-4357/ac7aa3
  * Wang et al. (2022) Wang, Y., Zhai, Z., Alavi, A., et al. 2022, The Astrophysical Journal, 928, 1, doi: 10.3847/1538-4357/ac4973
  * Watson et al. (2012a) Watson, D. F., Berlind, A. A., McBride, C. K., Hogg, D. W., & Jiang, T. 2012a, ApJ, 749, 83, doi: 10.1088/0004-637X/749/1/83
  * Watson et al. (2012b) Watson, D. F., Berlind, A. A., & Zentner, A. R. 2012b, ApJ, 754, 90, doi: 10.1088/0004-637X/754/2/90
  * Wechsler & Tinker (2018) Wechsler, R. H., & Tinker, J. L. 2018, ARA&A, 56, 435, doi: 10.1146/annurev-astro-081817-051756
  * Weinberger et al. (2017) Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291, doi: 10.1093/mnras/stw2944
  * White et al. (2007) White, M., Zheng, Z., Brown, M. J. I., Dey, A., & Jannuzi, B. T. 2007, ApJ, 655, L69, doi: 10.1086/512015
  * White et al. (2011) White, M., Blanton, M., Bolton, A., et al. 2011, ApJ, 728, 126, doi: 10.1088/0004-637X/728/2/126
  * Wright et al. (2022) Wright, R. J., Lagos, C. d. P., Power, C., et al. 2022, MNRAS, 516, 2891, doi: 10.1093/mnras/stac2042
  * Wu & Kragh Jespersen (2023) Wu, J. F., & Kragh Jespersen, C. 2023, arXiv e-prints, arXiv:2306.12327, doi: 10.48550/arXiv.2306.12327
  * Xu et al. (2021) Xu, X., Zehavi, I., & Contreras, S. 2021, MNRAS, 502, 3242, doi: 10.1093/mnras/stab100
  * Yang et al. (2003) Yang, X., Mo, H. J., & van den Bosch, F. C. 2003, MNRAS, 339, 1057, doi: 10.1046/j.1365-8711.2003.06254.x
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579, doi: 10.1086/301513
  * Yuan et al. (2022a) Yuan, S., Garrison, L. H., Eisenstein, D. J., & Wechsler, R. H. 2022a, MNRAS, 515, 871, doi: 10.1093/mnras/stac1830
  * Yuan et al. (2022b) —. 2022b, MNRAS, 515, 871, doi: 10.1093/mnras/stac1830
  * Yuan et al. (2022c) Yuan, S., Hadzhiyska, B., Bose, S., & Eisenstein, D. J. 2022c, MNRAS, 512, 5793, doi: 10.1093/mnras/stac830
  * Zehavi et al. (2018) Zehavi, I., Contreras, S., Padilla, N., et al. 2018, ApJ, 853, 84, doi: 10.3847/1538-4357/aaa54a
  * Zehavi et al. (2011) Zehavi, I., Zheng, Z., Weinberg, D. H., et al. 2011, ApJ, 736, 59, doi: 10.1088/0004-637X/736/1/59
  * Zentner et al. (2014) Zentner, A. R., Hearin, A. P., & van den Bosch, F. C. 2014, MNRAS, 443, 3044, doi: 10.1093/mnras/stu1383
  * Zentner et al. (2005) Zentner, A. R., Kravtsov, A. V., Gnedin, O. Y., & Klypin, A. A. 2005, ApJ, 629, 219, doi: 10.1086/431355
  * Zhai et al. (2019) Zhai, Z., Tinker, J. L., Becker, M. R., et al. 2019, ApJ, 874, 95, doi: 10.3847/1538-4357/ab0d7b
  * Zhai et al. (2023) Zhai, Z., Tinker, J. L., Banerjee, A., et al. 2023, ApJ, 948, 99, doi: 10.3847/1538-4357/acc65b
  * Zhang et al. (2019) Zhang, T., Liao, S., Li, M., & Gao, L. 2019, MNRAS, 487, 1227, doi: 10.1093/mnras/stz1370
  * Zheng et al. (2007) Zheng, Z., Coil, A. L., & Zehavi, I. 2007, ApJ, 667, 760, doi: 10.1086/521074
  * Zheng & Weinberg (2007) Zheng, Z., & Weinberg, D. H. 2007, ApJ, 659, 1, doi: 10.1086/512151
  * Zheng et al. (2009) Zheng, Z., Zehavi, I., Eisenstein, D. J., Weinberg, D. H., & Jing, Y. P. 2009, ApJ, 707, 554, doi: 10.1088/0004-637X/707/1/554
  * Zheng et al. (2005) Zheng, Z., Berlind, A. A., Weinberg, D. H., et al. 2005, ApJ, 633, 791, doi: 10.1086/466510
  * Zhou et al. (2023) Zhou, R., Dey, B., Newman, J. A., et al. 2023, AJ, 165, 58, doi: 10.3847/1538-3881/aca5fb
  * Zitrin et al. (2012) Zitrin, A., Bartelmann, M., Umetsu, K., Oguri, M., & Broadhurst, T. 2012, MNRAS, 426, 2944, doi: 10.1111/j.1365-2966.2012.21886.x
  * Zu & Weinberg (2013) Zu, Y., & Weinberg, D. H. 2013, MNRAS, 431, 3319, doi: 10.1093/mnras/stt411
  * Zu et al. (2014) Zu, Y., Weinberg, D. H., Jennings, E., Li, B., & Wyman, M. 2014, MNRAS, 445, 1885, doi: 10.1093/mnras/stu1739