# Generating mock galaxy catalogues for flux-limited samples like the DESI
Bright Galaxy Survey

A. Smith,1 C. Grove,1 S. Cole,1 P. Norberg,1,2 P. Zarrouk,3 S. Yuan,4 J.
Aguilar,5 S. Ahlen,6 D. Brooks,7 T. Claybaugh,5 A. de la Macorra,8 P. Doel,7
J. E. Forero-Romero,9,10 E. Gaztañaga,11,12,13 S. Gontcho A Gontcho,5 C.
Hahn,14 R. Kehoe,15 A. Kremin,5 M. E. Levi,5 M. Manera,16,17 A. Meisner,18 R.
Miquel,19,17 J. Moustakas,20 J. Nie,21 W. J. Percival,22,23,24 M. Rezaie,25 G.
Rossi,26 E. Sanchez,27 H. Seo,28 G. Tarlé,29 and Z. Zhou21
1 Institute for Computational Cosmology, Department of Physics, Durham
University, South Road, Durham DH1 3LE, UK
2 Centre for Extragalactic Astronomy, Department of Physics, Durham
University, South Road, Durham, DH1 3LE, UK
3 Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE), CNRS/IN2P3
& Sorbonne Université, 4 place Jussieu, 75005 Paris, France
4 SLAC National Accelerator Laboratory, Menlo Park, CA 94305, USA
5 Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720,
USA
6 Physics Dept., Boston University, 590 Commonwealth Avenue, Boston, MA 02215,
USA
7 Department of Physics & Astronomy, University College London, Gower Street,
London, WC1E 6BT, UK
8 Instituto de Física, Universidad Nacional Autónoma de México, Cd. de México
C.P. 04510, México
9 Departamento de Física, Universidad de los Andes, Cra. 1 No. 18A-10,
Edificio Ip, CP 111711, Bogotá, Colombia
10 Observatorio Astronómico, Universidad de los Andes, Cra. 1 No. 18A-10,
Edificio H, CP 111711 Bogotá, Colombia
11 Institut d’Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain
12 Institute of Cosmology & Gravitation, University of Portsmouth, Dennis
Sciama Building, Portsmouth, PO1 3FX, UK
13 Institute of Space Sciences, ICE-CSIC, Campus UAB, Carrer de Can Magrans
s/n, 08913 Bellaterra, Barcelona, Spain
14 Department of Astrophysical Sciences, Princeton University, Princeton NJ
08544, USA
15 Department of Physics, Southern Methodist University, 3215 Daniel Avenue,
Dallas, TX 75275, USA
16 Departament de Física, Serra Húnter, Universitat Autònoma de Barcelona,
08193 Bellaterra (Barcelona), Spain
17 Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of
Science and Technology, Campus UAB, 08193 Bellaterra Barcelona, Spain
18 NSF’s NOIRLab, 950 N. Cherry Ave., Tucson, AZ 85719, USA
19 Institució Catalana de Recerca i Estudis Avançats, Passeig de Lluís
Companys, 23, 08010 Barcelona, Spain
20 Department of Physics and Astronomy, Siena College, 515 Loudon Road,
Loudonville, NY 12211, USA
21 National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun
Rd., Chaoyang District, Beijing, 100012, P.R. China
22 Department of Physics and Astronomy, University of Waterloo, 200 University
Ave W, Waterloo, ON N2L 3G1, Canada
23 Perimeter Institute for Theoretical Physics, 31 Caroline St. North,
Waterloo, ON N2L 2Y5, Canada
24 Waterloo Centre for Astrophysics, University of Waterloo, 200 University
Ave W, Waterloo, ON N2L 3G1, Canada
25 Department of Physics, Kansas State University, 116 Cardwell Hall,
Manhattan, KS 66506, USA
26 Department of Physics and Astronomy, Sejong University, Seoul, 143-747,
Korea
27 CIEMAT, Avenida Complutense 40, E-28040 Madrid, Spain
28 Department of Physics & Astronomy, Ohio University, Athens, OH 45701, USA
29 University of Michigan, Ann Arbor, MI 48109, USA

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Accurate mock galaxy catalogues are crucial to validate analysis pipelines
used to constrain dark energy models. We present a fast HOD-fitting method
which we apply to the AbacusSummit simulations to create a set of mock
catalogues for the DESI Bright Galaxy Survey, which contain $r$-band
magnitudes and $(g-r)$ colours. The halo tabulation method fits HODs for
different absolute magnitude threshold samples simultaneously, preventing
unphysical HOD crossing between samples. We validate the HOD fitting procedure
by fitting to real-space clustering measurements and galaxy number densities
from the MXXL BGS mock, which was tuned to the SDSS and GAMA surveys. The
best-fitting clustering measurements and number densities are mostly within
the assumed errors, but the clustering for the faint samples is low on large
scales. The best-fitting HOD parameters are robust when fitting to simulations
with different realisations of the initial conditions. When varying the
cosmology, trends are seen as a function of each cosmological parameter. We
use the best-fitting HOD parameters to create cubic box and cut sky mocks from
the AbacusSummit simulations, in a range of cosmologies. As an illustration,
we compare the ${}^{0.1}M_{r}<-20$ sample of galaxies in the mock with BGS
measurements from the DESI one-percent survey. We find good agreement in the
number densities, and the projected correlation function is reasonable, with
differences that can be improved in the future by fitting directly to BGS
clustering measurements. The cubic box and cut-sky mocks in different
cosmologies are made publicly available.

###### keywords:

methods: analytical – large-scale structure of Universe – catalogues –
galaxies: statistics

††pubyear: 2023††pagerange: Generating mock galaxy catalogues for flux-limited
samples like the DESI Bright Galaxy Survey–Generating mock galaxy catalogues
for flux-limited samples like the DESI Bright Galaxy Survey

##  1 Introduction

The $\Lambda$CDM cosmological model has been very successful at describing the
formation and evolution of structure in the Universe (Planck Collaboration et
al., 2020). However, recent tensions have emerged, for example in measurements
of the Hubble parameter, between those derived from measurements of the cosmic
microwave background and measurements from supernovae in the local Universe
(e.g. Verde et al., 2019; Di Valentino et al., 2021; Freedman, 2021). In
addition, the nature of dark energy, which drives the accelerated expansion
(Riess et al., 1998; Perlmutter et al., 1999), and makes up the majority of
the energy density of the Universe, remains poorly understood.

These questions can be probed using the two-point clustering statistics of
galaxies in large galaxy surveys. The large-scale structure of the Universe
was seeded by primordial fluctuations at early times, which evolved to produce
the distribution of galaxies we observe today. Baryon acoustic oscillations
(BAO), which propagated through the early Universe, were frozen at the epoch
of recombination, leading to a characteristic distance scale in the clustering
of galaxies (Cole et al., 2005; Eisenstein et al., 2005). This can be used as
a standard ruler to measure the expansion history of the Universe. In
addition, the peculiar velocity of a galaxy along the line of sight has the
effect of shifting the observed redshift of the galaxy (Kaiser, 1987).
Measuring these redshift space distortions (RSD) in two-point galaxy
clustering statistics provides a test of general relativity, and constrains
modified gravity models (Guzzo et al., 2008).

The Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration et al.,
2016a, b, 2022) is currently undertaking a 5-year survey that will measure the
spectra of approximately 40 million galaxies and quasars between $0<z<3.5$,
and which aims to place our best constraints on models of dark energy. DESI is
targeting several galaxy tracers over this redshift range, including luminous
red galaxies (LRGs), emission line galaxies (ELGs) and quasars (QSOs). During
bright time, DESI is conducting the Bright Galaxy Survey (BGS), in addition to
the Milky Way Survey (MWS) of stars within the Milky Way galaxy. The BGS is a
flux-limited survey of low redshift galaxies ($z\sim 0.2$). The BGS-BRIGHT
sample will contain over 10 million galaxies brighter than $r=19.5$. The
secondary BGS-FAINT sample will extend this magnitude limit to $r=20.175$,
with additional cuts based on fibre magnitude and colour to ensure a high
redshift success rate (Hahn et al., 2023). The BGS will provide a highly
complete galaxy catalogue not only for cosmological measurements, but also for
studies of galaxy formation and evolution. The first DESI data from the survey
validation program (DESI Collaboration et al., 2023a) has recently been made
publicly available (DESI Collaboration et al., 2023b).

The use of accurate simulated mock galaxy catalogues are critical for large
surveys like DESI, to aid in survey design and assess survey strategies
(Looser et al., 2021). In addition, mocks are needed to test and optimise the
key data analysis pipelines that are designed to handle the large volumes of
data. Synthetic datasets can be used to test the recovery of key statistics
from survey data, investigating potential systematic errors and ensuring the
unbiased measurement of cosmological parameters. Since the volume probed by
modern galaxy surveys is so large, these simulations also need to cover
extremely large cosmological volumes. For covariance matrices, approximate or
low resolution simulations are used in order to generate many thousands of
mocks.

To create realistic mock galaxy catalogues, it is necessary to accurately
model the link between galaxies and their host dark matter halos. The galaxy-
halo connection has been a subject of research for several decades (Cole &
Kaiser, 1989; Mo & White, 1996; Wechsler & Tinker, 2018), with many
formulations from simple Halo Occupation Distributions (HODs) to more complex
treatments that aim to model the impact of various physical processes on
galaxy formation in hydrodynamic and semi-analytic simulations.

Hydrodynamic galaxy formation simulations model the galaxy-halo connection
directly by including both dark-matter and baryonic components interacting
over time directly through gravity and other forces (e.g. Vogelsberger et al.,
2014; Crain et al., 2015; McCarthy et al., 2017; Lee et al., 2021). However,
these simulations are very computationally expensive, limiting the volume that
can be simulated in a reasonable amount of time. While it is recently becoming
possible to run hydrodynamical simulations with Gpc box sizes (e.g. the
FLAMINGO simulations, Schaye et al., 2023) for the simulation volumes required
for large surveys like DESI, it is typically only feasible to run dark-matter-
only simulations, using methods to paint galaxies onto the dark matter halos.

Semi-analytic galaxy formation models can model the formation and evolution of
galaxies in an existing dark-matter-only $N$-body simulation. Physically-
informed models simulate a variety of processes such as star formation,
feedback, and radiative heating and cooling (e.g. Croton et al., 2006; Lacey
et al., 2016; Henriques et al., 2020). There are a large number of degrees of
freedom and assumptions underlying the physical models which are constrained
to match observations. While this can be applied to a large dark-matter-only
simulation, high resolution halo merger trees are required.

Sub-halo abundance matching (SHAM) techniques can be used to link galaxies to
sub-halos, based on a ranking of sub-halo and galaxy properties (e.g. Conroy
et al., 2006; Reddick et al., 2013; Prada et al., 2023). To reproduce galaxy
clustering statistics, the relationship between the sub-halo and galaxy
properties must include some scatter (Tasitsiomi et al., 2004). SHAM
techniques require high simulation resolution in order to resolve sub-halos
sufficiently, otherwise halos may be over-merged and systematic errors can be
produced in the clustering (Guo & White, 2014). The SHAM algorithm has been
extended to e.g. include the satellite fraction (Favole et al., 2016),
assembly bias (Contreras et al., 2021) and observational systematics (Yu et
al., 2022).

The HOD is a more empirical method for modelling the galaxy-halo connection.
Instead of assuming that galaxies are formed by specific physical processes
which affect their abundance and properties, HOD methods simply assume that
the abundance of galaxies is informed by the host halo properties without
reference to a physical model. This abstracts away much of the underlying
physics but is sufficient if one only cares about knowing what the galaxy-halo
connection is, rather than why it takes that form (Berlind & Weinberg, 2002).

