An emulator-based halo model in modified gravity – I. The halo concentration-mass relation and density profile

The $f(R)$ gravity is a generalisation of Einstein’s general relativity. In $f(R)$ gravity, the Einstein-Hilbert action in GR has an additional term, which is a function of the Ricci scalar $R$,

where $M_{\rm Pl}=(8\pi G)^{-1/2}$ is the reduced Planck mass, $G$ is Newton’s constant, $g$ is the determinant of the metric $g_{\mu\nu}$ and $\mathcal{L}_{m}$ is the Lagrangian density for matter fields. Varying the action with respect to the metric $g_{\mu\nu}$ gives the modified Einstein equation,

where

is the Einstein tensor, $f_{R}\equiv\differential f(R)/\differential R$, $\nabla_{\mu}$ is the covariant derivative corresponding to the metric $g_{\mu\nu}$, $\square\equiv\nabla^{\alpha}\nabla_{\alpha}$ and $T^{{\text{m}}}_{\mu\nu}$ is the energy momentum tensor for matter.

Equation (2) is a fourth-order partial differential equation in $g_{\mu\nu}$. This equation can also be considered as the standard Einstein equation in GR with a new dynamical degree of freedom, $f_{R}$, which is dubbed the scalaron (e.g., Zhao et al., 2011). The equation of motion of $f_{R}$ can be obtained by taking the trace of Equation (2):

where $\rho_{\mathrm{m}}$ is the matter density.

For cosmological simulations in standard gravity, the Newtonian limit is commonly adopted. This includes the approximations that the gravitational and scalar fields are weak (such that their higher-order terms can be neglected) and quasi-static (so that the time derivatives of the fields can be neglected compared to their spatial derivatives). Most modified gravity simulations (including the ones used in this work) adopt this assumption. In the context of $f(R)$ gravity and the Newtonian limit, the modified Einstein equation (2) becomes

where $\Phi$ is the Newtonian potential, ${\nabla}$ is the 3-dimensional gradient operator, and an overbar denotes the cosmic mean of a quantity.

An $f(R)$ gravity model is fully specified by the functional form of $f(R)$. Here, we adopt the well-studied Hu-Sawicki model (Hu & Sawicki, 2007), which is given by

where $m^{2}\equiv\Omega_{{\text{m}}0}H_{0}^{2}$, and $c_{1},c_{2}$ and $n$ are free parameters. The parameter $n$ is a positive number, which is set to $n=1$ in this work as in most previous studies of this model (however see e.g., Li & Hu, 2011; Ramachandra et al., 2021; Ruan et al., 2022, for some examples of $n\neq 1$). With this functional form, we have

where $\bar{R}_{0}$ and $\bar{f}_{R0}$ are, respectively, the present-day values of the background Ricci scalar and $\bar{f}_{R}$. For brevity, we will adopt the following nomenclature to label models: the model with $-\log_{10}\left(|\bar{f}_{R0}|\right)=5.5$ will be called F5.5, and so on.

The remaining free parameter of the theory is the background value of the scalar field $f_{R}$ at redshift $z=0$, $\bar{f}_{R0}$. With a suitable choice of this parameter, $f(R)$ gravity reverts to GR in high-density regions – this is necessary to be consistent with solar system tests through the associated chameleon mechanism (Khoury & Weltman, 2004b, a). A larger value of $|\bar{f}_{R0}|$ means a stronger deviation from standard gravity. The F5 model is in slight tension with small-scale tests, see, e.g., Lombriser (2014) for a recent review of current cosmological constraints on $\bar{f}_{R0}$. But since we aim to test gravity on cosmic scales, models with such strength of MG are nevertheless still valuable to study: given their stronger deviations from GR compared to models with smaller $|\bar{f}_{R0}|$, they can lead to important insights into how the deviations from GR can affect large-scale cosmological observables such as weak lensing and galaxy clustering statistics. In order to fully explore the gravity testing capacities of upcoming cosmological observations, it is important to gain a detailed understanding of how these measures are influenced by possible modifications to gravity.

In the Dvali-Gabadadze-Porrati braneworld model (Dvali et al., 2000b), the universe is a four-dimensional brane embedded in a five-dimensional space-time (called the bulk). The gravitational action in this model is given by

