# The BACCO simulation project: biased tracers in real space

Matteo Zennaro,1 Raul E. Angulo,1,2 Marcos Pellejero-Ibáñez,1 Jens Stücker,1
Sergio Contreras,1 and Giovanni Aricò1,3
1Donostia International Physics Center (DIPC), Paseo Manuel de Lardizabal, 4,
20018, Donostia-San Sebastián, Guipuzkoa, Spain.
2IKERBASQUE, Basque Foundation for Science, 48013, Bilbao, Spain.
3Universidad de Zaragoza, Pedro Cerbuna 12, 50009 Zaragoza, Spain
E-mail:matteo_zennaro001<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present an emulator for the two-point clustering of biased tracers in real
space. We construct this emulator using neural networks calibrated with more
than $400$ cosmological models in a 8-dimensional cosmological parameter space
that includes massive neutrinos an dynamical dark energy. The properties of
biased tracers are described via a Lagrangian perturbative bias expansion
which is advected to Eulerian space using the displacement field of numerical
simulations. The cosmology-dependence is captured thanks to a cosmology-
rescaling algorithm. We show that our emulator is capable of describing the
power spectrum of galaxy formation simulations for a sample mimicking that of
a typical Emission-Line survey at $z\sim 1$ with an accuracy of $1-2\%$ up to
nonlinear scales $k\sim 0.7h\,{\rm Mpc}^{-1}$.

###### keywords:

cosmology: theory – large-scale structure of Universe – methods: statistical –
methods: computational

††pubyear: 2020††pagerange: The BACCO simulation project: biased tracers in
real space–B

##  1 Introduction

The observed spatial distribution of galaxies and quasars offers an extremely
valuable window to the physics of the universe. For instance, their clustering
as a function of scale depends on the early universe physics, where properties
of the primordial fluctuation as well as the parameters of the cosmological
model leave distinctive signatures. Similarly, the relation between cosmic
velocities and densities, encoded in the so-called redshift space distortions,
offers a pathway to constrain the nature of gravity and the growth of
structure (e.g. Kaiser, 1987; Guzzo et al., 2008). Finally, baryonic acoustic
oscillations (e.g. Eisenstein & Hu, 1998) can be employed as a standard ruler
to measure the expansion history of the universe, offering an opportunity to
better constrain the properties of dark energy (see e.g. Dodelson & Schneider,
2013; Taylor et al., 2013; Paz & Sánchez, 2015).

There is a large number of ongoing observational campaigns that will take on
the wealth of information encoded in clustering (e.g. Euclid, see the review
by Amendola et al. 2018, DESI, Levi et al. 2013; DESI Collaboration et al.
2016, and J-PAS Bonoli et al. 2020). These observations will map the position
and shapes of hundreds of millions of galaxies up to tens of thousands of
square degrees (e.g. Bartelmann & Schneider, 2001; Takada et al., 2014; Ivezić
et al., 2019; Dore et al., 2019; Mandelbaum, 2018). Although statistical
uncertainties and observational systematic errors will be crucial to keep
under control, the actual limitation in exploiting the surveys will arise from
the precision with which the observed distribution of the luminous tracers can
be modelled.

The challenge of modelling galaxy clustering arises from the nonlinearity of
the physics involved. Although the early seeds of structure formation can be
predicted accurately and efficiently using linearised Bolztmann-Einstein
equations (Lewis et al., 2000; Lesgourgues, 2011), the subsequent
gravitational evolution will create nonlinearities which will dominate on
small scales. Furthermore, galaxies and quasars are expected to form in
particular regions of the universe with an efficiency that depends on the
details of the galaxy formation physics involved.

The most accurate way to follow all these processes is provided by numerical
simulations (see, e.g. Kuhlen et al., 2012, for a review). In these,
fluctuations in the early universe are represented by a set of $N$-body
particles which are then evolved solving their equations of motion. Recent
advances in the field have different codes and numerical approaches to
converge in the predictions for the nonlinear matter power spectrum to better
than 2% down to scales of $k\sim 10h\,{\rm Mpc}^{-1}$ (Schneider et al., 2016;
Springel et al., 2020; Garrison et al., 2018; Angulo et al., 2020). On the
other hand, simulations are computationally expensive, thus it is typically
not possible to carry them out for more than a handful of different choices of
cosmological parameters. Nevertheless, several approaches have been suggested
in the literature to speed up their production, (e.g. Monaco et al., 2002;
Tassev et al., 2013; Izard et al., 2016), or to interpolate among the outputs
of simulations (Heitmann et al., 2014; Liu et al., 2018; Nishimichi et al.,
2019; DeRose et al., 2019; Giblin et al., 2019; Euclid Collaboration et al.,
2019; Wibking et al., 2019; Winther et al., 2019; Angulo et al., 2020; Euclid
Collaboration et al., 2020).

Another problem that arises from a potential modelling from numerical
simulation refers to galaxy formation. Although there has been great progress
in understanding the formation of galaxies, correlations with halo properties
and the impact of various processes, predictions for galaxy properties are
still uncertain and model dependent. This has the drawback that, when compared
to observations, small uncertainties in galaxy modelling could heavily bias
cosmological constraints. An alternative could be marginalization over galaxy
formation. Though this is in principle possible, it would add much higher
computational demands, and a poor or incomplete modelling of galaxy properties
could still lead to biased cosmological inferences.

Perturbation theory offers a very attractive alternative. By solving
analytically the relevant fluid equations, predictions for the distribution of
matter down to quasilinear scales can be obtained efficiently (e.g. Bernardeau
et al., 2002) and with control of theoretical uncertainties (see e.g.
Chudaykin et al., 2020). The relation between biased tracers and the
underlying matter field can also be treated perturbatively (McDonald, 2009;
Desjacques et al., 2018a; Fujita & Vlah, 2020). By including all possible
dependences to a given order allowed by symmetries, it is possible to extend
the predictions for any biased tracers, without making an explicit connection
with particular galaxy formation physics. Recent advances in these fields have
significantly improved the accuracy of these predictions, which can typically
reach scales of $k\sim 0.2h\,{\rm Mpc}^{-1}$ (Baumann et al. 2012, Baldauf et
al. 2016, Vlah et al. 2016, Ivanov et al. 2020, d'Amico et al. 2020, Colas et
al. 2020, Nishimichi et al. 2020, Chen et al. 2020 ). Although this is a
remarkable achievement, these predictions still fall short with respect to the
observations and accuracy of future galaxy surveys (Blas et al. 2014, McQuinn
& White 2016). This approach has been extensively used under the assumption of
perturbed dynamics in cosmological parameter estimation (see e.g. Chuang et
al. 2017 and Pellejero-Ibanez et al. 2017).

In this paper we combine both approaches – numerical simulations and
perturbation theory – to create a framework that inherits the accuracy of
numerical simulations with the flexibility of a perturbative bias expansion,
which can reach even small nonlinear scales while being agnostic to galaxy
formation physics and details of any given particular observational surveys.
As shown by Modi et al. (2020), this approach has the potential to accurately
describe galaxy clustering down to smaller scales than those typically reached
by perturbation theory.

We are able to do this by combining the benefits of several recent
developments. On the simulation side, we employ large $N$-body simulations,
with accurate force and mass resolution, which have significant noise
suppression by using special initial conditions, and have been carefully
designed to be used in combinations with cosmology-rescaling algorithms. These
allow to densely cover a target parameter space with the equivalent of
hundreds to thousands of simulations. We combine these with feed forward
neural networks which allow us to quickly predict non-linear fields while
varying any cosmological parameter – including neutrinos and dynamical dark
energy – within defined regions in cosmological parameter space. On the
perturbative side, we employ a Lagrangian bias expansion up to second order
which captures dependences with the local density and tidal fields, and that
includes a higher-order derivative bias parameters, capturing the non-locality
of the galaxy formation process. Using the non-linear displacement field of
the rescaled simulations the Lagrangian bias descriptions are advected to
Eulerian space where the different fields can be combined to make predictions
for the non-linear power spectrum of biased tracers.

This paper is structured as follows. In §2 we describe the variety of
numerical tools we employ, including numerical simulations and the power
spectrum calculation, and we validate the cosmology-rescaling approach. We
also describe the implementation of the Lagrangian bias expansion. In §3 we
describe a perturbative solution for the matter fields we consider, which we
will employ to complement our numerical approach. In §4 we provide details of
our emulation via neural networks, including its validation and estimation of
its accuracy. In §5 we present an application of our emulator by fitting the
power spectrum of galaxies mimicking a Star Formation Rate selected (SFR)
sample at $z=1$. We conclude in §6.

##  2 Numerical Methods

Here we present the main numerical methods underlying this work. Specifically,
we describe our simulations in §2.1 and how we model biased tracers in §2.2.
We finish by validating the cosmology-rescaling approach for our purposes in
§2.3.

###  2.1 Simulations

We will employ two sets of simulations. The first one, referred to as the
BACCO simulations will be the core of our biased-tracers emulators. The second
suite will be employed to test the accuracy of our predictions.

Cosmology | $\Omega_{\rm cdm}$ | $\Omega_{\rm b}$ | $h$ | $n_{\rm s}$
---|---|---|---|---
Nenya | 0.265 | 0.050 | 0.60 | 1.01
Narya | 0.310 | 0.050 | 0.70 | 1.01
Vilya | 0.210 | 0.060 | 0.65 | 0.92
TheOne | 0.259 | 0.048 | 0.68 | 0.96
Table 1: Cosmological parameters of the four cosmologies simulated in the
BACCO project. All the cosmologies assume a flat geometry, no massive
neutrinos ($M_{\nu}=0$ eV), a dark energy equation of state with $w_{0}=-1$
and $w_{a}=0$, an amplitude of cold matter fluctuations $\sigma_{8}=0.9$, and
optical depth at recombination $\tau=0.0952$.

####  2.1.1 BACCO Simulations

Our main suite of simulations corresponds to the core of the BACCO simulation
projects. These corresponds to 8 large $N$-body simulations specially designed
to cover a wide range of cosmological parameters when combined with cosmology-
rescaling. These simulations were presented in Angulo et al. (2020), we here
simply provide a recap of their main characteristics.

Specifically, the BACCO simulations are a set of gravity-only simulations with
a box of sidelength $L=1440\,h^{-1}{\rm Mpc}$ resolved with $4320^{3}$
particles, which implies a mass resolution of $m_{p}\sim 3\times
10^{9}\,h^{-1}{\rm M_{\odot}}$. These simulations are carried out in pairs at
4 distinct sets of cosmological parameters. The parameter values, provided in
Table 1, were chosen so that, in combination, these simulations can
efficiently cover a region of approximately $10\sigma$ around the best fit
values obtained by the analysis of the Planck satellite. We refer to Contreras
et al. (2020a) for details on how these parameters were chosen. For each of
these cosmologies we carried out two simulations whose initial Fourier
amplitudes have been fixed but their initial phases are inverted. As shown by
Angulo & Pontzen (2016), this configuration allows for a dramatic reduction of
the noise in the power spectrum due to cosmic variance – by more than two
orders of magnitude for $k<0.1h\,{\rm Mpc}^{-1}$.

The gravitational evolution was carried out with an updated version of
L-Gadget3 (Springel, 2005; Angulo et al., 2020), and employing a Plummer-
equivalent softening length of $5h^{-1}{\rm kpc}$. Numerical parameters were
chosen so that the power spectrum displays a convergence better than 2% at
$k\sim 10\,h\,{\rm Mpc}^{-1}$. In fact, as shown in Angulo et al. (2020), such
configurations agree to better than 2% with other state-of-the-art codes in a
realization of the "Euclid Simulation Challenge" (Schneider et al., 2016).

