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

Cheng-Zong Ruan1,2, Carolina Cuesta-Lazaro1,3, Alexander Eggemeier4, Baojiu
Li1, Carlton M. Baugh1,3, Christian Arnold1, Sownak Bose1, César Hernández-
Aguayo5,6, Pauline Zarrouk7 and Christopher T. Davies8
1Institute for Computational Cosmology, Department of Physics, Durham
University, South Road, Durham DH1 3LE, UK
2Institute of Theoretical Astrophysics, University of Oslo, 0315 Oslo, Norway
3Institute for Data Science, Durham University, South Road, Durham DH1 3LE, UK
4Argelander-Institut für Astronomie, Auf dem Hügel 71, D-53121 Bonn, Germany
5Max-Planck-Institut fur Astrophysik, Karl-Schwarzschild-Str 1, D-85748
Garching, Germany
6Excellence Cluster ORIGINS, Boltzmannstrasse 2, D-85748 Garching, Germany
7Sorbonne Université, Université Paris Diderot, Sorbonne Paris Cité,
CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE),
74 place Jussieu, F-75252, Paris Cedex 5, France
8Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, 81679
Munich, Germany
Argelander FellowE-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form )

###### Abstract

In this series of papers we present an emulator-based halo model for the non-
linear clustering of galaxies in modified gravity cosmologies. In the first
paper, we present emulators for the following halo properties: the halo mass
function, concentration-mass relation and halo-matter cross-correlation
function. The emulators are trained on data extracted from the FORGE and
BRIDGE suites of $N$-body simulations, respectively for two modified gravity
(MG) theories: $f(R)$ gravity and the DGP model, varying three standard
cosmological parameters $\Omega_{\mathrm{m0}},H_{0},\sigma_{8}$, and one MG
parameter, either $\bar{f}_{R0}$ or $r_{\mathrm{c}}$. Our halo property
emulators achieve an accuracy of $\lesssim 1\%$ on independent test data sets.
We demonstrate that the emulators can be combined with a galaxy-halo
connection prescription to accurately predict the galaxy-galaxy and galaxy-
matter correlation functions using the halo model framework.

###### keywords:

dark energy – large-scale structure of Universe – cosmology: miscellaneous –
cosmology: theory.

††pubyear: 2023††pagerange: An emulator-based halo model in modified gravity –
I. The halo concentration-mass relation and density profile–An emulator-based
halo model in modified gravity – I. The halo concentration-mass relation and
density profile

##  1 Introduction

Ongoing and upcoming galaxy surveys, such as those that will be made with the
Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration et al., 2016),
the Vera Rubin Observatory (LSST Science Collaboration et al., 2009) and
Euclid (Laureijs et al., 2011; Amendola et al., 2013; Troja et al., 2022) will
map the large scale structure (LSS) of the Universe with unprecedented
statistical precision. Measurements of the large-scale structure can
potentially be used to unveil the nature of the dark matter and dark energy,
and to look for any deviation from the predictions of general relativity (GR).
Theories of gravity beyond GR – modified gravity (MG) models – can explain the
observed accelerated expansion of the Universe without invoking a cosmological
constant (e.g. Koyama, 2018; Ferreira, 2019). Studies of such models will not
only shed light on the nature of the cosmic acceleration, but also serve as
useful tests of GR on cosmic scales.

The impact of modifications to GR has been well studied in terms of the cosmic
expansion history and the large scale structure, i.e., at the background and
linear perturbation levels. In the late universe, the growth of LSS eventually
enters the non-linear regime over a wide range of length scales, and linear
theory predictions cease to be valid. This point becomes even more acute in
the context of MG, given that such models have intrinsically non-linear
features, such as screening mechanisms and non-linear field equations for new
degrees of freedom, which cannot be captured by linear theories (e.g. Li et
al., 2013). Often a choice is made to exclude small scale data, thereby losing
a wealth of information from high signal-to-noise ratio measurements. Such
nonlinearities must be properly incorporated into theoretical modelling if one
wishes to make the best use of the current and next generation cosmological
surveys to constrain cosmological parameters and test gravity theories.

A fully non-linear treatment – $N$-body simulations – is essential to
accurately solve the non-linear dynamics of cosmic structure formation (see
e.g. Kuhlen et al., 2012; Angulo & Hahn, 2022, for recent reivews). The main
hurdle of $N$-body simulations is their expensive computational cost. A Monte
Carlo Markov chain (MCMC) analysis, usually used to confront theoretical
predictions with data, requires sampling at least $10^{4}$-$10^{5}$ models in
the cosmological parameter space. Probing such a large number of models using
simulations is computationally prohibitive. The situation is even worse for MG
models, which usually involve partial differential equations governing the new
physics. Current MG simulations can take between $2$ to $\mathcal{O}(10)$
times longer than standard $\Lambda$CDM simulations with the same
specifications (e.g. Li et al., 2012; Arnold et al., 2019a).

There are several approaches to dealing with the non-linear regime in addition
to simulations. $N$-body simulation results can be used to develop
phenomenological or semi-analytical fitting formulae to describe the
statistical properties of matter and dark matter haloes, such as the halo mass
function calibrated by Tinker et al. (2008), and the halofit prescription for
the matter power spectrum (Smith et al., 2003; Takahashi et al., 2012; Smith &
Angulo, 2019; Mead et al., 2021) and bispectrum (Takahashi et al., 2020). The
most up-to-date version of halofit implemented by Mead et al. (2021) achieves
an accuracy of $5\%$ down to deeply non-linear scales. However, such
parametric fits may no longer be fit for purpose with the advent of next-
generation surveys that promise to reach measurements with per cent level
precision.

The halo model (Neyman & Scott, 1952; Ma & Fry, 2000; Peacock & Smith, 2000;
Seljak, 2000; Cooray & Sheth, 2002; Schmidt, 2016; Philcox et al., 2020) is a
successful analytical description of the LSS in the non-linear regime. In this
framework, all matter, including galaxies and any other tracers, is assumed to
reside within haloes. Then, the problem of predicting the clustering can be
split into the following steps:

  * •

the abundance of haloes as a function of halo mass, i.e. the halo mass
function (HMF);

  * •

the distribution of tracers around the halo centre, i.e. the halo density
profile, usually assumed to be a Navarro-Frenk-White (NFW) profile (Navarro et
al., 1996, 1997) specified by a halo mass-concentration relation; and

  * •

the clustering of the haloes themselves, e.g., the halo two-point correlation
function (TPCF) in configuration space and the halo auto power spectrum in
Fourier space.

These basic properties of dark matter haloes constitute the halo model
ingredients. The halo model provides a physically motivated description of the
clustering statistics and is flexible enough to be extended to incorporate new
physics such as massive neutrinos and baryonic feedback (e.g. Massara et al.,
2014; Bose et al., 2021; Carrilho et al., 2022), as well as models beyond
$\Lambda$CDM and GR (e.g. Barreira et al., 2014; Lombriser et al., 2014; Hu et
al., 2018; Cataneo et al., 2019).

The halo model ingredients can be derived from analytic methods. For example,
the HMF can be predicted using the spherical collapse model of the linear
matter density field (Press & Schechter, 1974) and the excursion set formalism
(Bond et al., 1991; Sheth et al., 2001), whereby the dependence of the HMF on
redshift and cosmology is expressed in terms of the root-mean-square
fluctuations in the linear matter power spectrum. Jenkins et al. (2001) found
that this universality of the HMF holds at an approximate level. As simulation
predictions improved further, it was discovered that the redshift evolution of
the mass function, even for $\Lambda$CDM, deviates from the universal
prediction at the $5\text{-}10$ per cent level, and several new fitting
formulae were proposed (e.g. Tinker et al., 2008; Courtin et al., 2011).
Moreover, the universality of the HMF is broken further in extensions of
$\Lambda$CDM (e.g. Bhattacharya et al., 2011, for $w$CDM) and modified gravity
models (e.g. Schmidt et al., 2010; Lam & Li, 2012; Li & Efstathiou, 2012;
Lombriser et al., 2013; Cataneo et al., 2016; Gupta et al., 2022). In order to
obtain even tighter constraints on cosmological parameters and to test gravity
theories, one therefore needs to proceed beyond the traditional approaches,
given their limitations in accuracy and coverage of parameter space.

In this series of papers, we develop simulation-based theoretical templates
called emulators, to obtain accurate predictions for basic halo properties as
a function of halo mass and (modified gravity) cosmology, and to construct
accurate predictions for clustering observables in preparation for ongoing and
future galaxy surveys. There have been several previous works on the emulation
of cosmological quantities in the $\Lambda$CDM model and its extensions, such
as Heitmann et al. (2006); Habib et al. (2007), the Coyote Universe (Heitmann
et al., 2010; Heitmann et al., 2009; Lawrence et al., 2010), PkANN (Agarwal et
al., 2012, 2014), the Mira-Titan Universe (Heitmann et al., 2016; Lawrence et
al., 2017; Bocquet et al., 2020), Kwan et al. (2013); Kwan et al. (2015),
Aemulus (DeRose et al., 2019; McClintock et al., 2019; Zhai et al., 2019),
Dark Quest (Nishimichi et al., 2019; Kobayashi et al., 2020; Miyatake et al.,
2020; Cuesta-Lazaro et al., 2022), matryoshka (Donald-McCann et al., 2022) and
AbacusSummit (Maksimova et al., 2021; Yuan et al., 2022), as well as in non-
standard cosmologies, such as Winther et al. (2019); Ramachandra et al.
(2021); Arnold et al. (2022); Brando et al. (2022); Harnois-Déraps et al.
(2022).

