Massive Neutrinos and the Non-linear Matter Power Spectrum

We perform an extensive suite of N-body simulations of the matter power spectrum, incorporating massive neutrinos in the range M = 0.15-0.6 eV, probing the non-linear regime at scales k<10 hMpc-1 at z<3. We extend the widely used HALOFIT approximation to account for the effect of massive neutrinos on the power spectrum. In the strongly non-linear regime HALOFIT systematically over-predicts the suppression due to the free-streaming of the neutrinos. The maximal discrepancy occurs at k ~ 1 hMpc-1, and is at the level of 10% of the total suppression. Most published constraints on neutrino masses based on HALOFIT are not affected, as they rely on data probing the matter power spectrum in the linear or mildly non-linear regime. However, predictions for future galaxy, Lyman-alpha forest and weak lensing surveys extending to more non-linear scales will benefit from the improved approximation to the non-linear matter power spectrum we provide. Our approximation reproduces the induced neutrino suppression over the targeted scales and redshifts significantly better. We test its robustness with regard to changing cosmological parameters and a variety of modelling effects.


INTRODUCTION
A variety of solar, atmospheric, reactor and accelerator neutrino experiments have firmly established the existence of neutrino flavour oscillations, implying that neutrinos have at least two non-zero mass eigenstates (Becker-Szendy & et al. 1992;Fukuda & et al. 1998;Ahmed & et al. 2004). These experiments yield a lower limit on the sum of the neutrino masses; Mν > 0.05eV . Classical β-decay experiments provide an upper limit on the mass of the electron neutrino of 2.2 eV (Lobashev 2003), but may in future reach constraints as low as 0.2 eV (Wolf & KATRIN collaboration 2010).
Sensitivity to well below 1 eV is necessary to distinguish the neutrino mass hierarchy , which is of particular interest when formulating theories of particle physics beyond the standard model. Fortunately, such strong constraints can be obtained from the cosmic neutrino background. Simply requiring Ων ΩM ≈ 0.3 and E-mail: spb41@ast.cam.ac.uk † E-mail: viel@oats.inaf.it ‡ E-mail: haehnelt@ast.cam.ac.uk fixing the Hubble constant implies Mν 15 eV (Gershtein & Zel'Dovich 1966;Lesgourgues & Pastor 2006). The cosmic microwave background (CMB; e.g. Komatsu & et al. (2011)) combined with measurements of the matter power spectrum from galaxy surveys provide more stringent constraints, presently around Mν 0.6 eV (Elgaroy et al. 2002;Seljak et al. 2006;Tegmark et al. 2006;Gratton et al. 2008;Mantz et al. 2010;Thomas et al. 2010;Reid et al. 2010). Seljak et al. (2006) combined the Lyman-α forest flux power spectrum with the CMB to claim the tightest constraint so far; Mν < 0.17 eV. This strong limit may, however, be partly caused by a slight tension between the power spectrum amplitude preferred by the Lyman-α forest data and that preferred by the CMB. In any case, it should be tested and improved by upcoming data (Slosar & et al. 2011). The Lyman-α data alone provides a limit of Mν < 1 eV (Viel 2005;Viel et al. 2010), slightly stronger than that provided by the CMB data alone (Larson & et al. 2011).
The sensitivity of the matter power spectrum to neutrino mass is due to the free-streaming of the neutrinos which suppresses the growth of structure. On large scales, where matter perturbations can be described by linear the-ory, this effect is relatively straightforward to quantify with numerical simulations or approximate schemes (e.g. Bond et al. (1980);Zeldovich et al. (1980); Ma & Bertschinger (1994); Valdarnini et al. (1998)) using standard Boltzmann codes (Lewis et al. 2000;Lesgourgues & Tram 2011). Probing the effect of neutrinos on smaller, more non-linear scales is more challenging, but has been attempted by, for example, Saito et al. (2008, 2009), Wong (2008 and Lesgourgues et al. (2009), with extensions to perturbation theory. Also of interest here is the approach of Abazajian et al. (2005); Hannestad et al. (2005), who described the fully non-linear regime by perturbatively evaluating the effect of massive neutrinos on halo profiles.
However, accurately probing the effect of massive neutrinos on the growth of structure arbitrarily far into the nonlinear regime requires fitting formulae precisely calibrated from high dynamic range N-body simulations. The formula we present here builds on the widely-used HALOFIT (Smith et al. 2003). HALOFIT models the non-linear matter power spectrum by assuming that the growth of structure is due to a combination of the density structure of individual DM haloes and the clustering of these DM haloes relative to each other. It has been calibrated against a wide range of CDM simulations, superseding earlier prescriptions based on the stable clustering hypothesis (Hamilton et al. 1991;Peacock & Dodds 1996). At low and moderate redshifts HALOFIT models the non-linear growth of structure in a ΛCDM cosmology with an accuracy of 5 − 10% (Heitmann et al. 2010) for k < 1 hMpc −1 , but it has not been calibrated for non-zero neutrino masses.
Fortunately, the quality of the numerical simulations has advanced to the point where it is possible to model the impact of neutrinos to the level of a few percent, despite the difficulties presented by the large thermal dispersion of light neutrinos (Brandbyge et al. 2008;Viel et al. 2010). A quantitative numerical study of the suppression of the matter power spectrum due to free-streaming of neutrinos appears timely, as future matter power spectrum surveys aim to make measurements accurate to a few percent. Realising the potential of these measurements to constrain the neutrino mass requires tight control of systematic error from theoretical predictions.
We present results from simulations with the recently developed version of the Tree-SPH code GADGET -3 (Springel 2005;Viel et al. 2010). This suite of simulations was designed to evaluate the performance of HALOFIT predictions when compared to full numerical simulations. We use our simulations to calibrate an improved fitting formula, capable of accurately reproducing the effect of massive neutrinos on the non-linear matter power spectrum. We discuss our methods in Section 2, present our results in Section 3, and conclude in Section 4.

MODELLING THE NON-LINEAR MATTER POWER SPECTRUM
Here we briefly describe the numerical methods we use to simulate structure growth in the presence of massive neutrinos, as well as the approximations underlying HALOFIT . In Section 2.1 we review the predicted effect of neutrinos in linear theory. In Section 2.2 we discuss the strategies used for 10 -1 10 0 10 1 k /(h Mpc-1) 10 0 10 1 10 2 Matter Power Spectra at z =0,1 Figure 1. Matter power spectra from a simulation with massive neutrinos (dashed lines), and Mν = 0.6 eV, compared to the corresponding power spectra from a ΛCDM simulation (solid lines). Dotted lines show the linear theory prediction for massive neutrinos, while dot-dashed lines show the linear theory prediction for ΛCDM. Upper (green) lines show spectra at z = 0, while lower (black) lines show them at z = 1. The simulations shown are S60 (small scales) and L60 (large scales), whose parameters are described in Table 1. S60 shows only modes where k > 4 × 2π/ Box-size, so sample variance is small, and L60 is cut off at k = 0.4 hMpc −1 , for clarity. Note that the absolute power spectrum is slightly less converged with respect to box size than the ratio of the massive neutrino power spectrum with and without massive neutrinos.
simulating the effect of massive neutrinos, describing methods based on both particle and fourier-space representations of neutrinos. In Section 2.4, we summarise HALOFIT . Finally, Section 2.3 gives details of our numerical simulations.

Massive neutrinos in linear theory
Neutrinos in the very early Universe are relativistic and tightly coupled to the primordial plasma. They decouple while still relativistic and then redshift adiabatically. The transition to non-relativistic behaviour takes place when Mν = 3Tν , at a redshift of (1+z) = 2×10 3 (Mν /1 eV). However, it is some time later before relativistic corrections are negligible, especially in the high-velocity tail of the Fermi-Dirac distribution. Once non-relativistic, massive, weaklyinteracting neutrinos behave qualitatively as a species of warm/hot dark matter, suppressing fluctuations on scales smaller than their thermal free-streaming length.
This suppression has a distinctive shape in the linear regime. Neutrinos of mass Mν have energy density Ων , and make up a fraction, fν , of the dark matter, where Ων = Mν 93.14 h 2 eV , fν = Ων ΩM .
ΩM is the present day matter density and h is the Hubble constant. The residual comoving peculiar thermal velocity of the neutrinos at redshift z is By analogy with the Jeans length, we can define a co-moving c 2011 RAS, MNRAS 000, 1-12 free-streaming scale Neutrinos cannot cluster on scales smaller than kFS, and so the matter power spectrum, defined as where a tilde denotes a Fourier transformed quantity, is suppressed by a constant factor proportional to the neutrino energy density, fν . The Boltzmann equation can be solved numerically to show that the fractional change in the power spectrum for fν < 0.07 is δP P ≈ −8fν .
On larger scales, k < kFS, neutrinos have free-streamed only for a portion of the cosmic history, inducing less suppression. The largest scale on which neutrinos alter the matter power spectrum is given by the wavenumber for which they freestreamed only immediately after becoming non-relativistic (Lesgourgues & Pastor 2006) knr ∼ 0.018 ΩM Mν 1 eV hMpc −1 .
Today, neutrinos in the observationally acceptable range have knr ∼ 10 −3 , placing the free-streaming effect at the scales probed by galaxy and Lyman-α measurements of the matter power spectrum, as illustrated by Figure 1, which shows non-linear matter power spectra estimated from our simulations. For further information on the effect of neutrinos in linear theory see, e.g. Kolb & Turner (1990);Bond et al. (1980); Eisenstein & Hu (1999) and references therein.

Incorporating massive neutrinos into N-body simulations
The first N-body simulations incorporating massive neutrinos were aimed at evaluating their potential as a dark matter candidate (White et al. 1983;Klypin et al. 1993;Ma & Bertschinger 1994). Recently attention has shifted to their perturbative effect as a sub-dominant component in a mainly ΛCDM universe. We build on the work of Viel et al. (2010), who incorporated massive neutrinos into GADGET -3 (Springel 2005) in two ways; as a new particle species, and using a Fourier-space representation (Brandbyge & Hannestad 2009). Viel et al. (2010) have described the implementation in detail and performed extensive validation tests. Here we briefly review the two implementations of neutrinos. We refer the reader to Brandbyge et al. (2008); Brandbyge & Hannestad (2009) and Viel et al. (2010) for further details. Both the particle and the Fourier space implementations incorporate massive neutrinos directly into the simulations and so are significantly more accurate than simulations altering only the initial conditions of a ΛCDM simulation, as e.g. in Agarwal & Feldman (2011). All our simulations are inherently Newtonian, ignoring relativistic corrections. It is therefore important to set up the initial conditions of the simulation at sufficiently low redshift that the vast majority of neutrinos are non-relativistic. Otherwise, neutrinos will be treated as massive while still effectively massless, over-estimating their effect by a constant factor in the linear regime, at the highest redshifts.

Particle-based Neutrinos
Massive neutrinos can be incorporated into cold dark matter simulations as a separate low-mass collisionless particle species. Structure growth is primarily driven by the cold dark matter. The neutrino species has significant thermal dispersion, preventing it from clustering on small scales. This allows us to avoid following the small-scale evolution of the neutrino component, first by ignoring the shortrange tree force for the neutrinos and second by setting the timestep purely from the dark matter species, effectively relaxing the Courant conditions for the neutrino component. Brandbyge et al. (2008) and Viel et al. (2010) found that these approximations significantly speed up computation while having a negligible impact on results.
One potential problem with particle based simulations is thermal shot noise from the finite sampling of the neutrino distribution (Wang & White 2007), especially for the higher thermal velocities encountered at high redshift or for particularly low mass neutrinos. We deal with this by explicitly checking the effect of shot noise, and increasing the neutrino particle number should it be non-negligible.

Fourier-space Neutrinos
The effect of neutrinos on the dark matter can also be approximated by adding their gravitational force to that of the dark matter in Fourier space (Brandbyge & Hannestad 2009), as part of the particle mesh computation. The massive neutrino component is assumed to have the power spectrum predicted by linear theory, in our case computed using the Boltzmann code CAMB 1 (Lewis et al. 2000). Since the neutrinos are not followed directly, we can avoid shot noise and the need to deal with large thermal velocities when computing their effect on the overall growth of structure. However, because the neutrinos are evolved using linear theory, Fourier-space simulations neglect any non-linear evolution in the neutrino component itself, including back-reaction from non-linearities in the dark matter. As we will discuss in Section 3 these effects are important, and so we shall focus mainly on simulations involving particles.

The numerical simulations
We primarily used particle-based simulations, summarised in Table 1. We also ran several simulations using the Fourierspace implementation of neutrinos using a grid with 512 3 elements, but as these were not used in the final results we omit them from the table for clarity.
Our fiducial simulations have 512 3 dark matter particles in a box of 150 Mpc h −1 . We tested Mν = 0.15 eV, 0.3 eV and 0.6 eV, running 6 simulations for each value of Mν with different realisations of the initial conditions, to suppress cosmic variance. To check we had correctly sampled the large scale modes, to explore the linear regime and to compare directly to Viel et al. (2010), we ran also simulations with box sizes of 512 Mpc h −1 . To facilitate comparison with Brandbyge et al. (2008) and Viel et al. (2010), our fiducial cosmology had ΩM = 0.3, a spectral index ns = 1.0, the  Table 1. Summary of simulation parameters. σ 8 has been derived from the other parameters. For each simulation with massive neutrinos we ran a ΛCDM simulation with identical parameters, but without massive neutrinos, in order to compute their ratio. The value of σ 8 given is that for the simulation with massive neutrinos.
Hubble parameter h = 0.7, and amplitude of the primordial power spectrum As = 2.27 × 10 −9 . This implies that a simulation without neutrinos has σ8 = 0.87 at z = 0, but a simulation with Mν = 0.3 eV has σ8 = 0.8. We tested our formula using a variety of simulations with different cosmological parameters, as summarised in Table 1.
The initial redshifts of our simulations were zi = 99 for Mν = 0.6 eV, zi = 49 for Mν = 0.3 eV and zi = 24 for Mν = 0.15 eV, which we found to be the earliest redshift where the neutrinos were fully non-relativistic. Our initial conditions were generated using the Zel'dovich approximation, and the linear theory power spectrum was calculated by CAMB. We acknowledge that second order Lagrangian perturbation theory is in principle more accurate, but our simulations are started at sufficiently high redshift that the effect should be small. The transfer function used for our cold dark matter (CDM) particles was a weighted average of the linear theory transfer functions for CDM and baryons, to account for the slight difference between them (Somogyi & Smith 2010). Neutrinos and CDM had identical initial random phase information, to simulate adiabatic initial conditions.
The numerical convergence of our results and the impact of neutrino shot noise were checked with simulations S60ND, S15NU and S60NU. We further checked the effect of baryonic physics with simulation S60BA, which included 512 3 baryon particles. S60BA modelled radiative cooling and star formation using the default multiphase model of GAD-GET -3 (Springel & Hernquist 2003). Differences much larger than 10% have been found for k > 20 hMpc −1 between simulations implementing different radiative processes (e.g. metal cooling) or feedback recipes in ΛCDM models (see e.g. Rudd et al. (2008); Guillet et al. (2010); van Daalen et al. (2011)). Furthermore, the model of AGN feedback suggested by van Daalen et al. (2011) in order to reproduce observations of groups of galaxies produces a power spectrum suppression of roughly 10% at k = 1 hMpc −1 and 1% at k = 0.3 hMpc −1 . We did not include this effect; however, van Daalen et al. (2011) found it to be independent of cosmological parame-ters, suggesting that it is also independent of the neutrino mass, and thus should not affect our results.

The HALOFIT formula
HALOFIT (Smith et al. 2003) estimates non-linear growth from a given linear theory power spectrum, based on the principles of the halo model. These assert that the growth of halos should depend only on the local physics at the scale of the halo, and not the details of the pre-collapse material, nor the larger-scale distribution of matter. Thus growth deep in the non-linear regime becomes independent of the details of the linear power spectrum, depending only on the non-linear scale and the slope and curvature of the power spectrum. In the mildly non-linear regime, the non-linear correction is assumed to be a function of the instantaneous linear power spectrum, without explicit dependence on the historical growth rate. HALOFIT is accurate to 5 − 10% at z < 3 (Heitmann et al. 2010), suggesting that these approximations are reasonably accurate for cold dark matter.
HALOFIT realises the above by splitting the non-linear power spectrum into two components. A quasilinear term ∆ 2 Q , and a non-linear term ∆ 2 H . This means the dimensionless power spectrum can be written as, Both terms depend on the non-linear scale, kσ, defined as the wave-number where the perturbation amplitude, filtered by a Gaussian, becomes greater than unity and linear theory breaks down. The quasilinear term, ∆ 2 Q , is given as an enhancement to the linear power spectrum, and incorporates an exponential cut-off ensuring it is negligible on small scales. The form of the quasilinear term is, Here so that n is the effective spectral index of the power spectrum. β and α are functions of n and f (y) = y/4 + y 2 /8. We refer readers interested in the exact values of these functions to Appendix C of Smith et al. (2003). ∆ 2 L is the dimensionless linear theory power spectrum.
Non-linear growth, ∆ 2 H , dominates when y ≡ k/kσ 1. The form of ∆H used by HALOFIT is, the spectral curvature of the power spectrum. a, b, c, µ and ν are functions of n and C, and ν(n) is not related to neutrinos. f1,2,3 are functions of ΩM , the total matter density. For the values of these functions, we again refer to Appendix C of Smith et al. (2003). Note that the non-linear growth depends only on y, n, C and the matter fraction.
Massive neutrinos affect this prescription because they suppress the linear power spectrum, increasing the nonlinear scale and retarding growth. HALOFIT estimates the effect of this suppression on non-linear power, but neglects any effect of the neutrinos on the process of halo collapse, as well as the effects of non-linear growth on the neutrinos themselves. Figure 2 shows ratios between the matter power spectrum for massive neutrinos, and ΛCDM. Plotted are the predictions of linear theory, of simulations S60 and L60, and of grid-based simulations with identical parameters. The maximal power spectrum suppression in linear theory is proportional to −8fν , compared to a maximal non-linear suppression of approximately −10fν for z 1, and −10.5fν for z = 2. As discussed in Section 2.1, linear theory predicts an increasing suppression approaching an almost flat plateau on small scales. In the numerical simulations, the flat plateau is replaced by a minimum. The suppression increases up to k ∼ 1 hMpc −1 , after which it gradually decreases. This distinctive spoon shaped pre-virialization feature was also found by Brandbyge et al. (2008) and Viel et al. (2010), and in CDM simulations by Lokas et al. (1996);Smith et al. (2007).

RESULTS
We can understand this shape by considering the effect of neutrinos on non-linear growth. The suppression caused by massive neutrinos delays the onset of non-linear growth and increases the wave-number of the non-linear scale. Between the non-linear scale for the simulation without neutrinos and the non-linear scale for the simulation incorporating neutrinos, the suppression will be enhanced by the absence of non-linear growth. On the smallest scales, the non-linear growth lost as a consequence of the massive neutrinos becomes an ever-smaller fraction of the total, and the extra suppression decreases. Eventually non-linear growth comes to dominate the linear effect, and the suppression is less than predicted by linear theory.
As mentioned in Section 2.3, we performed checks for numerical convergence. A simulation with increased resolution, S60ND (see Table 1), had a relative matter power spectrum which differed from the fiducal simulation, S60, by 1 − 2% for k < 7 hMpc −1 . Furthermore, the relative matter power spectra for our two different box sizes, 512 and 150 Mpc h −1 , were always in good agreement on the largest scales probed by the smaller box, showing that our smaller box size was sufficient for an accurate estimate of the suppression of matter power spectrum due to the free-streaming of the neutrinos. The larger box showed a slightly increased suppression on small scales, probably due to limited resolution. We found that the effect of baryons was, at all redshifts, less than 1% at k < 8 hMpc −1 , gradually increasing on smaller scales, in agreement with the results of e.g. Jing et al. (2006); van Daalen et al. (2011).
Simulations S15NU and S60NU were run with a higher number of neutrino particles, to check explicitly the effect of shot noise. We found that, although the neutrino power spectrum is indeed shot noise dominated on small scales, the power on these scales is so small that their impact on the growth of the dark matter power spectrum is negligible. Hence shot noise has only a minor impact for simulations with a box size of 150 Mpc h −1 , at the level of 1% at the relevant scales and redshifts. Figure 2 shows that particle and Fourier-space implementations of neutrinos give similar results for the matter power spectrum, and both agree with linear theory on large scales. The particle implementation, however, shows a smaller suppression, which scales roughly linearly with Mν , particularly near the trough of the dip. This discrepancy has again a simple physical explanation; structure growth in the dark matter induces structure growth in the neutrino component, as the less energetic neutrinos fall into the gravitational wells created by the dark matter. This increases the power in the neutrinos, drawing them closer to the non-linear dark matter power spectrum and leading to a reduced neutrino suppression. Essentially this is a backreaction effect, where the dark matter drags the neutrinos with it; since it results from the non-linear growth, it is not being fully accounted for by linear theory neutrinos. Note that structure growth in the neutrino component itself remains linear; the effect is due to non-linear growth in the dark matter. Because this effect is important at the accuracy we wish to achieve, our main results (Figures 3, 4 and 5) are based on simulations with the particle implementation of neutrinos.   Table 1. For the smaller box sizes, more than one realisation of structure is available, and hence we show (dotted grey lines) the one σ error due to sample variance, as estimated numerically. Particle and Fourier-space methods use the same initial structure realisations. The Fourier-space method has more power in the largest scale mode, probably because the phase structure of the neutrino component has been prevented from evolving, enhancing the effect of sample variance.   HALOFIT clearly over-predicts the suppression of the matter power spectrum due to massive neutrinos in the nonlinear regime. The cause is similar to that discussed in Section 3.1; HALOFIT includes massive neutrinos only through the linear theory neutrinos suppression on the non-linear scale, and thus neglects any back-reaction from the dark matter. This interpretation is supported by the good agreement of HALOFIT with the Fourier-space simulations.

Comparison to the HALOFIT model
The largest discrepancy, around 10% of the total suppression, occurs at k ∼ 1 hMpc −1 . The location of the maximal suppression in the numerical simulation moves to larger scales at lower redshifts, an effect which is again not captured by HALOFIT , although it was present to some extent in our Fourier-space simulations. Furthermore, the amplitude of the suppression decreases with redshift, which is to be expected from our discussion in 3.1; as non-linear growth occurs in the dark matter more neutrinos fall into the gravity wells. A cross-over redshift occurs at z = 1; on very small scales the suppression here is equal to that of linear theory, while at lower redshifts it is less. Note also that in the quasilinear regime, 0.05 < k < 0.2 hMpc −1 , the simulations clearly agree much better with HALOFIT rather than linear theory The dependence of our results on Mν is well described by a linear relation, with the maximal suppression for a given redshift being proportional to fν , although the redshift dependence is more complicated.
HALOFIT over-predicts the effect of neutrinos on the smallest scales. This is not due to neutrino physics, but is a discrepancy induced because, as found by Hilbert et al. (2009), HALOFIT under-predicts the growth of non-linear power for k > 2 hMpc −1 in a ΛCDM universe by up to a factor of two. We corrected this by re-fitting the HALOFIT parameter that controls the small-scale power using our ΛCDM simulations, after which we could reproduce the asymptotic neutrino effect. We detail our modifications in Appendix A. Note that HALOFIT also has reduced accuracy at z = 3 − 4, with errors of 15 − 20%, compared to 5 − 10% at z = 0 − 2. Because of this, it fails to accurately predict the location of the peak non-linear suppression at z > 2. We did not correct this effect, as our attempts were found to negatively impact accuracy at lower redshifts.
We have also performed some simulations varying the cosmological parameters from our fiducial values, as described in Section 2.3. The results of these simulations were similar to those for our fiducial cosmology, and agree with the results discussed in Viel et al. (2010). The dependence of the non-linear suppression on cosmology is largely captured by HALOFIT , with the exception of a weak dependence on ΩM , which we include in our fitting formula.

Towards an improved fitting formula
In the previous Section, we compared our simulations to HALOFIT . While there was some discrepancy, it seemed to capture many of the important physical features. In light of this, we have designed an improved fitting formula by modifying HALOFIT to explicitly account for the effect of massive neutrinos. Our fitting formula was derived by performing a series of least-squares minimisations of L = ∆ 2 sim /∆ 2 fit − 1, where ∆ 2 sim is the matter power spectrum from the simulation and ∆ 2 fit is the equivalent from the fitting formula. Both spectra have been rebinned to be evenly spaced in log k. To find a minimum, we used Powell's method as implemented in SciPy 2 . To avoid any contamination from sample variance, we omitted scales with k < 6 2π Box , where Box is the box size. The redshift bins used in our fits were z = 0, 0.2, 0.5, 1, 2, 3. HALOFIT models the effect of non-linear growth at z = 4 rather poorly, and there is some residual shot noise present at these redshifts in our neutrino simulations. Furthermore, the only current probe of the matter power spectrum at these redshifts is the Lyman-α forest, which should be examined by modelling the flux power directly with full hydro simulations, as in Viel et al. (2010). Hence we do not attempt to provide a good fit for z > 3. We checked that our fit was not significantly altered by excluding the simulation results at z = 3.
Before addressing the effect of neutrinos, we attempted to correct the lack of power shown by HALOFIT on the smallest scales. Eq. 12 shows that the asymptotic behaviour of the non-linear term in HALOFIT is given by Thus, to correct the behaviour on small scales we altered γ(n), fitting it to the ΛCDM simulations labelled in Table 1 as S60, S60AS, S60NS, S60H and S30OL. We were interested only in the asymptotic behaviour, so we omitted the results from the large simulation boxes from the fit. We thereby included extremely small scales into the fit; down to k = 30 hMpc −1 . The numerical minimisation algorithm proved unstable when constrained to k < 6 hMpc −1 . We checked that the resulting fit is a good match to larger scales. The largest remaining contribution to the discrepancy between HALOFIT and our simulations is the decrease in the maximal neutrino suppression due to neutrinos falling into dark matter halos. This is proportional to fν and roughly constant on large scales. To account for it we modified the nonlinear term with the ansatz with l and m found by minimisation. We considered other powers of k in the denominator of Qν , as well as additional terms and powers of the non-linear scale kσ, but none of these significantly enhanced the fit.
To correctly match the location of the maximal suppression at higher redshift, z 2, we needed to alter the quasilinear term, reflecting a leakage of non-linear power from small scales to larger scales caused by the high thermal dispersion of the neutrinos. Therefore, we modified Eq. 9 as follows, where p, r and s were found again by minimisation. These coefficients, as well as l and m were found by simultaneously fitting to the simulations S15, S30, S60, L15, L30 and L60. For the larger box sizes, we used only scales with k < 0.5 hMpc −1 , and for the smaller boxes, we used k < 6 hMpc −1 . Finally, there was some residual dependence on ΩM . To account for this, we modified l as l → l −t(ΩM −0.3), and fit t separately, based on a minimisation with simulations S30, S30OM and S30OL. The numerical values of the coefficients are shown in Appendix A, together with a summary of the fitting formula. The reduced error in the fit can be defined as where N d is the number of data points used in the fit. We find E = 0.007 with our coefficients, compared to E = 0.04 with the unmodifed HALOFIT . A rough estimate of the error on our fitting formula, including realisation noise at the 1 − σ level, is given by where kσ(z) is again the non-linear scale. Differences between the fitting formula and the simulations change sign between z = 0 and z = 0.2, suggesting the dominant contribution is realisation noise. The quoted errors are valid for k < 7 hMpc −1 , and z 3. At z > 3 our fitting formula is somewhat less accurate. The maximal error over all scales is Emax = fν ; for Mν = 0.3, Emax ∼ 2%.
Although they were not included in the fit, we checked that our formula accurately reproduced simulations with different cosmologies, chosen to be at the bounds of the range allowed by current data. The fitting formula errors for S60AS, S60NS, and S60H were very similar to those found when compared to the fiducial cosmology, although there were one or two outliers caused by realisation noise. Checking against S30OM, which has ΩM = 0.25, produced somewhat larger errors at the highest redshift bin, with z = 3. The lower value of ΩM means that z = 3 is deeper into the quasilinear regime, where HALOFIT is less accurate. Figure 6 shows the performance of our improved fitting formula, in a similar format to that of Figure 3. It is clear that our improved formula captures the suppression due to neutrinos significantly more accurately than HALOFIT , correctly predicting the location and depth of the maximal suppression to within the quoted errors.

DISCUSSION AND CONCLUSIONS
We have performed a suite of detailed simulations of the matter power spectrum, incorporating the effects of massive neutrinos using both particle and Fourier-space methods, but focussing on simulations with neutrino particles. We have compared our results to those predicted by the widely used HALOFIT approximation to the non-linear matter power spectrum, and presented an improved fitting formula.
Observational constraints on the neutrino mass prior to this work have either used linear theory (Reid et al. 2010) or HALOFIT to estimate the effect of massive neutrinos on the matter power spectrum. Thomas et al. (2010), for example, used photometric redshifts at z = 0.45 − 0.65, on scales k 0.2 hMpc −1 . This is in the weakly non-linear regime, where we find HALOFIT to be accurate. Hence our results should not significantly change their published constraints and that of studies probing similar scales.
However, while HALOFIT reproduces the suppresion of the matter power spectrum due to the free-streaming of neutrinos well on some scales, in the fully non-linear regime it systematically over-predicts the level of the effect. Furthermore, HALOFIT fails to capture the redshift dependence of the amplitude and scale of the peak suppression. Our quantitative numerical study of the non-linear regime will thus certainly be of importance for future galaxy (Wang et al. 2005;Schlegel et al. 2009;Carbone et al. 2011), CMB lensing (Cooray 1999Kaplinghat et al. 2003), Lyman-α forest (Gratton et al. 2008) and weak lensing (Ichiki et al. 2009;Beaulieu & et al. 2010) surveys. Abazajian & et al. (2011) have published a review of the potential constraints on the neutrino mass from a variety of observational measurements, many of which reach the Mν ∼ 0.1 eV range. Realising these limits will require theoretical predictions of the power spectrum accurate to the percent level, significantly exceeding the accuracy of current prescriptions.
We have proposed an improved fitting formula which modifies HALOFIT to take account of the effects of massive neutrinos for k < 7 hMpc −1 at z 3. Errors are proportional to the amplitude of the suppression, and thus to fν ; for Mν = 0.3 eV, the largest error is roughly 2% of the matter power spectrum. The maximal non-linear suppression we find is around −10fν , significantly greater than the −8fν expected from linear theory. Hence constraints obtained using our formula should be tighter than predicted by forecasts using linear theory. Saito et al. (2008Saito et al. ( , 2009Saito et al. ( , 2011 and Wong (2008)   proposed an approximation to the non-linear matter power spectrum based on one-loop perturbation theory, while Lesgourgues et al. (2009) have modelled the effects of neutrinos with an approach based on a time renormalisation group flow. These approaches are to some extent complementary to our efforts; while these formulae are designed for the quasilinear regime, for k < 0.3 hMpc −1 or z > 1, our formula should be useful even in the fully non-linear regime, at the lowest redshifts.
Our fitting formula should be sufficiently accurate to make forecasts for the accuracy of neutrino mass measurements of future surveys. We should stress, however that we have concentrated here on the relative impact of the neutrinos; HALOFIT itself only predicts the non-linear power spectrum to an accuracy of 5 − 10% (Heitmann et al. 2010) with larger discrepancies at very small scales and high redshifts (z > 3). Detailed quantitative analysis of actual data reaching into the non-linear regime will thus still need careful calibration with accurate numerical simulations. These will have to be tailored to the redshift and scales relevant for the particular data set at hand as the simulations still struggle to provide the required dynamical range and the fitting formula are not yet accurate at the percent level over the full range of scales, cosmological parameters, neutrino masses and redshifts of interest.
We have implemented our fitting formula in the form of a series of patches to CAMB, which have been included in the publicly released version.

ACKNOWLEDGMENTS
Calculations for this paper were performed on the COS-MOS Consortium supercomputer within the DiRAC Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS and the University of Cambridge, as well as the Darwin Supercomputer of the University of Cambridge High Performance Computing Service (http://www.hpc.cam.ac.uk/), provided by Dell Inc. using Strategic Research Infrastructure Funding from the Higher Education Funding Council for England. We thank Volker Springel for giving us permission to use GADGET -3. SB and MH are supported by STFC. The HALOFIT non-linear term is given by ∆ 2 H (k) = ∆ 2 H (k) 1 + µ(n)/y + ν(n)/y 2 (A1) ∆ 2 H (k) = a(n, C)y 3f 1 (Ω) 1 + b(n, C)y f 2 (Ω) + [c(n, C)f3(Ω)y] 3−γ(n,C) .