####  2.1.2 Test Suite of Simulations

A key aspect of our framework is the ability of employing the 4 BACCO
cosmologies to sample hundreds of different cosmologies using a cosmology-
rescaling.

To test the accuracy of such rescaling, we will employ a suite of $35$
independent $N$-body simulations. Each of these simulations consists of
$1536^{3}$ particles in a cubic volume of $512h^{-1}{\rm Mpc}$ a side.
Naturally, these simulations feature a much smaller volume than our main BACCO
simulations, however they have identical numerical parameters, mass
resolution, and force resolution. As for our main suite, these were carried
out with the latest version of L-Gadget3 and were initialised in pairs using
the approach of Angulo & Pontzen (2016) and employing 2LPT at $z=49$.

To further improve the accuracy of our comparison, we have carried out a
version of the BACCO simulations but with a volume and initial phase field
matching those of each of our test simulations. In this way, we reduce
significantly the role of cosmic variance and allow for an accurate testing.

The cosmologies of our test suite vary systematically one of eight
cosmological parameters $\boldsymbol{\vartheta}=\\{\Omega_{\rm m},\Omega_{\rm
b},\sigma_{8},n_{s},h,M_{\nu},w_{0},w_{a}\\}$, over the same range in which we
will build our emulator. When testing the accuracy of our predictions we will
compare the results of rescaling these BACCO simulations against simulations
carried out directly with the target cosmology.

Figure 1: Visualization of the different Lagrangian fields that have been
advected to Eulerian space. Each panel corresponds to a projection of the same
$100\times 281\times 25\,h^{-1}{\rm Mpc}$ volume. The first panel corresponds
to a uniform weighting – thereby recovering the dark matter density field,
whereas the other panels are weighted by different components of the
Lagrangian linear density field. One can see how the different Lagrangian bias
components emphasize different aspects of the large-scale structure – for
example the linear density weighting seems to emphasize filaments and clusters
in the cosmic web and the density squared weighting emphasizes the regions
associated with the most massive clusters (small red points).  Figure 2: Ratio
showing the comparison of the power spectrum of linear lagrangian fields
advected to Eulerian coordinates as predicted by $N$-body simulation and by
cosmology-rescaled simulations. Each panel displays the results for the cross-
spectra of different linear fields, $P_{ij}$, as indicated by the legend;
whereas lines of different colours shows different cosmological models in our
test suite, which includes dynamical dark energy an massive neutrinos.

###  2.2 Modelling biased tracers

To model biased tracers, we will consider a Lagrangian bias expansion. We
refer to Desjacques et al. (2018b) for a review of the perturbative bias
formalism. As discussed in the introduction, the use of a general biasing
formalism will allow us to model the clustering of any tracer of the
underlying matter field, without making any strong assumptions about galaxy
formation physics.

Specifically, we will describe the overdensity of objects in Lagrangian space,
$\delta_{\rm g}(\boldsymbol{q})$, as a second-order expansion in the linear
matter overdensity $\delta(\boldsymbol{q})$ and include potentially non-local
dependences via a higher-order derivative $\nabla^{2}\delta(\boldsymbol{q})$.
Hence:

$\begin{split}\delta_{\rm
g}(\boldsymbol{q})&=1+b_{1}\delta(\boldsymbol{q})+b_{2}\left[\delta^{2}(\boldsymbol{q})-\left\langle\delta^{2}\right\rangle\right]\\\
&+b_{s^{2}}\left[s^{2}(\boldsymbol{q})-\left\langle
s^{2}\right\rangle\right]+b_{\nabla^{2}\delta}\nabla^{2}\delta(\boldsymbol{q}),\end{split}$
(1)

where $1$ is a homogeneous field;
$s^{2}(\boldsymbol{q})=s_{ij}(\boldsymbol{q})s_{ij}(\boldsymbol{q})$ is the
shear field, defined by the tidal tensor
$s_{ij}(\boldsymbol{q})=\partial_{i}\partial_{j}\Phi(\boldsymbol{q})-\delta_{ij}\delta(\boldsymbol{q})$,
with $\Phi(\boldsymbol{q})$ the local matter gravitational potential. Here,
$b_{1},b_{2},b_{s^{2}},$ and $b_{\nabla^{2}\delta}$ are the Lagrangian bias
parameters, which we assume to be scale-independent in Lagrangian space. Note
that it is possible to include higher order bias terms and derivatives (e.g.
Fujita et al., 2020), albeit they are expected to have a small contribution
for low and intermediate mass haloes (Abidi & Baldauf 2018, Lazeyras & Schmidt
2019).

To describe the tracer overdensity field in Eulerian space, $\delta_{\rm
g}(\boldsymbol{x})$, the coordinate system needs to be advected
$\boldsymbol{x}=\boldsymbol{q}+\boldsymbol{\Psi}(\boldsymbol{q})$, where
$\boldsymbol{\Psi}(\boldsymbol{q})$ is referred to as the displacement field.
In terms of the advected fields $\delta_{i}(\boldsymbol{x})$, the power
spectrum can then be expressed in Eulerian space as

$P_{\rm
gg}=\sum_{i,j\in\\{1,\delta,\delta^{2},\nabla^{2}\delta\\}}(2-\delta_{ij})b_{i}b_{j}\,\langle|\delta_{i}(\boldsymbol{k})\delta_{j}^{*}(\boldsymbol{k})|\rangle,\\\
$ (2)

where $\delta_{i}(\boldsymbol{k})$ corresponds to the Fourier transform of the
advected field $\delta_{i}(\boldsymbol{x})$, and it is completely determined
by 15 cross-spectra,
$P_{ij}\equiv\langle|\delta_{i}(\boldsymbol{k})\delta_{j}^{*}(\boldsymbol{k})|\rangle$.
Usually, these spectra are computed perturbatively up to, at least, the same
order as the bias expansion. This has the advantage that predictions can be
computed quickly and accurately as a function of cosmology, but only on
relatively large scales.

In this work, we follow a different strategy and directly compute the 15
cross-spectra relevant for the 2nd-order bias expansion using the results of
numerical simulations. As shown by Modi et al. (2020), this approach improves
significantly the reach of scales accessible to analytic bias expressions and
specifically, was able to accurately describe the clustering of mock HOD
galaxies down to scales of $0.6h\,{\rm Mpc}^{-1}$.

Operationally, to implement this approach we create the 5 relevant fields –
$1$, $\delta(q)$, $\delta^{2}$, $s^{2}(q)$, and $\nabla^{2}\delta(q)$ – in
Lagrangian coordinates using a grid of $1080^{3}$ points, using Fourier
amplitudes and phases matching those of our $N$-body simulations. We then
apply a Gaussian smoothing of size $0.75h^{-1}{\rm Mpc}$, mimicking possible
exclusion effects and the non-locality of structure formation. We then compute
a non-linear displacement field on the same Lagrangian field by following
simulation particles initially located at the same grid points. Finally, we
compute the corresponding advected fields by using a cloud-in-cell assignment
in Eulerian space where the position is only given by the cosmology-dependent
displacement field and where the weight is given by the value of the
respective bias expansion field.

In Figure 1 we show a slice through each of these 5 fields in Eulerian space:
the ‘homogeneous’ is the field weighted by a constant $1$, the ‘density’ is
weighted by $\delta(\boldsymbol{q})$, the ‘quadratic density’ field is
weighted by
$\left[\delta^{2}(\boldsymbol{q})-\left\langle\delta^{2}\right\rangle\right]$,
the ‘tidal field’ by $\left[s^{2}(\boldsymbol{q})-\left\langle
s^{2}\right\rangle\right]$, and the ‘Laplacian’ by
$\nabla^{2}\delta(\boldsymbol{q})$.

A key aspect of our work is that the displacement field connecting Lagrangian
and Eulerian spaces can be easily and accurately evaluated for different
cosmologies with a cosmology rescaling algorithm. We will explore this in the
following subsection.

Figure 3:  Comparison of the power spectra of linear lagrangian field advected
to Eulerian coordinates at $z=0$, as predicted by analytic LPT calculations,
$P_{\rm LPT}$, and various kinds of simulations. Each panel shows results for
a different cross-spectra as indicated by the legend. In each panel coloured
lines show the results of using an ensemble of $1000$ realisations of
$L=1500h^{-1}{\rm Mpc}$ advected using 1st, 2nd, or 3rd-order LPT
displacements; whereas red symbols show the measurements in of our of $N$-body
simulations.

###  2.3 Cosmology Rescaling

To capture the cosmology dependence of the Lagrangian fields in Eulerian
coordinates, we will employ the so-called cosmology-rescaling algorithms.

The main idea of such approach is that the nonlinear structure in a given
cosmology can be mimicked by rescaling the outputs of a simulation carried out
in a nearby cosmology. The rescaling is performed by transforming the length
unit and by considering a different output time such that the linear variance
of the field as a function of scale matches that of the rescaled simulation.
Additionally, large-scale modes are modified by additionally displacing a
given simulation particle with the difference of the 2LPT displacements in the
original and target cosmologies. Velocities are modified in an analogous
manner.

Cosmology-rescaling was originally proposed by Angulo & White (2010) and has
been extensively tested in multiple subsequent studies. In particular, in
Zennaro et al. (2019) we extended the approach to massive neutrinos; in
Contreras et al. (2020a) we showed that the matter power spectrum of dark
matter, haloes, and subhaloes can be retrieved to better than 3% up to $k\sim
5h\,{\rm Mpc}^{-1}$; and in Ondaro et al (in prep.) that the halo mass
function can be obtained to better than 2% over the range
$10^{12}-10^{15}\,h^{-1}{\rm M_{\odot}}$. Additionally, this approach was
employed by Angulo et al. (2020) to construct an emulator for the nonlinear
matter power spectrum, and by Aricò et al. (2020) to incorporate the effects
of baryonic physics such as star formation, gas cooling and feedback.

Here, we will quantify the performance of cosmology-rescaling in predicting
the cross-spectra of our Eulerian linear fields. For this, we will compare the
measurements in our suite of test simulations (c.f. §2.1) with those obtained
by rescaling one of our BACCO simulations.

We present our results in Fig. 2 for the 15 possible combinations of our 5
Lagrangian fields at $z=0$, as indicated by the legend in each panel. Each of
the lines displays the results for one of our test cosmologies. We recall that
these cover a parameter space roughly set by a $10\sigma$ region around
Planck’s best fit values. Firstly, we can see that cosmology rescaling
retrieves highly-accurate predictions, agreeing with those of direct $N$-body
simulations to better than 1-2% in most cases. Although not shown here, we
have checked that the precision of the method is even slightly better at
higher redshifts.

The large-scale clustering is expected to be dominated by the terms involving
the linear density and the homogeneous field – $\langle 11\rangle$, $\langle
1\delta\rangle$, and $\langle\delta\delta\rangle$ – these combinations are
particularly well recovered, to better than 1%, even on the smallest scales we
consider. On smaller scales, we expect fields involving $\delta^{2}$ – the
combinations $1\delta^{2}$, $\delta^{2}\delta^{2}$, and $\delta^{2}$ – to also
have relevant contributions. Similarly, for such fields, we also obtain highly
accurate predictions, with deviations being typically below $\sim 2\%$ at
$k\sim 1h\,{\rm Mpc}^{-1}$. Note that these fields have particularly low
amplitudes on large scales, which is why the comparison becomes progressively
noisier as we consider smaller wavenumbers.

