# Improved Treatment of Cosmic Microwave Background Fluctuations Induced by a Late-decaying Massive Neutrino

## I Introduction

Anisotropies in the cosmic microwave background (CMB) contain an enormous amount of information about the universe. Data presently available has been used- to constrain from two up to eight cosmological parameters. With the promise of ever more precise measurements of these anisotropies, it has become possible to envision CMB fluctuations as a tool to go beyond this minimal set and constrain other areas of physics. Recent proposed constraints include limits on Brans-Dicke theories , constraints on time-variation in the fine-structure constant , tests of finite-temperature QED , and limits on various models for both stable and unstable  massive neutrinos. All of these additional constraints, with one exception, are based on the high-precision fluctuation spectra expected from the MAP and PLANCK satellites. The sole exception is reference , in which Lopez et al. pointed out that the radiation from a neutrino decaying into relativistic decay products could produce such a large integrated Sachs-Wolfe (ISW) effect, that a fairly large mass-lifetime range can be ruled out from current observations. Lopez et al. argued that a neutrino with a mass greater than 10 eV and a lifetime between $`10^{13}`$ and $`10^{17}`$ sec could be ruled out. (Although this calculation assumes nothing about the nature of the decay products other than that they are relativistic, this limit is most useful when applied to decay modes into “sterile” particles such as a light neutrino and a Majoron, since other, more restrictive limits apply to photon-producing decays). Hannestad  showed that the MAP and PLANCK experiments should produce an even larger excluded region in the neutrino mass-lifetime plane.

In this paper, we improve on a major approximation of references . In these papers, the relativistic decay products were simply added to the background neutrino energy density in the program CMBFAST . However, when the massive neutrinos decay, the spatial distribution of the decay products is determined by the distribution of the non-relativistic decaying particles; it is not identical to the distribution of the background massless neutrinos. In fact, the approach of references  violates energy-momentum conservation. Although in this approach energy and momentum are explicity conserved at zeroth order (the mean), the first order perturbations violate energy-momentum conservation. This may seem like a small effect, but it actually has significant consequences for the CMB fluctuation spectrum.

In the next section, we discuss the formalism for the Sachs-Wolfe effect in the presence of neutrinos decaying after recombination. In section III, we present our results, showing the effects of correctly incorporating the spatial distribution of the decay products, and provide a simple physical explanation of these effects. In section IV, we show how our revised calculation affects the excluded region in the neutrino mass-lifetime plane, and in section V we briefly summarize our conclusions. A comparison of our new results with current data leads to the excluded region $`m_h>100`$ eV, $`\tau >10^{12}`$ sec, although smaller masses can also be excluded for a smaller range of $`\tau `$.

## II The ISW Effect with an Unstable Neutrino: Formalism

To calculate the CMB fluctuations in the presence of a decaying massive neutrino, we first review the basic precepts of the pertinent linear perturbation theory. The perturbed homogeneous, isotropic FRW metric can be parametrized as

$$ds^2=a(\tau )^2\left[d\tau ^2(1+2\psi )d\stackrel{}{x}^2(12\varphi )\right],$$
(1)   
where $`a`$ is the scale factor normalized to unity today and $`\tau `$is the conformal time defined by $`d\tau =dt/a`$, $`t`$ being the proper time of a comoving observer. This particular gauge is referred to as the conformal Newtonian gauge because the behavior of the potentials ($`\varphi `$, $`\psi `$) is akin, loosely speaking, to that of the Newtonian potential. These potentials determine the large scale CMB behavior. In particular, the photon temperature perturbation decomposed into its Fourier and angular modes can be shown to be 

$$\mathrm{\Delta }_{\mathrm{}}(k)=_0^{\tau _0}𝑑\tau (\dot{\varphi }(k,\tau )+\dot{\psi }(k,\tau ))\mathrm{exp}(\kappa (\tau ))j_{\mathrm{}}(k\tau _0k\tau ),$$
(2)   
where the subscript ‘0’ refers throughout to the present time and $`\kappa `$ is the optical depth from the present to some conformal time $`\tau `$in the past. For the purpose of clarity, all the sources contributing to the anisotropy from inside the last scattering surface have been set to zero in Eq. 2. The effect of the sources contributing to the anisotropy between last scattering and the present (as given in Eq. 2) is called the Integrated Sachs-Wolfe (ISW) effect. The power in the $`\mathrm{}^{th}`$ multipole is normally defined as $`\mathrm{}(\mathrm{}+1)C_{\mathrm{}}`$ with 

