Euclid : Modelling massive neutrinos in cosmology – a code comparison

. 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 eﬀort 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 ﬁnd that the diﬀerent numerical implementations produce fully consistent results. We can therefore be conﬁdent 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.


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 eV.Using the well-known relation between neutrino mass and cosmological energy density, m ν ≈ Ω ν h 2 × 94 eV, this lower bound on the sum of neutrino masses corresponds to a lower bound of Ω ν h 2 > 6 × 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 Ω X , which is its present-day energy density in units of the critical density, and physical density parameters are then denoted as ω X = Ω 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-α 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 ≈ 7 h 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 ≈ 0.3 h 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 phasespace 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 "δf method" introduced in Elbers et al. [31] is another hybrid approach that uses a particle ensemble to estimate perturbations δ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.

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] ∆m 2  21 = 7.42 +0.21 −0.20 × 10 −5 eV 2 , ( ∆m 2 31 = 2.517 +0.026 −0.028 × 10 −3 eV 2 (NO), ( ∆m 2  32 = −2.498+0.028 −0.028 × 10 −3 eV 2 (IO), ( where "NO" denotes the normal mass ordering and "IO" the inverted one.This leads to lower bounds on the sum of neutrino masses of m ν ≳ 0.06 eV (NO) and m ν ≳ 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 m ν ≲ 2.4 eV.
However, cosmology already provides much more stringent bounds, typically around m ν ≲ 0.1 − 0.2 eV.Assuming a minimal ΛCDM cosmology with massive neutrinos, a joint analysis of cosmological probes currently obtains m ν < 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 cb , or total matter (which includes massive neutrinos), P m , with respect to a model with massless neutrinos, can be quantified as [18,39,40] respectively, where f ν = Ω ν /Ω m is the neutrino mass fraction.In the nonlinear regime, this suppression is even more prominent and exhibits a dip at k ≈ 1 h 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.

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

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 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 highdensity regions.Methods featured here include Tree-PM, adaptive mesh refinement (AMR), fast multipole method (FMM), particle-particle/particle-mesh (P 3 M), and moving mesh.No adaptive force computation is used in gevolution, COLA, and PINOCCHIO, which all work with a uniform mesh.[66,67] 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 δf method [31].This method works by decomposing the distribution function f into an analytical background component f and a perturbation δ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 δ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,30,68].
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 [69] where v is the peculiar velocity, p is the canonical momentum, m is the particle's rest-mass, a is the scale factor, and Ψ 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 p 2 ≪ m 2 c 2 a 2 in which the equations simplify to 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 ∝ a −1 in this case.
Even though the evolution of high-momentum particles suffers from severe errors, including the breakdown of causality for 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].

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 where ζ(k) is the comoving curvature perturbation of the mode 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.Ali-Haïmoud & Bird [24] solve for the transfer function of the neutrino perturbations using the nonlinear matter power spectrum of the simulation to construct an effective source term, assuming that the phase correlation between neutrino and matter perturbations remain largely intact even at nonlinear scales.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.

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.

COLA
The COLA approach by Tassev et al. [70] produces fast, approximate simulations of cosmological structure formation.Essentially, instead of solving for a full particle trajectory 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) δx(t) = x(t) − x 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. [61], which also included an implementation in the MG-PICOLA simulation code1 by Winther et al. [60].This implementation was carried over to the COLA solver within FML2 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 Λ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. [71].Therefore, we use a slightly modified3 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. [72,73] 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. [71].

PINOCCHIO
The PINOCCHIO code4 [62,74,75] 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 [76].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. [63] and is based on the result of Castorina et al. [18,77] 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. [63], which was limited to 2LPT, we extend here the computation to third order: as shown by Munari et al. [75] 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 [78] 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 ν,m (k), can be approximated by P ν,m (k) = P L ν (k)P 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.
We do not expect PINOCCHIO to be competitive with N -body codes in predicting the matter power spectrum: taking it as a sophisticated implementation of 3LPT, we expect it to lose power on scales smaller than k = 0.3 h Mpc −1 for the halo power spectrum and k = 0.2 h Mpc −1 for the matter power spectrum.It will not be competitive with COLA as well which, being a PM code, can converge to the solution (on scales larger than the mesh) if a sufficient number of time steps is used.This better accuracy comes at a higher cost in computing time, by approximately a factor of eight in the configuration used in this paper (see also Blot et al. [79] for a similar benchmark), as well as in memory since our COLA runs use a grid three times finer than the mean particle separation.The PINOCCHIO code is widely used especially to characterise the covariance of galaxy clustering measurements, thanks to its low computational cost and its ability to generate halo catalogues.

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 Nbody simulations.We would like to stress that our method is exact in the weak-field limit of general relativity (see Fidler et al. [80]) 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 [80]).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].