Finally, we expect fields involving the tidal field – $1s^{2}$, $\delta
s^{2}$, $s^{2}s^{2}$ – to have subdominant contributions on all scales.
Nevertheless, our predictions are still relatively accurate with typical
discrepancies of less than 5%. The least accurate predictions among all the
cases where we can reliably perform the comparison is $\delta^{2}s^{2}$, where
the cosmology-rescaling can over- or under-predict the results of our test
simulations by $5-10\%$.

We note that, even though these simulations have matching Fourier phases, some
cross-spectra are particularly noisy on large scales. As we will see in the
next sections, these are inherently noisy. These fields also display no
benefit from the “Paired-&-Fixed” method.

In summary, we have shown that, using a cosmology-rescaling, it is possible to
accurately predict all the 15 cross-spectra involved in the perturbative bias
expansion. This is remarkable as it opens up the possibility to compute
predictions for a densely-sampled cosmological parameter space which should
yield to an accurate emulation. Nevertheless, numerical predictions for some
combination of fields are particularly noisy on large-scales (e.g.
$\delta^{2}\nabla^{2}$, which is zero at linear order even on large scales),
which could pose a problem for a direct emulation. For this reason we will
combine our numerical results with analytic calculations using Lagrangian
perturbation theory, which will be the subject of our next section.

##  3 Lagrangian Perturbation Theory

In the previous section, we have introduced the Lagrangian bias expansion, and
in particular, the 15 cross-spectra of the different Lagrangian fields that
enter the biased-tracers power spectrum. In this section, we model these 15
terms theoretically, with the aim of using these LPT predictions to reduce the
noise present in the measurements from rescaled simulations, and provide us
with a way of reducing the dynamical range of the quantities we want to
emulate.

###  3.1 LPT predictions

To compare with quantities measured in simulations, we want each of the
components of the field presented in Eq. (1) to be evaluated at a given
Eulerian position $\boldsymbol{x}$. The field at that position will receive
contributions from all Lagrangian positions that, after being displaced, end
up in the same Eulerian position (e.g. Matsubara, 2008). Therefore,

$1+\delta_{F}(\boldsymbol{x})=\int\mathrm{d}^{3}\boldsymbol{q}F(\boldsymbol{q})\delta_{\rm
D}\left(\boldsymbol{x}-\boldsymbol{q}-\boldsymbol{\Psi}(\boldsymbol{q})\right),$
(3)

where $F(\boldsymbol{q})$ can be any of the components of the field.

We obtain the overdensities created in Eulerian space from displacing a
uniform distribution $F(\boldsymbol{q})=1$, a linear density field
$F(\boldsymbol{q})=\delta_{\rm L}(\boldsymbol{q})$, its square
$F(\boldsymbol{q})=\delta^{2}_{\rm L}(\boldsymbol{q})$, the tidal field
$F(\boldsymbol{q})=s^{2}(\boldsymbol{q})$ and the laplacian of the linear
density field $F(\boldsymbol{q})=\nabla^{2}\delta_{\rm L}(\boldsymbol{q})$.
Note that we expand the displacement field
$\boldsymbol{\Psi}\simeq\boldsymbol{\Psi}^{(1)}+\boldsymbol{\Psi}^{(2)}+\dots$,
retaining only its first order expansion up to second order in power, so that
in Fourier space

$\begin{split}\delta_{F}(\boldsymbol{k})=\int\mathrm{d}^{3}\boldsymbol{q}e^{-i\boldsymbol{k}\cdot\boldsymbol{q}}F(\boldsymbol{q})\left\\{1-\boldsymbol{k}\cdot\boldsymbol{\Psi}^{(1)}(\boldsymbol{q})-\dfrac{1}{2}[\boldsymbol{k}\cdot\boldsymbol{\Psi}^{(1)}(\boldsymbol{q})]^{2}\right\\},\end{split}$
(4)

ignoring an additional Dirac’s delta $\delta_{D}(\boldsymbol{k})$ that would
only contribute to the wavemode $\boldsymbol{k}=0$. Note that the first order
displacement field can be easily connected to the linear overdensity in
Fourier space, as

$\boldsymbol{\Psi}^{(1)}(\boldsymbol{k})=-i\dfrac{\boldsymbol{k}}{k^{2}}\delta(\boldsymbol{k}).$
(5)

After advecting these fields with Eq. (4) we combine them into the respective
cross spectra, retaining only terms of order $(11)$ and $(22)$, obtaining
contributions of the form

$\begin{split}P_{ij}(k)&=A_{ij}(k)P(k)\\\
&+B_{ij}(k)\int\dfrac{\mathrm{d}^{3}\boldsymbol{p}}{(2\pi)^{3}}P(\boldsymbol{p})P(|\boldsymbol{k}-\boldsymbol{p}|)C_{ij}(\boldsymbol{p},\boldsymbol{k}),\end{split}$
(6)

where the prefactors $A_{ij}(k)$ and $B_{ij}(k)$, and the kernels
$C_{ij}(\boldsymbol{p},\boldsymbol{k})$ are different for each combination of
the fields are are reported in appendix A.

Numerically, we compute all LTP integrals using a doubly-adaptive quadrature
integration, as implemented in the Gnu Scientific Library, whose accuracy we
have tested against a multi-dimensional numerical integration. Each integral
takes approximately $0.2$ seconds. In Appendix B we compare our calculations
against two publicly-available code FastPT (McEwen et al., 2016; Fang et al.,
2017) and velocileptors (Chen et al., 2020).

###  3.2 LPT validation

To validate our LPT expressions and their numerical evaluation, we will
compare them against a suite of realizations of an initially homogeneous
distribution advected to Eulerian space using Lagrangian perturbation theory
at different orders.

Specifically, we construct a suite of $1000$ realizations of the relevant
Lagrangian fields on $384^{3}$ grids in $L=1500h^{-1}{\rm Mpc}$ boxes, and
using a cosmology compatible with current observational constraints. These
grids are then advected to $z=0$ using either first, second, or third-order
Lagrangian perturbation theory, following the implementation of Michaux et al.
(2021).

We then measure the resulting cross-spectra using FFTs with $384^{3}$ grid
points. We show our results in Fig. 3 where we display the measured spectra
over our LPT analytic calculations. For comparison, we also display the
corresponding measurements from one of our main BACCO simulation as red
symbols.

Firstly, we can see that for most 15 power spectra, our analytic calculations
and the 1LPT simulations agree remarkably well on large-scale – the ratio
approaches to unity for low wavenumers. This is an important validation of our
calculations. Furthermore, for terms involving, $1$ and $\delta$, there is
also agreement with higher-order version of LPT and with the results of our
full $N$-body simulation.

In other cases, there is a disagreement between 1 and 2LPT simulations with
the later agreeing with the $N$-body results. This suggest that to predict
accurately these quantities, higher-orders in the displacement need to be
included. Yet in other instances, mostly those involving $\delta^{2}$, the LPT
simulations appear to converge to a different large-scale limit than the
$N$-body result, which points to shell-crossing and nonlinearities on small
scales being important to correctly predict the amplitude on large-scales. We
leave further investigation of this for future work.

Nevertheless, the discrepancies on large scales appear in fields which are
already highly subdominant on large scales. Thus we conclude that our analytic
LPT calculations are sufficiently accurate to complement our numerical
measurements on large scales.

##  4 Emulating biased tracers

Figure 4:  The cross spectra of various Lagrangian fields advected to Eulerian
space, $P_{ij}$ where $i,j=\\{1,\delta,\delta^{2},s^{2},\nabla^{2}\delta\\}$.
Symbols display the measurements employing a randomly-selected cosmology from
our training data at $z\sim 0$. We note that these measurements corresponds to
cosmology-rescaled high-resolution simulations, which have been created with
the “Paired-&-Fixed” method. For comparison, solid lines show the predictions
of our analytic LPT calculation for each respective cross-spectrum.

We now describe the procedure we follow to construct an emulator for real-
space biased tracers. We first define our target cosmological parameter space
in §4.1, and then describe our core power spectrum data in §4.2. We continue
by discussing our neural network emulation strategy in §4.3 and finish by
quantifying its accuracy in §4.4.

###  4.1 Parameter Space

We will construct and validate our emulator over a hyper-volume in
cosmological parameter space defined by the ranges:

$\displaystyle\sigma_{8}$ $\displaystyle\in$ $\displaystyle[0.73,0.9]$
$\displaystyle\Omega_{\rm m}$ $\displaystyle\in$ $\displaystyle[0.23,0.4]$
$\displaystyle\Omega_{\rm b}$ $\displaystyle\in$ $\displaystyle[0.04,0.06]$
$\displaystyle n_{\rm s}$ $\displaystyle\in$ $\displaystyle[0.92,1.01]$ (7)
$\displaystyle h$ $\displaystyle\in$ $\displaystyle[0.6,0.8]$ $\displaystyle
M_{\rm\nu}\,[{\rm eV}]$ $\displaystyle\in$ $\displaystyle[0.0,0.4]$
$\displaystyle w_{0}$ $\displaystyle\in$ $\displaystyle[-1.15,-0.85]$
$\displaystyle w_{\rm a}$ $\displaystyle\in$ $\displaystyle[-0.3,0.3]$

where $\sigma_{8}$ is the cold mass linear mass variance in $8\,h^{-1}{\rm
Mpc}$ spheres; $\Omega_{\rm m}$ and $\Omega_{\rm b}$ are the density of cold
matter and baryons in units of the critical density of the Universe; $n_{\rm
s}$ is the primordial spectral index; $h$ is the dimensionless Hubble
parameter $h=H_{0}/(100\,{\rm km}\,{\rm s^{-1}}{\rm Mpc^{-1}})$; $M_{\rm\nu}$
is the mass of neutrinos in units of eV; and $w_{0}$ and $w_{\rm a}$ are
parameters describing the time-evolving dark energy equation of state via
$w(z)=w_{0}+(1-a)\,w_{\rm a}$. We further consider a flat geometry and neglect
the impact of radiation in the background evolution.

We note that this parameter volume coincides with that of the non-linear
matter power spectrum and baryonic effects emulators of Angulo et al. (2020)
and Aricò et al. (2020). We also recall that our suite of test simulations –
with which we validated the cosmology-rescaling approach, c.f. §4.4 – span the
same range of parameters.

It is worth highlighting that although this represents a restricted parameter
space, it is significantly larger than that typically covered by emulation of
$N$-body simulations. This is because the cosmology-rescaling allows us to
densely sample the volume, thus keeping emulation errors under control. Also,
outside this parameter range, it could be possible to gracefully extrapolate
using less precise methods. This could be acceptable as the region outside
these boundaries is expected to be disfavored by large-scale-structure data,
which would justify the use of less precise methods.

###  4.2 Data

We sample our cosmology hyper-space with $400$ points generated with a latin-
hypercube algorithm. For each point, we first select the simulation that would
yield most accurate rescaling results (see Contreras et al., 2020a, for a
details), then apply our rescaling approach at $10$ different output times
over $0<z<1.5$, and finally advect the corresponding Lagrangian fields.

For each rescaled output we measure the power spectra and compute the
quantities to be emulated as we will describe in the next two subsections.
This whole procedure takes approximately 2 CPU hours, with the largest
fraction spent in the power spectra calculation and in the creation of the
linear fields.

####  4.2.1 Power Spectrum Measurements

We compute 15 cross-spectra using a Fast Fourier Transform with $1080^{3}$
grid points and a Cloud-in-Cells mass assignment scheme. To reduce the impact
of aliasing, we employ two interlaced grids, following Sefusatti et al.
(2016). We consider $50$ logarithmically-spaced bins over the range of
wavenumbers $k\in[10^{-2},1]\,h\,{\rm Mpc}^{-1}$.