$$C_{\mathrm{}}=(4\pi )^2_0^{\mathrm{}}𝑑kk^2|\mathrm{\Delta }_{\mathrm{}}(k)|^2.$$
(3)   
The sources mentioned in connection with the ISW effect can be varied. At any time, the modes that are important to the ISW effect correspond to those scales which are smaller than the sound horizon of the whole fluid (matter+radiation) at that time. For these modes, the potentials can decay if there is radiation pressure or if the universe expands rapidly. In models with no cosmological constant, the main contribution to the ISW effect comes from just after recombination (since radiation redshifts faster than matter). Inclusion of a cosmological constant leading to a rapid expansion of the universe late in its history would boost the power on larger scales (small $`\mathrm{}`$). Any other astrophysical process which contributes to the radiation content of the universe between last scattering and the present will lead to an increase in the total ISW effect. One such scenario is that of a massive particle decaying around or after last scattering. We will consider the case of a massive neutrino decaying non-relativistically into (effectively) massless particles. The details of the daughter particles turn out to be irrelevant.

To quantify the evolution of the massive neutrino density, we will consider the Boltzmann equation for its distribution. In a homogeneous and isotropic universe, the distribution of the collisionless massive neutrino decaying non-relativistically into two massless particles follows 

$$\frac{}{\tau }f_h^0(q_h,\tau )=\frac{a^2m_h}{t_dϵ_h}f_h^0(q_h,\tau ),$$
(4)   
where $`t_d`$ is the mean lifetime of the neutrino, and $`ϵ_h`$ and $`q_h`$ are the comoving energy and momentum: $`ϵ_h^2=q_h^2+m_h^2a^2`$, and a superscript ‘0’ will be used throughout to denote unperturbed quantities. We make the following simplifications throughout our treatment: (1) neglect inverse decays, (2) neglect spontaneous emission, (3) neglect Pauli blocking factor. The solution approaches the familiar $`\mathrm{exp}(t/t_d)`$ behavior as the neutrino becomes non-relativistic. The evolution equation for the energy density of the unstable neutrino is the integral of Eq. 4. It reads

$$\dot{\rho }_h^0+3\frac{\dot{a}}{a}(\rho _h^0+P_h^0)=\frac{am_hn_h^0}{t_d},$$
(5)   
where overdots represent differentiation with respect to conformal time. It should be noted (as it is important if the decay is not completely non-relativistic) that the right hand side contains the product of $`m_h`$ and $`n_h^0`$ (number density) and not $`\rho _h^0`$.

We now turn on the perturbations in the metric. Although the conformal Newtonian gauge is the most useful in which to understand the ISW effect, for computational purposes<sup>*</sup><sup>*</sup>*The main advantage is that CMBFAST  is written in synchronous gauge. we will define all our variables in the synchronous gauge. Thus, we will express the integrand in Eq. 2 in terms of perturbations in the synchronous gauge. The synchronous gauge has the property that the coordinate time and the proper time of a freely falling observer coincide. All the perturbations are in the spatial part of the metric ($`g_{ij}=a^2\delta _{ij}+a^2h_{ij}`$) in this gauge. The perturbation $`h_{ij}`$ can be Fourier transformed and broken up into its trace and a traceless part as 

$$h_{ij}(\stackrel{}{x},\tau )=d^3k\frac{\mathrm{exp}(i\stackrel{}{k}\stackrel{}{x})}{k^2}\left[h(\stackrel{}{k},\tau )k_ik_j+6\eta (\stackrel{}{k},\tau )\left(k_ik_j\frac{k^2}{3}\delta _{ij}\right)\right].$$
(6)   
Instead of working with the conjugate momentum in the perturbed space-time, we will use $`q_h`$ and $`ϵ_h`$ as defined above  and in keeping with that, we will write out the perturbed massive neutrino distribution as

$$f_h(\stackrel{}{x},\stackrel{}{q}_h,\tau )=f_h^0(q_h,\tau )\left[1+\mathrm{\Psi }_h(\stackrel{}{x},\stackrel{}{q}_h,\tau )\right].$$
(7)   
Due to the fact that the decay term is linear in $`f_h`$, the form of the equation for the evolution of $`\mathrm{\Psi }_h`$ is identical to that of the stable massive neutrino but with $`f_h^0`$ now given by Eq. 4. The stable massive neutrino case has been clearly worked out in Ref. .

The decay radiation rises exponentially from being negligible in the past to some maximum value at $`\tau \tau _d`$ and then drops off as $`a^4`$ like normal radiation. It is more informative therefore to follow the quantity $`r_{rd}=\rho _{rd}^0/\rho _\nu ^0`$ where ‘rd’ denotes the decay radiation and $`\rho _\nu ^0`$ is the cosmological density in a massless neutrino. The evolution equation for $`r_{rd}`$ is

$$\dot{r}_{rd}=\frac{m_hn_h^0}{\rho _\nu ^0}\frac{a}{t_d}.$$
(8)   
The treatment of the perturbations in the decay radiation will be analogous to that of the massless neutrino as worked out in Ref. . To evolve the perturbations in the decay radiation, we will integrate out the momentum dependence in the distribution function by defining (in Fourier space)