To build emulators we use the machine-learning interpolation technique of
neural networks, which allows us to predict halo properties for any given
cosmology within the range of parameters covered by the training data set. We
use the FORGE (F-Of-R Gravity Emulator) and BRIDGE (BRaneworld-Inspired Dgp
Gravity Emulator) modified gravity $N$-body simulations described in Arnold et
al. (2022), which together cover a very broad range of parameters in two MG
theories: $f(R)$ gravity and the DGP model. The emulated halo properties
incorporate all the complicated effects on non-linear scales, such as the non-
linear halo bias, the halo exclusion effect and the screening mechanism.

Following the spirit of the Dark Quest project (Nishimichi et al., 2019;
Cuesta-Lazaro et al., 2022), we do not perform an end-to-end emulation of
galaxy clustering statistics in the joint parameter space of cosmological and
galaxy-halo connection models. Instead, we develop emulators for each halo
property separately, and assemble these ingredients within the halo model
framework to construct analytical predictions of galaxy clustering statistics.
This emulator-based halo model gives us the flexibility to insert different
prescriptions of galaxy-halo connection for different galaxy samples.
Moreover, emulators for basic halo properties themselves are very useful. For
example, calibrating the cosmology dependence of the HMF is crucial to control
the systematic uncertainty in galaxy cluster abundance studies (e.g.
McClintock et al., 2019; Bocquet et al., 2020).

The layout of this paper is as follows. In Section 2, we present a short
description of the modified gravity theories studied here and a brief overview
of the FORGE and BRIDGE $N$-body simulation suites. In Section 3, we outline
the measurement and post-processing of the halo properties from the
simulations. Section 4 describes the construction of the halo property
emulators using neural networks, and Section 5 shows their performance in
reproducing the simulation results. In Section 6, we demonstrate how these
emulators can be combined with a galaxy-halo connection prescription to
predict galaxy statistics.

Throughout this paper, we use $\log$ to denote the base-$10$ logarithm,
$\log\equiv\log_{10}$, and $\ln$ to indicate the natural logarithm. Unless
otherwise stated, we use a subscript $0$ to denote the present-day value of a
physical quantity and an overbar for the background value of a quantity.

##  2 Modified Gravity Theories and Simulations

We briefly describe the two modified gravity models analysed in this work,
$f(R)$ gravity (Hu & Sawicki, 2007) and the Dvali-Gabadadze-Porrati (DGP)
brane-world models (Dvali et al., 2000a). These are two of the most widely
studied MG models and, as we discuss below, are representative examples of the
two main classes of screening mechanisms, which make them good test-beds for
generic MG models. For more detailed descriptions of these models, we refer
the reader to Sotiriou & Faraoni (2010) and De Felice & Tsujikawa (2010) for
$f(R)$ gravity, and Sahni & Shtanov (2003) and Maartens & Koyama (2010) for
DGP models.

###  2.1 $f(R)$ gravity

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$,

$\displaystyle S=\int\differential^{4}x\sqrt{-g}\,\quantity{\frac{M^{2}_{\rm
Pl}}{2}\quantity[R+f(R)]+\mathcal{L}_{m}}\ ,$ (1)

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,

$\displaystyle G_{\mu\nu}+f_{R}R_{\mu\nu}-\quantity(\frac{1}{2}f-\square
f_{R})g_{\mu\nu}-\nabla_{\mu}\nabla_{\nu}f_{R}=8\pi GT^{{\text{m}}}_{\mu\nu},$
(2)

where

$G_{\mu\nu}\equiv R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R,$ (3)

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):

$\displaystyle\square f_{R}=\frac{1}{3}\quantity(R-f_{R}R+2f+8\pi
G\rho_{\mathrm{m}})\ ,$ (4)

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

$\displaystyle{\nabla}^{2}\Phi$ $\displaystyle\approx\frac{16\pi
G}{3}a^{2}(\rho_{\mathrm{m}}-\bar{\rho}_{\mathrm{m}})+\frac{1}{6}a^{2}\quantity[R(f_{R})-\bar{R}]\
,$ (5) and the equation of motion of the scalaron reduces to
$\displaystyle{\nabla}^{2}f_{R}$
$\displaystyle\approx-\frac{1}{3}a^{2}\quantity[R(f_{R})-\bar{R}+8\pi
G(\rho_{\mathrm{m}}-\bar{\rho}_{\mathrm{m}})]\ ,$ (6)

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

$\displaystyle f(R)=-m^{2}\frac{c_{1}(-R/m^{2})^{n}}{c_{2}(-R/m^{2})^{n}+1}\
,$ (7)

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

$\displaystyle
f_{R}=-\left|\bar{f}_{R0}\right|\left(\frac{\bar{R}_{0}}{R}\right)^{n+1},$ (8)

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.

###  2.2 The Dvali-Gabadadze-Porrati (DGP) model

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

$\displaystyle
S=\int_{\text{brane}}\differential[4]{x}\sqrt{-g}\quantity(\frac{R}{16\pi
G})+\int_{\text{bulk}}\differential[5]{x}\sqrt{-g^{(5)}}\quantity(\frac{R^{(5)}}{16\pi
G^{(5)}}),$ (9)

where a superscript (5) 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,

$\displaystyle r_{c}\equiv\frac{1}{2}\frac{G^{(5)}}{G}\,.$ (10)

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

$\displaystyle\frac{H(a)}{H_{0}}=\sqrt{\Omega_{{\text{m}}0}a^{-3}+\Omega_{\rm
DE0}(a)+\Omega_{\rm rc}}-\sqrt{\Omega_{\rm rc}}\ ,$ (11)

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),

$\nabla^{2}\Phi=4\pi Ga^{2}\delta\rho_{\rm m}+\frac{1}{2}\nabla^{2}\varphi\,,$
(12)

and

$\nabla^{2}\varphi+\frac{r_{c}^{2}}{3\beta\,a^{2}c^{2}}\left[(\nabla^{2}\varphi)^{2}-(\nabla_{i}\nabla_{j}\varphi)^{2}\right]=\frac{8\pi\,G\,a^{2}}{3\beta}\delta\rho_{\rm
m}\,,$ (13)

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

$\displaystyle\beta(a)$ $\displaystyle\equiv
1+2H\,r_{c}\left(1+\frac{\dot{H}}{3H^{2}}\right)$
$\displaystyle=1+\frac{\Omega_{{\text{m}}0}a^{-3}+2\Omega_{\Lambda}0}{2\sqrt{\Omega_{\rm
rc}(\Omega_{{\text{m}}0}a^{-3}+\Omega_{\Lambda 0})}}.$ (14)

###  2.3 Modified gravity $N$-body simulations

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.

Figure 1:  Visualisation of the cosmological parameters
$\boldsymbol{\mathcal{C}}$ (Equations (16)-(18)) covered in the FORGE and
BRIDGE simulations. The $44$ training cosmologies, $5$ validation and $2$ test
cosmologies are shown as black dots, purple triangles and red stars,
respectively. $\Lambda$CDM models corresponding to
$\log{|\bar{f}_{R0}|}=-\infty$ and $\log{(H_{0}r_{\mathrm{c}})=\infty}$ are
not shown in the last two rows. (source code)

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

$\displaystyle S_{8}\equiv\sigma_{8}\sqrt{\frac{\Omega_{{\text{m}}0}}{0.3}},$
(15)

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,

$\displaystyle\boldsymbol{\mathcal{C}}$
$\displaystyle=\quantity{\Omega_{\mathrm{m0}},h,S_{8}},\ $
$\displaystyle\text{for $\Lambda$CDM},$ (16)
$\displaystyle\boldsymbol{\mathcal{C}}$
$\displaystyle=\quantity{\Omega_{\mathrm{m0}},h,S_{8},\log_{10}|\bar{f}_{R0}|},\
$ $\displaystyle\text{for $f(R)$},$ (17)
$\displaystyle\boldsymbol{\mathcal{C}}$
$\displaystyle=\quantity{\Omega_{\mathrm{m0}},h,S_{8},\log_{10}(H_{0}r_{\mathrm{c}})},\
$ $\displaystyle\text{for DGP},$ (18)

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

$\begin{split}0.11&<\Omega_{{\text{m}}0}<0.54\\\ 0.61&<h<0.81\\\
0.6&<S_{8}<0.9\\\ -6.2&<\log_{10}|\bar{f}_{R0}|<-4.6\\\
-0.45&<\log_{10}(H_{0}r_{\mathrm{c}})<1.0.\end{split}$ (19)

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

$\displaystyle\Omega_{\Lambda 0}=1-\Omega_{{\text{m}}0}.$ (20)

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),

$\displaystyle\Omega_{\mathrm{m0}}=0.31315,\quad\Omega_{\Lambda
0}=0.68685,\quad\Omega_{\mathrm{b0}}=0.049199,$ $\displaystyle
h=0.6737,\quad\sigma_{8}=0.82172,\quad n_{{\text{s}}}=0.9652.$ (21)

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

$M_{\mathrm{200c}}\equiv\frac{4\pi}{3}(R_{\mathrm{200c}})^{3}\times
200\rho_{\mathrm{crit}},$

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

$\displaystyle z=0.00,0.25,0.50,0.75,1.00,1.25,1.50,1.75\ \text{and}\ 2.00$

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.