On large scales, the regularity of the Fourier grid imprints features in the
measured spectra. To correct for this, we compute the average of Fourier modes
over the LPT expectation evaluated at the same wavelength, and then multiply
by the LPT expectation at the bin centre, i.e.:

$\hat{P}_{ij}(k)=P_{{\rm LPT},ij}(k)\sum_{|\vec{k}_{n}|\in[k-\Delta k,k+\Delta
k]}\frac{\delta_{i}(\vec{k}_{n})\delta_{j}(\vec{k}_{n})}{P^{\rm
LPT}_{ij}(\vec{k}_{n})}$ (8)

we have found this procedure to be particularly important for the terms
involving the $1$ and $\delta$ fields on large scales.

We perform the mass assignment and FFT using shared-memory parallelization
over $24$ cores. The mass deposit steps typically take 8 seconds per field,
and the corresponding FFT and power spectrum calculation approximately 17
seconds.

####  4.2.2 Emulated quantities

Figure 5:  Ratio between the cross spectra in our validation set over the
predictions of our neural network emulation. As in previous plots, each panel
display results for a different combination of Lagrangian fields in Eulerian
coordinates. Lines of different colour show specific test cosmologies, whereas
shaded regions denote the are enclosing 95% of the distribution with the mean
indicated by the tick black line. Note the predictions are percent accurate
over most cases, and the largest fractional discrepancies occur when the
emulated quantities crosses zero.

In total, we have computed a set of 15 spectra for approximately $4,000$
combinations of cosmology and redshifts.

In Fig. 4 we show as an example one set of measurements for a randomly-
selected cosmology at $z\sim 0$. We also display our analytic LPT predictions
as solid lines.

We can see that our measurements are in good agreement with LPT on large
scales, as it was already shown in Fig. 3, especially for fields involving
combinations of $1$, $\delta$, and $\delta^{2}$. In some spectra, there seems
to be a small disagreement, which in §3 we attribute to either higher-order
LPT contributions and/or shell-crossing. However, these later spectra seem to
contribute subdominantly.

It is worth noting that there does not seem to be any significant noise in the
fields that dominate the large-scale clustering $\langle 11\rangle$, $\langle
1\delta\rangle$ and $\langle\delta\delta\rangle$. This owns to the
“Paired-&-Fixed”method, which, as shown by Angulo & Pontzen (2016) (see also
Pontzen et al., 2016; Chuang et al., 2018; Villaescusa-Navarro et al., 2018),
results in <$\sim 0.1\%$ accurate density power spectrum for $k<0.1h\,{\rm
Mpc}^{-1}$ at $z=0$ using simulations of comparable volume to those in the
BACCO suite.

To reduce the dynamic range of our measurements, which would yield to a more
accurate subsequent emulation, we consider the ratio of our measured cross-
spectra over the corresponding analytic LPT expectation. We further consider
the logarithm of the numerical and analytic results whenever that ratio is
always positive over our range of wavenumbers, i.e:

$S_{ij}(k,z)\equiv\begin{cases}P_{ij}/P^{\rm LPT}_{ij},&\text{if
}\mathrm{min}[P_{ij}]\leq 0\\\ \log(P_{ij}/P^{\rm LPT}_{ij}),&\text{otherwise
}\end{cases}$ (9)

It is well known that large-scale coherent motions dilute the amplitude of the
baryonic acoustic oscillations (BAO Eisenstein et al., 2007). This effect can
be captured in perturbation theory via the so-called infrared-resummation,
however, this is not included in our LPT calculations. In principle, this is
not a problem for our approach since we simply use LPT to reduce the dynamic
range of our data. In practice, on the other hand, this means a small
oscillatory feature is present in $S$ which we ought to emulate, which could
introduce small uncertainties.

To avoid this and improve the accuracy our our predictions, in the definition
of $\langle 11\rangle$, $\langle 1\delta\rangle$ and
$\langle\delta\delta\rangle$, we replace the corresponding LPT predictions by
a linear theory power spectrum including an infra-red resummation, i.e:

$P_{\rm linear}^{\rm smeared-BAO}\equiv P_{\rm linear}\,G(k)+P_{\rm
linear}^{\rm no-BAO}\,[1-G(k)]$ (10)

where $G(k)\equiv\exp[-k^{2}/k_{*}]$, with
$k_{*}^{-1}\equiv(6\pi^{2})^{-1}\int{\rm d}k\,P_{\rm linear}$, and $P_{\rm
linear}^{\rm no-BAO}$ is a version of the linear theory power spectrum without
any BAO signal. Operationally, this is obtained by performing a discrete sine
transform, smoothing the result, and returning to Fourier space by an inverse
transform (Baumann et al., 2018; Giblin et al., 2019). We have checked that
this in fact results into $S(k)$ devoid of any residual oscillations, however,
we emphasise that our predictions are still essentially provided by our
simulation results.

As discussed previously, simulations results have a significant amount of
noise for some cross-spectra. Most notably, $1s^{2}$, $1\delta^{2}$,
$\delta\delta^{2}$, and $\delta s^{2}$, in such cases emulation could be
inaccurate as we would be trying to predict specific noise features. To avoid
this, for the aforementioned spectra, we have replaced our measurements for
$k<0.1\,h\,{\rm Mpc}^{-1}$ with our analytic LPT spectra.

To further reduce the dimensionality and improve our emulation, we have
performed a principal component decomposition over our dataset. We have found
that keeping the first $6$ vectors with the largest eigenvalues is enough to
explain most of the variance in all of our dataset.

###  4.3 Neural Network

Following Angulo et al. (2020) and Aricò et al. (2020), we construct our
emulators using a feed-forward Neural Network. We use a relatively simple
architecture with two fully-connected hidden layers with 200 neurons each a
Rectified Linear Unit activation function. Note that a traditional alternative
for building emulators are Gaussian Processes, which have the advantage of
naturally providing with an estimate of the emulator uncertainty. These
methods are, however, computationally expensive in high dimensions and for
abundant training data. Thus, we choose to adopt neural networks instead to
allow future expansions of the training set by including additional points in
our parameter space, which could in principle reduce arbitrarily the emulation
error.

We build a separate neural network for each cross-spectrum using the Keras
front-end of the Tensor-flow library (Abadi et al., 2015). We select a
standard Adam optimization algorithm with a $10^{-3}$ learning rate with a
loss function given by the mean squared error.

We split our dataset, consisting of $4000$ sets of cross-spectra, into
disjoint groups for training and validation. The training set comprises 95% of
the data with which the training of each 15 of the emulators takes
approximately 30 minutes in a single Nvidia Quadro RTX 8000 GPU card; the
evaluation of each emulator takes approximately $0.05$ seconds on the same
hardware.

###  4.4 Validation

We have estimated the accuracy of our neural network by comparing the
prediction of our emulator with the measurements of the spectra in our test
suite comprised of $\sim 200$ different cosmologies and redshifts.

Our results are presented in Fig. 5, where we display the ratio of the values
of $S(k)$ predicted by our neural network and those directly measured in our
test suite. Each panel shows the result for one of our 15 combinations of
Lagrangian fields. We recall that in this testing set there are simulations
with very dissimilar cosmologies including massive neutrinos and dark energy
with time-dependent equation of state. Individual results are indicated by the
coloured lines whereas the mean and the 95% region of the distribution are
indicated by thick solid and shaded regions.

We can see that our neural networks (NN) perform remarkably well overall. For
the main fields for large-scale predictions – i.e. combinations of $\delta$,
$1$ – the NN shows a typical uncertainty below 1% at scales $k<0.1h\,{\rm
Mpc}^{-1}$. On smaller scales, the accuracy somewhat degrades but it is still
below 2% for all but 1 case (pink line).

For fields further including $\delta^{2}$ terms, the accuracy is similar, with
differences being less than 5% on small scales. Recall that in the case of
$1s^{2}$, $1\delta^{2}$, $\delta\delta^{2}$, and $\delta s^{2}$, we have
replaced our numerical results on $k<0.1h\,{\rm Mpc}^{-1}$ by our analytic LPT
expressions, which yields to an almost perfect emulation and to the apparent
discontinuity of the emulator precision at that transition scale.

Note that the largest fractional errors occur at $k\sim 0.2h\,{\rm Mpc}^{-1}$
for $\delta^{2}\nabla^{2}\delta$ and at $k\sim 0.7h\,{\rm Mpc}^{-1}$ for
$\delta^{2}s^{2}$, which is where these cross-spectra cross zero (in many
cosmologies), thus the ratio can easily become unbounded. However, since these
fields are close to zero, their absolute contribution to the total power
spectrum is negligible.

##  5 Application: Fitting the galaxy power spectrum

Figure 6: The galaxy auto power spectrum and matter-galaxy cross spectrum for
SFR-selected mock galaxies with $\bar{n}=10^{-3}\,h^{3}\,\mathrm{Mpc}^{-3}$ at
$z=1$. Symbols show the measurements after a galaxy-formation model in one of
our BACCO simulations, whereas blue lines denote the best fitting model based
in the emulator for biased tracers developed in our work.

As an illustration of our emulator of the real-space clustering of biased
tracers, in this section we will constrain the bias parameters that best fit
the power spectrum of a simulated galaxy catalogue. While in principle we
could directly jointly constrain bias parameters and cosmological parameters,
for the sake of this paper we will fix our cosmology to the fiducial one. We
do so to provide a quick test of the flexibility of the bias model that can be
build with our emulator. We defer the study of possible degeneracies between
cosmology and Lagrangian bias parameters to a future work.

We consider a mock galaxy catalogue mimicking a sample of Emission Line
galaxies at $z\sim 1$ with a number density of
$\bar{n}=10^{-3}\,h^{3}\,\mathrm{Mpc}^{-3}$, which can be regarded as
analogous to the sample to be observed by EUCLID or DESI. To build this mock
catalogue, we employ the extended subhalo abundance matching model of
Contreras et al. (2020c). These authors showed that this particular
implementation – which includes models for tidal stripping and disruption of
satellite galaxies, as well as for the star formation rate of galaxies –
accurately reproduces the real and reshift-space clustering of galaxies in the
state-of-the-art hydrodynamical simulation TNG-300 (Nelson et al., 2018;
Springel et al., 2018; Marinacci et al., 2018; Pillepich et al., 2018; Naiman
et al., 2018).

For this paper, we applied the aforementioned model to one of our main BACCO
simulations (nenya). We employ the galaxy-formation parameters that best
described SFR-selected TNG-300 galaxies, as provided by Contreras et al.
(2020c). Specifically we set
$\beta=2.490,\gamma=5.127,\log_{1}0M_{1}=12.770,\tau_{0}=4.931,$ and
$\tau_{\rm S}=-0.363$. We then selected the galaxies with the highest SFR
values and kept a sample with a number density of
$\bar{n}=10^{-3}\,h^{3}\,\mathrm{Mpc}^{-3}$. Please note that galaxies
obtained with this method contain some galaxy assembly bias signal, since they
inherit correlations with secondary halo properties (i.e. $v_{peak}$) and
local environment (see Contreras et al., 2020b, for a more comprehensive
treatment of GAB within SHAM methods).

Figure 7: Marginalised 1$\sigma$ and $2\sigma$ credibility regions in the
Lagrangian bias parameters $b_{1}$, $b_{2}$, $b_{s^{2}},$ and
$b_{\nabla^{2}\delta}$ for a sample of mock galaxies with
$\bar{n}=10^{-3}\,h^{3}\,\mathrm{Mpc}^{-3}$ selected according to their star
formation rate at $z=1$. These constraints were derived by simultaneously
fitting the galaxy power spectrum and galaxy-matter cross-spectrum from
$k=10^{-2}h\,{\rm Mpc}^{-1}$ until three different maximum wavenumbers,
$k_{\rm max}$, $0.1$, $0.3$, and $0.7h\,{\rm Mpc}^{-1}$, as indicated by the
legend.