$$F_{rd}(\stackrel{}{k},\widehat{n},\tau )=\frac{𝑑qq^3f_{rd}^0(q,\tau )\mathrm{\Psi }_{rd}(\stackrel{}{k},q,\widehat{n},\tau )}{𝑑qq^3f_{rd}^0(q,\tau )}r_{rd},$$
(9)   
where $`\stackrel{}{q}=q\widehat{n}`$ and $`\mathrm{\Psi }_{rd}`$ is defined analagously to Eq. 7. The equation governing the evolution of $`F_{rd}`$ can be worked out to give

$`\dot{F}_{rd}+ik\mu F_{rd}+4\left({\displaystyle \frac{\dot{h}}{6}}+{\displaystyle \frac{\dot{h}+6\dot{\eta }}{3}}P_2(\mu )\right)r_{rd}`$ $`=`$ $`\dot{r}_{rd}\left(N_03iN_1P_1(\mu ){\displaystyle \frac{2}{3}}N_2P_2(\mu )+\mathrm{}\right),`$ (10)   
$`N_0(k,\tau )`$ $`=`$ $`{\displaystyle \frac{𝑑q_hq_h^2f_h^0(q_h,\tau )\mathrm{\Psi }_h(k,q_h,\tau )\left(1\frac{8}{3}\left(\frac{q_h}{am_h}\right)^2+\mathrm{}\right)}{𝑑q_hq_h^2f_h^0(q_h,\tau )}},`$ (11)   
where $`\mu =\widehat{k}\widehat{n}`$ and $`P_n(\mu )`$ are the Legendre polynomials of order $`n`$. The series of terms in these equations arises because the perturbed quantities depend on the direction of momentum and to get the contribution to a daughter particle with momentum $`\stackrel{}{q}`$, we need to integrate over all possible $`\stackrel{}{q}_h`$. Thus Eq. 10 depends on both $`\mu `$ and $`\stackrel{}{q}\stackrel{}{q}_h`$. The situation simplifies enormously for non-relativistic decays because each term $`N_p`$, which contributes to the $`p^{th}`$ multipole progressively, is of $`𝒪(q_h^p/a^pm_h^p)`$ or higher. In Eqs. 10 and 11, the series has been truncated by only keeping terms up to $`𝒪(q_h^2/a^2m_h^2)`$ in the integrand. Similar equations for the evolution of perturbations in the decay radiation can be found in references . Apart from $`N_0`$, the terms on the right-hand side of Eq. (10) are completely negligible for non-relativistic decays.

The use of Eq. 10 is our only difference from the treatment in ref. . In the latter paper, the relativistic decay products were simply added to the neutrino background in CMBFAST. This is equivalent to setting the right hand side of Eq. 10 to zero. Since the perturbations in the decay products are determined by the perturbations in both the metric and the decaying massive particles, they are correctly described by Eq. 10. Although this may seem like a minor difference, it produces very large effects, as we now show.

## III The ISW Effect with an Unstable Neutrino: Results

The formalism outlined above for the evolution of an unstable neutrino and its decay products was integrated into the CMBFAST code . We investigated a range of masses from 10 eV to $`10^4`$ eV and lifetimes from $`10^{12}`$ to $`10^{18}`$ seconds. The underlying cosmology was taken to be a standard ($`\mathrm{\Omega }=1`$) CDM model with $`h=0.5`$ (with $`H_0=100h`$ km sec<sup>-1</sup> Mpc<sup>-1</sup>); baryon density $`\mathrm{\Omega }_Bh^2=.02`$ and scale invariant isentropic initial conditions (the same model was used in ref. ). Our results are shown in Fig. 1 for several masses and lifetimes, along with the results obtained by simply adding the decay products to the relativistic background. As pointed out in Ref.  there is indeed an enhancement in the spectra at relatively large scales due to the ISW effect produced by the decaying neutrino. We will see in section IV that for many values of neutrino mass and lifetime, the spectrum produced is far from that observed today, and therefore a large region of parameter space is ruled out due to this effect.

The location of this ISW induced bump is determined by the lifetime of the neutrino. For lifetimes shorter than the age of the universe, inhomogeneities on scales $`k`$ project onto angular scales $`\mathrm{}k\tau _0`$ where $`\tau _0`$ is the conformal time today, and we assume a flat universe. The potentials vary in time (and hence cause the ISW effect) most significantly at the time of decays on scales of order the sound horizon: $`k_{sh}^23/(4\tau _d^2w)`$ where $`w=P/\rho `$. Therefore, the bump in the spectrum is produced at $`\mathrm{}k_{sh}\tau _0(\tau _0/\tau _d)(4w/3)^{1/2}`$. At these late times, the dominant contribution to $`w`$ comes from the decay radiation; hence $`w\mathrm{\Omega }_{rd}/3`$ where $`\mathrm{\Omega }_{rd}`$ is the fraction of critical density in decay radiation. Therefore, the ISW bump should be roughly at