##  3 Data Set

In this section we describe the measurement and post-processing of the halo
properties from the simulations, including the halo mass function,
concentration-mass relation and halo-matter cross-correlation function.

###  3.1 Halo Mass Functions

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

$\displaystyle\derivative{n(M;z,\boldsymbol{\mathcal{C}})}{M}\quad\text{or}\quad\derivative{n(\log
M;z,\boldsymbol{\mathcal{C}})}{\log M}.$ (22)

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

$\displaystyle
n(>M;z,\boldsymbol{\mathcal{C}})=\int_{M}^{\infty}\differential{m}\derivative{n(m;z,\boldsymbol{\mathcal{C}})}{m}.$
(23)

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 formulae111Numerically
implemented in the Python module hmf (Murray et al., 2013, 2021).,

$\displaystyle
r(M;\boldsymbol{\mathcal{C}},z)\equiv\frac{n_{{\text{h}}}^{\mathrm{sim}}(>M;\boldsymbol{\mathcal{C}},z)}{n_{{\text{h}}}^{\mathrm{T08}}(>M;\boldsymbol{\mathcal{C}},z)}.$
(24)

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.

Figure 2:  Flow chart illustrating the process of emulating halo mass
functions. We start from measuring the cumulative HMF from the simulations in
the training data set. For each simulation, we calculate the cHMF predicted by
the fitting formula of Tinker et al. (2008) for the same cosmology
$\boldsymbol{\mathcal{C}}$. We then interpolate the ratios between the two
cHMFs across parameter space using neural networks. To obtain the commonly
used differential HMF for any given cosmology
$\boldsymbol{\mathcal{C}}_{\mathrm{any}}$, we can calculate the ratio,
multiply the ratio by the Tinker et al. (2008) function and take the
derivative.

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.

Figure 3: Cumulative halo mass function (cHMF) ratios between simulation
measurements for $49$ $f(R)$ gravity models and the fitting formula calibrated
by Tinker et al. (2008) for the same cosmology (except for the MG parameter
$\bar{f}_{R0}$, which is set to zero), at redshift $0$. The dynamic range of
the ratios are significantly decreased compared with cHMFs themselves,
therefore increasing the emulation accuracy. (source code)

###  3.2 The individual halo density profile and the concentration-mass
relation

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,

$\displaystyle\rho_{\mathrm{NFW}}(r|r_{-2},\rho_{-2})=\frac{4\rho_{-2}}{\displaystyle\quantity(\frac{r}{r_{-2}})\quantity(1+\frac{r}{r_{-2}})^{2}},$
(25)

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$,

$\displaystyle\left.\frac{\differential{\log{\quantity[\rho_{\mathrm{NFW}}(r)]}}}{\differential{\log{r}}}\right|_{r=r_{-2}}=-2,$
(26)

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}$,

$\displaystyle c\equiv\frac{R_{200\mathrm{c}}}{r_{-2}}.$ (27)

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

$\displaystyle\rho_{-2}=\frac{\rho_{\mathrm{crit}}(a)}{4}\frac{200}{3\Omega_{{\text{m}}}(a)}\frac{c^{3}}{f(c)},$
(28)

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}$,

$\displaystyle\chi^{2}=\sum_{i}\quantity[\log{\rho_{\mathrm{sim}}(r_{i})}-\log{\rho_{\mathrm{NFW}}(r_{i}|c,\rho_{-2})}]^{2}.$
(29)

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,

$\displaystyle c(M)=c_{0}M^{\alpha},$ (30)

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.

Figure 4: Convergence test of halo concentration measurements. We fit the NFW
profile and measure the halo concentration from three $N$-body simulations
with the same particle number and different box sizes, $L200$ (the highest
resolution run, thick solid lines), $L500$ (the medium run, thin dashed lines)
and $L1000$ (the low resolution run, thin solid lines). Colours represent
results obtained for different minimum radial ranges in the fitting, as shown
by the key, with the maximum radial range fixed at $R_{\mathrm{200c}}$ in each
case. The vertical lines present the critical mass scales lower than which the
concentration measurements are unreliable due to insufficient resolution,
which are $\sim 800$ times the mass resolutions $m_{\text{particle}}$, as long
as $r_{\text{min}}>0.05R_{\text{200c}}$. The colour-shaded bands show the
power fitting ($c(M)=c_{0}M^{\alpha}$) to the measured $c(M)$ relations.
(source code)

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.

###  3.3 The averaged halo profile estimated from the halo-mass correlation
function

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

$\displaystyle\rho_{\text{m}}(\boldsymbol{x})$
$\displaystyle=\sum_{i}M_{i}\,u\quantity(|\boldsymbol{x}-\boldsymbol{x}_{i}|\Big{|}M_{i})$
(31)
$\displaystyle={\small\int\differential[3]{\boldsymbol{x}^{\prime}}\int\differential{M}\bigg{[}\sum_{i}\delta^{\mathrm{(D)}}(M-M_{i})\,\delta^{\mathrm{(D)}}(\boldsymbol{x}^{\prime}-\boldsymbol{x}_{i})}$
$\displaystyle\phantom{=\small\int\differential[3]{\boldsymbol{x}^{\prime}}\int\differential{M}}\quad
M\,u\quantity(|\boldsymbol{x}-\boldsymbol{x}^{\prime}|\Big{|}M)\bigg{]},$ (32)

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

$\displaystyle\xi_{ab}({\boldsymbol{r}})\equiv\expectationvalue{\delta_{a}({\boldsymbol{x}})\,\delta_{b}({\boldsymbol{x}}+{\boldsymbol{r}})},$
(37)

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

$\displaystyle\xi_{{\text{hm}}}(r)$
$\displaystyle\equiv\expectationvalue{\delta_{{\text{h}}}({\boldsymbol{x}})\,\delta_{{\text{m}}}({\boldsymbol{x}}+{\boldsymbol{r}})}$
(38)
$\displaystyle=\frac{1}{\bar{n}_{{\text{h}}}\bar{\rho}_{{\text{m}}}}\expectationvalue{n_{{\text{h}}}({\boldsymbol{x}})\,\rho_{{\text{m}}}({\boldsymbol{x}}^{\prime})}-1,$
(39)

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

$\displaystyle\delta_{{\text{h}}}({\boldsymbol{x}})\equiv\frac{n_{{\text{h}}}({\boldsymbol{x}})-\bar{n}_{{\text{h}}}}{\bar{n}_{{\text{h}}}}.$
(40)

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

$\displaystyle\phi(m|\mathrm{MR})\equiv\begin{cases}1,&\text{if}\
m\in\text{MR},\\\ 0,&\text{otherwise}.\end{cases}$ (41)

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

$\displaystyle
n_{{\text{h}}}(\boldsymbol{x})=\sum_{i}\delta^{\mathrm{(D)}}(\boldsymbol{x}-\boldsymbol{x}_{i})\,\phi(m_{i}|\mathrm{MR}).$
(42)

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

$\displaystyle\expectationvalue{n_{{\text{h}}}({\boldsymbol{x}})\,\rho_{{\text{m}}}({\boldsymbol{x}}^{\prime})}_{\mathrm{1h}}=\int\differential{m}\frac{\differential{n}}{\differential{m}}m\,\phi(m|\mathrm{MR})\,u(r|m)$
(43)
$\displaystyle=\begin{cases}\differential{n(M)}M\,u(r|M),&\mathrm{MR}=[M,M+\differential{M}),\\\
\displaystyle\int_{M}^{\infty}\differential{m}{\frac{\differential{n}}{\differential{m}}}\,m\,u(r|m),&\mathrm{MR}=[M,\infty),\end{cases}$
(44)

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

$\displaystyle\xi_{{\text{hm}}}^{\mathrm{1h}}\quantity(r|M)$
$\displaystyle=\frac{M\,u(r|M)}{\bar{\rho}_{{\text{m}}}}-1=\frac{\rho_{\mathrm{h}}(r|M)}{\bar{\rho}_{{\text{m}}}}-1,$
(45) $\displaystyle\xi_{{\text{hm}}}^{\mathrm{1h}}\quantity(r|\\!>\\!M)$
$\displaystyle=\frac{1}{\bar{\rho}_{{\text{m}}}\,\bar{n}_{{\text{h}}}(\\!>\\!M)}\int_{M}^{\infty}\differential{m}{\frac{\differential{n}}{\differential{m}}}\,m\,u(r|m)-1.$
(46)

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,

$\displaystyle u(r|M)$
$\displaystyle=-\frac{\bar{\rho}_{\mathrm{m}}}{M}\quantity[\frac{\differential{n(M)}}{\differential{M}}]^{-1}$
$\displaystyle\quad\times\frac{\partial}{\partial
M}\quantity{\bar{n}_{{\text{h}}}(>M)\,\Big{[}1+\xi_{\text{hm}}^{\mathrm{1h}}\quantity(r|>M)\Big{]}}.$
(47)

##  4 Using Neural Networks for Regression Problems

Given a data set comprised of independent variables (also called features) and
dependent variables (labels), there exists an unknown underlying function
mapping the inputs to the outputs. We can use supervised learning algorithms
to approximate this function, which falls in the category of a regression
problem. In the context of structure formation, the features consist of
cosmologies $\boldsymbol{\mathcal{C}}$, redshifts $z$, halo masses $M$ or
number densities $n_{{\text{h}}}$, etc. The labels of interest include basic
properties of haloes such as halo mass functions, concentration-mass relations
and correlation functions. In the case of structure formation, the “functions”
between the features and labels are known but expensive: we can run $N$-body
simulations given a set of cosmological parameters $\boldsymbol{\mathcal{C}}$,
save the snapshot at a given redshift $z$, identify haloes and measure the
properties of haloes. But it is computationally intractable to perform
$\gtrsim\mathcal{O}(10^{4})$ cosmological simulations to explore the parameter
space in a typical MCMC analysis.

As shown in previous works on cosmic emulation, such as Nishimichi et al.
(2019); DeRose et al. (2019); Cuesta-Lazaro et al. (2022), and as we will
report in the following sections, it is possible to construct emulators for
halo properties by running affordable numbers (e.g., $50\text{-}100$) of
simulations and interpolating in a high-dimensional parameter space. Emulators
can give accurate predictions for halo properties for new models without
running additional simulations. Thanks to the development of algorithms and
computing power, statistics and machine learning provide us with a wealth of
tools to solve such regression problems.

In a regression analysis, the typical progress of approximating a function can
be summarised as:

  * •

Define a functional form
$\boldsymbol{y}=f(\boldsymbol{x}|\boldsymbol{\theta})$ with adjustable or
trainable parameters $\boldsymbol{\theta}$. For example, in the simplest and
most common form — linear regression — the labels are assumed to be linear
combinations of the features and the coefficients are the trainable
parameters.

  * •

Define a loss function on the training data set $D$ to quantify the difference
between the real and predicted values of the target, e.g., the sum of the
absolute differences (which reduces the weight of outliers),

$\displaystyle
L(\boldsymbol{\theta}|D)\equiv\sum_{(\boldsymbol{x}_{i},\boldsymbol{y}_{i})\in
D}\big{|}\boldsymbol{y}_{i}-f(\boldsymbol{x}_{i}|\boldsymbol{\theta})\big{|}.$

  * •

Find the optimal parameters ${\boldsymbol{\theta}}_{\star}$ to minimise the
loss function (train the model).

Gaussian process (GP) regression (e.g. Williams & Rasmussen, 2006) has been
widely adopted in cosmological emulation projects, such as Dark Quest
(Nishimichi et al., 2019), Aemulus (DeRose et al., 2019), the Coyote Universe
(Heitmann et al., 2010), the Mira-Titan Universe (Bocquet et al., 2020),
Cosmic Emulators (Kwan et al., 2013) and FORGE (Arnold et al., 2022). GP
regression is non-parametric, i.e., no specific functional form is assumed.
However, GPs are not easy to apply to large data sets because of their
$\mathcal{O}(N^{3})$ scaling where $N$ is the size of the training data.
Therefore, the aforementioned projects usually emulate matter, halo and/or
galaxy properties using a combination of principal component analysis, to
reduce the dimensionality of the data vector, and GP, to fit the dependence of
the principal component coefficients on cosmology.

The machine learning algorithm we adopt is a fully connected neural network
(e.g. Bishop et al., 1995; Alom et al., 2019; Zhang et al., 2021). This is a
type of artificial neural network in which all neurons in one layer are
connected to the neurons in the next layer. The neural network algorithm has
been widely applied to cosmology (e.g., Agarwal et al., 2012, 2014; Jennings
et al., 2019; Kobayashi et al., 2020; Cuesta-Lazaro et al., 2022), and its
performance has been shown to be competitive or sometimes better than other
methods. Neural networks have a strong fitting ability, as reflected by the
universal approximation theorem (Cybenko, 1989; Hornik et al., 1989;
Goodfellow et al., 2016). However, this theorem does not provide a means to
construct optimised neural networks, but merely guarantees their existence.
Also, this strong fitting ability also makes neural networks more susceptible
to the overfitting problem compared to GPs. Therefore, we need to carefully
check the emulator performance using independent test data sets and tune the
network architecture so that the generalisation error is successfully
suppressed. Moreover, we show dimensionality reduction is not necessary when
using neural networks, which also improves the accuracy of the emulator
predictions.

A neural network is an interconnection of neurons arranged in a series of
layers, with each neuron in a layer connected to all other neurons in adjacent
layers with different weightings. One can impart values on to the neurons of
the first layer (called the input layer), have a number of hidden layers and
finally obtain the output from the last layer. For example, in this work, the
neural networks emulating the halo mass function at fixed redshift have five
nodes in the input layer, corresponding to the halo mass and four cosmological
parameters, and one node in the last layer outputting the HMF,

$\displaystyle
n\big{(}>M|\Omega_{\mathrm{m0}},h,S_{8},\log_{10}|\bar{f}_{R0}|\
\text{or}\,\log_{10}|H_{0}r_{{\text{c}}}|\big{)}.$ (48)

Neural networks use activation functions to impart non-linearities into the
fitting. Rectified Linear Unit (ReLU) (Agarap, 2018) is the most commonly used
activation function in current neural networks used to add non-linearities in
the mapping between inputs and outputs, which is defined as

$\mathrm{ReLU}(x)=\max(0,x),$ (49)

where $x$ is the output of the previous layer of the neural network. Note that
the activations of ReLU are not differentiable at $x=0$. Here, however, we are
interested in functions that are differentiable with respect to their inputs
and, in particular, with respect to the cosmological parameters. Therefore,
throughout, we use Gaussian error linear unit (GELU) (Hendrycks & Gimpel,
2016) as the activation function instead,

$\mathrm{GELU}(x)=0.5x\quantity[1+\erf\left(\frac{x}{\sqrt{2}}\right)].$ (50)

To find the optimal parameters $\boldsymbol{\theta}_{\star}$ that reproduce
the halo properties measured in the $N$-body simulations, we minimise the
L1-norm loss function,

$\mathcal{L}=\frac{1}{N}\sum_{i=0}^{N}|y^{i}_{\mathrm{true}}-y^{i}_{\mathrm{predicted}}|$
(51)

using the Adam optimiser (Kingma & Ba, 2014). Compared to the mean squared
error, the loss of L1 reduces the importance given to outlier errors. To avoid
fine-tuning the learning rate, we adopt a learning rate scheduler that reduces
the learning rate by a factor of $10$ every time the validation loss does not
improve after $30$ epochs. We also stop the training process when the
validation loss does not improve after $100$ epochs. This iterative reduction
of the learning rate allows the model to quickly learn the broad
characteristics of the data and later reduce the errors with a smaller
learning rate. The initial learning rate is always set to $0.01$.

###  4.1 Ideal emulation tests

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

$\displaystyle w(a)=w_{0}+w_{a}(1-a).$ (52)

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,

$\displaystyle\boldsymbol{\mathcal{C}}=\Big{\\{}\Omega_{{\text{m}}0},\Omega_{\mathrm{b0}},h,\sigma_{8},n_{{\text{s}}},w_{0},w_{a}\Big{\\}},$
(53)

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

$\begin{split}0.1&<\Omega_{{\text{m}}0}<0.7,\quad
0.02<\Omega_{\mathrm{b0}}<0.06,\\\ 0.5&<h<0.9,\quad 0.5<\sigma_{8}<1.2,\\\
0.92&<n_{{\text{s}}}<0.99,\\\
-1.3&<w_{0}<-0.7,\quad-0.1<w_{a}<0.1,\end{split}$ (54)

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.

Figure 5: Visualisation of the seven-parameter $(w_{0}\text{-}w_{a})$CDM
cosmologies studied in the ideal emulation tests. Grey dots show the training
set including $50$ nodes. Green dots represent full-range test set consisting
of $500$ nodes covering the same range as that of the training set. Blue dots
show the half-range test set including $500$ nodes in the inner half region of
the parameter space. (source code)

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

Figure 6: Ideal emulation tests. Top panels: The halo concentration-mass
relation (left) and halo mass function (right) of the training (grey), full-
range test (green) and half-range test (blue) sets, from the analytical
methods (dots) and emulators (lines). Bottom panels: The absolute value of the
relative difference between the models and emulators. The solid and dashed
lines present the median and 90-th percentile of the emulation fractional
errors among the training (black), test (green) and half-range test (blue)
sets, which include $50$, $500$ and $500$ cosmologies, respectively. (source
code 1, 2)

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.

Table 1: Summary of the neural network configurations for emulating the halo properties. The architecture of a neural network is specified from the input to output layer as $(N_{\text{input}},N_{\text{hidden1}},N_{\text{hidden2}},\cdots,N_{\text{output}})$ with $N$ the number of neurons in each layer.  Halo Property | Feature | Label |  Neural Network Architecture | Activation
---|---|---|---|---
HMF | $\boldsymbol{\mathcal{C}},M_{\mathrm{200c}}$ | Equation (24) | $(5,64,32,1)$ | GELU
Concentration | $\boldsymbol{\mathcal{C}},M_{\mathrm{200c}}$ | $c_{200\mathrm{c}}$ | $(5,32,16,1)$ | GELU
$\xi_{\text{hm}}$ | $\boldsymbol{\mathcal{C}},n_{\mathrm{h}}$ | $\\{r_{i}^{2}\xi_{\text{hm}}(r_{i})\\}_{i=1}^{N=30}$ | $(5,128,32,30)$ | GELU

##  5 Results

In this section we demonstrate the ability of a fully connected neural network
to reproduce the halo properties measured from FORGE and BRIDGE simulations.
We train different emulators for each gravity theory: $\Lambda$CDM, $f(R)$
gravity and DGP. The configurations of the neural networks for emulating three
halo properties are summarised in Table 1.