It is relatively simple to use a HOD model for a single galaxy sample to add
galaxies to the dark matter halos of a simulation. However, for the BGS, we
also want to assign galaxy properties, like magnitudes. Previously, this has
been done in Smith et al. (2017, 2022b), using the MXXL simulation. A method
was developed to use a set of ‘nested’ HODs for different magnitude thresholds
to assign each mock galaxy a SDSS $r$-band magnitude. However, these SDSS HODs
were measured for each magnitude-threshold sample independently, and needed to
be modified to prevent the HODs from crossing over each other. If one galaxy
sample is a subset of another, it is unphysical for the two HOD curves to
cross. The average number of galaxies in halos of mass $M$ that are brighter
than a certain magnitude threshold must always be a monotonic function of
magnitude.

In this work, we present a method for simultaneously fitting HODs for multiple
magnitude-threshold samples, which prevents unphysical HOD crossing, and uses
both the clustering measurements and number densities of each sample as
constraints. As a proof of concept, we fit HODs to the AbacusSummit
simulations, fitting to clustering measurements obtained from the MXXL mock.
We describe the HOD-fitting methodology, and show the results from fitting to
MXXL clustering measurements. The best-fitting HODs are used to produce a set
of mock catalogues for the BGS, which we make publicly available. We also
discuss the limitations of the method, and the modifications and improvements
that are needed in future work to apply the method directly to BGS data. The
outline of this paper is as follows. In Section 2 we describe our HOD model
for linking BGS galaxies to dark matter halos. The tabulation method for fast
HOD fitting is described in Section 3. In Section 4, we fit HODs to the
AbacusSummit simulations, testing the robustness of the method and explore how
the HOD parameters vary for simulations with different cosmologies. The method
for creating cubic box and cut-sky mocks is described in Section 5. Finally,
we summarise our conclusions in Section 6.

##  2 Linking galaxies and dark matter

###  2.1 Halo Occupation Distribution

Galaxies are biased tracers of the underlying matter density field of the
Universe. In the halo model, it is assumed that all galaxies reside within
larger dark matter halos. The link between galaxies and halos can be modelled
using a halo occupation distribution, which specifies the average number of
galaxies within each halo. In its simplest form, the HOD is purely a function
of the halo mass, and the form of the HOD is often split up into two
components, to model the abundance of central and satellite galaxies. The
specific choice of form for the HOD may depend on the selection of the galaxy
sample being modelled. For samples such as luminous red galaxies (LRGs), where
there is a correlation (with scatter) between halo mass and stellar mass (or
galaxy luminosity), the HOD of central galaxies is modelled with a smoothed
step function that approaches one at high masses (Zheng et al., 2005). For
emission line galaxies (ELGs), which are star forming, the central HOD is
instead modelled with a Gaussian-like distribution for the central galaxy
occupation, decreasing at high masses since these halos often host red
elliptical central galaxies, with low star formation rates (e.g. Avila et al.,
2020). The occupation number of satellite galaxies is often modelled as a
power law, with a cut-off at low masses.

More complex versions of HOD methods exist to account for the relationship
between other halo properties and galaxy abundance, known as assembly bias.
Decorated HODs have been explored as methods to improve HODs, which add
dependence of the HOD on several more parameters such as formation time and
concentration (Hearin et al., 2016; Yuan et al., 2022a).

The central galaxy is usually assumed to be at the centre of the halo, and
there are different methods of modelling the placement of satellite galaxies.
The satellites are often assumed to follow a Navarro-Frenk-White (NFW) profile
(Navarro et al., 1996, 1997). Alternatively, in a N-body simulation, the dark
matter particles within each halo can be used as tracers for the satellites.

###  2.2 HOD parametrisation

In this work, we use a HOD model of the same form as Smith et al. (2017) to
model the HOD of DESI BGS galaxies, where there is a correlation between
galaxy luminosity and halo mass. This is closely related to Zheng et al.
(2005), and is described by five parameters in total: two for the central
galaxy component and three for the satellites. The occupation number for
central galaxies brighter than luminosity $L$ is modelled using a smoothed
step function as a function of halo mass, $M$,

$\left\langle N_{\mathrm{cen}}(>L\mid
M)\right\rangle=\frac{1}{2}\left[1+F\left(\frac{\log M-\log
M_{\min}(L)}{\sigma_{\log M}(L)}\right)\right],$ (1)

where $M_{\textrm{min}}$ is the position of the step and $\sigma_{\log M}$ the
width.

In the standard HOD formalism of Zheng et al. (2005), the shape of the smooth
step function is modelled as an error function. We replace this with an
equivalent function, $F(x)$, which is defined as
$F(x)=2\int_{0}^{x}S\left(x^{\prime}\right)dx^{\prime}$, where $S(x^{\prime})$
is a pseudo-Gaussian function spline kernel function,

$\displaystyle S(x)$
$\displaystyle=1-6\mathopen{|}x\mathclose{|}^{2}+6\mathopen{|}x\mathclose{|}^{3}$
$\displaystyle\mathopen{|}x\mathclose{|}\leq 0.5$
$\displaystyle=2(1-\mathopen{|}x\mathclose{|})^{3}$ $\displaystyle
0.5<\mathopen{|}x\mathclose{|}\leq 1$ (2) $\displaystyle=0$ $\displaystyle
1<\mathopen{|}x\mathclose{|}.$

This has the advantage that the tails of the pseudo-Gaussian function are
truncated to exactly zero, helping to prevent unphysical crossing of the HODs
of different absolute magnitude threshold samples.

The satellite HOD is a power law weighted by the central HOD,

$\left\langle N_{\mathrm{sat}}(>L\mid M)\right\rangle=\left\langle
N_{\mathrm{cen}}(>L\mid
M)\right\rangle\left(\frac{M-M_{0}(L)}{M_{1}^{\prime}(L)}\right)^{\alpha(L)},$
(3)

where $M_{0}$ is the low mass cutoff mass scale, $M_{1}^{\prime}$ the
normalisation, and $\alpha$ the power law slope.

###  2.3 Varying HOD parameters with magnitude

A central assumption in this implementation is that HOD parameters are defined
to smoothly vary with an absolute magnitude limit. This is a reasonable
assumption because we are not modelling a galaxy population that is sensitive
to an on/off state such as star formation. This allows us to define the HOD at
any magnitude and therefore to populate the mock galaxy catalogue with
galaxies that have realistic luminosities, rather than populating a mock that
represents a fixed magnitude cut. Fitting HOD parameters to smoothly varying
curves presents some new challenges and there are several ways in which this
could be implemented.

One option is to first independently fit HOD parameters for a set of different
magnitude threshold samples. Smooth curves can then be fit to each parameter,
as a function of magnitude. This method has the benefit of being relatively
quick and simple, as one only has to repeat the standard HOD fitting procedure
several times, once for each magnitude limit. However, there is no guarantee
that the best-fitting HOD parameter values can be well approximated by a
smooth curve as a function of magnitude. The HOD parameters produced by these
smooth functions also may not preserve the target galaxy clustering which has
also been fitted. This is the approach that was taken in Smith et al. (2017),
but adjustments were required to prevent unphysical crossing between HODs of
bright and faint magnitude threshold samples.

The fitting method chosen in this paper is to first parametrise how each HOD
parameter varies with magnitude before fitting. The meta-parameters describing
these functions are then fit, meaning that no further adjustment needs to be
made to the fitted HOD parameters. Each HOD parameter, and therefore the shape
of the total HOD, can be obtained for any magnitude limit, using this
parametrisation. One potential pitfall of this method is that if the
parameterised curves are too constrained, then a good fit to the target
clustering may not be found and therefore different forms for the curves may
need to be tested in order to produce robust results. Another issue is that
fitting all of the meta-parameters at once is a higher dimensional fit than
fitting parameters for individual magnitudes separately. This takes up more
computing time to complete the fit and the large number of parameters means
that the fitting procedure is more likely to end up in a local minimum than be
able to locate the global minimum. We attempted to overcome this by running
the fitting procedure several times from different starting points in the
parameter space. This provides more confidence that the best fit HOD
parameters are stable without needing to allocate vastly greater computational
resources to the problem.

The parametrised forms of the 5 HOD parameters are shown below, where each HOD
parameter is written as a function of absolute magnitude. In total, there are
17 parameters. The functional forms were chosen as they approximately match
the measured HOD parameters when HOD fits are done independently for each
magnitude threshold. $\log M_{0}$ is described by a linear function of
magnitude, $\sigma_{\log M}$ is described by a smoothly varying step function
and $\alpha$ is described by a constant value with an exponential component at
bright magnitudes. We describe the mass parameters $\log M_{\mathrm{min}}$ and
$\log M_{1}$ with cubic functions. This is a different functional form then
used by e.g. Zehavi et al. (2011), but we find that it produces better fits.
The five functions are given by

$\displaystyle\textrm{log}M_{\textrm{min}}$
$\displaystyle=12+A_{\mathrm{min}}+B_{\mathrm{min}}M_{r}^{\prime}+C_{\mathrm{min}}M_{r}^{\prime
2}+D_{\mathrm{min}}M_{r}^{\prime 3}$ (4) $\displaystyle\sigma_{\textrm{log}M}$
$\displaystyle=A_{\sigma}+\frac{B_{\sigma}-A_{\sigma}}{1+\exp({C_{\sigma}(M_{r}^{\prime}+D_{\sigma}))}}$
(5) $\displaystyle\textrm{log}M_{0}$
$\displaystyle=11+A_{0}+B_{0}M_{r}^{\prime}$ (6)
$\displaystyle\textrm{log}M_{1}$
$\displaystyle=12+A_{1}+B_{1}M_{r}^{\prime}+C_{1}M_{r}^{\prime
2}+D_{1}M_{r}^{\prime 3}$ (7) $\displaystyle\alpha$
$\displaystyle=A_{\alpha}+B_{\alpha}^{-M_{r}^{\prime}+C_{\alpha}},$ (8)

For conciseness we define $M_{r}^{\prime}={}^{0.1}M_{r}+20$, where
${}^{0.1}M_{r}$ is the $r$-band absolute magnitude, $k$-corrected to a
reference redshift $z=0.1$ (see Section 5.2). The HOD parameters as a function
of magnitude are illustrated in Figure 1. These 5 functions are similar to how
the SDSS HODs were parametrised in order to construct the MXXL mock catalogue
of Smith et al. (2017). A different functional form was used for
$M_{\mathrm{min}}$ and $M_{1}$, which only had 3 free parameters (Zehavi et
al., 2011). However, in this work we find that cubic functions with an extra
free parameter are able to produce better HOD fits.

Figure 1: The best-fitting HOD parameters, as a function of magnitude, for the
AbacusSummit Planck ‘c000’ cosmology box with initial conditions ‘ph000’,
illustrating the functional forms of Eqs. 4-8. The upper panel shows the mass
parameters $M_{\textrm{min}}$, $M_{0}$, and $M_{1}$ on a logarithmic scale, in
blue, orange and green, respectively. The lower panel shows the parameters
$\sigma_{\textrm{log}M}$ in red and $\alpha$ in purple. We plot $\alpha-1$ so
that both parameters fall within the same y-axis range.

###  2.4 Preventing HOD crossing

