# The galaxy–halo connection of DESI luminous red galaxies with subhalo
abundance matching

Angela M. Berti11affiliation: Department of Physics & Astronomy, University of
Utah, Salt Lake City, UT 84112, USA , Kyle S. Dawson11affiliation: Department
of Physics & Astronomy, University of Utah, Salt Lake City, UT 84112, USA ,
Wilber Dominguez22affiliation: Department of Physics & Astronomy, Swarthmore
College, 500 College Ave, Swarthmore, PA 19081, USA

###### Abstract

We use subhalo abundance and age distribution matching to create magnitude-
limited mock galaxy catalogs at $z\sim 0.43$, 0.52, and 0.63 with $z$-band and
3.4 micron $W1$-band absolute magnitudes and ${r-z}$ and ${r-W1}$ colors. From
these magnitude-limited mocks we select mock luminous red galaxy (LRG) samples
according to the $(r-z)$-based (optical) and $(r-W1)$-based (infrared)
selection criteria for the LRG sample of the Dark Energy Spectroscopic
Instrument (DESI) Survey. Our models reproduce the number densities,
luminosity functions, color distributions, and projected clustering of the
DESI Legacy Surveys that are the basis for DESI LRG target selection. We
predict the halo occupation statistics of both optical and IR DESI LRGs at
fixed cosmology, and assess the differences between the two LRG samples. We
find that IR-based SHAM modeling represents the differences between the
optical and IR LRG populations better than using the $z$-band, and that age
distribution matching overpredicts the clustering of LRGs, implying that
galaxy color is uncorrelated with halo age in the LRG regime. Both the optical
and IR DESI LRG target selections exclude some of the most luminous galaxies
that would appear to be LRGs based on their position on the red sequence in
optical color–magnitude space. Both selections also yield populations with a
non-trivial LRG–halo connection that does not reach unity for the most massive
halos. We find the IR selection achieves greater completeness ($\gtrsim 90\%$)
than the optical selection across all redshift bins studied.

##  1\. Introduction

The Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration et al.
2016a, b) survey is a spectroscopic galaxy redshift survey of unprecedented
scale that will classify tens of millions of galaxies in four target classes
over $\sim 14,000$ square degrees at the lowest redshifts of this program.
DESI’s bright galaxy survey (BGS; Hahn et al. 2022) extends to $z\sim 0.4$,
while the the luminous red galaxy (LRG; Zhou et al. 2023) sample will reach to
$z\sim 1$ and cover 20 times the volume of the BGS111DESI will observe
emission line galaxies (ELGs) to $z\sim 1.6$, and quasars to $z\sim 3.5$..

The DESI target samples are optimized for precision measurements of
cosmological parameters. However, DESI also offers novel opportunities to
study galaxy evolution and the high-mass end of the stellar-to-halo mass
relation (SHMR), provided that sample selection effects are well understood.

Unlike the magnitude-limited BGS sample of relatively nearby galaxies, the
DESI LRG sample is selected with a comparatively complex set of magnitude and
color cuts, creating incompleteness that may depend on any combination of
galaxy color, stellar mass, and redshift. However, the LRG sample covers a
volume 20 times larger than that of the BGS sample, and will contain a much
higher number density of spectroscopic redshifts at ${0.4<z<1}$ than any
previous spectroscopic galaxy redshift survey, making it a good sample for
statistical studies with negligible sample and cosmic variance uncertainties.
The rarity and associated low number density of massive galaxies $(\sim
2\times 10^{-5}\ {\rm Mpc}^{-3}$ for $\log(M_{*}/{\rm M}_{\odot})>11.5)$ means
that such large volumes are essential for obtaining sample sizes large enough
for statistically significant measurements of the high-mass end of the galaxy
stellar mass function (SMF) and the SHMR of LRGs.

Existing studies of the galaxy–halo relationship for “DESI-like” LRG samples
include halo occupation distribution (HOD) modeling (e.g., Seljak 2000;
Berlind & Weinberg 2002; Bullock et al. 2002; Berlind et al. 2003; Zheng et
al. 2007; Zheng & Weinberg 2007) in Zhou et al. (2021), and semi-analytic
modeling in Hernández-Aguayo et al. (2021) with Galform (Cole et al. 2000),
which predicts absolute magnitudes with dust attenuation.

In its simplest form, the HOD model provides a statistical description for how
galaxies occupy dark matter halos solely as a function of halo mass. Zhou et
al. (2021) fit the projected clustering of “DESI-like” LRGs measured in five
redshift bins from $0.4<z<0.9$ with a five-parameter HOD model plus a sixth
nuisance parameter to account for photometric redshift uncertainties. They
find similar HOD parameters at $0.4<z<0.8$, and statistically significant
differences in model parameters for only the highest redshift bin
$(0.8<z<0.9)$. While Zhou et al. (2021) demonstrate that clustering
measurements using photometric redshifts are sufficient to constrain the HOD
parameters of “DESI-like” LRGs, the standard HOD framework they use does not
accommodate galaxy populations with halo occupation fractions that do not
reach unity for the most massive halos.

Hernández-Aguayo et al. (2021) apply the Galform semi-analytic model (SAM) to
the Planck-Millennium cosmological $N$-body simulation (Baugh et al. 2018).
They then convert predicted absolute magnitudes to apparent magnitudes at
various redshift snapshots, enabling the DESI LRG target selection to be
applied directly to Galform mock galaxy catalogs to select mock LRG samples.
Hernández-Aguayo et al. (2021) find that the DESI LRG selection criteria
exclude a small but important fraction of the most massive galaxies
($\log(M_{*}/{\rm M}_{\odot})>11.15$). Consequently, their model predicts that
the halo occupation fraction of LRGs does _not_ reach unity for the most
massive halos, and actually drops with increasing mass, indicative of a non-
trivial LRG–halo connection that is not modeled well with a standard HOD. By
comparing the HOD and subhalo mass functions of stellar mass-selected mock
galaxies against those of mock LRG samples, Hernández-Aguayo et al. (2021)
show that the DESI LRG selection cuts likely affect the selection of subhalos
populated by LRGs, i.e., (sub)halo mass is by itself insufficient to determine
whether a subhalo hosts a LRG.

Subhalo abundance matching (SHAM; e.g., Vale & Ostriker 2004; Conroy et al.
2006; Behroozi et al. 2010; Reddick et al. 2013) is an empirical technique for
assigning galaxies to dark matter halos in numerical $N$-body simulations by
assuming a correlation between a galaxy property—usually luminosity or stellar
mass—and halo property such as mass or circular velocity. In the simplest
application of SHAM, mock galaxy catalogs are constructed to reproduce the
number density and luminosity or stellar mass function (SMF) of a target
dataset with a single free parameter to allow scatter in the galaxy
property–halo property correlation.

Extensions to the SHAM framework to incorporate dependencies between
additional galaxy and halo properties are broadly referred to as conditional
abundance matching (CAM; e.g., Hearin et al. 2014; Zentner et al. 2014). In
addition to a primary galaxy property–halo property correlation, CAM models
assume a correlation between secondary galaxy and halo properties at a fixed
value of the primary property. Age distribution matching (Hearin & Watson
2013) is a form of CAM that equates galaxy color or a similar property with a
proxy for the age of dark matter halos. Unlike the standard HOD framework, in
which the statistical relationship between galaxies and dark matter halos is
solely a function of stellar and halo mass, CAM can naturally accommodate
galaxy assembly bias, the dependence of galaxy properties (besides stellar
mass) on the mass accretion history of their host halos (e.g., Zentner et al.
2014; Wechsler & Tinker 2018, and references therein).

The SHAM framework has been used to study the dependence of sample
completeness on redshift and stellar mass for LRGs from the Baryon Oscillation
Spectroscopic Survey (BOSS; Dawson et al. 2013), which obtained spectroscopic
redshifts of 1.5 million galaxies with $\log(M_{*}/{\rm M}_{\odot})>11$ to
$z\sim 0.7$. BOSS contains two color- and magnitude-selected samples of
massive galaxies ($\log(M_{*}/{\rm M}_{\odot})>11$; Reid et al. 2015): the
LOWZ sample of LRGs at $0.15<z<0.43$, and the approximately stellar mass-
limited constant mass (CMASS) sample, which includes galaxies of all colors at
$0.43<z<0.8$.

Saito et al. (2016) use SHAM to construct $z\sim 0.5$ mock galaxy catalogs for
the BOSS CMASS sample with and without added assembly bias effects. They use
the Stripe 82 Massive Galaxy Catalog (S82-MGC; Bundy et al. 2015) to replicate
the total galaxy SMF above $\log(M_{*}/{\rm M}_{\odot})>10.5$ over
$0.43<z<0.7$ and assign galaxies to halos in the MultiDark simulation (Riebe
et al. 2013). Saito et al. (2016) find that assembly bias does factor into the
galaxy–halo connection for high-mass galaxies, i.e., the SHMR for these
galaxies has some dependence on galaxy color, and should not be inferred from
the clustering signal without any consideration of color.

Yu et al. (2022) model LRGs from BOSS and eBOSS (Dawson et al. 2016) at
${0.2<z<1.0}$ with a SHAM framework that includes two additional free
parameters to account for redshift uncertainty and sample incompleteness.

SHAM has also been used to model “DESI-like” samples at low redshift. Safonova
et al. (2021) use SHAM and age distribution matching to create $z\sim 0.1$
mock galaxy catalogs representative of the DESI BGS sample with $r$-band
luminosities and ${g-r}$ colors. The low redshift range of the BGS sample
allows them to utilize spectroscopic redshifts from the Sloan Digital Sky
Survey (SDSS; York et al. 2000) and Galaxy and Mass Assembly (GAMA; Loveday et
al. 2012) project.

In this work, we use SHAM and age distribution matching to create magnitude-
limited, mock galaxy catalogs at multiple redshifts within the redshift range
of the DESI LRG sample. Two distinct target selection algorithms were
considered for the DESI LRG sample: an optical selection that uses $r-z$
color, and an infrared selection that uses $r-W1$222$W1$ is the 3.4 micron
band of the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010)
color. We select mock LRG samples from our magnitude-limited mocks based on
both the optical and IR DESI target selections that match the number density
and two-dimensional color–magnitude space distribution of each DESI LRG target
sample. We then predict the clustering signal and halo occupation statistics
of these samples as a function of redshift.

Our method is novel approach to modeling DESI LRGs that complements existing
SAM and HOD models by utilizing the full photometric samples that are the
basis of DESI LRG target selection. We offset the precision of training data
lost to photometric redshift errors by driving down cosmic variance
uncertainties with complete photometric samples from an unprecedented survey
volume. Finally, this work offers a comparative study of the samples selected
by the optical and IR selection algorithms considered for DESI LRGs.

The structure of this paper is as follows. In §2 describe the cosmological
simulation and photometric galaxy samples used in this work. §3 describes our
modeling procedure and the two-point statistics used to constrain our models.
In §4, we present the predicted properties of LRG samples, and in §5 we
summarize the conclusions of this work. Where applicable, we assume the
cosmological parameter values of Planck Collaboration et al. (2016):
$h=0.6777$ and $\Omega_{\rm m}=0.307115$, under the assumption of a flat,
$\Lambda$CDM cosmological model.

##  2\. Simulations and data

In this section we describe the cosmological simulation, halo finder, and
associated halo properties used for our models. We also present the data from
which we select parent galaxy samples for training our models, as well as the
DESI LRG target selection functions.

###  2.1. Simulations

We use halo catalogs and merger histories obtained with the publicly available
ROCKSTAR phase-space temporal halo finder (Behroozi et al. 2012b) for the
MultiDark Planck 2 (MDPL2) simulation333www.cosmosim.org (Klypin et al. 2016).
MDPL2 assumes Planck cosmology ($h=0.6777$, $\Omega_{\rm m}=0.307115$; Planck
Collaboration et al. 2016) and evolves $3840^{3}$ dark matter particles in a 1
$h^{-1}\ {\rm Gpc}$ cubic volume, beginning at $z=120$. The particle
resolution is $1.51\times 10^{9}\ h^{-1}\ {\rm M}_{\odot}$. In total 126
snapshots are available between $z\sim 15$ and $z=0$. For this work we use
three snapshots at $z=0.425$, 0.523, and 0.628.

The ROCKSTAR halo finder is designed to preserve particle–halo membership and
identify accurate halo merger trees across multiple time steps of a
simulation. For MDPL2, ROCKSTAR halo catalogs are mass-complete for halos
(including subhalos) above ${M_{\rm vir}\gtrsim 11.4\times 10^{11}\ h^{-1}\
{\rm M}_{\odot}}$.