where a superscript ⁽⁵⁾ denotes the quantity in the five-dimensional bulk. This model has a self-accelerating branch of solution (sDGP), which gives a natural explanation for the cosmic acceleration. However, the sDGP branch suffers from the ghost problems (Koyama, 2007) and cannot be considered as a physical model. Moreover, its predictions have been found to be inconsistent with observations such as cosmic microwave background (CMB) and local measurements of the Hubble parameter $H_{0}$ (e.g., Song et al., 2007; Fang et al., 2008).

The so-called normal branch DGP (nDGP) gravity (Koyama, 2007) cannot accelerate cosmic expansion by itself, so in order to explain the cosmological observations it has to introduce additional dark energy components. This model is nevertheless still of interest as a useful toy model that features the Vainshtein screening mechanism (Vainshtein, 1972). Here, we assume that there is an additional non-clustering dark energy component in this model, which results in its expansion history being identical to that of $\Lambda$CDM. The nDGP model provides an explanation of why gravity is much weaker than the other fundamental forces (Maartens & Koyama, 2010): all matter species are assumed to be confined to the brane, while gravity could propagate through (leak into) the extra spatial dimensions. There is one new free parameter in the nDGP model, which can be defined as the ratio of $G^{(5)}$ and $G$, and is known as the crossover scale,

The modified Friedmann equation in the normal branch DGP model is given by

where $\Omega_{\rm rc}\equiv 1/(4H_{0}^{2}r_{c}^{2})$, and $\Omega_{\rm DE0}$ is the density parameter of the additional dark energy component. The dimensionless quantity $H_{0}r_{c}$ is used to quantify departures from the standard gravity. If $H_{0}r_{c}\to\infty$ then Eqn. (11) returns to the $\Lambda$CDM case. A larger value of $H_{0}r_{c}$ means a weaker deviation from GR. Hereafter, an nDGP model with $H_{0}r_{c}=X$ will be referred to as N$X$. For example, a model with $H_{0}r_{c}=1$ is called N1.

The modified Poisson equation and the scalar field equation are given by (Koyama, 2007),

and

where $\varphi$ is a new scalar degree of freedom, $\delta\rho_{\rm m}=\rho_{\rm m}-\bar{\rho}_{\rm m}$ and

To construct emulators for dark matter halo properties, we use the FORGE and BRIDGE suites of $N$-body simulations (Arnold et al., 2022), covering $49$ $f(R)$ gravity and $49$ DGP models, along with $49$ $\Lambda$CDM counterparts. The simulations were performed using $1024^{3}$ dark matter particles in a cube of side $500\,h^{-1}\,\mathrm{Mpc}$ (hereafter the high-resolution runs, denoted HR) or $1500\,h^{-1}\,\mathrm{Mpc}$ (low-resolution runs, labelled LR), using the modified gravity version of the Arepo cosmological simulation code (Springel, 2010; Arnold et al., 2019b; Weinberger et al., 2020). The mass resolutions of the HR and LR runs are $9.1\times 10^{9}$ and $1.5\times 10^{12}(\Omega_{{\text{m}}0}/0.3)\,h^{-1}\,M_{\odot}$, respectively. The gravitational softening lengths of simulations are $15$ (HR) and $75\,h^{-1}\mathrm{kpc}$ (LR). The initial conditions were generated using second-order Lagrangian perturbation theory (2LPT, Crocce et al. (2006)) at $z_{\mathrm{init}}=127$. Each cosmology (also called “node”) has two independent realisations with the pairs of initial conditions selected to minimise the sample variance on large scales over the realisations. All nodes were initialised with the same (two) random seeds. See Section 3.2 of Arnold et al. (2022) for a detailed description.

The cosmological parameters were drawn from a Latin hypercube designed to efficiently sample a 4-dimensional parameter space (or a 3-dimensional space in the case of $\Lambda$CDM), as shown in Fig. 1, following a similar approach to that used in the cosmo-SLICS project (Harnois-Déraps et al., 2019). Since FORGE and BRIDGE were partly designed to emulate weak lensing statistics, they sampled directly in the composite structure growth parameter

instead of the physical matter fluctuation amplitude parameter $\sigma_{8}$. The use of $S_{8}$ accounts better for the degeneracy between $\Omega_{{\text{m}}0}$ and $\sigma_{8}$ in the cosmic shear analysis. The $49$ nodes form the training data set for each gravity model. It is common practice to use a small portion (called validation set) of the training set to determine whether the process has finished. We choose the nodes 11, 13, 22, 34, 36 as the validation set.