###  5.1 The emulator for the halo mass function

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

$\displaystyle\left.\derivative{n(M)}{\log{M}}\right|_{\log{M_{i}}}\approx\frac{n(>\log{M_{i}}-\frac{\mathrm{bw}}{2})-n(>\log{M_{i}}+\frac{\mathrm{bw}}{2})}{\mathrm{bw}},$
(55)

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,

$\displaystyle\frac{{\mathrm{HMF}}_{\text{emu}}-{\mathrm{HMF}}_{\text{sim}}}{{\mathrm{HMF}}_{\text{sim}}}.$
(56)

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.

Figure 7:  Cumulative (the first row) and differential (the second row) halo
mass functions from simulation measurements and emulator predictions at $z=0$,
for $\Lambda$CDM (left), $f(R)$ gravity (middle) and DGP (right). In each
panel, the thin lines show the results of the $49$ cosmologies in the training
data set, and the thick lines represent those of the test models which were
not used in the construction of the emulators. In the lower sub-panel, we
compare the relative differences between simulations and emulators. The dark
and light grey bands denote $\pm 1\%$\- and $\pm 2\%$-level errors,
respectively. The differential halo mass functions are obtained by finite
difference of the cumulative HMFs according to Eqn. (55). (source code 1, 2)

###  5.2 The emulator for the concentration-mass relation

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.

Figure 8:  Comparison between the halo concentration-mass relations from
simulations (points) and emulators (lines), for the $49$ $f(R)$ gravity
cosmologies in the training set (grey) and test models F5 (green), F6
(orange), N1 (blue), N5 (purple) and the fiducial Planck cosmology (denoted as
GR, in red). The results for the $\Lambda$CDM and DGP training sets are
similar to those for $f(R)$ gravity and are therefore not shown here. The
lower sub-panel shows the relative difference of the halo concentration
between simulations and emulators. The sub-percent differences between the
emulator and simulation results are much smaller than the differences between
the results for different cosmologies. (source code)

###  5.3 The emulator for the halo-mass cross-correlation function

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

$\displaystyle\log_{10}\quantity[\frac{\bar{n}_{{\text{h}}}}{(h^{-1}\mathrm{Mpc})^{-3}}]=[-5.1,-2.9],$
(57)

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

$\displaystyle
M_{\mathrm{200c}}(z)=\frac{4\pi}{3}(R_{\mathrm{200c}})^{3}200\rho_{\mathrm{crit}}(z).$
(58)

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

Figure 9: Left panel: Halo-mass cross-correlation functions from simulations
and emulators at $z=0$, for the $49$ $\Lambda$CDM gravity cosmologies in the
training set (grey), five test models F5 (cyan), F6 (orange), N1 (blue), N5
(purple) and the fiducial Planck $\Lambda$CDM model (red). In the lower sub-
panel, we compare the relative difference between the simulations and
emulators. Right panel: Comparison of the normalised halo density profiles,
$\rho(r|M)/M$, truncated at $r=R_{\mathrm{200c}}$, for the fiducial Planck
$\Lambda$CDM model (node 0) at $z=0$. The average halo profiles estimated from
$\xi_{\text{hm}}(r)$ according to Equation (47) are represented by points. The
solid lines show the fits to an NFW profile. The dashed, dotted and dash-
dotted lines show the NFW profiles with three different concentration-mass
relations: this work, Diemer & Joyce (2019) and Klypin et al. (2016). Colours
denote different halo masses. The lower sub-panel shows the relative
difference between the NFW profiles (lines) and the average profiles (points).
(source code)

##  6 Emulator Applications

With the emulators for the halo mass function and density profile as
ingredients, we are able to predict galaxy clustering statistics using the
halo model framework (Cooray & Sheth, 2002), and therefore fit galaxy
clustering measurements in the joint parameter space of cosmology and a
galaxy-halo connection model. In this section, we demonstrate that the
emulator-based halo model reproduces the signals measured from the mock HOD
catalogues generated with the same specifications, such as the cosmology, HOD
prescription, satellite profile and/or concentration-mass relation.

###  6.1 Galaxy two-point correlation function

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,

$\displaystyle\expectationvalue{N_{{\text{g}}}}(M)=\langle
N_{{\text{c}}}\rangle(M)+\langle N_{{\text{s}}}\rangle(M).$ (59)

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

$\displaystyle\bar{n}_{{\text{g}}}=\int\differential{M}\derivative{n(M)}{M}\big{[}\expectationvalue{N_{{\text{c}}}}\\!(M)\,+\expectationvalue{N_{{\text{s}}}}\\!(M)\big{]}.$
(60)

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

$\displaystyle\boldsymbol{\mathcal{G}}=\Big{\\{}\underbrace{M_{\text{min}},\sigma_{\log{M}}}_{\displaystyle\boldsymbol{\mathcal{G}}_{\text{cen}}},\underbrace{M_{1},\kappa,\alpha}_{\displaystyle\boldsymbol{\mathcal{G}}_{\text{sat}}}\Big{\\}}.$
(61)

The mean occupation number for central galaxies is given by

$\displaystyle\expectationvalue{N_{{\text{c}}}}\\!(M|{\boldsymbol{\mathcal{G}}})=\frac{1}{2}\quantity[1+\erf\quantity(\frac{\log{M}-\log{M_{\text{min}}}}{\sigma_{\log{M}}})].$
(62)

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

$\expectationvalue{N_{{\text{s}}}}\\!(M|{\boldsymbol{\mathcal{G}}})=\expectationvalue{N_{{\text{c}}}}\\!(M|{\boldsymbol{\mathcal{G}}})\lambda_{{\text{s}}}(M),$
(63)

where

$\lambda_{{\text{s}}}(M)=\quantity(\frac{M-\kappa
M_{\text{min}}}{M_{1}})^{\alpha}.$ (64)

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)$:

$\displaystyle P(N_{{\text{s}}}|N_{{\text{c}}}=1)$
$\displaystyle=\frac{\quantity[\lambda_{{\text{s}}}(M)]^{N_{{\text{s}}}}\mathrm{e}^{-\lambda_{{\text{s}}}(M)}}{N_{{\text{s}}}!},$
(65) and $\displaystyle P(N_{{\text{s}}}|N_{{\text{c}}}=0)$
$\displaystyle=\delta^{\mathrm{Kr}}_{N_{{\text{s}}},0},$ (66)

where $\delta^{\mathrm{Kr}}_{ij}$ stands for the Kronecker delta.

Given the HOD model, we populate dark matter haloes in $8$ test simulations of
the fiducial Planck cosmology with mock galaxies and measure the galaxy
clustering signals, using the following randomly selected HOD parameters:

$\begin{split}&\log{[M_{\text{min}}/(h^{-1}M_{\odot})]}=12.5,\
\sigma_{\log{M}}=0.6915,\\\
&\kappa=0.51,\log{[M_{1}/(h^{-1}M_{\odot})]}=12.9,\ \text{and}\
\alpha=0.9168.\end{split}$ (67)

We can also express the galaxy two-point correlation function (TPCF) in terms
of dark matter halo properties in the halo model framework. First, we split
the one- and two-halo terms into correlations of central and satellite
galaxies as

$\displaystyle\xi_{{\text{gg}}}(r)$
$\displaystyle=\xi_{\mathrm{cs}}^{\mathrm{1h}}(r)+\xi_{\mathrm{ss}}^{\mathrm{1h}}(r)$
$\displaystyle\phantom{=}+\xi_{\mathrm{cc}}^{\mathrm{2h}}(r)+\xi_{\mathrm{cs}}^{\mathrm{2h}}(r)+\xi_{\mathrm{ss}}^{\mathrm{2h}}(r).$
(68)

The terms involving both centrals and satellites lead to a convolution of the
halo profiles and/or the halo TPCF, following
$\displaystyle\int\differential[3]{{\boldsymbol{x}}}u(x|M)\,u\quantity(|{\boldsymbol{x}}+\boldsymbol{r}|\Big{|}M)$.
It is therefore more convenient to compute these terms in Fourier space, where
convolutions in coordinate space become simple products of the Fourier modes.
Here, we focus on the one-halo term only. The two-halo term involving the
emulation of halo clustering will be the topic of the subsequent papers in
this series. The expressions for the 1-halo term of the galaxy TPCF after the
central-satellite split are given by

$\displaystyle\xi(r)$
$\displaystyle=\int_{0}^{\infty}\frac{\differential{k}}{(2\pi)^{3}}4\pi
k^{2}\,\frac{\sin(kr)}{kr}\,P(k),$ (69) $\displaystyle
P_{\mathrm{cs}}^{\mathrm{1h}}(k)$
$\displaystyle=\frac{1}{\bar{n}_{{\text{g}}}^{2}}\int\differential{M}{\color[rgb]{0,0,1}\derivative{n(M)}{M}}\expectationvalue{N_{{\text{c}}}}\\!(M)\,\lambda_{{\text{s}}}(M)\,{\color[rgb]{0,0,1}u_{\mathrm{s}}(k|M)},$
(70) $\displaystyle P_{\mathrm{ss}}^{\mathrm{1h}}(k)$
$\displaystyle=\frac{1}{\bar{n}_{{\text{g}}}^{2}}\int\differential{M}{\color[rgb]{0,0,1}\derivative{n(M)}{M}}\expectationvalue{N_{{\text{c}}}}\\!(M)\,\big{[}\lambda_{{\text{s}}}(M)\big{]}^{2}\,\big{[}{\color[rgb]{0,0,1}u_{\mathrm{s}}(k|M)}\big{]}^{2},$
(71)

