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

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 Mr<-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.


INTRODUCTION
The Λ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. 2016aCollaboration et al. ,b, 2022) ) is currently undertaking a 5-year survey that will measure the spectra of approximately 40 million galaxies and quasars between 0 <  < 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 ( ∼ 0.2).The BGS-BRIGHT sample will contain over 10 million galaxies brighter than  = 19.5.The secondary BGS-FAINT sample will extend this magnitude limit to  = 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;Cooray & Sheth 2002;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 galaxyhalo 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 -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).
Another method is the conditional luminosity function (CLF), which was introduced by Yang et al. (2003).The CLF describes the halo occupation statistics in terms of galaxy luminosity, by modelling the luminosity function of galaxies residing in halos of mass .The CLF has previously been applied to SDSS data in order to obtain cosmological constraints (van den Bosch et al. 2013;Cacciato et al. 2013;More et al. 2013).
The HOD is a more empirical method for modelling the galaxyhalo 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).The HOD technique was applied to SDSS data in Zehavi et al. (2011).
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.This can be done using a set of 'nested' HODs, as has been done previously in Skibba et al. (2006); Smith et al. (2017); Paul et al. (2019); Smith et al. (2022b).In Smith et al. (2017Smith et al. ( , 2022b) ) a method based on Skibba et al. (2006) was developed to assign each mock galaxy a SDSS -band magnitude using a set of 'nested' HODs for different magnitude thresholds.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  that are brighter than a certain magnitude threshold must always be a monotonic function of magnitude.In Paul et al. (2019), nested HOD for different magnitude thresholds were fitted to SDSS clustering measurements, which were used to construct SDSS mocks with multi-band luminosities (Paranjape et al. 2021).
In the analysis of the DESI BGS data, a set of accurate mock galaxy catalogues are required that reproduce the BGS luminosity function and magnitude-dependent clustering.The mocks we produce here are a first step towards achieving this aim, by developing a mock creation pipeline based on existing datasets, which can then be extended to the real DESI data as it becomes available.In addition to validating the models used in the standard two-point clustering analyses, these mocks can be used for assessing new statistics that can be measured in the BGS, e.g.multi-tracer and three-point statistics.The mocks we make are created from AbacusSummit, which is a set of simulations that includes boxes run in a range of different cosmologies.This allows us to study how the galaxy halo connection varies with cosmology, in addition to testing the impact on the cosmological analyses.These mocks are made publicly available, and have been used as part of the DESI year-1 BAO (DESI Collaboration et al. 2024b) and RSD analyses (DESI Collaboration et al. 2024a).
To create the mocks, we simultaneously fit HODs to galaxy clustering and number density constraints from multiple magnitudethreshold samples, preventing any unphysical crossing of the HODs between samples.We fit to measurements from the previous MXXL BGS mock catalogue, which provides an accurate estimate of the expected clustering of the BGS survey.The volume of the 1% survey is too small and affected by cosmic variance for HOD fitting, but in future work we will apply this method to the larger DESI year-1 dataset.Our goal is to measure a set of HODs which can be used to create mocks that reproduce the clustering and number density of the BGS.We describe the HOD-fitting methodology, and show the results from fitting to MXXL clustering measurements.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.

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 (Paranjape et al. 2015;Hearin et al. 2016;Yuan et al. 2022a).For example, (Alam et al. 2024) found some evidence of assembly bias in the GAMA data, using an extended HOD model developed in Alam et al. (2021).
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(Navarro et al. , 1997)).Alternatively, in a N-body simulation, the dark matter particles within each halo can be used as tracers for the satellites.

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  is modelled using a smoothed step function as a function of halo mass, , where  min is the position of the step and  log  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,  (), which is defined as which has been normalised and rescaled to mean  = 0 and variance  2 = 1/2.This is done following 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, where  0 is the low mass cutoff mass scale,  ′ 1 the normalisation, and  the power law slope.

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 bestfitting 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 metaparameters 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  0 is described by a linear function of magnitude,  log  is described by a smoothly varying step function and  is described by a constant value with an exponential component at bright magnitudes.We describe the mass parameters log  min and log  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.5-9.The upper panel shows the mass parameters  min ,  0 , and  1 on a logarithmic scale, in blue, orange and green, respectively.The lower panel shows the parameters  log in red and  in purple.We plot  − 1 so that both parameters fall within the same y-axis range.
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 For conciseness we define  ′  = 0.1   + 20, where 0.1   is the band absolute magnitude, -corrected to a reference redshift  = 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  min and  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.

