#  Euclid: Modelling massive neutrinos in cosmology – a code comparison***This
paper is published on behalf of the Euclid Consortium

J. Adamek *  O R. E. Angulo *  O C. Arnold  M. Baldi *  O M. Biagetti *  O B.
Bose *  O C. Carbone  T. Castro *  O J. Dakin *  O K. Dolag  W. Elbers *  O C.
Fidler  C. Giocoli *  O S. Hannestad  F. Hassani *  O C. Hernández-Aguayo *  O
K. Koyama *  O B. Li *  O R. Mauland *  O P. Monaco *  O C. Moretti *  O D. F.
Mota  C. Partmann  G. Parimbelli *  O D. Potter *  O A. Schneider *  O S.
Schulz *  O R. E. Smith *  O V. Springel *  O J. Stadel *  O T. Tram *  O M.
Viel *  O F. Villaescusa-Navarro *  O H. A. Winther *  O B. S. Wright *  O M.
Zennaro *  O N. Aghanim  L. Amendola  N. Auricchio *  O D. Bonino  E.
Branchini *  O M. Brescia *  O S. Camera *  O V. Capobianco *  O V. F. Cardone
J. Carretero *  O F. J. Castander  M. Castellano *  O S. Cavuoti *  O A.
Cimatti  R. Cledassou *  O G. Congedo *  O L. Conversi *  O Y. Copin *  O A.
Da Silva *  O H. Degaudenzi *  O M. Douspis  F. Dubath *  O C. A. J. Duncan
X. Dupac  S. Dusini *  O S. Farrens *  O S. Ferriol  P. Fosalba *  O M.
Frailis *  O E. Franceschi *  O S. Galeotta *  O B. Garilli *  O W. Gillard *
O B. Gillis *  O A. Grazian *  O S. V. Haugan *  O W. Holmes  A. Hornstrup *
O K. Jahnke *  O S. Kermiche *  O A. Kiessling *  O M. Kilbinger *  O T.
Kitching *  O M. Kunz *  O H. Kurki-Suonio *  O P. B. Lilje *  O I. Lloro  O.
Mansutti *  O O. Marggraf *  O F. Marulli  R. Massey *  O E. Medinaceli *  O
M. Meneghetti *  O G. Meylan  M. Moresco *  O L. Moscardini *  O E. Munari *
O S.-M. Niemi  C. Padilla *  O S. Paltani  F. Pasian  K. Pedersen  W. J.
Percival *  O V. Pettorino  G. Polenta *  O M. Poncet  L. A. Popa  F. Raison *
O R. Rebolo  A. Renzi *  O J. Rhodes  G. Riccio  E. Romelli *  O M. Roncarelli
R. Saglia *  O D. Sapone *  O B. Sartoris  P. Schneider *  O T. Schrabback *
O A. Secroun *  O G. Seidel *  O C. Sirignano *  O G. Sirri *  O L. Stanco *
O J.-L. Starck  P. Tallada-Crespí  A. N. Taylor  I. Tereno  R. Toledo-Moreo *
O F. Torradeflot *  O I. Tutusaus *  O L. Valenziano *  O T. Vassallo *  O Y.
Wang *  O J. Weller *  O A. Zacchei *  O G. Zamorani *  O J. Zoubian  G.
Fabbian *  O and V. Scottez

###### Abstract

The measurement of the absolute neutrino mass scale from cosmological large-
scale clustering data is one of the key science goals of the Euclid mission.
Such a measurement relies on precise modelling of the impact of neutrinos on
structure formation, which can be studied with $N$-body simulations. Here we
present the results from a major code comparison effort to establish the
maturity and reliability of numerical methods for treating massive neutrinos.
The comparison includes eleven full $N$-body implementations (not all of them
independent), two $N$-body schemes with approximate time integration, and four
additional codes that directly predict or emulate the matter power spectrum.
Using a common set of initial data we quantify the relative agreement on the
nonlinear power spectrum of cold dark matter and baryons and, for the $N$-body
codes, also the relative agreement on the bispectrum, halo mass function, and
halo bias. We find that the different numerical implementations produce fully
consistent results. We can therefore be confident that we can model the impact
of massive neutrinos at the sub-percent level in the most common summary
statistics. We also provide a code validation pipeline for future reference.

##  1 Introduction

The upcoming Euclid mission [1] will provide very detailed observations of the
large-scale structure of our Universe, making it possible to probe physics
related to dark energy and neutrinos at an unprecedented level of precision.
The analysis and interpretation of these data require a very accurate
modelling of the process of structure formation. This is of particular
relevance since a precise modelling of the mass-dependent effect neutrinos
have on various summary statistics will allow a cosmological measurement of
the absolute neutrino mass scale, which is one of the key science objectives
of Euclid.

Here we focus on the treatment of massive neutrinos in cosmological $N$-body
simulations and investigate the convergence of a number of different codes
over a variety of different scales and redshifts. Oscillation experiments have
established a firm lower bound on the sum of neutrino masses of around
$0.06\,\mathrm{eV}$. Using the well-known relation between neutrino mass and
cosmological energy density, $\sum m_{\nu}\approx\Omega_{\nu}h^{2}\times
94\,\text{eV}$, this lower bound on the sum of neutrino masses corresponds to
a lower bound of $\Omega_{\nu}h^{2}>6\times 10^{-4}$, or approximately 0.5% of
the total matter density. As usual, the cosmological density of any component
X can be given in terms of $\Omega_{\text{X}}$, which is its present-day
energy density in units of the critical density, and physical density
parameters are then denoted as $\omega_{\text{X}}=\Omega_{\text{X}}h^{2}$,
where $h$ is the reduced Hubble parameter.

The sum of the neutrino masses is already constrained using information from
the cosmic microwave background (CMB) combined with observations of large-
scale structure like baryon acoustic oscillations [2], redshift-space
distortions [3], and the Lyman-$\alpha$ forest [4]. While currently providing
only upper bounds, these constraints are expected to improve significantly
with upcoming surveys like Euclid which will be able to measure the neutrino
mass fraction even if it is close to the lower bound. The matter power
spectrum is still affected at the level of $4\,\%$ in this scenario, which is
well within the sensitivity of the Euclid main probes.

For instance, the weak-lensing signal probed by Euclid is sensitive to the
matter power spectrum up to $k\approx 7\,h\,\mathrm{Mpc}^{-1}$ [5]. A linear
model would be completely inadequate at such short scales and we therefore
need robust nonlinear models. The forecast for galaxy clustering in Euclid
typically assumes that the matter power spectrum and the galaxy bias is well
understood at least up to $k\approx 0.3\,h\,\mathrm{Mpc}^{-1}$ which also
requires some nonlinear prescription [6]. Many codes considered here are used
in Euclid preparation papers and reference simulations. For instance, PKDGRAV3
ran Euclid’s flagship simulations and the simulations used by Knabenhans et
al. [7, 8] to train the EuclidEmulator2; openGADGET3 was used to calibrate the
halo mass function for the Euclid cluster abundance analysis by Castro et al.
[9], and PINOCCHIO was used to create the synthetic catalogues for the
validation of the covariance matrix of cluster abundance and the clustering of
clusters by Fumagalli et al. [10, 11]. Our objective is therefore to establish
a reliable calibration baseline for the measurement of the neutrino mass scale
within the cosmological analysis of Euclid data. We expect that our results
are also relevant in the context of other so-called “stage IV surveys” like
the Vera C. Rubin Observatory or the Nancy Grace Roman Space Telescope.

Over the past decade, a variety of different methods have been employed to
incorporate massive neutrinos in $N$-body simulations [12]. As discussed in
more detail later, they broadly fall into two categories which we may refer to
as “particle based” and “mesh based,” respectively. They follow different
philosophies of keeping track of the evolving neutrino phase-space
distribution function. Particle-based methods sample the six-dimensional
phase-space directly, see e.g. Refs. [13, 14, 15, 16, 17, 18, 19, 20, 21, 22].
Mesh-based methods, on the other hand, work under the approximation that
neutrino perturbations remain small and can be treated with perturbation
theory [23, 24, 25, 26]. In the simplest case that is sufficient for low
neutrino masses one works with a linear realisation of the neutrino density
field on a grid [23]. The linear theory for neutrinos may also be solved using
the full nonlinear gravitational potential calculated in the simulation [24,
25, 26]. A different approach is to treat the massive neutrinos as a fluid and
then solve the corresponding fluid equations, employing some approximation
scheme to close the set of equations [27]. One can also attempt to integrate
the Vlasov–Poisson equations on a six-dimensional phase-space grid [28].

Generally speaking, the mesh-based schemes work best for relatively small
neutrino masses where neutrino perturbations remain linear or quasi-linear at
all times. There are also hybrid schemes that use elements from both
approaches. For instance, one may use a linear mesh-based representation at
early times which is then converted to a particle representation at late time
[29, 30]. The so-called “$\delta f$ method” introduced in Elbers et al. [31]
is another hybrid approach that uses a particle ensemble to estimate
perturbations $\delta f$ to the smooth background phase-space distribution
function $f$. Finally, a coordinate (gauge) transformation can be used to
include linear neutrino perturbations without modifications to the $N$-body
simulation code [32]. Cosmology-rescaling algorithms that are applied in post-
processing have been shown to provide accurate results as well [33]. In this
work, we aim to compare these different numerical approaches by employing them
to run cosmological $N$-body simulations starting from the same initial
conditions, and comparing the properties of the resulting matter and halo
distributions using a controlled post-processing pipeline.

This paper is structured as follows. We begin with a brief review of neutrino
physics and its impact on cosmology in Section 2. In Section 3, we describe
the numerical methods that can be used to account for the cosmological effects
of neutrinos. Our simulations are described in Section 4 and in Section 5 we
present our numerical results. We conclude with a discussion in Section 6.

##  2 Neutrino physics

From oscillation experiments it is firmly established that at least two of the
standard-model neutrino mass states have non-zero mass, but the absolute mass
scale is unknown and two mass orderings remain possible: normal and inverted.
The current best-fit values for the mass-square differences measured in
oscillation experiments are given by [34]

$\displaystyle\Delta m_{21}^{2}$ $\displaystyle=$ $\displaystyle
7.42^{+0.21}_{-0.20}\times 10^{-5}\,{\rm eV}^{2},$ (2.1) $\displaystyle\Delta
m_{31}^{2}$ $\displaystyle=$ $\displaystyle 2.517^{+0.026}_{-0.028}\times
10^{-3}\,{\rm eV}^{2}\quad(\textrm{NO}),$ (2.2) $\displaystyle\Delta
m_{32}^{2}$ $\displaystyle=$ $\displaystyle-2.498^{+0.028}_{-0.028}\times
10^{-3}\,{\rm eV}^{2}\quad(\textrm{IO}),$ (2.3)

where “NO” denotes the normal mass ordering and “IO” the inverted one. This
leads to lower bounds on the sum of neutrino masses of $\sum m_{\nu}\gtrsim
0.06$ eV (NO) and $\sum m_{\nu}\gtrsim 0.1$ eV (IO), respectively. The best
current experimental upper bound on the neutrino masses comes from the KATRIN
experiment, which measures an incoherent sum of mass states using beta decay
of tritium [35, 36]. This bound approximately translates to $\sum
m_{\nu}\lesssim 2.4$ eV.

However, cosmology already provides much more stringent bounds, typically
around $\sum m_{\nu}\lesssim 0.1-0.2$ eV. Assuming a minimal $\Lambda$CDM
cosmology with massive neutrinos, a joint analysis of cosmological probes
currently obtains $\sum m_{\nu}<0.09$ eV at $95\%$ confidence, already putting
pressure on the inverted mass ordering scenario [37]. These constraints are
based on several physical effects affecting the CMB and large-scale structure
in different ways, see Archidiacono et al. [38] for a review. Some of the
effects are simply related to the change in the expansion history, others to
the explicit coupling of neutrinos to cosmological perturbations. More
specifically, massive neutrinos modify the shape of the matter power spectrum
both in the linear and the nonlinear regimes. First, as neutrinos behave like
radiation in the early Universe, they move the radiation-matter equality to
slightly later times, therefore shifting the peak of the power spectrum
towards smaller wavenumbers. Second, after the non-relativistic transition,
they slow down the linear growth of perturbations at scales smaller than the
free-streaming length, leading to a scale-dependent growth rate. The small-
scale suppression in the linear power spectra of cold dark matter (CDM) and
baryons, $P_{\mathrm{cb}}$, or total matter (which includes massive
neutrinos), $P_{\mathrm{m}}$, with respect to a model with massless neutrinos,
can be quantified as [18, 39, 40]

$\frac{\Delta P_{\mathrm{cb}}}{P_{\mathrm{cb}}}\approx
6\,f_{\nu}\,,\qquad\frac{\Delta P_{\mathrm{m}}}{P_{\mathrm{m}}}\approx
8\,f_{\nu}\,,$ (2.4)

respectively, where $f_{\nu}=\Omega_{\nu}/\Omega_{\mathrm{m}}$ is the neutrino
mass fraction. In the nonlinear regime, this suppression is even more
prominent and exhibits a dip at $k\approx 1\ h\ \mathrm{Mpc}^{-1}$ for low
redshift, giving rise to the well-known “spoon-like” feature [41]. In the
context of the halo model, this feature appears at the transition region of
the two-halo and the one-halo terms. In particular, the dip is caused by the
small-scale suppression of two-halo clustering that is induced by free-
streaming, while the subsequent rise reflects the fact that the number of the
most massive halos is rather independent of the neutrino masses. All these
effects can be accurately predicted by modelling the neutrino component in
cosmological $N$-body simulations. Assessing the relative accuracy and
convergence of such modelling over a range of different numerical methods and
simulation codes is the main goal of this paper.

##  3 Numerical methods

In this section, we give an overview of the various methods that have been
developed for the treatment of massive neutrinos in $N$-body simulations and
other numerical models. The most accurate results are expected when the local
density of neutrinos in configuration space is accounted for within a
simulation itself. This is technically challenging because of the large phase-
space volume that is populated by neutrino particles. The methods to deal with
this broadly fall into two categories that shall be discussed in turn:
particle-based and mesh-based. Hybrid methods that use concepts from both
categories have also been developed. Apart from full $N$-body simulations that
may try to incorporate (as much as possible) the neutrino physics, there also
exist approximate methods to generate realisations of large-scale structure.
These can be augmented with recipes to account for the effect of massive
neutrinos. Finally, if one is only interested in summary statistics like the
power spectrum, emulators are a powerful tool that can be calibrated to
include the sum of the neutrino masses as a free parameter. An overview of the
various numerical codes employed in this work is given in Table 1.

Table 1: Overview of the numerical codes employed in this code comparison.
$N$-body codes typically use a particle-mesh (PM) method coupled to some
scheme to increase the force resolution in high-density regions. Methods
featured here include Tree-PM, adaptive mesh refinement (AMR), fast multipole
method (FMM), particle-particle/particle-mesh (P3M), and moving mesh. No
adaptive force computation is used in gevolution, COLA, and PINOCCHIO, which
all work with a uniform mesh.
Code | type | neutrino method(s) | reference(s)
---|---|---|---
GADGET-3 |  $N$-body (Tree-PM) | particle | [42, 43]
L-GADGET3 |  $N$-body (Tree-PM) | mesh | [43, 44, 45]
openGADGET3 |  $N$-body (Tree-PM) | particle | [46, 47]
GADGET-4 |  $N$-body (Tree-PM) | particle | [48]
NM-GADGET4 |  $N$-body (Tree-PM) | Newtonian motion gauge | [32, 49]
AREPO |  $N$-body (moving mesh) | particle | [50, 51]
CONCEPT |  $N$-body (P3M) | mesh | [52, 27]
PKDGRAV3 |  $N$-body (Tree + FMM) | mesh | [53]
SWIFT |  $N$-body (PM + FMM) | particle / $\delta f$ | [54, 31, 55]
ANUBIS |  $N$-body (PM + AMR) | particle | [56]
gevolution |  $N$-body (uniform PM) | particle / mesh | [20, 57, 58]
COLA |  $N$-body surrogate | mesh | [59, 60]
PINOCCHIO |  $N$-body surrogate | linear growth factor | [61, 62]
ReACT |  $P(k)$ prediction | halo-model reaction | [63, 64]
BACCOemulator |  $P(k)$ prediction | emulation | [45]
EuclidEmulator2 |  $P(k)$ prediction | emulation | [8]
Cosmic Emu |  $P(k)$ prediction | emulation | [65, 66]