The following three standard cosmological parameters and one MG parameter are varied,

where $h\equiv H_{0}/(100\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1})$ and $H_{0}$ is the present day Hubble constant. The range of the parameters explored is

The density parameter of baryons was fixed to $\Omega_{\mathrm{b0}}=0.049199$ and massive neutrinos are ignored. The dark energy density is given by

The remaining cosmological parameter is the slope of the primordial curvature power spectrum normalised at $0.05\,\mathrm{Mpc}^{-1}$, which is $n_{{\text{s}}}=0.9652$ (Arnold et al., 2022).

We also need test data set to assess the performance of the emulators independently. In both cases, the test set consists of two models. The MG test cases are F5 and F6 for $f(R)$ gravity, and for DGP they are N1 and N5, and they share the same cosmological parameters, given by the fiducial Planck cosmology (Planck Collaboration et al., 2020),

Each test model has $8$ realisations.

The dark matter halo catalogues were obtained using the Subfind halo finder (Springel et al., 2001). The haloes were first identified using a fast parallel friends-of-friends (FOF) algorithm with link length set to $b=0.168$ times the mean interparticle separation. Spherical overdensity halo catalogues are then built out from the potential minimum of each FOF halo. The halo mass definition adopted is

where $\rho_{\mathrm{crit}}(z)\equiv 3H^{2}(z)/(8\pi G)$ is the critical density of the Universe, and $R_{\mathrm{200c}}$ is the spherical halo radius within which the spherically averaged mass density equals $200\rho_{\mathrm{crit}}$. Only main haloes with masses above $10^{12}\,h^{-1}M_{\odot}$ from the HR simulations are considered in this work. The halo catalogues at

are available for all nodes. Besides these common redshifts, Arnold et al. saved particle snapshots and halo catalogues at pre-selected redshifts to construct past light-cones for weak-lensing analysis, therefore enabling the emulation of halo properties as a function of redshift. In the main text, we present the performance of the emulators at redshift zero, since the measurement and emulation at other redshifts are performed in the same way.

The differential halo mass function (dHMF) quantifies the number density of haloes as a function of halo mass for a given cosmology $\mathcal{C}$ and redshift. It is denoted as

In the cumulative form, the cumulative HMF (cHMF) gives the number density of haloes above a given mass threshold $M$,

The HMF is measured by creating a histogram of the halo mass, which is affected by shot noise and sample variance, especially for massive haloes. Also, HMFs span many orders of magnitude in abundance, typically from $10^{-3}$ to $10^{-8}\,(h^{-1}\mathrm{Mpc})^{-3}$. Taking the logarithm of the HMF to reduce the dynamic range does not help much, since the interpolation errors in the logarithmic quantity would be exponentially amplified. To overcome these problems, a commonly used approach (e.g. followed by McClintock et al. 2019; Nishimichi et al. 2019; Cuesta-Lazaro et al. 2022) is to fit the measured HMF using fitting formulae like those proposed by Jenkins et al. (2001); Tinker et al. (2008), and then to emulate the mass and cosmology dependence of the fitting parameters. However, the performance of such fitting functions in MG simulations is not guaranteed (e.g. Schmidt et al., 2009; Gupta et al., 2022) and therefore may cause systematic errors.

We adopt an alternative method to emulate the HMF. The main ideas include: (1) Emulating the ratio between the simulation result for the HMF and a realistic fitting formula to reduce the dynamic range. (2) Considering the cumulative HMF instead of the differential one to allow for smaller steps in halo mass, thereby providing more training data.

Tinker et al. (2008) (hereafter T08) derived fitting functions for the HMF applicable over a wide range of halo masses and halo definitions, with a precision of $\lesssim 5\%$. We work with the ratio of the cumulative HMFs between the simulation measurements and the T08 fitting formulae¹¹1Numerically implemented in the Python module hmf (Murray et al., 2013, 2021).,

The HMF ratios are then interpolated across parameter space using neural networks to construct the emulator. To obtain the differential HMF for any given cosmology $\boldsymbol{\mathcal{C}}_{\text{any}}$, one can calculate the ratio for this cosmology using the trained emulator, multiply it with the T08 HMF and take the derivative. The emulation process is sketched in Fig. 2.