For two galaxy samples with lower luminosity limits/thresholds $L_{1}$ and
$L_{2}$ (with $L_{2}>L_{1}$), it must always be true that the number of
galaxies $N_{\mathrm{gal}}(>L_{1})\geq N_{\mathrm{gal}}(>L_{2})$, since the
galaxies with luminosity $L>L_{2}$ are a subset of those with $L>L_{1}$. It is
therefore unphysical for the HODs to cross, as this condition would no longer
be true, and would require a negative number of galaxies in the range
$L_{1}<L<L_{2}$. As already discussed, this motivates the choice of HOD
parametrisation, which uses a pseudo-Gaussian function to truncate the tail of
the central HODs at low masses.

This consideration can also be built into the HOD fitting procedure by adding
a high penalty to the likelihood if the HODs cross.

##  3 Halo tabulation

In order to efficiently explore the parameter space of possible HOD forms, a
fast method to evaluate the expected galaxy clustering is needed. We have
implemented a version of the halo paircounting method introduced by Neistein
et al. (2011) and Zheng & Guo (2016). This method tabulates the dark matter
halo paircounts by mass and separation, then the average effect of using a
particular HOD can be estimated by reweighting elements of these tables and
summing the contributions to the clustering.

Using this method means that the computationally expensive paircounting
routine needs only to be run once in total, instead of once for each set of
HOD parameters. The paircounting is run on the halos instead of the galaxies,
and the larger number of halos means that there is a greater fixed cost to
using this method. However when many evaluations are needed in the fitting
procedure, the much quicker time for each evaluation leads to a quicker
runtime in total.

Alternatives to explicitly calculating the clustering produced by a set of HOD
parameters can include using emulators or analytic methods (Kwan et al., 2015;
Zhai et al., 2019; Yuan et al., 2022b), but these are not guaranteed to
produce accurate results for a given simulation.

This tabulation method assumes that occupancy of a halo depends on only one
variable, the halo mass. We can bin halo paircounts by different combinations
of halo mass and then reweight these paircounts to account for changes in the
HOD. It is possible to include other halo parameters by mixing them with the
halo mass to create an _effective mass_ which maintains the HOD function as
dependent on a single variable (Yuan et al., 2018), however this method is not
tested here. HODs with two or more input variables can be used, but this
increases the dimension of the tabulated paircounts and leads to an
exponential increase in code runtime.

###  3.1 Tabulation Method

Here we give an overview of the tabulation method, but a more detailed
description can be found in Grove (2023), which includes additional run time
and accuracy tests.

We aim to estimate the real-space two-point correlation function, $\xi(r)$,
which can be written as

$\xi(r)=\frac{GG(r)}{RR(r)}-1,$ (9)

where $GG(r)$ and $RR(r)$ are the normalized galaxy-galaxy and random-random
pair counts respectively, in bins of separation $r$. For our periodic
simulation boxes, $RR(r)$ can be calculated analytically. We work in real
space since this simplifies the halo tabulation.

The galaxy pair counts can be split into components from central and satellite
galaxies,

$GG(r)=CC(r)+CS(r)+SS_{\mathrm{2halo}}(r)+SS_{\mathrm{1halo}}(r),$ (10)

where $CC$ and $CS$ represent the central-central, central-satellite galaxy
pair counts, respectively. The satellite-satellite pair counts are further
split into the 2-halo term, where the satellite occupy different halos,
$SS_{\mathrm{2halo}}$, and a 1-halo term of satellites that reside in the same
halo, $SS_{\mathrm{2halo}}$.

Each of these terms can be expressed as a sum over galaxy pairs which reside
in halos of different masses. For example, for the central-central pairs,

$CC(r)=\sum_{M_{ij}}{CC(M_{ij},r)},$ (11)

with similar expressions for the other terms in Eq. 10. Here, $CC(M_{ij},r)$
is the number of central-central pairs which live in halos of masses $M_{i}$
and $M_{j}$ (which we write as $M_{ij}$). We use evenly spaced logarithmic
mass bins, and pairs in the sum must not be double counted. Each of these sums
can be evaluated from the halo pair counts and the HOD model.

The HOD model we use, which is described in Section 2.2 and 2.3, provides the
mean number of central galaxies in halos of mass $M$, $N_{\mathrm{cen}}(M)$
(where $0\leq N_{\mathrm{cen}}(M)\leq 1$), and the mean number of satellite
galaxies, $N_{\mathrm{sat}}(M)$. We assume that there does not have to be a
central galaxy in a halo in order for it to host satellite galaxies. This
assumption simplifies the tabulation calculations significantly as otherwise
one would have to consider a larger number of correlations from halos
explicitly with and without central galaxies. In practice the majority of
satellite galaxies are placed in halos containing a central galaxy as they are
preferentially placed in high-mass halos. The low mass end of the HOD is
dominated by the central term, making it very unlikely to have satellites
without a central. Central galaxies are positioned at the centre of the halo
and satellite galaxies are positioned randomly according to an NFW profile
around the halo. The number of satellite galaxies per halo is chosen according
to a Poisson distribution with mean $N_{\rm sat}(M)$.

The central galaxies will share positions with the halos as they are placed at
the centre, therefore the distribution of potential central galaxy positions
is sampled by the halo centres themselves. For the satellite galaxies we place
a fixed number of tracer particles around each halo to represent where
satellite galaxies would be placed. In this work, we use 3 satellite tracers
per halo.

We represent the number of pairs of halos in mass bins $M_{i}$ and $M_{j}$ as
$W^{\rm CC}_{ij}(r)$. Using the satellite tracers for each halo, we similarly
define $W^{\rm CS}_{ij}(r)$, $W^{\rm SS1}_{ij}(r)$ and $W^{\rm SS2}_{ij}(r)$,
which are the halo-satellite pair counts, and the satellite-satellite counts
divided into 1-halo and 2-halo terms.

The galaxy central-central pair counts can then be evaluated from the sum

$CC(r)=\sum_{M_{ij}}N_{\rm cen}(M_{i})N_{\rm cen}(M_{j})W^{\rm CC}_{ij}(r).$
(12)

Similar expressions can be written for the central-satellite and satellite-
satellite 2-halo terms,

$CS(r)=\frac{1}{T}\sum_{M_{ij}}N_{\rm cen}(M_{i})N_{\rm sat}(M_{j})W^{\rm
CS}_{ij}(r)$ (13) $SS_{\mathrm{2halo}}(r)=\frac{1}{T^{2}}\sum_{M_{ij}}N_{\rm
sat}(M_{i})N_{\rm sat}(M_{j})W^{\rm SS2}_{ij}(r),$ (14)

where $T$ is the number of satellite tracers per halo. Note the different
dependence on the number of satellite tracer particles between Equation 13 and
14.

In practice, $T=1$ is sufficient for accurate calculation of both the central-
satellite term and the satellite-satellite two-halo term. For the one-halo
term, we must use more than one tracer particle per halo to sample the
possible galaxy pairs. The number of one-halo satellite galaxy pairs sampled
per halo when using $T$ tracer particles per halo is $T(T-1)/2$. Therefore in
this case we relate the galaxy and halo pair counts using

$SS_{\mathrm{1halo}}(r)=\frac{1}{T(T-1)}\sum_{M_{ij}}N_{\rm sat}(M_{i})N_{\rm
sat}(M_{j})W^{\rm SS1}_{ij}(r).$ (15)

Here all the off-diagonal terms are zero as there are no pairs of satellites
between halos of different masses in the one-halo term. In our calculations,
we use $T=3$ satellite tracers per halo.

In each case the halo paircounts are static and only depend on the halo
catalogue, meanwhile the HOD factors can vary allowing us to fit HOD
parameters rapidly by evaluating the sums shown above.

These estimated galaxy paircounts can be combined with analytic randoms to
produce the expected correlation function as shown in Eq. 9. The analytic
expression for the number of random pairs is

$RR(r)=\frac{\textrm{d}V(r)}{V}N(N-1),$ (16)

where $\textrm{d}V(r)$ is the volume of a spherical shell with minimum and
maximum radii corresponding to the edges of the radial bin in the clustering
calculation, $V$ is the total volume of the simulation box, and $N$ is the
expected number of galaxies in the whole simulation volume which can be
estimated by applying the HOD to the halo mass function. The number of pairs
of randoms vary as the HOD parameters change because the expected total number
of galaxies that will be populated is altered.

###  3.2 Halo Subsampling

Figure 2: Fraction of halos that are subsampled, $f$, as a function of halo
mass. Halos are subsamples in order to speed up the computation of the halo
pair counts. Subsampling is only done for halos with mass
$M<10^{12}~{}h^{-1}\mathrm{M_{\odot}}$; all halos are used above this. At
$M=10^{11}~{}h^{-1}\mathrm{M_{\odot}}$, the subsampled fraction is 1%.

The halos used in the paircounting were subsampled in order to speed up
computation of the tables of the binned paircounts. There are vastly more low
mass halos than high mass halos in the halo catalogue and low mass halos are
proportionally occupied by fewer galaxies. This means that HOD fits can be
robust to a subsampling of low mass halos. By testing of different subsampling
forms and mass cutoffs, it was found that the runtime could be improved by a
factor of ten with negligible impacts to the measured number density and
clustering. The form of the subsampling was a function of halo mass, with no
subsampling at masses above $10^{12}h^{-1}\textrm{M}_{\odot}$. The subsample
fraction $f$ as a function of halo mass $M$ is set as
$\log_{10}f=\mathrm{min}\\{0,2(\log_{10}M-12)\\}$, with $M$ in units
$h^{-1}\mathrm{M_{\odot}}$. This function is shown in Figure 2. We checked
with one of the simulations that changes to the number densities and
clustering of HOD fits introduced by the subsampling were smaller than 0.1%.
This is at a level much smaller than the precision we aim to achieve in the
HOD fits, so halo subsampling will not bias our results.

###  3.3 Unresolved Halos

Figure 3: The real-space 2-point clustering of field particles along with
halos of different masses from the base AbacusSummit simulation. Correlation
functions are truncated at small scales where there are insufficient pairs for
it to be measured. The clustering of the field particles is similar to that of
halos at the low mass limit of the halos, therefore field particles can
reasonably be used to extend the HOD below the mass resolution limit of the
halos.

A halo mass cut at $10^{11}~{}h^{-1}\mathrm{M_{\odot}}$ was applied to the
AbacusSummit halo catalogues because halos are not sufficiently resolved below
this limit. Randomly selected field particles were used as the locations of
halos below this mass cut and these were not given satellite tracers because
the satellite contribution from the HOD used in this work is negligible at low
masses. Figure 3 shows that using the field particles as tracers provides a
good approximation to the low mass, unresolved halos. The field particles have
a similar clustering amplitude as the low mass halos, and at low masses the
halo bias becomes approximately constant.

##  4 HOD fitting procedure

In this section we describe the HOD fitting procedure, which we apply to the
AbacusSummit simulations. The simulations and data used in this work are
described in Section 4.1, and the HOD fitting procedure in Section 4.2. The
best-fitting HODs in Planck cosmology are discussed in Section 4.3, and the
dependence of these results on cosmology in Section 4.4.

###  4.1 Simulations and data

The simulations and data used in this work are described below. We fit HODs
using the AbacusSummit simulations (Section 4.1.1) to clustering measurements
obtained from the MXXL mock catalogue (Section 4.1.2). The mocks constructed
from AbacusSummit are then compared to the DESI BGS One-percent survey data
(Section 4.1.3). The HODs are fit to MXXL, which in turn was fit to SDSS and
GAMA data (Smith et al., 2017), rather than the BGS data directly because the
one-percent survey is somewhat smaller, shallower and less complete than GAMA.
Direct comparisons between the mocks and data are enabled by using the same
magnitude selection, using the same $k$\- and $E$-corrections.