###  3.1 Particle-based methods

Conceptually, the most straightforward method of accounting for cosmic
neutrinos in a simulation is to represent them by a separate $N$-body
ensemble. However, this method faces several challenges that need to be
addressed carefully.

The first challenge is posed by the aforementioned phase-space volume that
needs to be sampled. The phase-space distribution of neutrino particles
typically has a very large velocity dispersion that is orders of magnitude
larger than the bulk velocity. Therefore, representing the neutrinos with a
small number of $N$-body particles that simply track the bulk velocities,
which is essentially the method of choice for CDM or baryons, would completely
miss the fact that most neutrinos are unbound and easily free stream out of
gravity wells. Thus, the common practice is to sample the $N$-body particles
from the true phase-space distribution, effectively performing a Monte-Carlo
integration of the evolution equations. The main drawback of this method is
that a poor sampling usually introduces significant shot noise, while high
sampling rates quickly become very expensive as both the memory requirement
and computations become completely dominated by the neutrino particle load.
This is undesirable since neutrinos are just a tiny fraction of the matter
after all.

Shot noise affects all moments of the distribution function, and in particular
the density. This means that shot noise will also propagate into the
gravitational field computed from the neutrino perturbations. While this could
in principle severely degrade the accuracy of the gravitational evolution as a
whole, the impact is actually mitigated by the fact that neutrinos only
account for a very small fraction of the total matter, and the gravitational
fields are therefore dominated by cold species. Nevertheless, some shot noise
does propagate into the other matter species, particularly on large scales
where the contribution of neutrino perturbations is largest. Various
strategies have been developed to reduce the impact of shot noise, e.g. by
filtering small-scale fluctuations. One effective strategy is to implement a
statistical weighting of the neutrino particles, as is done in the $\delta f$
method [31]. This method works by decomposing the distribution function $f$
into an analytical background component $\bar{f}$ and a perturbation $\delta
f$ computed from the particle ensemble. The weights are given at each time
step by requiring phase-space density conservation. These weights are
typically negligible, except for particles that are significantly perturbed,
such as those captured by halos. Shot noise is thereby minimised as particles
only contribute to the gravitational potential when needed. The $\delta f$
method has been implemented in a number of codes, but is only used by SWIFT in
this comparison (see Table 1). Other strategies aimed at reducing shot noise
while still using particles include alternate sampling of neutrino momenta
[21] and various hybrid methods [29, 67, 30].

The second challenge pertains to the kinematics of neutrino $N$-body
particles. If one were to apply a similar time-stepping criterion as for cold
matter species, the high velocities would typically result in extremely small
integration time steps, making the simulations considerably more expensive.
This is often solved by relaxing the time-stepping criterion and allowing the
neutrino particles to travel a larger distance in each integration step than
what would be considered “safe” for cold species. High-velocity neutrino
particles then respond poorly to small-scale fluctuations in the gravitational
forces. Given that the neutrino distribution at small scales is plagued by the
shot-noise problem anyway, this additional problem is often considered to be
of little concern.

Still regarding kinematics, further issues arise due to the relativistic
nature of neutrinos. In the weak-field limit, the propagation of a
collisionless massive particle is governed by the Hamiltonian equations of
motion [68]

$\displaystyle\bm{v}$ $\displaystyle=$
$\displaystyle\frac{\bm{p}c}{\sqrt{\bm{p}^{2}+m^{2}c^{2}a^{2}}}\,,$ (3.1)
$\displaystyle\bm{p}^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{2\bm{p}^{2}+m^{2}c^{2}a^{2}}{c\sqrt{\bm{p}^{2}+m^{2}c^{2}a^{2}}}\nabla\Psi\,,$
(3.2)

where $\bm{v}$ is the peculiar velocity, $\bm{p}$ is the canonical momentum,
$m$ is the particle’s rest-mass, $a$ is the scale factor, and $\Psi$ is the
gravitational potential, assuming as usual that gravitational slip can be
neglected, i.e. that non-relativistic and ultra-relativistic particles
essentially see the same potential. Here, a prime denotes the derivative with
respect to conformal time and $c$ denotes the speed of light. The canonical
momentum is conserved in the absence of a gravitational force. For non-
relativistic particles one usually considers the limit $\bm{p}^{2}\ll
m^{2}c^{2}a^{2}$ in which the equations simplify to

$\displaystyle\bm{v}$ $\displaystyle=$ $\displaystyle\frac{\bm{p}}{ma}\,,$
(3.3) $\displaystyle\bm{p}^{\prime}$ $\displaystyle=$ $\displaystyle-
ma\nabla\Psi\,.$ (3.4)

This simpler set of equations has the advantage that it is easy to find
integration methods that are symplectic, i.e. that preserve the phase-space
volume exactly as demanded by Hamiltonian time evolution. Note also that, in
the absence of a gravitational force, the peculiar velocity scales exactly as
$\propto a^{-1}$ in this case.

Even though the evolution of high-momentum particles suffers from severe
errors, including the breakdown of causality for $\bm{p}^{2}>m^{2}c^{2}a^{2}$,
some implementations might still use the simplified equations. The propagation
of these errors into the clustering amplitude of matter is limited by the fact
that high-momentum particles barely contribute to clustering in the first
place. However, using the more accurate Eqs. (3.1) and (3.2) is a more common
choice. The issue of symplectic time integration in this case is discussed
e.g. in Appendix A of Adamek et al. [20] and Appendix D of Elbers et al. [31].

###  3.2 Mesh-based methods

The main alternative to particle-based methods is to represent the
distribution function on a spatial mesh. Yoshikawa et al. [28] discretise the
distribution function on a six-dimensional mesh in phase-space, however, a
brute-force approach like this is rather expensive and cannot easily be
applied to very large simulations where memory requirements are a particular
concern. A possible way out is to take moments of the distribution function
(density, bulk velocity, and so on) where the momentum coordinates are
integrated out so that a discretisation in the three spatial dimensions is
sufficient.

Experience from solving the hierarchy of moments (the so-called Boltzmann
hierarchy) in linear perturbation theory shows that a considerable number of
moments must be taken into account in order to reach good accuracy for the
evolution of the lowest moments. This concerns in particular the density that
also affects the clustering of other matter components through gravitational
coupling. On the other hand, the numerical solutions are readily available in
the linear regime where they can be expressed in terms of linear transfer
functions. Throughout this paper, we follow the convention from standard
linear cosmological perturbation theory where the transfer function $T_{X}$ of
any perturbation variable $X$ in Fourier space is defined through the relation

$X(\bm{k},z)=T_{X}(k,z)\zeta(\bm{k})\,,$ (3.5)

where $\zeta(\bm{k})$ is the comoving curvature perturbation of the mode
$\bm{k}$ before it enters the horizon. Those transfer functions provide a
deterministic factor by which any given initial random perturbation mode needs
to be multiplied in order to obtain, for instance, the density perturbation at
any given time. Using these transfer functions that can be calculated at the
outset, a simulation code can therefore construct the linear density field of
neutrinos at any point in time for any given realisation of the random initial
conditions. This is precisely what basic mesh-based methods do: they use the
density field of neutrinos extrapolated from linear perturbation theory, which
is often a reasonable approximation because neutrinos do not cluster strongly.
The method is free of shot noise and is by construction exact in the limit of
linear perturbations. However, it obviously lacks any response to nonlinear
gravitational potentials that develop in a simulation.

While the linear method produces results that are sufficiently accurate for
many purposes, some more advanced approaches have been developed in attempts
to address the shortcomings. Dakin et al. [27] employ the coupled evolution
equations for the lowest moments in their nonlinear form. Then, to avoid
having to calculate a large Boltzmann hierarchy for every wavevector
represented in the simulation, the hierarchy is truncated by assuming that a
“scaling” holds approximately for ratios of higher moments, where the scaling
coefficients are taken from linear theory. This approximation is then used to
close the system of equations using only a small number of nonlinear moments.
By construction this method agrees with the simpler method in the limit of
linear perturbations. While the resulting nonlinear neutrino density is
somewhat more realistic, the distribution function still has some residual
errors that cannot easily be reduced without including further moments in the
nonlinear computation.

###  3.3 Approximations and other methods

For some purposes, such as the computation of covariance matrices for
different cosmological probes, it is useful to have methods for making
cosmological predictions that are faster, although less accurate, than the
traditional $N$-body methods discussed so far. Here we present two such
methods that can be used as surrogates for $N$-body simulations: the COmoving
Lagrangian Acceleration (COLA) approach and the PINpointing Orbit-Crossing
Collapsed HIerarchical Objects (PINOCCHIO) approach. In both cases, a speed-up
is achieved by drastically simplifying the time integration in the particle
evolution. Finally, we also present a method that avoids the need to include
any neutrino physics in the actual $N$-body simulation altogether, apart from
in the background solution. This method employs the so-called Newtonian motion
gauge and can be used with virtually any numerical scheme that solves the
Newtonian gravity problem.

####  3.3.1 COLA

The COLA approach by Tassev et al. [69] produces fast, approximate simulations
of cosmological structure formation. Essentially, instead of solving for a
full particle trajectory $\bm{x}(t)$, in this approach we solve for the
deviations of the full trajectory about the trajectory predicted by second-
order Lagrangian perturbation theory (2LPT)
$\delta\bm{x}(t)=\bm{x}(t)-\bm{x}_{\rm 2LPT}(t)$. Since the evolution of the
particles on large scales will be very close to that predicted by 2LPT, we can
decrease the number of time steps of the simulation to trade accuracy at small
scales for overall simulation speed while maintaining good accuracy at large
scales. For a large number of time steps, the method effectively converges to
a standard PM $N$-body method.

Adding massive neutrinos to the COLA method was described by Wright et al.
[60], which also included an implementation in the MG-PICOLA simulation
code222https://github.com/HAWinther/MG-PICOLA-PUBLIC. by Winther et al. [59].
This implementation was carried over to the COLA solver within
FML333https://github.com/HAWinther/FML/tree/master/FML/COLASolver. which
succeeded MG-PICOLA. It is this implementation of the COLA solver within FML
that we use in this paper. These implementations rely on the linear mesh-based
method described above, i.e. we use the density field of neutrinos
extrapolated from linear perturbation theory on a mesh for the PM part. For
the 2LPT part of the COLA code, we make a further approximation to the 2LPT
equation and use the $\Lambda$CDM kernel to speed up the computation.

To demonstrate the key advantage of COLA over traditional $N$-body codes, we
use only $50$ time steps linearly distributed in scale factor for the COLA
simulations in this paper. However, the COLA method does not work well for
simulations starting from high initial redshifts when using a relatively small
number of time steps; for a discussion on how to optimise initial redshift and
number of time steps in COLA simulations see Sections 4.1 and 4.3 of Izard et
al. [70]. Therefore, we use a slightly modified444The modifications are to the
order in which (pseudo-)random numbers are drawn for the phases and
amplitudes, such that we now do the same as the version of N-GenIC used for
the other methods in this paper, see Section 4.1. version of FML’s built-in
generator of initial conditions to generate initial particle data at $z=19$
instead of $z=127$ as is described in Section 4.1 and used for the other
methods in this paper. In addition, we use the CAMB Boltzmann solver by Lewis
et al. [71, 72] to generate the density transfer function for massive
neutrinos. Finally, we note that for all COLA simulations in this paper we use
a force grid that is a factor of three finer than the mean inter-particle
distance; for a thorough investigation of the impact of varying this factor in
COLA simulations see Section 4.4 of Izard et al. [70].

####  3.3.2 PINOCCHIO

The PINOCCHIO code555https://github.com/pigimonaco/Pinocchio. [73, 61, 74] is
an approximate method to generate halo catalogues in a very small fraction (of
the order of $1/1000$) of the time taken by an equivalent $N$-body simulation.
Starting from a linear density field generated in Lagrangian space over a
regular grid, its main goal is to construct halo catalogues by predicting
which particles will end up in dark matter halos. To achieve this goal the
algorithm first smoothes the linear density on a grid of smoothing radii, then
uses ellipsoidal collapse to compute the collapse time of each particle. In
the second step, it proceeds to group the collapsed particles into halos,
using an algorithm that mimics their hierarchical clustering and distinguishes
between halos and filaments.

As opposed to the other $N$-body methods employed in this work, PINOCCHIO does
not integrate particle orbits but places halos at their final position using a
single 2LPT or 3LPT displacement. Indeed, once the displacement fields are
averaged over the multi-stream region that corresponds to a dark matter halo,
LPT is very effective in predicting halo positions [75]. Another difference is
that PINOCCHIO does not start from a set of displaced particles but generates
the linear density field internally. Having implemented the same sequence for
populating modes in $k$-space as the N-GenIC code that is used for the purpose
of generating initial data in this work (see Section 4.1), it can reproduce
the same large-scale structure if the same seed for random numbers is
provided.

The extension of PINOCCHIO to massive neutrinos is presented by Rizzo et al.
[62] and is based on the result of Castorina et al. [76, 18] that the
nonlinear clustering of massive neutrinos is negligible. We use CAMB to
compute linear power spectra in massive neutrino cosmologies, and compute the
scale-dependent growth rate of matter by taking ratios of power spectra of CDM
and baryons (i.e. without neutrinos) at different times. We also adapt the
code to incorporate a scale-dependent growth rate. With respect to the
original implementation of Rizzo et al. [62], which was limited to 2LPT, we
extend here the computation to third order: as shown by Munari et al. [74]
this results in a significant improvement at mildly nonlinear scales.

Although the code has been conceived to predict the properties of dark matter
halos, it can produce a full nonlinear density field as follows: particles
that do not belong to halos are moved to their final position using 3LPT, halo
particles are distributed around their halo center of mass following a
Navarro–Frenk–White (NFW) profile [77] with Maxwellian velocity distributions.
This allows us to construct snapshots like an $N$-body simulation,
representing density fields that are far more accurate than a straight LPT
implementation. Because we have only one type of particle, to compute the
power spectrum of CDM and baryons needed below (Section 5) we subtract the
linear neutrino contribution from the total matter power spectrum obtained
from the snapshot. To this end we also assume that the neutrino-matter cross-
power spectrum, $P_{\nu,\mathrm{m}}(k)$, can be approximated by
$P_{\nu,\mathrm{m}}(k)=\sqrt{P_{\nu}^{\mathrm{L}}(k)P_{\mathrm{m}}(k)}$, where
the superscript “L” denotes a power spectrum from linear theory. This
approximation is strictly only true in the linear regime but we apply it at
all scales. While this is not the ideal configuration for which the code has
been conceived, the merit of showing results for such a density field is to
test how accurate a sophisticated version of 3LPT can be. In any case, we do
not expect PINOCCHIO to give accurate clustering on scales smaller than
$k=0.3\,h\,\rm{Mpc}^{-1}$ for the halo power spectrum, or
$k=0.2\,h\,\rm{Mpc}^{-1}$ for the matter power spectrum.

####  3.3.3 Newtonian motion gauge

The Newtonian motion gauge approach for massive neutrinos was developed by
Partmann et al. [32] and Heuschling et al. [49]. It allows for a simulation of
nonlinear CDM in an ordinary Newtonian $N$-body simulation while accounting
for the impact of linear neutrinos via a modification of the dark matter
initial conditions and by employing a dynamically evolving coordinate system.
The method is agnostic towards the implementation of the $N$-body simulation
and for this paper we choose to employ GADGET-4 [48]. However, any method
solving Newtonian nonlinear gravity is compatible, even methods other than
$N$-body simulations. We would like to stress that our method is exact in the
weak-field limit of general relativity (see Fidler et al. [78]) and therefore
captures the full effect of linear neutrino perturbations on the nonlinear
matter clustering.