$$\mathrm{}_{ISW}\frac{\tau _0}{\tau _d}\sqrt{\frac{9}{4\mathrm{\Omega }_{rd}(\tau _d)}}.$$
(12)   
For a matter dominated universe the conformal time and time are related as follows: $`\tau t^{1/3}`$. For a $`m_h=10`$ eV, $`t_d=10^{15}`$ sec neutrino, $`\mathrm{\Omega }_{rd}0.15`$ and $`\tau _0/\tau _d(4\times 10^{17}\mathrm{sec}/10^{15}\mathrm{sec})^{1/3}7.4`$. Therefore, in this case we expect $`\mathrm{}_{ISW}29`$. The actual peak occurs at a larger value of $`l`$, due to entropy fluctuations which decrease $`w`$, thereby increasing $`k_{sh}`$ and, finally, $`\mathrm{}_{ISW}`$.

Notice from Figure 1 that we find quantitative disagreement with the results of Lopez et al.  (dashed curves). The new results show that a more accurate treatment of the spatial distribution of the decay products produces a surprisingly large change in the CMB fluctuation spectrum compared to the results of reference . This difference is larger for smaller masses as can be seen in the figure.

At least for low masses, the most obvious difference between the the old and new spectra is the smaller size of the ISW effect for the new case. This difference has a physical explanation: by not properly treating the perturbations in the decay radiation we overestimate an important source of the potential decays that drive the ISW effect. To see this we first expand the Boltzmann equation for decay radiation perturbations, Eq. 9, in multipole moments, $`F_{rd}=_lF_{rd,l}P_l`$, to obtain the following hierarchy, shown here for $`l2`$:

$`\dot{\delta }_{rd}+{\displaystyle \frac{2}{3}}\left(\dot{h}+2\theta _{rd}\right)`$ $`=`$ $`{\displaystyle \frac{\dot{r}_{rd}}{r_{rd}}}\left(\delta _h\delta _{rd}\right),`$ (13)   
$`\dot{\theta }_{rd}k^2\left({\displaystyle \frac{\delta _{rd}}{4}}\sigma _{rd}\right)`$ $`=`$ $`{\displaystyle \frac{\dot{r}_{rd}}{r_{rd}}}\theta _{rd},`$ (14)   
$`\dot{\sigma }_{rd}{\displaystyle \frac{2}{15}}\left(2\theta _{rd}+\dot{h}+6\dot{\eta }\right)`$ $`=`$ $`{\displaystyle \frac{\dot{r}_{rd}}{r_{rd}}}\sigma _{rd},`$ (15)   
where $`\delta _{rd}=F_{rd,0}/r_{rd}`$, $`\theta _{rd}=3kF_{rd,1}/4r_{rd}`$ and $`\sigma _{rd}=F_{rd,2}/2r_{rd}`$. The treatment of Lopez et al.  is equivalent to neglecting the right hand sides of the equations above. This simplification breaks down near $`\tau \tau _d`$, where $`\dot{r}_{rd}/r_{rd}`$ is not negligible.

Neglecting the $`\dot{r}_{rd}/r_{rd}`$ terms in the Boltzmann equations for the decay radiation perturbations results in errors in the perturbations. Let us focus on $`\theta _{rd}`$, which turns out to be primarily responsible for the big difference. Consider Eq. 14 for modes above the horizon at $`\tau \tau _d`$, since for these modes the approximate treatment of Ref.  gives wrong results. For these modes the $`k^2`$ terms calculated in the approximate scheme can be shown (see Appendix) to be roughly similar to its exact value. Then the exact solution $`\theta _{rd}`$ is related to the approximate solution $`\theta _{rd}^a`$ by

$$\dot{\theta }_{rd}\dot{\theta }_{rd}^a\frac{\dot{r}_{rd}}{r_{rd}}\theta _{rd}.$$
(16)   
where the superscript here and in what follows denotes the solution to the set of equations 13-15 obtained in the approximate scheme by neglecting the feedback terms on the right hand side. The exact solution for $`\theta _{rd}`$ is therefore much smaller than the approximate one. Examples for several different modes are shown in Figure 2.