Preventing HOD crossing
For two galaxy samples with lower luminosity limits/thresholds  1 and  2 (with  2 >  1 ), it must always be true that the number of galaxies  gal (>  1 ) ≥  gal (>  2 ), since the galaxies with luminosity  >  2 are a subset of those with  >  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  1 <  <  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.
An alternative approach was adopted in (Paul et al. 2019), where HODs were fit to magnitude-threshold samples from the SDSS data.Here, rather than fitting directly to the clustering of magnitudethreshold samples, the fits were done to clustering measurements in magnitude bins.A joint likelihood was constructed which enforces the monotonicity of the HODs statistically.

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.The Zheng & Guo (2016) tabulation method was also used in the HOD fitting of Paul et al. (2019).
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.

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,  (), which can be written as where  () and () are the normalized galaxy-galaxy and random-random pair counts respectively, in bins of separation .For our periodic simulation boxes, () 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, where  and  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,  2halo , and a 1-halo term of satellites that reside in the same halo,  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 centralcentral pairs, with similar expressions for the other terms in Eq. 11.Here,  (   , ) is the number of central-central pairs which live in halos of masses   and   (which we write as    ).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 ,  cen () (where 0 ≤  cen () ≤ 1), and the mean number of satellite galaxies,  sat ().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  sat ().
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   and   as  CC   ().Using the satellite tracers for each halo, we similarly define  CS   (),  SS1   () and  SS2   (), 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 Similar expressions can be written for the central-satellite and satellite-satellite 2-halo terms, where  is the number of satellite tracers per halo.Note the different dependence on the number of satellite tracer particles between Equation 14 and 15.
In practice,  = 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  tracer particles per halo is  ( − 1)/2.Therefore in this case we relate the galaxy and halo pair counts using 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  = 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. 10.The analytic expression for the number of random pairs is where d () is the volume of a spherical shell with minimum and maximum radii corresponding to the edges of the radial bin in the clustering calculation,  is the total volume of the simulation box, and  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.See section 4.5.3 of Grove (2023) for tests of the tabulation accuracy.

Halo Subsampling
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 ℎ −1 M ⊙ .The subsample fraction  as a function of halo mass  is set as log 10  = min{0, 2(log 10  − 12)}, with  in units ℎ −1 M ⊙ .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.

Unresolved Halos
A halo mass cut at 10 11 ℎ −1 M ⊙ was applied to the AbacusSummit halo catalogues we use in this work (see Section 4.1.1)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 the clustering of the field particles, compared to the resolved halos of the simulation.The field particles have a similar clustering amplitude as the 10 11 ℎ −1 Mpc halos, showing that using the field particles is a reasonable approximation for extending the HOD below the mass resolution limit of the halos.While the clustering amplitude decreases further for halos less massive than 10 11 ℎ −1 Mpc, only a small fraction of BGS galaxies reside in such low mass halos.

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.

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 -and -corrections.