In each snapshot, we measure the number densities of haloes more massive than a series of masses, beginning at $10^{12}\,h^{-1}M_{\odot}$ and increasing in steps with a bin width of $\Delta\log{[M/(h^{-1}M_{\odot})]}=0.02$. The maximum mass varies across redshifts. Fig. 3 presents the ratios of HMFs for $49$ $f(R)$ gravity simulations at $z=0$. The ratios are gently varying functions over a lower dynamic range than HMFs themselves, therefore resulting in a higher emulation accuracy.

One of the most remarkable discoveries from cosmological $N$-body simulations was that dark matter haloes display a universal density profile (Navarro et al., 1996, 1997; Wang et al., 2020), from the host haloes of dwarf galaxies to those of massive galaxy clusters. Specifically, it was shown that the spherically averaged density profile of individual relaxed haloes can be described by the well-known Navarro, Frenk & White (NFW) profile. The NFW profile is described by two parameters, the characteristic density and scale radius of a halo, or equivalently the halo mass and concentration,

where $r_{-2}$ is the characteristic radius (also denoted as $r_{\mathrm{s}}$) of a halo at which the logarithmic slope of the density profile equal to $-2$,

and $\rho_{-2}\equiv\rho_{\mathrm{NFW}}(r=r_{-2})$. The halo concentration $c$ is defined as the ratio between the halo radius (which is adopted as $R_{200\mathrm{c}}$ in this work) and $r_{-2}$,

The other NFW parameter $\rho_{-2}$ is related to the concentration as

where $\Omega_{{\text{m}}}(a)\equiv\bar{\rho}_{{\text{m}}}(a)/\rho_{\mathrm{crit}}(a)$ is the matter density parameter at a given scale factor $a$; $f(c)\equiv\ln(1+c)-c/(1+c)$.

Previously, attention was focused on measuring halo concentrations for the best-fitting WMAP or Planck cosmologies, or similar models close by in parameter space (e.g. Gao et al., 2008; Prada et al., 2012; Diemer & Kravtsov, 2015; Klypin et al., 2016; Child et al., 2018). Such calibrated fitting functions cannot be simply extended beyond the cosmological and gravity models for which they have been tested. In order to overcome potential problems associated with the extrapolation of fitting functions to a wider range of cosmologies, we build emulators for the halo concentration-mass relation and halo density profiles.

We study the concentration-mass relation for haloes in the cosmologies covered by the FORGE and BRIDGE  simulation suites. We consider only the haloes containing more than $1\,000$ particles, corresponding to a mass of $10^{13}\,h^{-1}M_{\odot}$ for the fiducial Planck cosmology. For several reasons set out below we do not exclude unrelaxed haloes that contain a large amount of substructure as was done in some previous work (e.g. Neto et al., 2007; Prada et al., 2012; Klypin et al., 2016). First, our aim is to predict the halo profile as an ingredient of the halo model, instead of studying the formation and evolution of relaxed haloes. Galaxies are expected to reside in all haloes, regardless of their dynamic state. Second, excluding unrelaxed haloes would bias the concentration high because haloes in the rapid mass accretion stage tend to have low concentrations (Child et al., 2018). Third, such a cut removes typically $30\text{-}50\%$ haloes, which would make the measurement of halo properties less reliable, by introducing larger statistical errors.

To measure halo concentrations, we follow the approach taken by Mitchell et al. (2019) and briefly review the main aspects below. As mentioned in Section 2.3, we use $M_{\mathrm{200c}}$ and $R_{\mathrm{200c}}$ as the definitions of the halo mass and radius, respectively. The halo centre is adopted as the gravitational potential minimum. The halo particles are split into $20$ logarithmically spaced radial bins from the halo centre, covering the range $[0.05,1.00]R_{\mathrm{200c}}$. We then fit the NFW profile (Eqn. (25)) to the density in the radial bins of each halo, treating the characteristic density and concentration as free variables, by minimising an unweighted $\chi^{2}$,

We have checked that weighting the $\chi^{2}$ function by the number of particles in each radial bin has a negligible impact on the best-fitting concentration values recovered.

The mean value of the best-fitting concentration depends on technical details such as the halo finder used, and the number and range of the radial bins, which would have a non-negligible impact if the NFW profile is not a good fit to the halo profile (Meneghetti & Rasia, 2013; Dooley et al., 2014). It has been argued that using the median instead of the mean concentration would avoid such non-convergence (Diemer & Kravtsov, 2015). Neto et al. (2007) found that the median of the concentration depends only weakly on the radial range used in the fit, provided that the unrelaxed haloes are removed and $r_{\text{min}}\geq 0.05R_{\text{vir}}$ in the fitting.

