A Lightcone Catalogue from the Millennium-XXL Simulation

TL;DR

This paper presents a comprehensive lightcone galaxy catalogue from the Millennium-XXL simulation, matching observed galaxy distributions and clustering, useful for future cosmological surveys.

Contribution

It introduces a new halo lightcone catalogue and HOD galaxy catalogue that accurately reproduce observed galaxy properties and clustering across redshifts.

Findings

Galaxy clustering matches SDSS and GAMA measurements

Luminosity function aligns with SDSS and GAMA data

BAO and RSD features are recoverable in the mock catalogue

Abstract

Future galaxy surveys require realistic mock catalogues to understand and quantify systematics in order to make precise cosmological measurements. We present a halo lightcone catalogue and halo occupation distribution (HOD) galaxy catalogue built using the Millennium-XXL (MXXL) simulation. The halo catalogue covers the full sky, extending to z = 2 with a mass resolution of ~1e11 Msun/h . We use this to build a galaxy catalogue, which has an r-band magnitude limit of r < 20.0, with a median redshift of z~0.2. A Monte Carlo HOD method is used to assign galaxies to the halo lightcone catalogue, and we evolve the HODs to reproduce a target luminosity function; by construction, the luminosity function of galaxies in the mock is in agreement with the Sloan Digital Sky Survey (SDSS) at low redshifts and the Galaxy and Mass Assembly (GAMA) survey at high redshifts. A Monte Carlo method is used to assign a 0.1(g-r) colour to each galaxy, and the colour distribution of galaxies at different redshifts agrees with measurements from GAMA. The clustering of galaxies in the mock for galaxies in different magnitude and redshift bins is in good agreement with measurements from SDSS and GAMA, and the colour-dependent clustering is in reasonable agreement. We show that the baryon acoustic oscillation (BAO) can be measured in the mock catalogue, and the redshift space distortions (RSDs) are in agreement with measurements from SDSS, illustrating that this catalogue will be useful for upcoming surveys.

Tables

Table 1. Table 1: Polynomial coefficients of the median k -corrections of galaxies in GAMA in equally spaced bins of ( g − r ) 0.1 {}^{0.1}(g-r) colour, as defined in Eq.  38 . ( g − r ) med 0.1 {}^{0.1}(g-r)_{\rm med} is the median colour in each bin, and A i subscript 𝐴 𝑖 A_{i} are the polynomial coefficients. The constant term A 4 = − 2.5 ​ log 10 ⁡ ( 1.1 ) ≈ − 0.103 subscript 𝐴 4 2.5 subscript 10 1.1 0.103 A_{4}=-2.5\log_{10}(1.1)\approx-0.103 , as described in the text.

${}_{}{}^{0.1}(g-r)_{\mathrm{med}}^{}$	$A_0$	$A_1$	$A_2$	$A_3$
0.131	-45.33	35.28	-6.604	-0.4805
0.298	-20.08	20.14	-4.620	-0.04824
0.443	-10.98	14.36	-3.676	0.3395
0.603	-3.428	9.478	-2.703	0.7646
0.785	6.717	3.250	-1.176	1.113
0.933	16.76	-2.514	0.3513	1.307
1.067	20.30	-4.189	0.5619	1.494

Equations

$$\frac{dN}{dz}={\int }_{M_{\min }(z)}^{M_{\max }}{n}_{\mathrm{unres}}(M,z)\frac{dV}{dz}dM,$$

$$n_{\mathrm{unres}}(M,z)=[1-\mathrm{cut}(M,z)]n_{\mathrm{fit}}(M,z),$$

$$\mathrm{cut}(M,z)=10^{(-z-2)({\log }_{10}(M/h^{-1}{\mathrm{M}}_{\odot })-11{)}^{0.6}}.$$

$$⟨N_{\mathrm{gal}}(>L|M)⟩=⟨N_{\mathrm{cen}}(>L|M)⟩+⟨N_{\mathrm{sat}}(>L|M)⟩.$$

$$⟨N_{\mathrm{cen}}(>L|M)⟩=\frac{1}{2}\left[1+\mathrm{erf}\left(\frac{\log M-\log M_{\min }(L)}{{\sigma }_{\log M}(L)}\right)\right],$$

$$⟨N_{\mathrm{sat}}(>L|M)⟩=⟨N_{\mathrm{cen}}(>L|M)⟩{\left(\frac{M-M_0(L)}{M_1^{\prime }(L)}\right)}^{\alpha (L)},$$

$$L/L_{\ast }=A{\left(\frac{M}{M_t}\right)}^{{\alpha }_M}\exp \left(-\frac{M_t}{M}+1\right),$$

$$\operatorname{spline}(x)=\begin{array}[]{lr}1-6|x|^{2}+6|x|^{3}&|x|\leq 0.5\\ 2(1-|x|)^{3}&0.5<|x|\leq 1\\ 0&|x|>1,\end{array}$$

$$S(x)=\frac{4/3}{\sigma \sqrt{12}}\mathrm{spline}\left(\frac{x-\mu }{\sigma \sqrt{12}}\right).$$

$$⟨N_{\mathrm{cen}}(>L|M)⟩=\frac{1}{2}\left[1+F\left(\frac{\log M-\log M_{\min }(L)}{{\sigma }_{\log M}(L)}\right)\right],$$

$$n_{\mathrm{gal}}(>L,z)=\int n_{\mathrm{halo}}(M,z)⟨N(>L|M,z)⟩dM,$$

$$M_{\mathrm{HOD}}(n_{\mathrm{gal}},z)=f(n_{\mathrm{gal}},z)M_{\mathrm{HOD}}(n_{\mathrm{gal}}),$$

$$n_{\mathrm{gal}}(>L,f(z))-n_{\mathrm{gal}}^{\mathrm{target}}(>L,z)=0.$$

$$\varphi (M)=0.4\ln 10{\varphi }^{\ast }(10^{0.4(M^{\ast }-M)}{)}^{1+\alpha }\exp (-10^{0.4(M^{\ast }-M)}),$$

$$\alpha (z)$$

$$M^{\ast }(z)$$

$${\varphi }^{\ast }(z)$$

$${\varphi }_{\mathrm{SDSS}}(>L)=\int n_{\mathrm{halo}}(M)⟨N(>L|M)⟩dM,$$

$${\varphi }_{\mathrm{target}}(M,z)$$

$$w(z)=(1+e^{-100(z-0.15)}{)}^{-1}.$$

$${\sigma }^2(M)=\frac{G{M}_{200m}}{2R_{200m}},$$

$$⟨g-r|M_r{⟩}_{\mathrm{red}}^{\mathrm{Skibba}}$$

$$\mathrm{rms}(g-r|M_r{)}_{\mathrm{red}}^{\mathrm{Skibba}}$$

$$⟨g-r|M_r{⟩}_{\mathrm{blue}}^{\mathrm{Skibba}}$$

$$\mathrm{rms}(g-r|M_r{)}_{\mathrm{blue}}^{\mathrm{Skibba}}$$

$$f_{\mathrm{blue}}^{\mathrm{Skibba}}(M_r)=0.46+0.07(M_r+20).$$

$$f_{\mathrm{blue}}^{\mathrm{faint}}(M_r)=0.4+0.2(M_r+20),$$

$$f_{\mathrm{blue}}(M_r)=\max \{f_{\mathrm{blue}}^{\mathrm{faint}}(M_r),f_{\mathrm{blue}}^{\mathrm{Skibba}}(M_r)\};$$

$$⟨g-r|M_r{⟩}_{\mathrm{blue}}^{\mathrm{faint}}=0.4-0.03(M_r+16).$$

$$⟨g-r|M_r{⟩}_{\mathrm{sat}}$$

$$f_{\mathrm{sat}}^{\mathrm{blue}}(M_r)=\frac{⟨g-r|M_r{⟩}_{\mathrm{sat}}-⟨g-r|M_r{⟩}_{\mathrm{red}}}{⟨g-r|M_r{⟩}_{\mathrm{blue}}-⟨g-r|M_r{⟩}_{\mathrm{red}}},$$

$$⟨g-r|M_r{⟩}_{\mathrm{sat}}^{\mathrm{Skibba}}=0.83-0.08(M_r+20).$$

Full text

††pubyear: 2016††pagerange: A Lightcone Catalogue from the Millennium-XXL Simulation–A.2

1 Introduction

Upcoming galaxy surveys, such as the Dark Energy Spectroscopic Instrument (DESI) survey (DESI Collaboration et al., 2016a, b) and Euclid (Laureijs et al., 2011), aim to measure the expansion history of the Universe and the growth of cosmic structures. Measurements of galaxy clustering, redshift space distortions and weak lensing will test general relativity, constrain theories of dark energy, and give us precise cosmological constraints.