####  4.1.1 AbacusSummit simulations

AbacusSummit111https://abacussummit.readthedocs.io/(Maksimova et al., 2021) is
a suite of $N$-body simulations run using the abacus code (Garrison et al.,
2021) on the Summit supercomputer. There are 150 simulations that have been
run with different cosmologies, resolutions and initial conditions. In this
work we use the ‘base’ resolution simulations, which contain $6912^{3}$
particles in a cubic box of side length $2~{}h^{-1}\mathrm{Gpc}$, with a
particle mass of $2.11\times 10^{9}~{}h^{-1}\mathrm{M_{\odot}}$.

The AbacusSummit simulations were run from 2LPT initial conditions at $z=99$.
The primary ‘c000’ cosmology corresponds to Planck 2018 $\Lambda$CDM results
(Planck Collaboration et al., 2020), with cosmological parameters $h=0.6736$,
$\Omega_{\mathrm{cdm}}h^{2}=0.1200$, $\Omega_{\mathrm{b}}h^{2}=0.02237$,
$\sigma_{8}=0.8114$ and $n_{s}=0.9649$. A total of 25 boxes were produced in
the primary cosmology, with different initial conditions. The secondary
cosmologies are selected to match existing flagship $N$-body simulation
projects, in addition to a wide emulator grid around the primary cosmology. In
the emulator grid, the parameters $\sigma_{8}$, $\Omega_{\mathrm{cdm}}h^{2}$,
$n_{s}$, $\Omega_{\mathrm{b}}h^{2}$, $w_{0}$, $w_{a}$, $N_{ur}$ and
$\alpha_{s}$ are all varied around the Planck values, where $\alpha_{s}$ is a
running spectral index, $N_{ur}\sim 2$ is the effective number of massless
neutrinos, and $w_{0}$ and $w_{a}$ set an evolving dark energy equation of
state of the form $w(z)=w_{0}+w_{a}(1-a)$. These cosmologies include a single
species of massive 60 meV neutrinos, with
$\Omega_{\mathrm{ncdm}}h^{2}=6.4420\times 10^{-4}$. In the primary Planck
cosmology, $\alpha_{s}=0$, $N_{ur}=2.0328$, $w_{0}=-1$ and $w_{a}=0$
(Maksimova et al., 2021).

Halo catalogues are the primary data product from AbacusSummit, with the full
set of all particles being too large to store efficiently. Halo finding in
AbacusSummit was performed on the fly using CompaSO (Hadzhiyska et al.,
2022a), a spherical overdensity method. A friends-of-friends (Davis et al.,
1985) algorithm is first applied to identify ‘L0’ halos. ‘L1’ halos are found
using a a spherical overdensity algorithm, which finds particles within a
density threshold $\Delta=200$, while ‘L2’ subhalos are found using
$\Delta=800$. For each halo, we use use a halo mass which is defined as the
total number of particles within the L1 halo (corresponding to
$M_{\mathrm{200m}}$). The L2 cores of each halo are used to find the positions
and velocities. The CompaSO algorithm can deblend halos, producing unphysical
objects. We use the cleaned CompaSO halo catalogues which are processed using
merger tree information to fix these issues (see Hadzhiyska et al., 2022a).

Subsets of the particles are available at certain snapshots. This is split
into an ‘A’ sample (3% of the particles) and a ‘B’ sample (7% of the
particles). We use the random subset of field particles (which do not exist in
halos) to extend the halo catalogue below the mass resolution of the
simulation (see Section 3.3).

####  4.1.2 MXXL BGS mock

The Millennium-XXL (MXXL) simulation (Angulo et al., 2012) is a dark-matter-
only simulation that is a successor to the Millennium simulation (Springel et
al., 2005), with a much larger box size of $3~{}h^{-1}\mathrm{Gpc}$. The
particle mass is $6.17\times 10^{9}~{}h^{-1}\mathrm{Mpc}$, and the simulation
was run in the same WMAP1 cosmology as the original Millennium simulation,
with $\Omega_{\mathrm{m}}=0.25$, $\Omega_{\mathrm{b}}=0.045$
$\Omega_{\Lambda}=0.75$, $\sigma_{8}=0.9$, $h=0.73$ and $n_{s}=1$ (Spergel et
al., 2003). At each simulation snapshot, dark matter halos are first
identified using a friends-of-friends algorithm, with bound substructures
identified using subfind. To construct a halo merger tree, the 15 most bound
particles of each halo were found, and the descendent is the halo at the next
snapshot containing the majority of these particles.

The MXXL simulation was previously used to create a mock galaxy catalogue for
the BGS, as described in Smith et al. (2017, 2022b). Halos were first
interpolated between simulation snapshots to build a halo lightcone, and were
then populated with galaxies using a set of ‘nested’ HODs for different
$r$-band magnitude thresholds, which were measured from the SDSS survey. The
HODs were evolved with redshift to reproduce an evolving $r$-band target
luminosity function from SDSS and GAMA. A Monte Carlo method is used to assign
each galaxy a luminosity from the nested HODs, and satellites are positioned
following a NFW profile. The galaxies are also randomly assigned $g-r$ colours
from a parametrisation of the GAMA colour-magnitude diagram. In addition to
creating a lightcone mock, the same HOD methodology is applied to several
simulation snapshots to cubic box mocks at different redshifts.

In this work, we obtain our target galaxy clustering measurements from the
MXXL BGS mock. We measure the real-space correlation function from the cubic
box mock made from the $z=0.2$ snapshot. This snapshot is chosen as this is
the median redshift of the DESI BGS. Since we use the full cubic box, the
volume is large, and we can obtain precise measurements that are not affected
by cosmic variance.

####  4.1.3 DESI BGS One-percent survey

The DESI One-percent survey was conducted at the end of survey validation, and
before the start of the main 5-year survey. The one-percent survey observed a
footprint composed of 20 circular ‘rosettes’, covering an area of $\sim
140~{}\mathrm{deg}^{2}$ to very high completeness. The same target selection
as the main survey was used, with similar exposure times (DESI Collaboration
et al., 2023a).

Fibre assignment completeness is corrected for in the catalogue by applying a
completeness weight, $w_{\mathrm{comp}}$. This weight is determined from 128
alternative realisations of the fibre assignment algorithm, in addition to the
real survey. For each galaxy, this weight is given by

$w_{\mathrm{comp}}=\frac{129}{N_{\mathrm{assigned}}+1},$ (17)

where $N_{\mathrm{assigned}}$ is the number of alternative realizations in
which the target was assigned a fibre (DESI Collaboration et al., 2023b). Note
that this individual inverse probability weighting is not unbiased on very
small scales. However, the one-percent survey data is highly complete, so
these weights are typically close to 1. A pairwise inverse probability (PIP)
weighting can in principle provide an unbiased correction (Bianchi & Percival,
2017; Smith et al., 2019; Bianchi & Verde, 2020; Mohammad et al., 2020), and
will be provided in future DESI data releases (Lasker et al., 2023). The BGS
one-percent catalogues also provide FKP weights, $w_{\mathrm{FKP}}$, for the
complete BGS-BRIGHT sample, which reduce the variance in the correlation
function measurements where there are variations in the galaxy density
(Feldman et al., 1994). We do not apply any FKP weighting to our BGS
clustering measurements, since we use magnitude threshold galaxy samples with
a constant number density.

In Section 5.3, we compare the clustering of our AbacusSummit mock with
measurements from the One-Percent survey. In this work, we do not fit HODs
directly to this dataset, since the volume is small, and there are large
fluctuations in the clustering due to cosmic variance. In future work, we will
adapt the HOD fitting method to apply it directly to the larger DESI Year 1
dataset. This is discussed further in Section 5.4.

###  4.2 Fitting method

HOD fitting was performed using the emcee code, a Markov chain Monte Carlo
sampler (Foreman-Mackey et al., 2013). All 17 HOD meta-parameters were fitted
at once with target clustering and number densities taken from the BGS mock
catalogue described in Smith et al. (2017) at absolute $r$-band magnitude
limits in intervals of 0.5 from -18 to -22. These clustering and number
density values were themselves tuned to fit the results from the Sloan Digital
Sky Survey (SDSS) (Abazajian et al., 2009) & the Galaxy and Mass Assembly
(GAMA) survey (Liske et al., 2015) during creation of the original mock
catalogue (Smith et al., 2017). The halo catalogue corresponding to the
$z=0.2$ snapshot was used in all cases.

Clustering bias becomes scale-independent on large scales; this causes a
problem if the shape of the correlation function of the simulation and target
data have different shapes, due to the difference in cosmology. This can never
be fit well and therefore throws off the full HOD parameter fit. We adjust the
target data to match the shape of the simulation large-scale correlation
function by using the Zel’dovich clustering prediction at scales above
$8~{}h^{-1}\textrm{Mpc}$. This is evaluated from the Zel’dovich matter power
spectrum using the nbodykit package (Hand et al., 2018). The rescaled
correlation function is given by

$\xi_{\mathrm{rescaled}}(r)=\xi_{\mathrm{MXXL}}(r)\frac{\xi^{\mathrm{Zel}}_{\mathrm{Abacus}}(r)}{\xi^{\mathrm{Zel}}_{\mathrm{MXXL}}(r)}\frac{\xi^{\mathrm{Zel}}_{\mathrm{MXXL}}(r=8)}{\xi^{\mathrm{Zel}}_{\mathrm{Abacus}}(r=8)},$
(18)

where the superscript ‘Zel’ indicates the Zel’dovich correlation function, and
subscripts ‘MXXL’ and ‘Abacus’ indicate the cosmologies. On small scales
$r<8~{}h^{-1}\mathrm{Mpc}$, no rescaling is applied.

In addition, a volume correction was used to ensure that the target luminosity
function was reproduced despite the changing size of volume elements in
different cosmologies. The magnitude thresholds are first shifted to

${}^{0.1}M_{r}^{\mathrm{rescaled}}={}^{0.1}M_{r}+5\log_{10}\left(\frac{D_{L}^{\mathrm{Abacus}}}{D_{L}^{\mathrm{MXXL}}}\right),$
(19)

where $D_{L}$ is the luminosity distance to $z=0.2$, computed in the two
cosmologies. The target luminosity function (which is the same luminosity
function used in the construction of the MXXL mock) is then used to obtain the
number density of galaxies brighter than each rescaled magnitude threshold,
$\Phi(>{}^{0.1}M_{r}^{\mathrm{rescaled}})$. This number density is then
rescaled by the ratio of the volume element of the two cosmologies,

$\Phi^{\mathrm{rescaled}}=\Phi\cdot\frac{\mathrm{d}V_{\mathrm{MXXL}}}{\mathrm{d}V_{\mathrm{Abacus}}}$
(20)

The volume element in each cosmology is

$\mathrm{d}V=\frac{c(1+z)^{2}D_{A}^{2}(z)}{H(z)}\mathrm{d}\Omega\mathrm{d}z$
(21)

The clustering fit was limited to a maximum scale of
$50~{}h^{-1}\textrm{Mpc}$. At scales larger than this, BAO features become
important and this shape cannot be fit to by modifying HOD parameters as it is
an intrinsic feature of the cosmology.

The likelihood function used in the emcee fitting contained a contribution
from both the clustering and the number density expected from the HOD
parameters,

