An optimal nonlinear method for simulating relic neutrinos

Cosmology places the strongest current limits on the sum of neutrino masses. Future observations will further improve the sensitivity and this will require accurate cosmological simulations to quantify possible systematic uncertainties and to make predictions for nonlinear scales, where much information resides. However, shot noise arising from neutrino thermal motions limits the accuracy of simulations. In this paper, we introduce a new method for simulating large-scale structure formation with neutrinos that accurately resolves the neutrinos down to small scales and significantly reduces the shot noise. The method works by tracking perturbations to the neutrino phase-space distribution with particles and reduces shot noise in the power spectrum by a factor of $\mathcal{O}\left(10^2\right)$ at $z=0$ for minimal neutrino masses and significantly more at higher redshifts, without neglecting the back-reaction caused by neutrino clustering. We prove that the method is part of a family of optimal methods that minimize shot noise subject to a maximum deviation from the nonlinear solution. Compared to other methods we find permille level agreement in the matter power spectrum and percent level agreement in the large-scale neutrino bias, but large differences in the neutrino component on small scales. A basic version of the method can easily be implemented in existing N-body codes and makes it possible to run neutrino simulations with significantly reduced particle load. Further gains are possible by constructing background models based on perturbation theory. A major advantage of this technique is that it works well for all masses, enabling a consistent exploration of the full neutrino parameter space.