####  2.1.1 Subhalo abundance matching

In SHAM modeling, mock galaxies are assigned to dark matter halos by
exploiting the correlation between some galaxy property—usually stellar mass
or luminosity—and a halo property such as virial mass or circular velocity.
Circular velocity is defined as ${v_{\rm circ}(r,z)\equiv\sqrt{GM(<r,z)/r}}$,
where $M(<r,z)$ is the enclosed mass within radius $r$ at redshift $z$. SHAM
has been tested with several versions of halo circular velocity. In this work
we use $v_{\rm peak}$, the peak value of maximum circular velocity (${v_{\rm
max}(z)={\rm max}\\{v_{\rm circ}(r,z)\\}}$) achieved throughout a halo’s
entire assembly history.

####  2.1.2 Age distribution matching

Figure 1.—  Mass accretion histories and $z_{\rm starve}$ values of four
randomly selected halos from the $z_{\rm sim}=0.425$ snapshot of the MDPL2
simulation. The most massive halo (purple dash-dotted line) has the earliest
starvation redshift ($z_{\rm starve}=4.63$), corresponding to when its mass
first reaches the characteristic value ${M_{\rm char}=10^{12}\ h^{-1}{\rm
M}_{\odot}}$. The other three halos also reach $M_{\rm char}$ (at later
redshifts), but in accordance with Eq. 2 the redshift at which this occurs for
each halo is not necessarily the same as its $z_{\rm starve}$ value, e.g., the
least massive halo (orange solid line) has $z_{\rm starve}=1.36$ although its
mass doesn’t exceed $M_{\rm char}$ until $z\lesssim 1$. The shaded region
below each curve corresponds a halo’s mass accretion history before it reaches
$z_{\rm starve}$.

Age distribution matching assumes a correlation at fixed luminosity between
galaxy color and some proxy (at fixed model luminosity) for the age of the
halo in which each mock galaxy resides. Redder colors (i.e., older, quenched
galaxies) are generally assigned to older halos. Model colors are assigned at
fixed luminosity (in practice in narrow luminosity bins) because the galaxy
color distribution is highly dependent on luminosity.

We equate the cumulative distribution $\mathcal{D}_{\rm gal}$ of galaxies of
$K$-corrected color $\mathcal{C}$ at fixed absolute magnitude $M_{X}$ to the
cumulative distribution $\mathcal{D}_{\rm halo}$ of halo age proxy $A$ at
fixed model absolute magnitude:

$\mathcal{D}_{\rm gal}(<\mathcal{C}\ |\ M_{X})=\mathcal{D}_{\rm halo}(<A\ |\
{\rm model}\ M_{X}).$ (1)

In Equation 1 $(\mathcal{C},M_{X})={(r-z,M_{z})}$ OR ${(r-W1,M_{W1})}$ for
this work.

Implementations of age distribution matching at $z\sim 0$ using spectroscopic
galaxy redshifts from SDSS have used halo starvation redshift, $z_{\rm
starve}$, for the halo age proxy $A$ in Eq. 1 (Hearin & Watson 2013; Hearin et
al. 2014; Safonova et al. 2021). In general, $z_{\rm starve}$ represents the
redshift at which a galaxy loses its supply of cold gas, which leads to the
quenching of star formation and the reddening of the galaxy. Multiple physical
processes relevant to a halo’s assembly history can affect the value of
$z_{\rm starve}$ for a given halo, which Hearin & Watson (2013) incorporate
into the following definition:

$z_{\rm starve}\equiv\max\\{z_{\rm char},\,z_{\rm acc},\,z_{\rm form}\\}.$ (2)

In Equation 2:

  * •

$z_{\rm char}$ is either the redshift at which a halo’s mass first exceeds
some characteristic value, $M_{\rm char}$, or the redshift of the relevant
simulation snapshot ($z_{\rm sim}$) for halos that never achieve $M_{\rm
char}$;

  * •

$z_{\rm acc}$ is the redshift at which a subhalo accretes onto a parent halo
(for host halos ${z_{\rm acc}=z_{\rm sim}}$). Hearin & Watson (2013) follow
Behroozi et al. (2012c), which defines $z_{\rm acc}$ as the snapshot after
which a subhalo always remains a subhalo. They note that alternative
definitions, such as that of Wetzel et al. (2014), where $z_{\rm acc}$ is the
snapshot at which a subhalo has been identified as such for two consecutive
snapshots, have little impact on their results.

  * •

$z_{\rm form}$ is the “formation” redshift at which a halo transitions from
the fast to slow accretion regime.

We use same definition of $z_{\rm form}$ as Hearin & Watson (2013), motivated
by Wechsler et al. (2002):

$z_{\rm form}\equiv\frac{c_{\rm vir}}{4.1a_{0}}-1,$ (3)

where ${c_{\rm vir}=R_{\rm vir}/R_{\rm s}}$ is a halo’s concentration at the
time of observation, indicated by $a_{0}$. For host halos $a_{0}$ is the scale
factor of the relevant simulation snapshot, while for subhalos $a_{0}$ is the
scale factor at the time of accretion: ${z_{\rm acc}=1/{a_{0}}-1}$. $R_{\rm
vir}$ is the virial radius of a halo, and $R_{\rm s}$ is the NFW scale radius
(Navarro et al. 1997).

We adopt the value of ${M_{\rm char}=10^{12}\ h^{-1}{\rm M}_{\odot}}$ used in
Hearin & Watson (2013), who note that their results are insensitive to the
precise value of $M_{\rm char}$ used. The empirical and physical motivation
for $M_{\rm char}$ are described in detail in §6.3 of Hearin & Watson (2013).
Briefly, there is empirical support for a characteristic halo mass above which
star formation is highly inefficient: ${\sim 10^{12}\ h^{-1}{\rm M}_{\odot}}$
is the halo mass at which the SHMR peaks, falling off rapidly at higher halo
masses (Behroozi et al. 2013; Yang et al. 2012, 2013; Moster et al. 2013;
Watson & Conroy 2013), and Behroozi et al. (2012a) have shown that this mass
remains essentially constant throughout much of cosmic history.

We compute $z_{\rm starve}$ for all halos in our model from the publicly
available ROCKSTAR halo merger trees for MDPL2. Figure 1 shows sample halo
mass accretion histories and corresponding $z_{\rm starve}$ values for four
randomly selected halos from the $z_{\rm sim}=0.425$ snapshot of MDPL2.

###  2.2. Photometry and redshift estimates

We use publicly available catalogs from the ninth data release (DR9) of the
DESI Legacy Imaging Surveys444www.legacysurvey.org/dr9 (Dey et al. 2019). The
Legacy Surveys provide optical imaging in the $g$, $r$, and $z$ bands from a
combination of three public surveys: the DECam Legacy Survey (DECaLS; Flaugher
et al. 2015; Blum et al. 2016), the Beijing-Arizona Sky Survey (BASS; Zou et
al. 2017), and the Mayall $z$-band Legacy Survey (MzLS; Silva et al. 2016).
The Legacy Surveys also include four mid-infrared bands from the Wide-field
Infrared Survey Explorer (WISE; Wright et al. 2010), although only the 3.4
micron $W1$-band is relevant for DESI LRG target selection.

In total the Legacy Surveys cover $14,000$ square degrees visible from the
northern hemisphere, comprised of two contiguous regions within the northern
and southern galactic caps. To avoid effects from systematic differences among
data from the three component optical surveys, we limit our study to the
approximately 9,000 square degrees covered by DECaLS.

Zhou et al. (2021) compute photometric redshifts for the full catalog of
DECaLS DR7 objects using the random forest regression machine learning
algorithm in Scikit-Learn (Pedregosa et al. 2011) and a “truth” dataset of
spectroscopic and many-band photometric redshifts for objects within DR7. They
quantify the accuracy of their photometric redshifts with the normalized
median absolute deviation (NMAD; Dahlen et al. 2013): $\sigma_{\rm
NMAD}=1.48\times{\rm median}(|\Delta z|/(1+z_{\rm spec}))$, where $\Delta
z=z_{\rm phot}-z_{\rm spec}$ and $z_{\rm spec}$ are the redshift truth values
used to train the random forest algorithm, and report $\sigma_{\rm
NMAD}=0.021$ for LRGs. Their outlier rate for LRGs is 1%, where outliers are
objects with $|\Delta z|>0.1\times(1+z_{\rm spec})$. Additionally, Zhou et al.
(2021) estimate their redshifts are accurate for objects with apparent
$z$-band magnitude $z<21$, well beyond the $z<20.7$ cut we use to select our
target galaxy samples, described in §2.4 below.

###  2.3. DESI LRG target selection

The DESI LRG target sample is intended to serve as a cosmological tracer
spanning the redshift range $\sim 0.4<z\lesssim 1.0$. The sample lies between
the low-redshift Bright Galaxy Survey (BGS; Hahn et al. 2022) tracer sample at
$z\lesssim 0.4$, and the Emission Line Galaxy (ELG; Raichoor et al. 2022)
sample, optimized to trace the density field over the approximate range
${1.0<z<1.6}$.

Two different selection algorithms were considered for the DESI LRG sample: an
optical selection function based on $z$-band magnitude and ${r-z}$ color, and
an infrared selection function based on $W1$-band magnitude and ${r-W1}$
color. Both selections are tuned to yield a constant LRG target density of
$\sim 600$ objects per square degree and a comoving number density around
${5\times 10^{-4}\ h^{3}\ {\rm Mpc}^{-3}}$ at ${0.4<z<0.8}$. Both the optical
and IR selections were tested in DESI’s Survey Validation (SV; DESI et al., in
prep.) observations. Based largely on the calibration of $W1$-band imaging,
the DESI Main Survey uses exclusively the IR selection. A complete description
of DESI LRG target selection is given in Zhou et al. (2023). Here we cover the
details most relevant for this work.

Due to slight differences in photometry among BASS, MzLS, and DECaLS, the
optical and IR DESI LRG target selections use slightly different cuts for the
north (BASS and MzLS) and south (DECaLS) galactic caps. As this work uses $g$,
$r$, and $z$ magnitudes from DECaLS, we use the corresponding optical LRG
target selection cuts (Zhou et al. 2020):

$\displaystyle z-W1>0.8\times(r-z)-0.6,$ (4a)
$\displaystyle\begin{split}&((g-W1>2.6)\ {\rm AND}\ (g-r>1.4))\ {\rm OR}\\\
&\phantom{00}(r-W1>1.8),\end{split}$ (4b)
$\displaystyle\begin{split}&(r-z>0.45\times(z-16.83))\ {\rm AND}\ (r-z>0.7)\\\
&\phantom{00}{\rm AND}\ (r-z>0.19\times(z-13.80)),\end{split}$ (4c)
$\displaystyle z_{\rm fiber}<21.5.$ (4d)

The relevant IR LRG target selection cuts for this work are (Zhou et al.
2023):

$\displaystyle z-W1>0.8\times(r-z)-0.6,$ (5a) $\displaystyle(g-W1>2.9)\ {\rm
OR}\ (r-W1>1.8),$ (5b)
$\displaystyle\begin{split}&((r-W1>1.8\times(W1-17.14))\ {\rm AND}\\\
&\phantom{00}(r-W1>W1-16.33))\ {\rm OR}\ (r-W1>3.3),\end{split}$ (5c)
$\displaystyle z_{\rm fiber}<21.6.$ (5d)

Equations 4a and 5a are designed to reject stars, Eqs. 4b and 5b remove blue
and low-redshift objects, Eqs. 4c and 5c are color-dependent magnitude limits
that select only the most luminous objects at a given redshift, and $z_{\rm
fiber}$ in Eqs. 4d and 5d is the expected $z$-band flux within a DESI fiber.
All magnitudes in Eqs. 4 and 5 use the AB system, and are corrected for
galactic extinction using the relevant MW_TRANSMISSION values from the Legacy
Surveys DR9.

###  2.4. Parent galaxy samples

