Sound Speed and Viscosity of Semi-Relativistic Relic Neutrinos
PPrepared for submission to JCAP
Sound Speed and Viscosity ofSemi-Relativistic Relic Neutrinos
Lawrence Krauss a,b and Andrew J. Long c a Physics Department and School of Earth and Space Exploration, Arizona State Uni-versity, Tempe, Arizona 85287, USA. b Research School of Astronomy and Astrophysics, Mt. Stromlo Observatory, AustraliaNational University, Canberra, Australia 2611. c Kavli Institute for Cosmological Physics, The University of Chicago, Chicago, Illinois60637, USA.E-mail: [email protected], [email protected]
Abstract.
Generalized fluid equations, using sound speed c and viscosity c as effective parameters, provide a convenient phenomenological formalism for testingthe relic neutrino “null hypothesis,” i.e. that that neutrinos are relativistic and free-streaming prior to recombination. In this work, we relax the relativistic assumption andask “to what extent can the generalized fluid equations accommodate finite neutrinomass?” We consider both the mass of active neutrinos, which are largely still relativisticat recombination m /T ∼ .
2, and the effect of a semi-relativistic sterile component.While there is no one-to-one mapping between mass/mixing parameters and c and c , we demonstrate that the existence of a neutrino mass could induce a bias tomeasurements of c and c at the level of 0 . m /T ∼ − . Keywords: relic neutrinos, sterile neutrino, sound speed, viscosity, cosmic mi-crowave background a r X i v : . [ a s t r o - ph . C O ] J u l ontents Precision measurements of the cosmic microwave background (CMB) and large scalestructure (LSS) are providing a wealth of information about the early universe andits constituents. This information is particularly valuable in the neutrino sector wherea number of fundamental questions have yet to be answered: What is the absoluteneutrino mass scale? Are some neutrinos sterile? Do neutrinos self-interact througha long range force? The next-generation of CMB and LSS experiments will bringdramatic improvements in sensitivity and the promise of new insight into the physicsof neutrinos [1].To address the questions listed above in a model-independent way, it is customaryto use phenomenological parameters. These parameters are introduced “by hand” intothe equations of motion (Einstein or Boltzmann equations). They are not defined byany underlying fundamental parameters, such as Lagrangian couplings or masses.The most familiar phenomenological parameters are the effective number of neu-trino species N eff and the total neutrino mass (cid:80) m ν . Since the relic neutrinos aredecoupled at the time of recombination and structure formation, their effect on theCMB and LSS are only gravitational. Thus, the phenomenological parameters encodehow much the neutrinos contribute to the energy densities (see [1, 2] for notation) ρ rad = ρ γ + N eff (cid:18) (cid:19) / ρ γ and Ω ν h = (cid:80) m ν . . (1.1)Since N eff and (cid:80) m ν are not defined from fundamental parameters, there does notnecessarily exist a one-to-one mapping from any specific microphysical model onto theparameters ( N eff , (cid:80) m ν ). Rather, the phenomenological parameters are most useful– 1 –s a test of the “null hypothesis.” A combination of the concordance cosmology andStandard Model of particle physics predicts N eff = 3 .
046 and (cid:80) m ν = m + m + m > .
05 eV where m i are the three neutrino mass eigenvalues. Measurements compiled bythe Planck collaboration [3] ( Planck , TT + lensing + ext), N eff (cid:39) . ± .
40 and (cid:88) m ν < .
234 eV at 95% CL , (1.2)are consistent with the null hypothesis.Two additional phenomenological parameters affect the evolution of neutrino den-sity inhomogeneities. These are the effective sound speed c and viscosity c [4, 5].The effective sound speed sets the sound horizon, which in turn controls the growthof neutrino density perturbations, and the viscosity parameter leads to an anisotropicstress and the damping of neutrino density perturbations. (See Refs. [6] for a discussionof these effects on the CMB.) Once again, the phenomenological parameters providea model-independent formalism to test the null hypothesis: if the relic neutrinos arerelativistic and free-streaming then one expects c = 1 / c = 1 /
3. The Planckcollaboration furnishes the measurements [3] (
Planck , TT, TE, EE + lowP + BAO) c (cid:39) . ± . c (cid:39) . ± . , (1.3)which are consistent with the null hypothesis.As measurements of the four phenomenological parameters improve with the nextgeneration of CMB and LSS experiments, we must be mindful of any deviation fromthe null hypothesis, as this would indicate the presence of new physics. In order toprobe the nature of the new physics, we must understand how a specific microphysicalmodel maps onto the phenomenological parameters. For instance, many studies haveinvestigated how eV-scale sterile neutrinos (motivated in part by the short baseline andreactor anomalies [7, 8]) manifest themselves in the CMB and LSS (for one such recentpaper see Ref. [9]). This provides a mapping from the sterile mass and abundance to N eff and (cid:80) m ν . We seek to extend that correspondence to the perturbation parameters c and c .In Sec. 2 we study the formalism (generalized fluid equations) in which the phe-nomenological parameters c and c arise. While this formalism is convenient fortesting the null hypothesis, we will see that it cannot generally accommodate realisticdeviations from the null hypothesis. Specifically, if the neutrinos are assumed to befree-streaming but allowed to be semi-relativistic (such is the case for sterile neutrinos)then the fluid equations describing their evolution cannot be mapped onto the gener-alized fluid equations. In Sec. 3 we estimate the dependence of c and c on neutrinomass, and we calculate the predicted deviations from the null hypothesis, 1 / − c and 1 / − c , for a model of sterile neutrinos that saturates the Planck limits inEq. (1.2). We summarize our results and discuss directions for future work in Sec. 4.The main paper is accompanied by Appendix A, where we derive the fluid equations fora free-streaming species from the Boltzmann hierarchy. Appendix B contains formulasrelevant to a semi-relativistic Fermi-Dirac phase space distribution.– 2 – .1 0.2 0.5 1.0 2.0 5.0 10.0 20.00.000.050.100.150.200.250.30 Mass (cid:45) to (cid:45) Temperature Ratio: m (cid:144) T recombinationm (cid:61) (cid:206) (cid:72) (cid:76) eV recombinationm (cid:61) (cid:206) (cid:72) (cid:76) eV adi. sound speed: c adi2 equation of state: w 0.1 0.2 0.5 1.0 2.0 5.0 10.0 20.00.0010.010.11 Mass (cid:45) to (cid:45) Temperature Ratio: m (cid:144) T recombinationm (cid:61) (cid:206) (cid:72) (cid:76) eV recombinationm (cid:61) (cid:206) (cid:72) (cid:76) eV adi. sound speed: (cid:68) c adi2 (cid:61) (cid:144) (cid:45) c adi2 equation of state: (cid:68) w (cid:61) (cid:144) (cid:45) w Figure 1 . The equation of state w (solid) and squared adiabatic sound speed c (dashed)for a Fermi-Dirac distribution of relic neutrinos with mass m and temperature T . The shadedregions indicates the epoch of recombination, T rec ≈ . − . m = 0 .
08 eV and 0 . We are interested in the background of relic neutrinos at temperatures T (cid:46) ρ ( τ ) and pres-sure ¯ P ( τ ). The corresponding equation of state and adiabatic sound speed are w = ¯ P / ¯ ρ and c = ˙¯ P/ ˙¯ ρ , where the dot indicates differentiation with respect to conformal time τ . After decoupling the neutrino background maintains its Fermi-Dirac distributionwith temperature T . Using the notation established in Appendix A we calculate w and c in Appendix B and present the result in Fig. 1. For simplicity we assumethat the neutrino spectrum is in the degenerate regime, and the common neutrinomass is m ≈ (1 / (cid:80) m ν . Initially the neutrino temperature is high, m (cid:28) T , and w, c ≈ / w = 1 / − w and ∆ c = 1 / − c start to grow as the neutrinos becomesemi-relativistic. For a Fermi-Dirac distribution we find∆ c ≈ ∆ w ≈ π m T (2.1)for small m/T . The anomalously small prefactor, 5 / π (cid:39) .
02, invalidates the naivedimensional analysis prediction ∆ c ∼ m /T .During recombination, the photon temperature is T γ ≈ . − . (cid:15) q q /(cid:15) q /(cid:15) q /(cid:15) · · · δ Π ˜Π · · · θ ˜ θ · · · σ ˜ σ · · · χ · · · ... ... ... ... ... ... . . . Table 1 . The perturbation variables discussed in the text can be organized according themultipole moment (cid:96) of the phase space distribution function from which they were calcu-lated, and the factors of energy (cid:15) and momentum q that were included in the momentumintegral. Numerical factors of 1 /
3, etc., are not shown; see Eqs. (A.12) and (A.13) for detailedexpressions. neutrino temperature is smaller by a factor of (4 / / , which corresponds to T rec ≈ . − . m = 0 .
08 eV, which saturates thePlanck bound in Eq. (1.2), the deviations fall into the range 0 . (cid:46) ∆ c (cid:46)
2% atthe time of recombination. For a heavier eV-scale sterile neutrino, the deviation is(10 − w/w (cid:28) c /c (cid:28) w/w and ∆ c /c . Additionally it suggeststhat the effects of finite neutrino mass will be at most ∼ . m /T in magnitude.Let us now consider perturbations to the homogenous Fermi-Dirac distribution.The details of this calculation appear in Appendix A. Since the inhomogenous phasespace distribution function depends on both momentum q and position (or wavevec-tor k in Fourier space), it is convenient to organize the perturbations into multipolemoments with index (cid:96) . Each moment of the phase space distribution function can beintegrated over momentum q = | q | . It is possible to include additional factors of themomentum-to-energy ratio q/(cid:15) in the integrand. For the lowest order multipole mo-ments ( (cid:96) = 0 , ,
2) one obtains the energy density contrast δ ( k , τ ), energy flux θ ( k , τ ),and anisotropic stress σ ( k , τ ). Other combinations of (cid:96) and q/(cid:15) lead to a doubly-infinite tower of perturbation variables, shown in Table 1; see Eqs. (A.12) and (A.13)for detailed expressions. Specifically, Π( k , τ ) is the pressure perturbation. In the non-relativistic limit, T /m (cid:28)
1, the perturbation variables constructed from additionalfactors of q/(cid:15) are suppressed by powers of
T /m . In the relativistic limit,
T /m (cid:29) q/(cid:15) ≈ ≈ Π ≈ δ/ , ˜ θ ≈ θ , and ˜ σ ≈ σ . (2.2)The evolution of perturbations in a system of freely streaming particles is describedby the collisionless Boltzmann equation [10]. Upon performing the multipole expansiondescribed above, the Boltzmann equation yields a hierarchy of coupled first order– 4 –ifferential equations describing the evolution of each moment. See Appendix A fordetails of this calculation. Focusing on the first few multipole moments, we performthe momentum integrals to obtain˙ δ = − (1 + w ) (cid:16) θ + 12 ˙ h (cid:17) + 3 ˙ aa (cid:16) wδ − Π (cid:17) (2.3a)˙ θ = − aa (cid:16) − c (cid:17) θ + k Π1 + w − k σ (2.3b)˙ σ = 415 θ + (3 c ) 215 (cid:0) ˙ h + 6 ˙ η (cid:1) − kχ − aa (cid:16) − c (cid:17) σ + ˙ aa (cid:0) ˜ σ − σ (cid:1) + 415 (cid:0) ˜ θ − θ (cid:1) (2.3c)˙Π = − aa (cid:16) − w (cid:17) Π + ˙ aa (cid:0) ˜Π − Π (cid:1) −
13 (1 + w )˜ θ −
16 ˙ h (cid:0) w − ˜ w (cid:1) , (2.3d)which we call the collisionless fluid equations . We are working in the synchronousgauge where the metric perturbations are denoted as h ( k , τ ) and η ( k , τ ), and theirevolution is given by Einstein’s equations. The equations for ˙ δ and ˙ θ are the familiarcontinuity and Euler equations . Note that k = | k | and a ( τ ) is the FRW scale factor.The parameter ˜ w is the pseudo-equation of state, defined in the appendix. The equationfor ˙ σ depends on the next moment ( (cid:96) = 3) in the multipole expansion χ ( k , τ ). This isthe familiar result for the Boltzmann hierarchy: the evolution of lower-order multipolemoments depends on the higher-order moments. In addition, the equations for ˙ σ and˙Π also depend on the tilde’d variables ˜ θ , ˜ σ , and ˜Π. Consequently the equations shownin Eq. (2.3) do not form a closed system. However, we are only interested in comparingthe form of these equations with the generalized fluid equations below, and for thatpurpose we do not require the rest of the hierarchy.In the ultra-relativistic regime, m (cid:28) T , we can approximate w ≈ c ≈ / δ = − (cid:16) θ + 12 ˙ h (cid:17) (2.4a)˙ θ = 14 k δ − k σ (2.4b)˙ σ = 415 θ + 215 (cid:0) ˙ h + 6 ˙ η (cid:1) − kχ , (2.4c)and Π = δ/
3, which are the fluid equations for free-streaming, relativistic particles.A phenomenological generalization of the fluid equations was proposed in Ref. [4, 5].By introducing the sound speed and the viscosity parameters, c and c , one can write This name is something of an oxymoron. In a perfect fluid, collisions occur frequently and tend toisotropize the perturbations. This enforces a vanishing of the anisotropic stress σ and higher multipolemoments. One should view Eq. (2.3) as the analog of the fluid equations for a free-streaming species. These can also be derived from the conservation of stress-energy [10]. – 5 –11] ˙ δ = − (cid:16) θ + 12 ˙ h (cid:17) + 3 ˙ aa (cid:16) − c (cid:17) δ + 12 (cid:16) ˙ aa (cid:17) (cid:16) − c (cid:17) θk (2.5a)˙ θ = − aa (cid:16) − c (cid:17) θ + 14 (3 c ) k δ − k σ (2.5b)˙ σ = (3 c ) 415 θ + (3 c ) 215 (cid:0) ˙ h + 6 ˙ η (cid:1) − kχ , (2.5c)which we call the generalized fluid equations (GFE). The rest of the Boltzmann hier-archy, e.g. the equation for ˙ χ , is unmodified. Relativistic and free-streaming neutrinosobey the fluid equations in Eq. (2.4), which corresponds to the limit c = c = 1 / / i.e. that the relic neutrinos are relativistic and free-streaming. A numberof studies have investigated the effects of c and c on the cosmic microwave back-ground [6, 12–22], and recently the Planck collaboration reported the measurements inEq. (1.3) using a combination of CMB and BAO data. These measurements illustratethe utility of the generalized fluid equations for testing – and thus far confirming – thenull hypothesis of relativistic and free-streaming neutrinos.However, it is not clear the extent to which the GFE is able to capture specific mod-els when we relax the assumptions of relativistic free-streaming particles. For instance,it is often said that ( c , c ) = (1 / ,
0) corresponds to a relativistic perfect fluid, andtherefore this limit has been used to model the effect of neutrino self-interactions [23–28] (see also [29]). However, while c = 0 allows for solutions in which the anisotropicstress and higher moments vanish as in a perfect fluid, it also allows for solutions wherethey are nonzero and static, which is not the case for a perfect fluid. These criticismswere recently raised by Refs. [21, 28, 30].In this work, we consider the effect of finite neutrino mass either arising from theactive neutrinos themselves or a heavier sterile neutrino component. This problem hasbeen investigated recently in Ref. [6] by numerically solving the Boltzmann hierarchy,and it was found that there is no clear degeneracy between neutrino mass and the soundspeed parameters. Our goal is to develop an analytic understanding of this result whilealso deriving a parametric relationship between the parameters of the GFE and theneutrino mass.If we relax the assumption of relativistic neutrinos but maintain the assumption offree-streaming neutrinos, then the density perturbations satisfy the collisionless fluidequations of Eq. (2.3). Clearly it is not possible to put the GFE of Eq. (2.5) intothe form of Eq. (2.3) even with a judicious choice of the parameters ( c , c ); theequations have different structures. However, the neutrinos are still semi-relativisticat the time of recombination, see Fig. 1, and this observation motivates us to expand In comparing with Eqs. (2–4) of Ref. [11], note that q ν ( k , τ ) = 4 θ ( k , τ ) / (3 k ) and π ν ( k , τ ) =2 σ ( k , τ ) and F ν, ( k , τ ) = 2 χ ( k , τ ) in the massless limit. – 6 –round the relativistic limit. Using the results of Appendix B the equation of state,pseudo-equation of state, and sound speed are written as w = 13 − c , ˜ w = 13 − c , and c = 13 − ∆ c (2.6)where all of the perturbations are proportional to m /T and we have used Eq. (2.1).We can similarly expand the perturbation variables around Eq. (2.2) as Π = δ/ − ∆Π , ˜Π = δ/ − , ˜ θ = θ − ∆˜ θ , and ˜ σ = σ − ∆˜ σ (2.7)where the deviations are O ( m /T ). Making these replacements the collisionless fluidequations Eq. (2.3) become˙ δ = − (cid:16) θ + 12 ˙ h (cid:17) + 2∆ c (cid:16) θ + 12 ˙ h (cid:17) − aa (cid:16) c δ − ∆Π (cid:17) (2.8a)˙ θ = − aa ∆ c θ + 14 k δ − k σ + 316 k (cid:0) c δ − (cid:1) (2.8b)˙ σ = 415 θ + (cid:0) − c (cid:1) (cid:0) ˙ h + 6 ˙ η (cid:1) − kχ − aa ∆ c σ − ˙ aa ∆˜ σ −
415 ∆˜ θ . (2.8c)Here we keep only terms up to linear order in the deviations. In summary, a systemof free-streaming particle obeys the collisionless fluid equations of Eq. (2.3), and if theparticles are semi-relativistic these equations can be approximated as in Eq. (2.8).Now we seek to compare Eq. (2.8) with the generalized fluid equations of Eq. (2.5).To facilitate the comparison we difference the two sets of equations to obtain˙ δ : 3 ˙ aa (cid:16) − c + 2∆ c (cid:17) δ + 12 (cid:16) ˙ aa (cid:17) (cid:16) − c (cid:17) θk − c (cid:16) θ + 12 ˙ h (cid:17) − aa ∆Π(2.9a)˙ θ : − aa (cid:16) − c − ∆ c (cid:17) θ + 34 (cid:16) c − − ∆ c (cid:17) k δ + 34 k ∆Π (2.9b)˙ σ : (cid:16) c − (cid:17) θ + (cid:16) c −
13 + ∆ c (cid:17) (cid:0) ˙ h + 6 ˙ η (cid:1) + 3 ˙ aa ∆ c σ + ˙ aa ∆˜ σ + 415 ∆˜ θ . (2.9c)Evidently, there is no choice of c and c that brings the two expressions into thesame form, i.e. causes the three lines of Eq. (2.9) to vanish. The generalized fluidequations thus fail to capture even this minor deviation from the null hypothesis. While we have shown that there is no choice of c and c for which the generalizedfluid equations reduce to the collisionless fluid equations, nevertheless, it is reasonable Explicit calculation using Eqs. (A.12) and (A.13) reveals that δ/ − ˜Π ≈ δ/ − Π) to leadingorder in m /T . – 7 –o ask the following question. Suppose that the neutrinos have a small mass and aresemi-relativistic at the time of recombination. This affects the evolution of their densityperturbations according to Eq. (2.8) and ultimately impacts the CMB temperatureanisotropies. However, suppose one (naively) analyzes the observed CMB data usingthe generalized fluid equations, Eq. (2.5), which do not capture the physics of thesemi-relativistic neutrinos. How will the best fit parameters c and c depend on theneutrino mass?Inspecting Eq. (2.9), we ask what choice of the phenomenological sound speed andviscosity parameters would give the best agreement between the GFE and collisionlessfluid equations. Taking c = 1 / − ∆ c causes the gravitational source term tovanish from the equation for ˙ σ , and taking c = 1 / − ∆ c causes a number of otherterms to exactly or partially cancel. This observation suggests that as the neutrinosstart becoming semi-relativistic, the sound speed and viscosity will deviate from thenull hypothesis, ( c , c ) = (1 / , / c ≈ c and c ≈ c . (3.1)While the identification of c and c with the adiabatic sound speed is not rigorous,we propose here that it quantitatively reflects the correct parametric behavior andorder of magnitude of the effect.It is interesting to note that both c and c begin to deviate from 1 / c , c ) = (1 / ,
0) in which only c deviates from its value in the nullhypothesis.One additional comment is in order. Whereas c is temperature dependent, seeFig. 1, the phenomenological parameters c and c are assumed to be static. Thuswe should interpret Eq. (3.1) to mean that c and c are derived from a weightedtime average of c between the epoch of neutrino decoupling and recombination.Our analytic approximation does not determine which function will appear in thetime averaging. However, since c decreases monotonically from 1 /
3, any arbitrarilyweighted time average must satisfy∆ c , ∆ c ≤ ∆ c ( T rec ) , (3.2)where the deviation in the adiabatic sound speed is evaluated at the time of recombi-nation when the neutrino temperature was T rec (cid:39) . A similar identification was employed in the mixed dark matter scenario of Ref. [4]. – 8 –sing the analytic expression for ∆ c from Eq. (2.1) we estimate∆ c ≈ ∆ c (cid:46) . m T . (3.3)The sum of the relic neutrino masses is constrained as in Eq. (1.2) using Planck data.If the limit is saturated, the neutrinos are in the degenerate regime, and we can take m ∼ .
08 eV as a reference point. For this mass, the anticipated deviation in the soundspeed and viscosity parameters at the time of recombination are∆ c ≈ ∆ c (cid:46) . (cid:16) m .
08 eV (cid:17) (cid:18) T rec . (cid:19) − . (3.4)Comparing with Eq. (1.3), we see that the expected deviation is smaller than Planck’ssensitivity to c and c . If the sensitivity to c improves by an order of magni-tude, the estimate of Eq. (3.4) suggests that the effect of finite neutrino mass couldbecome relevant. In that case, a more detailed numerical analysis would be necessaryto determine actual constraints..Next we consider the possibility that the relic neutrino background contains a sub-dominant component of eV-scale sterile neutrinos. The fact that the neutrinos aresterile, i.e. not weakly interacting, will not actually be relevant for this discussion.Rather, it only matters that they are semi-relativistic and free-streaming at the timeof recombination. Once again we ask the question: suppose that the CMB sky gener-ated in this model is studied (naively) using the generalized fluid equations, which donot explicitly account for the sterile neutrino component. How will the best fit phe-nomenological parameters, c and c , depend on the sterile mass and abundance?To study this model, one writes down two sets of collisionless fluid equations witheach taking the form of Eq. (2.3) but labeled by subscripts “a” for active and “s”for sterile. This significantly complicates the analysis, but we now proceed to arguethat one can reduce the system to a single dynamical degree of freedom in the limitwhere both active and sterile neutrinos are relativistic. Since the neutrinos are free-streaming, they only influence the densities of other species ( e.g. , photons) throughtheir gravitational effect on the metric perturbations. Einstein’s equations, whichgovern the evolution of the metric perturbations, only depend on the diagonal linearcombinations, e.g. ¯ ρ a + ¯ ρ s and δρ a + δρ s (see Ref. [10] for complete expressions). Thus,as far as Einstein’s equations are concerned, we do not need to know the separateevolution of the active and sterile neutrino perturbation variables, but only their sumsare relevant: ¯ ρ ν = ¯ ρ a + ¯ ρ s , ¯ P ν = ¯ P a + ¯ P s , δρ ν = δρ a + δρ s (3.5)and so on for the other perturbation variables, θ ν , σ ν , etc. In this way, we can modelthe combined active and sterile neutrino background as a two-component fluid. The– 9 –orresponding adiabatic sound speed is given by c ,ν = ˙¯ P ν ˙¯ ρ ν = c , a (1 + w a ) ¯ ρ a + c , s (1 + w s ) ¯ ρ s (1 + w a ) ¯ ρ a + (1 + w s ) ¯ ρ s (3.6)where we have used Eq. (A.8). In the subsequent analysis we will assume, as above,that the effective sound speed and viscosity will follow the adiabatic sound speed asthe neutrinos become semi-relativistic.Note that while the sound speed formula bears a similarity to the baryon-photonfluid, the physics is very different. Before recombination the baryons and photons aretightly coupled due to frequent Thompson scattering [10]. Consequently, the baryonperturbation variables tend to track the photon perturbation variables, e.g. θ γ ≈ θ b and σ γ ≈ σ b ≈
0, and the single coupled fluid evolves as if it had an adiabatic soundspeed given by the analog of Eq. (3.6). In the case of free-streaming neutrinos, on theother hand, the active and sterile perturbations are not directly coupled. However, inthe relativistic regime, m a , m s (cid:28) T , the two sets of Boltzmann equations describingthe evolution of the active and sterile neutrinos are reduced to the same form, i.e. w a ≈ w s ≈ / c , a ≈ c , s ≈ /
3. If the isocurvature modes vanish initially, e.g. θ a ≈ θ s , then they remain vanishing as long at the both species are ultra-relativistic.Consequently the active and sterile neutrino perturbations evolve in the same way, eventhrough they are not directly coupled, and they can be modeled as a single fluid . Oncethe sterile neutrinos become non-relativistic, the isocurvature modes will grow, andthe two species will start evolving differently. Until that time, in the semi-relativisticregime, the sound speed given in Eq. (3.6) is appropriate.Further, we assume that the sterile neutrinos have a phase space distribution func-tion of the Fermi-Dirac form with the same temperature as the active neutrinos but adifferent overall normalization: f , a ( q ) = g (2 π ) e q/aT + 1 and f , s ( q ) = α f , a ( q ) . (3.7)The proportionality constant α controls the relative number densities, ¯ n s = α ¯ n a . Inthe relativistic limit this proportionality implies ¯ ρ s / ¯ ρ a ≈ α , and Eq. (1.1) gives∆ N eff ≈ (cid:18) (cid:19) / ¯ ρ s ρ γ ≈ α . (3.8)In the non-relativistic limit the proportionality implies∆ (cid:88) m ν ≈ m s ¯ n s ¯ n a ≈ α m s . (3.9)The effective number of neutrinos is measured with an error of δN eff ≈ .
4, see Eq. (1.2),which implies ¯ ρ s / ¯ ρ a = α (cid:46) δN eff / ≈ .
1. Similarly, imposing the bound on (cid:80) m ν implies m s (cid:46) (0 . /α ≈ α ≈ . One makes a similar reduction when modeling the Standard Model relic neutrino background asa single fluid, even though it is composed of three non-interacting components, corresponding to thethree neutrino mass eigenstates. – 10 – .1 0.2 0.5 1.0 2.0 5.0 10.0 20.010 (cid:45) (cid:45) (cid:45) Sterile Mass (cid:45) to (cid:45) Temperature Ratio: m s (cid:144) T S ound S peed : (cid:144) (cid:45) c ad i , Ν recombinationm (cid:61) (cid:206) (cid:72) (cid:76) eV recombinationm (cid:61) (cid:206) (cid:72) (cid:76) eV Ρ s (cid:144) Ρ a (cid:61) Ρ s (cid:144) Ρ a (cid:61) (cid:45) Ρ s (cid:144) Ρ a (cid:61) (cid:45) Ρ s (cid:144) Ρ a (cid:61) (cid:45) Figure 2 . The relic neutrino adiabatic sound speed from Eq. (3.6). Varying the sterileneutrino mass, m s , affects the equation of state and sound speed, w s and c , s , which appearin Eq. (3.6). The three lines correspond to different values of the sterile-to-active energyratios, ¯ ρ s / ¯ ρ a = 1 , − , − , and 10 − from top to bottom, which is a proxy for ∆ N eff / We evaluate the relic neutrino adiabatic sound speed using Eq. (3.6). The activeneutrinos are still relativistic at recombination, and we can set w a ≈ c , a ≈ /
3. Thesterile equation of state and sound speed, w s and c , s , are calculated from Eq. (3.7);they vary with the sterile neutrino mass m s as shown in Fig. 1. Figures 2 and 3 showhow the sound speed deviation ∆ c ,ν = 1 / − c ,ν varies with the sterile neutrino mass m s and relative abundance ¯ ρ s / ¯ ρ a . If the sterile neutrino is sufficiently light, then it isstill relativistic at recombination, and its effect on the sound speed is small. Similarly,if the relative sterile abundance is small, α = ¯ ρ s / ¯ ρ a (cid:28)
1, then it also has a suppressedimpact on the sound speed. As the sterile neutrino mass is increased, the sound speedbegins to deviate further from the null hypothesis value, c ,ν = 1 /
3. However, if themass is so large that the sterile neutrino is non-relativistic at recombination, then theapproximations used in our calculation are no longer valid, which is why we cut offFigures 2 at m s /T ≈
20. Fig. 3 also indicates the parameter space that is excluded bybounds on (cid:80) m ν and ∆ N eff as the red and blue shaded regions, respectively. Focusingon m s ≈ T ν (cid:39) . O (10 − ) before running into the bound on N eff . Generalized fluid equations (GFEs) provide a phenomenological formalism for testingthe relic neutrino “null hypothesis,” i.e. that the neutrinos are both relativistic andfree-streaming in the epoch prior to recombination. This formalism has two key advan-– 11 – .0 0.1 0.2 0.3 0.4 0.510 (cid:45) Sterile Mass: m s (cid:72) eV (cid:76) S t e r il e A bundan c e : Ρ s (cid:144) Ρ a (cid:68) c adi, Ν (cid:61) (cid:45) (cid:68) c adi, Ν (cid:61) (cid:45) (cid:68) c adi, Ν (cid:61) (cid:45) (cid:68) c adi, Ν (cid:61) (cid:45) (cid:68) c adi, Ν (cid:61) (cid:45) Figure 3 . The relic neutrino adiabatic sound speed from Eq. (3.6). The figure is calculated inthe same way as Fig. 2, but we hold T = 0 . c ,ν = 1 / − c ,ν are shown. The red shaded region is excluded by the bound on (cid:80) m ν , and the blue shadedregion is excluded by the bound on ∆ N eff . tages: it requires only minimal modifications to the fluid equations, which are easilyimplemented in a numerical Boltzmann solver; and these modifications are capturedby just two parameters, the effective sound speed c and the viscosity c , whichquantify deviations from the null hypothesis. However, as we have demonstrated inSec. 2, where we consider the effects of finite neutrino masses, specific microphysicalmodels that deviate from the null hypothesis cannot always be accommodated into theGFE formalism.One can nevertheless investigate how the presence of a finite neutrino mass wouldaffect the best fit values of c and c if the CMB sky were analyzed using a GFEanalysis. In Sec. 3 we propose that one can use the adiabatic sound speed at recom-bination, c ( T rec ), to gauge the magnitude of deviations in c and c from 1 /
3, thenull hypothesis prediction. Taking this as our measure, we estimate that a value ofthe neutrino mass saturating the Planck limit, m (cid:39) .
08 eV, could induce a deviationin the effective sound speed and viscosity by as much as 0 . c ≈ ∆ c ≈ − if Plank’s limits on ∆ N eff and (cid:80) m ν aresaturated. Since Planck’s error bars on the phenomenological parameters are relativelylarge, δc (cid:39) . δc (cid:39) .
037 respectively, the effect of finite neutrino massis currently imperceptible. However, if the next generation of CMB telescopes achievean order of magnitude improvement in sensitivity to the GFE parameters, then ourestimates suggest that the effect of finite neutrino mass cannot be neglected, and ana-lytical and phenomenological approximations will need to be supplemented by detailednumerical estimates. – 12 –n this regard we conclude by noting that our numerical estimates of Sec. 3 relyon the plausible argument that the best fit values of c and c will begin to deviatefrom 1 / c . If experimental sensitivities improve sufficientlyone could test this ansatz in detail using a numerical Boltzmann solver and the fol-lowing algorithm : for a particular neutrino mass solve the full Boltzmann hierarchy,Eq. (A.4), to generate realizations of the CMB sky; then for a particular ( c , c ) solvethe generalized fluid equations Eq. (2.5); using MCMC techniques find the values of( c , c ) that best fit the sky generated from the Boltzmann hierarchy. Acknowledgments
We are grateful for discussions with Sean Bryan, Daniel Grin, Wayne Hu, and MariusMillea, and we are grateful to Cecilia Lunardini for collaboration in the initial stages ofthis project. This work was supported at Arizona State University by the Departmentof Energy under Grant No. DE-SC0008016 and the National Science Foundation underGrant No. PHY-1205745. This work was supported at the University of Chicago bythe Kavli Institute for Cosmological Physics through grant NSF PHY-1125897 and anendowment from the Kavli Foundation and its founder Fred Kavli.
A Derivation of Fluid Equations
In this appendix we derive the fluid equations from the Boltzmann hierarchy. Ournotation mostly follows Ma & Bertschinger [10]: a ( τ ) is the scale factor, d τ = d t/a ( τ )is the conformal time, d x = d r /a ( τ ) is the comoving coordinate, k is the correspondingwave vector, q = a ( τ ) p is the comoving momentum, and (cid:15) ( τ ) = (cid:112) q + a ( τ ) m is a ( τ )times the proper energy measured by a comoving observer.The phase space distribution function is written as f ( k , q , τ ) = f ( q, τ ) (cid:0) k , q , τ ) (cid:1) . (A.1)For freely streaming particles, f satisfies the collisionless Boltzmann equation. For thehomogenous term this is simply ∂f ( q, τ ) /∂τ = 0, and the perturbations satisfy ∂ Ψ ∂τ + i qk(cid:15) (ˆ k · ˆ q )Ψ + d ln f d ln q (cid:32) ˙ η − ˙ h + 6 ˙ η k · ˆ q ) (cid:33) = 0 , (A.2) The algorithm we outline here is similar in spirit to the approach taken by Refs. [31, 32] to inferthe effect of neutrino mass on N eff . The analysis in this appendix does not assume any specific form for f ( q, τ ). – 13 –hich has been written here in synchronous gauge (with η ( k , τ ) and h ( k , τ ) the metricperturbations). The perturbation is decomposed onto the Legendre polynomials asΨ( k , q , τ ) = ∞ (cid:88) l =0 ( − i ) l (2 l + 1) Ψ l ( k , q, τ ) P l ( µ ) (A.3)with µ = ˆ k · ˆ q . The Boltzmann equation resolves to the set of coupled, first-orderdifferential equations˙Ψ = − qk(cid:15) Ψ + 16 d ln f d ln q ˙ h (A.4a)˙Ψ = qk (cid:15) (cid:16) Ψ − (cid:17) (A.4b)˙Ψ = qk (cid:15) (cid:16) − (cid:17) − d ln f d ln q (cid:0) ˙ h + 6 ˙ η (cid:1) (A.4c)˙Ψ l = 12 l + 1 qk(cid:15) (cid:16) l Ψ l − − ( l + 1)Ψ l +1 (cid:17) for l ≥ , (A.4d)which are collectively known as the Boltzmann hierarchy.If one is not interested in the momentum dependence of the perturbations, itwould seem that the problem is simplified by integrating Eq. (A.4) over q . In the caseof massless particles ( (cid:15) = q ) one can identify a new dynamical variable F l ( k , τ ) ∝ (cid:82) ∞ q d q qf ( q, τ )Ψ l ( k , q, τ ) for each original Ψ l , and in fact, the problem is simplified.However in the massive case ( (cid:15) (cid:54) = q ) the number of dynamical variables increases. Forinstance, Eq. (A.4a) gives the evolution of A = (cid:82) (cid:15) Ψ in terms of B = (cid:82) q Ψ (writtenschematically), but Eq. (A.4b) gives the evolution of B in terms of C = (cid:82) ( q /(cid:15) )Ψ .This second moment of Ψ requires its own evolution equation, and thus it is typicallyeasier to solve Eq. (A.4) directly. Nevertheless, the first few equations obtained byintegrating Eq. (A.4) correspond to the familiar fluid equations, and we now proceedto derive them.First we define the spatially averaged energy density, pressure, and pseudo-pressure:¯ ρ ( τ ) = 4 πa ( τ ) − (cid:90) ∞ q d q f ( q, τ ) (cid:15) ( q, τ ) (A.5)¯ P ( τ ) = 4 πa ( τ ) − (cid:90) ∞ q d q f ( q, τ ) q (cid:15) ( q, τ ) (A.6)˜ P ( τ ) = 4 πa ( τ ) − (cid:90) ∞ q d q f ( q, τ ) q (cid:15) ( q, τ ) . (A.7)For freely streaming particles ( ∂f /∂τ = 0), the energy density satisfies the homoge-nous continuity equation ˙¯ ρ ( τ ) = − aa (cid:0) ¯ ρ + ¯ P (cid:1) , (A.8)– 14 –nd the pressure satisfies ˙¯ P ( τ ) = − ˙ aa (cid:0) P − ˜ P (cid:1) . (A.9)We define the equation of state, pseudo-equation of state, and adiabatic sound speed, w ( τ ) = ¯ P ( τ )¯ ρ ( τ ) , ˜ w ( τ ) = ˜ P ( τ )¯ ρ ( τ ) , and c ( τ ) = ˙¯ P ( τ )˙¯ ρ ( τ ) . (A.10)They obey the useful relations˙ w w = 3 ˙ aa (cid:0) w − c (cid:1) and c = 5 w − ˜ w w ) . (A.11)Next we define a few of the lower order moments of Ψ l as δρ ( k , τ ) = 4 πa ( τ ) − (cid:90) ∞ q d q f ( q, τ )Ψ ( k , q, τ ) (cid:15) ( q, τ ) (A.12a) δP ( k , τ ) = 4 πa ( τ ) − (cid:90) ∞ q d q f ( q, τ )Ψ ( k , q, τ ) q (cid:15) ( q, τ ) (A.12b) δ ˜ P ( k , τ ) = 4 πa ( τ ) − (cid:90) ∞ q d q f ( q, τ )Ψ ( k , q, τ ) q (cid:15) ( q, τ ) (A.12c) δQ ( k , τ ) = 4 πa ( τ ) − (cid:90) ∞ q d q f ( q, τ )Ψ ( k , q, τ ) qk (A.12d) δ ˜ Q ( k , τ ) = 4 πa ( τ ) − (cid:90) ∞ q d q f ( q, τ )Ψ ( k , q, τ ) q k(cid:15) ( q, τ ) (A.12e) δS ( k , τ ) = 4 πa ( τ ) − (cid:90) ∞ q d q f ( q, τ )Ψ ( k , q, τ ) 2 q (cid:15) ( q, τ ) (A.12f) δ ˜ S ( k , τ ) = 4 πa ( τ ) − (cid:90) ∞ q d q f ( q, τ )Ψ ( k , q, τ ) 2 q (cid:15) ( q, τ ) (A.12g) δC ( k , τ ) = 4 πa ( τ ) − (cid:90) ∞ q d q f ( q, τ )Ψ ( k , q, τ ) 2 q (cid:15) ( q, τ ) . (A.12h)These correspond to perturbations in the energy density δρ , pressure δP , pseudo-pressure δ ˜ P , energy flux δQ , anisotropic stress δS , etc. We can also write δρ ( k , τ ) = ¯ ρ ( τ ) δ ( k , τ ) (A.13a) δP ( k , τ ) = ¯ ρ ( τ ) Π( k , τ ) (A.13b) δ ˜ P ( k , τ ) = ¯ ρ ( τ ) ˜Π( k , τ ) (A.13c) δQ ( k , τ ) = (cid:0) w ( τ ) (cid:1) ¯ ρ ( τ ) θ ( k , τ ) (A.13d) δ ˜ Q ( k , τ ) = (cid:0) w ( τ ) (cid:1) ¯ ρ ( τ ) ˜ θ ( k , τ ) (A.13e) δS ( k , τ ) = (cid:0) w ( τ ) (cid:1) ¯ ρ ( τ ) σ ( k , τ ) (A.13f) δ ˜ S ( k , τ ) = (cid:0) w ( τ ) (cid:1) ¯ ρ ( τ ) ˜ σ ( k , τ ) (A.13g) δC ( k , τ ) = (cid:0) w ( τ ) (cid:1) ¯ ρ ( τ ) χ ( k , τ ) , (A.13h) Here our notation diverges from that of Ma & Bertschinger [10]. – 15 –hich defines the dimensionless perturbation variables δ , Π, etc. For massless particles( (cid:15) = q ) we have w = ˜ w = c = 1 /
3, and the higher order moments simplify to thelower order ones, e.g. δ ˜ Q = δQ , ˜ θ = θ , ˜ P = P , and so on.Finally we are prepared to derive the fluid equations from the Boltzmann hierar-chy. Taking the appropriately-weighted momentum integral of Eq. (A.4a) leads to theinhomogenous continuity equation, which can be written in three equivalent forms:˙ δρ = − δQ −
12 ˙ h (cid:0) ¯ ρ + ¯ P (cid:1) − aa (cid:0) δρ + δP (cid:1) (A.14a)˙ δ = − (1 + w ) (cid:16) θ + 12 ˙ h (cid:17) + 3 ˙ aa (cid:16) wδ − Π (cid:17) (A.14b) (cid:18) δ w (cid:19) · = − (cid:16) θ + 12 ˙ h (cid:17) + 3 ˙ aa (cid:18) c δ w − Π1 + w (cid:19) . (A.14c)The relation Π = ( δP/δρ ) δ puts the second equation into a more familiar form. Usinga different weighting in the momentum integral yields,˙ δP = − aa δP + ˙ aa δ ˜ P − δ ˜ Q −
16 ˙ h (cid:0) P − ˜ P (cid:1) (A.15a)˙Π = − aa (cid:16) − w (cid:17) Π + ˙ aa (cid:0) ˜Π − Π (cid:1) −
13 (1 + w )˜ θ −
16 ˙ h (cid:0) w − ˜ w (cid:1) (A.15b) (cid:18) Π1 + w (cid:19) · = − aa (cid:16) − c (cid:17) Π1 + w + ˙ aa ˜Π − Π1 + w −
13 ˜ θ − c ˙ h , (A.15c)which gives the evolution of the momentum perturbation. Integrating Eq. (A.4b) leadsto the Euler equation, ˙ δQ = − aa δQ + k δP − k δS (A.16a)˙ θ = − aa (cid:16) − c (cid:17) θ + k Π1 + w − k σ , (A.16b)and integrating Eq. (A.4c) gives the shear equation,˙ δS = − aa δS + ˙ aa δ ˜ S + 415 δ ˜ Q + 215 ( ˙ h + 6 ˙ η ) (cid:0) P − ˜ P (cid:1) + 35 kδC (A.17a)˙ σ = − aa (cid:16) − c (cid:17) σ + ˙ aa (cid:0) ˜ σ − σ (cid:1) + 415 ˜ θ + 25 c ( ˙ h + 6 ˙ η ) − kχ . (A.17b)Equations (A.14)–(A.17) do not form a closed system, since the evolution of ˜Π, ˜ θ , ˜ σ ,and χ are undetermined. – 16 – Fermi-Dirac Distribution
For the relic neutrinos, which decoupled while they were relativistic, f ( q, τ ) maintainsthe Fermi-Dirac distribution f = g (2 π ) e q/aT + 1 (B.1)where aT = a T is independent of τ and g = 6 counts the two spin and three flavordegrees of freedom. The energy density, pressure, and pseudo-pressure are calculatedfrom Eqs. (A.5), (A.6), and (A.7) with (cid:15) = (cid:112) q + a ( τ ) m . In the limit m /T (cid:28) ρ ( τ ) ≈ gπ T + g m T (B.2)¯ P ( τ ) ≈ gπ T − g m T (B.3)˜ P ( τ ) ≈ gπ T − g m T (B.4)up to terms that are O ( m ). The equation of state, pseudo-equation of state, and adi-abatic sound speed are calculated using Eq. (A.10). The exact expressions, determinednumerically, are shown in Fig. 1. In the limit m /T (cid:28) w ≈ − π m T (B.5)˜ w ≈ − π m T (B.6) c ≈ − π m T , (B.7)up to terms of order O ( m /T ). These expressions give Eq. (2.1). References [1]
Topical Conveners: K.N. Abazajian, J.E. Carlstrom, A.T. Lee
Collaboration,K. N. Abazajian et al.,
Neutrino Physics from the Cosmic Microwave Background andLarge Scale Structure , Astropart. Phys. (2015) 66–80, [ arXiv:1309.5383 ].[2] G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti, and P. D. Serpico, Relicneutrino decoupling including flavor oscillations , Nucl. Phys.
B729 (2005) 221–234,[ hep-ph/0506164 ].[3]
Planck
Collaboration, P. A. R. Ade et al.,
Planck 2015 results. XIII. Cosmologicalparameters , arXiv:1502.0158 .[4] W. Hu, Structure formation with generalized dark matter , Astrophys. J. (1998)485–494, [ astro-ph/9801234 ]. – 17 –
5] W. Hu, D. J. Eisenstein, M. Tegmark, and M. J. White,
Observationally determiningthe properties of dark matter , Phys. Rev.
D59 (1999) 023512, [ astro-ph/9806362 ].[6] B. Audren et al.,
Robustness of cosmic neutrino background detection in the cosmicmicrowave background , JCAP (2015) 036, [ arXiv:1412.5948 ].[7] K. N. Abazajian et al.,
Light Sterile Neutrinos: A White Paper , arXiv:1204.5379 .[8] C. Giunti, Light Sterile Neutrinos , 2015. arXiv:1512.0475 .[9] M. Costanzi, B. Sartoris, M. Viel, and S. Borgani,
Neutrino constraints: whatlarge-scale structure and CMB data are telling us? , JCAP (2014), no. 10 081,[ arXiv:1407.8338 ].[10] C.-P. Ma and E. Bertschinger,
Cosmological perturbation theory in the synchronousand conformal Newtonian gauges , Astrophys. J. (1995) 7–25, [ astro-ph/9506072 ].[11] M. Archidiacono, E. Calabrese, and A. Melchiorri,
The Case for Dark Radiation ,Phys. Rev.
D84 (2011) 123008, [ arXiv:1109.2767 ].[12] R. Trotta and A. Melchiorri,
Indication for primordial anisotropies in the neutrinobackground from WMAP and SDSS , Phys. Rev. Lett. (2005) 011305,[ astro-ph/0412066 ].[13] F. De Bernardis, L. Pagano, P. Serra, A. Melchiorri, and A. Cooray, Anisotropies inthe Cosmic Neutrino Background after WMAP 5-year Data , JCAP (2008) 013,[ arXiv:0804.1925 ].[14] J. Lesgourgues and T. Tram,
The Cosmic Linear Anisotropy Solving System (CLASS)IV: efficient implementation of non-cold relics , JCAP (2011) 032,[ arXiv:1104.2935 ].[15] T. L. Smith, S. Das, and O. Zahn,
Constraints on neutrino and dark radiationinteractions using cosmological observations , Phys. Rev.
D85 (2012) 023001,[ arXiv:1105.3246 ].[16] R. Diamanti, E. Giusarma, O. Mena, M. Archidiacono, and A. Melchiorri,
DarkRadiation and interacting scenarios , Phys. Rev.
D87 (2013), no. 6 063509,[ arXiv:1212.6007 ].[17] M. Archidiacono, E. Giusarma, A. Melchiorri, and O. Mena,
Dark Radiation inextended cosmological scenarios , Phys. Rev.
D86 (2012) 043509, [ arXiv:1206.0109 ].[18] M. Gerbino, E. Di Valentino, and N. Said,
Neutrino Anisotropies after Planck , Phys.Rev.
D88 (2013), no. 6 063538, [ arXiv:1304.7400 ].[19] M. Archidiacono, E. Giusarma, S. Hannestad, and O. Mena,
Cosmic dark radiationand neutrinos , Adv. High Energy Phys. (2013) 191047, [ arXiv:1307.0637 ].[20] M. Archidiacono, E. Giusarma, A. Melchiorri, and O. Mena,
Neutrino and darkradiation properties in light of recent CMB observations , Phys. Rev.
D87 (2013),no. 10 103519, [ arXiv:1303.0143 ].[21] E. Sellentin and R. Durrer,
Detecting the cosmological neutrino background in theCMB , Phys. Rev.
D92 (2015), no. 6 063012, [ arXiv:1412.6427 ].[22] E. Di Valentino, S. Gariazzo, M. Gerbino, E. Giusarma, and O. Mena,
Dark Radiationand Inflationary Freedom after Planck 2015 , arXiv:1601.0755 . – 18 –
23] J. F. Beacom, N. F. Bell, and S. Dodelson,
Neutrinoless universe , Phys. Rev. Lett. (2004) 121302, [ astro-ph/0404585 ].[24] S. Hannestad, Structure formation with strongly interacting neutrinos - Implicationsfor the cosmological neutrino mass bound , JCAP (2005) 011,[ astro-ph/0411475 ].[25] N. F. Bell, E. Pierpaoli, and K. Sigurdson,
Cosmological signatures of interactingneutrinos , Phys. Rev.
D73 (2006) 063523, [ astro-ph/0511410 ].[26] R. F. Sawyer,
Bulk viscosity of a gas of neutrinos and coupled scalar particles, in theera of recombination , Phys. Rev.
D74 (2006) 043527, [ astro-ph/0601525 ].[27] A. Basboll, O. E. Bjaelde, S. Hannestad, and G. G. Raffelt,
Are cosmological neutrinosfree-streaming? , Phys. Rev.
D79 (2009) 043512, [ arXiv:0806.1735 ].[28] I. M. Oldengott, C. Rampf, and Y. Y. Y. Wong,
Boltzmann hierarchy for interactingneutrinos I: formalism , JCAP (2015), no. 04 016, [ arXiv:1409.1577 ].[29] M. Archidiacono and S. Hannestad,
Updated constraints on non-standard neutrinointeractions from Planck , JCAP (2014) 046, [ arXiv:1311.3873 ].[30] M. Shoji and E. Komatsu,
Massive Neutrinos in Cosmology: Analytic Solutions andFluid Approximation , Phys. Rev.
D81 (2010) 123516, [ arXiv:1003.0942 ]. [Erratum:Phys. Rev.D82,089901(2010)].[31] T. D. Jacques, L. M. Krauss, and C. Lunardini,
Additional Light Sterile Neutrinos andCosmology , Phys.Rev.
D87 (2013), no. 8 083515, [ arXiv:1301.3119 ].[32] E. Grohs, G. M. Fuller, C. T. Kishimoto, and M. W. Paris,
Effects of neutrino restmass on N eff and ionization equilibrium freeze-out , arXiv:1412.6875 ..