AbacusSummit simulations
AbacusSummit 1 (Maksimova et al. 2021) is a suite of -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 ℎ −1 Gpc, with a particle mass of 2.11 × 10 9 ℎ −1 M ⊙ .
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 friendsof-friends (Davis et al. 1985) algorithm is first applied to identify 'L0' halos.'L1' halos are found using a a spherical overdensity algorithm, 1 https://abacussummit.readthedocs.io/which finds particles within a density threshold Δ = 200, while 'L2' subhalos are found using Δ = 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  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).

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 ℎ −1 Gpc.The particle mass is 6.17 × 10 9 ℎ −1 Mpc, and the simulation was run in the same WMAP1 cosmology as the original Millennium simulation, with Ω m = 0.25, Ω b = 0.045 Ω Λ = 0.75,  8 = 0.9, ℎ = 0.73 and   = 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. (2017Smith et al. ( , 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 -band magnitude thresholds, which were measured from the SDSS survey.The HODs were evolved with redshift to reproduce an evolving -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  −  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  = 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.

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 ∼ 140 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,  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 where  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. 2024).The BGS one-percent catalogues also provide FKP weights,  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.

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 -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  = 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 ℎ −1 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 where the superscript 'Zel' indicates the Zel'dovich correlation function, and subscripts 'MXXL' and 'Abacus' indicate the cosmologies.On small scales  < 8 ℎ −1 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  rescaled  = 0.1   + 5 log 10 where   is the luminosity distance to  = 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, Φ(> 0.1  rescaled

𝑟
).This number density is then rescaled by the ratio of the volume element of the two cosmologies, The volume element in each cosmology is The clustering fit was limited to a maximum scale of 50 ℎ −1 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, where L () is the likelihood for a particular magnitude limited sample .The subscript t is used to represent the target values for the clustering, , and number density, .The total length of the data vector is 677, where for each of the 9 magnitude thresholds (  < −22.0, −21.5, ..., −18.0) there is a number density, and 74 correlation function bins between 0.01 <  < 50 ℎ −1 Mpc.In addition, to improve the luminosity function at faint magnitudes, number density constraints are also included for   < −17.5 and   < −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   () and   which produced robust fits with the correct balance between fitting the clustering and the number density.  () and   were the same for all nine of the magnitude fits.We choose to use a constant fractional error for   (), which is 7 %.This value was chosen to avoid overfitting to noise at very small scales.For   , we also use a constant fractional error of 1 %.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.

Best-fitting Meta-parameters
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 5-9 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 ex-ist 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   , which sets the minimum value of  log  at faint magnitudes, and   and   , which describe the exponential increase of  at bright magnitudes.Reference to Equations 6 and 9 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  0 and  0 .This is because  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.
Note that we are only fitting the HODs using diagonal errors.Our aim is to produce mocks that reproduce the galaxy clustering, to be useful for the analyses within DESI, and not to be used to determine the intrinsic errors on the HOD parameters.The errors illustrated by the contours of Figure 4 might change if a full covariance matrix is used.

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.