Figure 2.— Distribution of $K$-corrected $z$-band (left column) and $W1$-band (right column) absolute magnitudes for the three redshift bins in our study: ${0.4<z_{\rm phot}<0.5}$ (top), ${0.5<z_{\rm phot}<0.6}$ (middle), and ${0.6<z_{\rm phot}<0.7}$ (bottom). Each panel shows the distribution of all galaxies with $z<20.7$ (solid gray), as well as the distributions of optical (hatched purple; see Eq. 4) and IR (solid orange; see Eq. 5) DESI LRG targets. The leftmost dashed black line in each panel is the absolute magnitude cut that defines each absolute magnitude-limited parent galaxy sample for our models. At fainter magnitudes (to the left of the cut in each panel) the DECaLS sample is incomplete. The three dotted black lines in each panel denote the luminosity bins used to constrain model parameters (see §3.2 and Tables 2 and 3).  |
---|---
|
Figure 3.—  Color–magnitude diagrams of the parent galaxy samples and DESI LRG target samples (overlaid white contours) in the $0.4<z_{\rm phot}<0.5$ redshift bin. The left column shows optical ($K$-corrected $r-z$ color versus $M_{z}$ magnitude), while the right column shown IR ($K$-corrected $r-W1$ color versus $M_{W1}$ magnitude). The top (bottom) row shows the distribution of optical (IR) LRGs in each color–magnitude space. The boxed regions labeled “A” and “B” are referenced in §4.1 below.  Table 1 Bitmasks applied to photometry for parent galaxy sample selection. Additional details at legacysurvey.org/dr9/bitmasks.  MASKBIT | Description
---|---
5, 6, 7 | bad pixel in all of a set of overlapping $g$, $r$, or $z$-band
| images
8, 9 | bad pixel in a WISE $W1$ or $W2$ bright star mask
11 | pixel within locus of a radius-magnitude relation for
| GaiaaaGaia Collaboration (2018) DR2 stars to $G<16$
12 | pixel in a Siena Galaxy AtlasbbMoustakas et al. (2021) large galaxy
13 | pixel in a globular cluster
FITBIT | Description
6 | source is a medium-bright star
7 | Gaia sourceccGaia Collaboration et al. (2016)
8 | Tycho-2 starddHøg et al. (2000)
Table 2Properties of parent galaxy samples and corresponding simulation snapshot redshifts ($z_{\rm sim}$).  Redshift bin | $z_{\rm sim}$ | Luminosity cutaa$K$-correction redshifts are rounded to two decimal places for clarity, e.g., ${}^{0.43}M_{z}$ indicates absolute $z$-band magnitudes are $K$-corrected to $z_{\rm sim}=0.425$. | $N_{\rm gal}$ | $n_{\rm eff}$ | Included fraction of DESI LRG targets
---|---|---|---|---|---
$[\times 10^{-3}\ h^{3}\ {\rm Mpc}^{-3}]$ | Optical selection | IR selection
$0.4<z_{\rm phot}<0.5$ | 0.425 | ${}^{0.43}M_{z}<-21.60$ | 8,314,309 | 5.60 | 0.992 | 0.982
${}^{0.43}M_{W1}<-22.25$ | 7,565,153 | 4.33 | 0.992 | 0.992
$0.5<z_{\rm phot}<0.6$ | 0.523 | ${}^{0.52}M_{z}<-21.60$ | 7,804,346 | 2.86 | 0.992 | 0.988
${}^{0.52}M_{W1}<-22.85$ | 4,909,857 | 2.01 | 0.992 | 0.992
$0.6<z_{\rm phot}<0.7$ | 0.628 | ${}^{0.63}M_{z}<-21.85$ | 6,548,126 | 1.46 | 0.994 | 0.990
${}^{0.63}M_{W1}<-23.15$ | 4,758,470 | 1.12 | 0.993 | 0.994
Table 3Luminosity bins used to constrain the magnitude dependence of model parameters.  $0.4<z_{\rm phot}<0.5$ | $0.5<z_{\rm phot}<0.6$ | $0.6<z_{\rm phot}<0.7$
---|---|---
Luminosity bin | $N_{\rm gal}$ | Luminosity bin | $N_{\rm gal}$ | Luminosity bin | $N_{\rm gal}$
$z$-band
$-21.60>{{}^{0.43}M_{z}}>-21.85$ | 2,458,920 | $-21.60>{{}^{0.52}M_{z}}>-21.85$ | 2,304,853 | $-21.85>{{}^{0.63}M_{z}}>-22.10$ | 2,040,336
$-21.85>{{}^{0.43}M_{z}}>-22.10$ | 1,880,295 | $-21.85>{{}^{0.52}M_{z}}>-22.10$ | 1,732,820 | $-22.10>{{}^{0.63}M_{z}}>-22.35$ | 1,512,940
$-22.10>{{}^{0.43}M_{z}}>-22.35$ | 1,317,269 | $-22.10>{{}^{0.52}M_{z}}>-22.35$ | 1,205,545 | $-22.35>{{}^{0.63}M_{z}}>-22.60$ | 1,039,573
${{}^{0.43}M_{z}}<-22.35$ | 1,919,594 | ${{}^{0.52}M_{z}}<-22.35$ | 1,855,014 | ${{}^{0.63}M_{z}}<-22.60$ | 1,401,662
$W1$-band
$-22.25>{{}^{0.43}M_{W1}}>-22.55$ | 2,324,816 | $-22.85>{{}^{0.52}M_{W1}}>-23.15$ | 1,771,554 | $-23.15>{{}^{0.63}M_{W1}}>-23.45$ | 1,634,860
$-22.55>{{}^{0.43}M_{W1}}>-22.85$ | 1,837,859 | $-23.15>{{}^{0.52}M_{W1}}>-23.45$ | 1,245,107 | $-23.45>{{}^{0.63}M_{W1}}>-23.75$ | 1,273,814
$-22.85>{{}^{0.43}M_{W1}}>-23.15$ | 1,277,192 | $-23.45>{{}^{0.52}M_{W1}}>-23.75$ | 754,154 | $-23.75>{{}^{0.63}M_{W1}}>-24.05$ | 798,012
${{}^{0.43}M_{W1}}<-23.15$ | 1,565,914 | ${{}^{0.52}M_{W1}}<-23.75$ | 693,491 | ${{}^{0.63}M_{W1}}<-24.05$ | 685,799

A primary goal of this work is to create mock galaxy catalogs that are both
statistically complete and represent a superset of the color–magnitude space
occupied by DESI LRG targets. Selection of DESI LRG targets is based entirely
on $g$, $r$, $z$, and $W1$ apparent magnitudes (see §2.3), so we would ideally
create mock catalogs where every mock galaxy has an apparent magnitude in each
of these bands. SHAM, however, exploits the correlation between some physical
halo property (e.g., circular velocity) and a physical galaxy property
independent of redshift (e.g., luminosity). We therefore train our mock
catalogs on galaxy samples that are complete to an _absolute_ magnitude
threshold that includes all DESI LRG targets (in the relevant redshift bin;
see below).

To select suitable parent galaxy samples, we first apply an apparent $z$-band
magnitude cut of $z<20.7$, and take an additional step to remove stars by
excluding catalog sources with ${\rm TYPE}=\texttt{PSF}$. We also apply the
masks described in Table 1, which are provided with DECaLS DR9, to remove
sources affected by bad pixels or contamination from bright stars. Finally, a
geometric mask is applied to ensure complete angular coverage by the catalogs
of random points provided with DR9 (see §3.1). The resulting sample contains
$\mathcal{O}(10^{8})$ galaxies, sufficient to divide it into redshift bins and
maintain low statistical error.

We initially tested six redshift bins of width $\Delta z_{\rm phot}=0.1$
between $z_{\rm phot}=0.4$ and $z_{\rm phot}=1.0$, but found that DECaLS
photometry is only deep enough to apply our model up to $z_{\rm phot}\sim
0.7$. At $z_{\rm phot}\gtrsim 0.7$ the data are incomplete above the absolute
magnitude threshold that encompasses DESI LRG targets. We therefore limit our
study to three redshift bins of $\Delta z_{\rm phot}=0.1$ within $0.4<z_{\rm
phot}<0.7$.

For each redshift bin we compute $z$\- and $W1$-band absolute magnitudes using
photometric redshifts to obtain distance moduli. We $K$-correct absolute
magnitudes to the redshift of the relevant simulation snapshot, $z_{\rm sim}$
(see Table 2) with the IDL package kcorrect (Blanton & Roweis 2007) using
DECam $g$, $r$, and $z$ and WISE $W1$ and $W2$ filter responses. For each
redshift bin we use the simulation snapshot closest to the median $z_{\rm
phot}$ of the data, e.g., galaxies in the ${0.4<z_{\rm phot}<0.5}$ bin are
$K$-corrected to $z_{\rm sim}=0.425$.

The final step in selecting parent galaxy samples is to identify $z$\- and
$W1$-band absolute magnitude cuts in each redshift bin that yield complete
samples which also include the full absolute magnitude range of DESI LRG
targets in that bin. Figure 2 shows $K$-corrected absolute magnitude
distributions for all $z<20.7$ DECaLS galaxies in each redshift bin. For each
bandpass and redshift bin combination, we identify an absolute magnitude cut
(dashed black lines in Figure 2) that eliminates fainter galaxies where DECaLS
becomes incomplete, while preserving $\gtrsim 99\%$ of DESI LRG targets (solid
orange and hatched purple histograms).

Table 2 lists the details of each parent galaxy sample, including the relevant
absolute magnitude cut, sample size, effective number density, and the
included fractions of IR and optical DESI LRG targets. Besides enforcing
statistically complete parent samples, these magnitude cuts also eliminate
galaxies with larger photometric redshift errors, increasing the accuracy of
clustering measurements (see §3.1).

Figure 3 shows example optical ($r-z$ versus $M_{z}$) and IR ($r-W1$ versus
$M_{W1}$) color–magnitude diagrams of the magnitude-limited parent galaxy
samples for the $0.4<z_{\rm phot}<0.5$ redshift bin. Also shown in Figure 3
are the color–magnitude distributions of optical and IR DESI LRG targets.

##  3\. Modeling

In this section we describe our modeling procedure and the two-point
statistics we use to constrain the model parameters.

###  3.1. Projected correlation functions

One goal of this work is to exploit the completeness and enormous volume of
DECaLS data, which comes at the expense of the precision clustering
measurements achievable with spectroscopic redshifts. Zhou et al. (2021)
demonstrate the constraining power of the projected correlation function,
$\omega_{\rm p}(r_{\rm p})$, of DECaLS galaxies computed with line-of-sight
distances derived from photometric redshifts (they use this statistic to
compute HOD parameters of “DESI-like” LRGs selected from DECaLS DR7). The
projected correlation function conveys 3D correlation function, $\xi(r)$,
integrated along the line-of-sight, effectively eliminating the effects of
radial distance uncertainty due to photometric redshift errors:

$\omega_{\rm p}(r_{\rm p})\equiv\int_{-\pi_{\rm max}}^{\pi_{\rm
max}}\\!\xi(r_{\rm p},\pi)\,d\pi=2\\!\int_{0}^{\pi_{\rm max}}\\!\xi(r_{\rm
p},\pi)\,d\pi,$ (6)

where $r_{\rm p}$ is the projection of $r$ into the plane perpendicular to the
line-of-sight separation, $\pi$.

We use the Corrfunc package (Sinha & Garrison 2017, 2019) to calculate
$\omega_{\rm p}(r_{\rm p})$ for both our target galaxy samples and mock
catalogs in 19 logarithmic bins between ${r_{\rm p}>0.1\ h^{-1}\ {\rm Mpc}}$
and ${r_{\rm p}\lesssim 44\ h^{-1}\ {\rm Mpc}}$. We also compute $\omega_{\rm
p}(r_{\rm p})$ in additional bins at ${r_{\rm p}<0.1\ h^{-1}\ {\rm Mpc}}$ but
do not use these measurements for model fitting.

As our data samples are confined to narrow redshift bins of width ${\Delta
z_{\rm phot}=0.1}$, photometric redshift errors will cause some galaxies that
belong to a given redshift bin to scatter into an adjacent bin and be excluded
from the calculation of $\omega_{\rm p}(r_{\rm p})$ for their true bin. To
account for this we adopt the method used by Zhou et al. (2021) (see their
Figure 8): we use the Landy-Szalay estimator (Landy & Szalay 1993) for the
cross-correlation of two samples, $D_{1}$ and $D_{2}$:

$\omega_{\rm p}(r_{\rm p})=\\!\\!\sum_{-\pi_{\rm max}}^{\pi_{\rm
max}}\\!\\!\left(\frac{D_{1}D_{2}-D_{1}R_{2}-D_{2}R_{1}}{R_{1}R_{2}}+1\right),$
(7)

