#  CoLoRe: fast cosmological realisations over large volumes with multiple
tracers

César Ramírez-Pérez  Javier Sanchez  David Alonso  Andreu Font-Ribera

###### Abstract

We present CoLoRe, a public software package to efficiently generate synthetic
realisations of multiple cosmological surveys. CoLoRe can simulate the growth
of structure with different degrees of accuracy, with the current
implementation supporting lognormal fields, first, and second order Lagrangian
perturbation theory. CoLoRe simulates the density field on an all-sky light-
cone up to a desired maximum redshift, and uses it to generate multiple 2D and
3D maps: galaxy positions and velocities, lensing (shear, magnification,
convergence), integrated Sachs-Wolfe effect, line intensity mapping, and line
of sight skewers for simulations of the Lyman-$\alpha$ forest. We test the
accuracy of the simulated maps against analytical theoretical predictions, and
showcase its performance with a multi-survey simulation including DESI
galaxies and quasars, LSST galaxies and lensing, and SKA intensity mapping and
radio galaxies. We expect CoLoRe to be particularly useful in studies aiming
to characterise the impact of systematics in multi-experiment analyses,
quantify the covariance between different datasets, and test cross-correlation
pipelines for near-future surveys.

FERMILAB-PUB-21-562-PPD

##  1 Introduction

Ongoing and future cosmological surveys will explore large volumes with
multiple tracers to study dark energy, inflation and massive neutrinos.
Spectroscopic surveys such as the Dark Energy Spectroscopic Instrument (DESI,
[1]), Euclid [2] and Roman [3] will collect tens of millions of precise galaxy
redshifts. Meanwhile, the Vera Rubin Observatory’s Legacy Survey of Space and
Time (LSST, [4]) will photometrically observe billions of galaxies and provide
an exquisite weak lensing map over a large fraction of the sky. Lensing maps
will also be obtained from future experiments observing the cosmic microwave
background (CMB), such as the Simons Observatory (SO, [5]) and CMB-S4 [6].
Finally, large catalogs of radio galaxies and 21cm intensity maps will be
provided by experiments such as the Square Kilometer Array (SKA, [7]) or HIRAX
[8].

In order to obtain robust cosmological constraints from these large and
complex datasets, it is important to be able to efficiently generate mocks,
synthetic realisations of the data. For instance, mocks are often used to
compare survey strategies, to test the analysis pipeline, to study the impact
of astrophysical or instrumental contaminants, and to estimate the covariance
of the measurements.

When generating mocks one needs to trade off realism for computing costs. It
is just not feasible to generate hundreds or thousands of N-body simulations
covering volumes of tens of cubic gigaparsecs, let alone hydrodynamic
simulations that could model baryonic effects. On the other side of the
spectrum, lognormal realisations [9] offer an efficient way of obtaining
simplified mock catalogs with the correct distribution only on large, linear
scales (see for instance [10, 11, 12]). Lagrangian Perturbation Theory (LPT,
[13]) has inspired several approximated methods that can reproduce the
distribution of matter on intermediate, mildly non-linear scales. These
include PTHalos [14, 15, 16], Pinocchio [17], COLA [18], QPM [19], PATCHY [20]
ICE-COLA [21], L-PICOLA [22], HALOGEN [23] and EZMocks [24].

The best cosmological inference will come from joint analyses of multiple
cosmological probes, each providing independent and complementary information.
Two of the most important challenges in these joint analyses will be
characterising the effects of systematics affecting several experiments, and
estimating the cross-covariance between the different 2D and 3D datasets, with
partially overlapping area and redshift ranges. This publication addresses the
need of simultaneously simulating these surveys in a coherent and efficient
framework.

We present CoLoRe (Cosmological Lofty Realization), a parallel C code for
generating fast mock realizations of multiple cosmological surveys111The code
is publicly available at https://github.com/damonge/CoLoRe.. CoLoRe can
simulate the growth of structure using either a lognormal model or LPT (at 1st
or 2nd order), and it can simulate a plethora of cosmological tracers:
photometric and spectroscopic galaxies, weak lensing, intensity mapping,
Integrated Sachs-Wolfe effect and CMB lensing, or the Lyman-$\alpha$ forest in
the spectra of high-redshift quasars. It has been designed in a highly modular
fashion, making it easy to add new tracers or more complex models of growth of
structure. It uses both OpenMP and MPI parallelisation, and it is specially
suited to run with multiple nodes in high performance computing facilities.

The paper is organised as follows. In Section 2 we describe in detail CoLoRe,
its code structure, the cosmological assumptions made, and the list of tracers
already available. The validation of CoLoRe to generate reliable mocks for
intensity mapping and for Lyman-$\alpha$ forest studies was already presented
in previous work [25, 26]. In Section 3 we validate its ability to simulate
galaxy clustering and weak lensing statistics, and discuss the computing and
memory requirements to run large simulations. Finally, in Section 4 we draw
the conclusions.

##  2 Methods

###  2.1 Overall code structure

CoLoRe is written in a modular way that makes it relatively straightforward to
modify (e.g. to add a new non-linear structure formation model, or a new
tracer of the density fluctuations). In a standard run, CoLoRe goes through
the steps listed below, each of which is associated with a compartmentalised
piece of code:

  1. 1.

Initialisation. CoLoRe interprets the configuration file, allocates the
resources needed to carry out the requested simulation, and initialises a
number of cosmological quantities (redshift-distance relation, growth history,
linear matter power spectrum, background densities of all source tracers).

  2. 2.

Predictions: CoLoRe produces theoretical predictions for the three-dimensional
power spectrum of all biased matter tracers in the lognormal approximation.
This is mostly useful when using the lognormal structure formation model.

  3. 3.

Gaussian random fields. Two three-dimensional Cartesian grids are generated
containing the linear matter overdensity $\delta^{L}_{M}({\bf x})$ and the
Newtonian gravitational potential $\phi_{N}({\bf x})$ at redshift $z=0$. The
grid is sufficiently large to hold a sphere of comoving radius $\chi(z_{\rm
max})$, where $z_{\rm max}$ is the maximum redshift of the run. The spatial
resolution of the simulation is set by $N_{\rm grid}$, the number of grid
cells into which the box is divided in each dimension. The grid cell size is
therefore approximately $\Delta x=L_{\rm box}/N_{\rm grid}\simeq 2\chi(z_{\rm
max})/N_{\rm grid}$.

  4. 4.

Physical density field. The Gaussian overdensity $\delta^{L}_{M}({\bf x})$ is
transformed into a non-linear, physical overdensity field $\delta_{M}({\bf
x})$ through one of the structure formation models supported by CoLoRe. This
is done in the lightcone (i.e. the value of the field at comoving position
${\bf x}$ is $\delta_{M}(t(|{\bf x}|),{\bf x})$, where $t(\chi)$ is the cosmic
time at comoving distance $\chi$), with the observer located at the center of
the Cartesian box. The physical density field is such that $\delta\geq-1$
everywhere. The gravitational potential is also evolved in the lightcone
assuming linear growth.

  5. 5.

Density normalisation. CoLoRe uses non-linear transformations to generate
biased tracers of the matter overdensity. In general, these can be written as

$1+\delta_{k}=\frac{B_{k}(\delta_{M})}{\langle B_{k}(\delta_{M})\rangle},$
(2.1)

where $B_{k}$ is the non-linear biasing relation for tracer $k$. Due to the
non-linearity of these relations, the ergodic average in the denominator of
the previous equation is not necessarily equal to 1 (even if
$\langle\delta_{M}\rangle=0$), and therefore the normalising factor ensures
that $\langle\delta_{k}\rangle=0$ for all biased tracers. Since densities are
defined on the lightone, the normalisation factors are computed at this stage
independently for several redshift shells by averaging over grid cells.

  6. 6.

Get Cartesian information. At this stage the overdensity and Newtonian
potential grids are distributed across computer nodes as slabs of equal width
$N_{\rm slab}=N_{\rm grid}/N_{\rm nodes}$. Before proceeding further, CoLoRe
collects all information available in these slabs and needed by each of the
tracers requested for this simulation. This involves any data product not
involving line-of-sight integrals or interpolations, which are dealt with in
later stages. For instance, this is when source catalogs are generated by
Poisson-sampling the biased density field. All tracers are endowed with a
method tracer_set_cartesian that collects this information.

  7. 7.

Redistribution into beams. In order to carry out line-of-sight integrals and
interpolations, CoLoRe redistributes tracer data so each node has access to
all the data in a set of sky regions, labelled “beams”, covering the full
range of redshifts $0\leq z\leq z_{\rm max}$. Each beam is defined using the
HEALPix pixellation scheme [27], as the region of the celestial sphere covered
by a given low resolution pixel. The HEALPix $N_{\rm side}$ resolution
parameter used to define these beams is chosen to be large enough that the
full dataset is approximately evenly split between computer nodes. All tracers
have an associated method tracer_distribute in charge of distributing the
tracer data in each node’s slab to all other nodes whose beams intersect with
it. Note that, although the tracer information is now distributed across nodes
through beams, the density and Newtonian potential grids are still distributed
in slabs.

  8. 8.

Get beam information. Any calculation involving a line-of-sight integral (e.g.
gravitational lensing) or interpolation (e.g. Lyman-$\alpha$ skewers) is done
after the tracers have been redistributed into beams. The calculation is done
in three stages:

    1. (a)

Preprocessing. Initialisation of any necessary quantities (e.g. setting all
variables that eventually hold a gravitational lensing calculation to zero).
Each tracer has an associated function tracer_beams_preproc in charge of doing
this.

    2. (b)

Loop through slabs. All nodes send their current slab of the density and
Newtonian potential grids to the node on their right (assuming periodic
boundary conditions). Once the new slab is received, each node gathers the
necessary information from it (e.g. the contribution to a lensing convergence
integral from the section of the Newtonian potential held in that slab) and
adds it to each tracer222Note that all quantities calculated at this stage
(e.g. integrals and interpolations) are linear and additive on $\delta_{M}$
and $\phi_{N}$.. CoLoRe carries out this calculation through a method
tracer_get_beam_properties associated with each tracer. This is repeated
$N_{\rm nodes}$ times, at which point all nodes have had access to the full
density and potential grids.

    3. (c)

Post-processing. Each tracer finishes off any calculation still needed after
having gathered all information in the preceding step (e.g. multiplying maps
by an overall normalization factor). This is done by a method
tracer_beams_postproc associated with each tracer.

  9. 9.

Write output. Each tracer writes all its information (e.g. in the form of maps
or catalogs) to file through a method write_tracer. CoLoRe uses the FITS
standard in most cases, although it is also possible to save source catalog
data as ASCII or HDF5 files.

Modifying CoLoRe to support a new structure formation model would involve
implementing it as part of step 4 above, with no effect on the rest of the
code. Adding a new type of tracer would involve creating the corresponding
tracer methods for it enumerated above (_set_cartesian, _distribute,
_beams_preproc, _get_beam_properties, _beams_postproc, and write_).

The assumptions made by CoLoRe to compute the background cosmological
quantities (step 1) are described in Section 2.2. Section 2.3 describes the
Gaussian density fields and the different non-linear structure formation
models supported by CoLoRe (steps 3 and 4). Section 2.4 describes in detail
the calculations carried out for each of the tracers, including the bias
models supported (steps 5-8). Finally, the theory predictions computed by
CoLoRe for lognormal fields (step 2) are discussed in Section 2.5.

###  2.2 Cosmological assumptions

When generating simulated observations, CoLoRe makes a number of assumptions
about the underlying cosmological model. We describe these here.

CoLoRe assumes a flat $w$CDM cosmological background, characterised, at low
redshifts, by 3 cosmological parameters: the background matter density
$\Omega_{M}$, the current expansion rate $H_{0}$, and a constant dark energy
equation of state parameter $w$. The expansion rate is thus given by

$H(z)=H_{0}\,\left[\Omega_{M}(1+z)^{3}+(1-\Omega_{M})(1+z)^{3(1+w)}\right],$
(2.2)

in terms of which the comoving distance is333Note that we use units with $c=1$
throughout.

$\chi(z)=\int_{0}^{z}\frac{dz^{\prime}}{H(z^{\prime})}.$ (2.3)