where $u_{\mathrm{s}}(k|M)$ is (the Fourier transformation of) the radial
distribution of satellite galaxies within a halo, and we have highlighted the
emulated quantities in blue. In this section, we assume that the distribution
is given by an NFW profile with the concentration-mass relation from Diemer &
Joyce (2019).

The left panel of Fig. 10 compares the model predictions and the galaxy TPCF
measured from mock HOD catalogues at $z=0$. On the scales where the 1-halo
term dominates ($r\lesssim 1\,h^{-1}\mathrm{Mpc}$), the fractional difference
is within $1\%$. The colours in the plot represent different contributions:
the correlations of central-central, central-satellite and satellite-satellite
galaxy pairs.

Figure 10:  Emulator-based halo model predictions of galaxy clustering
statistics. Left panel: Galaxy two-point correlation functions from
simulations (marks) and emulator-based halo model predictions (solid lines),
for the fiducial Planck cosmology of FORGE at $z=0$. Colours represent
different terms of galaxy correlations after the central/satellite split:
galaxy-galaxy (black), central-central (red), central-satellite (orange) and
satellite-satellite (green). Only the one-halo terms of theory predictions are
shown. In the lower sub-panel, we show the relative difference between the
1-halo term and the full correlation function measured from simulations. Right
panel: Similar to the left panel but for the galaxy-matter cross power
spectrum for the same cosmology and HOD prescription. (source code)

As shown in Section 3.2 and Fig. 4, the halo concentrations measured from
simulations are sensitive to the minimum radius in the fitting, with a
relative difference of up to 10 per cent. To test the impact on galaxy
clustering, we calculate the one-halo terms of $\xi_{{\text{gg}}}$ (Eqns. (70)
and (71)), adopting the NFW profile with the concentration-mass relations
measured from this work and increasing or decreasing them by $10$ per cent.
Fig. 11 shows that a $10$ per cent change in the concentration-mass relation
will change the one-halo term by $5$ per cent. The impact on the two-halo term
and the degeneracy between the concentration amplitude and galaxy-halo
connection model parameters will be left for future work.

Figure 11: One-halo terms of the galaxy two-point correlation functions with
the NFW profiles combined with different concentration-mass relations. The
black line shows the fiducial case corresponding to the $c(M)$ relation
measured from this work. The other two coloured lines present the results for
increasing and decreasing the concentration by $10$ per cent. In the lower
sub-panel, we show the relative difference with respect to the fiducial
result. (source code)

###  6.2 Galaxy-matter cross-correlation function

In the galaxy-galaxy weak lensing observations, the excess surface mass
density profile around lensing galaxies, $\Delta\Sigma_{\text{gm}}(R)$ is
measured, which can be expressed in terms of the galaxy-matter cross-power
spectrum as (e.g., Murata et al., 2018; Nishimichi et al., 2019)

$\displaystyle\Delta\Sigma_{\text{gm}}(R)=\bar{\rho}_{{\text{m}}0}\int_{0}^{\infty}\frac{\differential{k}}{2\pi}kP_{\text{gm}}(k)J_{2}(kR),$
(72)

where $J_{2}(x)$ is the second-order Bessel function.

In the halo model framework, we can accurately predict $P_{\text{gm}}(k)$ with
the emulators providing the model ingredients. Under the same configurations
as in the last sub-section, $P_{\text{gm}}(k)$ is related to the halo
properties as

$\displaystyle P_{\text{gm}}(k)$
$\displaystyle=\frac{1}{\bar{n}_{{\text{g}}}}\int\differential{M}{\color[rgb]{0,0,1}\derivative{n(M)}{M}}\expectationvalue{N_{\mathrm{c}}}\\!(M)\times$
$\displaystyle\phantom{\frac{1}{\bar{n}_{{\text{g}}}}\int\differential{M}\quad}\Big{[}1+\lambda_{\mathrm{s}}(M)\,{\color[rgb]{0,0,1}u_{\mathrm{s}}(k|M)}\Big{]}{\color[rgb]{0,0,1}P_{\text{hm}}(k|M)}.$
(73)

$P_{\text{hm}}(k|M)$ is the cross-correlation between the matter overdensity
field with the halo sample in a narrow mass bin $[M,M+\differential{M}]$. The
quantity that our halo-mass cross-correlation emulators output is
$P_{\text{hm}}(k|n_{\text{h}}(>M))$, which can be converted as

$\displaystyle
P_{\text{hm}}(k|M)=-\quantity[\derivative{n_{\text{h}}(M)}{M}]^{-1}\frac{\partial}{\partial
M}\Big{[}n_{\text{h}}(>M)P_{\text{hm}}\quantity(k\big{|}n_{\text{h}}(>M))\Big{]}.$
(74)

We use the publicly available, open-source Python toolkit nbodykit (Hand et
al., 2018) to measure the cross power spectra between the mock galaxy
catalogues and the matter field, in linear $k$ bins from $0.3$ to
$6\,h\,\mathrm{Mpc}^{-1}$ with a width of $\Delta
k=0.02\,h\,\mathrm{Mpc}^{-1}$. These measurements are compared with the halo
model predictions in the right panel of Fig. 10. The relative difference shown
in the lower sub-panel is below 1 per cent except at low-$k$ bins, where the
cosmic variance dominates the error budget.

##  7 Discussions and Conclusions

We present accurate emulators for the halo mass function, concentration-mass
relation and halo-matter cross-correlation function, for $\Lambda$CDM and two
representative modified gravity theories, $f(R)$ gravity and DGP, using the
FORGE and BRIDGE suites of $N$-body simulations (Arnold et al., 2022). The
cosmological parameter space spans three non-MG parameters,
$\Omega_{\mathrm{m0}},h,S_{8}$, and one MG parameter, either $\bar{f}_{R0}$ or
$H_{0}r_{\mathrm{c}}$, depending on which modified gravity model we are using.

We construct emulators using fully connected neural networks implemented using
the open-source Python library PyTorch Lightning. We show the capabilities of
neural networks under noise-free conditions by emulating the existing fitting
formulae of halo properties, such as the fitting function for HMFs in Tinker
et al. (2008). The emulators mimic analytical models in a 7-D parameter space
with sub-percent accuracy, using only $50$ training cosmologies.

For realistic cases where the data come from $N$-body simulations and
therefore have noise, the accuracy of our halo property emulators is
summarised in Figs. 7-9. The emulation error is less than $1\%$ for most
cosmologies in both the training and the test data sets, in the halo mass
range of $10^{12}\leq M_{\text{200c}}/(h^{-1}M_{\odot})\leq 10^{14}$.

The primary purpose of this series of papers is to extend the modelling of
galaxy clustering to non-linear scales. We employ the halo model framework
(Cooray & Sheth, 2002) combined with an adopted galaxy-halo connection model
to predict galaxy clustering and other cosmological observables, following the
spirit of the Dark Quest project (Nishimichi et al., 2019; Kobayashi et al.,
2020; Cuesta-Lazaro et al., 2022). We demonstrate that the emulators can be
applied to the halo model framework combined with the HOD prescription to
predict the one-halo term of the galaxy clustering signal, achieving sub-
percent accuracy. The main advantages of this emulator-based halo model
approach can be summarised as follows.

  * •

Accuracy. The model ingredients provided by emulators incorporate the major
complicated effects in the non-linear regime of structure formation, such as
non-linear halo bias and the halo exclusion effect.

  * •

Versatility. The halo model approach enables a joint modelling of cosmological
observables, such as galaxy-galaxy and galaxy-matter correlation functions,
for a single population of galaxies. The combination of different probes can
mitigate the uncertainty of galaxy formation and evolution on cosmological
parameter inference.

  * •

Flexibility. Instead of making an end-to-end mapping between the cosmological
and HOD parameters to the final clustering statistics with the emulation
process, this “numerical” version of the halo model allows the flexibility of
combining with any specific HOD prescription for different types of galaxies,
without retraining the emulators.

To perform cosmological parameter inference by confronting the emulator-based
halo model prediction with galaxy survey observations, we plan to implement
the following improvements and extensions in the subsequent papers of this
series.

  * •

The excellent performance of the emulators is partly due to the smaller number
of parameters varied in the FORGE and BRIDGE simulations compared with other
emulation projects, as well as the limited halo mass range due to the
relatively small box size of the simulations. However, as indicated by the
ideal emulation tests, neural networks are capable of emulating halo
properties up to the halo masses of $10^{15.5}\,h^{-1}M_{\odot}$ in a higher-
dimensional parameter space with sub-percent accuracy. To extend the mass
range of the emulators, we plan to run additional simulations with different
specifications, such as box size and number of particles, to obtain halo
properties robustly and at low cost.

  * •

Galaxy clustering statistics are typically measured in redshift space from
surveys, which involves not only the information about galaxy positions but
also their peculiar velocities. We will build emulators for the halo peculiar
velocity statistics, such as pairwise velocity moments, and combine them using
a galaxy-halo connection model (e.g. Kobayashi et al., 2020; Cuesta-Lazaro et
al., 2022).

  * •

To resemble actual samples of galaxies, we need more realistic HOD
prescriptions and check the accuracy of the emulator-based halo model.