$\mathcal{L}(M)=-\frac{1}{2}\left(\sum_{r}{\left(\frac{\xi(r,M)-\xi_{\mathrm{t}}(r,M)}{\sigma_{\xi}(r)}\right)^{2}}+\left(\frac{n(M)-n_{\mathrm{t}}(M)}{\sigma_{n}}\right)^{2}\right)$
(22)

where $\mathcal{L}(M)$ is the likelihood for a particular magnitude limited
sample $M$. The subscript $\mathrm{t}$ is used to represent the target values
for the clustering, $\xi$, and number density, $n$. The total length of the
data vector is 677, where for each of the 9 magnitude thresholds
($M_{r}<-22.0,-21.5,...,-18.0$) there is a number density, and 74 correlation
function bins between $0.01<r<50~{}h^{-1}\mathrm{Mpc}$. In addition, to
improve the luminosity function at faint magnitudes, number density
constraints are also included for $M_{r}<-17.5$ and $M_{r}<-17.0$.

The total magnitude of the likelihood function is not important for the
position of the best fit HOD parameters, but it can affect how emcee explores
the parameter space, if the likelihood function is too steep then the walkers
can become stuck in local maxima and if it is too shallow then they may not
converge on the best fit solution. Through testing we were able to establish
values for $\sigma_{\xi}(r)$ and $\sigma_{n}$ which produced robust fits with
the correct balance between fitting the clustering and the number density.
$\sigma_{\xi}(r)$ and $\sigma_{n}$ were the same for all nine of the magnitude
fits. We choose to use a constant fractional error for $\sigma_{\xi}(r)$,
which is 6.5625 %. This value was chosen to avoid overfitting to noise at very
small scales. For $\sigma_{n}$, we also use a constant fractional error of
0.9375 %. This is chosen to avoid overfitting to either the number density or
the correlation function. We did not use a full error covariance matrix when
finding the best-fit HOD parameters. This is because we were not trying to
establish the likelihood surface for different HODs describing the data, but
instead trying to find one set of HOD parameters which produces a good fit to
the desired clustering measurements and luminosity function.

###  4.3 Best-fitting HODs

####  4.3.1 Best-fitting Meta-parameters

Figure 4: A corner plot from the emcee fitting chain of the 17 meta-parameters
describing the smoothly varying HOD curves. Best-fit values are shown by the
vertical and horizontal dashed lines. This example comes from the primary
AbacusSummit cosmology. This figure was generated using the pygtc package
(Bocquet & Carter, 2016).

The output of the fitting procedure is not the HOD parameters directly, but a
set of meta-parameters describing the evolution of smooth curves describing
the HOD parameters as a function of absolute magnitude limit. Figure 4 shows a
corner plot of all 17 parameters described in Equations 4-8 in the case of the
primary AbacusSummit Planck cosmology. The corner plot shows the parameter
space explored by the emcee fitting chain after an appropriate amount of burn-
in is discarded. The best-fit values are indicated by the dashed lines.
Correlations between parameters indicate that degeneracies exist as the
parameters vary along that axis. Note that since we did not use a full
covariace, caution should be used when interpreting these degeneracies.

The majority of the meta-parameters are well constrained, except for the
parameter $A_{\sigma}$, which sets the minimum value of $\sigma_{\log M}$ at
faint magnitudes, and $B_{\alpha}$ and $C_{\alpha}$, which describe the
exponential increase of $\alpha$ at bright magnitudes. Reference to Equations
5 and 8 explain the lack of constraint, because changes in these parameters
have very little effect on the shape of the HOD curve in certain regimes. In
addition, there is a long negative tail on $A_{0}$ and $B_{0}$. This is
because $M_{0}$ does not affect the shape of the HOD curve if it is at
sufficiently low mass that it applies when satellites have little impact.

####  4.3.2 Best-Fitting HOD Curves

The shape of the HOD curves that are produced from the best-fit meta-
parameters in the Planck cosmology case is shown in the upper panel of Figure
5, for galaxy samples with different absolute magnitude thresholds. The
contributions from both the central and satellite HOD components are combined
into a single curve. There are no unphysical overlaps between the HOD curves
as this is disallowed by the fitting procedure. The best-fitting HODs are
shown for each of the 25 Planck cosmology AbacusSummit boxes, showing that the
results are very consistent, with little scatter between them.

####  4.3.3 Comparison to Target Data

Figure 5: Top panel: best-fitting HODs of the 25 simulations with the same
c000 Planck cosmology, but with different initial conditions. Each line
represents a single HOD fit. The HOD curves are highly consistent with one
another for all the samples, with significant differences emerging only at the
low mass tails. Bottom panel: Variation of the best-fit HODs when using
AbacusSummit simulations with different cosmologies. Each line represents a
single HOD fit. The HOD variation is larger than the case where cosmology is
held constant as in the upper panel. We only plot the HODs for AbacusSummit
cosmologies in the emulator grid, to show the scatter around the primary
cosmology. Figure 6: Top panel: the ratio of the best-fitting correlation
function, $\xi(r)$, to the target MXXL clustering measurements, for the
primary ‘c000’ Planck cosmology. Each curve is for a different magnitude
threshold sample, as indicated in the legend. The grey shaded region is the
constant fractional error assumed during the HOD fitting procedure. Bottom
panel: The ratio of the predicted and target number densities of the same
magnitude threshold samples. The shaded regions shows the fractional
uncertainties assumed when fitting.

In addition to exploring the values of the best-fitting meta-parameters and
shape of the HOD curves produced from them, we have also compared the quality
of the best-fit with regards to the target data. Since we have fit to target
data for both the clustering and the number density simultaneously, the
quality of both of these aspects of the fit should be investigated.

The upper panel of Figure 6 shows the ratio of the predicted two-point
clustering for the c000 Planck cosmology to the target measurements for the
magnitude threshold samples used in the fitting procedure, where the shaded
region indicated the constant fractional error that was assumed. The
clustering measurements are mostly within this error. However, there is an
offset on large scales, where the best-fitting clustering is systematically
lower than the target clustering, by up to 10% for the faintest magnitude
threshold sample. The features above $\sim 50~{}h^{-1}\mathrm{Mpc}$ are
outside the range of the fitting because these are caused by BAO which cannot
be mitigated by varying HOD parameters. Different behaviour is seen in the
residuals on small scales, depending on the magnitude limit of the sample. The
brightest sample displays a peak relative to the target clustering at $\sim
2~{}h^{-1}\mathrm{Mpc}$, with a dip $\sim 0.4~{}h^{-1}\mathrm{Mpc}$, and this
behaviour is inverted for the faintest samples. Below
$0.1~{}h^{-1}\mathrm{Mpc}$, the slope of the best-fitting correlation function
is different to the target, leading to differences of up to 20% at
$0.01~{}h^{-1}\mathrm{Mpc}$.

In the lower panel of Figure 6, the ratio of the predicted number density to
the target value is shown as a function of the magnitude limit of the sample,
with the fractional error indicated by the shaded area. There is good
agreement for all the samples that were included in the fitting procedure,
where the number density for almost all the samples is within the error.

###  4.4 Varying the cosmology

Before describing the production of the mock galaxy catalogues, we shall
firstly test the robustness of the HOD fitting procedure, firstly by applying
it to multiple simulations in the same cosmology, and then explore the
variation of HOD parameters with cosmology. The HOD fitting procedure
described in Section 4.2 was applied to all base resolution AbacusSummit
boxes, including all 25 c000 Planck cosmology boxes, and the boxes in other
cosmologies. When applying the method to the different cosmologies, different
cosmology rescaling factors must be applied to the target number densities and
correlation functions.

An initial check that we performed was to find the level of sample variance
that exists for HOD fits on multiple independent simulation boxes in the same
cosmology. There are 25 AbacusSummit simulations in the base Planck
$\Lambda$CDM cosmology (see Section 4.1.1) that were produced from initial
conditions with different phase information. If there are large differences in
best-fit HODs produced by the fitting procedure on these boxes, then it
implies that there is a lack of robustness in the fitting procedure.

The best-fitting HODs are shown in the upper panel Figure 5. There is a high
degree of consistency between the HODs. Differences between HODs emerge at the
low mass tails where there are very few galaxies per halo. HOD differences at
these scales have a small effect on the properties of galaxy catalogues
because of the low occupation. Certain points have very little scatter. These
locations are where the amplitude of the HOD has the greatest effect on the
number density, which is constrained in the fitting process.

Creating mocks in different cosmologies allows us to test the impact of
cosmological parameters on mock observations. Cosmological parameter recovery
is an important application of mock catalogues that requires mocks in
different cosmologies. Exploring how the HOD parameters vary with cosmology
informs us how our models of the galaxy-halo connection may change.

The lower panel of Figure 5 shows the variation in HOD shapes using
simulations with different cosmologies. This variation is larger than the
sample variance when the same cosmology is used. The cosmology emulator grid
is sufficiently wide to lead to significant differences in the best-fit HODs.
This justifies our choice to refit the HODs for simulations that used
different cosmologies.

Figure 7: Left panel: best-fitting HOD parameter $\alpha$ vs $\sigma_{8}$,
from the AbacusSummit cosmologies in an emulator grid around the primary
Planck cosmology, for different magnitude threshold samples. A negative
correlation can be seen. Right panel: best-fitting $\alpha$ parameter vs
$\Omega_{\mathrm{m}}$. Here, a positive correlation is seen.

There is correlation between the best-fitting HOD parameters and the
cosmological parameters of the simulation. As an example, the left panel of
Figure 7 shows $\alpha$ plotted against $\sigma_{8}$ for different magnitude
threshold samples. A clear negative correlation can be seen, where $\alpha$ is
$\sim 1$ for the highest values of $\sigma_{8}$, but is larger for small
$\sigma_{8}$. In addition, this correlation becomes steeper for the brightest
samples, but with a larger scatter. $\sigma_{8}$ sets the amplitude of the
initial power spectrum at a fiducial scale of $8~{}h^{-1}\mathrm{Mpc}$.
Simulations run with higher $\sigma_{8}$ will produce a halo catalogue with
stronger clustering. Therefore, in order to reproduce the same target
clustering, the galaxy clustering bias must be reduced in these cases. This
can be achieved by reducing the number of satellite galaxies, because they are
located in the highest mass (and hence most biased) halos. Therefore it makes
intuitive sense that $\sigma_{8}$ should correlate with lower $\alpha$, if
$M_{0}$ and $M_{1}$ do not change significantly. For the brightest samples, it
is likely that the increase in the slope and scatter is due to the weak
constraints on the meta-parameters $B_{\alpha}$ and $C_{\alpha}$. Our HOD
model has too much freedom, since these are only constrained by the brightest
few samples. In the future, we will modify the HOD model to reduce the number
of meta-parameters for $\alpha$.

A second example in the right panel of Figure 7, which shows $\alpha$ plotted
against the matter density $\Omega_{\mathrm{m}}$. Here, a positive correlation
is seen, where in cosmologies with higher $\Omega_{\mathrm{m}}$, larger values
of $\alpha$ are obtained. There are many degeneracies between the different
HOD parameters, so it can be more revealing to see the effect of cosmology on
the full occupation functions, rather than an individual parameter.