Each term of Eq. 7 denotes pair counts between two samples, where $D$ and $R$
respectively indicate samples of data (i.e., galaxies) and random points with
the same angular and redshift distributions as the corresponding data sample.
Here $D_{1}$ is all galaxies within a given redshift bin: ${z_{\rm min}<z_{\rm
phot}<z_{\rm max}}$, where $z_{\rm min}$ and $z_{\rm max}$ are the limits of
the bin, while $D_{2}$ is all galaxies within a wider redshift range defined
by ${(z_{\rm min}-\pi_{\rm max})<z_{\rm phot}<(z_{\rm max}+\pi_{\rm max})}$,
where ${\pi_{\rm max}=150~{}h^{-1}\ {\rm Mpc}}$. We verify our implementation
of this method with Corrfunc by reproducing the projected correlation
functions of DECaLS LRGs from Zhou et al. (2021) (see their Figure 9) using
different clustering code.

DECaLS data include catalogs of random points with the same angular sky
coverage and mask information as the survey footprint, which we use to
construct our random samples. We use 20 times as many random as data points
for each galaxy sample, and draw redshifts for random points from the redshift
distribution of the corresponding data sample.

To measure $\omega_{\rm p}(r_{\rm p})$ of our mock catalogs we take advantage
of the Corrfunc theory module, which can quickly calculate the autocorrelation
function of a sample within a periodic volume using analytic random pair
counts. We confirmed that this method produces the expected result by
calculating $\omega_{\rm p}(r_{\rm p})$ of several mock catalogs directly from
pair counts between mock galaxies and catalogs of random points constructed
for the simulation volume.

###  3.2. Jackknife error estimation and goodness-of-fit

To estimate the uncertainty of the $\omega_{\rm p}(r_{\rm p})$ measurements of
our target galaxy samples we use healpy555healpix.sourceforge.net (Górski et
al. 2005; Zonca et al. 2019) with $N_{\rm side}=6$ to divide the angular sky
coverage of each galaxy sample into $N_{\rm jk}$ regions of roughly equal
area, suitable for jackknife resampling. We then measure $\omega_{\rm
p}(r_{\rm p})$ in each jackknife sample, where each jackknife sample consists
of the entire galaxy sample with one jackknife region removed, and compute the
covariance matrix as follows:

${\rm Cov}_{ij}=\frac{N_{\rm jk}-1}{N_{\rm jk}}\sum^{N_{\rm
jk}}_{\ell=1}\left(\omega^{\ell}_{i}-\overline{\omega}_{i}\right)\left(\omega^{\ell}_{j}-\overline{\omega}_{j}\right),$
(8)

where $\omega^{\ell}_{i}$ and $\omega^{\ell}_{j}$ are $\omega_{\rm p}(r_{\rm
p})$ of the $\ell$th jackknife region for the $i$th and $j$th $r_{\rm p}$
bins, respectively, and $\overline{\omega}_{i}$ and $\overline{\omega}_{j}$
are the mean values of $\omega_{\rm p}(r_{\rm p})$ across all jackknife
regions for the $i$th and $j$th $r_{\rm p}$ bins, respectively.

With the covariance matrix in hand we quantify how successful any instance of
our model is at fitting the projected correlation function of the data by
computing $\chi^{2}$ per degree of freedom $(\chi^{2}_{\nu})$:

$\chi^{2}_{\nu}=\frac{1}{\nu}\sum_{i=1}^{N_{r_{\rm p}}}\sum_{j=1}^{N_{r_{\rm
p}}}\left(\omega_{i}-\omega^{\rm mod}_{i}\right)\\!\left({\rm
Cov}^{-1}\right)_{ij}\\!\left(\omega_{j}-\omega^{\rm mod}_{j}\right)\\!,$ (9)

where $N_{r_{\rm p}}$ is the number of $r_{\rm p}$ bins used for fitting,
$\nu$ is equal to $N_{r_{\rm p}}$ minus the number of free model parameters,
$\omega_{i}$ and $\omega_{j}$ are the data $\omega_{\rm p}(r_{\rm p})$ values
in the $i$th and $j$th $r_{\rm p}$ bins, respectively, and $\omega^{\rm
mod}_{i}$ and $\omega^{\rm mod}_{j}$ are the $\omega_{\rm p}(r_{\rm p})$
values of the relevant mock catalog in the $i$th and $j$th $r_{\rm p}$ bins,
respectively.

###  3.3. Luminosity assignment

Figure 4.— Luminosity-binned projected clustering of the $W1$-band absolute
magnitude-limited parent galaxy samples and corresponding mock galaxy
catalogs, with different colors and line styles representing different
luminosity bins as denoted by the legend in each panel. Luminosity bins are
offset by 0.15 dex for clarity. Each panel shows a different redshift bin. The
shaded region at $r_{\rm p}<0.1\ h^{-1}\ {\rm Mpc}$ denotes measurements not
used for modeling.

Figure 5.— Same as Figure 4 but for the $z$-band parent galaxy samples and corresponding mock galaxy catalogs.  Figure 6.— Projected clustering of the $z$-band (top panel) and $W1$-band (bottom panel) absolute magnitude-limited parent galaxy samples and corresponding mock galaxy catalogs, with different colors and line styles representing different redshift bins as denoted by the legend in each panel. Redshift bins are offset by 0.15 dex for clarity.  Table 4Best-fit values for SHAM model parameters.  Redshift bin | Model | $\sigma_{\rm mag}$ | $s_{\rm LOS}$ | $\sigma_{{\rm LOS},0}$
---|---|---|---|---
$0.4<z_{\rm phot}<0.5$ |  $W1$-band | 0.66 | -22.3 | -439.2
$z$-band | 0.57 | -11.4 | -181.5
$0.5<z_{\rm phot}<0.6$ |  $W1$-band | 0.76 | 21.3 | 584.0
$z$-band | 0.76 | 60.9 | 1435.6
$0.6<z_{\rm phot}<0.7$ |  $W1$-band | 0.80 | 31.2 | 839.4
$z$-band | 0.77 | 50.8 | 1248.4

We use SHAM to create mock galaxy populations with the same number density and
luminosity distribution of the parent galaxy sample by assuming the following
relation:

$n_{\rm eff}(<M)=n_{\rm h}(>v_{\rm peak}),$ (10)

i.e., the (effective) number density of galaxies, $n_{\rm eff}$, of magnitude
$M$ or brighter equals the number density of halos, $n_{\rm h}$, with circular
velocity $v_{\rm peak}$ or greater.

To assign absolute magnitudes to halos in each redshift bin we implement the
following procedure:

  1. 1.

Compute for each parent galaxy sample the effective galaxy number density as a
function of absolute magnitude in band ${X\in\\{z,W1\\}}$, $n_{\rm
eff}(<M_{X})$. We compute $n_{\rm eff}(<M_{X})$ for each parent galaxy sample
as follows. For each galaxy in the sample we calculate an effective volume,
$V_{\rm eff}$:

${V_{\rm eff}(M_{X})=f_{\Omega}\left(V_{\rm max}(M_{X})-V_{\rm min}\right)},$
(11)

where $f_{\Omega}$ is the fractional solid angle covered by the parent galaxy
sample, and $V_{\rm min}$ is the comoving volume of the lower limit of the
redshift bin. $V_{\rm max}(M_{X})$ is the comoving volume of either the upper
limit of the redshift bin, or of the maximum possible redshift a galaxy of
magnitude $M_{X}$ could have and still be observed at the magnitude limit of
the sample, whichever is smaller. The effective galaxy number density is the
sum of the inverse of $V_{\rm eff}(M_{X})$ over all galaxies in the sample:

$n_{\rm eff}(<M_{X})=\sum_{i}\left[V^{i}_{\rm eff}(M_{X})\right]^{-1}.$ (12)

  2. 2.

Assign absolute magnitudes to mock galaxies with no scatter in the
luminosity–$v_{\rm peak}$ relation according to Eq. 10. For each halo we find
the cumulative number density $n_{\rm h}(>v_{\rm peak})$ corresponding to its
value of $v_{\rm peak}$. We then assign to each halo a mock galaxy with the
absolute magnitude $M_{X}$ at which the effective number density $n_{\rm
eff}(<M_{X})$ of the parent galaxy sample equals $n_{\rm h}(>v_{\rm peak})$.

  3. 3.

To incorporate magnitude-dependent scatter into the luminosity–$v_{\rm peak}$
relation, we assign to each mock galaxy a new absolute magnitude,
$M^{\prime}_{X}$, where $M^{\prime}_{X}$ is drawn from a Gaussian distribution
centered at $M_{X}$ with width $\sigma_{\rm mag}$. $\sigma_{\rm mag}$ is
proportional to the absolute value of $M_{X}$, and the constant of
proportionality is a free parameter in the model. We then rank order all mock
galaxies (including their $v_{\rm peak}$ values) by $M^{\prime}_{X}$, rank
order the original distribution of mock magnitudes $M_{X}$, and assign the
ordered original mock magnitude distribution to the mock galaxy catalog
ordered by $M^{\prime}_{X}$. This method incorporates scatter into the
luminosity–$v_{\rm peak}$ relation while replicating the target luminosity
function exactly in the resulting mock galaxy catalog.

  4. 4.

Scatter the positions of mock galaxies along one of the three axes of the
simulation volume to mimic the uncertainty in radial (line-of-sight) position
of our target galaxy samples due to photometric redshift errors. For each mock
galaxy we draw a “scattered” coordinate $x^{\prime}$ from a Gaussian
distribution of width $\sigma_{\rm LOS}$ centered at the galaxy’s original $x$
position. Mock galaxies that scatter out of the simulation volume of
$1~{}h^{-3}\ {\rm Gpc}^{3}$ are wrapped back in to preserve the periodic
boundary conditions, e.g., a mock galaxy at $x=25~{}h^{-1}\ {\rm Mpc}$ that
scatters to $x^{\prime}=-50~{}h^{-1}\ {\rm Mpc}$ is placed at $x=950~{}h^{-1}\
{\rm Mpc}$. We repeat this process for the other two simulation axes (i.e.,
$y$ and $z$).

  5. 5.

For each mock catalog, compute the projected correlation function,
$\omega_{\rm p}(r_{\rm p})$ (see §3.1), averaged over the three axes of the
simulation volume as described in the previous step. We then compute the
goodness-of-fit per degree of freedom, $\chi^{2}_{\nu}$ (see §3.2), of the
model fit to the mean $\omega_{\rm p}(r_{\rm p})$ of the relevant parent
galaxy sample. We measure $\chi^{2}_{\nu}$ for each mean $\omega_{\rm
p}(r_{\rm p})$ in bins of absolute magnitude ($M_{z}$ or $M_{W1}$).

  6. 6.

Repeat steps 3 through 5 for additional values of $\sigma_{\rm mag}$ and
$\sigma_{\rm LOS}$ as needed to minimize $\chi^{2}_{\nu}$ in each luminosity
bin. In practice we first coarsely sample a wide range of values of
$\sigma_{\rm mag}$ and $\sigma_{\rm LOS}$:

$\displaystyle 0\leq\frac{\sigma_{\rm mag}}{|M_{X}|}\leq 1.0,\ $
$\displaystyle\Delta\\!\\!\left[\frac{\sigma_{\rm mag}}{|M_{X}|}\right]=0.1$
(13a) $\displaystyle 0\leq\sigma_{\rm LOS}\leq 150~{}h^{-1}\ {\rm Mpc},\ $
$\displaystyle\Delta\sigma_{\rm LOS}=10~{}h^{-1}\ {\rm Mpc}.$ (13b)

We then more densely sample (using $\Delta(\sigma_{\rm mag}/|M_{X}|)=0.01$ and
${\Delta\sigma_{\rm LOS}=5~{}h^{-1}\ {\rm Mpc}}$) narrower ranges of both
parameters around the initial coarse-grained values that minimize
$\chi^{2}_{\nu}$.

  7. 7.

Finally, identify the value of $\sigma_{\rm mag}$ that coincides with the
region of $\sigma_{\rm mag}$–$\sigma_{\rm LOS}$ parameter space containing the
minimum $\chi^{2}_{\nu}$ across all luminosity bins, and parameterize the
dependence of $\sigma_{\rm LOS}$ on absolute magnitude, $M_{X}$, as follows:

$\sigma_{\rm LOS}(M_{X})=s_{\rm LOS}\,M_{X}+\sigma_{{\rm LOS},0}.$ (14)