Comparison to Target Data
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 top panel of Figure 6 shows the predicted two-point clustering for the c000 Planck cosmology compared 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 ratio is shown in the middle panel.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 ∼ 50 ℎ −1 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 ∼ 2 ℎ −1 Mpc, with a dip ∼ 0.4 ℎ −1 Mpc, and this behaviour is inverted for the faintest samples.Below 0.1 ℎ −1 Mpc, the slope of the best-fitting correlation function is different to the target, leading to differences of up to 20% at 0.01 ℎ −1 Mpc.
In the bottom 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.

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 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.
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 Λ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.
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  plotted against  8 for different magnitude threshold samples.A clear negative correlation can be seen, where  is ∼ 1 for the highest values of  8 , but is larger for small  8 .In addition, this correlation becomes steeper for the brightest samples, but with a larger scatter. 8 sets the amplitude of the initial power spectrum at a fiducial scale of 8 ℎ −1 Mpc.Simulations run with higher  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  8 should correlate with lower , if  0 and  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   and   .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 .We note that we have ignored the covariance in the correlation function estimates, but given the explanation above, we do not expect this correlation to strongly affected.
A second example in the right panel of Figure 7, which shows  plotted against the matter density Ω m .Here, a positive correlation is seen, where in cosmologies with higher Ω m , larger values of  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  8 and Ω m .The HODs are the same as in the lower panel of Figure 5.We only show the 0.1   < −22, -21 and -18 samples for clarity, but all samples show the same trends.For the parameter  8 , higher values are associated with fewer satellite galaxies, and have a broader central step function.As already discussed, halos in simulations with high  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 Ω m , the HODs are shifted to higher masses, with a sharper step function for the central galaxies.Even though , the satellite power law slope, as shown in Figure 7 is higher for large Ω 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.

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. (2017Smith et al. ( , 2022b) is used to assign each galaxy a rest frame ( − ) 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.

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  = 0.2 cubic box, assigning -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 ).We assume that the dispersion is constant, but in future, this could be modified to vary radially (equation A24 of Sheth et al. 2001).
The HOD parameters  min and  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 ) beyond this range.
In addition to assigning galaxy luminosities, we also assign rest frame ( − ) 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 of Smith et al. (2017), which itself was based on the method of Skibba & Sheth (2009).Other studies that add colours based on the prescription of Skibba & Sheth (2009) include Skibba et al. (2014); Carretero et al. (2015); Paranjape et al. (2015Paranjape et al. ( , 2021)).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 halos are more likely to be red.
The final cubic box mock is cut to an absolute magnitude limit of 0.1   < −18, which corresponds to a galaxy number density of 3.7 × 10 −2 (ℎ −1 Mpc) −3 for the c000 Planck cosmology simulations.