These large overestimates of $`\theta _{rd}`$ lead to correspondingly large overestimates of the ISW effect and are primarily responsible for the differences between our spectra and those generated in Ref. . The Appendix demonstrates precisely how the perturbations in the decay-produced radiation affect the potentials that govern the ISW effect, and how treating the decay products as identical to the massless neutrinos violates energy-momentum conservation. The bottom line is that the ISW effect depends significantly on the behavior of $`\theta _{rd}`$ and inaccuracies in it lead directly to inaccuracies in the $`C_l`$’s. Why does the approximation work better for higher mass neutrinos? The ISW effect is generated during times when the universe has appreciable radiation. For low-mass neutrinos whose decay radiation never dominates the energy density, the decay radiation redshifts away relative to the matter, and is only important near $`\tau \tau _d`$. Therefore, neglecting the $`\dot{r}_{rd}/r_{rd}`$ terms creates errors in the decay radiation perturbations at the crucial time when they are driving the ISW effect. If the neutrino is massive enough, then its decay products are important for a range of times with $`\tau \tau _d`$ when the approximation is good. So the approximate treatment works better for higher-mass neutrinos, like the $`m_h=10`$ keV, $`t_d=10^{12}`$ sec case.

There are other visible differences between the anisotropy spectra generated in reference  and our more accurate treatment. One difference, which exacerbates the rise in power at large scales, is a drop in the small-scale ISW effect. For modes which enter the horizon when there is significant radiation, the $`\delta _h`$ term in Eq. (13) is an important source term. This increases $`\delta _{rd}`$ relative to $`\delta _{rd}^a`$ and since $`\delta _{rd}`$ is a source for the evolution of $`\theta _{rd}`$, it implies that $`\theta _{rd}^a<\theta _{rd}`$. Thus there is a decrease in the ISW effect at small scales in the approximate scheme of ref . This is not visible for the 10 eV unstable neutrino (in Fig. 1) because of the comparatively large signature of the first peak, but it is readily apparent for the 10 keV neutrino because of the large ISW effect at small scales.

## IV Comparison with current CMB data

Since the detection of anisotropies in the CMB by COBE, there have been dozens of observations of anisotropies on a wide variety of angular scales (refs. -). We now use these observations to place more accurate limits on neutrino mass and lifetime.

In ref. , a very rough constraint was placed on decaying neutrino models: a model was excluded if the power at $`l=200`$ was greater than at $`l=10`$. As we have noted in the previous section, a more accurate treatment of the decaying neutrinos results in a much smaller distortion in the CMB spectrum for a certain range of neutrino masses and lifetimes. However, as we will see, consideration of all the data leads to constraints which are almost as stringent as the rough contours in ref. .

CMB experiments typically report an estimate of the band power

$$\widehat{C}_i=\frac{1}{4\pi }\frac{\underset{l}{}(2l+1)W_{i,l}C_l}{_lW_{i,l}/l}$$
(17)   
where $`W_{i,l}`$ is the window function which depends on beam size and chopping strategy of experiment $`i`$. Each of these comes with an error bar or, in the case of correlated measurements, an error matrix $`M^1`$. The naive way to constrain parameters in a theory then is to form

$$\chi ^2=\underset{i,i^{}}{}\left(\widehat{C}_iC_i(C_l)\right)M_{ii^{}}\left(\widehat{C}_i^{}C_i^{}(C_l)\right).$$
(18)   
Here we have explicitly written the dependence of $`C_i`$ on the theoretical $`C_l`$’s which in turn depend on the cosmological parameters. This naive statistic is useful only if the band power errors are Gaussian. In fact, the probability distribution is typically non-Gaussian, with a large tail at the high end and a sharp rise at the low end of the distribution. In recognition of this, and guided by some compelling theoretical arguments, Bond, Jaffe, and Knox  proposed forming an alternative statistic:

$$\chi ^2=\underset{i,i^{}}{}\left(\widehat{Z}_iZ_i(C_l)\right)M_{ii^{}}^Z\left(\widehat{Z}_i^{}Z_i^{}(C_l)\right)$$
(19)   
where

$$Z_i\mathrm{ln}\left(C_i+x_i\right)$$
(20)   
with $`x_i`$ an experiment dependent quantity, determined by the noise. The covariance matrix is now

$$M_{ij}^Z=\left(\widehat{C}_i+x_i\right)M_{ij}\left(\widehat{C}_i+x_i\right).$$
(21)   
Bond, Jaffe, and Knox  have tabulated and made available the relevant data from the experiments in refs. -. We use this information and formalismWe also account for calibration uncertainty in the manner set down in ref. . to constrain the mass and lifetime of unstable neutrinos.

The $`\chi ^2`$ in eq. 19 depends on the parameters of the cosmological model. In principle, it would be nice to allow as many parameters as possible to vary in addition to the mass and lifetime of the neutrino. This must be balanced against the constraints imposed by non-negligible time needed to run the modified version of CMBFASTThe modified version, accounting for decaying neutrinos, takes about ten times longer than the plain vanilla code.. Our strategy is to vary the mass and lifetime of the neutrino; the overall normalization of the $`C_l`$’s; the primordial spectral index (equal to one for Harrison-Zel’dovich fluctuations); and the calibration of each experiment. For the other cosmological parameters, we make “conservative” choices. That is, we choose values likely to make the power on small scales $`(l200)`$ as large as possible compared with the power on large scales. This acts against the effect of the decaying neutrino, which boosts up power on large scales, and therefore leads to more conservative limits. At each point in $`(m,\tau )`$ space, we use a Levenberg-Marquardt algorithm (see e.g. ) to find the values of normalization, spectral index, and calibration which minimize the $`\chi ^2`$ defined in Eq. 19. The contours in Fig. 3 show these best fit $`\chi ^2`$ in the $`(m,\tau )`$ plane.