The Newtonian motion approach allows for a very simple inclusion of massive
neutrinos, requiring only three additional steps applied to a simulation
without any neutrinos. First, we start from a set of “back-scaled” initial
conditions based on the present-day power spectrum of CDM and baryons in the
Newtonian motion coordinates, excluding neutrino perturbations. In contrast to
other neutrino methods presented in this paper, these initial conditions do
not assume a scale-dependent growth, i.e. the rescaling of the power spectrum
is done using the scale-independent growth factor $D_{+}$. The residual effect
of decaying modes due to neutrinos is included in the construction of the
present-day matter power spectrum in the Newtonian motion coordinate system.
This also has the added benefit that ordinary generators of initial conditions
can be used without modifications, provided that the correctly back-scaled
Newtonian motion gauge power spectrum is used. We then evolve the initial data
with the Newtonian solver, taking into account the impact of the massive
neutrinos on the background evolution via the Hubble rate. Finally, we obtain
the output in the Newtonian motion gauge. To make it comparable to the output
of other methods, we need to transform the result to the gauge employed
therein (usually the “$N$-boisson” gauge [78]). This step is realised by a
displacement field acting on the particle positions that is implemented in a
similar way to how the initial conditions are set. The transformation accounts
for the residual impact of neutrinos and other relativistic effects on the
evolution of the CDM and baryon particles. However, by construction, the
transformation vanishes exactly at $z=0$ (or another chosen target redshift)
while it is in general very small at late times, for small neutrino masses and
on small scales. Therefore, it can often be neglected as it will only lead to
small corrections in which case the output of the simulation can be used as-
is.

For this work, we include only the first two steps, while omitting the final
particle displacement. This leads to a small mismatch in the results shown for
$z=1$ at large scales for the cases of the highest neutrino masses. By leaving
this correction out, we demonstrate that for most neutrino masses, box sizes,
and redshifts the method is already sufficiently accurate in its simplest
form. For more details on the transformation we refer the reader to the
original work by Partmann et al. [32].

###  3.4 Halo-model reaction

The halo model reaction approach provides the nonlinear corrections caused by
massive neutrinos to a $\Lambda$CDM power spectrum through a ratio of halo-
model predictions. Following Cataneo et al. [79], the nonlinear power spectrum
is then given by

$P^{\rm NL}(k,z)=\mathcal{R}(k,z)P_{\rm pseudo}^{\rm NL}(k,z)\,,$ (3.6)

where $\mathcal{R}(k,z)$ is the halo-model reaction and $P_{\rm pseudo}^{\rm
NL}(k,z)$ is the nonlinear pseudo power spectrum.

The pseudo spectrum is a nonlinear $\Lambda$CDM power spectrum but with the
initial conditions tuned such that its linear clustering exactly matches the
linear clustering in the non-standard cosmology at the target redshift. For
example, if the non-standard physics introduces a simple rescaling of the
linear clustering amplitude, one could just rescale the amplitude of any
$\Lambda$CDM power spectrum to produce the pseudo spectrum. In the case of
scale-dependent modifications, this becomes a bit trickier in practice. We
approximate this quantity as in previous works by Cataneo et al. [79, 63] and
pass the modified linear spectrum as produced by CAMB to HMcode developed by
Mead et al. [80], with $\Lambda$CDM presets, i.e. no baryonic feedback nor
massive neutrinos. The benefit of using the pseudo rather than $\Lambda$CDM
cosmology is that it guarantees the mass functions in target and pseudo
cosmologies are similar as they have the same linear clustering. This produces
a smoother transition between the two-halo and one-halo terms. This transition
was one of the previous issues in calculating this nonlinear response using
the halo model [81, 82, 83].

Following Cataneo et al. [63], the halo-model reaction for massive neutrinos
is given by

$\mathcal{R}(k)=\frac{\left(1-f_{\nu}\right)^{2}P^{\mathrm{HM}}_{\text{cb}}(k)+2f_{\nu}\left(1-f_{\nu}\right)P^{\mathrm{HM}}_{\nu,\text{cb}}(k)+f_{\nu}^{2}P^{\mathrm{L}}_{\nu}(k)}{P^{\mathrm{L}}_{\mathrm{m}}(k)+P^{\mathrm{1h}}_{\mathrm{pseudo}}(k)}\,,$
(3.7)

with “cb” denoting the CDM and baryon component and “$\nu$” denoting massive
neutrinos. We include the effects of massive neutrinos at the linear level in
the numerator via the weighted sum of the nonlinear halo-model (cb) spectrum
and the massive neutrino linear spectrum [15]. The components of the reaction
are

$P^{\mathrm{HM}}_{\nu,\text{cb}}(k)\approx\sqrt{P^{\mathrm{HM}}_{\text{cb}}(k)P^{\mathrm{L}}_{\nu}(k)}\,,$
(3.8)
$P^{\mathrm{HM}}_{\text{cb}}(k)=P^{\mathrm{L}}_{\text{cb}}(k)+P^{\mathrm{1h}}_{\text{cb}}(k)\,.$
(3.9)

Explicitly, the one-halo terms are given as integrals over the Fourier space
halo density profile $u(k,M)$ and the halo mass function $n(M)$,

$P^{1\mathrm{h}}_{\text{cb}}(k)=\int\mathrm{d}\ln
M\,n_{\mathrm{cb}}(M)\left(\frac{M}{\bar{\rho}_{\mathrm{cb}}}\right)^{2}|u_{\mathrm{cb}}(k,M)|^{2}\,,$
(3.10) $P^{1\mathrm{h}}_{\mathrm{pseudo}}(k)=\int\mathrm{d}\ln
M\,n_{\mathrm{pseudo}}(M)\left(\frac{M}{\bar{\rho}_{\mathrm{m}}}\right)^{2}|u_{\mathrm{pseudo}}(k,M)|^{2}\,,$
(3.11)

where $\bar{\rho}$ is the background density for the specified matter species.
The halo mass functions are given as

$n_{\text{cb}}(M)=\frac{\bar{\rho}_{\text{cb}}}{M}[\nu^{\prime}f(\nu^{\prime})]\frac{\mathrm{d}\ln\nu^{\prime}}{\mathrm{d}\ln
M}\,,$ (3.12)
$n_{\mathrm{pseudo}}(M)=\frac{\bar{\rho}_{\mathrm{m}}}{M}[\nu^{\prime\prime}f(\nu^{\prime\prime})]\frac{\mathrm{d}\ln\nu^{\prime\prime}}{\mathrm{d}\ln
M}\,.$ (3.13)

The peak heights are defined as
$\nu^{\prime}=\delta_{\mathrm{sc,cb}}(M)/\sigma_{\mathrm{cb}}[R_{\mathrm{cb}}(M)]$
and
$\nu^{\prime\prime}=\delta_{\mathrm{sc,m}}/\sigma_{\mathrm{m}}[R_{\mathrm{m}}(M)]$,
where the subscript “sc” indicates this quantity is calculated by solving the
standard $\Lambda$CDM spherical-collapse equations but with the indicated
matter density. The mass fluctuation variances are given by

$[\sigma_{\mathrm{cb}}(R)]^{2}=\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}|\tilde{W}(kR)|^{2}P^{\mathrm{L}}_{\text{cb}}(k),$
(3.14)
$[\sigma_{\mathrm{m}}(R)]^{2}=\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}|\tilde{W}(kR)|^{2}P^{\mathrm{L}}_{\rm
m}(k).$ (3.15)

In all predictions from halo-model reaction, we employ a Sheth–Tormen halo
mass function [84, 85], a power-law concentration–mass relation (see for
example the work by Bullock et al. [86]), and an NFW halo density profile
[77]. The predictions are computed numerically using the public code ReACT
666https://github.com/nebblu/ReACT by Bose et al. [87, 64].

###  3.5 Power-spectrum emulation

Fast predictors of the matter power spectrum are an essential ingredient for
many inference pipelines in cosmology. Since numerical simulations are too
costly to be directly applied in this context, different approaches based on
approximate methods or elaborate fitting techniques have been used in the
past. Apart from the halo-model reaction discussed in the previous section,
well-known examples are the halofit predictor developed by Smith et al. [88]
and later improved by Takahashi et al. [89], and the HMcode predictor by Mead
et al. [80], the latter being based on the halo model. In terms of the power-
suppression signal of neutrinos, fitting routines by Bird et al. [16] have
been used in the past.

More recently, the emulation technique has become a popular alternative to
obtain fast predictions of the matter power spectrum within the cosmological
parameter space. Broadly speaking, emulators are regression routines based on
a suite of numerical simulations that sample the cosmological parameter space
and act as a training set. There are different surrogate techniques currently
used for cosmological emulators, such as Gaussian process regression [90, 65,
91, 92], polynomial chaos expansion [7, 8], or neural network approaches [93].

In addition to the predictions from halofit and HMcode mentioned above, we
focus in this paper on the Cosmic Emu [65], the EuclidEmulator2 [8], and the
BACCOemulator [45]. These emulators provide predictions of the matter power
spectrum and include a free parameter for the sum of the neutrino masses. We
will now summarise the particularities of these three emulators, specifically
focusing on the neutrino implementation.

  * •

The Cosmic Emu-2022 is built upon the Mira-Titan simulations [94, 66], a suite
of 111 simulations run with the HACC code [95]. The simulations are
distributed over an eight-dimensional cosmological parameter space comprising
($\omega_{\mathrm{m}}$, $\omega_{\mathrm{b}}$, $\omega_{\nu}$, $\sigma_{8}$,
$h$, $n_{\mathrm{s}}$, $w_{0}$, $w_{a}$), where $\sigma_{8}$ is the present-
day amplitude of linear matter density fluctuations at the scale of
$8\,${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc}$$, $n_{s}$ is the
scalar spectral index, and $w_{0}$, $w_{a}$ parameterise the effective
equation of state of dark energy in terms of the first two coefficients of the
Taylor series expansion around $a=1$. The emulator achieves an absolute
precision of about four percent for modes of $k<5\,h\,\mathrm{Mpc}^{-1}$
within the redshift range $z\in[0,2]$. Neutrinos are not incorporated in the
simulations and are effectively treated as a smooth background component. The
power spectra are then corrected on large scales for enhanced growth beyond
the neutrino free-streaming scale using the scale-dependent linear growth
factor.

  * •

The EuclidEmulator2 is trained on 200 paired and fixed simulations that were
run with PKDGRAV3. It emulates the nonlinear boost factor that is then
multiplied by the results of a linear Boltzmann solver. The emulator covers
eight cosmological parameters ($\Omega_{\mathrm{m}}$, $\Omega_{\mathrm{b}}$,
$\sum m_{\nu}$, $A_{\mathrm{s}}$, $h$, $n_{\mathrm{s}}$, $w_{0}$, $w_{a}$),
where $A_{s}$ is the amplitude of primordial perturbations at the scale
$k_{\text{p}}=0.05\,\mathrm{Mpc}^{-1}$, and includes redshifts of $z\in[0,3]$
and modes up to $k=10\,h\,\mathrm{Mpc}^{-1}$. It claims an error of below one
percent which is better than the other emulators discussed here. Note,
however, that the EuclidEmulator2 covers a somewhat narrower parameter space
motivated by the results of the Planck mission [96]. Within the training set
the neutrinos are modelled using the mesh-based method implemented in
PKDGRAV3.

  * •

The BACCOemulator is trained on a very large suite of simulations based on the
cosmology-rescaling technique [97]. More specifically, four high-resolution
simulations with judiciously chosen cosmologies are rescaled to more than 800
cosmologies at different redshifts. Whenever the target cosmologies included
massive neutrinos, their effect is added following the extension of the
cosmology-rescaling technique presented Zennaro et al. [33]. This emulator
varies eight cosmological parameters ($\Omega_{\text{cb}}$,
$\Omega_{\mathrm{b}}$, $\sum m_{\nu}$, $\sigma_{8}$, $h$, $n_{\mathrm{s}}$,
$w_{0}$, $w_{a}$) and covers a redshift range of $z\in[0,1.5]$ for modes up to
$k=5\,h\,\mathrm{Mpc}^{-1}$. The claimed precision is better than three
percent.

We refer to the original references for more information about the emulators.

##  4 Simulations

To compare different numerical methods, we carry out a large suite of $N$-body
simulations where we employ different codes to run the same set of ten
simulations summarised in Table 2. These ten simulations cover different
choices of total neutrino mass $\sum m_{\nu}$ (including the massless case),
different box sizes $L_{\mathrm{box}}$, and different mass resolutions to
check for numerical convergence with respect to finite-volume and
discretisation effects. $N_{\mathrm{part}}$ denotes the number of particles
used for CDM and baryons, as well as the number of particles for neutrinos if
a particle-based method is employed. For simplicity, we assume degenerate
neutrino mass eigenstates because cosmology is mainly sensitive to the total
neutrino mass scale [98]. We keep the total matter density (at redshift zero)
fixed at $\Omega_{\mathrm{m}}=0.319$ by adjusting the CDM density parameter
together with the neutrino mass. The baryon density is fixed at
$\Omega_{\mathrm{b}}=0.049$, and the remaining cosmological parameters are
$A_{\mathrm{s}}=2.215\times 10^{-9}$ at the pivot scale
$k_{\mathrm{p}}=0.05\,\mathrm{Mpc}^{-1}$, $n_{\mathrm{s}}=0.9619$, and
$h=0.67$, which are based on the Euclid Flagship 2 simulation. Dark energy is
modelled as a cosmological constant that provides a spatially flat background
solution, and the CMB temperature is set to $2.7255\,\rm{K}$ and the effect of
radiation is taken into account in the simulations at the linear level, either
by using carefully taylored initial conditions as detailed below or by
including the radiation component on the mesh if a mesh-based method is used.

Using identical initial data (see Section 4.1 below for details) in each case,
the ten simulations are run with each of the thirteen $N$-body methods listed
in Table 1 to produce particle snapshots at redshifts $z=1$ and $z=0$. For
AREPO, due to resource constraints, only the four simulations with
$N_{\mathrm{part}}=512^{3}$ are run, precluding the possibility of conducting
numerical convergence tests in this case. Therefore, a total of $248$
individual snapshots are analysed in this code comparison. In addition,
nonlinear power spectra are predicted for each distinct choice of neutrino
masses using the remaining methods in Table 1, as well as using the HMcode and
halofit fitting methods.

Table 2: Overview of the basic parameters used in our simulation suite. The
cases with $\sum m_{\nu}=0\,\mathrm{eV}$ and $\sum m_{\nu}=0.15\,\mathrm{eV}$
are the two main baselines for our comparison, but we include some cases with
larger masses, up to $\sum m_{\nu}=0.6\,\mathrm{eV}$, to probe more “extreme”
regions of parameter space.
Simulation | $L_{\mathrm{box}}$ | $N_{\mathrm{part}}$ | mass resolution | $\sum m_{\nu}$
---|---|---|---|---
0.0eV | 512 ${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc}$ | 5123 | 8.85$\times$10${}^{10}~{}h^{-1}$\mathrm{\text{${M}_{\odot}$}}$$ | 0.0 eV
0.15eV | 512 ${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc}$ | 5123 | 8.75$\times$10${}^{10}~{}h^{-1}$\mathrm{\text{${M}_{\odot}$}}$$ | 0.15 eV
0.3eV | 512 ${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc}$ | 5123 | 8.65$\times$10${}^{10}~{}h^{-1}$\mathrm{\text{${M}_{\odot}$}}$$ | 0.3 eV
0.6eV | 512 ${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc}$ | 5123 | 8.45$\times$10${}^{10}~{}h^{-1}$\mathrm{\text{${M}_{\odot}$}}$$ | 0.6 eV
0.0eV_HR | 512 ${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc}$ | 10243 | 1.11$\times$10${}^{10}~{}h^{-1}$\mathrm{\text{${M}_{\odot}$}}$$ | 0.0 eV
0.15eV_HR | 512 ${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc}$ | 10243 | 1.09$\times$10${}^{10}~{}h^{-1}$\mathrm{\text{${M}_{\odot}$}}$$ | 0.15 eV
0.0eV_1024Mpc | 1024 ${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc}$ | 10243 | 8.85$\times$10${}^{10}~{}h^{-1}$\mathrm{\text{${M}_{\odot}$}}$$ | 0.0 eV
0.15eV_1024Mpc | 1024 ${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc}$ | 10243 | 8.75$\times$10${}^{10}~{}h^{-1}$\mathrm{\text{${M}_{\odot}$}}$$ | 0.15 eV
0.3eV_1024Mpc | 1024 ${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc}$ | 10243 | 8.65$\times$10${}^{10}~{}h^{-1}$\mathrm{\text{${M}_{\odot}$}}$$ | 0.3 eV
0.6eV_1024Mpc | 1024 ${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc}$ | 10243 | 8.45$\times$10${}^{10}~{}h^{-1}$\mathrm{\text{${M}_{\odot}$}}$$ | 0.6 eV