We will consider as our data the average clustering measured in the two phase-
inverted Nenya simulations to reduce the impact of cosmic variance on large
scales. Specifically, we will employ the galaxy auto power spectrum, $P_{\rm
gg}$, and the cross-correlation, $P_{\rm gm}$, with the underlying nonlinear
matter field. As for the previous power spectrum calculations, we perform this
operation using $1080^{3}$ FFTs with a CiC assignment and interlacing. In Fig.
6 we display the measurements as red symbols. Note that $P_{\rm gg}$ and
$P_{\rm gm}$ involve different combinations of the Lagrangian fields, which
makes it important to prove that the model is able to jointly fit them.

We will now explore the values of the free parameters of our real-space model,
which is described by combinations of the cross spectra of Lagrangian fields
and one scale independent shot-noise contribution:

$\displaystyle P_{\rm gg}=$
$\displaystyle\sum_{i,j\in\\{1,\delta,\delta^{2},\nabla^{2}\delta\\}}(2-\delta_{ij})b_{i}b_{j}\,P_{ij}+\frac{A_{\rm
sn}}{\bar{n}}$ (11) $\displaystyle P_{\rm gm}=$
$\displaystyle\sum_{i=1,j\in\\{1,\delta,\delta^{2},\nabla^{2}\delta\\}}b_{j}\,P_{ij}$
(12)

where $b_{i=1}\equiv 1$, and $P_{ij}$ are the cross spectra predicted by our
neural network for the cosmological parameters of the target galaxy catalogue
(i.e. those of Nenya ). Note that in principle, there could be other
contributions to the noise (e.g. Eggemeier et al., 2020), which could be
included in our model but that we have neglected here. However, our model is
flexible enough for it to be extended with additional terms. We defer an
investigation of a minimal model for describing realistic galaxies to future
work. Our model is thus specified by 5 parameters:

$\boldsymbol{\vartheta}=\\{b_{1},b_{2},b_{s^{2}},b_{\nabla^{2}\delta},A_{\rm
sn}\\}$ (13)

where $b_{i}$ are the perturbative bias parameters of our sample. We assume
flat priors for these parameters over the range $b_{1}\in[-3,3]$,
$b_{2}\in[-5,5]$, $b_{s^{2}}\in[-10,20]$, $b_{\nabla^{2}\delta}\in[-10,20]$,
$A_{\rm sn}\in[0,2]$.

We describe the probability of observing a set of power spectra
$\tilde{P}(k)=\\{P_{\rm gg}(k),P_{\rm gm}(k)\\}$ as a multivariate normal
distribution:

$\ln
p[\tilde{P}(k),\boldsymbol{\vartheta}]=-\dfrac{1}{2}\sum_{i}\left[\dfrac{\tilde{P}(k_{i})-\tilde{P}_{\rm
model}(k_{i})}{\sigma_{i}}\right]^{2}.$ (14)

Where we have assumed a diagonal covariance matrix computed in the Gaussian
limit with a volume corresponding to $V=50h^{-3}\mathrm{Gpc}^{3}$, roughly
comparable to that to be surveyed by EUCLID. To account for the uncertainty
associated to our emulation process, we have added in quadrature a constant
fractional error of $2\%$.

We sample the posterior probability of $\boldsymbol{\vartheta}$,

$p[\boldsymbol{\vartheta},\tilde{P}(k)]=\dfrac{p[\tilde{P}(k),\boldsymbol{\vartheta}]p(\boldsymbol{\vartheta})}{p[\tilde{P}(k)]},$
(15)

using affine-invariant MonteCarlo Markov Chain code emcee (Foreman-Mackey et
al., 2013). At each step of the chain we evaluate our neural network emulator
to obtain $S_{ij}$ and evaluate the corresponding LPT expressions to recover
$P_{ij}$ and combine them as defined in Eqs. 11. We consider $20$ MCMC walkers
with 10,000 steps each and a burn-in phase of 500 steps. We checked that
doubling the number of steps does not change significantly our results and
therefore we consider our runs converged.

In Fig. 7 we present the two-dimensional marginalised constraints on the free
parameters of our model for the galaxy power spectra. Dark and light regions
denote the 1-$\sigma$ and 2-$\sigma$ two-dimensional confidence levels. We
display the results for when limiting our analysis to different maximum
wavenumbers, namely $k_{\rm max}=0.1,0.3,0.7h^{-1}{\rm Mpc}$. We note that, as
we include smaller scales, even though the width of the contours changes, the
best fitting parameters remain compatible.

Finally, in Fig. 6 we compare the measured galaxy statistics to those
predicted by our model at the values that maximise the posterior probability.
We can see that our model is an excellent description of our simulated galaxy
catalogues, over all the scales we consider and for both $P_{\rm gg}$ and
$P_{\rm gm}$. In particular, our model evaluated for the best-fitting set of
parameters is always within the error-bars of the data, and agrees with the
measured power spectra at the $1-2\%$ level. This is particularly remarkable
especially keeping in mind that, typically, perturbation-theory approaches are
limited to scales $k\sim 0.2h\,{\rm Mpc}^{-1}$, and that we compare against a
galaxy catalogue built directly in Eulerian space and following physically-
motivated recipes.

##  6 Conclusions

In this paper we have built and presented a set of 15 emulators for predicting
the cross power spectra necessary to describe the galaxy-galaxy and galaxy-
matter clustering in the context of a second order Lagrangian bias expansion.
Each emulator predicts one such term on a range of wavenumbers $0.01<k/h\,{\rm
Mpc}^{-1}<1$ and at redshifts $0<z<1.5$, in cosmologies spanning $\sim
10\sigma$ around the best-fitting parameters from the Planck collaboration,
including extensions of the standard $\Lambda$CDM paradigm, such as massive
neutrino and dynamical dark energy.

In order to build our emulator from a densely sampled parameter space, we have
used a cosmology-rescaling algorithm. For this reason, we have first used a
set of 35 intermediate volume simulations, with cosmologies spanning the same
parameter space as our emulator, to show that the Lagrangian fields measured
from simulations assuming different cosmologies can be accurately reproduced
rescaling the small set of strategically chosen BACCo simulations presented in
Contreras et al. (2020a). From this exercise (see Fig. 2) we conclude that the
rescaling technique can accurately reproduce the Lagrangian fields needed to
model Lagrangian bias expansion.

We have devoted particular care to obtaining theoretical predictions for the
power spectra of these Lagrangian fields. We have considered LPT terms up to
second order in density and advected them to Eulerian space. In Fig.3 we have
shown the accuracy of the expressions obtained this way, comparing them to
numerical simulations of these LPT fields. We have used out LPT predictions
for a two-fold purpose: one the one hand, we can combine the predictions to
the actual data to mitigate the large-scale noise affecting some of the terms;
on the other hand, we want to use them to reduce the dynamical range of the
quantities that we emulate.

We have built our 15 emulators rescaling the large BACCO simulations to
approximately 4000 different combinations of cosmology and redshifts (see Fig.
4) and we have then proceeded to validate them against a set of 200
combinations of cosmologies and redshifts (Fig. 5). We find that most of the
terms that principally contribute to the galaxy power spectrum (namely, the
terms involving combinations of the homogeneous and linear fields, $1,\delta$)
are emulated with percent accuracy. The terms involving the squared density
field $\delta^{2}$, the tidal field $s^{2}$ and the laplacian of the density
field $\nabla^{2}\delta$ are generally reproduced with an accuracy of a few
percent. However, two of these terms, the combination of $\delta^{2}s^{2}$ and
of $\delta^{2}\nabla^{2}\delta$, show particularly bad performances at the
scale where they cross zero.

Finally, we have used our emulator to constrain the bias parameters of a mock
galaxy sample of known cosmology, using the galaxy-galaxy and galaxy-matter
power spectra. We constructed our galaxy sample at redshift $z=1$, selecting
galaxies according to their Star Formation Rate until reaching a number
density of $0.001h^{3}\,\mathrm{Mpc}^{-3}$ and assigning errors that reflect
the cosmic variance expected for large upcoming surveys such as Euclid. This
proof of concept showed that this technique can be used to model the
clustering of realistic galaxy samples, obtaining constraints of the bias
parameters that are self consistent along a large range of scales (Fig. 6). In
particular we were able to fit our mock galaxy sample including scales down to
$k_{\rm max}=0.7h\,{\rm Mpc}^{-1}$.

We anticipate this approach will be valuable for studying Lagrangian biases
from realistic galaxy mock samples. Its applications could be extended also to
any problem in which the knowledge of the fully nonlinear clustering of each
of these lagrangian fields is required, as for example in the modelling of
galaxy-galaxy lensing. However non trivial, we also expect that an extension
of this emulator to redshift space would provide us with an unprecedented tool
to constrain cosmology from galaxy clustering.

In the latest stages of the preparation of this work, a similar approach to
ours has been made public in Kokron et al. (2021). The two works share many
similarities, although differ in many aspects of the implementation and
analysis. We defer to future works a comparison of the two approaches.

## Acknowledgments

The authors acknowledge the support of the ERC-StG number 716151 (BACCO). MPI
acknowledges the support of the “Juan de la Cierva Formación” fellowship
(FJC2019-040814-I). SC acknowledges the support of the “Juan de la Cierva
Formación” fellowship (FJCI-2017-33816). The authors acknowledge the computer
resources at MareNostrum and the technical support provided by Barcelona
Supercomputing Center (RES-AECT-2019-2-0012, RES-AECT-2020-3-0014).

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author. The Neural Network emulator will be made public at
http://www.dipc.org/bacco upon the publication of this article.