The best-fit values of $\sigma_{\rm mag}$, $s_{\rm LOS}$, and $\sigma_{{\rm
LOS},0}$ are given in Table 4. The values of $s_{\rm LOS}$ and $\sigma_{{\rm
LOS},0}$ are determined by linear fits to $\sigma_{\rm LOS}^{i}$ versus
$\langle M^{i}_{X}\rangle$, where $\sigma_{\rm LOS}^{i}$ are the values that
minimize $\chi^{2}_{\nu}$ of $\omega_{\rm p}(r_{\rm p})$ at fixed $\sigma_{\rm
mag}$ 666We tested using $\sigma_{\rm mag}$ with linear dependence on absolute
magnitude instead of a constant value across all luminosity bins, and found
that fixing $\sigma_{\rm mag}$ while allowing $\sigma_{\rm LOS}$ to scale
linearly with luminosity reproduces $\omega_{\rm p}(r_{\rm p})$ (§3.1) better
than if both parameters have linear luminosity dependence. in the $i$th
luminosity bin, $M^{i}_{X}$.

The luminosity assignment stage of our modeling procedure involves three free
parameters, $\sigma_{\rm mag}$, $s_{\rm LOS}$, and $\sigma_{{\rm LOS},0}$,
which account for scatter in the luminosity–$v_{\rm circ}$ relation and
photometric redshift errors of our target galaxy samples. We constrain these
parameters by fitting the projected correlation functions of mock galaxy
catalogs created from our model to those of the corresponding parent galaxy
samples.

Figures 4 and 5 show $\omega_{\rm p}(r_{\rm p})$ for the $W1$\- and $z$-band
absolute magnitude-limited parent galaxy samples, respectively, and
corresponding mock galaxy catalogs. The agreement between the data and model
increases with increasing luminosity within each redshift bin, and with
increasing redshift overall. Note that the shaded regions in Figures 4 and 5
at $r_{\rm p}<0.1\ h^{-1}\ {\rm Mpc}$ denotes measurements not used for model
fitting.

The clustering of the full magnitude-limited parent galaxy samples and
corresponding mock catalogs is shown for each redshift bin in Figure 6, with
$\omega_{\rm p}(r_{\rm p})$ for each redshift bin offset by 0.15 dex for
clarity. Agreement between the model and data increases with increasing
redshift. We emphasize that while the model $\omega_{\rm p}(r_{\rm p})$
deviates from that of the data for the full magnitude-limited parent samples
in the ${0.4<z_{\rm phot}<0.5}$ redshift bin, there is still good data–model
agreement within the highest luminosity bins, where the vast majority of LRGs
reside, for both the $W1$\- and $z$-band models (Figures 4 and 5,
respectively).

###  3.4. Color assignment

Figure 7.—  Illustration of the model color assignment procedure (§3.4), which
equates halo $z_{\rm starve}$ with galaxy color at fixed luminosity. The
purple dotted curve is the cumulative distribution of ${}^{0.43}(r-z)$ color
for galaxies in the $-23.0>{{}^{0.43}M_{z}}>-23.05$ luminosity bin, while the
gray solid curve is the cumulative distribution of halo $z_{\rm starve}$ for
mock galaxies in the same luminosity bin. The magenta arrows indicate that a
mock galaxy in this luminosity bin in a halo with $z_{\rm starve}\approx 2.6$
is assigned a ${}^{0.43}(r-z)$ color of $0.96$. The faint purple and gray
curves are the galaxy color and halo $z_{\rm starve}$ distributions of
additional luminosity bins.

Figure 7 shows an illustration of our color assignment algorithm (Eq. 1) with
age distribution matching for the ${0.4<z_{\rm phot}<0.5}$ redshift bin of the
$z$-band model. The dotted purple curve is the cumulative distribution of
${}^{0.43}(r-z)$ color for galaxies in the $-23.0>{{}^{0.43}M_{z}}>-23.05$
luminosity bin, and the solid gray curve is the halo $z_{\rm starve}$
distribution of $z_{\rm sim}=0.425$ mock galaxies in the same model luminosity
bin. The magenta arrows in Figure 7 indicate that a mock galaxy in this
luminosity bin in a halo with $z_{\rm starve}\approx 2.6$ is assigned a
${}^{0.43}(r-z)$ color of 0.96.

Implementation of the color assignment algorithm with no scatter in the
color–$z_{\rm starve}$ relationship yields mock LRG samples (see §3.5 below)
that generally overpredict the clustering amplitude compared to the data, with
the exception optical LRGs selected from the $W1$-band model (see Figure 10).
We tested the effect of introducing scatter into the color–$z_{\rm starve}$
relation on the clustering of mock LRGs using an analogous method to how
scatter is incorporated into the luminosity–$v_{\rm peak}$ relation (see
§3.3). After assigning model colors $\mathcal{C}$ without scatter according to
Eq. 1 as described above, we assign to each mock galaxy a new color
$\mathcal{C}^{\prime}$, where $\mathcal{C}^{\prime}$ is drawn from a Gaussian
distribution centered at $\mathcal{C}$ with width $\sigma_{\rm color}$. We
then rank order all mock galaxies by $\mathcal{C}^{\prime}$, rank order the
original distribution of mock colors $\mathcal{C}$, and assign the ordered
original color distribution to the mock galaxy catalog ordered by
$\mathcal{C}^{\prime}$.

For mock LRGs (both optical and IR) from the $z$-band model, as well as mock
IR LRGs from the $W1$-band model, clustering amplitude decreases with
increasing $\sigma_{\rm color}$ until $\sigma_{\rm color}$ is sufficiently
large that model colors are effectively assigned at random, at which point the
decrease in clustering amplitude levels off at a constant value that still
overpredicts $\omega_{\rm p}(r_{\rm p})$ relative to the data. The exception
to this is the predicted clustering of mock optical LRGs from the $W1$-band
model, which largely agrees with the data with _no_ scatter in the
color–$z_{\rm starve}$ relation. Additionally, $\sigma_{\rm color}>0$ in this
instance does _not_ affect the predicted clustering amplitude of mock optical
LRGs. This implies that color is uncorrelated with halo age for LRGs, which we
discuss further in §4.1, and motivates our decision to proceed with two
versions of our models (rather than introduce and constrain $\sigma_{\rm
color}$ as an additional parameter): (1) a “default” age distribution model,
in which galaxy color increases monotonically with increasing halo $z_{\rm
starve}$ at fixed luminosity (i.e., $\sigma_{\rm color}=0$), and (2) a “random
color” model, in which there is no correlation between galaxy color and halo
$z_{\rm starve}$. Both models yield magnitude-limited mock galaxy catalogs
with identical color–magnitude distributions that reproduce the target
distribution of the relevant data; the only difference is whether galaxy color
is maximally correlated (default model) or entirely uncorrelated (random color
model) with halo $z_{\rm starve}$.

###  3.5. Selecting mock LRGs

$z$-band model |  $W1$-band model
---|---
|
|
|
Figure 8.—  LRG fractions in color–magnitude space used to select optical and
IR LRG samples from our magnitude-limited mock galaxy catalogs. The first and
second columns show optical and IR LRG targets, respectively, in ${r-z}$ color
versus $M_{z}$ magnitude space. The third and fourth columns show optical and
IR LRG targets, respectively, in ${r-W1}$ color versus $M_{W1}$ magnitude
space. Each row shows a different redshift bin. Note that the distributions of
optical LRGs in the $z$-band model (first column) and IR LRGs in the $W1$-band
model (fourth column) are compact by design, i.e., the fraction of galaxies
that are LRGs is either $\sim 0$ or $\sim 1$ nearly everywhere in
color–magnitude space. In contrast, the distributions of IR LRGs in the
$z$-band model (second column) and optical LRGs in the $W1$-band model (third
column) have a broader gradient between 0 and 1, especially toward fainter
magnitudes. The boxed region in each panel marks the luminous end of the red
sequence, where essentially all galaxies are expected to be LRGs. The LRG
fraction in this boxed region is at least 98% for the $W1$-band model (third
and fourth columns), but only 92–94% for the $z$-band model (first and second
columns). Figure 9.—  The color-coded points in each panel are bright
(${M_{z}<-23.7}$), red-sequence objects with photometric redshifts between 0.4
and 0.7 that are _not_ DESI LRG targets, i.e., objects with ${M_{z}<-23.7}$ in
the first two columns of Figure 8 that cause the LRG fraction at the bright
end of the red sequence to be $<1$. The dashed and solid black lines in panels
(a), (b), (d), and (e) are the various color–color and color–magnitude cuts
that define the optical and IR DESI LRG target samples (Eqs. 4 and 5,
respectively). The light blue and gray contours show the distributions of
_all_ optical and IR DESI LRG targets with photometric redshifts between 0.4
and 0.7.

DESI LRG target selection (see §2.3 above) is entirely a function of apparent
$g$, $r$, $z$, and $W1$ (and $z_{\rm fiber}$) magnitudes. As such it would be
ideal to have mock galaxy catalogs with each of these model magnitudes for
every galaxy; mock LRG samples could then be obtained simply by applying the
DESI LRG target selection functions (Eqs. 4 and 5) to the full mock for each
redshift bin. However the SHAM and age distribution matching techniques we
employ create mock galaxy catalogs with absolute $z$\- or $W1$-band magnitudes
and corresponding ${r-z}$ or ${r-W1}$ colors.

To select mock DESI LRG targets using only these model quantities we first
identify where actual DESI LRG targets reside in $K$-corrected color–magnitude
space in each redshift bin for both optical and IR DESI LRG targets. In each
redshift bin we compute the fractions of galaxies that are optical and IR DESI
LRG targets in narrow two-dimensional bins in color–magnitude space. This is
shown in Figure 8.

We then use these optical and IR LRG fractions in color–magnitude space to
statistically select mock optical and IR LRGs from our mock galaxy catalogs in
the corresponding model color–magnitude space, i.e., in each color–magnitude
bin we draw mock galaxies and flag them as LRGs until the LRG fraction in that
bin matches the value from the data.

It is worth noting that the LRG fraction in most of color–magnitude space is
either $\sim 0$ or $\sim 1$ by design, especially for optical LRGs in ${r-z}$
versus $M_{z}$ (i.e., optical) space, and IR LRGs in ${r-W1}$ versus $M_{\rm
W1}$ (i.e., IR) space. This is clearly shown in Figure 8.

The first and second columns of Figure 8 show the fractions of optical and IR
LRG targets, respectively, in optical space, i.e., ${r-z}$ color versus
$M_{z}$ magnitude, and the red sequence is clearly visible in each panel. For
example, in the ${0.4<z_{\rm phot}<0.5}$ redshift bin the red sequence
corresponds almost entirely to the region where the LRG fraction is $\sim 1$,
i.e., where ${M_{z}\lesssim-22.5}$ and ${r-z}$ is between $\sim 0.9$ and $\sim
1.3$ (with a shift toward redder colors in higher redshift bins).

Interestingly, the LRG fraction begins to deviate from $\sim 1$ in the _most
luminous_ region of optical space ($M_{z}\lesssim-23.7$), and this is true for
both optical and IR LRG targets. Both selections (Eqs. 4 and 5) exclude some
of the most luminous red-sequence objects in optical space, but _not_ in IR
space (third and fourth columns of Fig. 8).

To investigate why these luminous, red-sequence objects in the first two
columns of Figure 8) are not selected as DESI LRG targets we looked at their
positions in color–color and color–magnitude space relative to the full set of
DESI LRG target selection cuts (Eqs. 4 and 5). This is shown in Figure 9,
where light blue and gray contours in each panel show the distributions of
_all_ optical and IR DESI LRG targets with photometric redshifts between 0.4
and 0.7.

The color-coded points in each panel are luminous (${M_{z}<-23.7}$), red-
sequence objects in the same redshift range that are _not_ DESI LRG targets
despite passing the stellar cut (Eqs. 4a and 5a; solid black line in panel
(a)) and the optical and IR color–magnitude cuts (Eq. 4c; dashed black line in
panel (e), and Eq. 5c; solid black line in panel (b)), as well as the $z_{\rm
fiber}$ cut (Eqs. 4d and 5d; not shown).

Panel (d) of Figure 9 shows that the low-redshift or blue color cut (Eqs. 4b
and 5b) excludes these objects from selection as DESI LRG targets.

##  4\. Predicted LRG properties

In this section we present the results of our $z$-band and $W1$-band models.
We report the HOD parameters of mock LRGs and compare the predicted clustering
to the data and to relevant studies in the literature.

###  4.1. LRG clustering

Figure 10.— Predicted clustering of IR (top row) and optical (bottom row) mock
LRGs compared to the clustering of the relevant LRG target samples from the
data (black points in each panel). Dotted purple lines in each panel show mock
LRGs from the default $z$-band model (monotonic correspondence between $r-z$
color and $z_{\rm starve}$ at fixed luminosity), while solid purple lines show
LRGs from the $z$-band model with randomized colors. Dashed orange lines in
each panel show mock LRGs from the default $W1$-band model, and solid orange
lines show LRGs from the $W1$-band model with randomized colors. Each column
shows a different redshift bin.  Figure 11.— IR color ($r-W1$) versus optical
color ($r-z$) for IR (left panel) and optical (right panel) LRGs at
$0.4<z_{\rm phot}<0.5$.

Figure 12.— Distributions of ${r-z}$, ${g-r}$, and ${r-W1}$ colors of optical
(shaded histograms) and IR (open histograms) LRGs from above (Box A; top three
panels) and along (Box B; bottom three panels) the red sequence in ${r-W1}$
versus $M_{W1}$ color–magnitude space. The regions of color–magnitude space
corresponding to Box A and Box B are shown in Figure 3. The mean, median, and
skewness of each distribution, as well as the number of LRGs in each
population, are given in Table 5. Optical and IR LRGs from along the red
sequence have similar distributions of each color, while optical LRGs above
the red sequence have redder ${r-z}$ and ${g-r}$ and bluer ${r-W1}$ colors
compared to IR LRGs in the same $W1$-band luminosity range.  Table 5Statistics
of distributions of ${r-z}$, ${g-r}$, and ${r-W1}$ colors of optical and IR
LRGs selected from above (Box A) and along (Box B) the red sequence in IR
color–magnitude space (${r-W1}$ versus $M_{W1}$). The boxed regions A and B
are shown for the $0.4<z_{\rm phot}<0.5$ redshift bin in Figure 3, and the
color distributions are shown in Figure 12. The mean, median, and skewness of
each distribution is shown for IR (optical) LRGs in bold (plain) text.
---

Our $z$-band and $W1$-band models yield magnitude-limited mock galaxy catalogs
from which we select mock LRGs. From these results we predict the clustering
of DESI LRG target samples.

Figure 10 shows the predicted clustering of IR and optical mock LRGs compared
to the relevant DESI LRG target sample. In addition to the data (black points
with error bars), each panel of Figure 10 shows the clustering of mock LRGs
according to the $z$-band and $W1$-band models, with and without the addition
of scatter to the model color–$z_{\rm starve}$ relation (§3.4).

As shown in the top row of Figure 10, the $z$-band and $W1$-band models
predict very similar clustering signals for IR LRGs. The predicted amplitude
of both the one-halo and two-halo terms exceeds that of the data in the two
lowest redshift bins, and is in better agreement in the highest redshift bin
($z_{\rm sim}=0.628$). Across all redshift bins the addition of scatter to the
color–$z_{\rm starve}$ relation brings the predicted clustering amplitude of
IR LRGs closer to the data, reducing the clustering amplitude of the one-halo
term by $\sim 15$–$20\%$, and the two-halo term by $\sim 10\%$, depending on
redshift bin. However, this reduction in the discrepancy between the model and
data is a ceiling achieved only by adding such a high degree of scatter that
color is uncorrelated with our proxy for halo assembly history777We also
tested using proxies for IR and optical LRGs by selecting the $N$ brightest
mock galaxies from the relevant magnitude-limited mock galaxy catalog to
recreate the number densities of the IR and optical DESI LRG target samples,
instead of applying the DESI LRG target selections to the magnitude-limited
mocks as described in §3.5. The result was a slightly larger predicted
clustering amplitude compared to the “random colors” models that do use the
DESI LRG target selections, shown in Figure 10. ($z_{\rm starve}$).

The bottom row of Figure 10 shows that the $z$-band and $W1$-band models
predict different clustering amplitudes for optical LRGs. The $z$-band model
overpredicts the clustering amplitude of optical LRGs to a similar degree as
the overprediction for IR LRGs. Assigning colors at random such that galaxy
color is uncorrelated with halo assembly history again reduces this
discrepancy of the one-halo (two-halo) term by $\sim 20$–$25\%$ ($\sim
10$–$15\%$).

On the other hand, the $W1$-band model prediction for the clustering of
optical LRGs is in very good agreement with the data, especially the two-halo
term (the one-halo term is actually slightly underpredicted). Additionally,
adding _any_ amount of scatter to the color–$z_{\rm starve}$ relation has no
effect on the predicted clustering amplitude in this case.

The distributions of optical and IR LRGs in color–magnitude space are helpful
for interpreting Figure 10. As Figure 3 shows, in optical space
(${M_{r}-M_{z}}$ versus $M_{z}$), both optical and IR LRGs are largely
confined to the red sequence, with the distribution extending below it toward
bluer color, especially at higher redshift.

In contrast, in IR space (${M_{r}-M_{W1}}$ versus $M_{W1}$) the galaxy
distribution does extend significantly above the red sequence. For example, in
the ${0.4<z_{\rm phot}<0.5}$ redshift bin, the red sequence is at
${M_{r}-M_{W1}\sim 1.6}$, but the distribution extends to
${M_{r}-M_{W1}\gtrsim 3.0}$. IR LRGs occupy the entire region of excess
${M_{r}-M_{W1}}$ color across the full range of $M_{W1}$, while optical LRGs
occupy only the part of this region that corresponds to more luminous
$M_{W1}$, i.e., the optical LRG selection excludes a population of galaxies
with very red ${M_{r}-M_{W1}}$ colors and moderately luminous $W1$-band
luminosities that are included by the IR selection.

The top panels of Figure 12 explore the differences in the color distributions
of optical and IR LRGs located above the red sequence (Box A in Figure 3) in
${r-W1}$ versus $M_{W1}$ color–magnitude space. The bottom three panels (Box B
in Figure 3) show the color distributions of optical and IR LRGs from the same
$M_{W1}$ range _along_ the red sequence.

Along the infrared red sequence, the SEDs of optical and IR LRGs are quite
similar. However, the two selections diverge above the red sequence. Table 5
lists the mean, median, and skewness of each color distribution, as well as
the number of LRGs in each population.

Above the red sequence IR LRGs have bluer ${r-z}$ and ${g-r}$ colors than
optical LRGs in the same region of color–magnitude space, while IR LRGs have
redder ${r-W1}$ colors than their optical counterparts from the same
color–magnitude region. The same general trend is seen along the red sequence,
but the color differences between the optical and IR LRG selections are
smaller, i.e., along the red sequence IR LRGs also have bluer ${r-z}$ and
${g-r}$ colors than optical LRGs, but by only half as much as above the red
sequence. Similarly, IR LRGs along the red sequence have redder ${r-W1}$
colors than their optical counterparts, but the color difference is two to
three times smaller than for IR and optical LRGs above the red sequence.

Activity from active galactic nuclei (AGN) may artificially inflate the
observed IR ($W1$-band) luminosities and artificially redden the $r-W1$ colors
of some galaxies that pass the IR LRG selection (e.g., Webster et al. 1995;
Georgakakis et al. 2009; Banerji et al. 2012; Glikman et al. 2012; Kim & Im
2018; Klindt et al. 2019; Rivera et al. 2021), relative to galaxies with
comparable IR luminosities and ${r-W1}$ colors but no AGN activity. The bluer
${r-z}$ colors of IR LRG targets cause these objects to be excluded by the
optical LRG selection (e.g., panel (e) of Figure 9). Our IR ($W1$-band) model
would assign mock IR LRGs with artificially red $r-W1$ colors and high IR
luminosities due to AGN activity to halos with higher bias than if their
$r-W1$ colors and $W1$-band magnitudes were not influenced by AGN activity.
This could explain why our $W1$-band model overpredicts the clustering
amplitude of IR LRGs relative to optical LRGs, as shown in Figure 10.

The influence of AGN activity may also explain the lack of correlation between
IR and optical color when assigning colors to mock galaxies based solely on a
proxy for halo age (here $z_{\rm starve}$). This modeling assumption
attributes galaxy color entirely to halo mass accretion history, and does not
account for how baryonic effects such as AGN activity may contribute to galaxy
color. In other words, optical color is not necessarily a reliable proxy for
IR color. This is also illustrated by Figure 11, which shows optical (${r-z}$)
versus IR (${r-W1}$) color for both IR and optical LRGs in the $0.4<z_{\rm
phot}<0.5$ redshift bin. The optical and IR colors of optical LRGs are more
closely correlated than for IR LRGs, but in both cases there is considerable
scatter in $r-W1$ at fixed $r-z$. Xu et al. (2022) also find a lack of
correlation between halo assembly properties and galaxy color using a semi-
analytic model.

###  4.2. The LRG–halo connection

Figure 13.— Predicted HODs of IR (top row) and optical (bottom row) mock LRGs from both the $z$-band (heavier purple lines) and $W1$-band (heavier orange lines) models. Also shown in fainter purple and orange lines are HODs of the $z$-band and $W1$-band magnitude-limited mock galaxy catalogs from which mock LRGs are selected. Solid (dotted) lines show results for central (satellite) galaxies and LRGs. Each column shows a different redshift bin. The best-fit values of a standard five-parameter HOD model, as well as the satellite fraction, for each mock LRG sample are given in Table 8.  Table 6Predicted halo occupation completeness for mock central LRGs by selection (IR or optical). Results for the $W1$-band ($z$-band) model are in bold (plain) text.  $z_{\rm sim}$ | LRG selection | peak $\langle N_{\rm cen}|M_{\rm vir}\rangle$aaHighest central occupation fraction reached across all halo $M_{\rm vir}$. | min $M_{\rm vir}$ of peak $\langle N_{\rm cen}|M_{\rm vir}\rangle$bbMinimum halo $M_{\rm vir}$ at which the highest central occupation fraction is achieved. |  $\langle N_{\rm cen}|(M_{\rm vir}>10^{14.75}\ h^{-1}{\rm M}_{\odot})\rangle$ccMean central occupation fraction at the largest halo masses ($M_{\rm vir}\gtrsim 10^{14.75}\ h^{-1}{\rm M}_{\odot}$), which is universally less or equal to than the peak central occupation fraction.
---|---|---|---|---
$0.425$ | IR |  0.92 0.91 |  14.45 14.85 |  0.74 0.70
optical |  0.85 0.90 |  14.35 14.55 |  0.62 0.71
$0.523$ | IR |  0.91 0.87 |  14.45 14.65 |  0.77 0.61
optical |  0.85 0.84 |  14.75 14.65 |  0.69 0.64
$0.628$ | IR |  0.96 0.82 |  14.85 14.55 |  0.96 0.31
optical |  0.92 0.84 |  14.95 14.55 |  0.92 0.38

We compute the mean central and satellite halo occupation statistics as a
function of halo mass directly from the mock galaxy catalogs and halo
abundances from the MDPL2 simulation. The results are shown in Figure 13. Each
panel of Figure 13 shows ${\langle N|M_{\rm vir}\rangle}$ for the full
magnitude-limited mock galaxy catalog in gray, while optical and IR LRGs are
shown in purple and orange, respectively.

The key conclusions of Figure 13 are summarized in Table 6, which gives the
peak value of ${\langle N_{\rm cen}\rangle}$ for optical and IR LRGs predicted
by both the $z$-band and $W1$-band models, as well as the mean value of
${\langle N_{\rm cen}\rangle}$ for the most massive halos, which is
universally less than the corresponding peak value.

The $z$-band model does not predict significant differences between the
populations selected by the IR and optical LRG selection functions. According
to this model the central LRG halo occupation fraction peaks at $\sim 90\%$ at
$z\sim 0.43$, and drops to $\sim 70\%$ at the largest halo masses, while by
$z\sim 0.63$ $\langle N_{\rm cen}\rangle$ peaks at 82–84% and falls to $\sim
35\%$ for the most massive halos, regardless of LRG selection function.

In contrast, the $W1$-band model _does_ predict a different LRG–halo
relationship for IR versus optical LRGs. In all redshift bins this model has
${\langle N_{\rm cen}\rangle}$ for IR LRGs peaking at around 92% and remaining
at $\gtrsim 75\%$ for the most massive halos, while ${\langle N_{\rm
cen}\rangle}$ for optical LRGs peaks at 85–92%, and falls at the largest halo
masses as low as 62% ($z\sim 0.43$) to 69% ($z\sim 0.52$), although it remains
at 92% for the highest redshift bin ($z\sim 0.63$).

