Effects of the neutrino mass splitting on the non-linear matter power spectrum

We have performed cosmological N-body simulations which include the effect of the masses of the individual neutrino species. The simulations were aimed at studying the effect of different neutrino hierarchies on the matter power spectrum. Compared to the linear theory predictions, we find that non-linearities enhance the effect of hierarchy on the matter power spectrum at mildly non-linear scales. The difference between the different hierarchies is about 0.5% for a sum of neutrino masses of 0.1eV. Albeit this is a small effect, it is potentially measurable from upcoming surveys. In combination with neutrinoless double-beta decay experiments, this opens up the possibility of using the sky to determine if neutrinos are Majorana or Dirac fermions.


INTRODUCTION
Neutrino oscillation experiments have demonstrated that neutrinos have non-zero mass. Measurements of the flavor changing oscillations have provided a difference in the squares of the masses of the lightest and heaviest mass eigenstates ∆m 2 (0.05eV) 2 , yielding therefore a lower limit on the total neutrino mass (Σ = m νi ). On-going and forthcoming ground based neutrino experiments are sensitive to neutrino flavor and to the nature of neutrino mass (Dirac or Majorana) but are only sensitive to the absolute mass scale for large masses. On the other hand, cosmological probes are blind to flavor but sensitive to the absolute neutrino mass scale and there has been recently significant progress in constraining the sum of neutrino masses from cosmological observations. Massive neutrinos affect the observed matter power spectrum: their free-streaming damps the smallscale power and modifies the shape of the matter power spectrum below the free-streaming length (Doroshkevich et al. 1980;Hu et al. 1998;Takada et al. 2006;Kiakotou et al. 2008). Up-to-date only upper limits have been obtained. However, the constraints are getting tighter -the current limits on the total mass are 0.3eV (e.g., Thomas et al. (2009); Reid et al. (2010a); Komatsu et al. (2011); Saito et al. (2011);Riemer-Sørensen et al. (2011);de Putter et al. (2012))-and closer to the experimental limits derived from accelerator, reactor, solar, and atmospheric neutrino oscillations (see the reviews by Lesgourgues & Pastor (2006); Gonzalez-Garcia & Maltoni (2008) and references therein). Detecting the effect of neutrino masses on cosmological structure and measuring the neutrino mass scale is well within the reach of upcoming cosmological surveys (e.g., Takada et al. (2006); Hannestad & Wong (2007); Kitching et al. (2008); LSST (2009); Hannestad (2010); Lahav et al. (2010); Reid et al. (2010a); Jimenez et al. (2010); Carbone et al. (2011) and references therein).
The neutrino mass splitting required to explain the oscillation results implies that for three neutrino species there are two possible hierarchies in the neutrino mass spectrum: normal hierarchy (NH) with two light states and one heavy one and a total mass Σ 0.05eV; inverted hierarchy (IH) with two heavy states and one light one with Σ 0.1eV. If the absolute mass scale is much higher than the mass splitting then the mass hierarchy does not matter and this case is referred to as degenerate mass spectrum. The degenerate hierarchy however is under pressure from observations (Thomas et al. 2009;Reid et al. 2010a;Gonzalez-Garcia et al. 2010); see To determine the possible hierarchy is very relevant as it can complement results from neutrinoless double-β decay to help determine the nature of the neutrino itself (Jimenez et al. 2010): is it its own anti-particle? that is, is it a Majorana fermion? As discussed in Bahcall et al. (2004) (see their Fig. 3 recalling that a light neutrino mass of 0.07eV corresponds to Σ ∼ 0.2 − 0.25eV), if the next generation of neutrinoless double-β decay experiments find a signal, then neutrinos are Majorana. If these experiments do not see a signal, it is important to discriminate if that is because the signal is below the detection threshold, or neutrinos are truly Dirac particles. Here is where cosmology enters (Jimenez et al. 2010): if Σ > 0.25eV, then the hierarchy is effectively degenerate and neutrinos are Dirac. The interesting region to determine the hierarchy is for 0.1eV < Σ < 0.25eV, where the absence of neutrinoless double-β decay indicates that neutrinos are Dirac only if the hierarchy is inverted. In Jimenez et al. (2010) we addressed the above question using linear theory to predict the effect of neutrinos on the matter power spectrum. Our conclusion was that an ambitious surveyà la stage IV dark energy task force (Albrecht et al. 2006), could in principle distinguish between the inverted and normal hierarchy if Σ < 0.15eV. How-ever, because the effect on the matter power spectrum of neutrinos extends into the non-linear regime, it is crucial to perform N-body simulations that include massive neutrinos to properly quantify the effect.
The main physical effect that distinguishes different hierarchies is the fact that neutrinos of different masses have different transition from relativistic to nonrelativistic thus influencing the shape of the matter power spectrum (Slosar 2006;de Bernardis et al. 2009). Note that neutrino oscillation experiments rule out large regions in the Σ − ∆ parameter space (see Fig. 1) and therefore it is worth investigating only the allowed region (gray swath in Fig. 1). Figure 1. Constraints on the mass splitting from neutrino oscillations (shaded regions) and total neutrino mass from cosmology (vertical dashed line) in the parameter space defined by the sum of neutrino masses Σ and the mass splitting parameter ∆ characterizing the hierarchy. The key region where it is interesting to determine ∆ is 0.1 eV < Σ < 0.3 eV. The triangles indicate the cases studied in this paper for normal and inverted hierarchies (at Σ = 0.1eV) and for the degenerate hierarchy (at Σ = 0.3eV). Plot adapted from Fig. 1 of Carbone et al. (2011).