To test the robustness of the halo concentration measurement, we use three $N$-body simulations from Mitchell et al. (2021) with the same cosmology and particle number, $N_{\text{p}}=1024^{3}$, and different box sizes, $L=200,500$ and $1000\,h^{-1}\mathrm{Mpc}$. We then bin the halo particles, fit the NFW profile and calculate the mean and median values of the concentrations in each mass bin, keeping the maximum radius fixed at $R_{\mathrm{200c}}$ and checking the results for four different $r_{\text{min}}$ values: $(0.05,0.07,0.10$ and $0.13)R_{\mathrm{200c}}$. The results are shown in Fig. 4. All concentration-mass relations are well fitted by a power law,

with $c_{0}$ the amplitude and $\alpha$ the index, as represented by the coloured bands in the figure. The median concentrations (as well as the mean) and the best-fitting amplitude are still sensitive to the minimum radius, with a relative difference of up to $10$ per cent. However, the power indices are practically the same for different fitting ranges.

The convergence test also confirms our choice of the minimum particle-number cut applied to define the halo sample used. The concentrations from the lower resolution simulations start deviating from those of the higher resolution runs around the halo mass corresponding to $\sim 800$ times the particle mass. This pivot mass also depends weakly on the radial range used to fit the density profile, with smaller minimum scales having larger pivot masses.

The normalised halo density profile, $u(r|M)$, that appears in the halo model can be estimated by measuring the halo-mass cross-correlation functions from an $N$-body simulation. As mentioned in Section 1, the halo model assumes that the matter density field consists of a superposition of haloes at locations ${\boldsymbol{x}}_{i}$ with masses $M_{i}$, so that the matter field can be written as

where

• •

  The summation is for all haloes; the same results can be derived by taking the summation over all microcells which are made to be so small that each cell contains no more than one halo centre (e.g. Peebles, 1980);

• •

  $\delta^{\mathrm{(D)}}(\cdots)$ is the Dirac delta function;

• •

  $\displaystyle\sum_{i}\delta^{\mathrm{(D)}}(M-M_{i})\,\delta^{\mathrm{(D)}}(\boldsymbol{x}^{\prime}-\boldsymbol{x}_{i})\equiv\frac{\differential{n({\boldsymbol{x}}^{\prime};M)}}{\differential{M}}$ is the local HMF, whose integral over the halo mass gives the halo number density field at the field point ${\boldsymbol{x}}^{\prime}$:

  $\displaystyle\int\differential{M}\frac{\differential{n({\boldsymbol{x}}^{\prime};M)}}{\differential{M}}=n_{{\text{h}}}({\boldsymbol{x}}^{\prime}),$

  (33)

  and its ensemble average gives the common dHMF,

  	$\displaystyle\expectationvalue{\frac{\differential{n({\boldsymbol{x}}^{\prime};M)}}{\differential{M}}}$	$\displaystyle=\expectationvalue{\sum_{i}\delta^{\mathrm{(D)}}(M-M_{i})\,\delta^{\mathrm{(D)}}(\boldsymbol{x}^{\prime}-\boldsymbol{x}_{i})}$
  		$\displaystyle=\derivative{n(M)}{M};$		(34)

• •

  $\displaystyle u\quantity(|\boldsymbol{x}-\boldsymbol{x}_{i}|\Big{|}M_{i})\equiv\frac{\rho_{\mathrm{h}}\quantity(|\boldsymbol{x}-\boldsymbol{x}_{i}|\Big{|}M_{i})}{M_{i}}$ denotes the density profile of a halo centred at $\boldsymbol{x}_{i}$, which is also assumed to be spherically symmetric and depends only on its mass,

  $\displaystyle u\quantity(|\boldsymbol{x}-\boldsymbol{x}_{i}|\big{|}M_{i})=\begin{cases}\rho_{\mathrm{h}}\quantity(|\boldsymbol{x}-\boldsymbol{x}_{i}|\Big{|}M_{i})/M_{i},&|\boldsymbol{x}-\boldsymbol{x}_{i}|\leq R_{\text{halo}};\\ 0,&|\boldsymbol{x}-\boldsymbol{x}_{i}|>R_{\text{halo}};\end{cases}$

  (35)

  and $u$ is normalised:

  $\displaystyle\int\differential[3]{\boldsymbol{x}}u\quantity(|\boldsymbol{x}-\boldsymbol{x}_{i}|\big{|}M_{i})$

  $\displaystyle=1.$

  (36)