Table 7Prior ranges for HOD fit parameters (Eqs. 15 and 16).  Parameter | Prior interval
---|---
$\log(M_{\rm min})$ | $(11.0,\ 14.0)$
$\sigma_{\log M}$ | $(0.001,\ 1.5)$
$\alpha$ | $(0.0,\ 2.0)$
$\log(M_{0})$ | $(11.0,\ 14.0)$
$\log(M_{1})$ | $(11.5,\ 15.5)$
Table 8Best-fit values of standard five-parameter HOD model (Eqs. 15 and 16), and the satellite fraction, $f_{\rm sat}$, of IR and optical mock LRG samples. Results for the $W1$-band ($z$-band) model are in bold (plain) text.  $z_{\rm sim}$ | LRG selection | $\log(M_{\rm min}/h^{-1}{\rm M}_{\odot})$ | $\sigma_{\log M}$ | $\alpha$ | $\log(M_{0}/h^{-1}{\rm M}_{\odot})$ | $\log(M_{1}/h^{-1}{\rm M}_{\odot})$ | $f_{\rm sat}$
---|---|---|---|---|---|---|---
0.425 | IR |  13.28 13.35 |  0.93 0.95 |  0.98 1.00 |  12.93 12.87 |  14.02 14.07 |  14.8 14.5
optical |  13.32 13.26 |  1.05 0.93 |  0.95 0.98 |  12.91 12.89 |  14.00 14.00 |  14.9 14.8
0.523 | IR |  13.37 13.52 |  0.99 1.14 |  0.94 0.91 |  12.94 12.98 |  14.05 14.09 |  14.4 14.5
optical |  13.45 13.44 |  1.15 1.11 |  0.95 0.96 |  12.84 12.85 |  14.05 14.05 |  14.7 14.6
0.628 | IR |  13.26 13.40 |  1.02 1.18 |  0.95 0.95 |  12.88 12.82 |  13.95 13.99 |  14.1 14.4
optical |  13.30 13.31 |  1.19 1.16 |  0.96 0.94 |  12.78 12.82 |  13.93 13.93 |  14.4 14.7

We fit a standard five-parameter HOD model to each optical and IR mock LRG
sample Table 8. For ease of comparison with their results we use the same
functional form as Zhou et al. (2021), who derive HOD parameters from
clustering measurements of DESI LRG targets using photometric redshifts in
bins that approximately match the redshift bins we use. This parameterization
is described in detail in Zheng et al. (2005) and Zheng et al. (2007).
Briefly, the probability $\langle N_{\rm cen}\rangle$ that a halo of virial
mass $M_{\rm vir}$ hosts a central galaxy is given by

$\langle N_{\rm cen}|M_{\rm vir}\rangle=\frac{1}{2}\\!\left(1+{\rm
erf}{\left[\frac{\log(M_{\rm vir})-\log(M_{\rm min})}{\sigma_{\log
M}}\right]}\right),$ (15)

where the parameter $\log(M_{\rm min})$ is the minimum halo mass for hosting a
central galaxy, and the parameter $\sigma_{\log M}$ defines the steepness of
the transition of $\langle N_{\rm cen}\rangle$ from $\sim 0$ at low $M_{\rm
vir}$ to $\sim 1$ at high $M_{\rm vir}$.

The mean number of satellite galaxies, $\langle N_{\rm sat}\rangle$, hosted by
a halo of virial mass $M_{\rm vir}$ is approximated by a power law given by

$\langle N_{\rm sat}|M_{\rm vir}\rangle={\left(\frac{M_{\rm
vir}-M_{0}}{M_{1}}\right)}^{\alpha},$ (16)

where $M_{0}$, $M_{1}$, and $\alpha$ are the remaining free parameters of the
HOD fit. The HOD fits of Zhou et al. (2021) incorporate a sixth nuisance
parameter to account for photometric redshift uncertainty, which in this work
is addressed by the magnitude-dependent line-of-sight scatter parameter,
$\sigma_{\rm LOS}$ (Eq. 14). Table 7 gives the prior interval for each HOD
parameter. The best-fit HOD parameters and satellite fraction of each mock LRG
sample are given in Table 8.

###  4.3. Comparison with other “DESI-like” LRG studies

Figure 14.— Comparison of LRG HOD parameters from this work at $z\sim 0.63$
(orange lines) with the analytic HOD fits of Zhou et al. (2021) (green lines)
and semi-analytic model of Hernández-Aguayo et al. (2021) (gray lines) in
comparable redshift bins. Solid, dotted, and dash-dotted lines show results
respectively for central LRGs, satellite LRGs, and the combination of both.

There are two main studies of the LRG–halo connection of “DESI-like” LRGs
suitable for comparison with our results. Zhou et al. (2021) (who provide the
photometric redshifts used here) measure the clustering of LRGs at
${0.41<z_{\rm phot}<0.93}$ selected from DECaLS DR7 photometry using the
optical target selection (Eq. 4), and fit a standard five-parameter analytic
HOD model in five redshift bins.

Zhou et al. (2021) find little to no evolution of HOD parameters for LRGs
across ${0.4\lesssim z\lesssim 0.8}$. Our results are consistent with a lack
of HOD parameter evolution across the subset of this redshift range that we
model, although some of our predicted parameter values themselves differ from
those of Zhou et al. (2021), particularly $\sigma_{\log M}$. Zhou et al.
(2021) find a very steep transition from a ${\langle N{\rm cen}\rangle\sim 0}$
to 1 ($\sigma_{\log M}\sim 0$–0.28; see their Figure 13), while our models
predict a more gradual transition, with ${\sigma_{\log M}\sim 1}$.

Figure 14 compares our predicted LRG halo occupation statistics at $z\sim
0.63$ with the best-fit HOD parameterization of Zhou et al. (2021) in the
relevant redshift bin from their study ($0.61<z_{\rm phot}<0.72$). Of the
redshift bins Zhou et al. (2021) use that overlap with this work, this is the
bin in which they find the greatest deviation of $\langle N{\rm cen}\rangle$
from a step function, although we also select this bin for Figure 14 to enable
comparison with an additional study, discussed later in this section.

Our HOD predictions for satellite LRGs agree with those of Zhou et al. (2021)
at halo masses below ${\log(M_{\rm vir}/h^{-1}{\rm M}_{\odot})\sim 14}$, and
our overall satellite fractions of $f_{\rm sat}\sim 0.14$–0.15 across both
model bands and LRG selections is consistent with Zhou et al. (2021), who find
$f_{\rm sat}\sim 0.13$–0.16, depending on redshift.

Zhou et al. (2021) obtain best-fit parameters for analytic forms of central
and satellite LRG HODs (Eqs. 15 and 16), which will necessarily correspond to
100% of massive halos containing a central LRG due to the functional form of
Eq. 15 (this is also the case for analytic HOD fits to our mock LRG samples,
shown in Table 8). Enforcing 100% halo occupation at the high-mass end by LRGs
also suppresses accounting for potential contamination by lower-mass halos.

As Figure 13 shows, our models do _not_ predict that central LRG halo
occupation reaches unity, instead peaking at 82–96% and falling in most cases
to $\sim 60$–77% for the most massive halos (${\log(M_{\rm vir}/h^{-1}{\rm
M}_{\odot})\gtrsim 14,75}$), although this depends on both model band ($z$ or
$W1$) and LRG selection (IR or optical; see Table 6).

The other study to which we can compare our results is Hernández-Aguayo et al.
(2021), who use the galform semi-analytic galaxy formation model (Cole et al.
2000) and nine snapshots from the Planck-Millennium $N$-body simulation (Baugh
et al. 2018) to study the galaxy–halo connection of “DESI-like” LRGs at
redshifts between 0.6 and 1.0. Their LRG HOD results for $z=0.64$ are overlaid
with our $z_{\rm sim}\sim 0.63$ results in Figure 14.

Hernández-Aguayo et al. (2021) predict central and satellite LRG abundances
below our model predictions, although they note that galform underpredicts the
abundance of LRGs across the entire redshift range they study, with the
greatest discrepancy at $z\sim 0.6$–0.7. The shape of their central LRG HOD is
notably similar to ours, consistent with a gradual transition to a peak
central LRG halo occupation less than unity that decreases at the largest halo
masses. Specifically, Hernández-Aguayo et al. (2021) find that ${\langle
N_{\rm cen}|M_{\rm vir}\rangle}$ for “DESI-like” LRGs reaches a maximum of
$\sim 70\%$ at ${\log(M_{\rm vir}/h^{-1}{\rm M}_{\odot})\sim 13.75}$, and is
as low as $\sim 40\%$ at ${\log(M_{\rm vir}/h^{-1}{\rm M}_{\odot})\sim 14.5}$.

##  5\. Summary and conclusion

We have used a SHAM and age distribution matching framework to construct
magnitude-limited mock galaxy catalogs at $z\sim 0.43$, 0.52, and 0.63 with
${r-z}$ and ${r-W1}$ colors. From these catalogs we select mock LRG samples
according to both the optical and IR DESI LRG target selection functions. This
work is the first application SHAM modeling in the infrared, complimenting the
few existing studies of the galaxy–halo connection of “DESI-like” LRGs.

Our models reproduce the number densities, luminosity functions, color
distributions, and magnitude-dependent projected clustering (in the luminosity
range of DESI LRG targets) of the parent galaxy samples from the Legacy
Surveys DR9 photometry that serve as the basis for DESI LRG target selection.
With the mock LRG samples selected from our magnitude-limited mock galaxy
catalogs, we predict the halo occupation statistics of both optical and IR
DESI LRGs at a fixed cosmology. We assess the differences between these two
LRG populations, as well as the effect of using the $z$-band versus the 3.4
micron $W1$-band for SHAM and age distribution matching.

The main results of this work are:

  1. 1.

Both the optical and IR DESI LRG target selections exclude some of the most
luminous galaxies that would appear to be LRGs based on their position on the
red sequence in optical color–magnitude space (Figure 8). This is a result of
the specific DESI LRG target selection cut intended to exclude blue, low-
redshift (${z\lesssim 0.4}$) galaxies (Figure 9, §3.5).

  2. 2.

Optical and IR LRGs occupy similar regions of optical color–magnitude space
(${r-z}$ versus $M_{z}$), where the red sequence corresponds to the very
reddest ${r-z}$ colors (Figure 3). In IR color–magnitude space (${r-W1}$
versus $M_{W1}$) there is a sizable galaxy population at significantly redder
colors than the red sequence. IR LRGs occupy most of this region of excess
${r-W1}$ color, but optical LRGs are largely excluded. This could be a result
of AGN activity artificially inflating the $W1$-band luminosities for some of
these objects (§4.1).

  3. 3.

There are clear distinctions between the LRG samples obtained from the optical
versus IR selections that are apparent from the data alone, namely among their
optical ($r-z$ versus $M_{z}$) and IR ($r-W1$ versus $M_{W1}$) color–magnitude
distributions, clustering, and color–color distributions. Our IR-based
($W1$-band) model predicts greater differences between the halo occupation
statistics of optical and IR LRGs than the $z$-band model (Figure 13, Table
6), and is therefore the preferred regime for this comparative study of the
two selections.

  4. 4.

Age distribution matching, which assumes a monotonic correlation between halo
age and galaxy color at fixed luminosity, tends to overpredict the clustering
amplitude of DESI LRGs (Figure 10). Introducing scatter into the assumed
age–color relation improves agreement between predictions and the data
somewhat, although the greatest improvement comes from increasing this scatter
to a level that is equivalent to assigning colors at random. Our models
therefore suggest that either galaxy color is uncorrelated with halo age in
the LRG regime, or there is some additional model parameter that has been
neglected and is possibly related to our use of photometric redshifts.

  5. 5.

Both DESI LRG target selections yield populations with a non-trivial LRG–halo
connection that does not reach unity ($\langle N_{\rm cen}\rangle\sim 1$) for
the most massive halos (Figure 13, §4.2). However, the IR selection achieves
greater completeness ($\langle N_{\rm cen}\rangle\geq 91\%$) than the optical
selection ($\langle N_{\rm cen}\rangle\sim 85$–$92\%$) across all redshift
bins studied. Our results of $\langle N_{\rm cen}\rangle<1$ is qualitatively
consistent with the SAM predictions of Hernández-Aguayo et al. (2021) (Figure
14, §4.3), although they find a lower maximum completeness of $\sim 70\%$ at
$z\sim 0.64$.