NUMERICAL SIMULATION METHOD
In this Letter we study the effect of the neutrino mass splitting on the non-linear matter power spectrum. In particular, we are interested in the question of how and if the neutrino hierarchy leaves an imprint on the matter power spectrum in the non-linear regime. To address this question we run cosmological N-body simulations for several different neutrino mass configurations: A) 3 massless neutrinos; B) 3 neutrinos of equal mass (degenerate case, Σ = 0.3eV); C) 1 massless neutrino and 2 massive neutrinos (inverted hierarchy, Σ = 0.1eV, ∆ = −0.50); D) 2 light neutrinos and 1 heavy neutrino (normal hierarchy, Σ = 0.1eV, ∆ = 0.32).
So far there have been mainly two approaches to incorporating neutrinos into cosmological N-body simulations: sampling the neutrino density with particles just like in the case of cold dark matter (White et al. 1983;Klypin et al. 1993;Brandbyge et al. 2008;Viel et al. 2010) or evolving the neutrino density on a grid (with a fixed size) using linear theory (Brandbyge & Hannestad 2009;Viel et al. 2010). These two approaches were combined into a hybrid method by Brandbyge & Hannestad (2010). The particle-based approach has the advantage of being able to capture the non-linear evolution of the neutrinos and their effect on the cold matter components, while the linear evolution of the neutrino density on a grid is by construction only able to model the linear gravitational effect of the neutrinos on the non-linear matter distribution. In the particle-based approach, however, the neutrinos are always treated as non-relativistic particles, since standard cosmological N-body codes do not include relativistic effects. Another problem of this approach is that the finite number of particles used to sample the neutrino density generate shot noise. As the neutrinos have high thermal velocities, they move quickly away from their initial positions and give rise to a Poisson-like shot noise. Both these problems can be somewhat circumvented by starting the simulation at a sufficiently low redshift, when the neutrinos are practically non-relativistic and their input power spectrum is large compared to the shot noise. On the other hand, one does not want to start too late, when non-linearities have already become significant on the scales of interest. Hence, in order to keep the shot noise sub-dominant already at the initial redshift of the simulation (ideally this would be the redshift at which the neutrinos become effectively non-relativistic), a large number of neutrino tracers is required.
Studies by Brandbyge & Hannestad (2010) and Bird et al. (2012) have shown that in spite of the shortcomings mentioned above the particle-based approach is more accurate in computing the effect of massive neutrinos on the non-linear matter power spectrum than the grid-based approach. This is especially the case, when looking at relative quantities like the ratio of the power spectrum with and without massive neutrinos, since in this case systematic effects due to the rather late start of the simulation cancel out to a large extent. Hence, in this paper, we present results of particle-based simulations only.
We focus on scales in the range of 0.01 h Mpc −1 < k < 1 h Mpc −1 . These scales will be probed to high accuracy by future surveys and are much less affected by baryonic physics (which we neglect in this paper) than smaller scales. Hence, we simulate a volume of 600 h −1 Mpc 3 with a particle load of 1 billion cold matter tracers and 2 billion neutrino tracers.
We adopt a flat ΛCDM cosmology with cosmological parameters compatible with current observational constraints. The primordial curvature power spectrum is specified by the scalar amplitude ∆ 2 R = 2.45 × 10 −9 at the pivot scale k p = 0.002 Mpc −1 with a spectral index n s = 0.97. We keep the present-day total matter fraction the same for all neutrino models: Ω m = Ω CDM + Ω ν +Ω b = 0.27 with the baryon fraction and the neutrino fraction given by Ω b = 0.046 and Ω ν = Σ/(93.8eVh 2 ), where Σ is the sum of neutrino masses in units of eV. We choose the Hubble parameter to be h = 0.7. This choice of cosmological parameters yields a present-day linear mass variance in spheres of 8 h −1 Mpc of σ 8 ≈ 0.8.
For setting up the initial conditions and for assessing the N-body simulation results, we need accurate linear predictions for the transfer functions for each component of the universe. To this end, we use the linear Einstein-Boltzmann solver CAMB (Version Jan 2012) (Lewis et al. 2000), which computes the individual linear transfer functions for cold dark matter, baryons, neutrinos, and photons. As we simulate the different neu- trino species separately with different neutrino tracers, we modified CAMB to store for each neutrino species a separate transfer function. To reduce transients due to the late start of the simulation (Scoccimarro 1998), we implement second-order Lagrangian perturbation theory (2LPT) for the cold matter components (cold dark matter and baryons). We compute numerically the (k-dependent) linear growth rate from two CAMB outputs around the initial redshift. We approximate the second-order growth function D 2 by D 2 ≈ −3/7 D 2 1 where D 1 is linear growth function (Bouchet et al. 1995). At the initial redshift (z i = 9) the neutrino perturbations are still very much in the linear regime. Hence, for the neutrino particles, it suffices to use the Zel'dovich approximation to displace the particles from the initial grid points and to assign the gravitational velocities. Then thermal velocities drawn from the appropriate Fermi-Dirac distribution are added to the neutrino velocities. We neglect higher-order multipoles of the neutrino phase space distribution (see Ma & Bertschinger (1994) for a treatment of the full neutrino phase space). Tests with CAMB where these multipoles were set to zero at z = z i , and were then evolved further to low redshift, have shown that the effect of these initial conditions on the linear matter power spectrum is negligible at z ≤ 2 (de Putter 2012).
The simulations were carried out with Gadget-2 (Springel 2005). We modified Gadget-2 slightly to take into account the effect of the massive neutrinos (and the radiation) on the evolution of the scale factor a.
Using smaller test runs, we performed convergence tests on the neutrino time stepping. We found that a maximum time step max(∆ ln a) = 0.025 in the longrange particle-mesh force computation is sufficient for an accurate computation of the effect of the neutrinos on the matter power spectrum. Note that we disabled the Courant condition for the neutrino tracers. Hence, the time step for the long-range force is determined from the velocity dispersion of the cold matter tracers alone. This speeds up the computations significantly and leaves the non-linear matter power spectrum virtually unchanged. We set the softening length of the short-range force to 20 kpc/h for both the cold matter and neutrino component.
Varying the initial redshift and the number of neutrino tracers, we confirmed that the measured ratio of the nonlinear matter power spectrum with and without massive neutrinos is robust against the neutrino shot noise and the residual transients due to the late start (Brandbyge et al. 2008). We compute the total matter density contrast δ by assigning the particles to a 2048 3 grid using the cloud-incell scheme. In this process cold matter and neutrino particles are weighted by the fraction they contribute to Ω m . Using Fast Fourier transforms and averaging |δ k | 2 over spherical shells with a bin width of ∆k = 0.01 h/Mpc, we obtain the power spectrum P (k) = |δ k | 2 . Note that we only consider the non-relativistic neutrino species sampled by particles, i.e. the fluctuations in the radiation (relativistic neutrinos and photons) are not taken into account. At the redshifts of interest (z 2), however, radiation contributes a negligible amount to the total energy budget. One possible uncertainty in our simulation results is the shot noise coming from the cold matter particles. As the particles start off from a grid, the shot noise is sub-Poissonian and scale-dependent at high redshifts. However, for z ≤ 1 and k < 1 h Mpc −1 , even a Poisson shot noise is negligible. Hence, we do not attempt to correct for shot noise.

RESULTS
First, we compare the results of the degenerate case with Σ = 0.3eV with previous works. We consider the fractional difference in the total matter power spectrum of models with and without massive neutrinos. The advantage of this relative quantity is that the sample variance present in the N-body simulations cancels out almost completely. The suppression of the total matter power spectrum due to the three degenerate massive neutrinos is shown in the left panel of Fig. 2. The simulation results show that non-linearities enhance the effect on mildly non-linear scales. This behavior is anticipated by perturbation theory (Saito et al. 2008;Wong 2008;Lesgourgues et al. 2009). In contrast to the linear theory predictions (solid lines), in the non-linear case there is a maximum suppression, whose depth and position in k depend on redshift. These numerical results are in excellent agreement with Brandbyge et al. (2008) and Viel et al. (2010) and can be reasonably well modeled with HALOFIT (Smith et al. 2003), which was extended by . Fractional difference in the total matter power spectrum (CDM+baryons+massive neutrinos) of the inverted hierarchy and normal hierarchy run (the sum of neutrino masses is kept fixed and only the mass splitting is varied). Note that also in this case, non-linearities enhance the effect on mildly non-linear scales compared to linear theory predictions (solid lines). Bird et al. (2012) to model the effect of massive neutrinos (black dotted lines).
In the middle and right panel of Fig. 2, we show the power spectrum suppression for the normal and inverted hierarchy with Σ = 0.1eV. We observe the same qualitative behavior as before. In order to make the small difference between the two models visible, the fractional difference in the total matter power spectrum of the two cases is shown in Fig. 3. Similar to the case of massive vs. massless neutrinos, the non-linear evolution enhances the difference between the two models on mildly non-linear scales. On even smaller scales which are in the stable clustering regime, the strong non-linearities eventually overcome the initial differences between the models and the remaining effect decreases with k and may eventually drop below the linear theory prediction at large k.
In order to test if the simulation results are affected by the realization of the initial Gaussian random field, we compare them with the results of another set of simulations with a different random realization. We find that both sets of simulation give virtually the same predictions.
For completeness, we show the neutrino power spectra for the normal and inverted hierarchy in Fig. 4. On large scales (k 0.1 h/Mpc), the neutrino power spectra from the simulations agree well with the linear predictions. As expected, even at high redshift we observe a Poisson-like shot noise due to the large thermal velocities of the neutrinos. Although this shot noise is large in the neutrino power spectrum, it is much smaller when one considers the total matter density, since neutrinos contribute less than 1% to the total matter fraction. Additionally, on scales in the non-linear regime, the matter power spectrum is dominated by the non-linear power sourced from large-scale modes, for which shot noise is negligible.

DISCUSSION AND CONCLUSIONS
The main result of our numerical experiments is that non-linearities enhance the dependence of the power spectrum on the different neutrino hierarchies, thus making the observational signature more pronounced. We estimate that, if all other cosmological parameters are known (including the sum of neutrino masses Σ), the two hierarchies can be distinguished with confidence, as illustrated in Fig. 5 as function of the maximum k considered, making the effect potentially measurable. We have assumed an effective volume of 1 (Gpc/h) 3 at z = 0 (red lines) and 10 (Gpc/h) 3 at z = 1 (blue lines). 4 Whether degeneracies with other cosmological parameters and systematic effects (galaxy bias, baryonic physics, observational limitations etc.) will cancel the detectability of the effect remains to be explored and will be considered elsewhere (Wagner et al. 2012).
Our findings indicate that cosmology has the potential of determining the neutrino hierarchy in the interesting window Σ 0.1eV. Signal-to-noise estimates done using linear theory predictions indicated that if Σ happens to be in the window 0.15eV < Σ < 0.25eV, cosmology could not help determine the hierarchy and thus the nature of neutrino masses (Jimenez et al. 2010), leaving an important gap in our knowledge of neutrino properties. The fact that non-linearities enhance the effect compared to the linear prediction, will potentially enable cosmology to determine the hierarchy in a wider Σ range and possibly close the gap.
As an aside and already noted in Jimenez et al. (2010), cosmology is more sensitive to |∆| than to its sign: a measurement of |∆| in agreement with that predicted by oscillations experiments for the measured Σ would provide a convincing consistency check for the total neutrino mass constraint from cosmology.
However, it is very important to keep in mind that one needs theoretical predictions of the absolute non-linear power spectrum at least accurate to the 0.1% level to actually be able to make any sensible measurement of the neutrino hierarchy. In this Letter we presented numerical predictions only for the relative non-linear power spectrum suppression. This relative quantity is much more robust against numerical errors. Even without massive neutrinos, it is challenging to compute the non-linear power spectrum to sub-percent precision (e.g., Heitmann et al. 2010). On small scales (k 1 h/Mpc), baryon physics, which is strongly model dependent and computationally very intensive, is non-negligible and makes it very hard to achieve this accuracy in the foreseeable future. In addition, although it is in principle much less demanding to compute the linear power spectrum to high accuracy, different linear Einstein-Boltzmann codes (e.g. CAMB (Lewis et al. 2000) and CLASS (Blas et al. 2011)) still do not agree to 0.1% precision on the relevant scales.
Despite these very challenging and open problems, precision measurements of the large-scale structure of the universe remain an interesting avenue to determine the neutrino hierarchy.
We acknowledge the participation of Carlos Peña-Garay at an early stage of this project and thank him for many interesting discussions. CW is grateful for in-  Figure 4. Neutrino power spectra for the normal (left panel: light neutrino, middle panel: heavy neutrino) and inverted hierarchy (right panel) at several redshifts. The shot noise effect described in the text can be easily seen. In our approximation and for the chosen value of Σ, the normal hierarchy has effectively two non-zero neutrino masses (M and m in the notation of Eq. 1), while the inverted hierarchy has effectively only one massive eigenstate M (the light neutrino is massless).