INTRODUCTION
The discovery of neutrino masses (Fukuda et al. 1998;Ahmad et al. 2002;Eguchi et al. 2003) calls for extensions of the Standard Model of particle physics and provides the only known form of dark matter. Measuring the masses is crucial for understanding their origin and for constraining cosmological parameters. While the neutrino mass squared differences are known to a few percent, the absolute masses are unknown and there remain two possible mass orderings: normal and inverted. A rich experimental programme is aimed at determining the mass ordering, measuring the mass scale set by the lightest neutrino and completing the overall picture of neutrino properties. Cosmology plays a vital role in this programme due its ability to provide an independent and complementary constraint on the sum of neutrino masses, (Bond et al. 1980;Hu et al. 1998) with a potential sensitivity below 0.02 eV (Font-Ribera et al. 2014;Chudaykin & Ivanov 2019;Sprenger et al. 2019).
Ongoing and planned neutrino experiments will establish the mass ordering with a discovery expected by the end of the decade. Although oscillation data have shown persistent hints of normal ordering, this preference has decreased to 1.6 over the past year (Esteban et al. 2020). The mass ordering can be established by exploiting matter effects in long baseline neutrino oscillation experiments, as in (Acciarri et al. 2015), and in the Earth for atmospheric neutrinos, as in approach is challenging, so information from multiple sources is essential. Single -decay is the experiment of choice for direct mass searches and provides a model-independent determination of their value. The experiment is ongoing and has put a bound of < 1.1 eV on the value of quasi-degenerate neutrino masses with the aim of reaching < 0.2 eV in the near future (Aker et al. 2019). Project 8 will have the potential to set a limit of < 0.04 eV (Esfahani et al. 2017). Neutrinoless double -decay can also provide information on neutrino masses (Bilenky et al. 2001;Pascoli & Petcov 2002;Nunokawa et al. 2002), albeit entangled with the value of the Majorana CP-violating phases and affected by uncertainty in the nuclear matrix elements (Vergados et al. 2016). For a recent review see e.g. Giuliani et al. (2019).
The complementarity between these different strategies is of great interest. A cosmological measurement of would provide a target for direct mass searches (Drexlin et al. 2013;Mertens 2016). An incompatibility between the two would indicate a non-standard cosmological evolution or new neutrino properties. A cosmological bound of < 100 meV would suggest a normal mass ordering, which should be confronted with evidence from neutrino experiments. Finally, there is a strong synergy with neutrinoless double -decay. Knowing the mass ordering and the sum of neutrino masses would narrow down the range of values for the effective Majorana mass parameter, providing a clear target for future experiments.
Measuring the mass scale, and potentially ruling out the inverted mass ordering, is therefore a major target of near-term cosmological surveys, including (Levi et al. 2013), (Laureĳs et al. 2011), and (Ivezic et al. 2009). In order to analyse these surveys and to extract a mass measurement, there has been a substantial effort to model precisely the effects of massive neutrinos on structure formation. From the analytical side, a swathe of new techniques such as time RG perturbation theory (Upadhye 2019) and effective field theories (Senatore & Zaldarriaga 2017;Colas et al. 2019), promise to push the validity of perturbation theory into the quasi-linear régime. In the nonlinear régime, -body simulations offer the most accurate picture of structure formation. Yet incorporating neutrinos intobody simulations has proved to be a challenge and some doubts remain about the validity of neutrino simulations on small scales.
The main obstacle to simulating neutrinos is that, in contrast to cold dark matter and baryons, neutrinos have a significant velocity dispersion. This effectively turns the 3-dimensional problem of structure formation, for which -body simulations are well suited, into a 6-dimensional phase-space problem. If no provisions are made, a far greater number of simulation particles is needed to sample properly the phase-space manifold. A further complication arises from the fact that neutrinos are relativistic at high redshifts, such that simulations need to handle both the régime where neutrinos are best described as radiation and the régime where neutrinos are better described as massive particles.
The first neutrino simulations were carried out by Klypin & Shandarin (1983) and Frenk et al. (1983), when neutrinos were thought to be much more massive and the velocity dispersion not as problematic. Modern simulations with sub-electronvolt neutrinos were pioneered by Brandbyge et al. (2008); Viel et al. (2010). Neutrinos are most commonly included in simulations as particles whose initial velocity is the sum of a peculiar gravitational component and a random component sampled from a Fermi-Dirac distribution (Brandbyge et al. 2008;Viel et al. 2010;Bird et al. 2012;Villaescusa-Navarro et al. 2014;Adamek et al. 2017;Banerjee et al. 2018;Emberson et al. 2017;Inman et al. 2015;Villaescusa-Navarro et al. 2015;Castorina et al. 2015;Villaescusa-Navarro et al. 2020). The main difficulty with particle simulations is shot noise caused by the velocity dispersion. This problem is more severe for the smallest neutrino masses, which are observationally most relevant. Because neutrinos are a subdominant component, the error in the total matter distribution is relatively small. However, shot noise obscures the small-scale behaviour of the neutrinos and is clearly undesirable if one is interested in the neutrino component and its effect on structure formation.
To overcome the problems with particle simulations, grid simulations evolve the neutrino distribution using a system of fluid equations, which requires a scheme to close the moment hierarchy at some low order (Brandbyge & Hannestad 2009;Viel et al. 2010;Hannestad et al. 2012;Archidiacono & Hannestad 2016;Banerjee & Dalal 2016;Dakin et al. 2019;Tram et al. 2019;Inman & Yu 2020), or as a linear response to the non-relativistic matter density (Ali-Haïmoud & Bird 2012;Liu et al. 2018;McCarthy et al. 2018). Even more efficiently, but in the same spirit of treating neutrinos perturbatively, the total effect of neutrinos has been included as a post-processing step in the form of a gauge transformation (Partmann et al. 2020). While these approaches do not suffer from shot noise, they are not able to capture the full nonlinear evolution of the neutrinos at late times. This problem becomes more severe for more massive neutrinos, but is present even for minimal neutrino masses. A number of hybrid simulations have therefore combined grid and particle methods (Brandbyge & Hannestad 2010;Banerjee & Dalal 2016;, typically transitioning from a fluid method to a particle method at some redshift when the neutrinos become nonlinear. Another interesting alternative is to integrate the Poisson-Boltzmann equations directly on the grid (Yoshikawa et al. 2020).
The method proposed in this paper can be considered as a type of hybrid method that integrates neutrino particles but only uses the information contained in the particles to the extent that it is necessary. This is accomplished by dynamically transitioning from a smooth background model to a nonlinear model at the individual particle level. It relies on the noiseless (but approximate) background model as much as possible, thereby producing the smallest amount of shot noise possible whilst solving the full nonlinear system. The main idea is to decompose the phase-space distribution function ( , , ) into a background model 0 ( , , ) which can be solved without noise, and a perturbation which is carried by the simulation particles: The choice of background model is arbitrary, but the method performs best whenever 0 ( , , ) is strongly correlated with ( , , ), in a way that will be made precise below. If the choice of background model is poor, the method performs no worse than an ordinarybody simulation, except for the small amount of overhead associated with evaluating 0 ( , , ). Note that the background model is just an approximation of and can itself be a perturbed Fermi-Dirac distribution.
This type of method has a long history in other fields and is variably known as the method of 'perturbation particles' or more commonly as the ' method', which is the name we shall adopt. Merritt (1987) and Leeuwin et al. (1993) discussed the method of perturbation particles in stellar dynamics. Around the same time, the method arose in plasma physics (Tajima & Perkins 1983;Parker & Lee 1993;Dimits & Lee 1993;Aydemir 1994). While the method of perturbation particles is not widely known today in astronomy, the method is standard fare in plasma physics. A major difficulty in astronomical applications is the absence of a background model that captures enough of the dynamics to be useful. In contrast, plasma physicists are often interested in turbulent phenomena arising in an otherwise stable system, with a natural candidate for a background model 0 at hand. Our work is motivated by the fact that there is also a natural background model for cosmic neutrinos, namely the phasespace density predicted by perturbation theory. There is a major synergy between -body simulations proposed here and work on improved perturbation theory methods. A better background model means a smaller dependence on the particles and therefore further reduced shot noise. We will show however that even the 0 th order approximation, which is just a homogeneous redshifted Fermi-Dirac distribution, provides a significant improvement over ordinarybody methods.
The remainder of the paper is structured as follows. In Section 2, we derive the method and describe its use as a variance reduction method for -body simulations. We also show that the method is part of a family of optimal hybrid methods. In Section 3, we illustrate the method with a one-dimensional test problem. In Section 4, we discuss how the method can be embedded in relativistic simulations. Our suite of simulations is then described in Section 5. The method is compared with commonly used alternatives in Section 6. We consider higher order background models based on perturbation theory in Section 7. Finally, we conclude in Section 8.

DERIVATION
The phase-space evolution of self-gravitating collisionless particles is described by the Poisson-Boltzmann equations, which in the single-fluid case read ≡ + · ∇ − ∇Φ · ∇ = 0, Here, Φ is the gravitational potential, the energy density, and the phase-space density. In general, the Liouville operator acts on each fluid separately and the potential should be summed over all fluid components. In relativistic perturbation theory, this system can be written as a hierarchy of moment equations for the neutrinos, which is solved to first order with Boltzmann codes such as (Lesgourgues 2011) or (Lewis & Challinor 2011). To extend our predictions to the nonlinear régime, we can use -body codes, which solve the Poisson-Boltzmann system by the method of characteristics. Characteristic curves satisfy By construction, one finds that d /d = = 0 along these curves. To infer statistics of the phase-space distribution, such as the number density, ( , ) = 1 , we simulate of these trajectories. Let ( , ) be the initial sampling distribution of Lagrangian marker particles in phase space. Then, a phase-space statistic is given by The usual approach is to set ( , ) = ( , , 0 ), in which case the sum reduces to a simple average over marker particles. The error in our estimate of is then / √ . Hence, if the distribution, ( , , ), has a large intrinsic scatter, we need a large to beat down the noise. Alternatively, we might construct an estimator with a smaller error. Let us therefore write the phase-space distribution function, , as a background model, 0 , together with some perturbation, : ( , , ) = 0 ( , , ) + ( , , ).
We can reduce the error by only using the particles to estimate the perturbed distribution, . We replace (1) with This is useful if can be computed efficiently and if and 0 are strongly correlated, so that the second term is small. The simplest choice of background model is a homogeneous Fermi-Dirac distribution 0 ( , , ) = (2 ) 3 1 /( with internal degrees of freedom. Here, = ( ) is the scale factor and the present-day momentum. Since the noise reduction scales with the correlation between 0 and , we can achieve further gains by adding more information to the background model. The obvious next step is to use perturbation theory to improve on (2). This option is considered in section 7.

Implementation
To implement the method in cosmological -body simulations, we simply replace the particle mass with a weighted mass: The weights are computed by comparing the true phase-space density with the background model. We know the true density for each particle, because = = 0 along characteristic curves. We record the two numbers and at the initial position of each particle and use these to compute the weights, , during each subsequent step. We note that any sampling distribution ( , ) is valid provided that ≠ 0 almost everywhere ≠ 0. We will continue to use the common choice, = , where is the (perturbed) Fermi-Dirac distribution. In general, the optimal choice of will depend on the phase-space statistic of interest. Choosing a distribution, , that oversamples slower particles can lead to an additional reduction in shot noise.
Given the homogeneous Fermi-Dirac background model (2), the neutrino density becomes Cosmological -body simulations only compute the perturbed potential, since the background density¯is accounted for in the background equations. The only change affecting the force calculation is therefore the weighting of the particles. The mean squared weight, = 1 2 2 , is a convenient statistic to quantify the importance of including the neutrino particles. We show the evolution of for a = 100 meV simulation with the homogeneous background model (2) in Fig. 1. At early times, particles deviate very little from their initial trajectory and the weights are negligible. We find that = 4×10 −7 at = 20, = 3×10 −6 at = 10, and = 2×10 −5 at = 5. This early reduction is important as shot noise at high redshifts inhibits the growth of physical structure and can seed additional fluctuations that grow by gravitational instability. At late times, when nonlinear effects become important, the weights increase to = 2 × 10 −4 at = 2, = 1 × 10 −3 at = 1, and = 6.7 × 10 −3 at = 0, independently of the starting redshift of the simulation. This translates to a reduction in shot noise, = 2 / , or an effective increase in particle number at = 0 by a factor (2 ) −1 = 75. Finally, we note that one can save computational resources by integrating only a fraction of the neutrino particles as long as remains small. We do not consider this possibility here.

Variance reduction
The method is an application of the much more general control variates method (Ross 2012;Aydemir 1994). This is a variance reduction technique commonly used in Monte Carlo simulations. See Chartier et al. (2020) for another recent application in cosmology. We briefly review the method here. Let be a random variable with an unknown expectation E[ ] = A. Given independent random samples , the standard estimator is given bŷ The error inÂ is

Let
be another random variable for which the expected value E[ ] = B is known. By adding and subtracting, we can construct a control variate estimator for :

for any constant . LikeÂ, this is an unbiased estimator of E[ ].
However, the error inÂ cv is given by Therefore, the error can be reduced if and are correlated. Differentiating, we see that the optimal value of is given by For the Fermi-Dirac model considered above, * is very close to unity and we simply set = 1 at all times. In general, the value of * could be estimated at runtime. This is useful if we add more information about the unknown variable and extend the method to a linear combination of multiple control variates (see section 7). Furthermore, the method can still be useful when a control variate is not exactly known but can be estimated more efficiently than .

Optimality
Let us consider how the method compares to other methods. To allow for the broadest possible comparison, we will write down an arbitrary hybrid method that involves some background model, 0 ( , , ), such as a fluid description or linear response, and a discrete sampling of the distribution with arbitrary particle weights, ( ): where ( ) is a weight function for the background. This parametrization captures virtually all existing methods. The ordinary -body particle method corresponds to ( , ) = (0, 1) at all times. Pure grid-based methods have ( , ) = (1, 0). Existing hybrid methods switch over from a grid method to a particle method after some time , which corresponds to ( , ) = (1 − , ) with ( ) = [ ≥ ] a step function. For simplicity, we consider only the case where all particles are switched on at the same time, but the argument extends readily to the more practical case where only some particles are switched on. Given a choice of weight function, ( ), for the background, what choice of particle weights is optimal?
Let ( , , ) be the nonlinear distribution and ( , , ) the sampling distribution of the markers. In the continuous limit, the expected error in the number density is given by Meanwhile, the shot noise term in the power spectrum grows as the square of the particle weights, so we want to minimize subject to the constraint ≤ for some maximum error . Assume that the bound is saturated. First, let us look for solutions that extremize the integral constraint. We find the unique solution This is the method introduced above, with optimal given by (3). Any further solution should extremize the Lagrangian, Writing down the Euler-Lagrange equations one finds a family of quadratic solutions The case = 0 corresponds to the trivial solution = 0. For ≠ 0, we obtain the minima These solutions correspond to small perturbations around the method that trade some accuracy for a possible reduction in shot noise. However, since the leading correction is ∝ ( ) 2 , this is only possible if the background model is skewed with respect to the nonlinear solution. Typically, the skewness and the additional reduction in shot noise is negligible. In fact, since the next-to-leading correction is positive, shot noise increases if the skewness is small. We have shown that within the broad class of hybrid methods described by equation (4), -type methods of the form (6) minimize the amount of shot noise, subject to the constraint that the error in the number density remains below a certain bound. The method given by (5), recovered from (6) in the limit → ∞, is the unique solution for which the expected error = 0. The optimal value of is given by (3), but will be close to 1 if 0 ≈ . This is the method we will use exclusively, with the choice = 1.
Density profiles for the 1-dimensional elliptical sine wave test problem. We used = 100 position bins to create the empirical density profiles. On the left, an ordinary -body simulation with = 10 6 particles was used. On the right, the -body simulation was extended with a step.

ONE-DIMENSIONAL EXAMPLE
We now illustrate the method by applying it to a one-dimensional test problem with a known solution. Readers that are satisfied with the mathematical derivation may skip ahead to section 4.

The elliptical sine wave
Consider the 1-dimensional collisionless Boltzmann equation where the particles move under a conservative force ( ) = −Φ ( ). Let us assume a periodic potential given by The steady-state solution can be found to be: in terms of the background density 0 and velocity dispersion . The corresponding density profile ( ) is given by To find the general time-dependent solution, we use the method of characteristics. The characteristic equations are These equations of motion can be solved in terms of the energy = 1 2 2 + sin 2 ( /2), which gives where 0 is an integration constant and sn( ) is the Jacobi elliptic sine function with elliptic modulus = 1/ √ (Weisstein 2002) 1 .
1 For → ∞, we have sn → sin , meaning that ∝ . The particle 'ignores' the potential. For = = 1, sn = tanh , meaning the particle asymptotically approaches a potential peak. For < 1, the particle is bounded and oscillates between peaks.
Assuming a homogeneous Gaussian distribution with mean¯for the initial momenta 0 , the general solution, ( , , ), at later times is a complicated expression involving elliptic sines and arcsines. The details are given in Appendix A. We replicate the problem using -body methods. A large number of particles are initialized on the interval ∈ [0, 4 ] with momenta drawn from the initial distribution (7). The particles are then integrated using In addition to the ordinary -body method, we use a method, where the background model is given by , and the weights are updated during each step via = / . The corresponding density profiles are shown in Fig. 2. The plots were created using = 10 6 particles and the model parameters are 0 = = 1 and¯= 10. The results show that both the ordinary -body simulation and the simulation with a step can reproduce the exact solution. However, the ordinary method is very noisy, whereas the method reproduces the expected profiles with remarkable accuracy. The reason for this is that while the distribution itself has a large dispersion, resulting in noisy results for the ordinary method, the perturbations from the steady solution are small, which allows the method to work. This is exactly analogous to the cosmic neutrino background.

RELATIVISTIC EFFECTS
Neutrinos constitute a relativistic fluid at early times, which introduces some subtleties when evolving such a fluid with a Newtonian code. Including relativistic effects is not necessary for the method, but we include them in our simulations to allow for a consistent comparison with recent works (Adamek et al. 2017;Tram et al. 2019;Partmann et al. 2020). Furthermore, the higher order methods discussed in section 7 provide a natural setting for including these effects without neglecting the nonlinear evolution of the neutrinos. We will work in the Newtonian motion framework of Fidler et al. (2017a) and make modifications to the initial conditions, long-range force calculation, and particle equations of motion as outlined below.

Initial Conditions
To generate initial conditions for massive neutrinos and to set up the higher order background models (section 7), accurate calculation of the linear theory neutrino distribution function is indispensable. This can be done with the Boltzmann codes (Lewis & Challinor 2011) and (Lesgourgues 2011). At their default settings, these codes produce accurate total matter and radiation power spectra (their intended purpose), but the neutrino related transfer functions (e.g. density and velocity) are not converged and can be very inaccurate (Dakin et al. 2019). To obtain converged results, we post-process perturbation vectors from by integrating source functions up to multipole ℓ max = 2000. This prevents the artificial reflection that can happen for low ℓ max . See appenfix B for more details.
Initial conditions are then created using the post-processed transfer functions from in -body gauge at = 100. We do not follow the usual approach of back-scaling the present-day power spectrum, but use the so-called forward Newtonian motion approach (Fidler et al. 2015(Fidler et al. , 2017a. To our knowledge, forward Newtonian motion initial conditions have always been set up with the Zel'dovich approximation. However, this approximation is known to be inadequate for precision simulations (Crocce et al. 2006). To go beyond Zel'dovich initial conditions, we determine the Lagrangian displacement vectors = − by solving the Monge-Ampère equation This equation is solved numerically with a fixed-point iterative algorithm that exploits the fact that the density perturbation is small. We note that this approach is not equivalent to Lagrangian perturbation theory, but merely provides a more accurate map from the Eulerian initial density field to a Lagrangian displacement field compared to the Zel'dovich approximation. A detailed analysis of this method will be presented elsewhere. Velocities were determined independently using the transfer functions for the velocity dispersion = · . Neutrino particles were displaced randomly in phase space according to the perturbed phase-space density function, PT ( , , ), including terms up to ℓ = 5.

Long-range forces
In the weak-field limit, working in -body gauge, the continuity and Euler equations for a collisionless fluid can be written as (Fidler et al. 2015(Fidler et al. , 2017bBrandbyge et al. 2017): where overdots denote conformal time derivatives, is the density contrast, v the peculiar velocity, and = / 2 . The scalar potential receives contributions from all fluid components: where the sum runs over cold dark matter, baryons, neutrinos, and photons. The -body gauge term, ∇ Nb , arises from the anisotropic stress of relativistic species. In the absence of such species, the continuity and Euler equations agree with the Newtonian equations solved in conventional -body codes. This is what makes -body gauge useful as it allows one to set up initial conditions in -body gauge, evolve them in a Newtonian simulation, and give the results a relativistic interpretation. The relativistic corrections become relevant at the 0.5% level on the largest scales in our Gpc simulations. We include the contribution from photons and massless neutrinos (and for some runs, the massive neutrinos 2 ) by realizing the corresponding transfer functions from on a grid as part of the long-range force calculation in our -body code .

Particle content
When simulating light neutrinos from high redshifts, we are evolving relativistic particles in a Newtonian simulation. Such particles can reach superluminal speeds when evolved using the ordinary equations of motion. Following Adamek et al. (2016), we address this issue by replacing the equations of motion with semi-relativistic equations that are valid to all orders in : Here, −1 u is the spatial part of the 4-velocity. Additionally, when computing the energy density, we replace the weighted mass of the particles with a weighted energy = √︃ 2 + 2 . We have verified that both changes result in virtually identical power spectra at = 0, but with the added benefit of slightly reducing the number of timesteps needed to integrate the particles.

SIMULATIONS
We now describe our neutrino simulations, which were run on the 6 computing facility in Durham. We have implemented the method in the cosmological hydrodynamics code (Schaller et al. 2016(Schaller et al. , 2018. uses a task-based parallelization paradigm to achieve strong scaling on large clusters and obtain significant speed-ups over competing -body codes. The main simulations presented in this paper use the basic version of the method with a homogeneous Fermi-Dirac distribution as background model. Our choice of cosmological parameters, based on Planck 2018 (Aghanim et al. 2018), is (ℎ, Ω c + Ω , Ω b , , ) = (0.6737, 0.265, 0.0492, 2.097 × 10 −9 , 0.9652). We run two sets of m ν = 100 meV

Linear theory
Particle method δ f method m ν = 500 meV Figure 3. Neutrino density plots of (256 Mpc) 3 cubes at = 0, simulated with two commonly-used methods and with the method. The particle and simulations used = 1024 3 particles. Shot noise is clearly visible for the particle method, although noticably less so for = 500 meV. The linear theory model fails to reproduce the small-scale behaviour. The method solves both problems. The inset zooms in on a neutrino halo, which should be compared with the linear response prediction in Fig. 4. simulations at different resolution to test the large-scale and smallscale behaviour of various methods. The cube sizes and particle numbers are listed in Table 1.
The best terrestrial constraint on the absolute mass scale comes from the detector, which places a bound of < 1.1 eV at the 90% C.L. on the effective neutrino mass (Aker et al. 2019). Assuming a quasi-degenerate mass spectrum, this corresponds to a neutrino mass sum of ≈ 3 eV. Recent cosmological limits are much stronger and are quoted below at the 95% C.L. Assuming a degenerate mass spectrum, the Planck temperature, polarization, and lensing likelihoods give a constraint of < 0.24 eV or < 0.26 eV, depending on the details of the high-ℓ polarization analysis (Aghanim et al. 2018). Adding BAO data from BOSS DR12, MGS, and 6dFGS, Choudhury & Hannestad (2020) found < 0.12 eV (degenerate), < 0.15 eV (normal), and < 0.17 eV (inverted). An analysis of the shape of the BOSS DR11 redshift-space power spectrum, combined with CMB data and Type 1a supernovae leads to < 0.18 eV (Upadhye 2019). Adding instead the latest SDSS DR14 BOSS and eBOSS Lyman-forest data to the Planck and BAO data leads to the strongest constraint: < 0.09 eV (Palanque-Delabrouille et al. 2020).
Given these limits, we consider three values for , keeping the present-day value Ω m,0 = Ω cb,0 + Ω ,0 fixed. Scenario one contains three massless neutrinos, scenario two corresponds to the inverted mass ordering with = 100 meV 3 , and scenario three to a degenerate spectrum with = 500 meV. The first two models bracket the most interesting range of values 0 < < 100 meV. The last model has surely been ruled out, but is included for several reasons. First of all, the method reduces to the ordinary particle method in the large mass limit at late times. Hence, the = 500 meV case provides a useful consistency check. Second, when simulations are used to emulate statistics for parameter extraction, we should allow for unlikely excursions in MCMC analyses without our simulation methods breaking down (Partmann et al. 2020). Finally, The two massive scenarios considered in this paper have degenerate neutrino masses (2 × 50 meV and 3 × 167 meV). However, the method can easily be extended to account for mass splittings. In that case, particles would be labelled with a given mass state, , and each state would have its own background model, 0, . The reduction in shot noise is largest for the smallest neutrino masses, placing different masses on a level footing. This allows for better load balancing between different neutrino masses.

RESULTS
We compare our neutrino method with three commonly used alternatives. The most common alternative is the ordinary -body particle method, which is the same in every respect as our method, but with the weighting step disabled. Next, we consider a linear theory method based on Tram et al. (2019) that does not evolve neutrino particles but instead realizes the linear theory neutrino perturbation in -body gauge on a grid. The neutrinos are then fully accounted for in the long-range forces. Finally, we consider the linear response method of Ali-Haïmoud & Bird (2012) in which the neutrino perturbation is calculated by applying the linear theory transfer function ratio, lin ( )/ lin cdm+b ( ), to the simulated cdm + baryon phases. A visual inspection of the neutrino density plots shown in Figs. 3 and 4 reveals the strengths and weaknesses of the four methods. Broadly, we see that the linear theory method does not suffer from shot noise, but fails to reproduce the small-scale behaviour resolved by the particle and methods. At the same time, shot noise is clearly visible in the particle simulation with = 100 meV, despite using = 1024 3 particles in a 256 Mpc cube. This is evidently cured in the plot. We also see that shot noise is much less of a problem for = 500 meV, but the plot is still less grainy than the corresponding particle plot. Finally, Fig. 4 shows that the linear response method greatly improves on the pure linear theory prediction, but still produces neutrino halos that are too diffuse compared to the particle and simulations.

Neutrino component
We start with an analysis of the probability density function of the neutrino density field, computed on a 1024 3 grid from the 256 Mpc simulations. Refer to the plots in Fig. 5, which bear out the basic picture sketched above. For the = 100 meV neutrinos, the particle method is plagued by shot noise, but agrees with the method in the high density tail where the particle number is sufficient to obtain a good signal-to-noise ratio. The linear prediction fails in the high and low density tails. Finally, the linear response method, which We compare the method with three commonly used alternatives. The particle and methods agree in the high density tail, because the largest overdensities have enough particles to achieve a high signal-to-noise ratio. Shot noise plagues the particle method, particularly in underdense regions. The linear methods fail in the high density tail. applies the linear theory ratio ( )/ cdm+b ( ) to the cdm+baryon phases, is an intermediate case between the linear theory and methods. For the more massive scenario, the situation is much the same, except that shot noise is much less of a problem for the particle method on these scales.
Next, we consider two-point statistics and show the neutrino power spectrum at = 0 in Fig. 6, combining the large and small simulations to show a wide range of scales. We use the Gpc simulations for < 0.1 Mpc −1 and the 256 Mpc simulations for ≥ 0.1 Mpc −1 . As expected, all methods agree on scales greater than = 0.1 Mpc −1 for both neutrino masses. On smaller scales, linear theory significantly underpredicts the amount of neutrino clustering. The linear response method also underpredicts the neutrino power spectrum, but not by as much. The relative difference between the nonlinear power spectrum and linear power spectrum is greater for neutrinos than for cdm and baryons. To account for this effect, we fit a nonlinear correction to the linear response power spectrum using the measured power spectrum up to = 1 Mpc −1 : and find = 0.006 ± 0.004 and = 0.90 ± 0.01 ( = 100 meV) and = −0.06 ± 0.03 and = 0.34 ± 0.09 ( = 500 meV). These are shown as the red curves in Fig. 6.
The particle simulations are clearly affected by shot noise, at the level of / = 1/64, obscuring the neutrino signal on scales smaller than = 0.2 Mpc −1 for the lightest scenario and on scales smaller than = 1 Mpc −1 for the more massive scenario. Using the method, shot noise is significantly reduced in the former case (factor of 87) and slightly reduced in the latter case (factor of 3.5), revealing a signal down to = 1 − 2 Mpc −1 . Hence, simulations can achieve a similar resolution independently of mass without adjusting the particle number.
We also show the cross-spectral coefficient which captures phase differences between the dark matter and neutri-nos. By definition, ,cb = 1 according to the linear response method. However, this does not hold on small scales as can be seen in the bottom panels. Up to the point where shot noise becomes a problem, the particle and methods agree, demonstrating that ,cb < 1. This is particularly clear for = 500 meV. Next, we consider how well the simulations can resolve the extended neutrino halos surrounding galaxies and clusters (Brandbyge et al. 2010;Villaescusa-Navarro et al. 2011). In Fig. 7, we show stacked neutrino profiles for halos with virial mass cdm+b in the range (5, 12) × 10 14 . The particle and methods agree almost perfectly, once again because of the high signal-to-noise ratio in the largest overdensities. In linear theory, the neutrino halos are completely absent as is evident also from the cross-sections in Fig. 3. Finally, the linear response method predicts neutrino halos that are too diffuse compared to the nonlinear simulations, and with too little dispersion from the mean profile. The larger dispersion found in the nonlinear simulations is not due to errors in individual profiles, but due to a stronger correlation between cdm+b and the local neutrino density.

Neutrino bias
On larger scales, the neutrino density field can be reconstructed from the density of halos for a given neutrino mass spectrum (Inman et al. 2015). We therefore construct the halo overdensity field, by calculating the number density, h ( ), of halos and the mean density,¯h, at = 0 in our Gpc simulations identified using the halo finder (Elahi et al. 2019). We restrict attention to halos with virial mass, cdm+b > 10 12 , and smooth h and with a tophat filter of comoving radius = 30ℎ −1 Mpc. Following Yoshikawa et al. (2020), we study the mean neutrino density at constant halo density¯( h ), defined in terms of the joint probability density function ( , h ) as (1 + 2 nl + 2 stoch ) −1/2 . This model is analogous to the nonlinear stochastic galaxy biasing model of Taruya & Suto (2000); Yoshikawa et al. (2001). We compute the four quantities ( , 2 nl , 2 stoch , ,h ) for each of the methods under consideration. The results are listed in Table 2 and the biasing relationship is shown in Fig. 8. As expected on these large scales, we find good agreement with differences of a few percent in the bias. The greater the level of neutrino clustering resolved by a given method, the greater the bias and correlation ,h . The stochasticity follows the opposite pattern. The nonlinearity follows no such pattern, but is very small in each case except (amusingly) for the linear theory runs. This is because linear theory does not resolve neutrino halos, causing thē ( h ) relation to level off in the high density tail.
The bias = 0.103 for the 100 meV scenario is in excellent agreement with the bias = 0.071 found by Yoshikawa et al. (2020), when the difference in mass ordering is factored in using the approximately linear relationship between neutrino mass and bias in their results. Yoshikawa et al. (2020) do not consider neutrino masses beyond 400 meV, but our finding of = 0.256 for 500 meV is slightly lower than expected when extrapolating from their results. We also find a larger stochasticity and smaller correlation than might be expected, although the small nonlinearities agree. Given the mutual agreement between the different runs in Table 2, these differences are unlikely to be due to our choice of neutrino method. Differences in the -body code or the identification of halos could also affect this comparison.

Matter power spectrum
The suppression of the total matter power spectrum at = 0, relative to a massless neutrino cosmology, is shown in Fig. 9. We see that all methods are in excellent agreement and reproduce the famous spoon-like feature, which has recently been explained in terms of the halo model (Hannestad et al. 2020). The linear theory method disagrees on the largest scales, because its neutrino component is not affected by cosmic variance unlike the other methods. The differences between the methods are otherwise most pronounced around = 0.6 Mpc −1 , where the suppression is largest. The inset graphs zoom in on these scales. For both neutrino masses, the method predicts a smaller suppression than the particle and linear methods. This is in line with expectation, as the additional small-scale neutrino clustering slightly offsets the suppression. Similarly, the pure linear theory method predicts the least neutrino clustering and the largest suppression. It is interesting to see that the particle and methods do not agree for Σ = 500 meV, despite having similar neutrino power spectra at = 0. This is most likely due to shot noise at high redshift, highlighting the importance of using a hybrid type method. The differences between the methods are at the permille level, corresponding to a shift in neutrino mass of several meV. In absolute terms, the differences are larger for = 500 meV, but less important overall.
The horizontal line corresponds to the empirical fitting formula, Δ / = −9.8 (Brandbyge et al. 2008). Compared to this formula, we find a slightly greater suppression in each case, regardless of the method used to model the neutrinos. For the 100 meV simulations, this can be attributed to our use of the inverted mass ordering. The = 500 meV case is perhaps more surprising, but seems to be in line with recent works. For example, Partmann et al. (2020) find increasingly larger differences with the fitting formula for increasing masses, although they do not consider models with > 300 meV.
Globally, the agreement between these very different methods is an encouraging sign and suggests that we have a good handle on the effects of massive neutrinos on the matter power spectrum. The differences, at most a few permille, may perhaps be relevant when trying to distinguish the effects of individual neutrino masses (Wagner et al. 2012).

HIGHER ORDER METHODS
The performance of the method scales with the correlation between the nonlinear solution ( , , ) and the background model 0 ( , , ), so it is worth investigating other background models. We can go beyond the 0 th order Fermi-Dirac model by including the linear theory prediction. In that case, the distribution function can be written as where the perturbation is decomposed into multipole moments (Ma & Bertschinger 1995), Here, ℓ (·) are Legendre polynomials and the coefficients Ψ ℓ satisfy an infinite hierarchy of moment equations. The Legendre representation yields simple expressions for the first few fluid moments, but is cumbersome for evaluating the distribution function itself. For our purposes, it is more convenient to use the following monomial representation where for a given ℓ max , the odd (even) Φ ℓ ( , , ) can be expressed in terms of all the odd (even) Ψ ( , , ) with ≤ ℓ. See Appendix C for details. With this choice of background model, the density integral becomes with particles weights = / and 0 given by (9). Here, ( ) is the linear neutrino overdensity, which is calculated using . The effect of the perturbation should now be included in the long-range force calculation.
As shown in Fig. 10, adding the multipoles Φ 0 and Φ 1 significantly improves the correlation and therefore reduces the shot noise by almost 50%. It is likely that higher order terms could contribute meaningfully too, as the multipole expansion converges only slowly. However, most of the gain is due to the 0 th order term, which on its own is much easier to implement.

DISCUSSION AND CONCLUSIONS
Shot noise in -body simulations is a major obstacle to modelling the nonlinear evolution of light relic neutrinos. In this paper, we demonstrate that the method, which decomposes the neutrino distribution into an analytically tractable background component, 0 , and a nonlinear perturbation, , carried by the simulation particles, is an effective variance reduction technique. The reduction in shot noise scales with the dynamic particle weights, parametrized by = 1 2 2 . Because the weights are negligible until very late times, the simulation is mostly immunized against the effects of shot noise. Furthermore, shot noise is greatly reduced even at = 0, which makes it possible to resolve neutrino clustering down to much smaller scales than is possible with conventional methods. Using higher order versions of the method, which incorporate additional information from perturbation theory, shot noise can be reduced by another factor of O (2), and possibly more if moments ℓ > 2 are included. Additional reduction in shot noise is possible by carefully tuning the sampling distribution of the marker particles.
The reduction in shot noise is more significant for smaller neutrino masses, because faster particles deviate less from their initial trajectory, resulting in smaller weights. This is fortunate as shot noise is most problematic for the fastest neutrinos. More generally, particles whose trajectories are not perturbed have negligible weights, whereas particles that are captured by halos have appreciable weights. This is again fortunate, because particles are needed in the vicinity of halos where grid methods tend to fail, while the unperturbed particles contain no information and contribute only noise. In between these extremes, particles will have intermediate weights. In this way, the method ensures an optimal combination of particles and background.
The method can in principle be combined with any grid or fluid background model to obtain an optimal hybrid method. Any simulation that evolves neutrino particles can be extended with a weighting step to minimize the shot noise as outlined in section 2.1. It is not necessary, as was done here, to evolve the neutrino particles from the beginning. The -statistic from a reference simulation can be used to gauge when the neutrinos become nonlinear and at what point they can safely be introduced (see Fig. 1).
We know from neutrino oscillations that at least one neutrino has a mass 0.05 eV. Our results indicate that even for masses close to that bound, neutrinos are not particularly well modelled by linear approximations. For instance, the linear response neutrino power spectrum is off by 10% (60%) at = 0.1 Mpc −1 ( = 1 Mpc −1 ) at = 0, and the pure linear theory prediction is off by 14% (96%). Because the neutrinos make up only a small fraction of the total mass, the effect on the matter power spectrum is at most a few permille. This is the level at which the mass splittings are important (Wagner et al. 2012). Other statistics may be affected at a greater level, particularly if they are more sensitive to neutrino effects. For example, we have shown that the neutrino bias relative to dark matter halos is affected at the percent level on 30ℎ −1 Mpc scales. In addition, some novel probes may require accurate modelling of the neutrino dynamics around halos, such as the neutrino-induced dynamical friction (Okoli et al. 2017) and torque (Yu et al. 2019) on halos. By reducing shot noise without neglecting nonlinear terms, the method makes it feasible to calculate these effects even for the lightest neutrinos. Bilenky S. M., Pascoli S., Petcov S., 2001, Phys. Rev. D, 64, 053010 and , which can be calculated accurately with much lower settings 4 . Therefore, we make the assumption that we can decouple the potential terms from most of the neutrino moments Ψ ℓ . We first evolve all source functions in at a reasonable precision setting. This gives the metric perturbations ℎ( , ) and ( , ), which we then take as given and use to integrate the multipole moments Ψ ℓ at high precision where they are needed.