Matter density perturbations are governed by a linear matter power spectrum at
$z=0$, $P_{0}(k)$, which must be provided to CoLoRe on input, and is then
normalised to the chosen value of $\sigma_{8}$. If the power spectrum is
needed on scales larger than those provided, it is extrapolated assuming a
power-law behaviour $P_{0}(k)\propto k^{n_{s}}$ on small $k$, where $n_{s}$ is
the scalar spectral index (also provided on input).

Finally, CoLoRe assumes a self-similar growth for the linear matter
overdensity: $\delta_{M}^{L}({\bf x},z)=\delta_{M}^{L}({\bf x},0)D(z)$, where
$D(z)$ is the linear growth factor. $D(z)$ is calculated from the cosmological
parameters by solving the differential equation

$\frac{d}{da}\left(a^{3}\,H(a)\frac{dD}{da}\right)=\frac{3}{2}\Omega_{M}(a)\,aH(a)\,D(a),$
(2.4)

where $a=1/(1+z)$ is the scale factor.

Although internally CoLoRe uses “$h$-inverse” units (i.e. distances are given
in units of ${\rm Mpc}\,h^{-1}$), all simulation outputs involve observable
quantities (redshift and angles), and therefore are insensitive to this
choice.

###  2.3 Matter box

The first step after initialising the cosmological model in a standard CoLoRe
run is the generation of a Gaussian realisation of the linear matter
inhomogeneities $\delta_{M}^{L}$ at $z=0$ on a Cartesian cubic grid. This is
done by drawing the Fourier coefficients of $\delta_{M}^{L}$ as independent
Gaussian random numbers from the input linear matter power spectrum using the
Box-Muller transform with variance:

$\sigma^{2}({\bf k})=\frac{P_{0}(k)}{(\Delta k)^{3}},$ (2.5)

where $\Delta k\equiv 2\pi/L_{\rm box}$ is the sampling rate in Fourier space.
CoLoRe can alternatively apply a Gaussian smoothing kernel with scale $R_{G}$
to the linear power spectrum when generating the linear Fourier coefficients.
This may be useful to control the behaviour of the non-linear overdensity
field (see discussion in Section 2.5).

At the same time, CoLoRe populates a similar cartesian grid with the values of
the Newtonian gravitational potential $\phi_{N}({\bf x})$, related to the
matter inhomogeneities in Fourier space via:

$\phi_{N}({\bf
k})=-\frac{3}{2}H_{0}^{2}\,\Omega_{M}\,\frac{\delta_{M}^{L}({\bf k})}{k^{2}}.$
(2.6)

The linear matter overdensity thus generated is then transformed into a
physical (i.e. positive-definite) non-linear matter overdensity in the
lightcone using one of the three structure formation models currently
supported by CoLoRe, which we describe below.

####  2.3.1 Lognormal fields

Lognormal fields were first proposed and analysed by [9] as a possible way to
describe the distribution of matter in the Universe. A lognormal random field
$x_{\rm LN}$ is defined in terms of a Gaussian random field $x_{\rm G}$
through the local transformation

$x_{\rm LN}=\exp x_{\rm G}.$ (2.7)

One of the nice properties of these fields is that, while the Gaussian
variable $x_{\rm G}$ is allowed to take any values in $(-\infty,+\infty)$,
$x_{\rm LN}$ can only take positive values by construction. Furthermore, as
discussed in [9], the density field evolved along Lagrangian trajectories
according to the linear velocity field along can be well described by a
lognormal distribution, which justifies the use of lognormal fields from a
physical point of view. In order to obtain a lognormal overdensity field with
zero mean from a Gaussian field, the transformation (2.7) must be slightly
varied as follows:

$1+\delta_{\rm LN}=\exp\left(\delta_{G}-\frac{\sigma_{G}^{2}}{2}\right),$
(2.8)

where $\sigma_{G}$ is the variance of the Gaussian overdensity field.

Lognormal density fields have been used in the past by different
collaborations in order to perform fast galaxy mock realisations [28, 29, 30,
31, 32], and are, therefore, a well established tool. Since the lognormal
transformation is a simple, local modification of the density field, it is by
far the fastest and most memory-efficient structure formation model
implemented in CoLoRe.

However, the simplicity of the lognormal transformation implies that lognormal
fields cannot be expected to describe all higher-order correlators of the
density field (e.g. bispectra), to give rise to filamentary structure, or to
reproduce the small-scale properties of the density field correctly [33], and
therefore this kind of mock realisations have a limited applicability.

Within this framework, the non-linear overdensity is generated in CoLoRe by
applying Eq. 2.8 to the linear overdensity field evolved in the lightcone
assuming linear growth (which is applied to both $\delta_{G}$ and $\sigma_{G}$
in this equation).

Caution must be exercised when making use of lognormal fields. Since the
lognormal transformation involves the exponentiation of a Gaussian field,
large $(\gg 1)$ values of $\delta_{G}$ lead to much larger fluctuations in
$\delta_{\rm LN}$. Thus, if the amplitude of $\delta_{G}$, characterised by
its standard deviation $\sigma_{G}$, is large, the resulting lognormal field
will exhibit a large degree of inhomogeneity, with an enormous variance
dominated by extreme fluctuations in a small number of voxels. This behaviour
can be avoided through the use of the Gaussian smoothing kernel described
above. This modifies the linear and non-linear power spectra in an
analytically predictable manner.

Note that other common implementations of the lognormal transformation to
generate mock cosmological realisations (e.g. [29]) have employed a different
method to avoid this problem. Instead of using the input power spectrum to
generate $\delta_{G}$ and then transform it into $\delta_{\rm LN}$, the input
power spectrum is taken to be that of the final $\delta_{\rm LN}$. The inverse
lognormal transformation is thus applied to the input power spectrum at the
level of the two-point correlator (see Section 2.5), to obtain the power
spectrum with which the Gaussian field is generated. CoLoRe does not
explicitly support this method, since it runs contrary to the idea of treating
the lognormal transformation as a non-linear structure formation model. It
would be, however, possible to pass as input to CoLoRe the “Gaussianised”
power spectrum in order to obtain the desired matter power spectrum at $z=0$
in the non-linear matter fluctuations.

####  2.3.2 Lagrangian perturbation theory

Lagrangian perturbation theory (LPT) [13] provides an alternative fast method
to generate non-linear physical matter overdensities, which has been used in
the past to generate mock galaxy catalogs [15, 24, 16]. The beauty of LPT lies
in its ability to capture non-linear aspects of the Eulerian matter
overdensity (e.g. the formation of filamentary structure) by carrying out low-
order perturbation theory calculations in the displacement field. The basic
premise is as follows: the linear overdensity is used to predict the
Lagrangian displacement of a set of massive test particles with respect to
their original unperturbed positions. The non-linear density field is then
given by the density of the displaced test particles which is, by definition,
positive-definite. CoLoRe supports the generation of Lagrangian displacements
at first and second order in perturbation theory. These first and second-order
displacements at $z=0$ are scaled with the corresponding growth factors before
interpolating the test particles onto a grid. We provide a brief overview of
LPT here, and direct the reader to [13] for further details.

Let ${\bf x}(t)=a(t)\,\left[{\bf q}+\bm{\Psi}({\bf q},t)\right]$ be the
physical position of a particle starting at comoving coordinates ${\bf q}$.
$\bm{\Psi}({\bf q},t)$ is the so-called Lagrangian displacement vector. In the
Newtonian approximation, the motion of these particles is governed by Newton’s
second law, which is sourced by the gravitational potential caused by the
particles themselves. This leads to two coupled equations that can be
summarized into a single equation for the divergence of $\bm{\Psi}$:

$J\left({\sf
J}^{-1}\nabla\right)\cdot\left(\bm{\Psi}^{\prime\prime}+a\,H\,\bm{\Psi}^{\prime}\right)=\frac{3}{2}a^{2}H^{2}\,\Omega_{M}(J-1),$
(2.9)

where all derivatives are taken with respect to conformal time $d\tau\equiv
dt/a$, ${\sf J}_{ij}\equiv\delta_{ij}+\partial_{i}\Psi_{j}$ is the Jacobian of
the Lagrangian flow, and $J\equiv{\rm det}({\sf J})$. Note that the matter
overdensity is given by $1+\delta=J^{-1}$.

At second order in the displacement field, and discarding all curl-like
components of $\bm{\Psi}$, the solution is given by

$\bm{\Psi}({\bf q},a)=D(a)\,\nabla\varphi_{\rm LPT}^{(1)}({\bf
q})+D^{(2)}(a)\,\nabla\varphi_{\rm LPT}^{(2)}({\bf q}).$ (2.10)

Here $D(a)$ is the linear growth factor, satisfying Eq. 2.4, $D^{(2)}$ is the
second-order growth factor, satisfying

$\frac{d}{da}\left(a^{3}\,H(a)\frac{dD^{(2)}}{da}\right)=\frac{3}{2}\Omega_{M}(a)\,aH(a)\,\left(D^{(2)}(a)-[D(a)]^{2}\right),$
(2.11)

and $\varphi_{\rm LPT}^{(1,2)}$ are the first- and second-order LPT
potentials, given by

$\displaystyle\nabla^{2}\varphi_{\rm LPT}^{(1)}=-\delta_{\rm M}^{L},$ (2.12)
$\displaystyle\nabla^{2}\varphi_{\rm
LPT}^{(2)}=\frac{1}{2}\sum_{ij}\left[\partial_{i}^{2}\varphi_{\rm
LPT}^{(1)}\partial_{j}^{2}\varphi_{\rm
LPT}^{(1)}-\left(\partial_{i}\partial_{j}\varphi_{\rm
LPT}^{(1)}\right)^{2}\right].$ (2.13)

CoLoRe solves the LPT Poisson equations (2.12 and 2.13) in Fourier space using
the Gaussian overdensity field as input, and computes the second-order growth
factor using the approximation [13]

$D_{2}(a)=-\frac{3}{7}\left[D(a)\right]^{2}\left[\Omega_{M}(a)\right]^{-1/143}.$
(2.14)

Once the displacement vector has been calculated and applied to a set of test
particles initially located at the centers of the Cartesian grid cells, the
density field is calculated by interpolating the displaced positions onto the
grid444To do this, CoLoRe supports three standard mass-assignment methods:
nearest-grid-point (NGP), cloud-in-cell (CIC), and triangular-shaped cloud
(TSC) [34]..

LPT is able to generate a more realistic non-linear density field than the
lognormal model at the cost of significantly higher memory requirements and
longer computation times. While generating the lognormal overdensity does not
require additional resources beyond those used to generate the
$\delta_{M}^{L}$ and $\phi_{N}$ grids, generating the first-order displacement
requires three additional Cartesian grids to hold the components of
$\bm{\Psi}$, and second-order LPT demands an additional 5 grids to store the
Hessian of the first-order LPT potential (the array holding the Gaussian
density field can be reused to store one of the 6 independent components of
the Hessian). The number of fast Fourier transforms needed, which dominate the
total computation time, is also different in each case: none in the case of
the lognormal model, 4 in the case of first-order LPT, and 13 for second-order
LPT.

###  2.4 Tracers

CoLoRe is able to generate simulated observations for a variety of
cosmological tracers of the same underlying matter fluctuations. The details
of the calculations involved in each of the supported tracer types, carried
out in steps 5 to 8 of the procedure outlined in Section 2.1, is discussed in
detail here.

####  2.4.1 Projected maps

Figure 1: Bias models implemented in CoLoRe. The exponential model preserves
the “lognormality” of the field if using the lognormal structure formation
model, but it can lead to numerically unstable results in the presence of
sufficiently large fluctuations in the Gaussian field. The exp-truncated model
can be used to curb this behaviour.

Figure 2: Simulated maps of the anisotropic stochastic gravitational wave
background from astrophysical sources at redshifts $z<0.4$ using the models of
[35]. The top and bottom plots show simulations using the lognormal and first-
order LPT sructure formation models respectively. The former is characterised
by strong positive fluctuations on a few regions, while the latter displays
the more physical filamentary structure of the cosmic web.

The simplest tracer supported by CoLoRe is the so-called “custom projected
tracer”. The associated observable is the overdensity in a biased tracer of
the matter fluctuations projected onto the celestial sphere after integrating
over an arbitrary radial kernel:

$\Delta_{W}(\hat{\bf n})=\int dz\,W(z)\,\delta_{W}(z,\chi(z)\hat{\bf n}),$
(2.15)

where the integral is over redshift $z$, $W(z)$ is the tracer’s radial kernel,
and $\delta_{W}$ are the three-dimensional fluctuations in the tracer, related
to the matter fluctuations via Eq. 2.1.

CoLoRe supports the following three local bias models, although more can be
easily added to the code:

Exponential: $\displaystyle B_{k}(\delta_{M})=(1+\delta_{M})^{b_{k}},$ (2.16)
Truncated: $\displaystyle B_{k}(\delta_{M})={\rm Max}(1+b_{k}\delta_{M},0),$
(2.17) Exp-truncated: $\displaystyle
B_{k}(\delta_{M})=\left\\{\begin{array}[]{cc}\exp\left[b_{k}\delta_{M}/(1+\delta_{M})\right]&\delta_{M}\leq
0\\\ 1+b_{k}\delta_{M}&\delta_{M}>0\end{array}\right..$ (2.20)

These three models are shown in Figure 1 for $b_{k}=2$, and are designed to be
positive definite ($B_{k}(\delta_{M})>0$ for $\delta_{M}\in[-1,\infty)$), and
to reduce to a linear biasing relation ($B_{k}(\delta_{M})\simeq
1+b_{k}\delta_{M}$) for small fluctuations. Given the modular nature of
CoLoRe, different biasing models could be easily added by the user.

The integral in Eq. 2.15 is calculated as follows: for each sky pixel, an
imaginary line of sight connecting it with the observer at the box center is
subdivided into intervals of constant comoving distance $\Delta\chi$,
commensurate with the Cartesian cell size $\Delta x$. The value of
$\delta_{M}$ at the center of each interval is calculated from the Cartesian
grid using trilinear interpolation, and is then translated into the
corresponding $\delta_{W}$. The integral is then calculated as a sum over all
intervals along the line of sight.

Custom projected tracers can be used in CoLoRe to make simulated maps of a
wide variety of two-dimensional cosmological anisotropic observables that
correlate with the large-scale structure. As an example, Figure 2 shows
simulated maps of the anisotropic stochastic gravitational wave background
from astrophysical sources at $z<0.4$ according to the models of [35]. The
figure shows results for the lognormal and first-order LPT structure formation
models, showcasing the morphological differences between them.

####  2.4.2 Integrated Sachs-Wolfe effect

The time evolution in the gravitational potential at late times due to the
accelerated background expansion causes an energy loss or gain in a background
of photons which correlates with the large-scale structure. This is the so-
called integrated Sachs-Wolfe effect (ISW, [36]).

The fluctuation in the temperature of a background of photons with black-body
spectrum emitted at redshift $z_{*}$ is given by

$\left.\frac{\Delta T}{T}\right|_{\rm ISW}(\hat{\bf
n})=2\int_{0}^{z_{*}}dz\frac{\dot{\phi}_{N}(z,\chi(z)\hat{\bf
n})}{(1+z)\,H(z)}.$ (2.21)

Assuming linear growth, appropriate on the large scales on which the ISW is
relevant, one can approximate $\dot{\phi}_{N}(z)=H(z)[f(z)-1]\phi_{N}(z)$,
where $f\equiv d\log D/d\log a$ is the growth rate.

The ISW tracer is therefore equivalent to the custom projected tracer
described in the previous section with the Newtonian potential $\phi_{N}$
taking the role of $\delta_{W}$, and with a kernel

$W_{\rm ISW}(z)=\frac{f(z)-1}{1+z}\Theta(z<z_{*}),$ (2.22)

where $\Theta$ is the Heaviside function. Thus, the same numerical methods
described in Section 2.4.1 are used by CoLoRe to generate simulated ISW maps.
The top panel of Figure 3 shows an example of a simulated map of the ISW
effect for $z_{*}=0.5$, characterised by features on very large scales due to
the $1/k^{2}$ relation between gravitational potential and matter overdensity.

####  2.4.3 Gravitational lensing

Figure 3: Top: simulated map of the low-redshift ISW effect for a source plane
at $z_{*}=0.5$. Bottom: map of the lensing convergence for a source plane at
$z_{*}=1.4$.

The fluctuations in the gravitational potential perturb the trajectories of
photons via gravitational lensing. The effects of gravitational lensing on
different cosmological observables are encoded in the so-called “lensing
potential”, defined as [37, 38]

$\psi(\hat{\bf
n})\equiv-2\int_{0}^{\chi_{*}}d\chi\frac{\chi_{*}-\chi}{\chi_{*}\chi}\phi_{N}(z(\chi),\chi\hat{\bf
n}).$ (2.23)

The trajectories of photons are deflected by an angle
$\bm{\alpha}\equiv\nabla_{\hat{\bf n}}\psi$, and the shapes of background
objects are distorted via the lensing distortion tensor
$\mathsf{\Gamma}\equiv-{\sf H}_{\hat{\bf n}}\psi$, where $\nabla_{\hat{\bf
n}}$ and ${\sf H}_{\hat{\bf n}}$ are the gradient and Hessian operators on the
sphere. The distortion tensor $\mathsf{\Gamma}$ is commonly decomposed into
its spin-0 and spin-2 components, the convergence $\kappa$ and shear
$(\gamma_{1},\gamma_{2})$:

$\mathsf{\Gamma}\equiv\left(\begin{array}[]{cc}\kappa+\gamma_{1}&\gamma_{2}\\\
-\gamma_{2}&\kappa-\gamma_{1}\end{array}\right).$ (2.24)

Computing $\bm{\alpha}$ and $\mathsf{\Gamma}$ from $\psi$ would require first
creating a map of $\psi$ to then differentiate. This is not feasible when
calculating the effects of lensing on a large number of sources at different
redshifts (see Section 2.4.4). Instead, we rewrite these quantities as

$\bm{\alpha}(\hat{\bf
n})=-2\int_{0}^{\chi_{*}}d\chi\frac{\chi_{*}-\chi}{\chi_{*}}\nabla_{\perp}\phi_{N}(z,\chi\hat{\bf
n}),\hskip 12.0pt\mathsf{\Gamma}(\hat{\bf
n})=-2\int_{0}^{\chi_{*}}d\chi\chi\frac{\chi_{*}-\chi}{\chi_{*}}{\sf
H}_{\perp}\phi_{N}(z,\chi\hat{\bf n}),$ (2.25)

where $\nabla_{\perp}$ and ${\sf H}_{\perp}$ are the gradient and Hessian
operators projected onto the plane perpendicular to $\hat{\bf n}$.

Explicitly, if $\hat{\bf
n}\equiv(\sin\theta\,\cos\varphi,\sin\theta\,\sin\varphi,\cos\theta)$,
defining the projector

${\sf P}_{\hat{\bf
n}}\equiv\left(\begin{array}[]{ccc}\cos\theta\,\cos\varphi&\cos\theta\,\sin\varphi&-\sin\theta\\\
-\sin\varphi&\cos\varphi&0\end{array}\right),$ (2.26)

the projected gradient and Hessian are

$\nabla_{\perp}\phi_{N}\equiv{\sf P}_{\hat{\bf n}}\nabla\phi_{N},\hskip
12.0pt{\sf H}_{\perp}\phi={\sf P}_{\hat{\bf n}}({\sf H}\phi_{N}){\sf
P}_{\hat{\bf n}}^{T},$ (2.27)

where $\nabla_{i}\equiv\partial/\partial x^{i}$ and ${\sf
H}_{ij}\equiv\partial^{2}/\partial x^{i}\partial x^{j}$.

The calculation of lensing-related quantities in CoLoRe is thus analogous to
the procedure outlined in Section 2.4.1: along a given line of sight $\hat{\bf
n}$ (corresponding to a map pixel or to the position of a given source), the
values of the first or second-order derivatives of $\phi_{N}$ are calculated
and interpolated from the Cartesian grid onto a set of equidistant points
along $\hat{\bf n}$. The corresponding quantities are then projected onto the
plane perpendicular to $\hat{\bf n}$, and the integrals in Eq. 2.25 are
computed as direct sums over the evaluated points.

Besides providing lensing information associated with its source catalogs (see
Section 2.4.4), CoLoRe returns maps of the lensing convergence
$\kappa(\hat{\bf n})$ for an arbitrary number of source planes at different
redshifts. An example at $z_{*}=1.4$ is shown in the bottom panel of Fig. 3.

####  2.4.4 Sources

Figure 4: Maps of source-related quantities for a simulated CoLoRe catalog in
the redshift range $z<0.3$. Top left: source overdensity. Top right: mean
redshift distortion. Middle: mean lensing shear. Bottom: mean lensing
displacement vector. Figure 5: Density (top) and radial velocity (bottom)
skewers for two arbitrary sources in a simulated CoLoRe catalog.

CoLoRe can also generate catalogs of discrete sources as biased tracers of the
matter distribution. The source distribution is modelled as a Cox process: the
number of sources of type $a$ in a given Cartesian cell $i$, $N_{i}^{a}$ is a
random Poisson variable with a stochastic mean given by

$\bar{N}_{i}^{a}=(\Delta x)^{3}\bar{n}_{a}(\chi)\frac{B_{a}(\delta_{M}({\bf
x}_{i}))}{\langle B_{a}(\delta_{M})\rangle},$ (2.28)

where ${\bf x}_{i}$ are the coordinates of cell $i$, $(\Delta x)^{3}$ is the
cell volume, $\bar{n}_{a}(\chi)$ is the redshift-dependent mean density of
sources, and $B_{a}$ is the biasing relation of type-$a$ sources (see Eqs.
2.16-2.20).

Once $N_{i}^{a}$ has been determined for each cell, the corresponding number
of sources are generated with three-dimensional positions randomly located
within their cell. The three-dimensional comoving position is then translated
into observable angular and redshift coordinates. In addition to the
cosmological redshift, the contribution from the source’s peculiar velocity,
necessary to simulate redshift-space distortions, is calculated as $\Delta
z_{\rm RSD}=v_{r}$, where the radial velocity is calculated from the gradient
of the gravitational potential as

$v_{r}(z,{\bf x})=-\frac{2}{3H_{0}^{2}\Omega_{M}}f(z)\,(\hat{\bf
n}\cdot\nabla)\phi_{N}(z,{\bf x}).$ (2.29)

Given this derivation of velocities they will only include linear effects.

If desired, the source catalogs can also contain gravitational lensing
information ($\alpha_{\theta}$, $\alpha_{\varphi}$, $\kappa$, $\gamma_{1}$,
$\gamma_{2}$) for each source. This can be used to construct a weak lensing
shear catalog with ellipticities $e_{i}=\gamma_{i}$, or to include the effects
of lensing magnifications by perturbing the source angular positions
$\theta\rightarrow\theta+\alpha_{\theta}$,
$\varphi\rightarrow\varphi+\alpha_{\varphi}$, and its flux by $F\rightarrow
F(1+2\kappa)$. In this case, all lensing quantities are calculated as
described in Section 2.4.3, by integrating the interpolated transverse
derivatives of the gravitational potential along each source’s line of sight.

Source catalogs can also be endowed with so-called “line-of-sight skewers”,
containing the matter overdensity and radial velocity interpolated from the
Cartesian box onto each source’s line of sight. As in the case of
gravitational lensing calculations, trilinear interpolation is used, and both
fields are sampled at radial comoving intervals equal to CoLoRe’s Cartesian
cell size. The use of skewers to produce simulated observations of the
Lyman-$\alpha$ forest was discussed in detail in [26].

It is worth noting that the computation of any quantity requiring full line-
of-sight information (lensing or skewers) requires redistributing sources from
“slabs” into “beams” across nodes, and for all nodes to loop through the full
Cartesian box to add the contribution of all slabs to their beams (i.e. points
7 and 8 in Section 2.1). This process requires significant communication
between MPI nodes and, depending on the particular case, can have a
significant impact on the total computing time. For this reason, if no line-
of-sight information is requested from CoLoRe (i.e. if only sources or
intensity maps are requested, with no associated lensing information or
skewers), steps 7 and 8 are skipped, often leading to a significant speed-up.

Furthermore, computing the lensing observables along each line-of-sight for
large catalogs containing billions of sources, as would be the case e.g. if
simulating the full LSST shear sample, is a computationally demanding task
that can dominate by far all other steps in a CoLoRe run. To avoid this, an
alternative approximate scheme to compute these quantities can be used. In the
“fast lensing” scheme, the five lensing observables (displacement, convergence
and shear) are precomputed on a set of spherical HEALPix maps at constant
radial comoving distance intervals, using the same method outlined in Section
2.4.3. The angular resolution parameter $N_{\rm side}$ of each map is
adaptively chosen so that the physical pixel size is smaller than the
Cartesian cell size, in order to avoid oversampling as well as degrading the
original three-dimensional resolution. The hierarchical nature of the HEALPix
scheme makes it possible to relate a pixel in a given map with pixels in all
lower-resolution maps with smaller radii. The lensing observables for a given
source are then calculated by interpolating between the values of that
observable along the line of pixels corresponding to the source’s angular
coordinates.

Figure 4 shows maps of the quantities described in this section for a small
source catalog simulated with CoLoRe. The simulated sample had a redshift
distribution

$\frac{dN}{dz}\propto\left(\frac{z}{z_{0}}\right)^{2}\exp\left[-\left(\frac{z}{z_{0}}\right)^{3/2}\right],$
(2.30)

with $z_{0}=0.07$, extending up to $z\lesssim 0.3$. The simulation was run
using the first-order LPT structure formation model. The different panels show
the source overdensity and mean redshift distortion (top panels), the mean
lensing shear (middle panels) and the mean lensing deflection vector (bottom
panels). Figure 5 shows density and velocity skewers (top and bottom panels
respectively) for two arbitrary sources in the same catalog.

####  2.4.5 Line intensity mapping

Figure 6: Slice through one of the 21cm intensity maps generated with CoLoRe.

Consider a species of gas emitting at a rest-frame frequency $\nu_{0}$ due to
some atomic or molecular transition. The intensity (flux per unit frequency)
measured in a patch around $\hat{\bf n}$ with solid angle $\delta\Omega$, and
in a frequency interval $\delta\nu$ around $\nu$ is given by [39]

$I(\nu,\hat{\bf n})=\frac{\hbar A_{21}\nu_{0}x_{2}}{2m_{a}}\frac{M_{\rm
em}}{(1+z)^{2}\chi^{2}\delta\Omega\,\delta\nu},$ (2.31)

where $A_{21}$ is the Einstein coefficient for the transition, $m_{a}$ is the
atomic mass of the emitting gas, $x_{2}$ is the fraction of the gas in the
excited state, $M_{\rm em}$ is the total mass of the emitting gas in the voxel
defined by $(\delta\nu,\delta\Omega)$. Associating this intensity with a
black-body temperature in the Rayleigh-Jeans regime
($T=c^{2}I/(2k_{B}\nu^{2})$), this can be rewritten as:

$T(\nu,\hat{\bf n})=\bar{T}(z)\left[1+\delta_{\rm em}(z,\hat{\bf n})\right],$
(2.32)

where $\delta_{\rm em}$ is the overdensity of the emitting gas in redshift
space (i.e. accounting for the effects of redshift-space distortion),
$z=\nu_{0}/\nu-1$, and the mean temperature $\bar{T}$ is

$\bar{T}(z)\equiv\frac{3\,\hbar\,A_{21}\,x_{2}\,x_{\rm
em}(z)\,\Omega_{b,0}\,H_{0}^{2}\,c^{2}\,(1+z)^{2}}{32\pi\,G\,k_{B}\,m_{a}\,\nu_{0}^{2}\,H(z)}.$
(2.33)

CoLoRe generates mock intensity mapping observations by modelling $\delta_{\rm
em}$ as a biased tracer of $\delta_{M}$ using one of the bias relations
described in Section 2.4.1. The code estimates the mean brightness temperature
in each Cartesian voxel in terms of the value of the matter overdensity, and a
redshift-dependent mean temperature and bias. In order to interpolate from the
Cartesian grid on to temperature maps, while accounting for the effects of
redshift space distortions, each cell is first sub-divided into smaller sub-
cells. The Cartesian coordinates of each sub-cell are translated into angular
coordinates $(\theta,\varphi)$ and emission frequency
$\nu=\nu_{0}/(1+z+v_{r})$, where $z$ is the redshift corresponding to the sub-
cell’s radial comoving distance, and $v_{r}$ is the value of the radial
comoving peculiar velocity field. Each sub-cell is then assigned to a given
pixel and frequency band based on $(\theta,\varphi,\nu)$. The final intensity
maps are then generated by averaging over the temperature of all sub-cells
thus assigned to each pixel.

The implementation used by CoLoRe is very similar to that used by CRIME [25].
CoLoRe does not emulate any other secondary effects, such as second-order
gravitational lensing [40, 41] or self-absorption. If desired, gravitational
lensing could be simulated by perturbing the angular positions of the
Cartesian sub-cells above with the deflection field described in Section
2.4.3. Figure 6 shows a slice through a set of intensity maps corresponding to
the 21cm line simulated using the model described in [25] in the range
$\nu\in[1015,646]\,\,{\rm MHz}$ ($0.4<z<1.2$).

###  2.5 Lognormal predictions

Although, as discussed in Section 2.3.1, the log-normal transformation is not
able to recover the right properties of the non-linear density fluctuations at
all orders, one of its key advantages is the possibility to produce exact
analytical predictions for the two-point correlators of a lognormal field.
This makes it possible to compare the output of a CoLoRe simulation against
precise predictions, making it straightforward to, for example, quantify the
impact of different observational systematic effects added in post-processing
on the main observables of a large-scale structure experiment.

The real-space two-point correlation function of a lognormal field is related
to that of its parent Gaussian field via:

$1+\xi_{\rm LN}(r)=\exp\left[\xi_{\rm G}(r)\right],$ (2.34)

where $\xi_{\rm LN}(r)\equiv\langle\delta_{\rm LN}({\bf x})\delta_{\rm
LN}({\bf x}+{\bf r})\rangle$ and $\xi_{\rm G}(r)\equiv\langle\delta_{\rm
G}({\bf x})\delta_{\rm G}({\bf x}+{\bf r})\rangle$. Its cross-correlation with
the Gaussian field itself is simply

$\langle\delta_{\rm G}({\bf x})\delta_{\rm LN}({\bf x}+{\bf
r})\rangle=\xi_{\rm G}(r).$ (2.35)

Computing the power spectrum of the lognormal field $P_{\rm LN}(k)$ from that
of the parent Gaussian field $P_{\rm G}(k)$ is thus a simple three-step
process:

  1. 1.

Compute the Gaussian correlation function from $P_{\rm G}(k)$ via

$\xi_{\rm G}(r)=\frac{1}{2\pi^{2}}\int_{0}^{\infty}dk\,k^{2}\,P_{\rm
G}(k)\,\frac{\sin(kr)}{kr}.$ (2.36)

  2. 2.

Compute $\xi_{\rm LN}$ from $\xi_{\rm G}$ through Eq. 2.34.

  3. 3.

Compute the $P_{\rm LN}(k)$ from $\xi_{\rm LN}(r)$ via

$P_{\rm LN}(k)=4\pi\int_{0}^{\infty}dr\,r^{2}\,\xi_{\rm
LN}(r)\,\frac{\sin(kr)}{kr}.$ (2.37)

CoLoRe automatically produces and outputs predictions for the power spectrum
and correlation function of any lognormal biased tracers simulated. These
predictions are valid as long as the simulation was run using the lognormal
structure formation model, and the exponential bias model (see Section 2.4.1),
which retains the lognormal nature of the biased fields.

It is worth noting that, in order to produce truly accurate predictions, it is
necessary to account for all sources of smoothing produced by the different
operations carried out by CoLoRe:

  * •

The chosen Gaussian smoothing scale.

  * •

The finite resolution of the Cartesian grid.

  * •

The effects of interpolating from the Cartesian grid onto beam-related
quantities (line-of-sight integrals, skewers etc.).

The latter two effects (finite grid resolution and interpolation), have an
associated Fourier-space window function given by

$W_{n}({\bf k})=\left[\left(2\frac{\sin(2k_{x}\Delta x)}{k_{x}\Delta
x}\right)\left(2\frac{\sin(2k_{y}\Delta x)}{k_{y}\Delta
x}\right)\left(2\frac{\sin(2k_{z}\Delta x)}{k_{z}\Delta x}\right)\right]^{n},$
(2.38)

where $\Delta x$ is the grid spacing, and $n=1$ and 2 for nearest-neighbour
and trilinear interpolation, respectively. For $k\ll 1/\Delta x$, we can
approximate $W_{n}$ as a Gaussian filter via

$W_{n}({\bf k})\simeq\left[1-\frac{(k_{x}^{2}+k_{y}^{2}+k_{z}^{2})(\Delta
x)^{2}}{24}\right]^{n}\simeq e^{-(kR_{G})^{2}/2}.$ (2.39)

where $R_{G}^{2}\equiv n(\Delta x)^{2}/12$. Thus, the effective smoothing
scale associated to the appropriate number of interpolation operations can be
simply added in quadrature to the chosen Gaussian smoothing scale to produce
theoretical predictions.

##  3 Results

###  3.1 Validation

Figure 7: Redshift distribution of the two tomographic bins used in the
analysis of the validation simulations (red and blue lines), as well as the
overall redshift distribution (black dashed line). The bins are defined by a
cut in photometric redshift space at $z_{\rm photo}=0.7$, where we assigned
each source a random photometric redshift error with standard deviation
$\sigma_{z}=0.03\,(1+z)$. Figure 8: Galaxy clustering, cosmic shear and
lensing displacement power spectra. The gray bands show the 68% scatter from
the 100 validation realizations, and the theoretical predictions, described in
Section 2.5, are shown as black solid lines. The lower-left corner shows all
auto- and cross-correlations between the two galaxy clustering and cosmic
shear bins. The upper-right corner shows the auto- and cross-correlations
between the lensing displacement vectors in both redshift bins. Figure 9:
Cross-correlations between ISW (top row) and convergence maps (bottom row) at
$z=1$ and the two high-redshift clustering and shear samples in the validation
simulations. The gray bands show the 68% scatter from the 100 validation
realizations, and the theoretical predictions, described in Section 2.5, are
shown as black solid lines. Figure 10: Shear power spectra for the first
redshift bin in the validation simulations. The blue and red points show the
results (mean and standard deviation of the 100 simulations) with no
correction for the effects of source clustering. The gray points show the
result of applying this correction by simply subtracting the $B$-mode power
spectrum from the $E$-mode one. The black line shows the theoretical
prediction for the latter. The orange band shows the 1$\sigma$ uncertainties
one would find in the presence of realistic shape noise, which would make the
source clustering effect undetectable in practice.

CoLoRe has been the basis of other analyses, and some of its functionality was
validated in previous work. The ability to produce large-scale intensity maps
for the 21cm neutral hydrogen line was presented, used, and validated in [25,
42], and [43] used the code to explore cross-correlations with galaxy
clustering data. [26] used CoLoRe’s line-of-sight skewers to produce mock
observations of the Lyman-$\alpha$ forest; these mocks were extensively used
to validate the final Lyman-$\alpha$ BAO analysis of the eBOSS collaboration
[44]. Our discussion here therefore focuses on presenting and validating
CoLoRe as a tool to produce fast simulations for wide galaxy surveys targeting
galaxy clustering (both spectroscopic and photometric), weak lensing shear,
and their cross-correlation with maps of the lensing convergence, and the ISW
effect.

To do so, we have run a set of 100 large CoLoRe realisations containing these
observables, and compared the relevant two-point correlations from different
pairs of tracers in the simulation with the corresponding theoretical
predictions.

####  3.1.1 Simulations

We generate 100 realisations covering the volume up to redshift $z=1.4$,
corresponding to a box size $L_{\rm box}=5{,}843\,{\rm Mpc}/h$. We use a grid
of size $N_{\rm grid}=2048$, with cell size $\Delta x=2.85\,{\rm Mpc}/h$. Each
realisation is populated with a single galaxy sample with the redshift
distribution shown in Fig. 7, and a total number density
$\bar{n}_{g}=2.9\,{\rm arcmin}^{-2}$. We used a lognormal structure formation
model, and an exponential bias model with bias $b(z)=1+0.65z+0.03z^{3}$,
compatible with a blue galaxy sample [45]. We also used CoLoRe to generate
maps of the lensing convergence and the ISW effect at $z=1$. We used
cosmological parameters $\Omega_{m}=0.3$, $\Omega_{b}=0.05$, $h=0.7$,
$n_{s}=0.96$, $\sigma_{8}=0.8$. The Gaussian overdensity field was smoothed
with a Gaussian kernel with width $R_{G}=2\,{\rm Mpc}/h$.

In order to simulate the effects of redshift uncertainties in photometric
redshift surveys, we assigned a random Gaussian error to each galaxy redshift
with standard deviation $\sigma_{z}=0.03\,(1+z)$. We then divided all galaxies
into two redshift bins, corresponding to sources above and below redshift
$z_{\rm thr}=0.7$. The redshift distributions of the resulting bins are shown
in Fig. 7. For each redshift bin, we created maps of the projected galaxy
overdensity, as well as the shear and lensing displacement vectors. We did not
make use of the fast lensing scheme described in Section 2.4.4 for these
simulations. The latter quantities are provided by CoLoRe at each source.
Finally, we computed all auto- and cross-power spectra between these maps, as
well as the corresponding correlations with the lensing convergence and ISW
maps at $z=1$ produced by CoLoRe. In order to optimally weight the shear and
displacement fields by the local number of sources when computing these power
spectra (and in order to account for their spin-2 and spin-1 nature), we make
use of NaMaster [46].

For the 3D clustering validation we used 10 out of the 100 realisations. Since
the aim in this case is simulating a spectroscopic survey, we do not add
photometric redshift uncertainties. We used the corrfunc [47] package to
obtain measurements of the monopole and the quadrupole correlation functions
for two different redshfit bins $(0.5,0.7)$ and $(0.7,0.9)$; probing
separations between $r=0.1\ {\rm Mpc}/h$ and $r=200\ {\rm Mpc}/h$ in 41
linearly spaced bins.

####  3.1.2 Validation of 2D observables

Before presenting the results from this validation exercise, it is worth
clarifying a technical point about the theoretical predictions used to compare
with the two-point functions estimated from the simulation. As described in
Section 2.5, besides the effects of the lognormal transformation, one must
account for the additional smoothing associated with the finite grid and the
different interpolation operations. We do so by adding an extra smoothing
scale in quadrature to the Gaussian smoothing scale used in the simulation,
with the form

$\Delta R_{G}^{2}=n_{\rm eff}\,\frac{(\Delta x)^{2}}{12},$ (3.1)

where the prefactor $n_{\rm eff}$ depends on the finite-resolution effects
that must be taken into account. For instance, a single nearest-neighbour
interpolation would correspond to $n_{\rm eff}\simeq 1$, while linear
interpolation would have $n_{\rm eff}\simeq 2$. Thus, since galaxies are
assigned to their nearest grid cell, and since two linear interpolations are
needed to assign lensing properties to sources, the galaxy clustering and
cosmic shear tracers produced by CoLoRe should take additional smoothing with
$n_{\rm eff}\simeq 1$ and $4$ respectively (not taking into account the
additional interpolations involved in the fast lensing scheme). Since this
additional Gaussian smoothing is not exact, we must adapt these prefactors
slightly, in order to improve the agreement between theory predictions and
simulations on small scales (high $\ell$s). For instance, we find that the
galaxy-galaxy, galaxy-matter, and matter-matter power spectra must be smoothed
by additional factors $n_{\rm eff}^{gg}=0.9$, $n_{\rm eff}^{gm}=4$, and
$n_{\rm eff}^{mm}=3.8$ to recover the galaxy clustering and cosmic shear power
spectra from the CoLoRe simulations. Therefore, care must be exercised when
interpreting the results of CoLoRe simulations, particularly on physical
scales smaller than, or comparable with, the grid resolution. The final users
are encouraged to tune these $n_{\rm eff}$ factors to their particular
analysis, depending on their use case and required level of accuracy.

The gray bands in Figure 8 show the 1$\sigma$ scatter of the power spectra
estimated from the 100 validation simulations, together with the corresponding
theoretical prediction in solid black. All auto- and cross-correlations
between the projected galaxy overdensity and the cosmic shear field in the two
redshift bins described above are shown in the lower-left panels in the
figure. The upper-right panels show the three auto- and cross-correlations
between the $E$-mode fields corresponding to the lensing displacement vectors
in both redshift bins. Similar results are shown in Figure 9 for correlations
involving the lensing convergence and ISW maps, and the projected overdensity
and cosmic shear fields in the second. A good agreement between theory and
simulations, at the $1\sigma$ level, is found in all cases.

One further complication must be noted. Cosmic shear is measured at the
positions of clustered sources, and this clustering is correlated with the
cosmic shear signal itself. This induces an additional contribution to the
observed cosmic shear statistics which gives rise to shear $B$-modes, in
addition to modifying the $E$-mode power spectrum. This effect is well known,
and should be unobservable in most cases [48]. However, in the absence of
shape noise (i.e. since we have access to the true weak lensing shear at each
source), the effect can be measured in the catalogs produced by CoLoRe. This
is illustrated in Fig. 10, which shows the $E$ and $B$-mode power spectra for
the first redshift bin in the validation simulations without any correction
for source clustering (blue and red respectively) in comparison with the
theory prediction without source clustering (black line) and the expected
uncertainties in the presence of realistic shape noise (orange band). We find
that, in practice, the effects of source clustering in the $E$-mode power
spectrum can be corrected by simply subtracting the $B$-mode power spectrum
from the $E$-mode power spectrum (gray points in the figure).

####  3.1.3 Validation of 3D clustering

We validated the 3D clustering of simulated galaxies both in real and in
redshift space. The top panel of Figure 11 shows the real-space monopole for
two spectroscopic redshift bins: (0.5-0.7) and (0.7-0.9). The measurement of
the correlation function was done separately in 48 different healpix pixels
defining the CoLoRe beams, and the error bands (estimated from the scatter of
these measurements) show the uncertainty for a single realization. Solid lines
show the prediction derived in Section 2.5 with the linear bias left as a free
parameter. The best-fit value of bias (using separations larger than 10 Mpc/h)
agrees with the input value at the 2% level.

Figure 11: Measurements of the correlation function from the stack of 10
realizations used to validate the 3D clustering. The low-z samples take
redshifts from 0.5 to 0.7, while high-z samples take redshifts from 0.7 to
0.9. The lines show the model, solid in the regions where the bias was fitted,
and the shaded bands show the error for a single realization. _Top_ :
Measurements of the monopole in real space. _Bottom_ : Measurements of the
monopole and quadrupole in redshift space.

The bottom panel in Figure 11 shows the redshift-space monopole and quadrupole
for the same redshift bins and linear bias. We do not have a prediction for
the clustering of galaxies in redshift space that is valid on all scales, but
we use a modified version of the lognormal model that includes linear
redshift-space distortions (Kaiser model [49]):

$\delta_{\rm LN}^{s}({\bf k})=\delta_{\rm LN}({\bf k})+f\mu^{2}\delta_{G}({\bf
k})$ (3.2)

where the Gaussian term stands due to the fact that velocities comes directly
from the gravitational potential (Eq. 2.29). The redshift-space power spectrum
is then:

$P_{\rm LN}^{s}(k,\mu)=P_{\rm LN}(k)+f^{2}\mu^{4}P_{\rm G}(k)+2bf\mu^{2}P_{\rm
G}(k)$ (3.3)

where we used Eq. 2.35 when computing the last term. Similarly to the 2D
clustering, we added a smoothing associated with the finite grid with $n_{\rm
eff}=1$. We also added an extra smoothing to correct for the binning of the
correlation function measurement.

Differences on small scales between the measured and the predicted quadrupoles
are due to inaccuracies in our RSD modelling, in particular higher-order terms
ignored in Eq. 3.2. We discuss the role of these higher-order terms in
Appendix A, and show that they are indeed the cause of this disagreement.

###  3.2 Performance at scale

In the context of existing and next-generation cosmological experiments,
CoLoRe should be able to generate mock observations covering large volumes
($z\lesssim 3$), with reasonable resolution ($\Delta x={\cal O}(1)\,{\rm
Mpc}$) for a wide range of observables (galaxy positions, shear, CMB
lensing555Note that CoLoRe can provide the contribution to the CMB lensing
convergence up to the highest redshift covered by the simulation box.
Contributions from higher redshifts can then be added as a correlated Gaussian
random field., intensity maps, etc.). This section quantifies the feasibility
of these simulations.

#### Large-volume multi-tracer simulations

Figure 12: Redshift distribution of the different tracers simulated in our two
flagship simulations. These include the LSST gold sample (black), the 4 DESI
samples (blue) and 4 different radio continuum samples observable by SKA
(red).

As an example of the type of mocks needed for Stage-IV surveys, we have used
CoLoRe to generate two simulations containing cosmological probes for three
major experiments:

  * •

DESI. We simulate the Bright Galaxy Survey (BGS), Large Red Galaxies (LRG),
Emission Line Galaxies (ELG), and quasi-stellar object (QSO) targets
(including Ly-$\alpha$ skewers for the QSO sample) using the nominal $N(z)$,
and bias functions from the DESI white paper [1]. The resulting catalog
contains $\sim 3.4\times 10^{7}$ BGSs, $\sim 1.5\times 10^{7}$ LRGs, $\sim
10^{8}$ ELGs, and $\sim 6.5\times 10^{6}$ QSOs over the full celestial sphere.

  * •

LSST. We simulate a sample similar to the LSST "Gold" sample ($i\lesssim
25.3$, $\bar{n}\sim 40$ galaxies/arcmin2), resulting in $\sim 6$ billion
sources across the full sky. We follow the redshift distribution of the DESC
Science Requirements Document [50], and assume a linear bias with redshift
dependence $b(z)=0.95/D(z)$ [50, 51]. The lensing shear, convergence and
displacement is calculated for all sources using the fast lensing scheme
described in Section 2.4.4.

  * •

SKA. We simulate a radio continuum catalog comprised of 4 different radio
galaxy samples: FRI radio-loud AGNs, radio-quiet AGNs (RQQs), normal star-
forming galaxies (SFGs), and starbursts (SBs). For this we follow the models
for the redshift distributions and linear bias described in [52]. The
resulting samples contain $\sim 4.5\times 10^{7}$ FRIs, $\sim 1.3\times
10^{8}$ RQQs, $\sim 8\times 10^{8}$ SFGs, and $\sim 8\times 10^{7}$ SBs over
the full sky. In addition to this, we generate simulated HI intensity mapping
observations for 490 frequency bands covering the range $473\,{\rm
MHz}<\nu<947\,{\rm MHz}$ (corresponding to redshifts $0.5<z<2$). The maps were
generated with an angular resolution $N_{\rm side}=256$, corresponding to
pixels about 4 times smaller than the SKA primary beam in single-dish mode at
the highest frequency ($\theta_{\rm FWHM}\sim 1.2^{\circ}$).

In addition to these, to showcase the ability of CoLoRe to generate CMB
lensing observations, we generate a convergence map at resolution $N_{\rm
side}=1024$ at $z=3$, caused by the same matter density field that serves as
seed for the tracers listed above. The redshift distributions of the different
galaxy samples simulated are shown in Fig. 12. For illustrative purposes, Fig.
13 shows the 3-dimensional representation of some of the quantities simulated
by CoLoRe in one of the beams used by the code.

Figure 13: Visual description of the multi-tracer products that can be
simulated with CoLoRe. The upper plot shows one of the beams used internally
by CoLoRe for domain decomposition, with the redshift and angular coordinates
of 3 different DESI samples, and maps of the density, radial velocity, and
lensing shear constructed from sources at two different redshifts. The lower
plot shows the density skewers calculated for three arbitrary DESI quasars
contained in the same beam, as a function of comoving distance.

The simulations were generated with a $\Lambda$CDM model compatible with the
best-fit Planck cosmological parameters [53]. The boxes span the redshift
range $(0<z<3)$, with $N_{\rm grid}=4096$, resulting in a spatial resolution
of $\approx 2$ Mpc/$h$.

Figure 14: Fraction of the total run time taken up by different stages of a typical CoLoRe simulation. The stages shown are, the generation of the initial Gaussian random fields in Fourier space ("GRF"), their transformation to real space ("FFT"), the structure formation model leading to a positive-definite matter overdensity ("LN" or "1LPT" for the lognormal and first-order LPT simulations), the generation of source catalogs via Poisson sampling ("Poisson"), the redistribution of these sources across different nodes before any line-of-sight calculations ("Source redist."), the calculation of all relevant line-of-sight quantities ("LOS tracing"), and the output of all final products to disk ("Write"). Note that the LOS-tracing stage is more sensitive to inter-node communication than other stages, and thus it takes up the same fraction of the total compute time in both simulations in spite of the additional time taken by the 1LPT stage, given the larger number of MPI nodes needed for that simulation. Tracer | $\Delta z_{\rm RSD}$ | $\bm{\alpha}$ | $\bm{\gamma}$ | $\kappa$ | $\delta(z),v_{r}(z)$ | Memory (GB) | Disk (GB)
---|---|---|---|---|---|---|---
LSST | ✓ | ✓ | ✓ | ✓ | ✗ | 308 | 155
DESI-BGS | ✓ | ✗ | ✗ | ✗ | ✗ | 1.7 | 0.6
DESI-LRG | ✓ | ✗ | ✗ | ✗ | ✗ | 0.8 | 0.3
DESI-ELG | ✓ | ✗ | ✗ | ✗ | ✗ | 5.2 | 1.9
DESI-QSO | ✓ | ✗ | ✗ | ✗ | ✓ | 99 | 99
SKA-FRI | ✓ | ✗ | ✗ | ✗ | ✗ | 2.3 | 0.8
SKA-RQQ | ✓ | ✗ | ✗ | ✗ | ✗ | 6.8 | 2.5
SKA-SFG | ✓ | ✗ | ✗ | ✗ | ✗ | 39 | 15
SKA-SB | ✓ | ✗ | ✗ | ✗ | ✗ | 4.2 | 1.5
SKA-21cm | ✓ | N.A. | N.A. | N.A. | N.A. | 2.8 | 1.4
Convergence map | N.A. | N.A. | N.A. | ✓ | N.A. | 0.1 | 0.05
$\delta$, $\phi_{N}$ grids | N.A. | N.A. | N.A. | N.A. | N.A. | 768 | N.A.
$\bm{\Psi}_{\rm 1LPT}$ grids | N.A. | N.A. | N.A. | N.A. | N.A. | 1229 | N.A.
Table 1: Simulation products generated by CoLoRe. The first 11 rows show the
different tracers generated for the large-volume simulations described in the
text. For each tracer we show the physical quantities simulated (RSDs, lensing
displacements, shear, convergence, and density/velocity skewers), as well as
the memory and disk space taken. The last two rows display the memory
requirements for the different Cartesian grids stored for lognormal and 1LPT
simulations. Although the memory requirements are dominated by these cartesian
grids (particularly for LPT simulations), the simulated tracers can take up a
non-negligible fraction of the available memory. This is patent for the LSST
sample, given its size, and the need to store lensing information, and for the
DESI quasar sample, since we save a full density/velocity skewer for each
source.

#### Run time, memory usage, and fast lensing

Both simulations were initialised with the same random seed, but using
different structure formation models, lognormal (LN) and first-order LPT
(1LPT) respectively. Both simulations were run at NERSC 666https://nersc.gov.
The LN simulation was generated using 40 MPI tasks, distributed across 20
Cori-haswell nodes, using 16 OMP threads per task. This simulation ran in 1.15
hours, using approximately 730 CPU-hours. The 1LPT simulation required a
larger number of nodes, given the additional memory needed to allocate the
three more Cartesian grids mentioned in Section 2.3.2. This simulation was run
using 72 MPI tasks distributed across 36 Cori-haswell nodes, using 16 OMP
threads per task. The simulation ran in 0.93 hours using a total of 1,075 CPU-
hours. A suite of 1000 such simulations could therefore be run using $\sim 1$
million CPU-hours.

Table 1 lists the memory and disk requirements associated with each of these
tracers777Note that, naively, we have simulated all galaxy tracers as disjoint
samples when, in reality, there would be significant overlap between them
(e.g. between SKA SFGs and the LSST sample). A more realistic setting should
therefore account for these overlaps when generating the different galaxy
catalogs.. Although the most memory-demanding task is the generation of the 3D
density and Newtonian potential fields, or the Lagrangian displacement
components if using LPT, the final data products also lead to a non-negligible
memory requirement, particularly in the case of high-density samples such as
LSST, or for a large number of density/velocity skewers ($\sim 27\%$ and $\sim
4\%$ respectively in the case of the lognormal simulation).

Figure 14 shows the time taken by the different stages in both runs. In both
cases the slowest stage is the collection of line-of-sight information, such
as the density and velocity skewers stored for all DESI QSOs, or the lensing
information associated with the LSST sources. Given the large number of
sources in the LSST sample, this stage would completely dominate the run time
if we had not used the “fast lensing” method described in Section 2.4.4.

To quantify this, we ran two additional identical simulations with a total of
$5.9$ billion galaxies, emulating a blue galaxy population of LSST 10-year
depth. Both simulations were generated using 32 MPI tasks distributed across
16 Cori Haswell nodes at NERSC. Each MPI task had 16 OMP threads. The
simulation using the fast lensing scheme took $\sim 35$ minutes to run (a
total of $173$ CPU-hours), while the one integrating along each galaxy’s line
of sight took $\sim 36.3$ hours (i.e. $\sim 18,600$ CPU-hours), more than 60
times longer. The additional time was completely taken by the line of sight
calculations. The fast-lensing scheme thus allows for a factor $\sim 60$
speed-up and, in the current implementation of CoLoRe, is absolutely necessary
for simulations with billions of sources.

Figure 15: Relative difference between the shear power spectrum of a CoLoRe
simulation run using the "fast lensing" scheme, and a simulation run without
this approximation. The effects caused by the various interpolations carried
out as part of the fast scheme should be carefully modelled if an accurate
theoretical description of the CoLoRe output is required.

The price to pay for this speed-up is the additional complexity associated
with the interpolation operations involved in the fast scheme (interpolation
from the cartesian grid to the fixed pixel lines-of-sight, and from those to
the galaxy positions). Fig. 15 shows the relative difference between the
cosmic shear power spectrum computed from a simulation run with the fast
lensing scheme, and one run without this approximation. The effects of the
fast lensing approximation are a decrease in the final lensing amplitude, and
a gradual loss of power at higher multipoles, on scales comparable with the
size of the adaptive pixels used in this approximation. These effects would
need to be accurately characterized if the application of the CoLoRe
realizations requires an accurate theoretical prediction of two-point
statistics involving weak lensing observables. The accuracy/speed-up trade-off
can be controlled by the user by changing the number of shells in which to
compute the lensing information, as well as the size of the adaptive pixels in
relation with the Cartesian cell size.

##  4 Conclusions

In this paper we have introduced CoLoRe, a public code to efficiently generate
synthetic realisations of multiple cosmological surveys. We started in Section
2 by describing the overall structure of the code, and the different methods
to simulate the density field. We have presented the available tracers in
CoLoRe, and how to add new ones using its highly modular structure. We
concluded this section discussing the accurate predictions available when
working with the lognormal model of structure formation.

In Section 3 we presented the validation results of some of the key summary
statistics from the simulated maps, using a large set of CoLoRe boxes. We
showed that the measured angular power spectra from simulated photometric
surveys agree with the theoretical predictions up to $\ell=1000$. This is true
for galaxy correlations, galaxy-shear cross-correlations and several lensing
statistics (shear, convergence, displacements). The 3D clustering in simulated
spectroscopic galaxy surveys was also validated in the absence of redshift
space distortions, where the lognormal model can be accurately predicted on
all scales. As discussed in Appendix A we do not have a good model for the
small-scales multipoles in redshift space, but the agreement is very good on
large, linear scales.

We have discussed the performance of CoLoRe at scale by presenting two joint
simulations of DESI, LSST and SKA. These large boxes cover the whole comoving
volume out to $z=3$, with a resolution of $\approx 2\rm{Mpc}/h$, and include
spectroscopic and photometric galaxies, lensing, intensity mapping and radio
galaxies. The more realistic simulation, using Lagrangian perturbation theory,
only used about 1,000 CPU-hours. Simulating hundreds of these boxes is
entirely feasible, and can be used to characterise systematic effects in
multi-experiment analyses, or estimate cross-survey covariances.

Finally, we discussed the differences between the two options used in CoLoRe
to compute weak lensing variables from source galaxies. The fast lensing
implementation provides a factor of $\sim$60 speed-up for an LSST-like sample,
while maintaining an accuracy better than a 2-3% bias in the amplitude of the
shear power spectrum on large scales ($\ell<1000$).

There are several features that could be added to CoLoRe without major changes
in the code structure:

  * •

Better structure formation: Currently CoLoRe can simulate the growth of
structure using a lognormal model or Lagrangian Perturbation Theory (LPT)
computed at first or second order. Future versions of the code could add new
modules to use more complex models of the growth of structure, such as
COmoving Lagrangian Acceleration (COLA, [18] as in [54], where it was used to
generate weak lensing maps and halo catalogues in the lightcone) or Fast
Particle-Mesh (FastPM, [55]) algorithms. The potential additional compute time
and memory requirements associated with these methods should be weighed
against the need for more accurate clustering statistics of a given
application.

  * •

Non-linear RSD: Minor modifications of the code could improve the level of
realism of the RSDs, by sourcing the velocities from the computed LPT fields.
This change could have a significant impact on the 3D clustering of galaxies
intermediate scales. However, CoLoRe does not currently simulate virialized
objects, and therefore we are not able to properly capture non-linear peculiar
velocities or Fingers of God. Random virial motions (or redshift errors) can
be added in post-processing, by assigning random shifts to the galaxy
redshifts, but a more realistic approach of non-linear RSD is beyond the scope
of this paper.

  * •

Halos: The models currently used by CoLoRe to generate different tracer
observations directly connect the latter with the underlying smooth density
field. The complexity and fidelity of these simulations could be improved if
this density field was endowed with a halo catalog. This could be done
directly at the level of the linear density field using Press-Schechter-
inspired methods (as in e.g. [56]), or from the LPT displacement field using a
modified friends-of-friends search [15], or peak-patch methods [57, 58]. A
halo catalog would allow us to improve the fidelity of the resulting density
field on small scales by including the expected density profile of these
halos, and would make it possible to use the halo model to generate simulated
observations of additional tracers (e.g. thermal or kinematic Sunyaev-
Zel’dovich effects, halo-occupation distributions for better galaxy catalogs
etc.).

  * •

Small areas and flat-skies: The current version of CoLoRe simulates the whole
Universe to a given redshift, with the observer placed in the center of a
large box. Users interested in simulating surveys with relatively small areas
might prefer to place the observer in one side of rectangular box and simulate
only a (literal) light-cone. For sufficiently small sky areas, these
simulations could be made faster making use of the flat-sky approximation.
Both of these features should be easy to implement in CoLoRe.

The addition of these features should not significantly impact the performance
of CoLoRe reported here. By making our code public we want to encourage the
different collaborations preparing the next generation of large cosmological
surveys to use CoLoRe to efficiently generate realistic synthetic simulations
of their datasets to be used in multi-survey analyses.

## Acknowledgments

We would like to thank James Farr, David Kirkby, Stephane Plaszczynski, and
Anẑe Slosar for contributions to the CoLoRe repository and for useful
discussions in early stages of this project. We would also like to thanks
Hanyu Zhang for his assistance with the measurement of the 3D correlations
using corrfunc. CRP is partially supported by the Spanish Ministry of Science
and Innovation (MICINN) under grants PGC-2018-094773-B-C31 and SEV-2016-0588.
DA is supported by the Science and Technology Facilities Council through an
Ernest Rutherford Fellowship, grant reference ST/P004474. AFR is supported by
MICINN with a Ramón y Cajal contract (RYC-2018-025210). IFAE is partially
funded by the CERCA program of the Generalitat de Catalunya.

We made extensive use of computational resources at the University of Oxford
Department of Physics, funded by the John Fell Oxford University Press
Research Fund. We also used resources of the National Energy Research
Scientific Computing Center (NERSC), a U.S. Department of Energy Office of
Science User Facility located at Lawrence Berkeley National Laboratory,
operated under Contract No. DE-AC02-05CH11231. This manuscript has been
authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359
with the U.S. Department of Energy, Office of Science, Office of High Energy
Physics.

We made extensive use of the numpy [59, 60], scipy [61], astropy [62, 63],
healpy [64], NaMaster [46], corrfunc [47], matplotlib [65] and lmfit [66]
python packages.

## References

  * [1] DESI Collaboration, A. Aghamousa, J. Aguilar, S. Ahlen, S. Alam, L. E. Allen et al., _The DESI Experiment Part I: Science,Targeting, and Survey Design_ , _arXiv e-prints_ (2016) arXiv:1611.00036 [1611.00036].
  * [2] R. Laureijs, J. Amiaux, S. Arduini, J. L. Auguères, J. Brinchmann, R. Cole et al., _Euclid Definition Study Report_ , _arXiv e-prints_ (2011) arXiv:1110.3193 [1110.3193].
  * [3] D. Spergel, N. Gehrels, C. Baltay, D. Bennett, J. Breckinridge, M. Donahue et al., _Wide-Field InfrarRed Survey Telescope-Astrophysics Focused Telescope Assets WFIRST-AFTA 2015 Report_ , _arXiv e-prints_ (2015) arXiv:1503.03757 [1503.03757].
  * [4] LSST Science Collaboration, P. A. Abell, J. Allison, S. F. Anderson, J. R. Andrew, J. R. P. Angel et al., _LSST Science Book, Version 2.0_ , _arXiv e-prints_ (2009) arXiv:0912.0201 [0912.0201].
  * [5] P. Ade, J. Aguirre, Z. Ahmed, S. Aiola, A. Ali, D. Alonso et al., _The Simons Observatory: science goals and forecasts_ , _J. Cosmology Astropart. Phys_ 2019 (2019) 056 [1808.07445].
  * [6] K. N. Abazajian, P. Adshead, Z. Ahmed, S. W. Allen, D. Alonso, K. S. Arnold et al., _CMB-S4 Science Book, First Edition_ , _arXiv e-prints_ (2016) arXiv:1610.02743 [1610.02743].
  * [7] Square Kilometre Array Cosmology Science Working Group, D. J. Bacon, R. A. Battye, P. Bull, S. Camera, P. G. Ferreira et al., _Cosmology with Phase 1 of the Square Kilometre Array Red Book 2018: Technical specifications and performance forecasts_ , _PASA_ 37 (2020) e007 [1811.02743].
  * [8] L. B. Newburgh, K. Bandura, M. A. Bucher, T. C. Chang, H. C. Chiang, J. F. Cliche et al., _HIRAX: a probe of dark energy and radio transients_ , in _Ground-based and Airborne Telescopes VI_ , H. J. Hall, R. Gilmozzi and H. K. Marshall, eds., vol. 9906 of _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , p. 99065X, Aug., 2016, 1607.02059, DOI.
  * [9] P. Coles and B. Jones, _A lognormal model for the cosmological mass distribution._ , _MNRAS_ 248 (1991) 1.
  * [10] H. S. Xavier, F. B. Abdalla and B. Joachimi, _Improving lognormal models for cosmological fields_ , _MNRAS_ 459 (2016) 3693 [1602.08503].
  * [11] A. Agrawal, R. Makiya, C.-T. Chiang, D. Jeong, S. Saito and E. Komatsu, _Generating log-normal mock catalog of galaxies in redshift space_ , _J. Cosmology Astropart. Phys_ 2017 (2017) 003 [1706.09195].
  * [12] R. Makiya, I. Kayo and E. Komatsu, _Ray-tracing log-normal simulation for weak gravitational lensing: application to the cross-correlation with galaxies_ , _J. Cosmology Astropart. Phys_ 2021 (2021) 095 [2008.13195].
  * [13] F. Bernardeau, S. Colombi, E. Gaztañaga and R. Scoccimarro, _Large-scale structure of the Universe and cosmological perturbation theory_ , _Phys. Rep._ 367 (2002) 1 [astro-ph/0112551].
  * [14] R. Scoccimarro and R. K. Sheth, _PTHALOS: a fast method for generating mock galaxy distributions_ , _MNRAS_ 329 (2002) 629 [astro-ph/0106120].
  * [15] M. Manera, R. Scoccimarro, W. J. Percival, L. Samushia, C. K. McBride, A. J. Ross et al., _The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: a large sample of mock galaxy catalogues_ , _MNRAS_ 428 (2013) 1036 [1203.6609].
  * [16] M. Manera, L. Samushia, R. Tojeiro, C. Howlett, A. J. Ross, W. J. Percival et al., _The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: mock galaxy catalogues for the low-redshift sample_ , _MNRAS_ 447 (2015) 437 [1401.4171].
  * [17] G. Taffoni, P. Monaco and T. Theuns, _PINOCCHIO and the hierarchical build-up of dark matter haloes_ , _MNRAS_ 333 (2002) 623 [astro-ph/0109324].
  * [18] S. Tassev, M. Zaldarriaga and D. J. Eisenstein, _Solving large scale structure in ten easy steps with COLA_ , _J. Cosmology Astropart. Phys_ 2013 (2013) 036 [1301.0322].
  * [19] M. White, J. L. Tinker and C. K. McBride, _Mock galaxy catalogues using the quick particle mesh method_ , _MNRAS_ 437 (2014) 2594 [1309.5532].
  * [20] F. S. Kitaura, G. Yepes and F. Prada, _Modelling baryon acoustic oscillations with perturbation theory and stochastic halo biasing._ , _MNRAS_ 439 (2014) L21 [1307.3285].
  * [21] A. Izard, M. Crocce and P. Fosalba, _Ice-cola: towards fast and accurate synthetic galaxy catalogues optimizing a quasi-n-body method_ , _Monthly Notices of the Royal Astronomical Society_ 459 (2016) 2327–2341.
  * [22] C. Howlett, M. Manera and W. J. Percival, _L-PICOLA: A parallel code for fast dark matter simulation_ , _Astron. Comput._ 12 (2015) 109 [1506.03737].
  * [23] S. Avila, S. G. Murray, A. Knebe, C. Power, A. S. G. Robotham and J. Garcia-Bellido, _Halogen: a tool for fast generation of mock halo catalogues_ , _Monthly Notices of the Royal Astronomical Society_ 450 (2015) 1856–1867.
  * [24] C.-H. Chuang, F.-S. Kitaura, F. Prada, C. Zhao and G. Yepes, _EZmocks: extending the Zel’dovich approximation to generate mock galaxy catalogues with accurate clustering statistics_ , _MNRAS_ 446 (2015) 2621 [1409.1124].
  * [25] D. Alonso, P. G. Ferreira and M. G. Santos, _Fast simulations for intensity mapping experiments_ , _MNRAS_ 444 (2014) 3183 [1405.1751].
  * [26] J. Farr, A. Font-Ribera, H. du Mas des Bourboux, A. Muñoz-Gutiérrez, F. J. Sánchez, A. Pontzen et al., _LyaCoLoRe: synthetic datasets for current and future Lyman- $\alpha$ forest BAO surveys_, _J. Cosmology Astropart. Phys_ 2020 (2020) 068 [1912.02763].
  * [27] K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke et al., _HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere_ , _ApJ_ 622 (2005) 759 [astro-ph/0409513].
  * [28] S. Cole, W. J. Percival, J. A. Peacock, P. Norberg, C. M. Baugh, C. S. Frenk et al., _The 2dF Galaxy Redshift Survey: power-spectrum analysis of the final data set and cosmological implications_ , _MNRAS_ 362 (2005) 505 [astro-ph/0501174].
  * [29] F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, L. Campbell et al., _The 6dF Galaxy Survey: baryon acoustic oscillations and the local Hubble constant_ , _MNRAS_ 416 (2011) 3017 [1106.3366].
  * [30] C. Blake, T. Davis, G. B. Poole, D. Parkinson, S. Brough, M. Colless et al., _The WiggleZ Dark Energy Survey: testing the cosmological model with baryon acoustic oscillations at z= 0.6_ , _MNRAS_ 415 (2011) 2892 [1105.2862].
  * [31] A. Font-Ribera, P. McDonald and J. Miralda-Escudé, _Generating mock data sets for large-scale Lyman- $\alpha$ forest correlation measurements_, _J. Cosmology Astropart. Phys_ 2012 (2012) 001 [1108.5606].
  * [32] J. M. Le Goff, C. Magneville, E. Rollinde, S. Peirani, P. Petitjean, C. Pichon et al., _Simulations of BAO reconstruction with a quasar Ly- $\alpha$ survey_, _A &A_ 534 (2011) A135 [1107.4233].
  * [33] F.-S. Kitaura, J. Jasche and R. B. Metcalf, _Recovering the non-linear density field from the galaxy distribution with a Poisson-lognormal filter_ , _MNRAS_ 403 (2010) 589 [0911.1407].
  * [34] R. W. Hockney and J. W. Eastwood, _Computer Simulation Using Particles_. 1981.
  * [35] G. Cusin, I. Dvorkin, C. Pitrou and J.-P. Uzan, _Stochastic gravitational wave background anisotropies in the mHz band: astrophysical dependencies_ , _MNRAS_ 493 (2020) L1 [1904.07757].
  * [36] R. K. Sachs and A. M. Wolfe, _Perturbations of a Cosmological Model and Angular Variations of the Microwave Background_ , _ApJ_ 147 (1967) 73.
  * [37] M. Bartelmann and P. Schneider, _Weak gravitational lensing_ , _Phys. Rep._ 340 (2001) 291 [astro-ph/9912508].
  * [38] A. Lewis and A. Challinor, _Weak gravitational lensing of the CMB_ , _Phys. Rep._ 429 (2006) 1 [astro-ph/0601594].
  * [39] F. B. Abdalla and S. Rawlings, _Probing dark energy with baryonic oscillations and future radio surveys of neutral hydrogen_ , _MNRAS_ 360 (2005) 27 [astro-ph/0411342].
  * [40] A. Pourtsidou and R. B. Metcalf, _Gravitational lensing of cosmological 21 cm emission_ , _MNRAS_ 448 (2015) 2368 [1410.2533].
  * [41] E. Schaan, S. Ferraro and D. N. Spergel, _Weak lensing of intensity mapping: The cosmic infrared background_ , _Phys. Rev. D_ 97 (2018) 123539 [1802.05706].
  * [42] F. Villaescusa-Navarro, D. Alonso and M. Viel, _Baryonic acoustic oscillations from 21 cm intensity mapping: the Square Kilometre Array case_ , _MNRAS_ 466 (2017) 2736 [1609.00019].
  * [43] A. Witzemann, D. Alonso, J. Fonseca and M. G. Santos, _Simulated multitracer analyses with H I intensity mapping_ , _MNRAS_ 485 (2019) 5519 [1808.03093].
  * [44] H. du Mas des Bourboux, J. Rich, A. Font-Ribera, V. de Sainte Agathe, J. Farr, T. Etourneau et al., _The Completed SDSS-IV Extended Baryon Oscillation Spectroscopic Survey: Baryon Acoustic Oscillations with Ly $\alpha$ Forests_, _ApJ_ 901 (2020) 153 [2007.08995].
  * [45] A. Gabasch, U. Hopp, G. Feulner, R. Bender, S. Seitz, R. P. Saglia et al., _The evolution of the luminosity functions in the FORS deep field from low to high redshift. II. The red bands_ , _A &A_ 448 (2006) 101 [astro-ph/0510339].
  * [46] D. Alonso, J. Sanchez, A. Slosar and LSST Dark Energy Science Collaboration, _A unified pseudo-C ℓ framework_, _MNRAS_ 484 (2019) 4127 [1809.09603].
  * [47] M. Sinha and L. H. Garrison, _CORRFUNC - a suite of blazing fast correlation functions on the CPU_ , _MNRAS_ 491 (2020) 3022 [1911.03545].
  * [48] P. Schneider, L. van Waerbeke and Y. Mellier, _B-modes in cosmic shear from source redshift clustering_ , _A &A_ 389 (2002) 729 [astro-ph/0112441].
  * [49] N. Kaiser, _Clustering in real space and in redshift space_ , _Monthly Notices of the Royal Astronomical Society_ 227 (1987) 1 [https://academic.oup.com/mnras/article-pdf/227/1/1/18522208/mnras227-0001.pdf].
  * [50] The LSST Dark Energy Science Collaboration, R. Mandelbaum, T. Eifler, R. Hložek, T. Collett, E. Gawiser et al., _The LSST Dark Energy Science Collaboration (DESC) Science Requirements Document_ , _arXiv e-prints_ (2018) arXiv:1809.01669 [1809.01669].
  * [51] A. Nicola, D. Alonso, J. Sánchez, A. Slosar, H. Awan, A. Broussard et al., _Tomographic galaxy clustering with the Subaru Hyper Suprime-Cam first year public data release_ , _J. Cosmology Astropart. Phys_ 2020 (2020) 044 [1912.08209].
  * [52] R. J. Wilman, L. Miller, M. J. Jarvis, T. Mauch, F. Levrier, F. B. Abdalla et al., _A semi-empirical simulation of the extragalactic radio continuum sky for next generation radio telescopes_ , _MNRAS_ 388 (2008) 1335 [0805.3413].
  * [53] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi et al., _Planck 2018 results. VI. Cosmological parameters_ , _A &A_ 641 (2020) A6 [1807.06209].
  * [54] A. Izard, P. Fosalba and M. Crocce, _Ice-cola: fast simulations for weak lensing observables_ , _Monthly Notices of the Royal Astronomical Society_ 473 (2017) 3051–3061.
  * [55] Y. Feng, M.-Y. Chu, U. Seljak and P. McDonald, _FASTPM: a new scheme for fast simulations of dark matter and haloes_ , _MNRAS_ 463 (2016) 2273 [1603.00476].
  * [56] M. G. Santos, L. Ferramacho, M. B. Silva, A. Amblard and A. Cooray, _Fast large volume simulations of the 21-cm signal from the reionization and pre-reionization epochs_ , _MNRAS_ 406 (2010) 2421 [0911.2219].
  * [57] G. Stein, M. A. Alvarez and J. R. Bond, _The mass-Peak Patch algorithm for fast generation of deep all-sky dark matter halo catalogues and its N-body validation_ , _MNRAS_ 483 (2019) 2236 [1810.07727].
  * [58] G. Stein, M. A. Alvarez, J. R. Bond, A. van Engelen and N. Battaglia, _The Websky extragalactic CMB simulations_ , _J. Cosmology Astropart. Phys_ 2020 (2020) 012 [2001.08787].
  * [59] T. E. Oliphant, _A guide to NumPy_ , vol. 1. Trelgol Publishing USA, 2006.
  * [60] S. Van Der Walt, S. C. Colbert and G. Varoquaux, _The numpy array: a structure for efficient numerical computation_ , _Computing in Science & Engineering_ 13 (2011) 22.
  * [61] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau et al., _SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python_ , _Nature Methods_ 17 (2020) 261.
  * [62] Astropy Collaboration, T. P. Robitaille, E. J. Tollerud, P. Greenfield, M. Droettboom, E. Bray et al., _Astropy: A community Python package for astronomy_ , _A &A_ 558 (2013) A33 [1307.6212].
  * [63] Astropy Collaboration, A. M. Price-Whelan, B. M. SipHocz, H. M. G"unther, P. L. Lim, S. M. Crawford et al., _The Astropy Project: Building an Open-science Project and Status of the v2.0 Core Package_ , _aj_ 156 (2018) 123 [1801.02634].
  * [64] A. Zonca, L. Singer, D. Lenz, M. Reinecke, C. Rosset, E. Hivon et al., _healpy: equal area pixelization and spherical harmonics transforms for data on the sphere in python_ , _Journal of Open Source Software_ 4 (2019) 1298.
  * [65] J. D. Hunter, _Matplotlib: A 2d graphics environment_ , _Computing in Science & Engineering_ 9 (2007) 90.
  * [66] M. Newville, T. Stensitzki, D. B. Allen and A. Ingargiola, _LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python_ , Sept., 2014\.  10.5281/zenodo.11813.

##  Appendix A Higher-order terms in the modelling of redshift-space
distortions

The redshift-space galaxy overdensity, $\delta^{s}_{\rm LN}=\delta_{\rm
LN}(\mathbf{s})$, is related to its real-space equivalent $\delta_{\rm
LN}(\mathbf{x})$ and the normalized gradient of line-of sight velocities
$\eta=-\partial_{z}v_{z}/H(z)$ via:

$1+\delta^{s}_{\rm LN}=\frac{1+\delta_{\rm LN}}{1-\eta}$ (A.1)

Assuming Gaussian RSDs are small, we can expand this in powers of $\eta$,

$\delta^{s}_{\rm LN}=\delta_{\rm LN}+\eta+\epsilon~{},$ (A.2)

where $\epsilon(\mathbf{x})\equiv\delta_{\rm
LN}(\mathbf{x})~{}\eta(\mathbf{x})$. This term would be order-2 assuming
$\delta_{\rm LN}$ is small, but we do not make this approximation here.
Ignoring this term one would recover Eq. 3.2, used for the theoretical
predictions in Section 3.1.3. If we keep this extra term, however, the model
for the redshift-space galaxy power spectrum will have new contributions with
respect to the model described in Eq. 3.3:

$\displaystyle P_{\rm LN}^{s}(k,\mu)=$ $\displaystyle~{}\left(\Delta
k\right)^{3}\langle|\delta^{s}_{\rm LN}(\mathbf{k})|^{2}\rangle$
$\displaystyle=$ $\displaystyle~{}P_{\rm LN}(k)+f^{2}\mu^{4}P_{\rm
G}(k)+2bf\mu^{2}P_{\rm G}(k)$ (A.3) $\displaystyle~{}+\left(\Delta
k\right)^{3}\left[2\langle\delta_{\rm
LN}(\mathbf{k})~{}\epsilon(\mathbf{k})\rangle+2\langle\eta(\mathbf{k})~{}\epsilon(\mathbf{k})\rangle+\langle\epsilon(\mathbf{k})~{}\epsilon(\mathbf{k})\rangle\right]~{}.$

Figure 16: Measurements and predictions for the different cross-correlations
discussed in A for a single simulation with two special tracers ($b_{A}=0.001$
and $b_{B}=2$). The predictions are from the simpler model discussed in 3.1.3
and are missing those terms involving $\epsilon(\mathbf{x})=\delta_{\rm
LN}(\mathbf{x})\eta(\mathbf{x})$. This explains the disagreement seen on the
small-scales quadrupole in panel (c), it also shows a small disagreement in
panels (e) and (f).

Even though it is possible to compute analytical predictions for the
correlation function including these new $\epsilon$ terms (as convolutions of
$\delta_{\rm LN}$ and $\eta$ in Fourier space), we leave this for future work.
Here we only quantify each of the terms by looking at cross-correlations of
test galaxy samples from a custom CoLoRe simulation designed for this study.
This simulation was similar to those described in 3.1.1 but with two
spectroscopic galaxy samples: a first sample $\delta_{A}$ with an extremely
low clustering amplitude ($b_{A}=0.001$), and a second sample $\delta_{B}$
with a redshift-independent bias of $b_{B}=2$.

While the clustering of sample $A$ in real space is negligible, its redshift-
space power spectrum can be modeled as
$\langle\delta_{A}^{s}(\mathbf{k})\delta_{A}^{s}(\mathbf{k})\rangle\propto
f^{2}\mu^{4}P_{\rm G}(k)$ (up to a normalisation factor $(\Delta k)^{3}$).
This can be clearly seen in the (a) panel of Figure 16, and validates the
simulation of RSDs in CoLoRe. The clustering of sample $B$ in real (redshift)
space is shown in the (b) ((c)) panels of the same figure, and are similar to
Fig. 11 discussed in 3.1.3. While its real space clustering is well described
by the lognormal model
$\langle\delta_{B}(\mathbf{k})\delta_{B}(\mathbf{k})\rangle\propto P_{\rm
LN}(k)$, the small-scales quadrupole of its redshift-space equivalent can not
be described by the simple model of Eq. 3.3. We are missing the contributions
from the terms $\langle\delta_{\rm LN}~{}\epsilon\rangle$,
$\langle\eta~{}\epsilon\rangle$ and $\langle\epsilon~{}\epsilon\rangle$
introduced in Eq. A.

In panel (d) we show the cross-correlation of $\delta_{A}^{s}$ and
$\delta_{B}$, compared to its prediction
$\langle\delta_{A}^{s}(\mathbf{k})\delta_{B}(\mathbf{k})\rangle\propto
b_{B}f\mu^{2}P_{\rm G}(k)$. Because there is no $\epsilon$ term involved in
this cross-correlation, the model describes the measurement very well on all
scales. In panel (e) we cross-correlate $\delta_{A}^{s}$ with $\delta_{B}^{s}$
instead, and compare it to a theoretical prediction that includes the two
Kaiser terms ($b_{B}f\mu^{2}P_{\rm G}(k)+f^{2}\mu^{4}P_{\rm G}(k)$) but is
missing an extra term $\langle\eta~{}\epsilon\rangle$. The minor disagreement
of the quadrupole on small scales allows us to estimate the magnitude and sign
of the missing term. Finally, in (f) we show the cross-correlations of
$\delta_{B}$ and $\delta_{B}^{s}$, i.e., the cross-correlation of the same B
sample in real and in redshift space. We plot its prediction from the simpler
model from section 3.1.3, including the terms $P_{\rm
LN}(k)+b_{B}f\mu^{2}P_{\rm G}(k)$. Again, following Eq. A this cross-
correlation should include an extra term $\langle\delta_{\rm
LN}~{}\epsilon\rangle$ that explains the small disagreement of the quadrupole
on small scales.

Figure 17: Contributions to the monopole (top panel/blue lines and bands) and
the quadrupole (bottom panel/red lines and bands) from the terms in Eqs.
A.4-A.6 using the same simulation as in Figure 16. The shaded bands show the
error in the measurement of each term from the scatter between 48 healpixels
of $N_{\rm{side}}=2$. The dashed lines show the full model prediction as in
section 3.1.3, where the extra terms are not included. Particularly the term
$\langle\epsilon~{}\epsilon\rangle$ (right panel) has an important impact on
the small-scales quadrupole.

Following Eq. A.1 it is clear that
$\epsilon=\delta_{B}^{s}-\delta_{B}-\delta_{A}^{s}$. Therefore, by combining
the different cross-correlations discussed above we can now isolate each of
the new terms in A:

$\displaystyle\langle\delta_{\rm LN}~{}\epsilon\rangle=$
$\displaystyle~{}\langle\delta_{B}~{}\delta_{B}^{s}\rangle-\langle\delta_{B}~{}\delta_{B}\rangle-\langle\delta_{B}~{}\delta_{A}^{s}\rangle$
(A.4) $\displaystyle\langle\eta~{}\epsilon\rangle=$
$\displaystyle~{}\langle\delta_{A}^{s}~{}\delta_{B}^{s}\rangle-\langle\delta_{A}^{s}~{}\delta_{B}\rangle-\langle\delta_{A}^{s}~{}\delta_{B}\rangle$
(A.5) $\displaystyle\langle\epsilon~{}\epsilon\rangle=$
$\displaystyle\langle\delta_{B}^{s}~{}\delta_{B}^{s}\rangle+\langle\delta_{B}~{}\delta_{B}\rangle+\langle\delta_{A}^{s}~{}\delta_{A}^{s}\rangle+2\langle\delta_{A}^{s}~{}\delta_{B}\rangle-2\langle\delta_{A}^{s}~{}\delta_{B}^{s}\rangle-2\langle\delta_{B}~{}\delta_{B}^{s}\rangle~{}.$
(A.6)

The multipoles corresponding to these three combinations of correlations are
shown in Figure 17, together with the model prediction from equation 3.3. One
can see that the largest correction to the model should come from the
$\langle\epsilon~{}\epsilon\rangle$ term, while the other two terms are small
and partially cancel each other.