Figure 3 shows the constraints on the neutrino mass and lifetime for a Hubble constant $`h=0.5`$ and $`\mathrm{\Omega }_Bh^2=.02`$ in a flat ($`\mathrm{\Omega }=1`$) matter dominated ($`\mathrm{\Omega }_\mathrm{\Lambda }=0`$) universe. The high baryon content is above the favored value of Tytler and Burles  and serves to raise the power on small scales. Masses greater than $`100`$ eV are ruled out for almost all lifetimes we have explored ($`\tau >10^{12}`$ sec). For lifetimes between $`10^{14}`$ and $`10^{15}`$ sec, masses as low as $`30`$ eV are excluded at the two-sigma level. These results are similar to those of ref. , but more reliable because of the improvements in the calculated spectra and the more careful treatment of the data.

We checked that the contours for a different set of $`(h,\mathrm{\Omega }_b)`$ were similar to the contours in Fig 3. Fig. 4 shows the results for a cosmological constant-dominated universe. Again, a sizable region is ruled out, reflecting the robustness of the constraint.

Hannestad  performed a similar calculation, using future CMB experiments to rule out decaying neutrino models, but he used the same approximation as in reference ; the decay products were added into the background neutrino density. We expect that his excluded-region contours for low masses should shrink since ISW effect is the main discriminator for these masses.

It has been noted that the decay products from a very massive neutrino could keep the universe substantially populated with radiation or even radiation-dominated for most of its history. The presence of radiation has the effect of stopping the growth of density perturbations, which in a matter-dominated universe would grow as $`\delta a`$. Since these density perturbations should (eventually) collapse into the structure we see today, it is clear that structure formation arguments can also provide constraints on the neutrino mass and lifetime. Very coarse constraints on the radiation density can placed by requiring that the scales relevant to structure formation are able to grow sufficiently (assuming of course, we know the initial perturbations), as is done in Ref. . In fact, for a scale-invariant initial spectrum, the structure formation arguments of Ref.  also rule out a region at the bottom-right of our excluded region. A more detailed analysis yields more stringent constraints . In light of this, it is important to understand that the constraints from CMB are most useful for low masses, i.e., for massive decaying neutrinos which do not affect the late-time growth of the density perturbations appreciably. Future experiments (MAP and PLANCK) have the potential to constrain neutrino masses as low as 1 eV and maybe even lower . In the end, CMB and large scale structure constraints on massive decaying neutrinos both overlap and complement each other.

## V Conclusions

Our results indicate that for calculations involving the effects of decaying particles on CMB fluctuations, exact conservation of energy-momentum (not just conservation of the mean energy-momentum) is crucial. When perturbations in the decay products are correctly treated as being determined by the perturbations in both the metric and the decaying massive particle, energy and momentum of the massive particle plus its decay products are conserved. The result is a much smaller change (when an unstable neutrino is added) in the CMB fluctuation spectrum than was noted in ref. . However, by using a comparison with current data, rather than a simple constraint on $`C_{200}/C_{10}`$, we have been able to obtain an excluded region only slightly less restrictive than that obtained in ref. . This excluded region will grow as more data becomes available, culminating potentially in very restrictive limits from MAP and PLANCK . Our results, of course, can be generalized to arbitrary decaying particles.

The CMB spectra used in this work were generated with a modified version of CMBFAST . We thank Lloyd Knox for providing the data used to generate the constraints in section IV. This work was supported by the DOE and the NASA grant NAG 5-7092 at Fermilab and by the DOE grant DE-FG02-91ER40690 at Ohio State.

## A Source of Decaying Potentials

Here we show that the ISW effect in the decaying neutrino model is primarily driven by the dipole of the decay-produced radiation, $`\theta _{rd}`$. The ISW effect is driven by time changes to the potentials (eq. 2), which are determined from Einstein’s equations. In synchronous gauge, the source of these time changes is