Cut-sky mocks
To create a mock that is more representative of what DESI will observe, we convert the cubic box mock from the  = 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  = 0.6.Around half of the volume of the final mock fits inside a single copy of the cubic box, and it is above  = 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  = 0.2, which used HODs that were fit to reproduce the target luminosity function at the same redshift  = 0.2.However, the BGS covers a wide redshift range of 0 <  < 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  > 0.15, the target luminosity function is a Schechter function fit to the GAMA luminosity function (Loveday et al. 2012).Below  = 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 * () =  * (0)10 0.4 , where  sets the evolution in number density, and  sets the magnitude evolution (Lin et al. 1999).We use the values  0 = 0.1,  = 1.8, and  = 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  = 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   and the rest-frame 0.1 ( −) colour, which we then convert to the observed quantities.The absolute magnitudes assigned to each galaxy are converted to an apparent magnitude following where   ( obs ) is the luminosity distance from the observer to the galaxy at observed redshift  obs , in units of ℎ −1 Mpc.  ( obs ) is the -band -correction, which accounts for the shift in the band pass with redshift.We -correct to a reference redshift of  ref = 0.1.The distance,   , to each galaxy should only depend on its cosmological redshift,  cos , which does not include the effect of peculiar velocities.However, we use  obs to be consistent with the data, where the cosmological redshift is not known.The exact form of the -correction depends on the filter being used and the type of object being observed.In this work, the -corrections we use come from the GAMA survey and are described by fourthorder polynomials.The -corrections are split into seven different rest frame  −  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 where   ( obs ) is the -band -correction.As with the -band, we use a set of colour-dependent polynomial -corrections.
Finally, an apparent magnitude cut of  < 20.2 is applied, which encompasses both the BGS-BRIGHT and BGS-FAINT samples.

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 ( < 19.5) galaxies with 0.1   < −20 from the c000 Planck cosmology Aba-cusSummit 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  () = ( −  0 ), where  = 0. galaxies brighter than 0.1   = −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   = −20.0, in the redshift range 0.05 <  < 0.25.upper panel of Figure 9, for galaxies with 0.1   < −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  ∼ 0.25, the density decreases since the sample become incomplete due to the  = 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 () 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.
We measure the projected correlation function for a volume limited sample of BGS-BRIGHT galaxies with 0.1   < −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 <  < 0.25, where the sample is complete (and the () is flat).The projected correlation function is defined as where  (  , ) is the correlation function in bins of   and , perpendicular and parallel to the line of sight, respectively.Mocks and data are integrated to  max = 40 ℎ −1 Mpc.The current AbacusSummit mocks are not tuned to reproduce the   (  ) 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  () measurements, but here we show the   (  ) 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 onepercent 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.

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  () 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  () 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 onepercent 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  () measurements.In the real data, this is not available, and we will instead fit to the projected correlation function   (  ).To do this, the HOD fitting method must be modified to compute halo pair counts in bins of   and , and the correlation function  (  , ) can be computed and integrated to obtain   (  ).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 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, assembly bias was detected in the GAMA data in Alam et al. (2024), andPearl et al. (2023) detected a 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 -corrections from the GAMA survey.These colour-dependent -corrections are used to convert the observer frame SDSS -and -bands to the rest frame bands at  ref .
We use these -corrections since they can easily be applied to the mock galaxies where we only have 0.1 ( − ) colours and do not have a full spectrum.However, these -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 -correction can be computed from the spectrum using the fastspecfit code (Moustakas et al. 2023).2Colour-dependent -corrections can then be determined, to create a set of -corrections that are appropriate for the DESI BGS survey.
The cut-sky mocks we have created are constructed from a single simulation snapshot at  = 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).

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 cre-ate 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  ∼ 0.2, consisting of the flux-limited BGS-BRIGHT sample, with -band apparent magnitude limit  = 19.5.The secondary BGS-FAINT sample extends this limit to  = 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 -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  = 0.2, and fitting to real-space correlation functions,  () 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  > 8 ℎ −1 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 metaparameters 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  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 Ω m , the best-fitting HODs in cosmologies with high Ω 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  = 0.2 with galaxies, using the nested set of HODs to assign absolute -band magnitudes following the same method used to create the MXXL mock.A Monte Carlo method is also used to assign 0.1 (−) 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 -band apparent magnitude and observer frame  −  colours are calculated using colour-dependent -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   < −20.The () 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,   (  ) 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   (  ) 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 -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.
Figure1.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.5-9.The upper panel shows the mass parameters  min ,  0 , and  1 on a logarithmic scale, in blue, orange and green, respectively.The lower panel shows the parameters  log in red and  in purple.We plot  − 1 so that both parameters fall within the same y-axis range.

Figure 2 .
Figure 2. Fraction of halos that are subsampled,  , 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  < 10 12 ℎ −1 M ⊙ ; all halos are used above this.At  = 10 11 ℎ −1 M ⊙ , the subsampled fraction is 1%.

Figure 4 .
Figure4.A corner plot from the emcee fitting chain of the 17 meta-parameters describing the smoothly varying HOD curves.This illustrates the correlations between HOD parameters even with our assumption of uncorrelated clustering errors.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).

Figure 5 .
Figure5.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 .
Figure6.Top panel: Best-fitting projected correlation function,  ( ), for the primary 'c000' Planck cosmology (solid curves).Each curve is for a different magnitude threshold sample, as indicated in the legend in the middle panel.The target MXXL clustering is shown by the dashed black curves, which have been rescaled above  = 8 ℎ −1 Mpc to correct for the difference in cosmology.The black shaded region indicates the constant fractional error assumed during the HOD fitting procedure.For clarity, the curves have been offset by 0.1 dex, relative to the 0.1   < −20.0 sample.Middle panel: the ratio of the best-fitting correlation functions to the target MXXL clustering measurements.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.

Figure 7 .
Figure 7. Left panel: best-fitting HOD parameter  vs  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  parameter vs Ω m .Here, a positive correlation is seen.

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

Figure 9 .
Figure 9. Upper panel: the number density of BGS-BRIGHT ( < 19.5) galaxies brighter than 0.1   = −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   = −20.0, in the redshift range 0.05 <  < 0.25.