###  4.1 Initial conditions

The initial conditions of all simulations are generated at redshift
$z=127$.777COLA exceptionally uses $z=19$ as initial redshift (see Section
3.3.1). The linear matter power spectra and transfer functions are obtained by
running either CAMB or the CLASS Boltzmann code by Blas et al. [99]. These
files are then used by the REPS888https://github.com/matteozennaro/reps code
to compute the rescaled power spectra and transfer functions at $z=127$ by
solving the multi-fluid linear equations as outlined by Zennaro et al. [100].
This procedure, known as rescaling, guarantees that the power spectrum of the
output of the simulation on linear scales at low redshift will match the
correct linear power spectra. A realisation of initial data is then generated
by drawing random phases for all perturbation modes and fixing their
amplitudes according to the initial transfer functions. This approach of
“fixing” the amplitudes effectively removes cosmic variance at linear scales
and has been shown to generally produce less noisy summary statistics [101].
It introduces a specific type of non-Gaussianity that is not expected to
affect any of our results. Given the density field of CDM and baryons, the
initial positions and velocities of the $N$-body particles are computed from
the Zeldovich approximation [102]. We employ a modified version of the N-GenIC
code999https://github.com/franciscovillaescusa/N-GenIC_growth that accounts
for the scale-dependence present in both the growth rate and growth factor in
cosmologies with massive neutrinos.

For neutrinos, the different implementations make use of distinct methods. For
particle-based implementations, the positions and velocities of the massive
neutrino particles are generated in a similar fashion as CDM. This means we
effectively make use of the Zeldovich approximation to set the first two
moments of the phase-space distribution function which correspond to the
density perturbation and bulk velocity, respectively. The initial velocities
of neutrino particles are then offset by a random thermal component drawn from
their Fermi–Dirac distribution [103, 104]. We assume the standard Big Bang
scenario where this distribution is set in equilibrium before the weak
interaction freezes out — when the Universe was about 1 second old — and after
freeze-out simply redshifts as the Universe expands. Note that this means that
the typical thermal velocities are much larger than those of a distribution
that is in thermal equilibrium at low redshift. For SWIFT, neutrino particles
are instead set up with the FastDF code,101010https://github.com/wullm/fastdf
using geodesic integration from high redshift [68, 20], to reproduce the full
distribution function and to prevent the initial perturbations from being
erased by thermal motions [55]. For mesh-based implementations, on the other
hand, the density field of each neutrino species is directly computed using
the phases from the random field realisation of the linear initial conditions.

###  4.2 Post-processing pipeline

In order to quantify differences in our numerical schemes as precisely as
possible, we analyse the snapshots of all our $N$-body simulations in a common
pipeline. We compute the power spectra and bispectra of the CDM and baryon
component, and produce halo catalogues from which we measure the halo mass
functions and halo bias. In cases where the simulations provide a neutrino
distribution, we also compute the cross-power spectra of neutrinos with the
CDM and baryon component, as well as the neutrino auto-power spectra.

####  4.2.1 Power spectra

The power spectra of the different snapshots have been estimated using
Pylians3.111111https://pylians3.readthedocs.io The routine first deposits
particle masses into a regular 3D grid with $N^{3}$ voxels using the cloud-in-
cell mass-assignment scheme. In this work, we always use a mesh with
$N^{3}=N_{\mathrm{part}}$ such that the Nyquist scales match between particles
and mesh. The constructed field is then Fourier transformed and the effects of
the mass-assignment scheme are corrected. Next, for each mode the square of
its amplitude is computed, $|\delta(\bm{k})|^{2}$. The modes are then binned
in intervals of width equal to the fundamental frequency $k_{\mathrm{F}}=2\pi
L_{\mathrm{box}}^{-1}$ and the power spectrum is finally estimated as

$P(k_{i})=\frac{1}{N_{i}}\sum_{\bm{k}\in B_{i}}|\delta(\bm{k})|^{2}\,,$ (4.1)

where $N_{i}$ is the number of independent modes in the considered bin
$B_{i}=\\{\bm{k}\,|\,i\,k_{\mathrm{F}}\leq|\bm{k}|<(i+1)k_{\mathrm{F}}\\}$ and
$k_{i}$ is computed as

$k_{i}=\frac{1}{N_{i}}\sum_{\bm{k}\in B_{i}}|\bm{k}|\,.$ (4.2)

To compute the cross-power spectrum of two fields instead, the estimator is
generalised in the most straightforward way,

$P_{\mathrm{X},\mathrm{Y}}(k_{i})=\frac{1}{2N_{i}}\sum_{\bm{k}\in
B_{i}}\Bigl{[}\delta_{\mathrm{X}}(\bm{k})\delta_{\mathrm{Y}}^{\ast}(\bm{k})+\delta_{\mathrm{X}}^{\ast}(\bm{k})\delta_{\mathrm{Y}}(\bm{k})\Bigr{]}\,,$
(4.3)

where the subscripts “X” and “Y” denote the two different fields. When
presenting the results, we combine the measurements into larger bins
logarithmically spaced in $k$. This reduces the noise at large $k$ and makes
our plots more readable.

Some codes do not produce snapshots at exactly the desired redshift, but at
redshifts that deviate by less than $\pm\,0.01$ from the target redshift.
These differences can be visible when comparing power spectra on a sub-percent
level. For those cases, we rescale the power spectra by the square of the
ratio of the linear growth factors at the respective redshift values. Such a
rescaling is applied to ANUBIS for all snapshots and to L-GADGET3 and PKDGRAV3
at $z=1$.

####  4.2.2 Bispectra

We measure the bispectrum of CDM and baryons (ccc) using the estimator

$\hat{B}_{\mathrm{ccc}}(k_{l},k_{m},k_{n})\equiv\frac{k_{\mathrm{F}}^{3}}{N_{\mathrm{tr}}(k_{l},k_{m},k_{n})}\sum_{\bm{q}_{1}\in
B_{l}}\sum_{\bm{q}_{2}\in B_{m}}\sum_{\bm{q}_{3}\in
B_{n}}\,\delta_{\mathrm{K}}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})\,\,\delta(\bm{q}_{1})\,\delta(\bm{q}_{2})\,\delta(\bm{q}_{3}),$
(4.4)

where $N_{\mathrm{tr}}$ is the number of “fundamental triangles”,

$N_{\mathrm{tr}}(k_{l},k_{m},k_{n})\equiv\sum_{\bm{q}_{1}\in
B_{l}}\sum_{\bm{q}_{2}\in B_{m}}\sum_{\bm{q}_{3}\in
B_{n}}\,\delta_{\mathrm{K}}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})\,,$ (4.5)

formed by the vectors $\bm{q}_{i}$ satisfying the triangle condition
$\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3}=0$ that are included within the “triangle
bin” defined by the triplet of centers $(k_{l},k_{m},k_{n})$ and corresponding
bins $B_{l},B_{m},B_{n}$.

We use a Python code implementing the fourth-order density interpolation and
the interlacing scheme described by Sefusatti et al. [105]. In order to
compare the large-volume simulations
($L_{\mathrm{box}}=1024\,h^{-1}\,\mathrm{Mpc}$) more easily with the small-
volume ones, we use the same $k$-space binning in both cases, fixing the bin
width to $k_{\mathrm{F}}$ of the small box. Just like for the power spectra,
to account for inaccuracies in the redshift of some snapshots, we rescale some
of the resulting bispectra by the cube of the ratio of the linear growth
factors at the respective redshift values. Such a rescaling is applied to
ANUBIS for all snapshots and to L-GADGET3 and PKDGRAV3 at $z=1$.

####  4.2.3 Halo catalogues

For the considered snapshots of the various simulations, we identify halos
with the code Denhf [106, 107, 108, 109] which uses a “spherical overdensity”
criterion. The algorithm does not rely on any pre-identification method. Only
CDM and baryon particles are considered in the characterisation of halos;
neutrino particles (if present at all in the simulation) are considered as a
background component [17, 76].

Denhf estimates the local density at the position of each $N$-body particle by
calculating the distance to its 10th-nearest neighbour $d_{10}$, and assigning
to each particle a density that is proportional to $d_{10}^{-3}$. Centered on
the particle with the highest density value, the algorithm grows a sphere and
stops when the mean density within the sphere falls below a desired
overdensity threshold, set to $200$ times the background density of CDM and
baryons for the purpose of this work. All particles assigned to this spherical
overdensity halo are then removed from the global list of particles, and the
algorithm proceeds recursively until none of the remaining particles has a
local density large enough to be the center of a 10-particle halo. Particles
not assigned to halos will be part of the field.

Only in the case of PINOCCHIO we use the halo catalogues as produced by the
code itself instead of Denhf. Because PINOCCHIO is calibrated on the friends-
of-friends halo mass function, we translate its masses to spherical
overdensity ones by measuring friends-of-friends and spherical overdensity
masses for halos of the set of GADGET-3 runs, and relating them through
abundance matching.

For the codes that do not produce snapshots at the exact values of the desired
redshifts, we apply no further corrections here. The error in the redshift is
less than $\pm\,0.01$ while our halo properties typically display
disagreements larger than 1% between different codes. We therefore assume that
the error due to mismatching redshift values is subdominant.

##  5 Results

###  5.1 Power spectra

A key prediction from neutrino simulations is a suppression of the matter
power spectrum that exceeds the maximum linear theory prediction of $\Delta
P_{\text{m}}/P_{\text{m}}\approx 8f_{\nu}$. The matter power spectrum can be
decomposed as follows

$P_{\text{m}}(k)=(1-f_{\nu})^{2}P_{\text{cb}}(k)+2f_{\nu}(1-f_{\nu})P_{\nu,\text{cb}}(k)+f_{\nu}^{2}P_{\nu}(k)\,,$
(5.1)

where $P_{\text{cb}}$ is the power spectrum of CDM and baryons,
$P_{\nu,\text{cb}}$ the cross-power spectrum of neutrinos with CDM and
baryons, and $P_{\nu}$ the neutrino power spectrum. Various methods treat
these components differently or make predictions for only some of them, so we
discuss each component in turn. Finally, we will also compare the total
$P_{\text{m}}(k)$ with various power spectrum emulators, including
EuclidEmulator2 which only predicts this quantity.

####  5.1.1 CDM and baryons

Figure 1: Power spectrum of CDM and baryons as measured from different codes,
relative to GADGET-3 at $z=0$ and $z=1$ for a neutrino mass of $\sum
m_{\nu}=0.15\rm{eV}$ in the simulation with
$L_{\mathrm{box}}=512\,h^{-1}\,\mathrm{Mpc}$ and $N_{\mathrm{part}}=512^{3}$
CDM and baryon particles. The corresponding particle Nyquist wavenumber is
indicated by the grey dash-dotted line. The grey bands highlight the interval
of $\pm\,0.01$.

The leading contribution to $P_{\text{m}}(k)$ is $P_{\text{cb}}(k)$, which is
suppressed in massive neutrino models. Figure 1 shows the ratio of
$P_{\text{cb}}(k)$ for models with a neutrino mass of $\sum
m_{\nu}=0.15\,\rm{eV}$ relative to the GADGET-3 simulation which we
arbitrarily pick as the reference. In all our figures, results from codes
where massive neutrinos are represented through an $N$-body ensemble are
plotted using solid lines, while other $N$-body methods, including surrogates,
use dashed lines. Any additional predictions use dotted lines. In Figure 1,
the linear prediction, computed with CLASS, is shown by the rose dotted line
that drops off sharply at $k\approx 0.1\,h\,\mathrm{Mpc}^{-1}$ beyond which
the error quickly exceeds $10\%$. At $z=0$, on the largest scales, all codes
deviate less than 1% — the fluctuations seen in the dotted lines are largely
due to the lack of cosmic variance in the codes that predict
$P_{\mathrm{cb}}(k)$ directly. On smaller scales, some of the codes start to
deviate from the GADGET-3 reference. PINOCCHIO is in agreement within 1% up to
$k\approx 0.1\,h\,\rm{Mpc}^{-1}$, ANUBIS up to $k\approx
0.3\,h\,\rm{Mpc}^{-1}$, gevolution up to $k\approx 0.5\,h\,\rm{Mpc}^{-1}$, and
COLA up to $k\approx 0.7\,h\,\rm{Mpc}^{-1}$. The other codes stay within a 1%
deviation from GADGET-3 for all scales down to the particle Nyquist scale of
$k_{\text{Nyq}}=\pi\,N_{\mathrm{part}}^{1/3}\,L_{\mathrm{box}}^{-1}\approx
3\,h\,\rm{Mpc}^{-1}$. This scale is indicated by a vertical dash-dotted line.
The ReACT and BACCOemulator codes stay accurate within 1% to 5% on all scales
down to the Nyquist scale. At $z=1$ the qualitative behaviour is similar, but
most notably the codes disagree more on large scales while still staying
within a 1% agreement. Here, PINOCCHIO is accurate up to $k\approx
0.2\,h\,\rm{Mpc}^{-1}$.

Figure 2: Suppression of the power spectrum of CDM and baryons at $z=0$ and
$z=1$ for three different neutrino masses, $\sum
m_{\nu}\in\\{0.15,\,0.3,\,0.6\\}\,\rm{eV}$, when compared to the massless
case. Results are from the simulations with
$L_{\mathrm{box}}=512\,h^{-1}\,\mathrm{Mpc}$ and $N_{\mathrm{part}}=512^{3}$.
The corresponding particle Nyquist wavenumber is indicated by the grey dash-
dotted line.

Figure 2 shows the ratio of $P_{\text{cb}}(k)$ for models with massive
neutrinos relative to the massless case for $\sum
m_{\nu}\in\\{0.15,\,0.3,\,0.6\\}\,\rm{eV}$. For later convenience, we denote
this suppression ratio by $S_{\text{cb}}(k)$, or in general

$S_{\mathrm{X}}(k)=\frac{P_{\mathrm{X}}^{\mathrm{(massive)}}(k)}{P_{\mathrm{X}}^{\mathrm{(massless)}}(k)}\,,$
(5.2)

where the subscript “X” denotes any component in question, here
$\mathrm{X}\to\text{cb}$. Linear calculations (taken from CLASS) predict that,
on the largest scales, the growth of structure is mostly unaffected so that
$S_{\text{cb}}(k)$ approaches unity, while on small scales $S_{\text{cb}}(k)$
reaches a plateau. Compared to this linear expectation, all codes in the
comparison predict slightly less suppression around $k=0.1\,h\,\rm{Mpc}^{-1}$
and a much greater suppression for $k>0.3\,h\,\rm{Mpc}^{-1}$, followed by an
upturn on nonlinear scales. This upturn has been repeatedly demonstrated and
results from the reduced sensitivity of the one-halo contribution [110, 41].
At $z=0$, we obtain excellent agreement between all simulations and most
approximate methods up to the scale of maximum suppression at
$k_{\rm{max}}\approx 1\,h\,\rm{Mpc}^{-1}$, where the suppression is $20\%$
greater than the linear prediction. At $z=1$, the scale of maximum suppression
shifts to $k_{\rm{max}}\approx 2\,h\,\rm{Mpc}^{-1}$ and the differences are
greater, both with the linear prediction and between the codes, with the
exception of PINOCCHIO 121212The reason for this disagreement is due to the
fact that the accuracy of LPT-based PINOCCHIO depends on the level of
nonlinearity that varies with neutrino mass. which fares significantly better
compared to $z=0$.

Figure 3: The suppression of the power spectrum of CDM and baryons relative to
the one measured in the GADGET-3 reference runs at $z=0$ and $z=1$ for $\sum
m_{\nu}=0.15\,\rm{eV}$ for the simulations with
$L_{\mathrm{box}}=512\,h^{-1}\,\mathrm{Mpc}$ and $N_{\mathrm{part}}=512^{3}$.
The corresponding particle Nyquist wavenumber is indicated by the grey dash-
dotted line. The grey bands highlight the interval of $\pm\,0.01$.