$$\dot{\varphi }+\dot{\psi }=t_1+t_2+t_3+t_4,$$
(A1)   
where   
$`t_1`$ $`=`$ $`\left[2+{\displaystyle \frac{3}{k^2}}\left({\displaystyle \frac{\dot{a}}{a}}\right)^2(5+3w)\right]\dot{\eta },`$    $`t_3`$ $`=`$ $`2{\displaystyle \frac{\dot{a}}{a}}\eta ,`$        $`t_2`$ $`=`$ $`{\displaystyle \frac{1}{2k^2}}\left({\displaystyle \frac{\dot{a}}{a}}\right)^2(5+3w)\dot{h},`$    $`t_4`$ $`=`$ $`{\displaystyle \frac{3}{k^2}}\left(2{\displaystyle \frac{\dot{a}}{a}}D_\sigma \dot{D}_\sigma \right).`$          
The quantity $`D_\sigma `$ is related to the anisotropic stress of the fluid: $`D_\sigma =(3/2)(\dot{a}/a)^2(1+w)\sigma `$, where $`w=P/\rho `$ is the equation of state of the universe. We will consider the behavior of superhorizon-scale perturbations, where $`k\tau 1`$. We assume that the neutrinos decay well into the matter-dominated phase of the universe, and that the decay radiation never dominates the energy density of the universe, but does come to dominate the standard radiation, i.e., photons and massless neutrinos. Then the equation of state takes a simple form near neutrino decay: $`w1/3\mathrm{\Omega }_{rd}`$. In addition, the total fluid perturbation sources $`\theta `$ and $`\sigma `$ are dominated by the decay radiation, so that we can write $`\theta 4w\theta _{rd}`$, and $`\sigma 4w\sigma _{rd}`$. These assumptions are well motivated for $`m_h=10`$ eV, $`t_d=10^{15}`$ sec neutrinos which decay well into the matter dominated era, with $`\mathrm{\Omega }_{rd}0.15`$ at decay.

We first examine the behavior of $`\dot{\varphi }+\dot{\psi }`$in the approximation where we neglect the $`\dot{r}_{rd}/r_{rd}`$ terms in the Boltzmann equations for the decay radiation perturbations, i.e., we treat the decay radiation as massless neutrinos as in Lopez, et al . We then consider the effect of relaxing the approximation and calculating the decay radiation perturbations correctly. We denote the use of the Lopez et al.  approximation in all quantities by the superscript-$`a`$.

The potentials do not decay in a completely matter-dominated universe; $`\dot{\varphi }+\dot{\psi }`$ is sourced by the decay radiation and is therefore first order in $`w`$. The term $`t_4^a`$ is directly related to $`\sigma `$ and so is of order $`w`$. The linearized Einstein equations imply that $`\dot{\eta }\theta w\theta _{rd}`$, so that $`t_1^a`$ is also of order $`w`$. However, $`t_2^a`$ and $`t_3^a`$ are each zeroth order in $`w`$, so their sum must cancel to lowest order. Using the linearized Einstein equations and the continuity equation we find that

$$t_2^a+t_3^a8\frac{w\eta }{\tau },$$
(A3)   
demonstrating the required cancellation. In our approximation, the decay radiation perturbations can be calculated from the Boltzmann equation for massless neutrinos, which admit analytic solutions for the superhorizon modes of $`\theta _{rd}^a`$ and $`\sigma _{rd}^a`$. Using these solutions, it can be shown that

$$t_1^a29\frac{w\eta }{\tau },t_4^a12\frac{w\eta }{\tau },$$
(A4)   
so that $`t_1^a`$, $`t_2^a+t_3^a`$ and $`t_4^a`$ each contribute roughly comparable amounts to $`\dot{\varphi }+\dot{\psi }`$. In calculating the effect of the approximation on $`\dot{\varphi }+\dot{\psi }`$we will therefore have to separately consider each term. The quantities $`D_\sigma `$ and $`\dot{\eta }`$ are very much affected by the approximation, since they directly depend on the decay radiation perturbations, and the error in the decay radiation perturbations is of order the quantities themselves. This implies that $`\delta \dot{\eta }|\dot{\eta }|`$ and $`\delta D_\sigma |D_\sigma |`$, where $`\delta x|xx^a|`$ is the absolute error in the variable $`x`$. The zeroth order quantities $`\eta `$ and $`\dot{h}`$ are much less affected by the approximation. For super-horizon modes, we do not expect $`\eta `$ to evolve much from its initial value, and so the error in it (determined by the error in $`\dot{\eta }`$) is naturally small. Following this line of reasoning, one can write

$`\delta \eta `$ $``$ $`\delta \dot{\eta }\tau k^2\tau ^2w\eta ,`$ (A5)   
$`\delta \dot{h}`$ $``$ $`k^2\tau \delta \eta k^4\tau ^3w\eta .`$ (A6)   
Using these relations we find that

$$\delta t_1\frac{w\eta }{\tau },\delta t_2(k\tau )^2\frac{w\eta }{\tau },\delta t_3(k\tau )^2\frac{w\eta }{\tau },\delta t_4\frac{w\eta }{\tau },$$
(A7)   
which makes it clear that for superhorizon modes, the errors in $`t_1`$ and $`t_4`$ dominate the error in $`\dot{\varphi }+\dot{\psi }`$. Numerically it is seen that the error in $`t_1`$ is the most important.