A natural extension of this work would be to utilize spectroscopic redshifts
from the DESI BGS sample (Hahn et al. 2022) to conduct a similar study at
lower redshift. The BGS888Specifically the BGS Bright sample; BGS also
contains BGS Faint and AGN samples with different selection criteria. is
nearing completion and is obtaining spectroscopic redshifts for a magnitude-
limited ($r<19.5$) sample of over 800 galaxies per $\mathrm{deg}^{2}$ at
$z<0.4$. While the LRG target selection algorithms used in this work are
designed to identify LRGs at $0.4\lesssim z\lesssim 1.0$, these cuts could
easily be modified to select LRGs in the same redshift range as the BGS. Mock
LRGs could then be selected from SHAM-based magnitude-limited mock galaxy
catalogs tuned to recreate the magnitude and color distributions of the BGS
sample.

## Acknowledgements

The authors thank Zheng Zheng for comments on an earlier draft of this work.

The Legacy Surveys consist of three individual and complementary projects: the
Dark Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David
Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID
#2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall $z$-band Legacy
Survey (MzLS; Prop. ID #2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS
together include data obtained, respectively, at the Blanco telescope, Cerro
Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward
Observatory, University of Arizona; and the Mayall telescope, Kitt Peak
National Observatory, NOIRLab. The Legacy Surveys project is honored to be
permitted to conduct astronomical research on Iolkam Du’ag (Kitt Peak), a
mountain with particular significance to the Tohono O’odham Nation.

NOIRLab is operated by the Association of Universities for Research in
Astronomy (AURA) under a cooperative agreement with the National Science
Foundation.

This project used data obtained with the Dark Energy Camera (DECam), which was
constructed by the Dark Energy Survey (DES) collaboration. Funding for the DES
Projects has been provided by the U.S. Department of Energy, the U.S. National
Science Foundation, the Ministry of Science and Education of Spain, the
Science and Technology Facilities Council of the United Kingdom, the Higher
Education Funding Council for England, the National Center for Supercomputing
Applications at the University of Illinois at Urbana-Champaign, the Kavli
Institute of Cosmological Physics at the University of Chicago, Center for
Cosmology and Astro-Particle Physics at the Ohio State University, the
Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M
University, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho
de Amparo, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de
Amparo a Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de
Desenvolvimento Cientifico e Tecnologico and the Ministerio da Ciencia,
Tecnologia e Inovacao, the Deutsche Forschungsgemeinschaft and the
Collaborating Institutions in the Dark Energy Survey. The Collaborating
Institutions are Argonne National Laboratory, the University of California at
Santa Cruz, the University of Cambridge, Centro de Investigaciones
Energeticas, Medioambientales y Tecnologicas-Madrid, the University of
Chicago, University College London, the DES-Brazil Consortium, the University
of Edinburgh, the Eidgenossische Technische Hochschule (ETH) Zurich, Fermi
National Accelerator Laboratory, the University of Illinois at Urbana-
Champaign, the Institut de Ciencies de l’Espai (IEEC/CSIC), the Institut de
Fisica d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig
Maximilians Universitat Munchen and the associated Excellence Cluster
Universe, the University of Michigan, NSF’s NOIRLab, the University of
Nottingham, the Ohio State University, the University of Pennsylvania, the
University of Portsmouth, SLAC National Accelerator Laboratory, Stanford
University, the University of Sussex, and Texas A&M University.

BASS is a key project of the Telescope Access Program (TAP), which has been
funded by the National Astronomical Observatories of China, the Chinese
Academy of Sciences (the Strategic Priority Research Program “The Emergence of
Cosmological Structures” Grant #XDB09000000), and the Special Fund for
Astronomy from the Ministry of Finance. The BASS is also supported by the
External Cooperation Program of Chinese Academy of Sciences (Grant
#114A11KYSB20160057), and Chinese National Natural Science Foundation (Grant
#11433005).

The Legacy Survey team makes use of data products from the Near-Earth Object
Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet
Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by
the National Aeronautics and Space Administration.

The Legacy Surveys imaging of the DESI footprint is supported by the Director,
Office of Science, Office of High Energy Physics of the U.S. Department of
Energy under Contract No. DE-AC02-05CH1123, by the National Energy Research
Scientific Computing Center, a DOE Office of Science User Facility under the
same contract; and by the U.S. National Science Foundation, Division of
Astronomical Sciences under Contract No. AST-0950945 to NOAO.

The Photometric Redshifts for the Legacy Surveys (PRLS) catalog used in this
paper was produced thanks to funding from the U.S. Department of Energy Office
of Science, Office of High Energy Physics via grant DE-SC0007914.

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V.
(www.gauss-centre.eu) and the Partnership for Advanced Supercomputing in
Europe (PRACE, www.prace-ri.eu) for funding the MultiDark simulation project
by providing computing time on the GCS Supercomputer SuperMUC at Leibniz
Supercomputing Centre (LRZ, www.lrz.de).

Some of the results in this paper have been derived using the healpy and
HEALPix package.

## References

  * Banerji et al. (2012) Banerji, M., et al. 2012, MNRAS, 427, 2275
  * Baugh et al. (2018) Baugh, C. M., et al. 2018, MNRAS, 483, 4922
  * Behroozi et al. (2010) Behroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, ApJ, 717, 379
  * Behroozi et al. (2012a) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2012a, ApJ, 762, L31
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57
  * Behroozi et al. (2012b) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2012b, ApJ, 762, 109
  * Behroozi et al. (2012c) Behroozi, P. S., et al. 2012c, ApJ, 763, 18
  * Berlind & Weinberg (2002) Berlind, A. A., & Weinberg, D. H. 2002, ApJ, 575, 587
  * Berlind et al. (2003) Berlind, A. A., et al. 2003, ApJ, 593, 1
  * Blanton & Roweis (2007) Blanton, M. R., & Roweis, S. 2007, AJ, 133, 734
  * Blum et al. (2016) Blum, R. D., et al. 2016, in AAS Meeting Abstracts # 228, Vol. 228, 317–01
  * Bullock et al. (2002) Bullock, J. S., Wechsler, R. H., & Somerville, R. S. 2002, MNRAS, 329, 246
  * Bundy et al. (2015) Bundy, K., et al. 2015, ApJS, 221, 15
  * Cole et al. (2000) Cole, S., et al. 2000, MNRAS, 319, 168
  * Conroy et al. (2006) Conroy, C., Wechsler, R. H., & Kravtsov, A. V. 2006, ApJ, 647, 201
  * Dahlen et al. (2013) Dahlen, T., et al. 2013, ApJ, 775, 93
  * Dawson et al. (2013) Dawson, K. S., et al. 2013, AJ, 145, 10
  * Dawson et al. (2016) Dawson, K. S., et al. 2016, AJ, 151, 44
  * DESI Collaboration et al. (2016a) DESI Collaboration et al. 2016a
  * DESI Collaboration et al. (2016b) DESI Collaboration et al. 2016b
  * Dey et al. (2019) Dey, A., et al. 2019, ApJ, 157, 168
  * Flaugher et al. (2015) Flaugher, B., et al. 2015, AJ, 150, 150
  * Gaia Collaboration (2018) Gaia Collaboration. 2018, VizieR Online Data Catalog, I
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al. 2016, A&A, 595, A1
  * Georgakakis et al. (2009) Georgakakis, A., et al. 2009, MNRAS, 394, 533
  * Glikman et al. (2012) Glikman, E., et al. 2012, ApJ, 757, 51
  * Górski et al. (2005) Górski, K. M., et al. 2005, ApJ, 622, 759
  * Hahn et al. (2022) Hahn, C., et al. 2022, DESI Bright Galaxy Survey: Final Target Selection, Design, and Validation
  * Hearin & Watson (2013) Hearin, A. P., & Watson, D. F. 2013, MNRAS, 435, 1313
  * Hearin et al. (2014) Hearin, A. P., et al. 2014, MNRAS, 444, 729
  * Hernández-Aguayo et al. (2021) Hernández-Aguayo, C., et al. 2021, MNRAS, 503, 2318
  * Høg et al. (2000) Høg, E., et al. 2000, A&A, 355, L27
  * Kim & Im (2018) Kim, D., & Im, M. 2018, A&A, 610, A31
  * Klindt et al. (2019) Klindt, L., et al. 2019, MNRAS, 488, 3109
  * Klypin et al. (2016) Klypin, A., et al. 2016, MNRAS, 457, 4340
  * Landy & Szalay (1993) Landy, S. D., & Szalay, A. S. 1993, ApJ, 412, 64
  * Loveday et al. (2012) Loveday, J., et al. 2012, MNRAS, 420, 1239
  * Moster et al. (2013) Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121
  * Moustakas et al. (2021) Moustakas, J., et al. 2021, Bulletin of the AAS, 53, https://baas.aas.org/pub/2021n1i527p04
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  * Pedregosa et al. (2011) Pedregosa, F., et al. 2011, the Journal of Machine Learning Research, 12, 2825
  * Planck Collaboration et al. (2016) Planck Collaboration et al. 2016, A&A, 594, A13
  * Raichoor et al. (2022) Raichoor, A., et al. 2022, arXiv preprint arXiv:2208.08513
  * Reddick et al. (2013) Reddick, R. M., et al. 2013, ApJ, 771, 30
  * Reid et al. (2015) Reid, B., et al. 2015, MNRAS, 455, 1553
  * Riebe et al. (2013) Riebe, K., et al. 2013, Astronomische Nachrichten, 334, 691
  * Rivera et al. (2021) Rivera, C. G., et al. 2021, A&A, 649, A102
  * Safonova et al. (2021) Safonova, S., Norberg, P., & Cole, S. 2021, MNRAS, 505, 325
  * Saito et al. (2016) Saito, S., et al. 2016, MNRAS, 460, 1457
  * Seljak (2000) Seljak, U. 2000, MNRAS, 318, 203
  * Silva et al. (2016) Silva, D. R., et al. 2016, in AAS Meeting Abstracts # 228, Vol. 228, 317–02
  * Sinha & Garrison (2017) Sinha, M., & Garrison, L. 2017, Corrfunc: Blazing fast correlation functions on the CPU, Astrophysics Source Code Library, record ascl:1703.003
  * Sinha & Garrison (2019) Sinha, M., & Garrison, L. 2019, in Software Challenges to Exascale Computing, ed. A. Majumdar & R. Arora (Singapore: Springer Singapore), 3–20
  * Vale & Ostriker (2004) Vale, A., & Ostriker, J. P. 2004, MNRAS, 353, 189
  * Watson & Conroy (2013) Watson, D. F., & Conroy, C. 2013, ApJ, 772, 139
  * Webster et al. (1995) Webster, R. L., et al. 1995, Nature, 375, 469
  * Wechsler & Tinker (2018) Wechsler, R. H., & Tinker, J. L. 2018, ARA&A, 56, 435
  * Wechsler et al. (2002) Wechsler, R. H., et al. 2002, ApJ, 568, 52
  * Wetzel et al. (2014) Wetzel, A. R., et al. 2014, MNRAS, 439, 2687
  * Wright et al. (2010) Wright, E. L., et al. 2010, AJ, 140, 1868
  * Xu et al. (2022) Xu, X., Zheng, Z., & Guo, Q. 2022, MNRAS, 516, 4276
  * Yang et al. (2012) Yang, X., et al. 2012, ApJ, 752, 41
  * Yang et al. (2013) Yang, X., et al. 2013, ApJ, 770, 115
  * York et al. (2000) York, D. G., et al. 2000, ApJ, 120, 1579
  * Yu et al. (2022) Yu, J., et al. 2022, MNRAS, 516, 57
  * Zentner et al. (2014) Zentner, A. R., Hearin, A. P., & van den Bosch, F. C. 2014, MNRAS, 443, 3044
  * Zheng et al. (2007) Zheng, Z., Coil, A. L., & Zehavi, I. 2007, ApJ, 667, 760
  * Zheng & Weinberg (2007) Zheng, Z., & Weinberg, D. H. 2007, ApJ, 659, 1
  * Zheng et al. (2005) Zheng, Z., et al. 2005, ApJ, 633, 791
  * Zhou et al. (2020) Zhou, R., et al. 2020, Research Notes of the AAS, 4, 181
  * Zhou et al. (2021) Zhou, R., et al. 2021, MNRAS, 501, 3309
  * Zhou et al. (2023) Zhou, R., et al. 2023, AJ, 165, 58
  * Zonca et al. (2019) Zonca, A., et al. 2019, Journal of Open Source Software, 4, 1298
  * Zou et al. (2017) Zou, H., et al. 2017, PASP, 129, 064101