To study these relative differences in greater detail, we show
$S_{\text{cb}}(k)$ for the smallest neutrino mass, $\sum
m_{\nu}=0.15\,\rm{eV}$, relative to the GADGET-3 prediction for this quantity
in Figure 3. With the exception of PINOCCHIO and CLASS, all codes agree to
better than $1\%$ at $z=0$ up to the particle Nyquist scale. Near
$k_{\text{Nyq}}$, the approximate COLA method and the BACCOemulator differ by
more than $1\%$ from the bulk of the simulations, while gevolution differs by
slightly less than 1%. Beyond this scale, the predictions diverge and should
be compared to higher-resolution runs since our estimator of the power
spectrum is computed on a mesh with a matching Nyquist scale. The measured
power spectra therefore cannot be used beyond $k_{\text{Nyq}}$ where they
become strongly biased. This particularly affects the comparison between
simulations — where the power spectrum estimator is employed — and other
methods to predict $P_{\mathrm{cb}}(k)$. At $z=1$, nonlinearities are smaller
and the agreement between the simulations is better. Here the snapshot
produced by PINOCCHIO achieves percent accuracy to $k\approx
0.3\,h\,\rm{Mpc}^{-1}$. Some of the approximate methods fare slightly worse at
this earlier time compared to $z=0$, with the difference between CLASS and the
simulations increasing by $50\%$ and COLA diverging beyond
$k=2\,h\,\rm{Mpc}^{-1}$. The NM-GADGET4 method requires an additional post-
processing coordinate transformation at any redshift except $z=0$. Because
this additional step is omitted in this work for simplicity, a small error at
low wavenumbers remains in the $z=1$ data. This explains why the error is
larger at that redshift than at the final time.

ANUBIS notably drops off for $k>1\,h\,\rm{Mpc}^{-1}$ at $z=1$ and also shows a
small excess on linear scales both at $z=0$ and $z=1$. This excess originates
from the massless case and is a result of ANUBIS being run with a coarser base
grid than the other codes ($512^{3}$ for all simulations) due to limited
resources. A finer base grid requires more memory but the lack of it can be
somewhat compensated by using a smaller time step. For the ANUBIS massive
neutrino runs, this is done automatically but for the massless case the time
step has been set to half of that originally calculated by the code. Tests
indicate that further reducing the time step or ideally using a finer base
grid should lessen the excess at large scales, but finding the optimal choice
of code settings is not the aim of this work. The drop-off observed for ANUBIS
at $z=1$ for $k>1\,h\,\rm{Mpc}^{-1}$ is due to differences in resolution
between the various simulations. As ANUBIS is an AMR-code, a modified version
of the RAMSES code originally written by Teyssier [56], the inclusion of
massive neutrinos, which suppresses clustering on small scales, also reduces
the number of refinements reached in the simulation compared to the massless
case. This effect increases with the neutrino mass and is more noticeable for
higher redshifts where there is also less refinement due to less clustering.
This can be solved by a higher particle density which automatically leads to
more refinement. Generally this is also necessary to find a better agreement
between ANUBIS and GADGET-3 for smaller scales, as can be observed from the
high-resolution runs and also in a previous comparison conducted by Schneider
et al. [111] between RAMSES, GADGET-3, and PKDGRAV3.

Figure 4: The suppression of the power spectrum of CDM and baryons relative to
the one measured in the GADGET-3 reference runs at $z=0$ for $\sum
m_{\nu}\in\\{0.3,0.6\\}\,\rm{eV}$ for the simulations with
$L_{\mathrm{box}}=512\,h^{-1}\,\mathrm{Mpc}$ and $N_{\mathrm{part}}=512^{3}$.
The corresponding particle Nyquist wavenumber is indicated by the grey dash-
dotted line. The grey bands highlight the interval of $\pm\,0.01$.

The observed differences are greater when the neutrino mass is increased,
especially for the approximate methods. Figure 4 shows $S_{\text{cb}}(k)$
relative to the GADGET-3 prediction for models with $\sum
m_{\nu}=0.3\,\rm{eV}$ (left panel) and $\sum m_{\nu}=0.6\,\rm{eV}$ (right
panel) at $z=0$. While the BACCOemulator agrees to better than $1\%$ with
GADGET-3 up to $k_{\text{Nyq}}$ for $0.3\,\rm{eV}$, it makes no prediction for
$0.6\,\rm{eV}$ because that value lies far outside the range covered by the
training set of the emulator. We deliberately include a case of such a large
neutrino mass to exacerbate the differences between the various numerical
implementations. The reaction method differs by slightly more than $1\%$ and
$2\%$ at $k=0.7\,h\,\rm{Mpc}^{-1}$ for total neutrino masses of $0.3\,\rm{eV}$
and $0.6\,\rm{eV}$, respectively. The differences with the linear prediction
(CLASS) and with PINOCCHIO are large, as noted above for $0.15\,\rm{eV}$.

The agreement is excellent for the other codes, but some subtle differences
can be discerned. On large scales, we observe that SWIFT, CONCEPT, COLA,
gevolution, PKDGRAV3, NM-GADGET4, and L-GADGET3 show the same coherent scatter
relative to GADGET-3, especially for a neutrino mass of $0.6\,\rm{eV}$. This
is due to the contamination from shot noise in the neutrino particle
implementation used by the GADGET-3 run. The other mentioned codes have
implementations that do not suffer from shot noise or take measures to limit
the contamination. For instance, in gevolution the neutrino $N$-body ensemble
is evolved throughout the simulation, but it is only used as source of
gravitational fields from redshift $z=7$ and below. At higher redshifts, the
code uses the linear grid-based density instead. The reasoning behind this
strategy is that shot noise is constant over time and hence more problematic
at high redshift where cosmological perturbations are smaller in comparison.
On the other hand, the linear prediction for neutrinos is expected to be very
accurate at high redshift. One can therefore reduce the total error by
judiciously choosing the time at which the code switches from linear to fully
nonlinear neutrino treatment. It is nonetheless reassuring that even without
mitigating against shot noise the scatter remains far below $1\%$.

On small scales the differences are also larger. The lines for COLA and
gevolution track each other closely, but differ by more than $1\%$ from the
other codes for $k>1\,h\,\rm{Mpc}^{-1}$. SWIFT, PKDGRAV3, and CONCEPT are low
compared to GADGET-3 for $k>1\,h\,\rm{Mpc}^{-1}$, unlike what was seen in
Figure 3. ANUBIS diverges from GADGET-3 by more than $1\%$ for
$k>1\,h\,\rm{Mpc}^{-1}$ for the $0.6\,\rm{eV}$ neutrino mass case. As
mentioned earlier, this is due to the fact that the AMR scheme has a lower
effective resolution as the neutrino mass increases, simply because refinement
is triggered by clustering. Finally, NM-GADGET4 is slightly higher than
GADGET-3. However, these differences remain below $1\%$ well beyond
$k_{\text{Nyq}}$.

####  5.1.2 Convergence tests

To study the numerical convergence of our results, we consider the effects of
finite box size and resolution. Figure 5 shows the relative suppression for
the runs with larger volume (left panel) and higher resolution (right panel).
When $L_{\mathrm{box}}$ is doubled at fixed resolution, the agreement remains
excellent on linear scales and is sometimes even slightly better around
$k_{\text{Nyq}}$, providing an important consistency check for most codes.
Increasing instead the mass resolution by doubling $k_{\text{Nyq}}$, the
agreement between the simulations improves significantly on nonlinear scales.
Including more scales in either direction, most codes remain within $1\%$ of
the reference runs done with GADGET-3. The excess on large scales for the case
of ANUBIS persists as a result of the coarse base grid. For the runs with
$N_{\mathrm{part}}=1024^{3}$, this base grid is even less suited and the time-
steppings for the massless cases are set to $0.1$ and $0.15$ times the
original time step for the larger-box and high-resolution runs, respectively.

Figure 5: The suppression of the power spectrum of CDM and baryons relative to
the one measured in the GADGET-3 reference runs at $z=0$ for $\sum
m_{\nu}=0.15\,\rm{eV}$ for the simulations with a larger volume,
$L_{\mathrm{box}}=1024\,h^{-1}\,\mathrm{Mpc}$ and $N_{\mathrm{part}}=1024^{3}$
(left panel), and at a higher resolution in the small volume,
$L_{\mathrm{box}}=512\,h^{-1}\,\mathrm{Mpc}$ and $N_{\mathrm{part}}=1024^{3}$
(right panel). The respective particle Nyquist wavenumbers are indicated by
the grey dash-dotted lines. The grey bands highlight the interval of
$\pm\,0.01$.

####  5.1.3 Contributions from neutrinos

The subdominant contributions to $P_{\text{m}}(k)$ are the cross-power
spectrum $P_{\nu,\text{cb}}(k)$ between neutrinos and the CDM and baryon
component, and the auto-power spectrum $P_{\nu}(k)$ of neutrinos. As can be
seen from Eq. (5.1), these are suppressed by the small factors $f_{\nu}$ and
$f_{\nu}^{2}$, respectively, and they are themselves additionally strongly
suppressed with respect to $P_{\text{cb}}(k)$ on scales smaller than the
neutrino free-streaming scale. While both of these contributions will be
exceedingly hard to constrain individually from observations, it is
nonetheless interesting to study them in the context of our code comparison in
order to highlight some more subtle differences in the numerical schemes.
Figure 6 (left panel) shows $P_{\nu,\text{cb}}(k)$ for various codes relative
to the result from the SWIFT code, for the smallest neutrino mass, $\sum
m_{\nu}=0.15\,\rm{eV}$, computed from the high-resolution simulations. We use
SWIFT as the reference here because it has a very low level of shot noise in
the neutrino component, yet is able to track the nonlinear evolution of
neutrinos very accurately. As is the case with CDM and baryons, linear theory
cannot describe neutrino clustering on nonlinear scales and therefore CLASS
significantly underestimates the cross-power spectrum for
$k>0.1\,h\,\rm{Mpc}^{-1}$. Most other codes use a particle implementation of
neutrinos and scatter about the SWIFT prediction partially due to shot noise.
The results from CONCEPT have no shot noise, but depart from the other codes
for $k>0.2\,h\,\rm{Mpc}^{-1}$. The relative difference to SWIFT is $8\%$ at
$k=1\,h\,\rm{Mpc}^{-1}$. Much the same applies to the neutrino auto-power
spectrum, $P_{\nu}(k)$, shown relative to the neutrino auto-power spectrum of
SWIFT in the right panel of Figure 6. Here, for the codes with a neutrino
particle ensemble, the dominant contribution to the shot noise is removed.
Except for SWIFT this is done by subtracting the inverse of the average
neutrino particle density, $\bar{n}^{-1}$, from the measured power spectrum,
where $\bar{n}^{-1}=0.125\,h^{-3}\,\mathrm{Mpc}^{3}$ in our high-resolution
simulations. For SWIFT the subtracted values are derived from the $\delta f$
method. The difference between SWIFT and CONCEPT is $19\%$ (off the chart) at
$k=1\,h\,\rm{Mpc}^{-1}$ and the effects of shot noise are even more evident in
the other codes. The agreement between gevolution and SWIFT is quite
remarkable, and on the largest scales these results are also more consistent
with CONCEPT than with the other $N$-body codes.

Figure 6: Relative cross-power spectrum of neutrinos with CDM and baryons
(left) and neutrino auto-power spectrum (right) with respect to SWIFT at $z=0$
for $\sum m_{\nu}=0.15\,\rm{eV}$ for the higher-resolution simulation with
$L_{\mathrm{box}}=512\,h^{-1}\,\mathrm{Mpc}$ and $N_{\mathrm{part}}=1024^{3}$.
The Nyquist wavenumber is indicated by the grey dash-dotted line. The grey
bands highlight the interval of $\pm\,0.01$.

####  5.1.4 Total matter

Some of the codes, in particular the EuclidEmulator2, only provide predictions
for the power spectrum of total matter, $P_{\mathrm{m}}(k)$. In Figure 7, we
show the relative agreement of different emulators and other codes predicting
$P_{\mathrm{m}}(k)$ to our high-resolution reference run with GADGET-3. The
left panel shows the result for the matter power spectrum itself,
$P_{\mathrm{m}}(k)$, at $\sum m_{\nu}=0.15$ eV. The correlated fluctuations on
large scales are due to sample variance which is only present in the reference
simulation and not in the predicted spectra. We note that emulators, HMcode,
and the halo-model reaction method perform slightly better than the fitting
recipe of halofit. The results of the EuclidEmulator2 are marginally
consistent with the claimed accuracy of $1\%$, but the neutrino mass lies
close to the boundary of parameter space the emulator was trained for. The
right panel of Figure 7 shows corresponding results for the power suppression
factor, $S_{\mathrm{m}}(k)$, with respect to the massless scenario. Cosmic Emu
shows a disagreement larger than $1\%$ around $k\approx
1\,h\,\mathrm{Mpc}^{-1}$ where the power suppression is largest. Also halofit
performs poorly, worse even than linear theory (CLASS). The other codes remain
within $1\%$ of the GADGET-3 result up to the particle Nyquist scale.

Figure 7: Total matter power spectrum $P_{\text{m}}(k)$ at $z=0$ for $\sum
m_{\nu}=0.15\,\rm{eV}$ from emulators and fitting methods (left panel), and
the respective suppression $S_{\text{m}}(k)$ with respect to the massless case
(right panel), compared to the higher-resolution reference GADGET-3 simulation
with $L_{\text{box}}=512\,h^{-1}\,\rm{Mpc}$ and $N_{\mathrm{part}}=1024^{3}$.
The corresponding particle Nyquist wavenumber of the GADGET-3 run is indicated
by the grey dash-dotted line and marks the limit beyond which the estimator of
the power spectrum becomes unreliable. Disagreements near and beyond this line
are therefore not indicative of errors in the emulators. The grey bands
highlight the interval of $\pm\,0.01$.

###  5.2 Bispectra

The nonlinear evolution of matter fluctuations generates a non-vanishing
bispectrum, the three-point correlation function of matter in Fourier space,
even if non-Gaussianity is negligibly small in the initial conditions. This
represents an opportunity for the measurement of neutrino masses, as we expect
that the suppression predicted by linear theory is enhanced at the nonlinear
level [112, 113, 114, 115]. The total matter bispectrum in the presence of
massive neutrinos can be schematically defined as

$B_{\mathrm{mmm}}(k_{1},k_{2},k_{3})=(1-f_{\nu})^{3}\,B_{\mathrm{ccc}}(k_{1},k_{2},k_{3})+f_{\nu}\,(1-f_{\nu})^{2}\,B^{(\mathrm{sym})}_{\mathrm{cc}\nu}(k_{1},k_{2},k_{3})\\\
+f_{\nu}^{2}\,(1-f_{\nu})\,B^{(\mathrm{sym})}_{\mathrm{c}\nu\nu}(k_{1},k_{2},k_{3})+f_{\nu}^{3}\,B_{\nu\nu\nu}(k_{1},k_{2},k_{3}),$
(5.3)

where “ccc” denotes the CDM and baryon bispectrum (we do not write “cb cb cb”
to avoid clutter) and we note the presence of cross cold-cold-neutrino
“$\mathrm{cc}\nu$” and cold-neutrino-neutrino “$\mathrm{c}\nu\nu$” terms that
are symmetrised, as indicated by the superscript “$(\mathrm{sym})$”. In this
work, we focus on the bispectrum of CDM and baryons only, but cross terms have
also been investigated in the literature, e.g. by Ruggeri et al. [116]. As for
the power spectrum case, the leading term is the one of CDM and baryons.

Figure 8: The coloured bands show the scatter of $B_{\mathrm{ccc}}$ for all
triangle configurations with respect to the measurement from GADGET-3 at
redshift $z=1$. The horizontal axis indicates the maximum wavenumber in each
triangle configuration, $k_{\mathrm{max}}=\max(k_{1},k_{2},k_{3})$, in order
to present the results in a simple plot. The four panels show different
neutrino masses. Each plot uses three subpanels to make the results of the
many different codes more visible. Figure 9: Squeezed bispectrum of CDM and
baryons as measured from different codes, relative to GADGET-3 at $z=0$ and
$z=1$ for a neutrino mass of $\sum m_{\nu}=0.15\rm{eV}$ in the simulation with
$L_{\mathrm{box}}=512\,h^{-1}\,\rm{Mpc}$ and $N_{\mathrm{part}}=512^{3}$ CDM
and baryon particles. Results shown are for the squeezed configuration where
$k_{1}=k_{2}\equiv k$, $k_{3}=0.012\,h\,\mathrm{Mpc}^{-1}$. The grey bands
highlight the interval of $\pm\,0.01$.