In order to reach the high precision required to meet these aims, it is necessary to understand and quantify the systematic uncertainties in measurements from surveys, which requires the use of accurate mock catalogues (Baugh, 2008). Since a mock catalogue has a known cosmology and the ‘true’ value of a statistic can be measured directly, they can be used to develop and test the analysis tools which will be used on real observations.

Mocks are also required to test observational strategies and quantify the resultant levels of sample incompleteness. It is often not possible to assign a fibre to every galaxy due to mechanical constraints on fibre positioning (e.g. Hawkins et al., 2003; Guo et al., 2012; Burden et al., 2016; Hahn et al., 2016; Pinol et al., 2016) and even if a fibre is assigned, a redshift measurement can fail if the galaxy has weak emission lines or low surface brightness. This incompleteness may have a significant effect on clustering measurements, and therefore in order to make precise baryon acoustic oscillation (BAO) and redshift space distortion measurements, it is important that this incompleteness is well understood, and that methods are developed and tested in order to mitigate these effects on the measured clustering. The differences in the clustering statistics expected in viable models is small, making it essential that systematics like these are understood.

Mock catalogues which have realistic galaxy clustering can be created from cosmological simulations. In order to see the BAO peak in clustering measurements, at a scale of the order of $100h^{-1}\rm Mpc$, these simulations need to have a very large box size of the order of a few Gpc. Running a hydrodynamical simulation that has both the large volume needed to model such scales, and the resolution to produce faint galaxies down to the flux limit of the survey is infeasible, due to the large computational expense. Dark matter only simulations are much less expensive. There are several schemes which can be used to populate haloes in a dark matter only simulation with galaxies. These include the halo occupation distribution (HOD) (e.g. Peacock & Smith, 2000; Seljak, 2000; Scoccimarro et al., 2001; Berlind & Weinberg, 2002; Kravtsov et al., 2004; Zheng et al., 2005), which describes the probability a halo with mass $M$ contains $N$ galaxies with some property; the closely related conditional luminosity function (CLF) (e.g. Yang et al., 2003), which specifies the luminosity function of galaxies at each halo mass; subhalo abundance matching (SHAM) (e.g. Vale & Ostriker, 2004; Conroy et al., 2006), which assumes a correlation between halo or subhalo properties (e.g. mass or circular velocity), and galaxy properties (e.g. luminosity or stellar mass); and semi-analytic models (SAMs) (e.g. Baugh, 2006; Benson, 2010; Somerville & Davé, 2015), which uses analytic prescriptions to model the formation and evolution of galaxies.

In order to apply a SAM to a simulation, high resolution merger trees are needed, and these are difficult to construct for large volume simulations. However, there are approaches which can augment the resolution of the simulation merger trees (e.g. de la Torre & Peacock, 2013; Angulo et al., 2014; Benson et al., 2016). The SHAM prescription assigns galaxies to subhaloes, requiring a complete subhalo catalogue. Since subhaloes are disrupted when they undergo mergers, this catalogue will only be complete for large subhaloes with thousands of particles, and so a very high resolution simulation is needed to resolve the low mass subhaloes that will be populated by faint galaxies. The HOD, on the other hand, can be applied to a lower resolution simulation, since satellite galaxies can be placed around the central galaxy following an analytic distribution, without knowledge of the subhaloes. The HOD method can also be applied to simulations in which the underlying cosmology has been rescaled (e.g. Angulo & White, 2010).

Ideally, these methods would be used to populate a halo lightcone that is the direct output from a simulation. However, most simulations do not output lightcones, but output snapshots at discrete times. Typically when a HOD method is used, it is applied to a single snapshot. However, this means that the halo bias is constant, and so the clustering of haloes does not evolve with redshift in the mock. Multiple snapshots can be joined together to create a lightcone, but this leads to discontinuities at the boundaries; the same halo could appear twice at either side of the boundary, or not at all (e.g. Fosalba et al., 2015).

The standard abundance matching and HOD schemes do not incorporate evolution. Attempts have been made to extend the abundance matching scheme, such as Moster et al. (2013), which reproduces the observed stellar mass function at different redshifts. There is currently no complete model for HOD evolution, as this evolution would depend on the galaxy sample under consideration. Contreras et al. (2016) use the HODs produced in SAMs to build a simple parametrisation of the evolution of the HOD parameters.

Here, we describe a HOD method which we use to populate haloes over a range of redshifts from the Millennium-XXL (MXXL) simulation with galaxies. We first create a halo lightcone catalogue from the simulation by interpolating the positions of haloes between snapshots, which is then populated with galaxies using HODs, reproducing the observed clustering from the Sloan Digital Sky Survey (SDSS) (Abazajian et al., 2009) and the Galaxy and Mass Assembly (GAMA) survey (Driver et al., 2009, 2011; Liske et al., 2015).

This paper is organised as follows: in Section 2 we describe the MXXL simulation and outline the method for generating the halo lightcone catalogue. In Section 3, we describe the halo occupation distribution model, and our method of evolving the HODs with redshift. In Section 4 we outline the method used to populate the halo lightcone with galaxies, and the method used to assign each galaxy a ${}^{0.1}(g-r)$ colour. In Section 5, we give examples of potential applications of the mock catalogue.

2 Halo Lightcone Catalogue

2.1 The MXXL Simulation

The Millennium-XXL (MXXL) simulation (Angulo et al., 2012b) is a large dark-matter only N-body simulation in the same family as the Millennium simulation (Springel et al., 2005). The volume of MXXL is 216 times larger than Millennium, with a box size of $3h^{-1}\rm Gpc$, and the particle mass is $6.17\times 10^{9}h^{-1}\rm M_{\odot}$, with a force softening of $13.7$ kpc. MXXL adopts a $\Lambda$CDM cosmology with the same WMAP-1 cosmological parameters as the Millennium simulation, $\Omega_{\rm m}=0.25$, $\Omega_{\Lambda}=0.75$, $\sigma_{8}=0.9$, $h=0.73$, and $n=1$ (Spergel et al., 2003). The initial conditions were set at a starting redshift of $z=63$, and the simulation was evolved to $z=0$ with 63 outputs. The large volume of the simulation means that it can be used to study features such as baryon acoustic oscillations (BAOs) and redshift space distortions with good statistics.

2.2 Merger Trees

We use the halo merger trees computed by Angulo et al. (2012b). Haloes were found using a friends-of-friends (FOF) algorithm (Davis et al., 1985), and bound subhaloes were identified using SUBFIND (Springel et al., 2001). Halo merger trees were built by identifying the unique descendant of each subhalo at the subsequent snapshot. For each subhalo, the 15 most bound particles were found, and the subhalo at the next snapshot which contains the greatest number of these particles was defined as the descendant. In the case that two subhaloes contain equal numbers of these particles, the subhalo with the greatest total binding energy was chosen (Angulo et al., 2012a).

2.3 Constructing the Halo Lightcone Catalogue

The full sky halo lightcone catalogue is created using the standard interpolation method (e.g. Merson et al., 2013, but applied to haloes rather than galaxies). An observer is firstly placed randomly inside the MXXL box. If the observer happened to be placed at the centre, haloes at the edge of the box would have a redshift $z\sim 0.5$; multiple periodic replications of the box must therefore be used in order to construct a catalogue that goes to redshifts higher than this. This replication is done without any artificial rotation or translation in order to prevent the introduction of discontinuities. The positions and velocities of each halo at each snapshot are used to interpolate their trajectories through the simulation. From the position and redshift of a halo at two adjacent snapshots, it can be determined whether the halo crossed the observer’s lightcone; if it has, a binary search algorithm is used to find the interpolated position (and velocity) at the redshift where it crosses.

The halo occupation distribution method of creating the galaxy catalogue (Section 3) assigns galaxies to FOF groups. Since the merger tree is defined for SUBFIND subhaloes, we need to infer the merger tree for the FOF haloes. To do this, we make the assumption that the position (and velocity) of the main subhalo (i.e. the most massive subhalo) in each FOF group is the same as that of the FOF group itself. The descendant FOF group is then found from the descendant of the main subhalo. To interpolate the position (and velocity) of each subhalo, we use cubic interpolation (i.e. use a cubic polynomial to describe the path of the halo in each dimension, using the positions and velocities at the previous and next snapshot as boundary conditions).

We use a halo mass definition of $M_{200\rm m}$ (the mass enclosed by a sphere, centred on the halo, in which the average density is 200 times the mean density of the Universe), as stored in the MXXL output for each FOF group. Since the number of galaxies in each halo depends on its mass, $M_{200\rm m}$ must be interpolated between snapshots. Below $z=2$, the simulation snapshots are approximately spaced linearly with expansion factor. We use the descendant of the most massive subhalo to find the descendant of each halo, and then interpolate linearly in mass between snapshots, finding the mass at the redshift at which it crosses the lightcone. In the case that two or more haloes merge between snapshots, the total mass of the haloes is interpolated linearly, and each halo is assigned a constant fraction of the total mass. If the halo is not the most massive progenitor, a random time between snapshots is chosen for the merger to take place. If the halo crosses the observer’s lightcone after this time, the merger has happened, and the interpolated mass of the halo is transferred to the most massive progenitor.

The mass function of the halo lightcone at low redshifts $z<0.1$ is shown in Fig. 1 and compared to the Sheth & Tormen (1999) and Jenkins et al. (2001) analytic mass functions. At high masses, the mass function of the lightcone catalogue is in reasonable agreement with Sheth & Tormen (1999) and Jenkins et al. (2001), but there is a lower abundance of less massive haloes. This difference is because haloes in the simulation are identified using a FOF algorithm, and not a spherical overdensity (SO) finder. However, a SO mass is calculated for each FOF halo ($M_{\rm 200m}$). Any small overdensities close to a large FOF group would be identified as part of the large FOF group, and therefore these would be missing from the halo catalogue.

In order to add haloes to the catalogue that are below the MXXL mass resolution (Section 2.5), and to evolve the HODs with redshift (Section 3.2), it is useful to have a smooth function which is in close agreement with the actual mass function of the catalogue. For this, we take a mass function with the same form as Sheth & Tormen (1999), but refit the parameters to the MXXL mass function. This fit is shown as the green curve in Fig. 1, which is in better agreement with the MXXL mass function at low masses, and is close to the fit given in equation 2 of Angulo et al. (2012b). The MXXL mass function peels away from this fit slightly at masses close to the resolution limit, but is complete for masses greater than $\sim 10^{12}h^{-1}\rm M_{\odot}$.

The number density of haloes as a function of redshift in the halo lightcone catalogue is shown in Fig. 2 for several mass thresholds. If halo masses are kept fixed between snapshots, step features can be seen, since the mass function is being kept frozen and will only change at the next snapshot. These features are most apparent for the highest mass threshold, for which the number density decreases more rapidly at high redshifts. Mass interpolation greatly reduces these features.

The large scale real space correlation function of the lightcone catalogue for FOF groups with masses $M_{200\rm m}>3\times 10^{12}h^{-1}\rm M_{\odot}$ and $z<0.5$ is shown in Fig. 3. This redshift limit avoids structures being repeated due to periodic replication of the box111Lightcones with a wide opening angle, or directed along the principle axes of the simulation, that extend beyond $z=0.5$ will contain repeated structures. This will result in clustering errors being underestimated.. The BAO peak can be seen clearly in the clustering of haloes.

2.4 Caveats

The interpolation scheme uses the position and velocity of a halo at two snapshots as boundary conditions in order to find the path the halo moved through in the simulation. If two haloes merge together, there is not enough information to determine when this occurs, so we assume they merge at a random time. If a new halo forms, or drops below the resolution limit, we assume this happens exactly on a snapshot.

To construct the merger trees, the descendant of a subhalo is defined as the subhalo which contains the majority of its 15 most bound particles (Angulo et al., 2012a). However, it is likely that some of the particles of the descendant subhalo were not in its progenitor, and vice versa. All of these particles are used to calculate the position and velocity of the subhalo, which can occasionally lead to jumps in the position of a subhalo that are inconsistent with its velocity.

Sometimes, a halo can be lost by the halo finder at one snapshot, but is then found again at a later snapshot. This can happen if a small halo passes very close to a more massive halo at one snapshot; the SUBFIND algorithm can fail to identify the small halo as the algorithm finds that its particles are bound to the massive halo. The MXXL merger trees we use do not make any attempt to add in these haloes lost by SUBFIND. However, since we use the most massive subhalo in a FOF group to trace the FOF merger trees, and use $M_{\rm 200m}$ as the mass definition, our results should not be affected much by small subhaloes being lost by the halo finder.

2.5 Haloes below the mass resolution

Populating the resolved haloes in the MXXL halo lightcone with galaxies will result in incompleteness in a magnitude limited galaxy catalogue at low redshifts. This is because intrinsically faint galaxies which are sufficiently close to the observer to be bright enough to be included in the catalogue occupy haloes which fall below the MXXL mass resolution. We use our fit to the MXXL halo mass function in order to add these haloes into the lightcone catalogue, and position them randomly in the catalogue so that they are unclustered. Other methods for augmenting the halo catalogue exist (e.g. de la Torre & Peacock, 2013; Angulo et al., 2014; Benson et al., 2016), but we find that this simple method is able to bring the $dN/dz$ of galaxies in the catalogue into better agreement with the measured $dN/dz$ from GAMA, while only having a very small effect on the measured clustering.

The redshift distribution of the haloes which need to be added to the lightcone catalogue can be calculated from the integral,

[TABLE]

where $n_{\rm unres}(M,z)$ is the number density of unresolved haloes, $dV/dz$ is the comoving volume per unit redshift, $M_{\rm min}(z)$ is the minimum halo mass that could host a galaxy brighter than the faintest observable galaxy in the survey at that redshift, and $M_{\rm max}=10^{12}h^{-1}M_{\odot}$ is the mass at which the MXXL mass function is judged to be complete. If the survey is flux limited, then the faintest observable galaxy at each redshift is set by an apparent magnitude threshold; for our mock catalogue we set this threshold to $r=20$, as this is the magnitude threshold for the DESI Bright Galaxy Survey (BGS) (DESI Collaboration et al., 2016a). The number density of unresolved haloes is given by $n_{\rm unres}(M,z)=n_{\rm fit}(M,z)-n_{\rm res}(M,z)$, where $n_{\rm fit}(M,z)$ is our fit to the number density of haloes in the lightcone, extrapolated to low masses, and $n_{\rm res}(M,z)$ is the number density of haloes resolved in MXXL. We model the mass function of resolved haloes by multiplying the fit to the mass function by a cutoff at the mass resolution limit of $M_{200\rm m}\sim 10^{11}h^{-1}{\rm M_{\odot}}$,

[TABLE]

where a good approximation to the cutoff is given by

[TABLE]

In order to add unresolved haloes to the catalogue, we first randomly draw a redshift for each unresolved halo from the $dN/dz$ distribution defined in Eq. 1. The mass of each halo is then randomly assigned using the mass distribution at the redshift of the halo defined by $n_{\rm unres}(M,z)$. The haloes are then randomly positioned uniformly on the sky. Since the unresolved haloes are randomly positioned so that they are unclustered, redshift space distortions do not affect their clustering, and so we set the velocity of each of these haloes to zero. A random concentration is also assigned from the mass-concentration relation of MXXL (with scatter), extrapolated to lower masses.

While the introduction of unclustered haloes only has a small effect on the two-point correlation function, other statistics might also change, for example, density estimators and void statistics. We have not checked the size of this effect, but the final galaxy catalogue includes a flag which indicates whether a galaxy lives inside one of these haloes, enabling these galaxies to be removed when calculating other statistics.

3 Halo Occupation Distribution

Galaxies are biased tracers of the underlying dark matter density field. The halo occupation distribution (HOD) describes this bias using the probability that a halo of mass $M$ contains $N$ galaxies with a certain property, $P(N|M)$, providing a physical interpretation of galaxy clustering measurements.

The mean number of galaxies in a halo of mass $M$ which are brighter than some luminosity threshold, $L$, can be written as a sum of central and satellite galaxies (e.g. Zheng et al., 2005),

[TABLE]

We use central and satellite occupation functions of the same form as Zehavi et al. (2011). The mean number of central galaxies brighter than $L$ is described by a smoothed step function,

[TABLE]

where $\operatorname{erf}(x)=2\pi^{-1/2}\int^{x}_{0}e^{-x^{\prime 2}}dx^{\prime}$ is the error function. The parameter $M_{\rm min}$ is the halo mass for which half of haloes contain a galaxy brighter than $L$, and $\sigma_{\log M}$ sets the width of the step. For $M\gg M_{\rm min}$, $\langle N_{\rm cen}(>L|M)\rangle=1$, while for $M\ll M_{\rm min}$, $\langle N_{\rm cen}(>L|M)\rangle=0$. The mean number of satellites per halo brighter than $L$ is given by a power law,

[TABLE]

where $M_{0}$ is the cutoff mass scale, $M^{\prime}_{1}$ the normalisation, and $\alpha$ the power law slope. $M^{\prime}_{1}$ is different to $M_{1}$, the mass of a halo which on average contains 1 satellite, although the two quantities are related222Since the satellite occupation function is modified by the centrals, the relation $M_{1}=M^{\prime}_{1}+M_{0}$ is not exact.. The power law is also multiplied by the central occupation function, which ensures that the brightest galaxy in the halo is the central; there cannot be a satellite brighter than $L$ without there first being a central galaxy brighter than $L$.

3.1 HODs at redshift 0.1

We use HOD parameters calculated from the SDSS using the procedure of Zehavi et al. (2011). These are calculated for different luminosity threshold galaxy samples, using an MCMC code to find the best fitting HOD parameters which reproduce the measured projected correlation functions to within the SDSS uncertainties. Since the cosmology of the MXXL simulation is different to that used by Zehavi et al. (2011), the parameter fitting was redone using the Millennium cosmology. The SDSS HODs use the absolute r-band magnitude of each galaxy, k-corrected to a reference redshift of $z_{\rm ref}=0.1$ (see Section 4.3), which is the median redshift of the survey. We denote absolute magnitudes k-corrected to this redshift as ${}^{0.1}\!M_{r}$. Absolute magnitudes written as ${}^{0.1}\!M_{r}$ assume $h=1$.

The best fitting HOD parameters, in Millennium cosmology, are shown by the points in Fig 4. We do not show the errors as they are misleading, due to the probability distributions being highly asymmetric. Projecting these asymmetric probability distributions to 1 dimensional errors can lead to the best fitting values of some of the parameters being outside the error bars, particularly for $\sigma_{\log M}$.

For each HOD parameter, a least squares routine is used to fit a function which describes the variation with luminosity, which are shown by the dashed lines in Fig. 4. The top panel shows $M_{\rm min}(L)$ and $M^{\prime}_{1}(L)$, for which we fit curves of the same functional form as Eq. 11 from Zehavi et al. (2011),

[TABLE]

where $A$, $M_{t}$ and $\alpha_{M}$ are free parameters. We fit a power law to $M_{0}(L)$ (second panel). This is a poor fit for the points at ${}^{0.1}M_{r}=-18$, which is over 3 orders of magnitude lower than the fit, and ${}^{0.1}M_{r}=-19$, which is 20 orders of magnitude lower. However, increasing the value of this parameter by many orders of magnitude has a very small effect on the shape of the HODs. This is because the occupation function of central galaxies adds a second cutoff to Eq. 6; if $M_{0}$ is below this cutoff, it will not affect the shape of the HODs. The parameter $\alpha(L)\sim 1$ at low luminosities, but increases for the highest luminosity samples (third panel). We fit a linear relation, which smoothly transitions to $\alpha=1$ at low luminosities. $\sigma_{\log M}(L)$ (bottom panel) is small at low luminosities, with a step up to $\sim 0.7$ for the brightest two samples. We fit a sigmoid function to $\sigma_{\log M}$, where the width of the step is set such that the HODs do not overlap.

The large step in $\sigma_{\log M}(L)$ means that, as the luminosity threshold is increased, there is a rapid jump in the amount of scatter in the luminosities of central galaxies. This results in overlapping HODs, as can be seen for the ${}^{0.1}M_{r}<-21$ and ${}^{0.1}M_{r}<-21.5$ samples in Fig. 5. For two luminosity thresholds $L_{1}$ and $L_{2}$, where $L_{1}<L_{2}$, it must be true that $\langle N_{\rm gal}(>L_{1}|M)\rangle\geq\langle N_{\rm gal}(>L_{2}|M)\rangle$ since all galaxies brighter than $L_{2}$ are also brighter than $L_{1}$. However, if the two occupation functions cross, then this condition is not satisfied for haloes below the mass at which they cross. We therefore must model the HODs such that there is no overlap. There exist HOD frameworks in which the occupation functions cannot overlap (see e.g. Leauthaud et al., 2011), but since we are using HOD parameters obtained using the standard HOD framework, we make a small modification to these HODs to prevent any crossing, as set out below.

Eq. 5 assumes that the scatter set by the parameter $\sigma_{\log M}(L)$ is Gaussian. Since a Gaussian function has a long tail which extends to infinity, there will always be an overlap between HODs if $\sigma_{\log M}(L_{2})>\sigma_{\log M}(L_{1})$. We instead approximate the Gaussian by using a spline kernel (Schoenberg, 1946),

[TABLE]

which has $\operatorname{spline}(0)=1$, mean = 0, variance = 1/12, and $\operatorname{spline}(x)=0$ for $|x|>1$. This function can be rescaled and normalised to approximate any Gaussian of mean $\mu$ and variance $\sigma^{2}$ as,

[TABLE]

The HOD for central galaxies can therefore be written as

[TABLE]

where $F(x)=2\int^{x}_{0}S(x^{\prime})dx^{\prime}$. The best fitting values of $\sigma_{\log M}$ (shown by the points in the bottom panel of Fig. 4), suggest a sharp step between ${}^{0.1}M_{r}=-21$ and ${}^{0.1}M_{r}=-21.5$. Even using Eq. 10, the HODs will still overlap with this abrupt step, but they will not overlap if the step is gradual, unlike Eq. 5. We make the step in $\sigma_{\log M}(L)$ as narrow as we can while preventing the HODs crossing (shown by the dashed curve).

The HODs using the SDSS HOD parameters, and our fits, are shown in Fig. 5. Our fits produce halo occupation functions which are in reasonable agreement with the SDSS HODs, with the exception of the ${}^{0.1}M_{r}<-21$ and ${}^{0.1}M_{r}<-21.5$ samples, where the width of the step set by the parameter $\sigma_{\log M}$ is too broad and narrow respectively. This is necessary to prevent the HODs from overlapping. The SDSS HOD for the ${}^{0.1}M_{r}<-21$ sample appears to have a sharp transition from centrals to satellites, which is due to a large value of $M_{0}$ compared to $M^{\prime}_{1}$.

3.2 Redshift Evolution

In order to evolve the occupation functions with redshift, we first choose a target luminosity function, $\phi_{\rm target}(L,z)$, that we would like the galaxies in the mock catalogue to reproduce as a function of redshift. This luminosity function defines a mapping between a luminosity threshold $L$ at redshift $z$, and the number density of galaxies brighter than this, $n_{\rm gal}^{\rm target}(>L,z)$.

For a given HOD, the number density of galaxies brighter than $L$ can be calculated from the integral

[TABLE]

where $n_{\rm halo}(M,z)$ is the number density of haloes of mass $M$ at redshift $z$, and $\langle N(>L|M,z)\rangle$ is the halo occupation function at redshift $z$. The HODs must evolve with redshift such that the condition $n_{\rm gal}(>L,z)=n_{\rm gal}^{\rm target}(>L,z)$ is satisfied.

Since the target luminosity function defines a mapping between a luminosity threshold and the number density of galaxies, the occupation functions can be rewritten as a function of number density, $n_{\rm gal}$: $\langle N(>L|M,z)\rangle\equiv\langle N(n_{\rm gal}|M,z)\rangle$. The shape of the HOD could evolve in a complex way, but for simplicity we keep the shape of the occupation function fixed for constant galaxy number density, but slide the HODs along the halo mass axis such that the target luminosity function is achieved. That is, the HOD parameters $\sigma_{\log M}(n_{\rm gal},z)$ and $\alpha(n_{\rm gal},z)$ are kept constant, but the 3 mass parameters $M_{\rm min}$, $M_{0}$ and $M_{1}$ are all multiplied by some factor $f$,

[TABLE]

where $M_{\rm HOD}$ is one of the HOD mass parameters. The value of $f$ required to achieve the target luminosity function is found by finding the root of the equation

[TABLE]

At high redshifts, the target luminosity function we use is the evolving Schechter function fit to the luminosity function estimated from the Galaxy and Mass Assembly (GAMA) survey. The Schechter function can be written in terms of magnitudes as

[TABLE]

where $\phi^{*}$ is the normalisation, $M^{*}$ is a characteristic magnitude and $\alpha$ is the faint end slope. For GAMA, Loveday et al. (2012, 2015) model the evolution of the Schechter parameters with redshift as

[TABLE]

where $Q$ parametrises the evolution in luminosity, $P$ parametrises the evolution in number density, and $z_{0}=0.1$ is the same reference redshift as used for the k-corrections (Section 4.3). The faint end slope is kept constant with redshift since there is not enough data to constrain it at high redshifts. We use the evolving Schechter function, $\phi_{\rm GAMA}$, from Loveday et al. (2012) with $P=1.8$ and $Q=0.7$.

However, the shape of the GAMA Schechter luminosity function is slightly different than the SDSS luminosity function. Using it as the target at all redshifts would result in the evolution parameter $f\neq 1$ at $z=0.1$, meaning that the HODs would change from the HODs measured from SDSS. In order to not change the HODs at $z=0.1$, we use the luminosity function from SDSS, $\phi_{\rm SDSS}$, as the target at low redshifts. The SDSS target luminosity function we use is the result of the integral

[TABLE]

where $n_{\rm halo}(M)$ is the number density of haloes at $z=0.1$, and $\langle N(>L|M)\rangle$ is the (unevolved) occupation function. The result of this integral is close to the Blanton et al. (2003) luminosity function for absolute magnitudes brighter than ${}^{0.1}M_{r}=-19$, and by definition $f=1$, so the HODs remain unchanged from SDSS at this redshift. However, at magnitudes fainter than ${}^{0.1}M_{r}=-19$, the result of this integral is very flat, while the Blanton et al. (2003) SDSS luminosity function is steeper; at ${}^{0.1}M_{r}=-17$ they differ by a factor of $\sim 2$. We therefore smoothly transition to the Blanton et al. (2003) luminosity function at ${}^{0.1}M_{r}=-19$. This is then extrapolated to fainter magnitudes with a power law.

We interpolate the target luminosity function from $\phi_{\rm SDSS}$ at low redshifts to $\phi_{\rm GAMA}$ at high redshifts,

[TABLE]

where the transition between $0.1<z<0.2$ is set by the sigmoid function

[TABLE]

The evolution parameter, $f$, for this target luminosity function, is shown in Fig. 6 as a function of magnitude for different redshifts. At $z=0.1$, $f$ is close to 1, by definition. However, it is not exactly 1 because the function $w(z)$, which sets the transition between the two target luminosity functions is close to, but not exactly 0 at $z=0.1$. At $z=0.1$, $f$ is equal to 1 to within 1%. Fainter than magnitude $-19$, $f(z=0.1)<1$. At these faint magnitudes, the target luminosity function is transitioning to the Blanton et al. (2003) luminosity function. Keeping $f(z=0.1)=1$ at all magnitudes produces a luminosity function which, while being close to SDSS at the bright end, is too flat at the faint end, so this transition is required to bring the luminosity function of the mock into better agreement with the data.

The evolution of the parameter $M_{1}$ implied by this evolution of $f$ is shown in Fig. 7 for galaxy samples of a fixed number density, up to $z=0.6$. Since the shape of the HODs are kept fixed for a fixed number density, but the HODs are evolved along the mass axis, the other mass parameters $M_{\rm min}$ and $M_{0}$ show the same evolution, while $\sigma_{\log M}$ and $\alpha$ are held constant. For comparison, we also show the evolution reported in Contreras et al. (2016) from their fit to the evolution found in the Gonzalez-Perez et al. (2014) version of the GALFORM semi-analytic galaxy formation model (Cole et al., 2000). We find that $M_{1}$ decreases slightly with redshift, in remarkably close agreement with what is found in Contreras et al. (2016), although the highest number density samples show slightly less evolution. By construction, the ratio of the parameters $M_{1}/M_{\rm min}$ is kept constant in the mock. This is in contrast to the behaviour found in Contreras et al. (2016), where they reported that this ratio decreases over the same redshift range. The amount by which this ratio decreases depends on the semi analytic model used, and on the number density of galaxies; at most it decreases by $\sim 50\%$. The evolution of the HOD model could be extended to include this change in $M_{1}/M_{\rm min}$, but we find that simply keeping the mass ratio fixed produces a good match to the measured clustering (see Fig. 12).

4 Mock Galaxy Catalogue

Now we describe in detail the HOD method used to populate the halo catalogue with galaxies, assign each galaxy a luminosity and ${}^{0.1}(g-r)$ colour and compare the resultant clustering in the mock with measurements from SDSS and GAMA. These details can be skipped by the reader, but we give a brief summary below.

Section 4.1 describes the HOD method for populating the halo lightcone with galaxies with luminosities. This Monte Carlo method is based on Skibba et al. (2006), but extended to an evolving 5 parameter HOD. The number and luminosity of galaxies in each halo are randomly generated such that the input HODs are reproduced. Central galaxies are assigned the position and velocity of the halo, and satellites are randomly positioned around the central, following an NFW density profile, and assigned a random virial velocity.

The method for assigning a ${}^{0.1}(g-r)$ colour to each galaxy is described in Section 4.2. This is based on Skibba & Sheth (2009), and randomly assigns a colour from a parametrisation of the SDSS colour magnitude diagram. Section 4.2 describes our modification to the parametrisation of the colour magnitude diagram, which includes evolution, and is in agreement with measurements from GAMA. The colour assigned to each galaxy depends only on its luminosity, its redshift, and whether it is a central or satellite galaxy; there is no explicit dependence on halo mass. Colour-dependent k-corrections derived from GAMA are described in Section 4.3. The colour-dependent clustering of galaxies in the mock is shown in Section 4.4.

4.1 Constructing the Galaxy Catalogue

We use a modified version of the method of Skibba et al. (2006) to populate the halo lightcone catalogue with galaxies, and to assign each galaxy an $r$-band absolute magnitude, k-corrected to $z=0.1$. Skibba et al. (2006) use a 3 parameter HOD in which the occupation function of central galaxies is simply a step function; we have extended this method in order to reproduce the 5 parameter HOD given by Eq. 6 & 10, which adds scatter to the luminosity of central galaxies, as required by the SDSS clustering data. We also use the fits to the HOD parameters as a function of luminosity as described in Section 3.1 and shown in Fig. 4. To be consistent with the mass definition used in Zehavi et al. (2011), we take the halo mass to be $M_{200\rm m}$; i.e. the mass enclosed by a sphere in which the average density is 200 times the mean density of the Universe.

For each halo, a number, $x$, is randomly drawn from the spline kernel probability distribution, $S(x)$ (Eq. 9), with $\mu=0$ and $\sigma=1$. This introduces the scatter in the luminosity of the central galaxy, relative to the average luminosity in a halo of this mass. The luminosity $L$ which is required to produce this scatter is found by solving $x\sigma_{\log M}(L)/\sqrt{2}=\log M-\log M_{\rm min}(L)$, where the factor of $\sqrt{2}$ comes from how $\sigma_{\log M}$ is defined. Finally the central galaxy is positioned at the centre of the halo, with the same velocity.

To populate a halo with satellite galaxies, a minimum luminosity, $L_{\min}$, must first be chosen. We vary $L_{\min}$ with redshift, choosing it to be slightly fainter than the luminosity corresponding to $r=20$. This ensures that the final mock catalogue is complete to $r=20$ at all redshifts, while preventing galaxies that are too faint to be observed being unnecessarily added to the catalogue. The number of satellite galaxies to be added to each halo is drawn from a Poisson distribution with mean $\langle N_{\rm sat}(>L_{\rm min}|M)\rangle$, which is given by Eq. 6. For each satellite, a uniform random number $0<u<1$ is drawn, and the luminosity is found such that $\langle N_{\rm sat}(>L|M)\rangle/\langle N_{\rm sat}(>L_{\rm min}|M)\rangle=u$. The satellite galaxies are assigned a random virial velocity, relative to the velocity of the central galaxy, which is drawn from a Maxwell-Boltzmann distribution with a line of sight velocity dispersion

[TABLE]

where $R_{\rm 200m}$ is the radius of the sphere, centred on the halo, in which the enclosed density is 200 times the mean density of the Universe. Finally, the satellite galaxies are positioned randomly around the centre of the halo such that they follow an NFW (Navarro et al., 1997) density profile, which is truncated at $R_{200\rm m}$. We find that using the same concentration, $c$, as the halo, calculated from $c=2.16R_{200\rm m}/R_{\rm Vmax}$, where $R_{\rm Vmax}$ is the radius at which the maximum circular velocity occurs, produces angular clustering which is too strong at small angular scales compared to SDSS (Wang et al., 2013). This can be improved by reducing the concentration of all haloes by a factor of 2 (see Fig. 8). We therefore use these reduced concentrations when positioning satellite galaxies inside each halo.

The HODs from Zehavi et al. (2011) were fit to the projected correlation functions, using the mass-concentration relation of Bullock et al. (2001), modified to be consistent with their mass definition. This mass-concentration relation is close to what is seen in MXXL. However, modifying the concentrations only has a small effect on the 1-halo term of the projected correlation functions. Down to separations of $0.1h^{-1}\rm Mpc$, the clustering in the mock catalogue only changes by a small amount. It is the change in the clustering at physical scales smaller than this which causes the small scale angular clustering to improve, and this is below the scale at which the projected correlation functions were measured in SDSS. The angular correlation function, $\omega(\theta)$, in the mock catalogue is shown in Fig. 8 for galaxies in bins of apparent magnitude, compared to the angular clustering measured in SDSS (Wang et al., 2013). Solid lines show the clustering in the mock with concentrations reduced by a factor of 2, which is in good agreement with SDSS down to a small angular separation of 20 arcsec, although for the faintest sample the clustering is a little low. Using unmodified concentrations results in $\omega(\theta)$ having a slope which is steeper than the SDSS measurements, shown by the dashed curves, resulting in clustering which is too strong at small angular scales. The introduction of unclustered haloes below the mass resolution has the effect of reducing the clustering in the mock, but as we show later in Section 4.1.3, this effect is small.

4.1.1 The luminosity function of the mock

The Petrosian r-band luminosity function of the galaxy catalogue is shown in Fig. 9 for galaxies in three redshift bins. The dashed lines show the target luminosity at the median redshift of each bin, showing that this evolving target luminosity function is reproduced in the mock catalogue. The smaller panel in Fig. 9 compares the luminosity function in the mock at low redshifts with the Blanton et al. (2003) luminosity function from SDSS. Brighter than ${}^{0.1}M_{r}=-19$, the luminosity function in the mock is in good agreement with SDSS, which indirectly shows that the mass function of the MXXL lightcone is close to the Jenkins et al. (2001) mass function assumed by Zehavi et al. (2011), and our fits to the HOD parameters as a function of luminosity are a good approximation to the actual values. Fainter than ${}^{0.1}M_{r}=-19$, the luminosity functions agree by construction.

4.1.2 The redshift distribution of the mock

The redshift distribution of galaxies brighter than an apparent magnitude limit of $r=19.8$ is shown in Fig. 10 (see Section 4.3 for the k-corrections used), and compared to the GAMA survey. The $dN/dz$ of the mock catalogue is in good agreement with GAMA, within 15% of the fitted curve at most redshifts. Without adding in the low mass, unresolved haloes at low redshifts (see Section 2.5), there is a deficit in the $dN/dz$ for $z\lesssim 0.1$ (dashed red curve); adding in these haloes increases the number of low redshift haloes, bringing the $dN/dz$ into better agreement with GAMA (solid red curve).

4.1.3 Clustering of the mock

Projected correlation functions of galaxies in the mock catalogue are shown by the solid curves in Fig. 11 for different luminosity threshold samples at $z\sim 0.1$, where we have calculated the two point correlation functions using the publicly available code CUTE (Alonso, 2012)333http://members.ift.uam-csic.es/dmonge/CUTE.html. These are compared to the measured clustering from SDSS (points with error bars), and the clustering predicted by the best fitting HODs (dashed lines). We use the same redshift ranges as the SDSS volume limited luminosity threshold samples (see table 2 in Zehavi et al., 2011). To be consistent with the definition of magnitude used in Zehavi et al. (2011), magnitudes are evolved to $z=0.1$ using the evolution model $E(z)=Q_{0}(1+Q_{1}(z-z_{0}))(z-z_{0})$, where $Q_{0}=2$, $Q_{1}=-1$ and $z_{0}=0.1$. The clustering in our galaxy catalogue is in reasonable agreement with the projected correlation functions measured from SDSS.

The small differences in the large scale clustering between the mock catalogue and the clustering predicted by the best fitting HODs can be understood by comparing the HODs in Fig. 5. For example, the ${}^{0.1}M_{r}<-19$ sample is slightly less clustered in the mock. The fit to the HOD has a smaller $M_{\rm min}$ than the best fitting SDSS HOD, meaning that this sample contains more low mass haloes. These haloes are less biased, and therefore the clustering is reduced compared to SDSS. The ${}^{0.1}M_{r}<-22$ sample contains more high mass haloes, and should therefore be more clustered than SDSS, but the opposite is seen. This is because the brightest samples cover a wider redshift range, and are affected more by the evolution of the HODs.

The clustering of galaxies in the mock catalogue is also affected by the introduction of haloes below the MXXL mass resolution, which are unclustered. Adding these haloes will therefore have the effect of reducing the measured galaxy clustering. The galaxies which reside in these haloes are faint, and have low redshifts, and so the faintest galaxy samples in Fig. 5 are affected by this more than the bright samples. For the ${}^{0.1}M_{r}<-18.5$ sample, we illustrate the size of this effect: the magenta dashed curve shows the projected correlation function with galaxies residing in unresolved haloes omitted. Including these galaxies reduces the clustering, but only by a very small amount.

We have checked that if we modify our fits to the HOD parameters to agree exactly with the best fitting SDSS parameters at one magnitude, and do not evolve the HODs, we reproduce the SDSS correlation functions very closely for that magnitude limit.

In Fig. 12 we show the projected correlation functions in the mock catalogue at high redshifts, compared to the clustering measured in GAMA by Farrow et al. (2015), which has a high completeness of galaxy pairs (Robotham et al., 2010). Here, the magnitude ranges are defined for magnitudes k-corrected to a reference redshift of $z=0$, denoted as ${}^{0.0}M_{r}$, which have also had evolutionary corrections applied. To be consistent with Farrow et al. (2015), we use the same evolutionary correction $E(z)=-Q(z-z_{\rm ref})=-1.45z$ for this comparison. We find that the clustering at high redshifts is in good agreement with the clustering seen in GAMA. The agreement is least good in the lowest redshift bin, but at low redshifts we find good agreement with SDSS, which covers a much larger area of the sky than GAMA.

4.2 Assigning Colours

We use the method of Skibba & Sheth (2009) to assign each galaxy a ${}^{0.1}(g-r)$ colour, where $g$ and $r$ are SDSS DR7 model magnitudes (Abazajian et al., 2009; Baldry et al., 2010). This method parametrises the red and blue sequence of the colour-magnitude diagram as two Gaussians with a mean and rms that are linear functions of magnitude. A galaxy is randomly chosen to be red or blue, then a colour is drawn from the appropriate Gaussian. We have modified the parametrisation of the colour-magnitude diagram given in Skibba & Sheth (2009) to bring the faint end into agreement with the colour-magnitude diagram from GAMA, and to add evolution. For those interested, this is described in detail below. For clarity, the 0.1 superscript has been omitted from the equations in the following subsections.

4.2.1 Low redshift

For redshifts $z<0.1$, we use the parametrisation used by Skibba & Sheth (2009) that produces a good approximation to the SDSS colour-magnitude diagram, which we summarise below. However, we make some slight modifications to the parametrisation to also bring it into agreement with GAMA at faint magnitudes.

The mean and rms of the red and blue sequences are given by

[TABLE]

and

[TABLE]

The total fraction of galaxies which are blue is also parametrised as a linear function of magnitude, given by

[TABLE]

These relations from Skibba & Sheth (2009) produce a colour-magnitude diagram which is in good agreement with SDSS at the bright end. However, the faint end does not agree with what is seen in GAMA (e.g. the first panel in figure 6 of Loveday et al., 2012). Firstly, at ${}^{0.1}M_{r}=-16$, all the galaxies should lie on the blue sequence, while the fraction of galaxies which are blue given by Eq. 22 is 0.74. At faint magnitudes, we instead use a blue fraction given by

[TABLE]

so the total fraction of blue galaxies is

[TABLE]

$f_{\rm blue}(M_{r})$ is capped so it is always in the range $0\leq f_{\rm blue}(M_{r})\leq 1$. Another issue with the parametrisation of Skibba & Sheth (2009) is that faint galaxies which lie on the blue sequence are too blue in comparison to the galaxies in GAMA. At ${}^{0.1}M_{r}=-18.7$, we transition to a flatter blue sequence, given by

[TABLE]

If the fraction of satellite galaxies that are blue, $f_{\rm sat}^{\rm blue}(M_{r})$, is specified, then the mean colour of satellite galaxies is given by

[TABLE]

Conversely, Eq. 4.2.1 can be rearranged, and the mean satellite colour can be used to specify the fraction of satellites that are blue,

[TABLE]

(equation 8 from Skibba & Sheth, 2009, but for blue galaxies). The average colour of a satellite galaxy is parametrised by Skibba & Sheth (2009) as

[TABLE]

Modifying the mean satellite colour has the effect of changing the strength of the colour dependent clustering. We find that we get a better agreement with the clustering in SDSS by modifying the mean satellite colour to

[TABLE]

At, for example, ${}^{0.1}M_{r}=-16$, the fraction of blue satellites given by Eq. 27 & 29 is less than 1, meaning that some satellite galaxies are red. However, all galaxies at this magnitude should lie on the blue sequence, as determined from Eq 24. In order to achieve the correct $f_{\rm blue}$ from Eq 24, the fraction of satellites which are blue must be $(f_{\rm blue}-f_{\rm cen})/f_{\rm sat}$ if all central galaxies are on the blue sequence. If the value of $f_{\rm sat}^{\rm blue}$ is greater than this, it is still possible to get the correct $f_{\rm blue}$ by making central galaxies red, but $f_{\rm sat}^{\rm blue}$ cannot be less than this. To ensure that at faint magnitudes we get the total fraction of blue galaxies given by Eq 24, we take the fraction of blue satellites to be

[TABLE]

The fraction of central galaxies that are blue can then be determined from $f_{\rm blue}(M_{r})$ and $f_{\rm sat}^{\rm blue}(M_{r})$. However, Skibba & Sheth (2009) erroneously state that the fraction of central galaxies which are blue is

[TABLE]

where $f_{\rm cen}(M_{r})$ is the fraction of galaxies which are centrals. Eq. 31 is only true if all satellite galaxies are red; since a significant fraction of faint satellites are blue, the fraction of blue central galaxies needs to be reduced to ensure we get the correct total fraction of blue galaxies given by Eq. 22. This is achieved by changing Eq. 31 to

[TABLE]

For each galaxy, a uniform random number $x$ is drawn in the interval $0<x<1$. For central galaxies, if $x<f_{\rm cen}^{\rm blue}(M_{r})$ (given by Eq. 32), the galaxy is blue, and a colour is drawn randomly from the Gaussian distribution defined by Eq. 4.2.1, otherwise it is red, and the colour is drawn from Eq. 4.2.1. Similarly, satellite galaxies are assigned to the blue sequence if $x<f_{\rm sat}^{\rm blue}(M_{r})$, and the red sequence otherwise.

4.2.2 Evolution of colours with redshift

The colour magnitude diagram evolves with redshift, as seen for example in figure 6 of Loveday et al. (2012) from GAMA. We therefore need to evolve the expressions given in Section 4.2.1 in order to produce a mock which has a realistic distribution of colours at these redshifts. In the GAMA data at high redshifts, only the brightest tip of the red and blue sequences can be seen, making it difficult to constrain their slopes. We therefore keep the slope of the red and blue sequence fixed with redshift.

We keep the red and blue sequences fixed at $z<0.1$, and evolve them with redshift as

[TABLE]

and

[TABLE]

respectively, where we stop evolving the mean of the sequences above $z=0.4$ in order to prevent too many high redshift galaxies being assigned as blue.

The mean satellite colour is also evolved as,

[TABLE]

and the fraction of blue galaxies is evolved as

[TABLE]

Fig. 14 shows the distribution of colours in the mock catalogue compared to GAMA for galaxies in different redshift and magnitude bins. Our parametrisation of the colour evolution is able to produce a good approximation to the GAMA colour distributions at all redshifts.

To evolve the luminosity function, we have assumed a fixed $Q$ parameter for all galaxies. We note that Loveday et al. (2012, 2015) hint that red and blue galaxies evolve differently, with a different $Q_{\rm red}$ and $Q_{\rm blue}$. However, the assumption of fixed $Q$ with this parametrisation is able to reproduce the observed colour-magnitude diagram.

4.3 Colour dependent k-corrections

In the mock catalogue, we use the HOD method to assign each galaxy an r-band absolute magnitude ${}^{0.1}M_{r}$, and the method outlined above to randomly generate a ${}^{0.1}(g-r)$ colour. However, the apparent magnitude, $r$, is the quantity which would be measured directly by the survey, and this is related to the absolute magnitude, ${}^{0.1}M_{r}$, through the equation

[TABLE]

where $d_{L}(z)$ is the luminosity distance in units of $h^{-1}\rm{Mpc}$, and ${}^{0.1}k(z)$ is the k-correction. The superscript $0.1$ denotes that the magnitude has been k-corrected to a reference redshift of $z_{\rm ref}=0.1$. In order to calculate an apparent magnitude for each galaxy in the mock catalogue, we use colour-dependent k-corrections derived from the GAMA survey, similar to those given in table 1 of McNaught-Roberts et al. (2014), except for k-correcting to $z_{\rm ref}=0.1$, rather than $z_{\rm ref}=0$.

The k-correction for each individual galaxy in GAMA is fit with a 4th order polynomial of the form

[TABLE]

The median k-correction is then found in 7 equally spaced bins of ${}^{0.1}(g-r)$ colour. Strictly speaking, the constant term in Eq. 38 should have the value $A_{4}=-2.5\log_{10}(1+z_{\rm ref})$ (Hogg et al., 2002); McNaught-Roberts et al. (2014) do not require this, but they end up with values of $A_{4}$ close to [math] at their $z_{\rm ref}=0$. We force our k-corrections to have the value $A_{4}=-2.5\log_{10}(1.1)$ at our $z_{\rm ref}=0.1$, but this only has a small effect on the k-corrections.

Using 7 distinct k-corrections based on colour leads to artificial features being added to the mock catalogue; for example step features can be seen in the colour-magnitude diagram at the boundaries between colour bins. In order to remove these features, we interpolate the k-corrections between the median colour in each bin.

Fig. 13 shows the polynomial fits to the k-corrections as a function of redshift. By definition, all the curves cross at ${}^{0.1}k(z=0.1)=-2.5\log_{10}(1.1)\approx-0.103$. The polynomial coefficients are shown in Table 1.

4.4 Colour dependent clustering in the mock

The projected correlation function of galaxies in the mock at low redshifts, split by red and blue galaxies, is shown in Fig. 15 for different bins in absolute magnitude and compared to the clustering in the corresponding volume limited samples from SDSS (Zehavi et al., 2011), where the red and blue samples are defined using the same colour cut as their equation 13. In the SDSS data, red galaxies are clustered more strongly than blue galaxies, since red elliptical galaxies are more likely to reside in more massive haloes, which are more strongly biased (Eisenstein et al., 2005). As the samples get fainter, the strength of the colour dependence becomes stronger. These trends are reproduced in the mock catalogue, using the modified satellite colour in Eq. 29.

Projected correlation functions for red and blue galaxies are also shown for higher redshift galaxies in Fig. 16, compared with the clustering seen in GAMA (Farrow et al., 2015). The galaxy samples are defined using the same ${}^{0.0}M_{r}$ magnitude ranges as figure 14 of Farrow et al. (2015), and using the same ${}^{0.0}(g-r)$ colour cut (their equation 4), where the superscript 0.0 denotes that these magnitudes are k-corrected to a reference redshift of $z_{\rm ref}=0$. The clustering of the red and blue galaxies in the mock is in reasonable agreement with the GAMA data.

5 Applications

As shown in Section 4, the galaxies in the mock catalogue have realistic clustering, which is in agreement with measurements from SDSS and GAMA. The galaxies also have a realistic distribution of ${}^{0.1}(g-r)$ colours at different redshifts. Future surveys, such as DESI and Euclid, aim to make measurements of the BAO and redshift space distortions, which probe larger scales than have been considered so far. Here, we show as an example some of the measurements that can be made using this mock catalogue.

5.1 BAO

As described in Section 2.1 and shown in Fig. 3, the large box size of the MXXL simulation enables the BAO feature to be seen clearly in the clustering of haloes. Here, we show that the BAO can also be seen in measurements of the redshift-space galaxy clustering. Fig. 17 shows the large scale redshift-space correlation function for several apparent magnitude threshold galaxy samples, using a redshift weighting $w(z)=1/(1+4\pi J_{3}\bar{n}(z))$ (Efstathiou et al., 1990), where $\bar{n}(z)$ is the number density of galaxies in the sample at redshift $z$, and $J_{3}=\int\xi r^{2}dr$, where we have assumed $4\pi J_{3}=3\times 10^{4}h^{3}\rm Mpc^{-3}$. The BAO peak can be seen in all samples, but the errors in the correlation function are largest for the $r<18.0$ sample. The $r<20.0$ sample contains fainter galaxies, and covers a larger volume, which greatly reduces the errors. For comparison, the crosses with error bars show measurements of the BAO from the Baryon Oscillation Spectroscopic Survey (BOSS) (Reid et al., 2016) for galaxies in the redshift range $0.2<z<0.5$ (Ross et al., 2017). The BAO scale in the mock catalogue is $\sim 7\%$ larger than is measured in BOSS. This is consistent with the difference in cosmology between that used in the MXXL simulation and the best fit to observations, including the BOSS results, and is mostly driven by the difference between $\Omega_{\rm m}=0.25$ in MXXL and $\Omega_{\rm m}=0.31$ in the Planck cosmology. The amplitude of the BAO peak in the mock catalogue also differs with the BOSS results, but we have not made any attempt to match the BOSS colour selection.

5.2 Redshift Space Distortions

The two-point correlation function $\xi(s,\mu)$, in bins of $s$ and $\mu$, can be decomposed into multipoles (Hamilton, 1992)

[TABLE]

where $s$ is the separation between a pair of galaxies in redshift space, $\mu=\cos\theta$ is the cosine of the angle between the vector $\mathbf{s}$ and the line of sight, and $P_{l}(\mu)$ is the lth order Legendre polynomial. The multipoles can be determined from the measured $\xi(s,\mu)$ by evaluating the integral

[TABLE]

Due to the symmetry $\xi(s,\mu)=\xi(s,-\mu)$, all odd-numbered terms are zero, and in linear theory, it is only the monopole, $\xi_{0}(s)$, quadrupole, $\xi_{2}(s)$, and hexadecapole, $\xi_{4}(s)$, that are non-zero. The amplitude of these multipoles depends on the strength of the redshift space distortions, and can provide a way to measure $f(z)\sigma_{8}(z)$ (Samushia et al., 2012).

The multipoles of the redshift-space correlation function of galaxies in the mock catalogue are shown in Fig. 18 for different volume limited magnitude threshold samples, and compared to measurements of clustering from SDSS (Guo et al., 2015). The monopole and quadrupole show reasonable agreement with the SDSS measurements, although the amplitude of the hexadecapole is a little high for some of the samples. Overall, the redshift space distortions in the mock catalogue look reasonably realistic, showing that the catalogue will be useful for future surveys that will take redshift space distortion measurements. We have extended the predictions beyond the range of the SDSS results, where they are easier to model and can be probed by surveys like DESI and Euclid.

6 Conclusions

For upcoming galaxy surveys, such as DESI and Euclid, it is important to have realistic mock catalogues in order to test and verify analysis tools, assess incompleteness and determine error covariances. The mock catalogues can also be used to make predictions and set expectations in advance of the first data from the survey.

We have outlined a method for creating a mock catalogue from the Millennium-XXL (MXXL) simulation. We first created a halo lightcone catalogue from the simulation, which we then populated with galaxies using a halo occupation distribution (HOD) scheme.

The halo lightcone catalogue is created from the simulation by finding the interpolated values of the position, velocity and mass of each halo at the redshift at which it crosses the observer’s lightcone. The halo lightcone catalogue covers the full sky, and extends to redshift $z=2.2$ with a mass resolution of $\sim 10^{11}h^{-1}\rm M_{\odot}$. Extending the catalogue to high redshifts requires multiple periodic replications of the MXXL box; these replications are only necessary to extend to redshifts greater than $z\sim 0.5$.

The halo catalogue is populated with galaxies using a Monte Carlo method, which randomly assigns galaxies with luminosities to dark matter haloes, following a HOD. This is an extension of the method outlined in Skibba et al. (2006) to a 5 parameter HOD. Galaxies in the mock catalogue are also assigned a ${}^{0.1}(g-r)$ colour, based on the Monte Carlo method of Skibba & Sheth (2009). A galaxy is assigned a colour based on its luminosity, redshift, and whether it is a central or a satellite galaxy, which is randomly drawn from a parametrisation of the SDSS colour-magnitude diagram.

The values of the HOD parameters we use are based on the best fitting HODs which reproduce the measured clustering from SDSS (Zehavi et al., 2011), but in Millennium cosmology. In the standard 5 parameter HOD model, the parameter $\sigma_{\log M}$ introduces Gaussian scatter in the luminosity of central galaxies for haloes of a fixed mass, which can lead to the unphysical crossing of HODs in different luminosity bins. We modify the HOD model so that this scatter follows a pseudo-Gaussian spline kernel, which prevents this crossing. The HODs are evolved with redshift such that they reproduce a target luminosity function, which is chosen to be the SDSS luminosity function at low redshift, and the luminosity function from GAMA at high redshifts. For a sample of galaxies with a fixed number density, the shape of the HOD is kept fixed with redshift, and the mass HOD parameters are all multiplied by the factor $f$ which is required to produce the correct number density. By construction, the mock catalogue reproduces the SDSS and GAMA luminosity functions, and the ratio of the HOD parameters $M_{1}/M_{\rm min}$ is constant with redshift.

We modify the parametrisation of the colour-magnitude diagram outlined in Skibba & Sheth (2009), and add evolution, such that the distribution of ${}^{0.1}(g-r)$ colours in the mock catalogue agrees with measurements from GAMA at different redshifts.

The galaxy catalogue is a flux limited mock galaxy catalogue, covering the full sky with an r-band magnitude limit of $r<20.0$ and median redshift $z\sim 0.2$. The angular and projected correlation functions of galaxies in the mock show good agreement with measurements from SDSS and GAMA, and the colour dependence of the clustering is reasonable. The BAO peak can be seen in the large scale clustering of galaxies in the mock, and galaxies show realistic redshift space distortions, making this mock useful for upcoming surveys which will measure these.

Here we have presented one mock galaxy catalogue, but to enable model inferences and place tight constraints on cosmological parameters, error covariances need to be determined. This requires the use of many mock catalogues, of the order of 100s to 10,000s. This could be achieved by coupling the HOD component of the mock with an approximate but fast method of generating halo catalogues (e.g. Manera et al., 2013; Monaco et al., 2013; Tassev et al., 2013; White et al., 2014; Avila et al., 2015; Chuang et al., 2015; Kitaura et al., 2015).

Acknowledgements

The authors would like to thank Sam Lyddon for making preliminary versions of several of the figures, and Jon Loveday for useful discussions. We thank the referee, Darren Croton, for several insightful comments. This work was supported by the Science and Technology Facilities Council (ST/M503472/1). AS, SMC, CMB & PN acknowledge the support of the Science and Technology Facilities Council (ST/L00075X/1). ZZ was partially supported by NSF grant AST-1208891 and NASA grant NNX14AC89G. REA acknowledges support from AYA2015-66211-C2-2. PN acknowledges the support of the Royal Society through the award of a University Research Fellowship and the European Research Council, through receipt of a Starting Grant (DEGAS-259586). IZ was supported by NSF grant AST-1612085, a CWRU Faculty Seed Grant and by ERC starting grant DEGAS-259586. IZ also acknowledges the hospitality of the ICC at Durham University.

GAMA is a joint European-Australasian project based around a spectroscopic campaign using the Anglo-Australian Telescope. The GAMA input catalogue is based on data taken from the Sloan Digital Sky Survey and the UKIRT Infrared Deep Sky Survey. Complementary imaging of the GAMA regions is being obtained by a number of independent survey programmes including GALEX MIS, VST KiDS, VISTA VIKING, WISE, Herschel-ATLAS, GMRT and ASKAP providing UV to radio coverage. GAMA is funded by the STFC (UK), the ARC (Australia), the AAO, and the participating institutions. The GAMA website is http://www.gama-survey.org/.

This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grants ST/H008519/1 and ST/K00087X/1, STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure.

Appendix A Database

The full sky halo and galaxy lightcone catalogues are made publicly available on the Theoretical Astrophysical Observatory database444https://tao.asvo.org.au/tao/ (Bernyk et al., 2016). The catalogues are also available at http://icc.dur.ac.uk/data/.

A.1 Halo catalogue

The halo catalogue contains a total of 5.1 billion haloes out to a redshift of $z=2.2$, and contains the following halo properties:

• •

$z_{\rm obs}$, the observed redshift, which takes into account the peculiar velocity of the halo

• •

$z_{\rm cos}$, the cosmological redshift, which ignores the effect of the peculiar velocity

• •

Right ascension, in degrees

• •

Declination, in degrees

• •

$M_{200\rm m}$, the mass enclosed by a sphere in which the average density is 200 times the mean density of the Universe, interpolated to the redshift at which the halo crosses the lightcone, in units of $10^{10}h^{-1}\rm M_{\odot}$

• •

$M_{200\rm c}$, the mass enclosed by a sphere in which the average density is 200 times the critical density of the Universe, interpolated to the redshift at which the halo crosses the lightcone, in units of $10^{10}h^{-1}\rm M_{\odot}$

• •

$V_{\rm max}$, the maximum circular velocity, in units of $\rm kms^{-1}$

• •

$R_{V_{\rm max}}$, the radius at which $V_{\rm max}$ occurs, in $h^{-1}\rm Mpc$

• •

$\sigma_{R_{\rm 200m}}$, velocity dispersion of particles within $R_{\rm 200m}$, in units of $\rm kms^{-1}$.

• •

Snapshot number in the MXXL simulation

• •

Halo id in the MXXL simulation

A.2 Galaxy catalogue

The galaxy catalogue contains 58.1 million galaxies with $r<20$, out to redshift $z=0.8$, and contains the following properties:

• •

$z_{\rm obs}$, the observed redshift, which takes into account the peculiar velocity of the galaxy

• •

$z_{\rm cos}$, the cosmological redshift, which ignores the effect of the peculiar velocity

• •

Right ascension, in degrees

• •

Declination, in degrees

• •

$M_{200\rm m}$ of the host halo, interpolated to the redshift at which the halo crosses the lightcone, in units of $10^{10}h^{-1}\rm M_{\odot}$

• •

Apparent r-band magnitude

• •

${}^{0.1}M_{r}-5\log h$, the rest frame absolute r-band magnitude, k-corrected to a reference redshift of $z_{\rm ref}=0.1$, with no evolutionary correction

• •

${}^{0.1}(g-r)$ colour, k-corrected to a reference redshift $z_{\rm ref}=0.1$

• •

A flag indicating whether the galaxy is a central or a satellite, and whether it is in a resolved or unresolved halo

• •

Snapshot number in the MXXL simulation

• •

Halo id in the MXXL simulation