## References

  * Abadi et al. (2015) Abadi M., et al., 2015, TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems, http://tensorflow.org/
  * Abidi & Baldauf (2018) Abidi M. M., Baldauf T., 2018, Journal of Cosmology and Astroparticle Physics, 2018, 029
  * Amendola et al. (2018) Amendola L., et al., 2018, Living Reviews in Relativity, 21, 2
  * Angulo & Pontzen (2016) Angulo R. E., Pontzen A., 2016, MNRAS, 462, L1
  * Angulo & White (2010) Angulo R. E., White S. D. M., 2010, MNRAS, 405, 143
  * Angulo et al. (2020) Angulo R. E., Zennaro M., Contreras S., Aricò G., Pellejero-Ibañez M., Stücker J., 2020, arXiv e-prints, p. arXiv:2004.06245
  * Aricò et al. (2020) Aricò G., Angulo R. E., Contreras S., Ondaro-Mallea L., Pellejero-Ibañez M., Zennaro M., 2020, arXiv e-prints, p. arXiv:2011.15018
  * Baldauf et al. (2016) Baldauf T., Schaan E., Zaldarriaga M., 2016, Journal of Cosmology and Astroparticle Physics, 2016, 017
  * Bartelmann & Schneider (2001) Bartelmann M., Schneider P., 2001, Physics Reports, 340, 291
  * Baumann et al. (2012) Baumann D., Nicolis A., Senatore L., Zaldarriaga M., 2012, Journal of Cosmology and Astroparticle Physics, 2012, 051
  * Baumann et al. (2018) Baumann D., Green D., Wallisch B., 2018, J. Cosmology Astropart. Phys., 2018, 029
  * Bernardeau et al. (2002) Bernardeau F., Colombi S., Gaztañaga E., Scoccimarro R., 2002, Phys. Rep., 367, 1
  * Blas et al. (2014) Blas D., Garny M., Konstandin T., 2014, Journal of Cosmology and Astroparticle Physics, 2014, 010
  * Bonoli et al. (2020) Bonoli S., et al., 2020, arXiv e-prints, p. arXiv:2007.01910
  * Chen et al. (2020) Chen S.-F., Vlah Z., White M., 2020, J. Cosmology Astropart. Phys., 2020, 062
  * Chuang et al. (2017) Chuang C.-H., et al., 2017, Monthly Notices of the Royal Astronomical Society, 471, 2370
  * Chuang et al. (2018) Chuang C.-H., et al., 2018, arXiv e-prints, p. arXiv:1811.02111
  * Chudaykin et al. (2020) Chudaykin A., Ivanov M. M., Simonović M., 2020, arXiv e-prints, p. arXiv:2009.10724
  * Colas et al. (2020) Colas T., d'Amico G., Senatore L., Zhang P., Beutler F., 2020, Journal of Cosmology and Astroparticle Physics, 2020, 001
  * Contreras et al. (2020a) Contreras S., Zennaro R. E. A. M., Aricó G., Pellejero-Ibañez M., 2020a, arXiv e-prints, p. arXiv:2001.03176
  * Contreras et al. (2020b) Contreras S., Angulo R., Zennaro M., 2020b, arXiv e-prints, p. arXiv:2005.03672
  * Contreras et al. (2020c) Contreras S., Angulo R., Zennaro M., 2020c, arXiv e-prints, p. arXiv:2012.06596
  * DESI Collaboration et al. (2016) DESI Collaboration et al., 2016, arXiv e-prints, p. arXiv:1611.00036
  * DeRose et al. (2019) DeRose J., et al., 2019, ApJ, 875, 69
  * Desjacques et al. (2018a) Desjacques V., Jeong D., Schmidt F., 2018a, Phys. Rep., 733, 1
  * Desjacques et al. (2018b) Desjacques V., Jeong D., Schmidt F., 2018b, J. Cosmology Astropart. Phys., 2018, 035
  * Dodelson & Schneider (2013) Dodelson S., Schneider M. D., 2013, Phys. Rev. D, 88, 063537
  * Dore et al. (2019) Dore O., et al., 2019, Bulletin of the AAS, 51
  * Eggemeier et al. (2020) Eggemeier A., Scoccimarro R., Crocce M., Pezzotta A., Sánchez A. G., 2020, Phys. Rev. D, 102, 103530
  * Eisenstein & Hu (1998) Eisenstein D. J., Hu W., 1998, ApJ, 496, 605
  * Eisenstein et al. (2007) Eisenstein D. J., Seo H.-J., Sirko E., Spergel D. N., 2007, ApJ, 664, 675
  * Euclid Collaboration et al. (2019) Euclid Collaboration et al., 2019, MNRAS, 484, 5509
  * Euclid Collaboration et al. (2020) Euclid Collaboration et al., 2020, arXiv e-prints, p. arXiv:2010.11288
  * Fang et al. (2017) Fang X., Blazek J. A., McEwen J. E., Hirata C. M., 2017, J. Cosmology Astropart. Phys., 2017, 030
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Fujita & Vlah (2020) Fujita T., Vlah Z., 2020, Journal of Cosmology and Astroparticle Physics, 2020, 059
  * Fujita et al. (2020) Fujita T., Mauerhofer V., Senatore L., Vlah Z., Angulo R., 2020, J. Cosmology Astropart. Phys., 2020, 009
  * Garrison et al. (2018) Garrison L. H., Eisenstein D. J., Ferrer D., Tinker J. L., Pinto P. A., Weinberg D. H., 2018, The Astrophysical Journal Supplement Series, 236, 43
  * Giblin et al. (2019) Giblin B., Cataneo M., Moews B., Heymans C., 2019, MNRAS, 490, 4826
  * Guzzo et al. (2008) Guzzo L., et al., 2008, Nature, 451, 541
  * Heitmann et al. (2014) Heitmann K., Lawrence E., Kwan J., Habib S., Higdon D., 2014, ApJ, 780, 111
  * Ivanov et al. (2020) Ivanov M. M., Simonović M., Zaldarriaga M., 2020, Journal of Cosmology and Astroparticle Physics, 2020, 042
  * Ivezić et al. (2019) Ivezić Ž., et al., 2019, The Astrophysical Journal, 873, 111
  * Izard et al. (2016) Izard A., Crocce M., Fosalba P., 2016, MNRAS, 459, 2327
  * Kaiser (1987) Kaiser N., 1987, MNRAS, 227, 1
  * Kokron et al. (2021) Kokron N., DeRose J., Chen S.-F., White M., Wechsler R. H., 2021, arXiv e-prints, p. arXiv:2101.11014
  * Kuhlen et al. (2012) Kuhlen M., Vogelsberger M., Angulo R., 2012, Physics of the Dark Universe, 1, 50
  * Lazeyras & Schmidt (2019) Lazeyras T., Schmidt F., 2019, Journal of Cosmology and Astroparticle Physics, 2019, 041
  * Lesgourgues (2011) Lesgourgues J., 2011, ArXiv e-prints 1104.2932,
  * Levi et al. (2013) Levi M., et al., 2013, preprint,  (arXiv:1308.0847)
  * Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473
  * Liu et al. (2018) Liu J., Bird S., Zorrilla Matilla J. M., Hill J. C., Haiman Z., Madhavacheril M. S., Petri A., Spergel D. N., 2018, Journal of Cosmology and Astro-Particle Physics, 2018, 049
  * Mandelbaum (2018) Mandelbaum R., 2018, Annual Review of Astronomy and Astrophysics, 56, 393
  * Marinacci et al. (2018) Marinacci F., et al., 2018, MNRAS, 480, 5113
  * Matsubara (2008) Matsubara T., 2008, Phys. Rev. D, 78, 083519
  * McDonald (2009) McDonald P., 2009, Journal of Cosmology and Astroparticle Physics, 2009, 026
  * McEwen et al. (2016) McEwen J. E., Fang X., Hirata C. M., Blazek J. A., 2016, J. Cosmology Astropart. Phys., 2016, 015
  * McQuinn & White (2016) McQuinn M., White M., 2016, Journal of Cosmology and Astroparticle Physics, 2016, 043
  * Michaux et al. (2021) Michaux M., Hahn O., Rampf C., Angulo R. E., 2021, MNRAS, 500, 663
  * Modi et al. (2020) Modi C., Chen S.-F., White M., 2020, MNRAS, 492, 5754
  * Monaco et al. (2002) Monaco P., Theuns T., Taffoni G., 2002, MNRAS, 331, 587
  * Naiman et al. (2018) Naiman J. P., et al., 2018, MNRAS, 477, 1206
  * Nelson et al. (2018) Nelson D., et al., 2018, MNRAS, 475, 624
  * Nishimichi et al. (2019) Nishimichi T., et al., 2019, ApJ, 884, 29
  * Nishimichi et al. (2020) Nishimichi T., D’Amico G., Ivanov M. M., Senatore L., Simonović M., Takada M., Zaldarriaga M., Zhang P., 2020, Phys. Rev. D, 102, 123541
  * Paz & Sánchez (2015) Paz D. J., Sánchez A. G., 2015, MNRAS, 454, 4326
  * Pellejero-Ibanez et al. (2017) Pellejero-Ibanez M., et al., 2017, Monthly Notices of the Royal Astronomical Society, 468, 4116
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 475, 648
  * Pontzen et al. (2016) Pontzen A., Slosar A., Roth N., Peiris H. V., 2016, Phys. Rev. D, 93, 103519
  * Schneider et al. (2016) Schneider A., et al., 2016, Journal of Cosmology and Astro-Particle Physics, 2016, 047
  * Sefusatti et al. (2016) Sefusatti E., Crocce M., Scoccimarro R., Couchman H. M. P., 2016, MNRAS, 460, 3624
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel et al. (2018) Springel V., et al., 2018, MNRAS, 475, 676
  * Springel et al. (2020) Springel V., Pakmor R., Zier O., Reinecke M., 2020, arXiv e-prints, p. arXiv:2010.03567
  * Takada et al. (2014) Takada M., et al., 2014, Publications of the Astronomical Society of Japan, 66
  * Tassev et al. (2013) Tassev S., Zaldarriaga M., Eisenstein D. J., 2013, JCAP, 6, 036
  * Taylor et al. (2013) Taylor A., Joachimi B., Kitching T., 2013, MNRAS, 432, 1928
  * Villaescusa-Navarro et al. (2018) Villaescusa-Navarro F., et al., 2018, ApJ, 867, 137
  * Vlah et al. (2016) Vlah Z., Castorina E., White M., 2016, Journal of Cosmology and Astroparticle Physics, 2016, 007
  * Wibking et al. (2019) Wibking B. D., et al., 2019, MNRAS, 484, 989
  * 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
  * Zennaro et al. (2019) Zennaro M., Angulo R. E., Aricò G., Contreras S., Pellejero-Ibáñez M., 2019, MNRAS, 489, 5938
  * d'Amico et al. (2020) d'Amico G., Gleyzes J., Kokron N., Markovic K., Senatore L., Zhang P., Beutler F., Gil-Marín H., 2020, Journal of Cosmology and Astroparticle Physics, 2020, 005

##  Appendix A LPT power spectra

In Eq. (6) we presented a general form expressing the cross power spectra
$P_{ij}$ of the different fields, with
$i,j=\\{1,\delta,\delta^{2},s^{2},\nabla^{2}\delta\\}$. Here, we report the
prefactors of the terms of order $(11)$ and $(22)$, and the different kernels
entering the $(22)$ terms.

The $(11)$ terms are multiplied by these prefactors (note the $k$-dependence
arising for terms including the Laplacian of the density field),

$A_{ij}(k)=\begin{bmatrix}1&1&0&0&-k^{2}\\\ 1&1&0&0&-k^{2}\\\ 0&0&0&0&0\\\
0&0&0&0&0\\\ -k^{2}&-k^{2}&0&0&k^{4}\\\ \end{bmatrix}.$ (16)

The $(22)$ terms are preceded by

$B_{ij}(k)=\begin{bmatrix}1&1&1&1&0\\\ 1&1/2&2&2&0\\\ 1&2&1&1&1\\\
1&2&1&1&1\\\ 0&0&1&1&1\\\ \end{bmatrix}.$ (17)

Finally, we report all the kernels for the $(22)$ terms. Please note that we
use

$F_{\rm
ZA}(\boldsymbol{k}_{1},\boldsymbol{k}_{2})=1+\dfrac{\boldsymbol{k}_{1}\cdot\boldsymbol{k}_{2}}{k_{1}k_{2}}\left(\dfrac{k_{1}}{k_{2}}+\dfrac{k_{2}}{k_{1}}\right)+\left(\dfrac{\boldsymbol{k}_{1}\cdot\boldsymbol{k}_{2}}{k_{1}k_{2}}\right)^{2},$
(18)

and

$S_{2}(\boldsymbol{k}_{1},\boldsymbol{k}_{2})=\left(\dfrac{\boldsymbol{k}_{1}\cdot\boldsymbol{k}_{2}}{k_{1}k_{2}}\right)^{2}-\dfrac{1}{3}.$
(19)