Figure 8 shows how the shape of the best-fitting HOD curves is affected by
$\sigma_{8}$ and $\Omega_{\mathrm{m}}$. The HODs are the same as in the lower
panel of Figure 5. We only show the ${}^{0.1}M_{r}<-22$, -21 and -18 samples
for clarity, but all samples show the same trends. For the parameter
$\sigma_{8}$, higher values are associated with fewer satellite galaxies, and
have a broader central step function. As already discussed, halos in
simulations with high $\sigma_{8}$ are more strongly clustered, and reducing
the number of satellites and placing more centrals in lower mass halos both
have the effect of reducing the clustering of the galaxy catalogue. For
simulations with large values of the parameter $\Omega_{\mathrm{m}}$, the HODs
are shifted to higher masses, with a sharper step function for the central
galaxies. Even though $\alpha$, the satellite power law slope, as shown in
Figure 7 is higher for large $\Omega_{\mathrm{m}}$, the shift in the HODs to
higher masses results in a total number of satellites which is smaller.

There are too many potential combinations of cosmological and HOD parameters
to investigate all the relationships in this way. These variations in HOD
parameters with cosmology may be useful if one wanted to build an emulator for
HOD parameters. Such an emulator could produce an estimated best-fit HOD for
any cosmology that lies within the region spanned by the AbacusSummit
cosmologies.

Figure 8: Top panel: best-fitting HOD curves from the AbacusSummit simulations
in different cosmologies, for the ${}^{0.1}M_{r}<-22$, -21 and -18 samples,
from right to left. Each curve is coloured by the simulation value of
$\sigma_{8}$, with low values in blue and high values in red. Bottom panel: as
above, but the curves are coloured by the value of $\Omega_{\mathrm{m}}$.

##  5 Mock creation

In this section we provide a brief overview of the methodology for creating
mocks from the AbacusSummit simulations. Galaxies are first positioned within
halos in the cubic box simulations. The nested set of HODs that we have
determined are used to assign each galaxy a luminosity, as in Smith et al.
(2017). In addition, the semi-empirical method of Smith et al. (2017, 2022b)
is used to assign each galaxy a rest frame $(g-r)$ colour. We then process the
cubic box mocks to produce ‘cut-sky’ mocks, which convert the positions of the
galaxies into sky coordinates, producing a catalogue that is more realistic
and similar to what DESI observes. The mocks we create are full sky, which can
then be cut to cover the DESI survey footprint. For each AbacusSummit
simulation, we create one mock catalogue, using a single random realisation of
the best-fitting HODs which reproduce the target clustering and luminosity
function.

###  5.1 Cubic Box Mock Catalogues

Following the Monte Carlo method of Smith et al. (2017), we use our best-
fitting HODs to add galaxies to halos in the $z=0.2$ cubic box, assigning
$r$-band absolute magnitudes. This is described in section 4.1 of Smith et al.
(2017). Since each HOD parameter is described as a smooth function of absolute
magnitude, we can easily evaluate the HOD for any magnitude threshold, which
is necessary for this method. The best-fitting HODs are also constrained to
prevent any unphysical crossing of the HODs of different magnitude thresholds,
which is also important when using the method to assign magnitudes. Each
central galaxy is then placed at the centre of the halo, and assigned the same
velocity, which is defined in AbacusSummit as the location of the centre of
mass of the largest sub-halo. The number of satellite galaxies above a minimum
luminosity threshold is drawn from a Poisson distribution. Magnitudes are
assigned following Smith et al. (2017), and the satellite galaxies are
positioned following a NFW profile, with a random virial velocity along each
direction drawn from a Gaussian with velocity dispersion
$\sigma^{2}=GM_{\mathrm{200m}}/({2R_{\mathrm{200m}}})$.

The HOD parameters $M_{\textrm{min}}$ and $M_{1}$ are modelled as cubic
functions, which means that they diverge rapidly at magnitudes beyond the
fitting range. To avoid this, these parameters are both extrapolated linearly
(in $\log M$) beyond this range.

In addition to assigning galaxy luminosities, we also assign rest frame
$(g-r)$ colours to the galaxies. This allows colour cuts to be applied to the
mock and colour dependent clustering to be investigated. We assign colours
following the method of Smith et al. (2022b), which uses a parametrisation of
the colour-magnitude diagram of the GAMA survey. This process for assigning
colours is an extension to the method of Skibba & Sheth (2009); Smith et al.
(2017). It is first randomly decided whether a galaxy should lie on the red or
blue sequence, which is different for central and satellite galaxies. A random
colour is then drawn from a Gaussian distribution, from the fit to the GAMA
colour-magnitude diagram. The assigned colour depends on the absolute
magnitude and redshift of each galaxy, and therefore galaxies in more massive
haloes are more likely to be red.

The final cubic box mock is cut to an absolute magnitude limit of
${}^{0.1}M_{r}<-18$, which corresponds to a galaxy number density of
$3.7\times 10^{-2}~{}(h^{-1}\mathrm{Mpc})^{-3}$ for the c000 Planck cosmology
simulations.

###  5.2 Cut-sky mocks

To create a mock that is more representative of what DESI will observe, we
convert the cubic box mock from the $z=0.2$ snapshot into a ‘cut-sky’ mock.

An observer is placed at the corner of the box, and the cubic box is
replicated to cover the volume required to make a mock to a maximum redshift
of $z=0.6$. Around half of the volume of the final mock fits inside a single
copy of the cubic box, and it is above $z=0.37$ where the same large scale
structure is repeated. Cartesian coordinates are then converted to an angular
position on the sky and redshift, where the observed redshift of each galaxy
includes the effect of the peculiar velocity along the line of sight.

The cubic box we use to construct the mock is at $z=0.2$, which used HODs that
were fit to reproduce the target luminosity function at the same redshift
$z=0.2$. However, the BGS covers a wide redshift range of $0<z<0.6$, so we
must model the evolution of the galaxy luminosity function in the cut-sky
catalogue.

The target luminosity function we aim to reproduce in the mock comes from
existing SDSS and GAMA survey measurements. For $z>0.15$, the target
luminosity function is a Schechter function fit to the GAMA luminosity
function (Loveday et al., 2012). Below $z=0.15$, the target luminosity
function smoothly interpolates to that from SDSS measurements (Blanton et al.,
2003). The evolution of the Schechter parameters as a function of redshift is
described as

$\displaystyle M^{*}(z)$ $\displaystyle=M^{*}(z_{0})-Q(z-z_{0})$ (23)
$\displaystyle\phi^{*}(z)$ $\displaystyle=\phi^{*}(0)10^{0.4Pz},$ (24)

where $P$ sets the evolution in number density, and $Q$ sets the magnitude
evolution (Lin et al., 1999). We use the values $z_{0}=0.1$, $P=1.8$, and
$Q=0.7$, which were measured from the GAMA survey (McNaught-Roberts et al.,
2014).

A rescaling is applied to the absolute magnitudes to reproduce the evolving
target luminosity function, as is done in Dong-Páez et al. (2022). In narrow
redshift bins, we measure the cumulative luminosity function in the mock, and
assign new magnitudes from the target cumulative luminosity function at that
redshift, which correspond to the same number density for each galaxy. We also
re-run the colour assignment algorithm on the cut-sky mock to ensure that the
colour distributions also evolve smoothly with redshift, reproducing the GAMA
measurements.

The full DESI BGS sample contains faint galaxies at low redshifts which live
within low mass halos that fall below the resolution of the AbacusSummit
simulations. We add these unresolved halos to the cut-sky mock by using the
field particles (which are not in halos) as tracers. Halo masses are generated
from extrapolating the measured AbacusSummit halo mass function to low masses
(assuming a power law), and these halos are assigned the same position and
velocity of randomly-selected field particles. The unresolved halos are
assigned galaxies with magnitudes at $z=0.2$ using the same methods as the
resolved halos, and the same rescaling is used to add redshift evolution.

The cut-sky mock we have created contains the rest-frame absolute magnitude,
${}^{0.1}M_{r}$ and the rest-frame ${}^{0.1}(g-r)$ colour, which we then
convert to the observed quantities. The absolute magnitudes assigned to each
galaxy are converted to an apparent magnitude following

$r=M_{r}+5\log_{10}D_{\textrm{L}}(z)+25+k_{r}(z),$ (25)

where $D_{L}(z)$ is the luminosity distance from the observer to the galaxy at
redshift, $z$, in units of $h^{-1}\mathrm{Mpc}$. $k_{r}(z)$ is the $r$-band
$k$-correction, which accounts for the shift in the band pass with redshift.
We $k$-correct to a reference redshift of $z_{\mathrm{ref}}=0.1$.

The exact form of the $k$-correction depends on the filter being used and the
type of object being observed. In this work, the $k$-corrections we use come
from the GAMA survey and are described by fourth-order polynomials. The
$k$-corrections are split into seven different rest frame $g-r$ colour bins
with interpolation between the bins (see section 4.3 of Smith et al., 2017).

Similarly, galaxy colours can be transformed into the observer frame using

$g-r={}^{0.1}(g-r)+{}^{0.1}k_{g}(z)-{}^{0.1}k_{r}(z),$ (26)

where $k_{g}(z)$ is the $g$-band $k$-correction. As with the $r$-band, we use
a set of colour-dependent polynomial $k$-corrections.

Finally, an apparent magnitude cut of $r<20.2$ is applied, which encompasses
both the BGS-BRIGHT and BGS-FAINT samples.

###  5.3 Illustration: comparison with the DESI one-percent survey

In this section, as an illustration, we compare number density and clustering
measurements for a sample of BGS-BRIGHT ($r<19.5$) galaxies with
${}^{0.1}M_{r}<-20$ from the c000 Planck cosmology AbacusSummit mock with
measurements from the DESI one-percent survey. In these comparisons, the same
evolutionary correction has been applied to the absolute magnitudes in the
data and mock, which is of the form $E(z)=Q(z-z_{0})$, where $Q=0.97$ and
$z_{0}=0.1$ (McNaught-Roberts et al., 2014).

The number density as a function of redshift is shown in the upper panel of
Figure 9, for galaxies with ${}^{0.1}M_{r}<-20$ in the AbacusSummit mock, MXXL
mock and DESI one-percent survey. At low redshifts, since an evolutionary
correction has been applied the number density is almost constant with
redshift. Above $z\sim 0.25$, the density decreases since the sample become
incomplete due to the $r=19.5$ cut in apparent magnitude. As expected, we see
good agreement in the number densities between the AbacusSummit mock and the
MXXL mock, which have the same luminosity function. A small offset is seen at
high redshifts, but these mocks have different cosmologies. Both mocks agree
well with the number $n(z)$ measured from the DESI one-percent survey, however
the data error bars are large since the volume is much smaller than the full-
sky mocks.

Figure 9: Upper panel: the number density of BGS-BRIGHT ($r<19.5$) galaxies
brighter than ${}^{0.1}M_{r}=-20.0$, as a function of redshift, for galaxies
in the full-sky AbacusSummit Planck cosmology mock (purple line), the MXXL
mock (dashed black line), and DESI one-percent survey measurements (points
with error bars). An evolutionary correction has been applied to all absolute
magnitudes, and error bars are jackknife errors with 20 jackknife regions. A
weight is applied to the data measurements to take into account systematics
and incompleteness. Lower panel: projected correlation function of a volume
limited sample of galaxies brighter than ${}^{0.1}M_{r}=-20.0$, in the
redshift range $0.05<z<0.25$.

We measure the projected correlation function for a volume limited sample of
BGS-BRIGHT galaxies with ${}^{0.1}M_{r}<-20$, for the same MXXL and
AbacusSummit mocks and DESI one-percent survey data. The same evolutionary
correction is applied, and we cut to the redshift range $0.05<z<0.25$, where
the sample is complete (and the $n(z)$ is flat). The projected correlation
function is defined as