Figure 8 shows a comparison of the bispectrum $B_{\mathrm{ccc}}$ as measured
at redshift $z=1$ by all simulation codes131313PINOCCHIO is not used in this
test that is mostly focused on nonlinear scales beyond its range of validity.
for all triangle configurations and all neutrino masses considered, for the
set of runs with $N_{\mathrm{part}}=512^{3}$. Triangles are plotted as a
function of $k_{\mathrm{max}}=\max(k_{1},k_{2},k_{3})$. The comparison shows
the maximum percentage difference with respect to GADGET-3 of all triangle
configurations at each value of $k_{\mathrm{max}}$. For L-GADGET3,
openGADGET3, GADGET-4, PKDGRAV3, CONCEPT, and NM-GADGET4, discrepancies are
mostly within $5\%$ for $\sum m_{\nu}=0.0$ eV, while steadily growing up to
around $10\%$ for massive neutrino cosmologies. Other simulation codes are
within $10\%$ already at $\sum m_{\nu}=0.0$ eV. ANUBIS and gevolution show
some of the strongest deviations at large $k_{\mathrm{max}}$, but this is
mainly a result of finite resolution as we find much better agreement in the
higher-resolution runs. At low $k_{\mathrm{max}}$, strong fluctuations can be
observed where measurements can cross zero because of sampling variance, which
in turn leads to numerical issues when taking ratios.

In Figure 9, we consider specifically a squeezed configuration for which
$k_{1}=k_{2}$ and $k_{3}=0.012\,h\,\mathrm{Mpc}^{-1}$ and show the agreement
between different codes for a total neutrino mass of $\sum m_{\nu}=0.15$ eV at
redshift $z=0$ (left panel) and $z=1$ (right panel). We use measurements from
our simulations with
$L_{\mathrm{box}}=512\,${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc}$$
and $N_{\mathrm{part}}=512^{3}$. The relative agreement is better at low
redshift, partially due to the fact that the signal amplitude is larger there.

Figure 10: The four panels show bispectrum measurements at redshift $z=1$ in
the simulations with $L_{\mathrm{box}}=512\,h^{-1}\,\rm{Mpc}$ and
$N_{\mathrm{part}}=512^{3}$ for different triangle configurations: squeezed
(top left panel), equilateral (top right panel), and scalene configurations A
and B (bottom panels). In each panel, the top subpanel shows the suppression
ratio of the bispectrum of CDM and baryons for three different neutrino masses
$\sum m_{\nu}\in\\{0.15,0.3,0.6\\}$ eV with respect to the massless case, and
the three bottom subpanels show the respective relative differences of the
various codes when compared to GADGET-3.

Figure 10 is a comparison of simulation codes for different triangle
configurations, for all choices of neutrino masses at redshift $z=1$. We
consider four triangle configuations: squeezed, equilateral and two different
scalene configurations. Squeezed configurations refer to triangles where one
of the side is much shorter than the other two (in Fourier space) which
corresponds to looking at the correlation of a distant point with two points
close to each other. Equilateral configurations, instead, refer to the
correlation of three points at equal distance. Scalene triangles do not have
any specific symmetry. In detail, we consider

  * •

squeezed configurations, for which $k_{1}=k_{2}$ and
$k_{3}=0.012\,h\,\mathrm{Mpc}^{-1}$, plotted as a function of $k_{1}$ as in
Figure 9;

  * •

equilateral configurations, for which $k_{1}=k_{2}=k_{3}$, plotted as a
function of $k_{1}$;

  * •

scalene configurations A, for which $k_{1}=0.79\,h\,\mathrm{Mpc}^{-1}$,
$k_{2}=0.56\,h\,\mathrm{Mpc}^{-1}$, plotted as a function of $k_{3}$;

  * •

scalene configurations B, for which $k_{1}=0.39\,h\,\mathrm{Mpc}^{-1}$,
$k_{3}=0.2\,h\,\mathrm{Mpc}^{-1}$, plotted as a function of $k_{2}$.

In each of the four panels of Figure 10, the top subpanel shows the ratio
between the CDM and baryon bispectrum for massive neutrino cosmologies over
massless ones, which we define in analogy to the case of the power spectrum as

$T_{\mathrm{ccc}}(k_{1},k_{2},k_{3})=\frac{B_{\mathrm{ccc}}^{\mathrm{(massive)}}(k_{1},k_{2},k_{3})}{B_{\mathrm{ccc}}^{\mathrm{(massless)}}(k_{1},k_{2},k_{3})}\,.$
(5.4)

The three bottom subpanels show the relative differences of the measurements
of the suppression ratio in the various codes with respect to GADGET-3 at each
of the three neutrino masses, i.e. $\sum m_{\nu}\in\\{0.15,0.3,0.6\\}$ eV
(from top to bottom). For all configurations considered, discrepancies fall
broadly within the $5\%$ range. As expected, massive neutrinos suppress the
bispectrum of CDM and baryons at all scales, with a stronger effect at smaller
scales. For comparison, we also show the tree-level prediction from
perturbation theory, using

$B_{\mathrm{ccc}}^{\text{(tree-
level)}}(k_{1},k_{2},k_{3})=2F_{2}(\bm{k}_{1},\bm{k}_{2})P_{\mathrm{cb}}^{\mathrm{L}}(k_{1})P_{\mathrm{cb}}^{\mathrm{L}}(k_{2})+2\,\textrm{permutations}\,,$
(5.5)

where