$C_{11}(\boldsymbol{p},\boldsymbol{k})=F_{\rm
ZA}^{2}(\boldsymbol{p},\boldsymbol{k}-\boldsymbol{p})$ (20)
$C_{1\delta}(\boldsymbol{p},\boldsymbol{k})=F_{\rm
ZA}(\boldsymbol{p},\boldsymbol{k}-\boldsymbol{p})\left[\dfrac{\boldsymbol{k}\cdot(\boldsymbol{k}-\boldsymbol{p})}{|\boldsymbol{k}-\boldsymbol{p}|^{2}}+\dfrac{\boldsymbol{k}\cdot\boldsymbol{p}}{p^{2}}\right]$
(21) $C_{1\delta^{2}}(\boldsymbol{p},\boldsymbol{k})=F_{\rm
ZA}(\boldsymbol{p},\boldsymbol{k}-\boldsymbol{p})$ (22)
$C_{1s^{2}}(\boldsymbol{p},\boldsymbol{k})=F_{\rm
ZA}(\boldsymbol{p},\boldsymbol{k}-\boldsymbol{p})S_{2}(\boldsymbol{p},\boldsymbol{k}-\boldsymbol{p})$
(23) $C_{1\nabla^{2}\delta}(\boldsymbol{p},\boldsymbol{k})=0$ (24)
$C_{\delta\delta}(\boldsymbol{p},\boldsymbol{k})=\dfrac{\boldsymbol{k}\cdot(\boldsymbol{k}-\boldsymbol{p})}{|\boldsymbol{k}-\boldsymbol{p}|^{2}}\left[\dfrac{\boldsymbol{k}\cdot(\boldsymbol{k}-\boldsymbol{p})}{|\boldsymbol{k}-\boldsymbol{p}|^{2}}+\dfrac{\boldsymbol{k}\cdot\boldsymbol{p}}{p^{2}}\right]$
(25)
$C_{\delta\delta^{2}}(\boldsymbol{p},\boldsymbol{k})=\dfrac{\boldsymbol{k}\cdot(\boldsymbol{k}-\boldsymbol{p})}{|\boldsymbol{k}-\boldsymbol{p}|^{2}}$
(26) $C_{\delta
s^{2}}(\boldsymbol{p},\boldsymbol{k})=\dfrac{\boldsymbol{k}\cdot(\boldsymbol{k}-\boldsymbol{p})}{|\boldsymbol{k}-\boldsymbol{p}|^{2}}S_{2}(\boldsymbol{p},\boldsymbol{k}-\boldsymbol{p})$
(27) $C_{\delta\nabla^{2}\delta}(\boldsymbol{p},\boldsymbol{k})=0$ (28)
$C_{\delta^{2}\delta^{2}}(\boldsymbol{p},\boldsymbol{k})=1$ (29)
$C_{\delta^{2}s^{2}}(\boldsymbol{p},\boldsymbol{k})=S_{2}(\boldsymbol{p},\boldsymbol{k}-\boldsymbol{p})$
(30)
$C_{\delta^{2}\nabla^{2}\delta}(\boldsymbol{p},\boldsymbol{k})=\left[p^{2}\dfrac{\boldsymbol{k}\cdot(\boldsymbol{k}-\boldsymbol{p})}{|\boldsymbol{k}-\boldsymbol{p}|^{2}}+|\boldsymbol{k}-\boldsymbol{p}|^{2}\dfrac{\boldsymbol{k}\cdot\boldsymbol{p}}{p^{2}}\right]$
(31)
$C_{s^{2}\nabla^{2}\delta}(\boldsymbol{p},\boldsymbol{k})=\left[p^{2}\dfrac{\boldsymbol{k}\cdot(\boldsymbol{k}-\boldsymbol{p})}{|\boldsymbol{k}-\boldsymbol{p}|^{2}}+|\boldsymbol{k}-\boldsymbol{p}|^{2}\dfrac{\boldsymbol{k}\cdot\boldsymbol{p}}{p^{2}}\right]S_{2}(\boldsymbol{p},\boldsymbol{k}-\boldsymbol{p})$
(32)
$C_{\nabla^{2}\delta\nabla^{2}\delta}(\boldsymbol{p},\boldsymbol{k})=\left[\dfrac{p^{4}}{|\boldsymbol{k}-\boldsymbol{p}|^{4}}(k^{2}-\boldsymbol{k}\cdot\boldsymbol{p})^{2}+\boldsymbol{k}\cdot\boldsymbol{p}(k^{2}-\boldsymbol{k}\cdot\boldsymbol{p})\right]$
(33)

##  Appendix B Comparison against public PT codes

In this appendix we compare our implementation of LPT integrals and the
resulting prediction for Lagrangian fields in Eulerian space against those
computed using two publicly available codes.

Specifically, we compare against FastPT, which implements standard
perturbation theory at 1 loop, and against velocileptors, which employs 1-loop
Lagrangian perturbation theory with effective-field theory counterterms. Note
that not all power spectrum combinations are available in those codes, thus we
restrict the comparison to a subset of integrals.

Figure 8: Various cross-spectra, $P_{ij}$ at $z=0$, as indicated by the
legend. Shaded areas represent the standard deviation of 1000 fast simulations
evolved with 2LPT. Solid lines show our results, whereas dashed lines and
crosses denote the values obtained using velocileptors and FastPT,
respectively.

We first compare our LPT solutions against velocileptors. We see that, on
large scales, there is a good agreement for $11$, $1\delta$ $\delta\delta$,
$\delta\delta^{2}$, $\delta s^{2}$, $\delta^{2}\delta^{2}$, $\delta^{2}s^{2}$,
and $s^{2}s^{2}$. We find some differences in the large scale behavior of the
terms $1\delta^{2}$ and $1s^{2}$ that we plan to ivestigate in the future.
When compared against FastPT, we also see a good agreement with our
predictions on large scales for $\delta\delta$, $\delta^{2}\delta^{2}$ and
$s^{2}s^{2}$, but systematic differences for $\delta\delta^{2}$, $\delta
s^{2}$.