$w_{p}(r_{p})=2\int_{0}^{\pi_{\mathrm{max}}}\xi(r_{p},\pi)\mathrm{d}\pi,$ (27)

where $\xi(r_{p},\pi)$ is the correlation function in bins of $r_{p}$ and
$\pi$, perpendicular and parallel to the line of sight, respectively. Mocks
and data are integrated to $\pi_{\mathrm{max}}=40~{}h^{-1}\mathrm{Mpc}$. The
current AbacusSummit mocks are not tuned to reproduce the $w_{p}(r_{p})$
measurements of the BGS data, so we do not expect a perfect agreement.
However, this comparison demonstrates the current level of agreement, and
gives insights into improvements that could be made to our HOD modelling for
future mocks tuned directly to the BGS.

The projected correlation function is shown in the lower panel of Figure 9. On
intermediate scales, there is good agreement between the MXXL and AbacusSummit
mocks. On large scales, there is an offset between them. The two mocks have
different cosmologies, and therefore the clustering on large scales is
different. This was taken into account when fitting the HODs to the real-space
$\xi(r)$ measurements, but here we show the $w_{p}(r_{p})$ measurements
without any cosmology rescaling. Agreement with the one-percent survey
measurements is reasonable, but not perfect. On large scales, the data is more
strongly clustered than the AbacusSummit mock, and shows better agreement with
MXXL. On small scales, both mocks show stronger clustering than the the one-
percent survey measurements. The DESI clustering measurements are affected by
fibre incompleteness on small scales, since it is not possible to place a
fibre on every galaxy, particularly in dense regions such as large galaxy
clusters. However, in the one-percent survey, each region is observed with
multiple passes, so the completeness is very high, and a weighting is applied
to correct for any small incompleteness. On these small scales, the data
clustering measurements are therefore accurate, and the MXXL mock, which was
tuned to SDSS and GAMA is more strongly clustered. Since AbacusSummit is fit
to the MXXL clustering measurements, it also shows similar clustering to MXXL.
This could potentially be improved by modifying the parameters of the NFW
profiles used to position the satellites in the mock. Alternatively, fitting
HODs directly to BGS clustering measurements could produce fits with fewer
satellites, which would also reduce the small-scale clustering.

###  5.4 Limitations and improvements for DESI

We have presented a set of HOD mock catalogues for the DESI BGS, produced from
the AbacusSummit simulations. The HODs have been fit to real-space $\xi(r)$
measurements from the MXXL mock, and to number densities from the SDSS and
GAMA surveys. The aim of this work is a proof of concept that the fast HOD
fitting method, using halo tabulation, can be used to fit HODs to multiple
magnitude threshold samples simultaneously. We have validated that the best-
fitting $\xi(r)$ and number densities are mostly within the errors assumed.
Using these HODs to create a mock, we find good agreement in the number
densities compared to the DESI one-percent survey, but there are differences
in the clustering measurements which we aim to improve in future work by
fitting the HODs directly to BGS data.

The DESI survey will soon be releasing data from the first full year of the
survey. This large dataset will be ideal for tuning our future mocks, covering
a much larger area on the sky than the one-percent survey, with smaller
uncertainties in the number density and clustering measurements. However,
currently there are limitations to our HOD fitting method, which we aim to
address in future work.

In the current mock, we fit to real-space $\xi(r)$ measurements. In the real
data, this is not available, and we will instead fit to the projected
correlation function $w_{p}(r_{p})$. To do this, the HOD fitting method must
be modified to compute halo pair counts in bins of $r_{p}$ and $\pi$, and the
correlation function $\xi(r_{p},\pi)$ can be computed and integrated to obtain
$w_{p}(r_{p})$. This adds an extra dimension to the arrays of halo pair
counts, slowing down each step in the MCMC chain. The effect of velocities
also need to be included in the $\pi$-direction. While the increase in
dimensionality reduces the speed, this can be compensated for by reducing the
number of mass bins. We have checked that our fits are robust when reducing
the total number of mass bins from 120 to 30.

Our HOD model assumes that the number of galaxies in each halo depends only on
the halo mass. Recently, Pearl et al. (2023) detected an assembly bias signal
in counts-in-cylinders measurements from the BGS one-percent survey data. In
future work, the HOD model could be extended to include assembly bias,
enabling a cross-check of these results.

Currently when fitting the HODs, we assume a constant fraction error, both in
the correlation function and number density measurements. This can be improved
by using the uncertainties in the DESI measurements, which can be estimated
e.g. using jackknife sub sampling. However, this does not take into account
the covariances, so the fitting could be improved further by using a
covariance matrix. For the larger DESI Y1 dataset, covariance matrices will be
produced using analytic and mock-based methods (e.g. Rashkovetskyi et al.,
2023; Trusov et al., 2023)

To compute absolute magnitude from the DESI BGS data, we use colour-dependent
polynomial $k$-corrections from the GAMA survey. These colour-dependent
$k$-corrections are used to convert the observer frame SDSS $r$\- and
$g$-bands to the rest frame bands at $z_{\mathrm{ref}}$. We use these
$k$-corrections since they can easily be applied to the mock galaxies where we
only have ${}^{0.1}(g-r)$ colours and do not have a full spectrum. However,
these $k$-corrections are not sufficient for DESI, since there are differences
between the DECam and SDSS bands, and the photometry is different in the north
and south (Dey et al., 2019; Zarrouk et al., 2022). For individual DESI
galaxies, the $k$-correction can be computed from the spectrum using the
fastspecfit code (Moustakas et al.,
2023).222https://github.com/desihub/fastspecfit Colour-dependent
$k$-corrections can then be determined, to create a set of $k$-corrections
that are appropriate for the DESI BGS survey.

The cut-sky mocks we have created are constructed from a single simulation
snapshot at $z=0.2$. While we are able to reproduce the right evolving
luminosity function by applying a rescaling to the magnitudes in the mock,
there is no evolution in the underlying dark matter halos. In the AbacusSummit
simulations, lightcone outputs were produced on the fly, which have been used
to create halo lightcones (Hadzhiyska et al., 2022b). These can be used in the
future to make lightcone mocks which avoid issues when using snapshots to
build approximate lightcones, such as joining together the snapshots in shells
(Smith et al., 2022a), or interpolating halo properties between snapshots
(Smith et al., 2022b).

##  6 Conclusions

For large galaxy surveys such as DESI, it is essential to use realistic mock
galaxy catalogues to ensure the analyses are able to make unbiased
cosmological measurements, by testing how well key statistics are recovered,
and assessing systematics. A common method to create mock galaxy catalogues
from large dark-matter-only simulations is to add galaxies to halos using a
HOD model.

The DESI BGS is a survey of low redshift galaxies with median redshift $z\sim
0.2$, consisting of the flux-limited BGS-BRIGHT sample, with $r$-band apparent
magnitude limit $r=19.5$. The secondary BGS-FAINT sample extends this limit to
$r=20.175$, with additional colour and fibre magnitude cuts to ensure a high
redshift success rate. As was previously done for the MXXL simulation, mock
galaxy catalogues with $r$-band magnitudes can be created using a set of
‘nested’ HODs for different absolute magnitude thresholds. However, when
fitting the HODs of sample independently to data clustering measurements, it
is possible for the HODs of different magnitude thresholds to cross
unphysically.

We have developed a fast HOD fitting method to simultaneously fit HODs for
multiple magnitude threshold samples. By fitting all the samples at once, we
constrain the HODs to prevent unphysical crossing of the HODs.. Halo pair
counts are first tabulated as a function of halo mass, which only needs to be
done once. The correlation function for a given HOD model can then be
evaluated quickly from a weighted sum of the halo pair counts. Each HOD
parameter is defined as a smooth function of absolute magnitude, and we fit
the meta-parameters defining these smooth curves. When fitting, we include
constraints on both the galaxy clustering and number densities for each of the
magnitude threshold samples.

As a proof of concept, we apply the HOD fitting procedure to the AbacusSummit
simulations, using the snapshot at $z=0.2$, and fitting to real-space
correlation functions, $\xi(r)$ from the previous MXXL mock, and number
densities from the SDSS and GAMA surveys. Since the MXXL mock is in a WMAP1
cosmology, which is different to the AbacusSummit simulation cosmologies, we
apply a rescaling to the clustering measurements on scales
$r>8~{}h^{-1}\mathrm{Mpc}$ so that the large-scale clustering of the two
simulations matches. A cosmology rescaling is also applied to the target
number densities. We first apply the HOD fitting to the 25 AbacusSummit
simulations in the primary Planck cosmology. We find that there is very little
scatter between the best-fitting HODs of the 25 simulations, which verifies
the robustness of the fitting procedure. Most of the meta-parameters are well
constrained, and the parameters with the weakest constraints only have little
effect on the shape of the HODs. The number densities achieved are within the
assumed errors, and the clustering measurements are mostly within the errors,
although the clustering is low on large scales for the faintest samples.

We also apply the HOD fitting procedure to the AbacusSummit simulations in a
range of different cosmologies. Varying the cosmology produces much more
variation in the HODs, and there are trends with the different cosmological
parameters. For example, increasing the parameter $\sigma_{8}$ increases the
amplitude of the initial power spectrum and hence the halo catalogues in these
cosmologies are more strongly clustered. To match the same target galaxy
correlation function, the best-fitting HODs have fewer satellites, and a
broader central step function. For the parameter $\Omega_{\mathrm{m}}$, the
best-fitting HODs in cosmologies with high $\Omega_{\mathrm{m}}$ are shifted
to higher masses, and have a sharper central step.

We use the best-fitting HODs to create AbacusSummit mock catalogues for the
DESI BGS. We first populate the cubic box at $z=0.2$ with galaxies, using the
nested set of HODs to assign absolute $r$-band magnitudes following the same
method used to create the MXXL mock. A Monte Carlo method is also used to
assign ${}^{0.1}(g-r)$ colours. The snapshot is then converted to a ‘cut-sky’
mock by replicating the box, and converting the galaxy coordinates to a sky
coordinate and redshift, where the redshift takes into account the velocity of
the galaxy along the observer’s line of sight. Field particles are used as
tracers of halos which fall below the mass resolution of the simulation. A
rescaling is applied to the magnitudes to achieve a smoothly evolving
luminosity function, and colours are re-assigned so that the colour
distributions also evolve smoothly. The $r$-band apparent magnitude and
observer frame $g-r$ colours are calculated using colour-dependent
$k$-corrections from the GAMA survey.

As an illustration, we compare the number density and clustering of one of the
base Planck cosmology cut-sky mocks with BGS measurements from the one-percent
survey, for a sample of galaxies with magnitude threshold ${}^{0.1}M_{r}<-20$.
The $n(z)$ shows good agreement, although the uncertainties in the BGS
measurements are large, due to the small volume of the one-percent survey. The
projected correlation function, $w_{p}(r_{p})$ measurement in the mock is
reasonable, but the data is more strongly clustered than the mock on large
scales, while the opposite is true on small scales. In the future, we aim to
improve these mocks by extending the tabulation method to fit directly to
$w_{p}(r_{p})$ measurements from the DESI year 1 dataset, which covers a much
larger area than the one-percent survey. We also aim for future mocks to
improve the errors and take into account covariances, and to use DESI
$k$-corrections.

The AbacusSummit mocks we have presented here are being used within the DESI
collaboration as part of the first generation of BGS mocks, which we make
publicly available.

## Acknowledgements

AS, CG, SC and PN acknowledge the support of STFC grants ST/T000244/1 and
ST/X001075/1. Significant work presented in this paper is part of Cameron
Grove’s PhD thesis ‘Gravity to Galaxies – N-body Simulations for the DESI
Survey’ (Grove, 2023).

This material is based upon work supported by the U.S. Department of Energy
(DOE), Office of Science, Office of High-Energy Physics, under Contract No.
DE–AC02–05CH11231, and by the National Energy Research Scientific Computing
Center, a DOE Office of Science User Facility under the same contract.
Additional support for DESI was provided by the U.S. National Science
Foundation (NSF), Division of Astronomical Sciences under Contract No.
AST-0950945 to the NSF’s National Optical-Infrared Astronomy Research
Laboratory; the Science and Technology Facilities Council of the United
Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation;
the French Alternative Energies and Atomic Energy Commission (CEA); the
National Council of Science and Technology of Mexico (CONACYT); the Ministry
of Science and Innovation of Spain (MICINN), and by the DESI Member
Institutions: https://www.desi.lbl.gov/collaborating-institutions. Any
opinions, findings, and conclusions or recommendations expressed in this
material are those of the author(s) and do not necessarily reflect the views
of the U. S. National Science Foundation, the U. S. Department of Energy, or
any of the listed funding agencies.

The authors are honored to be permitted to conduct scientific research on
Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the
Tohono O’odham Nation.

## Data Availability

The AbacusSummit cubic box and cut-sky mocks produced in this work will be
made available at https://icc.dur.ac.uk/data/.

The data points and Python scripts for reproducing the figures in this
publication will be made available on Zenodo.

## References

  * Abazajian et al. (2009) Abazajian K. N., et al., 2009, ApJS, 182, 543
  * Angulo et al. (2012) Angulo R. E., Springel V., White S. D. M., Jenkins A., Baugh C. M., Frenk C. S., 2012, MNRAS, 426, 2046
  * Avila et al. (2020) Avila S., et al., 2020, MNRAS, 499, 5486
  * Berlind & Weinberg (2002) Berlind A. A., Weinberg D. H., 2002, ApJ, 575, 587
  * Bianchi & Percival (2017) Bianchi D., Percival W. J., 2017, MNRAS, 472, 1106
  * Bianchi & Verde (2020) Bianchi D., Verde L., 2020, MNRAS, 495, 1511
  * Blanton et al. (2003) Blanton M. R., et al., 2003, ApJ, 592, 819
  * Bocquet & Carter (2016) Bocquet S., Carter F. W., 2016, The Journal of Open Source Software, 1
  * Cole & Kaiser (1989) Cole S., Kaiser N., 1989, MNRAS, 237, 1127
  * Cole et al. (2005) Cole S., et al., 2005, MNRAS, 362, 505
  * Conroy et al. (2006) Conroy C., Wechsler R. H., Kravtsov A. V., 2006, ApJ, 647, 201
  * Contreras et al. (2021) Contreras S., Angulo R. E., Zennaro M., 2021, MNRAS, 508, 175
  * Crain et al. (2015) Crain R. A., et al., 2015, MNRAS, 450, 1937
  * Croton et al. (2006) Croton D. J., et al., 2006, MNRAS, 365, 11
  * DESI Collaboration et al. (2016a) DESI Collaboration et al., 2016a, arXiv e-prints, p. arXiv:1611.00036
  * DESI Collaboration et al. (2016b) DESI Collaboration et al., 2016b, arXiv e-prints, p. arXiv:1611.00037
  * DESI Collaboration et al. (2022) DESI Collaboration et al., 2022, AJ, 164, 207
  * DESI Collaboration et al. (2023a) DESI Collaboration et al., 2023a, arXiv e-prints, p. arXiv:2306.06307
  * DESI Collaboration et al. (2023b) DESI Collaboration et al., 2023b, arXiv e-prints, p. arXiv:2306.06308
  * Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
  * Dey et al. (2019) Dey A., et al., 2019, AJ, 157, 168
  * Di Valentino et al. (2021) Di Valentino E., et al., 2021, Classical and Quantum Gravity, 38, 153001
  * Dong-Páez et al. (2022) Dong-Páez C. A., et al., 2022, arXiv e-prints, p. arXiv:2208.00540
  * Eisenstein et al. (2005) Eisenstein D. J., et al., 2005, ApJ, 633, 560
  * Favole et al. (2016) Favole G., et al., 2016, MNRAS, 461, 3421
  * Feldman et al. (1994) Feldman H. A., Kaiser N., Peacock J. A., 1994, ApJ, 426, 23
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Freedman (2021) Freedman W. L., 2021, ApJ, 919, 16
  * Garrison et al. (2021) Garrison L. H., Eisenstein D. J., Ferrer D., Maksimova N. A., Pinto P. A., 2021, MNRAS, 508, 575
  * Grove (2023) Grove C., 2023, PhD thesis, Durham University, http://etheses.dur.ac.uk/14821/
  * Guo & White (2014) Guo Q., White S., 2014, MNRAS, 437, 3228
  * Guzzo et al. (2008) Guzzo L., et al., 2008, Nature, 451, 541
  * Hadzhiyska et al. (2022a) Hadzhiyska B., Eisenstein D., Bose S., Garrison L. H., Maksimova N., 2022a, MNRAS, 509, 501
  * Hadzhiyska et al. (2022b) Hadzhiyska B., Garrison L. H., Eisenstein D., Bose S., 2022b, MNRAS, 509, 2194
  * Hahn et al. (2023) Hahn C., et al., 2023, AJ, 165, 253
  * Hand et al. (2018) Hand N., Feng Y., Beutler F., Li Y., Modi C., Seljak U., Slepian Z., 2018, AJ, 156, 160
  * Hearin et al. (2016) Hearin A. P., Zentner A. R., van den Bosch F. C., Campbell D., Tollerud E., 2016, MNRAS, 460, 2552
  * Henriques et al. (2020) Henriques B. M. B., Yates R. M., Fu J., Guo Q., Kauffmann G., Srisawat C., Thomas P. A., White S. D. M., 2020, MNRAS, 491, 5795
  * Kaiser (1987) Kaiser N., 1987, MNRAS, 227, 1
  * Kwan et al. (2015) Kwan J., Heitmann K., Habib S., Padmanabhan N., Lawrence E., Finkel H., Frontiere N., Pope A., 2015, ApJ, 810, 35
  * Lacey et al. (2016) Lacey C. G., et al., 2016, MNRAS, 462, 3854
  * Lasker et al. (2023) Lasker J., et al., 2023, in prep.
  * Lee et al. (2021) Lee J., et al., 2021, ApJ, 908, 11
  * Lin et al. (1999) Lin H., Yee H. K. C., Carlberg R. G., Morris S. L., Sawicki M., Patton D. R., Wirth G., Shepherd C. W., 1999, The Astrophysical Journal, 518, 533
  * Liske et al. (2015) Liske J., et al., 2015, MNRAS, 452, 2087
  * Looser et al. (2021) Looser T. J., Lilly S. J., Sin L. P. T., Henriques B. M. B., Maiolino R., Cirasuolo M., 2021, MNRAS, 504, 3029
  * Loveday et al. (2012) Loveday J., et al., 2012, MNRAS, 420, 1239
  * Maksimova et al. (2021) Maksimova N. A., Garrison L. H., Eisenstein D. J., Hadzhiyska B., Bose S., Satterthwaite T. P., 2021, MNRAS, 508, 4017
  * McCarthy et al. (2017) McCarthy I. G., Schaye J., Bird S., Le Brun A. M. C., 2017, MNRAS, 465, 2936
  * McNaught-Roberts et al. (2014) McNaught-Roberts T., et al., 2014, MNRAS, 445, 2125
  * Mo & White (1996) Mo H. J., White S. D. M., 1996, MNRAS, 282, 347
  * Mohammad et al. (2020) Mohammad F. G., et al., 2020, MNRAS, 498, 128
  * Moustakas et al. (2023) Moustakas J., et al., 2023, in prep.
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Neistein et al. (2011) Neistein E., Li C., Khochfar S., Weinmann S. M., Shankar F., Boylan-Kolchin M., 2011, MNRAS, 416, 1486
  * Pearl et al. (2023) Pearl A. N., et al., 2023, arXiv e-prints, p. arXiv:2309.08675
  * Perlmutter et al. (1999) Perlmutter S., et al., 1999, ApJ, 517, 565
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Prada et al. (2023) Prada F., et al., 2023, arXiv e-prints, p. arXiv:2306.06315
  * Rashkovetskyi et al. (2023) Rashkovetskyi M., et al., 2023, MNRAS, 524, 3894
  * Reddick et al. (2013) Reddick R. M., Wechsler R. H., Tinker J. L., Behroozi P. S., 2013, ApJ, 771, 30
  * Riess et al. (1998) Riess A. G., et al., 1998, AJ, 116, 1009
  * Schaye et al. (2023) Schaye J., et al., 2023, arXiv e-prints, p. arXiv:2306.04024
  * Skibba & Sheth (2009) Skibba R. A., Sheth R. K., 2009, MNRAS, 392, 1080
  * Smith et al. (2017) Smith A., Cole S., Baugh C., Zheng Z., Angulo R., Norberg P., Zehavi I., 2017, MNRAS, 470, 4646
  * Smith et al. (2019) Smith A., et al., 2019, MNRAS, 484, 1285
  * Smith et al. (2022a) Smith A., Cole S., Grove C., Norberg P., Zarrouk P., 2022a, MNRAS, 516, 1062
  * Smith et al. (2022b) Smith A., Cole S., Grove C., Norberg P., Zarrouk P., 2022b, MNRAS, 516, 4529
  * Spergel et al. (2003) Spergel D. N., et al., 2003, ApJS, 148, 175
  * Springel et al. (2005) Springel V., et al., 2005, Nature, 435, 629
  * Tasitsiomi et al. (2004) Tasitsiomi A., Kravtsov A. V., Wechsler R. H., Primack J. R., 2004, ApJ, 614, 533
  * Trusov et al. (2023) Trusov S., et al., 2023, arXiv e-prints, p. arXiv:2306.16332
  * Verde et al. (2019) Verde L., Treu T., Riess A. G., 2019, Nature Astronomy, 3, 891
  * Vogelsberger et al. (2014) Vogelsberger M., et al., 2014, MNRAS, 444, 1518
  * Wechsler & Tinker (2018) Wechsler R. H., Tinker J. L., 2018, ARA&A, 56, 435
  * Yu et al. (2022) Yu J., et al., 2022, MNRAS, 516, 57
  * Yuan et al. (2018) Yuan S., Eisenstein D. J., Garrison L. H., 2018, MNRAS, 478, 2019
  * Yuan et al. (2022a) Yuan S., Garrison L. H., Hadzhiyska B., Bose S., Eisenstein D. J., 2022a, MNRAS, 510, 3301
  * Yuan et al. (2022b) Yuan S., Garrison L. H., Eisenstein D. J., Wechsler R. H., 2022b, MNRAS, 515, 871
  * Zarrouk et al. (2022) Zarrouk P., et al., 2022, MNRAS, 509, 1478
  * Zehavi et al. (2011) Zehavi I., et al., 2011, ApJ, 736, 59
  * Zhai et al. (2019) Zhai Z., et al., 2019, ApJ, 874, 95
  * Zheng & Guo (2016) Zheng Z., Guo H., 2016, MNRAS, 458, 4015
  * Zheng et al. (2005) Zheng Z., et al., 2005, ApJ, 633, 791