In the future, we plan to use the neural network emulators on the upcoming
data from DESI and Euclid to constrain the cosmological and modified gravity
parameters. This requires that the models be trained on simulations with
higher resolution to meet the demand of the cutting-edge observational data
with unprecedented volume and much better controlled systematic errors.

## Acknowledgements

C-ZR thanks the Research Council of Norway for their support. C-ZR and BL are
supported by the European Research Council (ERC) through a starting Grant
(ERC-StG-716532 PUNCA). AE is supported at the AIfA by an Argelander
Fellowship. CMB acknowledges support from the Science Technology Facilities
Council through ST/T000244/1. BL and CMB are further supported by the UK
Science and Technology Funding Council (STFC) Consolidated Grant No.
ST/I00162X/1 and ST/P000541/1. CTD is funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy – EXC-2094 – 390783311.

This work used the DiRAC@Durham facility managed by the Institute for
Computational Cosmology on behalf of the STFC DiRAC HPC Facility
(www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC
capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1,
Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the
National e-Infrastructure.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author. The source code and data to generate the figures in this
manuscript are available at this GitHub repository.

## References

  * Agarap (2018) Agarap A. F., 2018, arXiv preprint arXiv:1803.08375
  * Agarwal et al. (2012) Agarwal S., Abdalla F. B., Feldman H. A., Lahav O., Thomas S. A., 2012, MNRAS, 424, 1409
  * Agarwal et al. (2014) Agarwal S., Abdalla F. B., Feldman H. A., Lahav O., Thomas S. A., 2014, MNRAS, 439, 2102
  * Alom et al. (2019) Alom M. Z., et al., 2019, Electronics, 8
  * Amendola et al. (2013) Amendola L., et al., 2013, Living Reviews in Relativity, 16, 6
  * Angulo & Hahn (2022) Angulo R. E., Hahn O., 2022, Living Reviews in Computational Astrophysics, 8, 1
  * Arnold et al. (2019a) Arnold C., Leo M., Li B., 2019a, Nature Astronomy, 3, 945
  * Arnold et al. (2019b) Arnold C., Leo M., Li B., 2019b, Nature Astron., 3, 945
  * Arnold et al. (2022) Arnold C., Li B., Giblin B., Harnois-Déraps J., Cai Y.-C., 2022, MNRAS, 515, 4161
  * Ba et al. (2015) Ba S., Myers W. R., Brenneman W. A., 2015, Technometrics, 57, 479
  * Barreira et al. (2014) Barreira A., Li B., Hellwing W. A., Lombriser L., Baugh C. M., Pascoli S., 2014, J. Cosmology Astropart. Phys., 2014, 029
  * Bhattacharya et al. (2011) Bhattacharya S., Heitmann K., White M., Lukić Z., Wagner C., Habib S., 2011, ApJ, 732, 122
  * Bishop et al. (1995) Bishop C. M., et al., 1995, Neural networks for pattern recognition. Oxford university press
  * Bocquet et al. (2020) Bocquet S., Heitmann K., Habib S., Lawrence E., Uram T., Frontiere N., Pope A., Finkel H., 2020, ApJ, 901, 5
  * Bond et al. (1991) Bond J. R., Cole S., Efstathiou G., Kaiser N., 1991, ApJ, 379, 440
  * Bose et al. (2021) Bose B., et al., 2021, MNRAS, 508, 2479
  * Brando et al. (2022) Brando G., Fiorini B., Koyama K., Winther H. A., 2022, J. Cosmology Astropart. Phys., 2022, 051
  * Carrilho et al. (2022) Carrilho P., Carrion K., Bose B., Pourtsidou A., Hidalgo J. C., Lombriser L., Baldi M., 2022, MNRAS, 512, 3691
  * Cataneo et al. (2016) Cataneo M., Rapetti D., Lombriser L., Li B., 2016, J. Cosmology Astropart. Phys., 2016, 024
  * Cataneo et al. (2019) Cataneo M., Lombriser L., Heymans C., Mead A. J., Barreira A., Bose S., Li B., 2019, MNRAS, 488, 2121
  * Child et al. (2018) Child H. L., Habib S., Heitmann K., Frontiere N., Finkel H., Pope A., Morozov V., 2018, ApJ, 859, 55
  * Cooray & Sheth (2002) Cooray A., Sheth R., 2002, Phys. Rep., 372, 1
  * Courtin et al. (2011) Courtin J., Rasera Y., Alimi J. M., Corasaniti P. S., Boucher V., Füzfa A., 2011, MNRAS, 410, 1911
  * Crocce et al. (2006) Crocce M., Pueblas S., Scoccimarro R., 2006, MNRAS, 373, 369
  * Cuesta-Lazaro et al. (2022) Cuesta-Lazaro C., et al., 2022, arXiv e-prints, p. arXiv:2208.05218
  * Cybenko (1989) Cybenko G., 1989, Mathematics of control, signals and systems, 2, 303
  * DESI Collaboration et al. (2016) DESI Collaboration et al., 2016, arXiv e-prints, p. arXiv:1611.00036
  * De Felice & Tsujikawa (2010) De Felice A., Tsujikawa S., 2010, Living Reviews in Relativity, 13, 3
  * DeRose et al. (2019) DeRose J., et al., 2019, ApJ, 875, 69
  * Diemer (2018) Diemer B., 2018, ApJS, 239, 35
  * Diemer & Joyce (2019) Diemer B., Joyce M., 2019, ApJ, 871, 168
  * Diemer & Kravtsov (2015) Diemer B., Kravtsov A. V., 2015, ApJ, 799, 108
  * Donald-McCann et al. (2022) Donald-McCann J., Beutler F., Koyama K., Karamanis M., 2022, MNRAS, 511, 3768
  * Dooley et al. (2014) Dooley G. A., Griffen B. F., Zukin P., Ji A. P., Vogelsberger M., Hernquist L. E., Frebel A., 2014, ApJ, 786, 50
  * Dvali et al. (2000a) Dvali G., Gabadadze G., Porrati M., 2000a, Physics Letters B, 485, 208
  * Dvali et al. (2000b) Dvali G., Gabadadze G., Porrati M., 2000b, Physics Letters B, 485, 208
  * Fang et al. (2008) Fang W., Wang S., Hu W., Haiman Z., Hui L., May M., 2008, Phys. Rev. D, 78, 103509
  * Ferreira (2019) Ferreira P. G., 2019, ARA&A, 57, 335
  * Gao et al. (2008) Gao L., Navarro J. F., Cole S., Frenk C. S., White S. D. M., Springel V., Jenkins A., Neto A. F., 2008, MNRAS, 387, 536
  * García et al. (2021) García R., Rozo E., Becker M. R., More S., 2021, MNRAS, 505, 1195
  * Goodfellow et al. (2016) Goodfellow I., Bengio Y., Courville A., 2016, Deep learning. MIT press
  * Gupta et al. (2022) Gupta S., Hellwing W. A., Bilicki M., García-Farieta J. E., 2022, Phys. Rev. D, 105, 043538
  * Habib et al. (2007) Habib S., Heitmann K., Higdon D., Nakhleh C., Williams B., 2007, Phys. Rev. D, 76, 083503
  * Hand et al. (2018) Hand N., Feng Y., Beutler F., Li Y., Modi C., Seljak U., Slepian Z., 2018, AJ, 156, 160
  * Harnois-Déraps et al. (2019) Harnois-Déraps J., Giblin B., Joachimi B., 2019, A&A, 631, A160
  * Harnois-Déraps et al. (2022) Harnois-Déraps J., Hernandez-Aguayo C., Cuesta-Lazaro C., Arnold C., Li B., Davies C. T., Cai Y.-C., 2022, arXiv e-prints, p. arXiv:2211.05779
  * Heitmann et al. (2006) Heitmann K., Higdon D., Nakhleh C., Habib S., 2006, ApJ, 646, L1
  * Heitmann et al. (2009) Heitmann K., Higdon D., White M., Habib S., Williams B. J., Lawrence E., Wagner C., 2009, ApJ, 705, 156
  * Heitmann et al. (2010) Heitmann K., White M., Wagner C., Habib S., Higdon D., 2010, ApJ, 715, 104
  * Heitmann et al. (2016) Heitmann K., et al., 2016, ApJ, 820, 108
  * Hendrycks & Gimpel (2016) Hendrycks D., Gimpel K., 2016, arXiv e-prints, p. arXiv:1606.08415
  * Hornik et al. (1989) Hornik K., Stinchcombe M., White H., 1989, Neural networks, 2, 359
  * Hu & Sawicki (2007) Hu W., Sawicki I., 2007, Phys. Rev. D, 76, 064004
  * Hu et al. (2018) Hu B., Liu X.-W., Cai R.-G., 2018, MNRAS, 476, L65
  * Jenkins et al. (2001) Jenkins A., Frenk C. S., White S. D. M., Colberg J. M., Cole S., Evrard A. E., Couchman H. M. P., Yoshida N., 2001, MNRAS, 321, 372
  * Jennings et al. (2019) Jennings W. D., Watkinson C. A., Abdalla F. B., McEwen J. D., 2019, MNRAS, 483, 2907
  * Khoury & Weltman (2004a) Khoury J., Weltman A., 2004a, Phys. Rev. D, 69, 044026
  * Khoury & Weltman (2004b) Khoury J., Weltman A., 2004b, Phys. Rev. Lett., 93, 171104
  * Kingma & Ba (2014) Kingma D. P., Ba J., 2014, arXiv e-prints, p. arXiv:1412.6980
  * Klypin et al. (2016) Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340
  * Kobayashi et al. (2020) Kobayashi Y., Nishimichi T., Takada M., Takahashi R., Osato K., 2020, Phys. Rev. D, 102, 063504
  * Koyama (2007) Koyama K., 2007, Classical and Quantum Gravity, 24, R231
  * Koyama (2018) Koyama K., 2018, International Journal of Modern Physics D, 27, 1848001
  * Kuhlen et al. (2012) Kuhlen M., Vogelsberger M., Angulo R., 2012, Physics of the Dark Universe, 1, 50
  * Kwan et al. (2013) Kwan J., Bhattacharya S., Heitmann K., Habib S., 2013, ApJ, 768, 123
  * Kwan et al. (2015) Kwan J., Heitmann K., Habib S., Padmanabhan N., Lawrence E., Finkel H., Frontiere N., Pope A., 2015, ApJ, 810, 35
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration et al., 2009, arXiv e-prints, p. arXiv:0912.0201
  * Lam & Li (2012) Lam T. Y., Li B., 2012, MNRAS, 426, 3260
  * Laureijs et al. (2011) Laureijs R., et al., 2011, arXiv e-prints, p. arXiv:1110.3193
  * Lawrence et al. (2010) Lawrence E., Heitmann K., White M., Higdon D., Wagner C., Habib S., Williams B., 2010, ApJ, 713, 1322
  * Lawrence et al. (2017) Lawrence E., et al., 2017, ApJ, 847, 50
  * Li & Efstathiou (2012) Li B., Efstathiou G., 2012, MNRAS, 421, 1431
  * Li & Hu (2011) Li Y., Hu W., 2011, Phys. Rev. D, 84, 084033
  * Li et al. (2012) Li B., Zhao G.-B., Teyssier R., Koyama K., 2012, J. Cosmology Astropart. Phys., 2012, 051
  * Li et al. (2013) Li B., Hellwing W. A., Koyama K., Zhao G.-B., Jennings E., Baugh C. M., 2013, MNRAS, 428, 743
  * Linder (2003) Linder E. V., 2003, Phys. Rev. Lett., 90, 091301
  * Lombriser (2014) Lombriser L., 2014, Annalen der Physik, 264, 259
  * Lombriser et al. (2013) Lombriser L., Li B., Koyama K., Zhao G.-B., 2013, Phys. Rev. D, 87, 123511
  * Lombriser et al. (2014) Lombriser L., Koyama K., Li B., 2014, J. Cosmology Astropart. Phys., 2014, 021
  * Ma & Fry (2000) Ma C.-P., Fry J. N., 2000, ApJ, 543, 503
  * Maartens & Koyama (2010) Maartens R., Koyama K., 2010, Living Reviews in Relativity, 13, 5
  * Maksimova et al. (2021) Maksimova N. A., Garrison L. H., Eisenstein D. J., Hadzhiyska B., Bose S., Satterthwaite T. P., 2021, MNRAS, 508, 4017
  * Massara et al. (2014) Massara E., Villaescusa-Navarro F., Viel M., 2014, J. Cosmology Astropart. Phys., 2014, 053
  * McClintock et al. (2019) McClintock T., et al., 2019, ApJ, 872, 53
  * Mead et al. (2021) Mead A. J., Brieden S., Tröster T., Heymans C., 2021, MNRAS, 502, 1401
  * Meneghetti & Rasia (2013) Meneghetti M., Rasia E., 2013, arXiv e-prints, p. arXiv:1303.6158
  * Mitchell et al. (2019) Mitchell M. A., Arnold C., He J.-h., Li B., 2019, MNRAS, 487, 1410
  * Mitchell et al. (2021) Mitchell M. A., Hernández-Aguayo C., Arnold C., Li B., 2021, MNRAS, 508, 4140
  * Miyatake et al. (2020) Miyatake H., et al., 2020, arXiv e-prints, p. arXiv:2101.00113
  * Murata et al. (2018) Murata R., Nishimichi T., Takada M., Miyatake H., Shirasaki M., More S., Takahashi R., Osato K., 2018, ApJ, 854, 120
  * Murray et al. (2013) Murray S. G., Power C., Robotham A. S. G., 2013, Astronomy and Computing, 3, 23
  * Murray et al. (2021) Murray S. G., Diemer B., Chen Z., Neuhold A. G., Schnapp M. A., Peruzzi T., Blevins D., Engelman T., 2021, Astronomy and Computing, 36, 100487
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Neto et al. (2007) Neto A. F., et al., 2007, MNRAS, 381, 1450
  * Neyman & Scott (1952) Neyman J., Scott E. L., 1952, ApJ, 116, 144
  * Nishimichi et al. (2019) Nishimichi T., et al., 2019, ApJ, 884, 29
  * Peacock & Smith (2000) Peacock J. A., Smith R. E., 2000, MNRAS, 318, 1144
  * Peebles (1980) Peebles P. J. E., 1980, The large-scale structure of the universe
  * Philcox et al. (2020) Philcox O. H. E., Spergel D. N., Villaescusa-Navarro F., 2020, Phys. Rev. D, 101, 123520
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Prada et al. (2012) Prada F., Klypin A. A., Cuesta A. J., Betancort-Rijo J. E., Primack J., 2012, MNRAS, 423, 3018
  * Press & Schechter (1974) Press W. H., Schechter P., 1974, ApJ, 187, 425
  * Ramachandra et al. (2021) Ramachandra N., Valogiannis G., Ishak M., Heitmann K., LSST Dark Energy Science Collaboration 2021, Phys. Rev. D, 103, 123525
  * Ruan et al. (2022) Ruan C.-Z., Hernández-Aguayo C., Li B., Arnold C., Baugh C. M., Klypin A., Prada F., 2022, J. Cosmology Astropart. Phys., 2022, 018
  * Sahni & Shtanov (2003) Sahni V., Shtanov Y., 2003, J. Cosmology Astropart. Phys., 2003, 014
  * Schmidt (2016) Schmidt F., 2016, Phys. Rev. D, 93, 063512
  * Schmidt et al. (2009) Schmidt F., Lima M., Oyaizu H., Hu W., 2009, Phys. Rev. D, 79, 083518
  * Schmidt et al. (2010) Schmidt F., Hu W., Lima M., 2010, Phys. Rev. D, 81, 063005
  * Seljak (2000) Seljak U., 2000, MNRAS, 318, 203
  * Sheth et al. (2001) Sheth R. K., Mo H. J., Tormen G., 2001, MNRAS, 323, 1
  * Sinha & Garrison (2020) Sinha M., Garrison L. H., 2020, MNRAS, 491, 3022
  * Smith & Angulo (2019) Smith R. E., Angulo R. E., 2019, MNRAS, 486, 1448
  * Smith et al. (2003) Smith R. E., et al., 2003, MNRAS, 341, 1311
  * Song et al. (2007) Song Y.-S., Sawicki I., Hu W., 2007, Phys. Rev. D, 75, 064003
  * Sotiriou & Faraoni (2010) Sotiriou T. P., Faraoni V., 2010, Reviews of Modern Physics, 82, 451
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
  * Takahashi et al. (2012) Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, ApJ, 761, 152
  * Takahashi et al. (2020) Takahashi R., Nishimichi T., Namikawa T., Taruya A., Kayo I., Osato K., Kobayashi Y., Shirasaki M., 2020, ApJ, 895, 113
  * Tinker et al. (2008) Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709
  * Troja et al. (2022) Troja A., Tutusaus I., Sorce J. G., 2022, arXiv e-prints, p. arXiv:2211.09668
  * Vainshtein (1972) Vainshtein A. I., 1972, Physics Letters B, 39, 393
  * Wang et al. (2020) Wang J., Bose S., Frenk C. S., Gao L., Jenkins A., Springel V., White S. D. M., 2020, Nature, 585, 39
  * Weinberger et al. (2020) Weinberger R., Springel V., Pakmor R., 2020, ApJS, 248, 32
  * Williams & Rasmussen (2006) Williams C. K., Rasmussen C. E., 2006, Gaussian processes for machine learning. MIT press Cambridge, MA
  * Winther et al. (2019) Winther H. A., Casas S., Baldi M., Koyama K., Li B., Lombriser L., Zhao G.-B., 2019, Phys. Rev. D, 100, 123540
  * Yuan et al. (2022) Yuan S., Garrison L. H., Eisenstein D. J., Wechsler R. H., 2022, MNRAS, 515, 871
  * Zhai et al. (2019) Zhai Z., et al., 2019, ApJ, 874, 95
  * Zhang et al. (2021) Zhang A., Lipton Z. C., Li M., Smola A. J., 2021, arXiv e-prints, p. arXiv:2106.11342
  * Zhao et al. (2011) Zhao G.-B., Li B., Koyama K., 2011, Phys. Rev. D, 83, 044007
  * Zheng et al. (2005) Zheng Z., et al., 2005, ApJ, 633, 791
  * Zheng et al. (2007) Zheng Z., Coil A. L., Zehavi I., 2007, ApJ, 667, 760

  *[DGP]: Dvali-Gabadadze-Porrati
  *[CMB]: cosmic microwave background
  *[ReLU]: Rectified Linear Unit
  *[GELU]: Gaussian error linear unit
  *[TPCF]: two-point correlation function