$F_{2}(\bm{k}_{1},\bm{k}_{2})=\frac{5}{7}+\frac{1}{2}\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{1}k_{2}}\left(\frac{k_{1}}{k_{2}}+\frac{k_{2}}{k_{1}}\right)+\frac{2}{7}\frac{(\bm{k}_{1}\cdot\bm{k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}\,,$
(5.6)

and $P_{\mathrm{cb}}^{\mathrm{L}}$ is the linear power spectrum of CDM and
baryons generated by CLASS. We can see in Figure 10 that the suppression ratio
is in good agreement with this prediction when all the three scales
$k_{1},k_{2},k_{3}$ have a small wavenumber, and that the measured suppression
is generally stronger if some of the wavenumbers are large. This is in line
with the results seen in the power spectrum, where strong nonlinearities lead
to additional suppression.

At the largest scales, i.e. when $k_{1},k_{2},k_{3}\lesssim
0.1\,h\,\mathrm{Mpc}^{-1}$, a distinct feature can be discerned which is
particularly prominent in the equilateral configuration and for larger
neutrino masses: the measurements from the various codes separate into two
groups, GADGET-3, openGADGET3, GADGET-4, AREPO, and ANUBIS on the one side,
and L-GADGET3, NM-GADGET4, CONCEPT, PKDGRAV3, SWIFT, gevolution, and COLA on
the other side. The latter group includes all the codes that employ a mesh-
based method, and all of them have means to mitigate against shot noise. We
therefore suspect that this dichotomy originates from shot noise in the
particle method. This could be tested, e.g. by increasing the number of
neutrino particles until convergence is achieved.

As with the case of the power spectrum, we conducted various checks concerning
numerical convergence with respect to finite-volume and resolution effects.
These show a consistent picture that is in line with what we discussed in
Section 5.1.2. As an example, Figure 11 presents the results for the case of
the squeezed configuration in simulations with $N_{\mathrm{part}}=1024^{3}$,
with a larger volume (left panel) or a higher resolution (right panel) than
our simulations with $N_{\mathrm{part}}=512^{3}$. In both cases, the agreement
on smaller scales is improved: for the larger volume this happens because more
independent triangles contribute to each measurement, while for the higher
resolution this is due to the better numerical convergence of the density
field on small scales.

Figure 11: The suppression of the squeezed bispectrum of CDM and baryons
relative to the one measured in the GADGET-3 reference runs at $z=1$ for $\sum
m_{\nu}=0.15\,\rm{eV}$ for the simulations with a larger volume,
$L_{\mathrm{box}}=1024\,h^{-1}\,\rm{Mpc}$ and $N_{\mathrm{part}}=1024^{3}$
(left panel), and at a higher resolution in the small volume,
$L_{\mathrm{box}}=512\,h^{-1}\,\rm{Mpc}$ and $N_{\mathrm{part}}=1024^{3}$
(right panel). Results shown are for the squeezed configuration where
$k_{1}=k_{2}\equiv k$, $k_{3}=0.012\,h\,\mathrm{Mpc}^{-1}$. The grey bands
highlight the interval of $\pm\,0.01$.

Overall we may conclude that the $N$-body methods which produce highly
consistent two-point statistics also tend to agree very well on the three-
point statistics presented here. NM-GADGET4 appears to be an outlier, showing
considerable deviations for squeezed configurations when the sum of the
neutrino masses is large. This is however not unexpected since the final gauge
transformation from Newtonian motion gauge has been neglected here. This
transformation mainly acts at large scales and would therefore affect the
squeezed configurations. At low neutrino mass, where the method works best,
this effect is almost negligible though. It also becomes minimal at redshift
$z=0$ which was set as the target redshift for this method.

###  5.3 Halo mass function

From the halo catalogues produced by Denhf, we estimate the halo mass
functions (considering only the contribution from CDM and baryons) and compare
them to the predictions by Tinker et al. [117], hereafter Tinker10, as well as
Despali et al. [109], hereafter Despali16. For these predictions, we use the
linear power spectra of CDM and baryons calculated by CLASS for the respective
neutrino cosmologies in the modelling of the theoretical halo mass functions.
It has been shown by Costanzi et al. [118] that this approach reproduces the
halo mass function well for neutrino cosmologies. Figure 12 shows the ratio of
the halo mass functions relative to the halo mass function of GADGET-3 at a
neutrino mass of $\sum m_{\nu}=0.15\,\rm{eV}$ for the runs with
$N_{\mathrm{part}}=1024^{3}$, in the large volume where
$L_{\mathrm{box}}=1024\,${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc}$$
(left panel), as well as for the higher-resolution setup with
$L_{\mathrm{box}}=512\,${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc}$$
(right panel). At the high-mass end, the agreement between different codes is
generally very good, and fluctuations are smaller in the larger volume due to
better statistics. At low masses
$M_{\mathrm{200b}}<10^{13}\,h^{-1}$\mathrm{\text{${M}_{\odot}$}}$$, the number
density of halos is underestimated by COLA, PINOCCHIO, gevolution, and ANUBIS
by up to $50\%$ for the lower resolution. At higher resolution the agreement
improves. The relatively poor performance of gevolution in predicting the halo
mass function can be understood from the fact that the code uses a uniform
mesh. This leads to a smoothing of small-scale structures and generally to a
mass estimate of halos that is poorly converged at the low-mass end. COLA
suffers from the same limitation, but the simulations used a mesh with
significantly higher resolution in this case. Specifically, in the runs with
$N_{\mathrm{part}}=1024^{3}$, COLA used a mesh of $3072^{3}$ grid points,
while gevolution used a mesh of $2048^{3}$ grid points.141414Running
gevolution on such a fine mesh, i.e. with less than one particle per cell on
average, was not supported in the public release of the code. In such a
situation, the evolution would become unstable due to the low order of the
finite-difference gradients used in the particle update. A second-order
gradient computation was therefore implemented for this work, a feature that
is made available in a recent patch of the code.

Figure 12: Halo mass functions relative to the one of GADGET-3, for
$N_{\mathrm{part}}=1024^{3}$ at $z=0$ and neutrino mass $\sum m_{\nu}=0.15$
eV. The result for the larger-volume runs with
$L_{\mathrm{box}}=1024\,h^{-1}\,\rm{Mpc}$ is shown in the left panel while the
right panel shows the result for the higher-resolution runs with
$L_{\mathrm{box}}=512\,h^{-1}\,\rm{Mpc}$.

Figure 13 shows the suppression of the halo mass function due to neutrinos
with masses $\sum m_{\nu}\in\\{0.15,\,0.3,\,0.6\\}\,\rm{eV}$ at redshift $z=0$
(left panel) and $z=1$ (right panel). At low halo masses, there is little
suppression, while going to higher masses the number density of halos is more
and more suppressed. The higher the neutrino masses, and the higher the
redshift, the stronger the suppression: at $z=0$ and $\sum
m_{\nu}=0.15\,\rm{eV}$ the suppression goes down to a factor of $0.9$ at halo
masses of $10^{14}\,h^{-1}$\mathrm{\text{${M}_{\odot}$}}$$, while at $z=1$ and
$\sum m_{\nu}=0.6\,\rm{eV}$ the suppression goes down to $0.4$ at the same
halo mass. In analogy to the case of the power spectrum, we define the
suppression ratio with respect to the massless case as

$R=\frac{~{}\frac{\mathrm{d}n^{\mathrm{(massive)}}}{\mathrm{d}\ln
M}~{}}{\frac{\mathrm{d}n^{\mathrm{(massless)}}}{\mathrm{d}\ln M}}\,.$ (5.7)

Figure 14 shows this suppression ratio relative to the one measured from
GADGET-3 for a sum of neutrino masses of $\sum m_{\nu}=0.15$ eV. The different
codes generally agree within $3\%$ on the suppression ratio, even in cases
where the halo mass function was poorly converged in Figure 12. Our data show
a trend that the suppression is about $1\%$ stronger than predicted by the
models of Tinker10 and Despali16. The COLA results are a slight outlier,
agreeing more with these models than with the other simulations.

Figure 13: The halo mass function for different neutrino masses, relative to
the case with massless neutrinos for each code, respectively. We show the
results for the simulations with $L_{\mathrm{box}}=512\,h^{-1}\,\rm{Mpc}$ and
$N_{\mathrm{part}}=512^{3}$ at $z=0$ (left panel) and at $z=1$ (right panel).

###  5.4 Halo bias

Finally, we study the halo bias for a fixed selection of halos defined by a
mass threshold of
$M_{\mathrm{200b}}>10^{13}\,h^{-1}$\mathrm{\text{${M}_{\odot}$}}$$. However,
since the halo mass function shows considerable differences between the
different simulations, sometimes due to the fact that the mass estimate is not
well converged at the low-mass end (certainly for gevolution, COLA, PINOCCHIO,
and ANUBIS), we apply the halo selection as follows. First, we select the
halos above the mass threshold in the reference runs done with GADGET-3. We
may call the size of the selected population $N_{\mathrm{h}}$. Then, for each
other code, we generate the sample by selecting the $N_{\mathrm{h}}$ most
massive halos. The reasoning for this approach is that, while the estimated
masses of individual halos may differ significantly between different codes,
we still expect there to be a tight correlation that largely preserves the
mass ordering. Another way to think about this is to consider a simple
abundance matching of $N_{\mathrm{h}}$ sources, assigned to the centers of the
most massive halos. We compare the bias measurements to the prediction by
Tinker et al. [117] (Tinker10). Here the large-scale bias is estimated from
the mass-dependent peak height of halos in the linear density field. Given the
GADGET-3 halo masses, we model the peak heights using the linear power
spectrum of CDM and baryons calculated by CLASS for the respective neutrino
cosmology. We then take the average of all biases obtained for each halo mass
to get the final prediction. Note that the prediction of Tinker10 is only
expected to work in the linear regime.

Figure 14: Suppression of the halo mass functions relative to the one measured
in GADGET-3 at redshift $z=0$ and for a total neutrino mass of $\sum
m_{\nu}=0.15$ eV. We show the results for the simulations with
$L_{\mathrm{box}}=1024\,h^{-1}\,\rm{Mpc}$ and $N_{\mathrm{part}}=1024^{3}$ in
the left panel, and $L_{\mathrm{box}}=512\,h^{-1}\,\rm{Mpc}$ and
$N_{\mathrm{part}}=1024^{3}$ in the right panel. The grey bands highlight the
interval of $\pm\,0.01$. Figure 15: Halo bias with respect to CDM and baryons
at $z=0$ and neutrino mass $\sum m_{\nu}=0.15$ eV for the simulations with
$L_{\mathrm{box}}=1024\,h^{-1}\,\rm{Mpc}$ and $N_{\mathrm{part}}=1024^{3}$. In
the GADGET-3 reference simulations, all halos with
$M_{\mathrm{200b}}>10^{13}\,h^{-1}$\mathrm{\text{${M}_{\odot}$}}$$ are
selected, providing a sample of size $N_{\mathrm{h}}$. For the other
simulations, we then select the most massive $N_{\mathrm{h}}$ halos. The grey
bands in the lower panels highlight the interval of $\pm\,0.01$.

Figure 15 shows the bias measurements from the simulations with larger volume,
$L_{\mathrm{box}}=1024\,h^{-1}\,\rm{Mpc}$. We define the scale-dependent halo
bias $b(k)$ as the ratio of the cross-power spectrum of halos with CDM and
baryons and the auto-power spectrum of CDM and baryons, i.e.

$b(k)=\frac{P_{\mathrm{h,cb}}(k)}{P_{\mathrm{cb}}(k)}\,.$ (5.8)

It has been shown by Castorina et al. [76] that defining the bias factor with
respect to cold species gives closer-to-universal and less scale-dependent
results than using the same definition with respect to total matter. As can be
seen in the left panel of Figure 15, the bias measurements agree reasonably
well on large scales except for gevolution where the bias is measured to be
about $4\%$ larger, and for PINOCCHIO where the bias is measured to be about
$10\%$ smaller. In analogy to the case of the power spectrum, we again define
a bias ratio with respect to the massless case, which in this situation will
quantify the increase (rather than suppression) of the bias in the presence of
massive neutrinos,

$Q(k)=\frac{b^{\mathrm{(massive)}}(k)}{b^{\mathrm{(massless)}}(k)}\,.$ (5.9)

Results for this bias ratio are shown in the right panel of Figure 15. Here
the agreement between different codes is excellent, well within $1\%$ over
almost the entire range of scales probed. For PINOCCHIO the bias ratio is
about $1\%$ accurate up to $k\simeq 0.3\,h\,\mathrm{Mpc}^{-1}$.

Figure 16: Halo bias with respect to CDM and baryons at $z=0$ and neutrino
mass $\sum m_{\nu}=0.15$ eV for the higher-resolution simulations with
$L_{\mathrm{box}}=512\,h^{-1}\,\rm{Mpc}$ and $N_{\mathrm{part}}=1024^{3}$. In
the GADGET-3 reference simulations, all halos with
$M_{\mathrm{200b}}>10^{13}\,h^{-1}$\mathrm{\text{${M}_{\odot}$}}$$ are
selected, providing a sample of size $N_{\mathrm{h}}$. For the other
simulations, we then select the most massive $N_{\mathrm{h}}$ halos. The grey
bands in the lower panels highlight the interval of $\pm\,0.01$.

To study the robustness of our results with respect to the mass resolution of
the simulations we repeat the bias measurements in the runs with higher
resolution, i.e. with
$L_{\mathrm{box}}=512\,${\mathrm{\mathnormal{h}}}^{-1}\text{\,}\mathrm{Mpc}$$
and $N_{\mathrm{part}}=1024^{3}$. The smaller simulation volume leads to a
higher level of shot noise in the halo counts, which incurs somewhat larger
fluctuations when compared to the larger volume. As can be seen in Figure 16,
the agreement between the different codes is improved significantly, in
particular for codes that have difficulties in predicting the halo mass
function accurately (gevolution, COLA, and ANUBIS). The results for PINOCCHIO
do not improve significantly, as the discrepancy is mainly due to the
approximate nature of the method rather than lack of resolution.

It is worth pointing out that the bias ratio $Q$ has a shape similar to the
inverted power spectrum ratio $S_{\mathrm{cb}}$. This is actually expected,
and was recently discussed by Hassani et al. [119]. It means that the power
spectrum of halos, when selected at fixed mass threshold, is much less
sensitive to the neutrino mass than the power spectrum of CDM and baryons. On
the other hand, synthetic catalogues are often created in such a way that the
observed abundance of a certain type of object is reproduced (abundance
matching). In such a situation, it may be more appropriate to study the
dependence of the bias on the neutrino mass at fixed number count. We
therefore repeat our measurements, but keeping the size of all samples fixed
at $N_{\mathrm{h}}^{\mathrm{(massless)}}$, which is the number of halos with
$M_{\mathrm{200b}}>10^{13}\,h^{-1}$\mathrm{\text{${M}_{\odot}$}}$$ in the
GADGET-3 reference simulation at zero neutrino mass. In other words, when
selecting the halo sample for non-zero neutrino mass, we still select the
$N_{\mathrm{h}}^{\mathrm{(massless)}}$ most massive halos in all simulations.
This effectively reduces the mass threshold of the selection for the massive
neutrino case when compared to the previous procedure.

Figure 17 shows the results of the bias measurement obtained through this
procedure, using the higher-resolution simulations. We observe that the bias
still increases with neutrino mass, but not as much as in Figure 16 where a
fixed halo mass threshold was used. The interplay between bias enhancement and
the suppression of the matter power spectrum therefore depends on the
situation, and needs to be studied carefully in the context of the specific
numerical recipes that are employed in the production of synthetic catalogues.

Figure 17: Halo bias with respect to CDM and baryons at $z=0$ and neutrino
mass $\sum m_{\nu}=0.15$ eV for the higher-resolution simulations with
$L_{\mathrm{box}}=512\,h^{-1}\,\rm{Mpc}$ and $N_{\mathrm{part}}=1024^{3}$. In
the GADGET-3 reference simulation for $\sum m_{\nu}=0$ eV, all halos with
$M_{\mathrm{200b}}>10^{13}\,h^{-1}$\mathrm{\text{${M}_{\odot}$}}$$ are
selected, providing a sample of size $N_{\mathrm{h}}^{\mathrm{(massless)}}$.
We then select the most massive $N_{\mathrm{h}}^{\mathrm{(massless)}}$ halos
in all the other simulations, even for those with non-zero neutrino mass. The
grey bands in the lower panels highlight the interval of $\pm\,0.01$.

##  6 Discussion

Accurate and reliable modelling of the signatures that the neutrino mass
imprints on the observables used to test the cosmological model is an
essential ingredient for the data analysis of all upcoming galaxy surveys, and
in particular for Euclid. Such modelling necessarily requires a self-
consistent description of the linear and possibly nonlinear clustering of
neutrinos along with the nonlinear evolution of dark matter and baryonic
structures. By comparing results across different implementations, including
eleven full $N$-body implementations, two $N$-body schemes with fast time
integration based on Lagrangian perturbation theory, and a further four codes
that predict the nonlinear matter power spectra directly, we establish that
current numerical techniques are in sub-percent agreement with regards to
modelling the impact of massive neutrinos on the most common summary
statistics of cosmological large-scale structure. We identify several specific
situations where larger modelling errors can occur, but such shortfalls are
generally well understood in terms of approximations or other compromises that
were made in these situations. Our results can therefore be used as a detailed
guide for choosing the preferred modelling techniques for any application
given its requirements in terms of resources, accuracy, and quantities that
need to be modelled. The validation presented here is the crucial first step
for building up confidence in the numerical tools employed in the data
analysis pipeline of Euclid. It is particularly vital when considering the
actual measurement of the neutrino mass scale from the data, which is one of
the key science goals of the mission.

The fastest method of predicting simple summary statistics like the nonlinear
matter power spectrum are emulators, and they will therefore play a crucial
role in the cosmological likelihood analysis of Euclid and other large-scale
structure surveys. They are of course many orders of magnitude faster than
simulations but tend to outperform even semi-analytic models which often have
some bottlenecks in their numerical evaluation. However, emulators can only be
as accurate as the simulations they are trained on, and it is therefore
important to understand the modelling errors of simulations too. Overall our
results demonstrate that we are in a fairly comfortable position, with several
independent numerical techniques at our disposal that produce consistent
results at the sub-percent level if we employ them diligently. The community
has also implemented such techniques in a large number of different $N$-body
codes, such that there is no shortage in choice of which code one wants to
use. Moreover, our detailed comparison of particle-based and mesh-based
techniques shows that the assumption of linear neutrinos is clearly sufficient
to reach percent accuracy, even up to scales of the order of $k\approx
7\,h\,\mathrm{Mpc}^{-1}$ relevant for predicting the weak-lensing signal in
Euclid. Here we do of course not attempt to address the additional challenge
of modelling baryonic effects, i.e. astrophysical processes, down to such
scales as this can be treated separately, see e.g. Martinelli et al. [120] for
a discussion. On mildly nonlinear scales and in particular at redshifts
$z\gtrsim 1$ relevant for Euclid, a “sophisticated 3LPT” realisation like the
one produced with PINOCCHIO, based on a scale-dependent linear growth rate
computed from CAMB and propagated to second- and third-order LPT with standard
techniques, can be useful in contexts where emulators are not available and
full simulations would be too expensive.

The summary statistics considered in our analysis include auto- and cross-
power spectra of the CDM and baryon component and neutrinos, bispectra of the
CDM and baryon component, halo mass functions, and halo bias. We present
results for redshifts $z=0$ and $z=1$, relevant for galaxy surveys like
Euclid. We do not consider redshift-space distortions or other effects that
occur due to taking observations on our past light cone, and leave a detailed
investigation of these effects to future work. However, we expect that no big
surprises would appear given our level of confidence in the modelling of the
summary statistics presented here.

In order to aid future code development, we make our reference simulations and
analysis pipelines available via a public repository (see data availability
statement below). This provides a reliable baseline against which further
numerical methods can be validated, and it showcases the current state-of-the-
art in modelling massive neutrinos in cosmology.

## Acknowledgments

JA acknowledges financial support from the Swiss National Science Foundation.
BB was supported by a UK Research and Innovation Stephen Hawking Fellowship
(EP/W005654/1). MBa acknowledges support by the project “Combining Cosmic
Microwave Background and Large Scale Structure data: an Integrated Approach
for Addressing Fundamental Questions in Cosmology”, funded by the MIUR
Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN) Bando 2017 - grant
2017YJYZAH. CH-A acknowledges support from the Excellence Cluster ORIGINS
which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy - EXC-2094 - 390783311. BSW is
supported by a Royal Society Enhancement Award (grant no.
RGF$\backslash$EA$\backslash$181023). RM, DFM, HAW, FH thank the Research
Council of Norway for their support and our computations were performed on
resources provided by UNINETT Sigma2 – the National Infrastructure for High
Performance Computing and Data Storage in Norway. CM acknowledges support from
a UK Research and Innovation Future Leaders Fellowship [grant MR/S016066/1].
The PINOCCHIO simulations were performed on Queen Mary’s Apocrita HPC
facility, supported by QMUL Research-IT. CG acknowledges the support from the
grant PRIN-MIUR 2017 WSCC32 ZOOMING and from the Italian National Institute of
Astrophysics under the grant “Bando PRIN 2019,” PI: Viola Allevato. TC is
supported by the INFN INDARK PD51 grant and the FARE MIUR grant
‘ClustersXEuclid’ R165SBKTMA. KD acknowledges support by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy - EXC-2094 - 390783311 and the COMPLEX project from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program grant agreement ERC-2019-AdG 882679. The
openGADGET3 simulations were carried out at the Leibniz Supercomputer Center
(LRZ) under the project pr86re. CA and BL are supported by an European
Research Council Starting Grant (ERC-StG-716532), and BL also acknowledges
support by the UK Science and Technology Facilities Council Consolidated
Grants No. ST/I00162X/1 and ST/P000541/1. KK is supported by the UK Science
and Technology Facilities Council (grant numbers ST/S000550/1 and
ST/W001225/1). MV and GP are supported by the INFN INDARK PD 51 grant and by
the ASI-INAF n. 2017-14-H.0 agreement. REA and MZ acknowledge the support of
the ERC-STG number 716151 (BACCO). This work used the DiRAC@Durham facility
managed by the Institute for Computational Cosmology on behalf of the STFC
DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS via
STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1,
Durham University and STFC operation grant ST/R000832/1. DiRAC is part of the
UK National e-Infrastructure.

The Euclid Consortium acknowledges the European Space Agency and a number of
agencies and institutes that have supported the development of Euclid, in
particular the Academy of Finland, the Agenzia Spaziale Italiana, the Belgian
Science Policy, the Canadian Euclid Consortium, the French Centre National
d’Etudes Spatiales, the Deutsches Zentrum für Luft- und Raumfahrt, the Danish
Space Research Institute, the Fundação para a Ciência e a Tecnologia, the
Ministerio de Ciencia e Innovación, the National Aeronautics and Space
Administration, the National Astronomical Observatory of Japan, the
Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the
Romanian Space Agency, the State Secretariat for Education, Research and
Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom
Space Agency. A complete and detailed list is available on the Euclid web site
(http://www.euclid-ec.org).

For the purpose of open access, the authors have applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising
from this work and upload the accepted version to arXiv with CC BY licence.

## Data availability

Raw data from our reference runs (GADGET-3), and all reduced data (summary
statistics) presented in this work are available at
https://doi.org/10.5281/zenodo.7297976. We also provide initial data,
parameter specifications for CLASS, and a documented analysis pipeline to
facilitate further validations against our results.

## References

  * [1] R. Laureijs et al., “Euclid Definition Study Report,” arXiv:1110.3193 [astro-ph.CO].
  * [2] E. Di Valentino, A. Melchiorri, and J. Silk, “Cosmological constraints in extended parameter space from the Planck 2018 Legacy release,” JCAP 01 (2020) 013.
  * [3] eBOSS Collaboration: S. Alam et al., “Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Cosmological implications from two decades of spectroscopic surveys at the Apache Point Observatory,” Phys. Rev. D 103 (2021) 083533.
  * [4] N. Palanque-Delabrouille, C. Yèche, N. Schöneberg, J. Lesgourgues, M. Walther, S. Chabanier, and E. Armengaud, “Hints, neutrino bounds and WDM constraints from SDSS DR14 Lyman-$\alpha$ and Planck full-survey data,” JCAP 04 (2020) 038.
  * [5] P. L. Taylor, T. D. Kitching, and J. D. McEwen, “Preparing for the Cosmic Shear Data Flood: Optimal Data Extraction and Simulation Requirements for Stage IV Dark Energy Experiments,” Phys. Rev. D 98 (2018) 043532.
  * [6] Euclid Collaboration: A. Blanchard et al., “Euclid preparation: VII. Forecast validation for Euclid cosmological probes,” Astron. Astrophys. 642 (2020) A191.
  * [7] Euclid Collaboration: M. Knabenhans et al., “Euclid preparation: II. The EuclidEmulator – A tool to compute the cosmology dependence of the nonlinear matter power spectrum,” Mon. Not. Roy. Astron. Soc. 484 (2019) 5509.
  * [8] Euclid Collaboration: M. Knabenhans et al., “Euclid preparation: IX. EuclidEmulator2 – power spectrum emulation with massive neutrinos and self-consistent dark energy perturbations,” Mon. Not. Roy. Astron. Soc. 505 (2021) 2840.
  * [9] Euclid Collaboration: T. Castro et al., “Euclid preparation. XXIV. Calibration of the halo mass function in $\Lambda(\nu)$CDM cosmologies,” arXiv:2208.02174 [astro-ph.CO].
  * [10] Euclid Collaboration: A. Fumagalli et al., “Euclid : Effects of sample covariance on the number counts of galaxy clusters,” Astron. Astrophys. 652 (2021) A21.
  * [11] Euclid Collaboration: A. Fumagalli et al., “Euclid preparation. TBD. Covariance model validation for the 2-point correlation function of galaxy clusters,” in preparation .
  * [12] R. E. Angulo and O. Hahn, “Large-scale dark matter simulations,” arXiv:2112.05165 [astro-ph.CO].
  * [13] J. Brandbyge, S. Hannestad, T. Haugbølle, and B. Thomsen, “The Effect of Thermal Neutrino Motion on the Non-linear Cosmological Matter Power Spectrum,” JCAP 08 (2008) 020.
  * [14] M. Viel, M. G. Haehnelt, and V. Springel, “The effect of neutrinos on the matter distribution as probed by the Intergalactic Medium,” JCAP 06 (2010) 015.
  * [15] S. Agarwal and H. A. Feldman, “The effect of massive neutrinos on the matter power spectrum,” Mon. Not. Roy. Astron. Soc. 410 (2011) 1647.
  * [16] S. Bird, M. Viel, and M. G. Haehnelt, “Massive Neutrinos and the Non-linear Matter Power Spectrum,” Mon. Not. Roy. Astron. Soc. 420 (2012) 2551.
  * [17] F. Villaescusa-Navarro, F. Marulli, M. Viel, E. Branchini, E. Castorina, E. Sefusatti, and S. Saito, “Cosmology with massive neutrinos I: towards a realistic modeling of the relation between matter, haloes and galaxies,” JCAP 03 (2014) 011.
  * [18] E. Castorina, C. Carbone, J. Bel, E. Sefusatti, and K. Dolag, “DEMNUni: The clustering of large-scale structures in the presence of massive neutrinos,” JCAP 07 (2015) 043.
  * [19] J. D. Emberson et al., “Cosmological neutrino simulations at extreme scale,” Res. Astron. Astrophys. 17 (2017) 085.
  * [20] J. Adamek, R. Durrer, and M. Kunz, “Relativistic N-body simulations with massive neutrinos,” JCAP 11 (2017) 004.
  * [21] A. Banerjee, D. Powell, T. Abel, and F. Villaescusa-Navarro, “Reducing Noise in Cosmological N-body Simulations with Neutrinos,” JCAP 09 (2018) 028.
  * [22] J. Brandbyge, S. Hannestad, and T. Tram, “Momentum space sampling of neutrinos in $N$-body simulations,” JCAP 03 (2019) 047.
  * [23] J. Brandbyge and S. Hannestad, “Grid Based Linear Neutrino Perturbations in Cosmological N-body Simulations,” JCAP 05 (2009) 002.
  * [24] Y. Ali-Haimoud and S. Bird, “An efficient implementation of massive neutrinos in non-linear structure formation simulations,” Mon. Not. Roy. Astron. Soc. 428 (2012) 3375.
  * [25] J. Liu, S. Bird, J. M. Z. Matilla, J. C. Hill, Z. Haiman, M. S. Madhavacheril, A. Petri, and D. N. Spergel, “MassiveNuS: Cosmological Massive Neutrino Simulations,” JCAP 03 (2018) 049.
  * [26] J. Z. Chen, A. Upadhye, and Y. Y. Y. Wong, “One line to run them all: SuperEasy massive neutrino linear response in $N$-body simulations,” JCAP 04 (2021) 078.
  * [27] J. Dakin, J. Brandbyge, S. Hannestad, T. Haugbølle, and T. Tram, “$\nu$CO$N$CEPT: Cosmological neutrino simulations from the non-linear Boltzmann hierarchy,” JCAP 02 (2019) 052.
  * [28] K. Yoshikawa, S. Tanaka, N. Yoshida, and S. Saito, “Cosmological Vlasov–Poisson Simulations of Structure Formation with Relic Neutrinos: Nonlinear Clustering and the Neutrino Mass,” Astrophys. J. 904 (2020) 159.
  * [29] J. Brandbyge and S. Hannestad, “Resolving Cosmic Neutrino Structure: A Hybrid Neutrino N-body Scheme,” JCAP 01 (2010) 021.
  * [30] S. Bird, Y. Ali-Haïmoud, Y. Feng, and J. Liu, “An Efficient and Accurate Hybrid Method for Simulating Non-Linear Neutrino Structure,” Mon. Not. Roy. Astron. Soc. 481 (2018) 1486.
  * [31] W. Elbers, C. S. Frenk, A. Jenkins, B. Li, and S. Pascoli, “An optimal non-linear method for simulating relic neutrinos,” Mon. Not. Roy. Astron. Soc. 507 (2021) 2614.
  * [32] C. Partmann, C. Fidler, C. Rampf, and O. Hahn, “Fast simulations of cosmic large-scale structure with massive neutrinos,” JCAP 09 (2020) 018.
  * [33] M. Zennaro, R. E. Angulo, G. Aricò, S. Contreras, and M. Pellejero-Ibáñez, “How to add massive neutrinos to your $\Lambda$CDM simulation - extending cosmology rescaling algorithms,” Mon. Not. Roy. Astron. Soc. 489 (2019) 5938.
  * [34] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, “The fate of hints: updated global analysis of three-flavor neutrino oscillations,” JHEP 09 (2020) 178.
  * [35] KATRIN Collaboration: M. Aker et al., “Improved Upper Limit on the Neutrino Mass from a Direct Kinematic Method by KATRIN,” Phys. Rev. Lett. 123 (2019) 221802.
  * [36] KATRIN Collaboration: M. Aker et al., “Direct neutrino-mass measurement with sub-electronvolt sensitivity,” Nature Phys. 18 (2022) 160.
  * [37] E. Di Valentino, S. Gariazzo, and O. Mena, “Most constraining cosmological neutrino mass bounds,” Phys. Rev. D 104 (2021) 083504.
  * [38] M. Archidiacono, T. Brinckmann, J. Lesgourgues, and V. Poulin, “Physical effects involved in the measurements of neutrino masses with future cosmological data,” JCAP 02 (2017) 052.
  * [39] J. R. Bond, G. Efstathiou, and J. Silk, “Massive Neutrinos and the Large Scale Structure of the Universe,” Phys. Rev. Lett. 45 (1980) 1980.
  * [40] W. Hu, D. J. Eisenstein, and M. Tegmark, “Weighing neutrinos with galaxy surveys,” Phys. Rev. Lett. 80 (1998) 5255.
  * [41] S. Hannestad, A. Upadhye, and Y. Y. Y. Wong, “Spoon or slide? The non-linear matter power spectrum in the presence of massive neutrinos,” JCAP 11 (2020) 062.
  * [42] V. Springel, “The Cosmological simulation code GADGET-2,” Mon. Not. Roy. Astron. Soc. 364 (2005) 1105.
  * [43] V. Springel, J. Wang, M. Vogelsberger, A. Ludlow, A. Jenkins, A. Helmi, J. F. Navarro, C. S. Frenk, and S. D. M. White, “The Aquarius Project: the subhalos of galactic halos,” Mon. Not. Roy. Astron. Soc. 391 (2008) 1685.
  * [44] R. E. Angulo, V. Springel, S. D. M. White, A. Jenkins, C. M. Baugh, and C. S. Frenk, “Scaling relations for galaxy clusters in the Millennium-XXL simulation,” Mon. Not. Roy. Astron. Soc. 426 (2012) 2046.
  * [45] R. E. Angulo, M. Zennaro, S. Contreras, G. Aricò, M. Pellejero-Ibañez, and J. Stücker, “The BACCO simulation project: exploiting the full power of large-scale structure for cosmology,” Mon. Not. Roy. Astron. Soc. 507 (2021) 5869.
  * [46] A. M. Beck et al., “An improved SPH scheme for cosmological simulations,” Mon. Not. Roy. Astron. Soc. 455 (2016) 2110.
  * [47] T. Marin-Gilabert, M. Valentini, U. P. Steinwandel, and K. Dolag, “The role of physical and numerical viscosity in hydrodynamical instabilities,” arXiv:2205.09135 [astro-ph.CO].
  * [48] V. Springel, R. Pakmor, O. Zier, and M. Reinecke, “Simulating cosmic structure formation with the gadget-4 code,” Mon. Not. Roy. Astron. Soc. 506 (2021) 2871.
  * [49] P. Heuschling, C. Partmann, and C. Fidler, “A minimal model for massive neutrinos in Newtonian N-body simulations,” JCAP 09 (2022) 068.
  * [50] V. Springel, “E pur si muove: Galiliean-invariant cosmological hydrodynamical simulations on a moving mesh,” Mon. Not. Roy. Astron. Soc. 401 (2010) 791.
  * [51] R. Weinberger, V. Springel, and R. Pakmor, “The Arepo public code release,” Astrophys. J. Suppl. 248 (2020) 32.
  * [52] J. Dakin, S. Hannestad, and T. Tram, “The cosmological simulation code CONCEPT 1.0,” Mon. Not. Roy. Astron. Soc. 513 (2022) 991.
  * [53] D. Potter, J. Stadel, and R. Teyssier, “PKDGRAV3: Beyond Trillion Particle Cosmological Simulations for the Next Era of Galaxy Surveys,” arXiv:1609.08621 [astro-ph.IM].
  * [54] M. Schaller et al., “SWIFT: SPH With Inter-dependent Fine-grained Tasking.” Astrophysics source code library, May 2018.  https://ascl.net/1805.020.
  * [55] W. Elbers, “Geodesic motion and phase-space evolution of massive neutrinos,” arXiv:2207.14256 [astro-ph.CO].
  * [56] R. Teyssier, “Cosmological hydrodynamics with adaptive mesh refinement: a new high resolution code called ramses,” Astron. Astrophys. 385 (2002) 337.
  * [57] J. Adamek, D. Daverio, R. Durrer, and M. Kunz, “General relativity and cosmic structure formation,” Nature Phys. 12 (2016) 346.
  * [58] J. Adamek, D. Daverio, R. Durrer, and M. Kunz, “gevolution: a cosmological N-body code based on General Relativity,” JCAP 07 (2016) 053.
  * [59] H. A. Winther, K. Koyama, M. Manera, B. S. Wright, and G.-B. Zhao, “COLA with scale-dependent growth: applications to screened modified gravity models,” JCAP 08 (2017) 006.
  * [60] B. S. Wright, H. A. Winther, and K. Koyama, “COLA with massive neutrinos,” JCAP 10 (2017) 054.
  * [61] P. Monaco, E. Sefusatti, S. Borgani, M. Crocce, P. Fosalba, R. K. Sheth, and T. Theuns, “An accurate tool for the fast generation of dark matter halo catalogs,” Mon. Not. Roy. Astron. Soc. 433 (2013) 2389.
  * [62] L. A. Rizzo, F. Villaescusa-Navarro, P. Monaco, E. Munari, S. Borgani, E. Castorina, and E. Sefusatti, “Simulating cosmologies beyond $\Lambda$CDM with PINOCCHIO,” JCAP 01 (2017) 008.
  * [63] M. Cataneo, J. D. Emberson, D. Inman, J. Harnois-Deraps, and C. Heymans, “On the road to per cent accuracy – III. Non-linear reaction of the matter power spectrum to massive neutrinos,” Mon. Not. Roy. Astron. Soc. 491 (2020) 3101.
  * [64] B. Bose, B. S. Wright, M. Cataneo, A. Pourtsidou, C. Giocoli, L. Lombriser, I. G. McCarthy, M. Baldi, S. Pfeifer, and Q. Xia., “On the road to per cent accuracy – V. The non-linear power spectrum beyond $\Lambda$CDM with massive neutrinos and baryonic feedback,” Mon. Not. Roy. Astron. Soc. 508 (2021) 2479.
  * [65] E. Lawrence, K. Heitmann, J. Kwan, A. Upadhye, D. Bingham, S. Habib, D. Higdon, A. Pope, H. Finkel, and N. Frontiere, “The Mira-Titan Universe II: Matter Power Spectrum Emulation,” Astrophys. J. 847 (2017) 50.
  * [66] K. R. Moran, K. Heitmann, E. Lawrence, S. Habib, D. Bingham, A. Upadhye, J. Kwan, D. Higdon, and R. Payne, “The Mira-Titan Universe IV. High Precision Power Spectrum Emulation,” arXiv:2207.12345 [astro-ph.CO].
  * [67] A. Banerjee and N. Dalal, “Simulating nonlinear cosmological structure formation with massive neutrinos,” JCAP 11 (2016) 015.
  * [68] C.-P. Ma and E. Bertschinger, “A Calculation of the full neutrino phase space in cold + hot dark matter models,” Astrophys. J. 429 (1994) 22.
  * [69] S. Tassev, M. Zaldarriaga, and D. Eisenstein, “Solving Large Scale Structure in Ten Easy Steps with COLA,” JCAP 06 (2013) 036.
  * [70] A. Izard, M. Crocce, and P. Fosalba, “ICE-COLA: Towards fast and accurate synthetic galaxy catalogues optimizing a quasi $N$-body method,” Mon. Not. Roy. Astron. Soc. 459 (2016) 2327.
  * [71] A. Lewis, A. Challinor, and A. Lasenby, “Efficient computation of CMB anisotropies in closed FRW models,” Astrophys. J. 538 (2000) 473.
  * [72] A. Lewis and S. Bridle, “Cosmological parameters from CMB and other data: A Monte Carlo approach,” Phys. Rev. D 66 (2002) 103511.
  * [73] P. Monaco, T. Theuns, and G. Taffoni, “Pinocchio: pinpointing orbit-crossing collapsed hierarchical objects in a linear density field,” Mon. Not. Roy. Astron. Soc. 331 (2002) 587.
  * [74] E. Munari, P. Monaco, E. Sefusatti, E. Castorina, F. G. Mohammad, S. Anselmi, and S. Borgani, “Improving fast generation of halo catalogues with higher order Lagrangian perturbation theory,” Mon. Not. Roy. Astron. Soc. 465 (2017) 4658.
  * [75] E. Munari, P. Monaco, J. Koda, F.-S. Kitaura, E. Sefusatti, and S. Borgani, “Testing approximate predictions of displacements of cosmological dark matter halos,” JCAP 07 (2017) 050.
  * [76] E. Castorina, E. Sefusatti, R. K. Sheth, F. Villaescusa-Navarro, and M. Viel, “Cosmology with massive neutrinos II: on the universality of the halo mass function and bias,” JCAP 02 (2014) 049.
  * [77] J. F. Navarro, C. S. Frenk, and S. D. M. White, “A Universal density profile from hierarchical clustering,” Astrophys. J. 490 (1997) 493.
  * [78] C. Fidler, A. Kleinjohann, T. Tram, C. Rampf, and K. Koyama, “A new approach to cosmological structure formation with massive neutrinos,” JCAP 01 (2019) 025.
  * [79] M. Cataneo, L. Lombriser, C. Heymans, A. Mead, A. Barreira, S. Bose, and B. Li, “On the road to percent accuracy: non-linear reaction of the matter power spectrum to dark energy and modified gravity,” Mon. Not. Roy. Astron. Soc. 488 (2019) 2121.
  * [80] A. Mead, S. Brieden, T. Tröster, and C. Heymans, “HMcode-2020: Improved modelling of non-linear cosmological power spectra with baryonic feedback,” arXiv:2009.01858 [astro-ph.CO].
  * [81] A. Cooray and R. K. Sheth, “Halo Models of Large Scale Structure,” Phys. Rept. 372 (2002) 1.
  * [82] M. Cacciato, F. C. v. d. Bosch, S. More, R. Li, H. J. Mo, and X. Yang, “Galaxy Clustering & Galaxy-Galaxy Lensing: A Promising Union to Constrain Cosmological Parameters,” Mon. Not. Roy. Astron. Soc. 394 (2009) 929.
  * [83] C. Giocoli, M. Bartelmann, R. K. Sheth, and M. Cacciato, “Halo model description of the non-linear dark matter power spectrum at $k\gg 1$ Mpc-1,” Mon. Not. Roy. Astron. Soc. 408 (2010) 300.
  * [84] R. K. Sheth and G. Tormen, “Large scale bias and the peak background split,” Mon. Not. Roy. Astron. Soc. 308 (1999) 119.
  * [85] R. K. Sheth and G. Tormen, “An Excursion Set Model of Hierarchical Clustering : Ellipsoidal Collapse and the Moving Barrier,” Mon. Not. Roy. Astron. Soc. 329 (2002) 61.
  * [86] J. S. Bullock, T. S. Kolatt, Y. Sigad, R. S. Somerville, A. V. Kravtsov, A. A. Klypin, J. R. Primack, and A. Dekel, “Profiles of dark haloes. Evolution, scatter, and environment,” Mon. Not. Roy. Astron. Soc. 321 (2001) 559.
  * [87] B. Bose, M. Cataneo, T. Tröster, Q. Xia, C. Heymans, and L. Lombriser, “On the road to per cent accuracy IV: ReACT – computing the non-linear power spectrum beyond $\Lambda$CDM,” Mon. Not. Roy. Astron. Soc. 498 (2020) 4650.
  * [88] VIRGO Consortium Collaboration: R. E. Smith, J. A. Peacock, A. Jenkins, S. D. M. White, C. S. Frenk, F. R. Pearce, P. A. Thomas, G. Efstathiou, and H. M. P. Couchmann, “Stable clustering, the halo model and nonlinear cosmological power spectra,” Mon. Not. Roy. Astron. Soc. 341 (2003) 1311.
  * [89] R. Takahashi, M. Sato, T. Nishimichi, A. Taruya, and M. Oguri, “Revising the Halofit Model for the Nonlinear Matter Power Spectrum,” Astrophys. J. 761 (2012) 152.