For two different populations of objects with overdensity fields $\delta_{a}({\boldsymbol{x}})$ and $\delta_{b}({\boldsymbol{x}})$, the two-point cross-correlation function is defined as

where $\expectationvalue{\cdots}$ denotes ensemble averaging. Assuming statistical isotropy reduces $\xi_{ab}({\boldsymbol{r}})$ to a function of separation only. In the case of the cross-correlation between halo centres and the matter field, the cross-correlation function is given by

where ${\boldsymbol{x}}^{\prime}\equiv{\boldsymbol{x}}+{\boldsymbol{r}}$, and the halo number density fluctuation field is defined as

Now we derive the expressions for the halo and matter fields in Eqn. (39) within the halo model framework. In realistic $N$-body simulations, halo catalogues always have a finite halo mass range, limited by the simulation force/mass resolution and the box size. We can define the halo selection function for a given mass range (MR) as

In practice, MR can be a narrow mass interval, $[M,M+\differential{M})\approx[M,M+\Delta M)$, or a mass threshold interval, $[M,\infty)$. For a given halo sample, the corresponding halo number density field can be expressed as

To obtain the expression for $\langle n_{{\text{h}}}({\boldsymbol{x}})\,\rho_{{\text{m}}}({\boldsymbol{x}}^{\prime})\rangle$ in the halo model, we can plug the expressions for the halo and mass fields, Eqns. (32) and (42), into Eqn. (39). The full expression can be found in Eqns. (7) and (8) of García et al. (2021). We focus only on the internal structure of the halo and, therefore, on the 1-halo term, which reads

where $\differential{n(M)}$ is the number density of halos in the mass range $[M,M+\differential{M})$. The 1-halo term of the halo-mass correlation function is

We are interested in the average density profile of haloes in a narrow mass interval $[M,M+\differential{M})$. However, the noise level of the simulation measurements for such halo properties and statistics is high because of the low number density. We bypass this problem by measuring $\xi_{{\text{hm}}}^{\mathrm{1h}}\quantity(r|>M)$ (which involves more haloes and therefore is a smoother function) and taking the partial derivative with respect to mass,

To gain a preliminary impression of the emulation process, and to guide the design of emulators in the future, we perform emulation tests under ideal conditions. The halo properties are generated for a limited number of randomly selected cosmologies using analytical methods or fitting formulae, which are noise-free mappings from cosmologies to halo properties. Then we use these data to train neural networks and emulate the “true” model. To evaluate the performance of the emulator, we compare the true values with the emulator predictions using independent test data sets.

The cosmologies of the training set cover $50$ flat geometry $(w_{0}\text{-}w_{a})$ CDM models (Linder, 2003), where the equation-of-state parameter for dark energy is parameterised in terms of the expansion factor, $a$, as

A key aspect of building emulators is an efficient sampling scheme. As the training dataset, the $50$ cosmologies were sampled using optimal minimax distance sliced Latin hypercube designs (Ba et al., 2015) in a seven-dimensional cosmological parameter space,

as shown by the grey dots in Fig. 5. The range of parameters is

while the upper and lower parameter limits depart significantly from the current best-fitting $\Lambda$CDM background cosmology from the Planck satellite (Planck Collaboration et al., 2020). We also generate two test data sets that were not used in the training: both consist of $500$ random cosmologies; one set covers the same parameter range as that of the training set, and the other one covers the inner half-region (in terms of the length per dimension, instead of volume) of the parameter space. The cosmologies in the full- and half-range test data sets are shown in green and blue dots in Fig. 5, respectively.

We choose to emulate two basic properties of haloes in the tests: the concentration-mass relation $c(M)$, and the cumulative HMF $\bar{n}_{{\text{h}}}(>M)$. For given cosmologies, we generate the $c(M)$ relation calibrated in Prada et al. (2012), using the publicly available Python toolkit COLOSSUS (Diemer, 2018), and compute the Tinker et al. (2008) HMF with the Python package hmf (Murray et al., 2013, 2021). As described in Section 3.1 and Fig. 2, we train the emulators directly on the ratio of the cumulative HMF between the target HMF and a fitting formula to reduce the dynamic range and improve the interpolation accuracy. We choose the HMF calibrated in Jenkins et al. (2001) as the reference.