We can see why the error in $`t_1`$ is the most important in a simple way. The dominant source term (see eq. 14) for $`\dot{\theta }_{rd}^a`$ is $`\delta _{rd}^a`$. Now $`\dot{\delta }_{rd}`$ and $`\dot{\delta }_{rd}^a`$ differ only by the terms on the right-hand side of eq. 13 (since $`\dot{h}`$ is not much affected by the approximation and $`\theta _{rd}\dot{h}`$ for super-horizon modes). But the right-hand side of eq. 13 contains the difference of $`\delta _h`$ and $`\delta _{rd}`$, and hence the fractional error in $`\delta _{rd}`$ is expected to be much smaller relative to that in $`\theta _{rd}`$ or $`\sigma _{rd}`$. So we can assume that $`\delta _{rd}`$ and $`\delta _{rd}^a`$ are roughly the same for the purpose of estimating the errors in $`\theta _{rd}`$ and $`\sigma _{rd}`$ in the approximate scheme. Therefore, for super-horizon modes (at $`\tau \tau _d`$), we can write to a good approximation

$$\dot{\theta }_{rd}\dot{\theta }_{rd}^a\frac{\dot{r}_{rd}}{r_{rd}}\theta _{rd}\text{and}\dot{\sigma }_{rd}\dot{\sigma }_{rd}^a\frac{\dot{r}_{rd}}{r_{rd}}\sigma _{rd}.$$
(A8)   
From this equation, we can gauge that the fractional errors in $`\theta _{rd}`$ and $`\sigma _{rd}`$ are roughly the same, and close to $`\tau \dot{r}_{rd}/r_{rd}`$. But since the coefficient for $`t_4`$ in $`\dot{\varphi }+\dot{\psi }`$is much less than that for $`t_1`$ (see eq. A4), we expect that the error in $`t_1`$ dominates which implies that $`\delta (\dot{\varphi }+\dot{\psi })\dot{\eta }\theta _{rd}`$. From eq. A8 we have that $`\theta _{rd}^a>\theta _{rd}`$ for super-horizon modes at $`\tau \tau _d`$. Therefore for these modes,

$$\left|\dot{\varphi }^a+\dot{\psi }^a\right|>\left|\dot{\varphi }+\dot{\psi }\right|.$$
(A9)   
This is the reason for the dramatic rise in power at large scales when we neglect the decay terms in the decay radiation Boltzmann equations.

Merely adding the decay radiation to the massless neutrino background causes large errors in the ISW effect. However, there exists a method of calculation that yields good results without introducing a separate Boltzmann hierarchy for the decay radiation. Since this method might be helpful for other late time processes which affect the CMB spectrum, and since it demonstrates that we really have isolated the source of our disagreement with Lopez et al. , we present it here.

The fix can be accomplished by explicitly evolving $`\alpha =1/(2k^2)(\dot{h}+6\dot{\eta })`$, a quantity that contains the problematic $`\dot{\eta }`$, within CMBFAST in the following differential equation:

$$\dot{\alpha }+2\frac{\dot{a}}{a}\alpha =\eta \frac{9}{2k^2}\left(\frac{\dot{a}}{a}\right)^2(1+w)\sigma $$
(A10)   
The dominant quantity on the right-hand side is $`\eta `$, which is quite unaffected by the approximation (recall that $`\delta \eta k^2\tau ^2w\eta `$. In contrast, when $`\alpha `$ is set by the equation (the default in CMBFAST),

$$\alpha =\frac{a}{\dot{a}}\eta +\frac{3}{2k^2}\frac{\dot{a}}{a}\delta +\frac{9}{2k^4}\left(\frac{\dot{a}}{a}\right)^2(1+w)\theta ,$$
(A11)   
the $`\theta `$ term is important, and inaccuracies in the large scale behavior are generated. The condition that $`\dot{\alpha }`$ given by Eq. (A10) should match that obtained from Eq. (A11) is conservation of momentum for the massive neutrino and its decay products. In the approximate scheme, the unperturbed quantities for the decay radiation are calculated correctly while the perturbations in it are set equal to that of the massless neutrino. This violates the energy-momentum conservation conditions for the system of the massive neutrino plus its decay products, and this is the reason behind the fact that different combinations of the Einstein equations lead to different potential decay rates. It may be noted that the CMBFAST code  used to calculate the fluctuation spectrum implicitly assumes energy-momentum conservation. When this condition is violated, the code cannot produce internally consistent results. In Fig. 5, we have plotted the result of using Eq. A10 (in place of Eq. A11) with the approximate scheme and as expected, good agreement with the actual curves is obtained. This exercise clearly shows that it is important to check for energy-momentum conservation when using approximate methods to model any part of the energy-momentum tensor.