Our ensemble of LPT simulations shows very good agreement with our predictions
(as well as with those measured in full $N$-body simulations (c.f. §3). Note
that, for the term $1\delta^{2}$, the LPT simulations appear to have less
power than the predictions. However, for this term, the 2LPT simulations also
exhibits less power than the $N$-body solution, that in turn is in better
agreement with our predictions, as shown in Fig. 3.


# Cortical thickness and functional networks modules by cortical lobes

Vesna Vuksanović

Aberdeen Biomedical Imaging Centre
University of Aberdeen

<EMAIL_ADDRESS>
##  1 Abstract

This study aims to investigate topological organization of cortical thickness
and functional networks by cortical lobes. First, I demonstrated modular
organization of these networks by the cortical surface frontal, temporal,
parietal and occipital divisions. Secondly, I mapped the overlapping edges of
cortical thickness and functional networks for positive and negative
correlations. Finally, I showed that overlapping positive edges map onto
within-lobe cortical interactions and negative onto between-lobes
interactions.

##  2 Introduction

Classical, localized characterization of human brain surface morphology has
moved in recent years towards maps of inter-regional co-variations of gray
matter volume or thickness measurements in structural magnetic resonance
imaging data (sMRI) [1]. To make such cortical maps or brain networks (graphs)
the edges (links) between the regions (nodes) are defined by the strength of
correlation between regional volume or cortical thickness. The physiological
mechanisms underlying the thickness correlation among cortical areas remain
unclear. Neuroimaging findings propose strong correlation of the cortical
thickness measurements between regions that are directly connected via axonal
connections [2, 3]. However, only a fraction of variations in regional
thickness correlations across the cortex has been explained by direct axonal
links [3]. In this context it is interesting to ask whether the regional
thickness (morphological) co-variations also reflect functional associations
(or parts of them). Recent studies suggest that functional network properties,
such as the strengths or the number of correlations, can impact structural
properties. The synchronization of neuronal activity in response to specific
functional demands might induce synchronized plastic changes among related
regions during brain development [4] or degenerative diseases [5]. More
pertinent to this study, such interplay between regional structural and
functional properties may contribute to strong thickness correlations observed
within the visual areas of the occipital lobe [6]. Here it is also worth
noting that, while it has been observed that functional modules mirror the
local brain anatomy a number of studies show deviations between these two
suggesting many-to-one function-structure mapping [7].

The natural divisions of the cortical surface – frontal, temporal, parietal
and occipital – have been defined by regional morphological characteristics,
which are also known to support different functions [8, 9]. In this context,
questions arise about the relationship between cortical surface inter-regional
correlations, i.e., how do variations in functional and thickness correlations
map onto each other? To what extent these cortical structural and functional
correlations share similar topological organization? If topological measures
of functional network organization are structurally relevant we might expect
them to be impacted at the lobe level of the cortical surface organization.

To address these questions I examined thickness and functional correlational
networks based on two different neuroimaging modalities – structural and
functional MRI. The networks were represented by $68\times 68$ correlation
matrices, where functional and thickness network maps inter-regional
correlations between brain activity and thickness measurements respectively.
The main focus of the study is to shed light on intrinsic topological
properties of the frontal, temporal, parietal and occipital divisions, which
may impact structural and functional interactions across the cortical surface.
To this end I studied how modular organization of the cortical surface
divisions affect overlap of the cortical thickness and functional correlations
across cortical regions.

##  3 Methods

The MRI data sets considered in this short communication are from a public
database ($http://fcon_{1}000.projects.nitrc.org/$) provided by the Max Planck
Institute (MPI) Leipzig. A study group included 37 (16M/21F) participants
between 20 and 42 years. MRI data was acquired on 3 Tesla Magnetom Tim Trio
scanner (Siemens, Erlangen, Germany) using a 32-channel head coil. T1-weighted
images were acquired using a MPRAGE sequence (TR = 1.3 s; TE=3.46 ms; flip
angle=10$\deg$; FOV=$256\times 240mm^{2}$; 176 sagittal slices; voxel size =
$1\times 1\times 1.5mm$). Functional MRI/EPI data were acquired on a 3T MRI
scanner (Siemens Tim Trio) using TR=2.3 sec, TE=30ms, $3\times 3$ in-plane
resolution, 3 mm slice thickness, 1 mm gap between slices. Each scanning
session was a task-absent (”resting state”) scan lasting 7.6 minutes during
which subjects were asked to fixate a fixation cross.

Processing the structural MRI data was done using FreeSurfer v5.3.0
(https://surfer.nmr.mgh.harvard.edu/) pipeline according to the procedure
described in more details in [6]. Cortical thickness was measured on $N$ = 68
cortical surface regions segmented according to the Desikan-Keliany Atlas
(DKA) [10]. FSL v5.0.11 (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/) was used to
process fMRI data and extract BOLD fMRI time-series on the same 68 cortical
regions of the DKA parcellations. The fMRI data were processed using the FSL
toolbox FEAT fMRI analysis, similarly to our previous work [11] and global
signal regression was applied in order to reduce motion artifacts [12].

Figure 1: Workflow diagram of the cortical networks’ construction on processed
structural and functional magnetic resonance images. (Upper panel) Structural
magnetic resonance images (sMRI) are processed using FreeSurfer v5.3.0
pipeline: cortical surface was reconstructed and regional cortical thickness
measures were averaged over N = 68 structures of the Desikan-Keliany Atlas
(DKA). Structural correlation network was constructed on partial correlations
between regional thickness across all participants (S = 37) while controlling
for age, sex and mean CT. (Lower panel) Functional MRI (fMRI) were parcellated
and time series were averaged over the same 68 DKA regions. Individual
functional correlations (FC) matrix was constructed on pair-wise correlations
between the fMRI time series, Fisher z-transformed, averaged over S subjects
and transformed back to linear Pearson coefficient. Group-wise CT and FC
matrices were reordered according to node affiliation with cortical lobes
(frontal, temporal, parietal and occipital) and their coupling topology was
assessed to capture positive and negative correlations maps.

###  3.1 Brain Graph Construction

Brain graphs considered in this study were constructed on correlations between
either regional cortical thickness or BOLD activity measured by structural or
functional MRI. In these correlational networks, a brain region represent a
network node and a pair-wise correlation between the regional measurements
represents network’s edge/link or connection between the nodes. See Fig. 1
that summarizes analysis pipeline and approach.

####  3.1.1 Functional Correlations

To obtain the FC matrices, time series of the 68 DKA regions were calculated
by averaging the respective BOLD signals over all regional voxels. The
correlation matrices were obtained for each subject in the study. Since voxels
are in MNI space and a given voxel has approximately the same anatomical
position in all subjects, the individual correlation matrices can be averaged
across subjects. In detail, each matrix is first Fisher z-transformed and
averaged across subjects and then transformed back into correlation
coefficients again as described in previous studies [11, 13]. The resulting
correlation matrix $f_{ij},(i,j=1,…,N)$, was used to create a group-wise
functional correlation matrix between the regions of interest. See Fig. 1
(lower panel) for analysis pipeline and approach.

####  3.1.2 Cortical Thickness Correlations

A group-wise structural network was constructed on inter-regional cortical
thickness correlations across all study participants. The correlations were
calculated between each pair of 68 regional thickness measurements while
controlling for age, gender and mean CT, similarly to previously described
methodology [6, 3]. See Fig. 1 (upper panel) for analysis pipeline and
approach.

###  3.2 Brain Graph Analysis

####  3.2.1 Network Density

Network topological properties depend on the network thresholding. Threshold
affects network density ($\kappa$), i.e., number of links relative to the
total number of all possible links in the network [14]. The threshold
considered here was set to yield a fully connected network, i.e., one network
component. In addition, the threshold was applied to ensure that both networks
CT and FC have the same density, i.e., to ensure that the networks are
analyzed and compared across the same number of links.

####  3.2.2 Modularity by Lobes

Each network was assessed at the scale of its lobar organization – the
frontal, temporal, parietal and occipital divisions of the cortical surface.
Modularity index (Q), was calculated to determine whether cortical lobes as
conventionally defined correspond to modules of the CT/FC network. This can be
done by calculating the Q according to lobe (by employing vector of nodal
affiliation with the particular lobe as initial community affiliation vector).
The modularity index quantifies the observed fraction of within-module degree
values relative to those expected if connections were randomly distributed
across the network. In the context of this short communication, cortical lobes
have been considered and analyzed as network modules. Since the constructed
CT/FC network contains both positive and negative edge strengths, I used the
asymmetric generalization of the modularity quality function introduced in
Rubinov and Sporns [15]. Since the algorithm returns varying results on
different runs, to control for the variance in the analysis, modularity index
Q was calculated 100 times for each network, i.e., from running the same
analysis 100 times on CT/FC network. In other words, index Q was estimated
from 100 runs on CT/FT network and on 100 randomized (CT/FC each) versions of
these networks. (See section 4.2) for details.) For networks of similar size,
it is accepted that a Q value above 0.3 is a good indicator of the existence
of significant modules in a network [16].

##  4 Results

###  4.1 Network Density

As described in the Methods, a single value for network density was obtained
by considering percolation threshold for each network. This yielded a network
density of $\kappa$ = 0.16, which includes both positive and negative edges.
This choice of the threshold insures that each network (either CT or FC) has
the same ($N=362$) total number of positive and negative links.

Figure 2: Cortical thickness and functional correlations networks. (A) Group-
wise cortical thickness (CT) and functional correlations (FC) networks
arranged (from top to bottom) according to nodal affiliations with frontal,
temporal, parietal and occipital lobes. (B) Histograms of modularity index (Q)
by lobe calculated on 100 randomized CT (magenta) and FC (cyan) networks and
compared against the modularity index calculated on group-wise CT (magenta)
and FC (cyan) networks (values of 0.446 and 0.442, respectively). (C)
Histogram of the cortical distances for the convergent positive (red) and
negative (blue) networks. Vertical, color-coded lines represent mean value of
the corresponding cortical distances (i.e., distances between regions with the
overlapping positive or negative correlations). (D) Overlapping edges of
positive and negative CT and FC networks; red – positive correlations; blue –
negative correlations.

###  4.2 Modularity by Lobes

Modularity index of the CT and FC networks by lobes were 0.446 and 0.442 (mean
values on 100 runs for each network) respectively, thus indicating modular
organization of CT/FC networks by cortical lobes. To get the confidence
interval for network modularity distribution by lobar (frontal, temporal,
parietal and occipital) divisions, 200 randomized networks were generated on
CT/FC networks (100 each), i.e., each null network model was tested to obtain
distribution of the modularity indexes (Qs) on randomly distributed
correlations of the CT/FC network. Fig. 2(B) shows how the distributions of
100 Q values on random CT/FC networks differ from those calculated on real
networks (i.e., on 100 calculations of the modularity index for each network).

###  4.3 Convergent Cortical Thickness and Functional Networks Maps

Overlapping thickness and functional networks were calculated by element-wise
multiplication of the two matrices. Fig. 2(D) shows products of this
multiplication for the coupling of positive and negative
correlations/networks. The survived network edges of positive correlations map
intra-lobar and inter-hemispheric homologous (off-diagonal matrix elements).
In contrast, the coupled negative networks map exclusively onto inter-lobar
connections. In particular, the frontal lobe regions – rostral anterior
cingulate, rostral middle frontal and superior frontal – are correlated with
regions of the occipital lobe. The mean distance between the coupled networks
nodes ($N_{agree}$) was ($72\pm 30)mm$ for positive and ($81\pm 30)mm$ for
negative network nodes ($p>0.05$), Fig. 2(C).

##  5 Discussion

This study demonstrates convergence of the cortical thickness and functional
correlations networks, with reference to modular organization of the
conventional lobar divisions. Noticeably different patterns of the overlapping
correlations were observed for positive and negative networks. In contrast to
the positive networks’ overlaps, which map almost exclusively within-lobe
interactions, negative networks show only between-lobes interactions.

Both cortical networks – inter-regional thickness and functional correlations
– demonstrated modular organizations that mirror conventional frontal,
temporal, parietal and occipital divisions of the cortical surface. Thus
allowing for inferences about overlapping interactions of these networks by
cortical lobes. While several previous studies reported modular topology of
the cortex by cortical morphological (thickness and volume) [17, 18] or
functional [19] measurements, no data are available on functional modules as
an instinct property of the cortical surface lobar divisions. Motivated by the
lobe-specific patterns of atrophy progression in two types of dementia, we
have recently reported similar findings on cortical thickness and surface area
networks in healthy elderly subjects [6]. This study represents first effort
to map patterns of positive and negative correlations in cortical thickness
and functional networks in healthy adults.

Several interesting results observed when mapping overlapping links between
the cortical thickness and functional networks by lobar divisions need further
attention. First, the approach adopted here was to compare the networks at a
single value of network density. The advantage of this approach is that is
maps only the most pronounced correlations within the both networks. This
means that the correlations statistically not different from zero were not
considered in the analysis. However, the obtained topologies may include
different numbers of correlations for a larger study-group size. Thus, the
calculation of the maps of overlapping correlations across a range of networks
densities (based also on a larger study group) may reveal some of the cortical
interactions that did not pass statistical significance test in this study.
Second, the definition of conventional lobar divisions used in this study was
based on the Desikan-Kelliany cortical parcellation. Although well
established, there exist alternative approaches for the cortical segmentation
by using different brain templates [20]. Since adoption of different brain
templates may have influence on the patterns of network correlations, further
studies could validate these results using different anatomical (or
functional) classifications. Finally, a study combining structural, functional
and diffusion-weighted MRI data could test to what extent variations and
strengths of positive and negative correlations can be inferred from
underlying direct axonal links or synchronous/asynchronous functional links.
These future studies could also resolve problems of estimating networks
individually for each subject.

### Acknowledgements

The author would like to acknowledge the support of the Maxwell compute
cluster funded by the University of Aberdeen.

## References

  * [1] Evans, A. C.  Networks of anatomical covariance.  _Neuroimage_ 80, 489–504 (2013).
  * [2] Lerch, J. P. _et al._ Mapping anatomical correlations across cerebral cortex (macacc) using cortical thickness from mri.  _Neuroimage_ 31, 993–1003 (2006).
  * [3] Gong, G., He, Y., Chen, Z. J. & Evans, A. C.  Convergence and divergence of thickness correlations with diffusion connections across the human cerebral cortex.  _Neuroimage_ 59, 1239–1248 (2012).
  * [4] Hyde, K. L. _et al._ The effects of musical training on structural brain development: a longitudinal study.  _Annals of the New York Academy of Sciences_ 1169, 182–186 (2009).
  * [5] Witiuk, K. _et al._ Cognitive deterioration and functional compensation in als measured with fmri using an inhibitory task.  _Journal of Neuroscience_ 34, 14260–14271 (2014).
  * [6] Vuksnaović, V., Staff, R., Ahearn, T., Murray, A. & Claude, W.  Cortical thickness and surface area networks in healthy aging, alzheimer’s disease and behavioral variant fronto-temporal dementia.  _Int J Neur Sys_ (in press).
  * [7] Park, H.-J. & Friston, K.  Structural and functional brain networks: from connections to cognition.  _Science_ 342, 1238411 (2013).
  * [8] Kandel, E. R. _et al._ _Principles of neural science_ , vol. 4 (McGraw-hill New York, 2000).
  * [9] Mesulam, M.-M.  From sensation to cognition.  _Brain: a journal of neurology_ 121, 1013–1052 (1998).
  * [10] Desikan, R. S. _et al._ An automated labeling system for subdividing the human cerebral cortex on mri scans into gyral based regions of interest.  _Neuroimage_ 31, 968–980 (2006).
  * [11] Vuksanović, V. & Hövel, P.  Functional connectivity of distant cortical regions: role of remote synchronization and symmetry in interactions.  _NeuroImage_ 97, 1–8 (2014).
  * [12] Ciric, R. _et al._ Benchmarking of participant-level confound regression strategies for the control of motion artifact in studies of functional connectivity.  _Neuroimage_ 154, 174–187 (2017).
  * [13] Vuksanović, V. & Hövel, P.  Dynamic changes in network synchrony reveal resting-state functional networks.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 25, 023116 (2015).
  * [14] van Wijk, B. C., Stam, C. J. & Daffertshofer, A.  Comparing brain networks of different size and connectivity density using graph theory.  _PLoS One_ 5, e13701 (2010).
  * [15] Rubinov, M. & Sporns, O.  Complex network measures of brain connectivity: uses and interpretations.  _NeuroImage_ 52, 1059–1069 (2010).
  * [16] Clauset, A., Newman, M. E. & Moore, C.  Finding community structure in very large networks.  _Physical review E_ 70, 066111 (2004).
  * [17] Chen, Z. J., He, Y., Rosa-Neto, P., Germann, J. & Evans, A. C.  Revealing modular architecture of human brain structural networks by using cortical thickness from mri.  _Cerebral cortex_ 18, 2374–2381 (2008).
  * [18] Bassett, D. S. _et al._ Hierarchical organization of human cortical networks in health and schizophrenia.  _Journal of Neuroscience_ 28, 9239–9248 (2008).
  * [19] Meunier, D., Achard, S., Morcom, A. & Bullmore, E.  Age-related changes in modular organization of human brain functional networks.  _Neuroimage_ 44, 715–723 (2009).
  * [20] Tzourio-Mazoyer, B. _et al._ Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain.  _NeuroImage_ 15, 273–289 (2002).