The upper panels of Fig. 6 show the halo properties calculated using the fitting formulae and emulators. The fractional errors are shown in the lower sub-panels. The results show that the emulator reproduces the analytical halo $c(M)$ relation with a sub-percent error in the 7D parameter space with only $50$ training models. The performance of the HMF emulator is even better than that of the concentration emulator, although the HMF data span $\sim 20$ orders of magnitude. The median absolute error of the emulator predictions is lower than $1$ per cent for halo masses $M\lesssim 10^{16}\,h^{-1}M_{\odot}$.

We also note that the emulator precision is generally different at the edge and centre of the parameter space, as revealed by the green and blue lines in the bottom panels of Fig. 6. This suggests that we should design the parameter space to be wider than the existing cosmological constraints. The parameter space in our ideal tests is designed to be wide enough that covers some extreme cosmologies, such as $\Omega_{\text{m0}}=0.7$ and $h=0.5$.

In general, ideal emulation tests show that under noise-free conditions, neural network emulators can provide accurate interpolations in high-dimensional parameter space, using $50$ efficiently sampled models as the training set. In the next section, we will present the cosmic emulation in the real situation: the data are measured from simulations and, therefore, are influenced by sample variance and noise.

As discussed in Section 3.1, we train the neural network emulators directly on the ratio of the cumulative HMF between simulation measurements and the fitting formula calibrated in Tinker et al. (2008) to reduce the dynamic range and improve the emulation accuracy. Fig. 7 compares the cumulative and the differential HMFs from simulations and emulators at $z=0$, for the $\Lambda$CDM, $f(R)$ gravity and DGP models. The differential HMF of a given mass bin centred on $\log{M_{i}}$ is obtained by

where $\mathrm{bw}\equiv\log{M_{i+1}}-\log{M_{i}}$ is the bin width. The lower sub-panels show the fractional difference between the emulator prediction and the measured HMFs in a given mass bin,

The performance of the emulator on the training data set is shown by the thin lines of Fig. 7. The emulator achieves sub-percent accuracy in reproducing most of the cumulative HMF for halo masses between $10^{12}$ and $10^{14}\,h^{-1}M_{\odot}$. The residuals of the differential HMF obtained using Eqn. (55) are slightly larger but still show $\lesssim 2$ per cent scatter around zero, with a mean consistent with zero. The thick lines in Fig. 7 show the emulator predictions for three test models which are not used in the training, as described in Section 2.3. We again find percent-level agreement between the emulator predictions and simulation results. Furthermore, the fluctuations of the residuals in the test set are much smaller than those of the training models, since each test cosmology has $8$ realisations and the sample variance of the measured dHMFs is suppressed. This confirms that the errors of the emulator predictions are mainly random instead of systematic.

As shown in Section 5.2 and Fig. 4, the halo concentrations measured from simulations are sensitive to the radial range used in the fitting, which indicates that this is not the optimal way to describe the halo density profiles. However, the power law index is not sensitive to the range adopted, which indicates that we can treat the amplitude of the concentration-mass relation as a free parameter to take into account this variation. We build an emulator for the $c(M)$ relation taking $r_{\text{min}}=0.10R_{\mathrm{200c}}$ as a representative value. The performance of the emulator at $z=0$ is shown in Fig. 8. The fractional errors are sub-percent for most of the cosmologies in the training and test data sets.

As discussed in Section 3.3, the average halo profile can be estimated from the halo-mass cross-correlation function. The halo profile $u(r|M)$ is directly related to the matter density field cross-correlated with the halo sample in a narrow mass range $[M,M+\Delta M]$. However, such correlation functions measured from simulations would be rather noisy because of the low halo number density. To feed the neural networks with smoother data, we measure the cross-correlation functions between the matter field and the halo samples with fixed number densities, $\xi_{{\text{hm}}}(r|\boldsymbol{\mathcal{C}},\bar{n}_{{\text{h}}})$. We then use the HMF emulator to translate $\xi_{{\text{hm}}}(r|\boldsymbol{\mathcal{C}},\bar{n}_{{\text{h}}})$ as a function of number density into $\xi_{{\text{hm}}}(r|\boldsymbol{\mathcal{C}},M)$ as a function of mass, according to Equation (47).