Halo-model reaction
The halo model reaction approach provides the nonlinear corrections caused by massive neutrinos to a ΛCDM power spectrum through a ratio of halo-model predictions.Following Cataneo et al. [81], the nonlinear power spectrum is then given by where R(k, z) is the halo-model reaction and P NL pseudo (k, z) is the nonlinear pseudo power spectrum.
The pseudo spectrum is a nonlinear Λ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 Λ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. [64,81] and pass the modified linear spectrum as produced by CAMB to HMcode developed by Mead et al. [82], with ΛCDM presets, i.e. no baryonic feedback nor massive neutrinos.The benefit of using the pseudo rather than Λ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 [83][84][85].
Following Cataneo et al. [64], the halo-model reaction for massive neutrinos is given by with "cb" denoting the CDM and baryon component and "ν" 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 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 ), where ρ is the background density for the specified matter species.The halo mass functions are given as ) The peak heights are defined as where the subscript "sc" indicates this quantity is calculated by solving the standard ΛCDM spherical-collapse equations but with the indicated matter density.The mass fluctuation variances are given by In all predictions from halo-model reaction, we employ a Sheth-Tormen halo mass function [86,87], a power-law concentration-mass relation (see for example the work by Bullock et al. [88]), and an NFW halo density profile [78].The predictions are computed numerically using the public code ReACT5 by Bose et al. [65,89].

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. [90] and later improved by Takahashi et al. [91], and the HMcode predictor by Mead et al. [82], 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 interpolation 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 [66,[92][93][94], polynomial chaos expansion [7,8], or neural network approaches [95].
In addition to the predictions from halofit and HMcode mentioned above, we focus in this paper on the Cosmic Emu [66], 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 [67,96], a suite of 111 simulations run with the HACC code [97].The simulations are distributed over an eightdimensional cosmological parameter space comprising (ω m , ω b , ω ν , σ 8 , h, n s , w 0 , w a ), where σ 8 is the present-day amplitude of linear matter density fluctuations at the scale of 8 h −1 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 Mpc −1 within the redshift range z ∈ [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 (Ω m , Ω b , m ν , A s , h, n s , w 0 , w a ), where A s is the amplitude of primordial perturbations at the scale k p = 0.05 Mpc −1 , and includes redshifts of z ∈ [0, 3] and modes up to k = 10 h 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 [98].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 [99].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 (Ω cb , Ω b , m ν , σ 8 , h, n s , w 0 , w a ) and covers a redshift range of z ∈ [0, 1.5] for modes up to k = 5 h Mpc −1 .The claimed precision is better than three percent.
We refer to the original references for more information about the emulators.

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 m ν (including the massless case), different box sizes L box , and different mass resolutions to check for numerical convergence with respect to finite-volume and discretisation effects.N 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 [100].We keep the total matter density (at redshift zero) fixed at Ω m = 0.319 by adjusting the CDM density parameter together with the neutrino mass.The baryon density is fixed at Ω b = 0.049, and the remaining cosmological parameters are A s = 2.215 × 10 −9 at the pivot scale k p = 0.05 Mpc −1 , n 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 K and the effect of radiation is taken into account in the simulations at the linear level, either by using carefully tailored 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 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.

Initial conditions
The initial conditions of all simulations are generated at redshift z = 127. 6The linear matter power spectra and transfer functions are obtained by running either CAMB or the CLASS Table 2. Overview of the basic parameters used in our simulation suite.The cases with m ν = 0 eV and m ν = 0.15 eV are the two main baselines for our comparison, but we include some cases with larger masses, up to m ν = 0.6 eV, to probe more "extreme" regions of parameter space. Simulation Boltzmann code by Blas et al. [101].These files are then used by the REPS7 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. [102].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 [103].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 [104].
We employ a modified version of the N-GenIC code8 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 particlebased 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 [105,106].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,9 using geodesic integration from high redshift [20,69], 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.

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.

Power spectra
The power spectra of the different snapshots have been estimated using Pylians3. 10The routine first deposits particle masses into a regular 3D grid with N 3 voxels using the cloudin-cell mass-assignment scheme.In this work, we always use a mesh with N 3 = N part such that the Nyquist scales match between particles and mesh.Although using larger grids may improve measurements on smaller scales, we recommend caution due to potential systematics and advise against relying on results near or beyond the Nyquist scale set by the mean particle separation.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, |δ(k)| 2 .The modes are then binned in intervals of width equal to the fundamental frequency k F = 2πL −1 box and the power spectrum is finally estimated as where N i is the number of independent modes in the considered bin To compute the cross-power spectrum of two fields instead, the estimator is generalised in the most straightforward way, 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 ± 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.

Bispectra
We measure the bispectrum of CDM and baryons (ccc) using the estimator l q 2 ∈Bm q 3 ∈Bn δ K (q 1 + q 2 + q 3 ) δ(q 1 ) δ(q 2 ) δ(q 3 ), (4.4) where N tr is the number of "fundamental triangles", formed by the vectors q i satisfying the triangle condition q 1 + q 2 + 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. [107].In order to compare the large-volume simulations (L box = 1024 h −1 Mpc) more easily with the small-volume ones, we use the same k-space binning in both cases, fixing the bin width to k 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.

Halo catalogues
For the considered snapshots of the various simulations, we identify halos with the code Denhf [108][109][110][111] 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,77].
Denhf estimates the local density at the position of each N -body particle by calculating the distance to its 10 th -nearest neighbour d 10 , and assigning to each particle a density that is proportional to d −3 10 .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 applying the rescaling of halo masses that translates the halo mass function of Watson et al. [112], that has been used to calibrate the code, to the one of Tinker et al. [113].Such a rescaling has been used, e.g., by Fumagalli et al. [10] to force the halo mass function averaged over 1000 realisations to follow a target one.We compute the rescaling only once, in the case of massless neutrinos, and use it for all neutrino masses.
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 ± 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.

Power spectra
A key prediction from neutrino simulations is a suppression of the matter power spectrum that exceeds the maximum linear theory prediction of ∆P m /P m ≈ 8f ν .The matter power spectrum can be decomposed as follows where P cb is the power spectrum of CDM and baryons, P ν,cb the cross-power spectrum of neutrinos with CDM and baryons, and P ν 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 m (k) with various power spectrum emulators, including EuclidEmulator2 which only predicts this quantity.

CDM and baryons
The leading contribution to P m (k) is P cb (k), which is suppressed in massive neutrino models.
Figure 1 shows the ratio of P cb (k) for models with a neutrino mass of m ν = 0.15 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 ≈ 0.1 h 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 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 ≈ 0.1 h Mpc −1 , ANUBIS up to k ≈ 0.3 h Mpc −1 , gevolution up to k ≈ 0.5 h Mpc −1 , and COLA up to k ≈ 0.7 h Mpc −1 .The other codes stay within a 1% deviation from GADGET-3 for all scales down to the particle Nyquist scale of 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 ≈ 0.2 h Mpc −1 .
Figure 2 shows the ratio of P cb (k) for models with massive neutrinos relative to the massless case for m ν ∈ {0.15, 0.3, 0.6} eV.For later convenience, we denote this suppression ratio by S cb (k), or in general where the subscript "X" denotes any component in question, here X → cb.Linear calculations (taken from CLASS) predict that, on the largest scales, the growth of structure is mostly unaffected so that S cb (k) approaches unity, while on small scales S cb (k) reaches a plateau.Compared to this linear expectation, all codes in the comparison predict slightly less suppression around k = 0.1 h Mpc −1 and a much greater suppression for k > 0.3 h 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 [41,114].At z = 0, we obtain excellent agreement between all simulations and most approximate methods up to the scale of maximum suppression at k max ≈ 1 h Mpc −1 , where the suppression is 20% greater than the linear prediction.At z = 1, the scale of maximum suppression shifts to k max ≈ 2 h Mpc −1 and the differences are greater, both with the linear prediction and between the codes, with the exception of PINOCCHIO11 which fares significantly better compared to z = 0.
To study these relative differences in greater detail, we show S cb (k) for the smallest neutrino mass, m ν = 0.15 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 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 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 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 ≈ 0.3 h 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 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 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 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. [115] between RAMSES, GADGET-3, and PKDGRAV3.
The observed differences are greater when the neutrino mass is increased, especially for the approximate methods.Figure 4 shows S cb (k) relative to the GADGET-3 prediction for models with m ν = 0.3 eV (left panel) and m ν = 0.6 eV (right panel) at z = 0.While the BACCOemulator agrees to better than 1% with GADGET-3 up to k Nyq for 0.3 eV, it makes no prediction for 0.6 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 Mpc −1 for total neutrino masses of 0.3 eV and 0.6 eV, respectively.The differences with the linear prediction (CLASS) and with PINOCCHIO are large, as noted above for 0.15 eV.
The agreement is excellent for the other codes, but some subtle differences can be discerned.On large scales, we observe that SWIFT, CON CEPT, 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 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 Mpc −1 .SWIFT, PKDGRAV3, and CON CEPT are low compared to GADGET-3 for k > 1 h Mpc −1 , unlike what was seen in Figure 3. ANUBIS diverges from GADGET-3 by more than 1% for k > 1 h Mpc −1 for the 0.6 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 Nyq .

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 box is doubled at fixed resolution, the agreement remains excellent on linear scales and is sometimes even slightly better around k Nyq , providing an important consistency check for most codes.Increasing instead the mass resolution by doubling k 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 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.

Contributions from neutrinos
The subdominant contributions to P m (k) are the cross-power spectrum P ν,cb (k) between neutrinos and the CDM and baryon component, and the auto-power spectrum P ν (k) of neutrinos.As can be seen from Eq. (5.1), these are suppressed by the small factors f ν and f 2 ν , respectively, and they are themselves additionally strongly suppressed with respect to P 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 ν,cb (k) for various codes relative to the result from the SWIFT code, for the smallest neutrino mass, m ν = 0.15 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 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 CON CEPT have no shot noise, but depart from the other codes for k > 0.2 h Mpc −1 .The relative difference to SWIFT is 8% at k = 1 h Mpc −1 .Much the same applies to the neutrino auto-power spectrum, P ν (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, n−1 , from the measured power spectrum, where n−1 = 0.125 h −3 Mpc 3 in our high-resolution simulations.For SWIFT the subtracted values are derived from the δf method.The difference between SWIFT and CON CEPT is 19% (off the chart) at k = 1 h 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 CON CEPT than with the other N -body codes.

Total matter
Some of the codes, in particular the EuclidEmulator2, only provide predictions for the power spectrum of total matter, P m (k).In Figure 7, we show the relative agreement of different emulators and other codes predicting P m (k) to our high-resolution reference run with GADGET-3.The left panel shows the result for the matter power spectrum itself, P m (k), at m ν = 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 m (k), with respect to the massless scenario.Cosmic Emu shows a disagreement larger than 1% around k ≈ 1 h 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.Overall our results are in fair agreement with the detailed comparison carried out by Parimbelli et al. [116].

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 [117][118][119][120].The total matter bispectrum in the presence of massive neutrinos can be schematically defined as 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 "ccν" and cold-neutrino-neutrino "cνν" terms that are symmetrised, as indicated by the superscript "(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. [121].As for the power spectrum case, the leading term is the one of CDM and baryons.Figure 8 shows a comparison of the bispectrum B ccc as measured at redshift z = 1 by all simulation codes12 for all triangle configurations and all neutrino masses considered, for the set of runs with N part = 512 3 .Triangles are plotted as a function of The comparison shows the maximum percentage difference with respect to GADGET-3 of all triangle configurations at each value of k max .For L-GADGET3, openGADGET3, GADGET-4, PKDGRAV3, CON CEPT, and NM-GADGET4, discrepancies are mostly within 5% for m ν = 0.0 eV, while steadily growing up to around 10% for massive neutrino cosmologies.Other simulation codes are within 10% already at m ν = 0.0 eV.ANUBIS and gevolution show some of the strongest deviations at large k max , but this is mainly a result of finite resolution as we find much better agreement in the higher-resolution runs.At low k 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 Mpc −1 and show the agreement between different codes for a total neutrino mass of m ν = 0.15 eV at redshift z = 0 (left panel) and z = 1 (right panel).We use measurements from our simulations with L box = 512 h −1 Mpc and N 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 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 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 Mpc −1 , k 2 = 0.56 h Mpc −1 , plotted as a function of k 3 ; • scalene configurations B, for which k 1 = 0.39 h Mpc −1 , k 3 = 0.2 h 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 . (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. m ν ∈ {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 where and P L cb 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 ≲ 0.1 h 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, CON CEPT, 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 part = 1024 3 , with a larger volume (left panel) or a higher resolution (right panel) than our simulations with N 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.
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.

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. [113], hereafter Tinker10, as well as Despali et al. [111], 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. [122] 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 m ν = 0.15 eV for the runs with N part = 1024 3 , in the large volume where L box = 1024 h −1 Mpc (left panel), as well as for the higher-resolution setup with L box = 512 h −1 Mpc (right panel).At the highmass 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 200b < 10 13 h −1 M ⊙ , the number density of halos is underestimated by COLA, 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 part = 1024 3 , COLA used a mesh of 3072 3 grid points, while gevolution used a mesh of 2048 3 grid points. 13igure 13 shows the suppression of the halo mass function due to neutrinos with masses m ν ∈ {0.15, 0.3, 0.6} 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 m ν = 0.15 eV the suppression goes down to a factor of 0.9 at halo masses of 10 14 h −1 M ⊙ , while at z = 1 and m ν = 0.6 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 = dn (massive)  d ln M dn (massless)  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 m ν = 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.

Halo bias
Finally, we study the halo bias for a fixed selection of halos defined by a mass threshold of M 200b > 10 13 h −1 M ⊙ .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, 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 h .Then, for each other code, we generate the sample by selecting the N 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 h sources, assigned to the centers of the most massive halos.We compare the bias measurements to the prediction by Tinker et al. [113] (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 15 shows the bias measurements from the simulations with larger volume, L box = 1024 h −1 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. 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,

b(k)
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 ≃ 0.3 h Mpc −1 .
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 box = 512 h −1 Mpc and N 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 cb .This is actually expected, and was recently discussed by Hassani et al. [123].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 (massless) h , which is the number of halos with M 200b > 10 13 h −1 M ⊙ 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 (massless) h 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 quite as much as in Figure 16 where a fixed halo mass threshold was used.The large-scale bias indicated by the dotted line (Tinker10) is about 0.01, or one percentage point, lower so that the corresponding bias ratio Q drops by 0.7%.The largest effect in Q appears around k ≈ 1 h Mpc −1 were the change can be as much as 2% due to the different halo selection.The interplay between bias enhancement and the suppression of the matter power spectrum has another layer of complexity due to the way in which the sample is selected.This needs to be studied carefully in the context of the specific numerical recipes that are employed in the production of synthetic catalogues.

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 .We then select the most massive N (massless) h 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 ± 0.01.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.We find that the current state of the art for emulators yields an absolute precision on the power spectrum of total matter better than 2% and can predict the relative change due to the neutrino mass to better than 1% on all scales considered in this work.Interestingly, the best semi-analytic fitting methods available, in particular ReACT and HMcode, can achieve similar performance.
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 ≈ 7 h Mpc −1 relevant for predicting the weak-lensing signal in Euclid.We note that some codes have inherent difficulties reaching such levels of absolute accuracy due to effects of finite resolution.This is obvious in the cases where a uniform mesh is employed in the computation of gravitational interactions (gevolution and COLA), but AMR does not solve the issue entirely as the example of ANUBIS illustrates.The relative impact of massive neutrinos can nonetheless be predicted very accurately with those codes.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. [124] for a discussion.On mildly nonlinear scales and in particular at redshifts z ≳ 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 to produce a large number of halo catalogues in a limited amount of computing time, see e.g.Fumagalli et al. [10].
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.

Figure 1 .
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 m ν = 0.15eV in the simulation with L box = 512 h −1 Mpc and N 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 ± 0.01.

Figure 2 .
Figure 2. Suppression of the power spectrum of CDM and baryons at z = 0 and z = 1 for three different neutrino masses, m ν ∈ {0.15, 0.3, 0.6} eV, when compared to the massless case.Results are from the simulations with L box = 512 h −1 Mpc and N part = 512 3 .The corresponding particle Nyquist wavenumber is indicated by the grey dash-dotted line.

Figure 3 .
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 m ν = 0.15 eV for the simulations with L box = 512 h −1 Mpc and N part = 512 3 .The corresponding particle Nyquist wavenumber is indicated by the grey dash-dotted line.The grey bands highlight the interval of ± 0.01.

Figure 4 .
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 m ν ∈ {0.3, 0.6} eV for the simulations with L box = 512 h −1 Mpc and N part = 512 3 .The corresponding particle Nyquist wavenumber is indicated by the grey dash-dotted line.The grey bands highlight the interval of ± 0.01.

Figure 5 .
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 m ν = 0.15 eV for the simulations with a larger volume, L box = 1024 h −1 Mpc and N part = 1024 3 (left panel), and at a higher resolution in the small volume, L box = 512 h −1 Mpc and N 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 ± 0.01.

Figure 6 .
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 m ν = 0.15 eV for the higherresolution simulation with L box = 512 h −1 Mpc and N part = 1024 3 .The Nyquist wavenumber is indicated by the grey dash-dotted line.The grey bands highlight the interval of ± 0.01.

Figure 7 .
Figure 7.Total matter power spectrum P m (k) at z = 0 for m ν = 0.15 eV from emulators and fitting methods (left panel), and the respective suppression S m (k) with respect to the massless case (right panel), compared to the higher-resolution reference GADGET-3 simulation with L box = 512 h −1 Mpc and N 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 ± 0.01.

Figure 8 .Figure 9 .
Figure8.The coloured bands show the scatter of B 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 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.

1 -Figure 10 .
Figure 10.The four panels show bispectrum measurements at redshift z = 1 in the simulations with L box = 512 h −1 Mpc and N 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 massesm ν ∈ {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 11 .
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 m ν = 0.15 eV for the simulations with a larger volume, L box = 1024 h −1 Mpc and N part = 1024 3 (left panel), and at a higher resolution in the small volume, L box = 512 h −1 Mpc and N part = 1024 3 (right panel).Results shown are for the squeezed configuration where k 1 = k 2 ≡ k, k 3 = 0.012 h Mpc −1 .The grey bands highlight the interval of ± 0.01.

14 M 14 MFigure 12 .
Figure 12.Halo mass functions relative to the one of GADGET-3, for N part = 1024 3 at z = 0 and neutrino mass m ν = 0.15 eV.The result for the larger-volume runs with L box = 1024 h −1 Mpc is shown in the left panel while the right panel shows the result for the higher-resolution runs with L box = 512 h −1 Mpc.

14 M 14 MFigure 13 .
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 box = 512 h −1 Mpc and N part = 512 3 at z = 0 (left panel) and at z = 1 (right panel).

. 8 ) 12 10 13 10 14 M 14 MFigure 14 .
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 m ν = 0.15 eV.We show the results for the simulations with L box = 1024 h −1 Mpc and N part = 1024 3 in the left panel, and L box = 512 h −1 Mpc and N part = 1024 3 in the right panel.The grey bands highlight the interval of ± 0.01.

Figure 15 .
Figure 15.Halo bias with respect to CDM and baryons at z = 0 and neutrino mass m ν = 0.15 eV for the simulations with L box = 1024 h −1 Mpc and N part = 10243 .In the GADGET-3 reference simulations, all halos with M 200b > 10 13 h −1 M ⊙ are selected, providing a sample of size N h .For the other simulations, we then select the most massive N h halos.The grey bands in the lower panels highlight the interval of ± 0.01.

k [h Mpc − 1 ]Figure 16 .
Figure16.Halo bias with respect to CDM and baryons at z = 0 and neutrino mass m ν = 0.15 eV for the higher-resolution simulations with L box = 512 h −1 Mpc and N part = 10243 .In the GADGET-3 reference simulations, all halos with M 200b > 10 13 h −1 M ⊙ are selected, providing a sample of size N h .For the other simulations, we then select the most massive N h halos.The grey bands in the lower panels highlight the interval of ± 0.01.

k [h Mpc − 1 ]Figure 17 .
Figure 17.Halo bias with respect to CDM and baryons at z = 0 and neutrino mass m ν = 0.15 eV for the higher-resolution simulations with L box = 512 h −1 Mpc and N part = 1024 3 .In the GADGET-3 reference simulation for m ν = 0 eV, all halos with M 200b > 10 13 h −1 M ⊙ are selected, providing a sample of size N (massless) h