We measure $\xi_{{\text{hm}}}(r|\boldsymbol{\mathcal{C}},\bar{n}_{{\text{h}}})$ using the high-performance code Corrfunc (Sinha & Garrison, 2020) for the halo number densities in logarithmically spaced bins over the range

using a bin width of $\Delta\log_{10}\bar{n}_{{\text{h}}}=0.05$. The separation $r$ is split into $30$ logarithmically-spaced bins from $0.05\,h^{-1}\mathrm{Mpc}$ (three times the force resolution) to $3\,h^{-1}\mathrm{Mpc}$. Furthermore, to reduce the dynamic range of the data vector, we opt to emulate $r^{2}\xi_{{\text{hm}}}(r)$ instead of $\xi_{{\text{hm}}}(r)$ itself. The upper-left panel of Fig. 9 shows $\xi_{{\text{hm}}}\big{(}r|\bar{n}_{{\text{h}}}=10^{-3.5}\,(h^{-1}\,\mathrm{Mpc})^{-3}\big{)}$ at $z=0$ for the $49$ $\Lambda$CDM gravity cosmologies along with the test models.

The average halo profile is only related to the 1-halo term of $\xi_{{\text{hm}}}$. The halo-mass correlation enters the transition between 1- and 2-halo terms as the scale increases. To estimate the range of the one-halo term, we only consider the scales below $R_{\mathrm{200c}}$, which is related to the adopted mass definition $M_{\mathrm{200c}}$ as

We then build emulators for $\xi_{{\text{hm}}}(r|\boldsymbol{\mathcal{C}},\bar{n}_{{\text{h}}})$ at each redshift, to reduce the number of features and minimise emulation errors. The lower left sub-panel of Fig. 9 shows the fractional difference between the simulation measurements and the emulator predictions, in the training set along with the test models. The emulator achieves sub-percent accuracy for the both the training and the test models.

The average halo density profile can be estimated from $\xi_{\text{hm}}$ using Eqn. (47). In the right panel of Fig. 9, we compare this type of profile with the NFW profiles combined with three concentration-mass relations from this work, Klypin et al. (2016) and Diemer & Joyce (2019), in five halo mass bins. We also fit the average profile with an NFW form, using the the data over the range of $[0.1,1.0]R_{\text{200c}}$.

In the right sub-panel of Fig. 9, we check the relative difference between the average profiles measured from the simulations and the NFW fits. There is a $\sim 5\%$ discrepancy between the two types of profiles, regardless of the concentration-mass relation, which shows that the differences between the NFW profiles with different $c(M)$ relations are small. This level of difference is consistent with the results discovered in Section 5.1.2 of Nishimichi et al. (2019).

We adopt the halo occupation distribution (HOD) (e.g. Zheng et al., 2005) prescription to model the average number of galaxies in a halo as a function of halo mass. The occupation of central galaxies is parameterised as a Bernoulli distribution, whereas that of satellites is a Poisson distribution (Zheng et al., 2005). Both distributions are described by their mean occupation number,

The galaxy number density $\bar{n}_{{\text{g}}}$ can then be obtained by integrating the HMF weighted by the mean occupation,

Following Cuesta-Lazaro et al. (2022), we adopt the HOD model in Zheng et al. (2007) by introducing the following HOD parameters,

The mean occupation number for central galaxies is given by

The mean central HOD, $\expectationvalue{N_{{\text{c}}}}\!(M)\,$, can be interpreted as the probability that a halo with mass $M$ hosts a central galaxy. The mean central HOD considered here has the asymptotic behaviour that $\expectationvalue{N_{{\text{c}}}}\to 0$ for haloes with $M\ll M_{\text{min}}$, while $\expectationvalue{N_{{\text{c}}}}\to 1$ for haloes with $M\gg M_{\text{min}}$.

The mean satellite HOD is parameterised as

where

Following the commonly-used prescription, we assume that satellite galaxies reside only in a halo that already hosts a central galaxy. Hence, in the above equation, satellite galaxies can only reside in haloes with $\expectationvalue{N_{{\text{c}}}}=1$. Then we assume that the number distribution of satellite galaxies in a given host halo follows the Poisson distribution with mean $\lambda_{{\text{s}}}(M)$: