# 1 Introduction

## 1 Introduction

Direct experimental data permits the tau-neutrino to be quite heavy, $`m_{\nu _\tau }<18`$ MeV , and if the mass indeed lays in the MeV region, the $`\nu _\tau `$ would have a very strong impact on primordial nucleosynthesis. The recent data from Super Kamiokande , however, indicates possible oscillations between $`\nu _\mu `$ and either $`\nu _\tau `$ or a new sterile neutrino, $`\nu _s`$, with a large mixing angle, $`\mathrm{sin}2\theta 1`$ and a small mass difference $`\delta m^2=10^210^3\text{eV}^2`$. If oscillations proceed into $`\nu _\tau `$ then both $`\nu _\mu `$ and $`\nu _\tau `$ are light and the impact of a $`\nu _\tau `$ mass on nucleosynthesis is negligible. On the other hand it is not excluded that the oscillations go into a sterile $`\nu _s`$ and, if this is the case, $`\nu _\tau `$ may be heavy. Of course a heavy tau-neutrino should be unstable because otherwise it violates Gerstein-Zeldovich bound , and would overclose the universe.

If the life-time of $`\nu _\tau `$ is larger than a few hundred seconds, it can be considered as effectively stable on the nucleosynthesis time scale. The role of massive “stable” $`\nu _\tau `$ in the primordial nucleosynthesis was considered in several papers with chronologically improving accuracy. In pioneering work  and a subsequent paper  the assumptions of kinetic equilibrium for all participating particles and of validity of Maxwell-Boltzmann statistics were made (see also the later paper ), hence the problem was reduced to the solution of an ordinary differential equation of Riccatti type. However, nonequilibrium corrections to the spectra of both $`\nu _\tau `$ and massless neutrinos in the case of $`m_{\nu _\tau }`$ in the MeV range happened to be quite significant  and a more refined treatment of the problem had to be developed. In ref.  the simplifying assumption of Maxwell-Boltzmann statistics was dropped in favor of the exact Fermi-Dirac one, but it was assumed that kinetic equilibrium is maintained for all the species. Nonequilibrium corrections have been treated by one of the authors , in the update to , who found that they do not strongly change their original results. Exact numerical solutions of the full system of kinetic equations for all neutrino species without any simplifications have been done in refs. . In the last work a somewhat better precision was achieved and in particular a higher cut-off in particle momenta was used. Also the expressions for some matrix elements of the weak interaction reactions with Majorana particles were corrected. It has indeed been proven, that nonequilibrium effects are quite significant, almost up to 50%.

Consideration of the impact of an unstable $`\nu _\tau `$ on primordial nucleosynthesis also has a long history. In the earlier papers \[13-18\] some relevant effects were approximately estimated, but the accuracy of the calculations was typically rather low. In the next generation of papers \[19-23\] the level of calculations was improved considerably but is still subject to some criticism. General shortcomings of these calculations are the assumption of kinetic equilibrium for some of the species or the use of Maxwell-Boltzmann statistics. The recent paper  makes use of similar simplifications. However, the nonequilibrium corrections, as we have seen already for the stable case and will see below for unstable $`\nu _\tau `$, are quite significant. The corrections related to quantum statistics are typically at the level of 10% . The complete set of kinetic equations for the case of the decay $`\nu _\tau \nu _e+\varphi `$, where $`\varphi `$ is a massless scalar was numerically solved without any simplifying assumptions in ref.  for $`\nu _\tau `$ mass in the interval 0.1-1 MeV and life-times larger than 0.1 sec.

In this paper we extend the analysis of ref.  to a much larger mass interval, $`0.1<m_{\nu _\tau }/\text{MeV}<20`$, and life-times between 0.001 and 400 sec both for the decays of Majorana type tau-neutrino through the channels $`\nu _\tau \nu _e+\varphi `$ and $`\nu _\tau \nu _\mu +\varphi `$. Our numerical procedure is somewhat more accurate than that of ref. , we have a larger cut-off in dimensionless momentum $`y=pa`$, 20 instead of $`13`$ ($`a`$ is the expansion parameter of the universe), and we correct some expressions for the matrix elements squared, as in our previous paper related to the stable $`\nu _\tau `$. As in the preceding papers on decaying $`\nu _\tau `$ we neglect all other interactions of the scalar $`\varphi `$, except for the decay and inverse decay. The exact treatment of the scattering $`\varphi \nu \varphi \nu `$ and the annihilation $`\overline{\nu }\nu 2\varphi `$ presents considerable technical difficulties because of a larger number of contributions into the collision integral and because of the non-polynomial form of the matrix elements squared of these processes. It will be considered in the subsequent publication .

The paper is organized as follows. After an overview in section 2 we describe in detail our numerical solution in section 3. The results are presented in section 4 followed by a general discussion and conclusion in section 5. Throughout the paper the natural system of units, $`\mathrm{}=c=k=1`$, is used.

## 2 Decays of tau neutrino

We assume that the $`\nu _\tau `$ is a Majorana type fermion which is coupled to a scalar boson $`\varphi `$, possibly a Majoron , which is light or even massless. The coupling of $`\varphi `$ to neutrinos may have diagonal terms as e.g. $`g_1\overline{\nu }_\tau \nu _\tau \varphi `$ which are important for elastic scattering $`\nu _\tau +\varphi \nu _\tau +\varphi `$ and annihilation $`\overline{\nu }_\tau +\nu _\tau 2\varphi `$. The non-diagonal coupling $`g_a\overline{\nu }_\tau \nu _a\varphi `$ is responsible for the decay of $`\nu _\tau `$ into lighter neutrinos, $`\nu _e`$ or $`\nu _\mu `$ (correspondingly $`a=e`$ or $`\mu `$). We assume that one of these two couplings dominates, i.e. that $`\nu _\tau `$ predominantly decays either into $`\nu _e\varphi `$ or $`\nu _\mu \varphi `$ and consider these two possibilities separately. It is implicitly assumed that both $`\nu _e`$ and $`\nu _\mu `$ are the usual active neutrinos. Since chirality is changed by the coupling to a scalar field, the corresponding light neutrinos should also be Majorana particles, otherwise new sterile states would be produced by the decay. We also assume, that $`\varphi `$ is a weak singlet, because the LEP measurements  of the total decay width of $`Z^0`$ do not leave room for any other light weakly interacting particles except for the already known ones.

There are several possible ways of production of $`\varphi `$ in the primeval plasma. The first and evident one is by the decay $`\nu _\tau \varphi +\nu _a`$. Another possibility is the annihilation $`\nu _\tau +\nu _\tau \varphi +\varphi `$ and the third one is a possible non-thermal production in the course of a phase transition similar to the production of axions at the QCD phase transition. We neglect the last possibility, assuming that even if (pseudo)goldstone bosons were created in the course of the phase transition, the phase transition took place early enough so that the created bosons were diluted by a subsequent entropy release in the course of the universe cooling down. The rate of $`\varphi `$-production in $`\nu _\tau `$-annihilation can be estimated as:

$`{\displaystyle \frac{\dot{n}_\varphi }{n_{\nu _\tau }}}=\sigma _{ann}vn_{\nu _\tau },`$ (1)   
where $`v`$ is the relative velocity and $`\sigma _{ann}`$ is the annihilation cross-section. In the limit of large energies, $`s=4E_{cm}^2m_{\nu _\tau }^2`$ it is equal to: $`\sigma _{ann}(g_1^4/32\pi s)\mathrm{ln}(s/m_{\nu _\tau }^2)`$ (see e.g. ). One can check that this rate is small in comparison with the universe expansion rate $`H=\dot{a}/a`$, if $`g_1<10^5`$. In this case the production of Majorons through annihilation can be neglected and they would dominantly be produced by the decay of $`\nu _\tau `$. The opposite case of dominant production of $`\varphi `$’s by $`\nu _\tau `$-annihilation and their influence on nucleosynthesis was approximately considered in ref. .

The life-time of $`\nu _\tau `$ with respect to the decay into massless particles $`\varphi `$ and $`\nu _a`$ is equal to:

$`\tau _{\nu _\tau }={\displaystyle \frac{8\pi }{g_a^2m_{\nu _\tau }}}.`$ (2)   
This is faster than the universe expansion rate if $`T<410^{10}g_a\sqrt{m}_{\nu _\tau }`$, where the temperature $`T`$ and $`m_{\nu _\tau }`$ are expressed in MeV. Hence for $`g_a>10^{10}`$ the decay is essential, i.e. $`\tau _{\nu _\tau }<H^1`$, while still $`T>m_{\nu _\tau }`$. The interval of life-times of $`\nu _\tau `$, that we consider below, $`\tau _{\nu _\tau }=0.001400`$ sec, corresponds to $`g_a=\left(410^96.510^{12}\right)\sqrt{m/\text{MeV}}`$. Thus there is a large range of parameters (coupling constants and masses) for which the decay is faster than the expansion while annihilation is effectively frozen.

## 3 Numerical solution

### 3.1 Kinetic equations

The basic equations governing the evolution of the distribution functions $`f_a`$ ($`a=\nu _e`$, $`\nu _\mu `$, $`\nu _\tau `$, and $`\varphi `$) are discussed in some detail in our previous works . The tables of matrix elements squared of the relevant weak interaction processes, $`A^2`$, are presented there. We have some disagreement of our expressions for $`A^2`$ with those of refs. . In addition to the processes considered for the case of a stable $`\nu _\tau `$ we should also take decay and inverse decay of $`\nu _\tau `$ into account and include an extra distribution function, $`f_\varphi `$. The distributions of photons and electrons/positrons are assumed to have the equilibrium form:

$`f_{e,\gamma }=\left[1\pm \mathrm{exp}\left(E/T_\gamma \right)\right]^1,`$ (3)   
because the electromagnetic interactions are very strong and keep the electromagnetic plasma in thermal equilibrium. The behaviour of the unknown function of one variable $`T(t)`$ is found from the law of covariant energy conservation:

$`\dot{\rho }=3H(\rho +P),`$ (4)   
where the Hubble parameter, $`H=\dot{a}/a`$, is expressed through the total energy density, $`\rho `$, as $`H^2=(8\pi \rho /3m_{Pl}^2)`$; here $`a(t)`$ is the universe expansion factor (scale factor), P is the total pressure and $`m_{Pl}=1.2210^{19}`$ GeV is the Planck mass.

The set of integro-differential kinetic equations has the form:

$$(_tHp_j_{p_j})f_j(p_j,t)=I_j^{scat}+I_j^{decay},$$
(5)   
where the collision integral for two-body reactions $`1+23+4`$ is given by the expression:

$`I_1^{scat}={\displaystyle \frac{1}{2E_1}}{\displaystyle \frac{d^3p_2}{2E_2(2\pi )^3}\frac{d^3p_3}{2E_3(2\pi )^3}\frac{d^3p_4}{2E_4(2\pi )^3}S|A|_{1234}^2}`$   
$`(2\pi )^4\delta ^{(4)}(p_1+p_2p_3p_4)F_{scat}(f_1,f_2,f_3,f_4),`$ (6)   
$`F_{scat}=f_3f_4(1f_1)(1f_2)f_1f_2(1f_3)(1f_4)`$, $`|A|^2`$ is the weak interaction amplitude squared summed over spins of all particles, and $`S`$ is the symmetrization factor which includes $`1/2`$ from the averaging over the spin of the first particle (if necessary), $`1/2!`$ for each pair of identical particles in initial and final states and the factor 2 if there are 2 identical particles in the initial state; the summation is done over all possible sets of leptons 2, 3, and 4.

The ”decay” part of the collision integrals enters the r.h.s. of the equations for $`f_{\nu _\tau }`$, $`f_{\nu _a}`$ ($`a=e`$ or $`\mu `$), and $`f_\varphi `$ and has respectively the form:

$`I_{\nu _\tau }^{decay}={\displaystyle \frac{m}{E_{\nu _\tau }p_{\nu _\tau }\tau _{\nu _\tau }}}{\displaystyle _{\left(E_{\nu _\tau }+p_{\nu _\tau }\right)/2}^{\left(E_{\nu _\tau }p_{\nu _\tau }\right)/2}}𝑑E_\varphi F_{dec}(E_{\nu _\tau },E_\varphi ,E_{\nu _\tau }E_\varphi ),`$ (7)   
$`I_{\nu _a}^{decay}={\displaystyle \frac{m}{E_{\nu _a}p_{\nu _a}\tau _{\nu _\tau }}}{\displaystyle _{\left(m^2/4p_{\nu _a}\right)p_{\nu _a}}^{\mathrm{}}}{\displaystyle \frac{dp_{\nu _\tau }p_{\nu _\tau }}{E_{\nu _\tau }}}F_{dec}(E_{\nu _\tau },E_{\nu _\tau }E_{\nu _a},E_{\nu _a}),`$ (8)   
$`I_\varphi ^{decay}={\displaystyle \frac{2m}{E_\varphi p_\varphi \tau _{\nu _\tau }}}{\displaystyle _{\left(m^2/4p_\varphi \right)p_\varphi }^{\mathrm{}}}{\displaystyle \frac{dp_{\nu _\tau }p_{\nu _\tau }}{E_{\nu _\tau }}}F_{dec}(E_{\nu _\tau },E_\varphi ,E_{\nu _\tau }E_\varphi ),`$ (9)   
where $`m`$ is the mass of $`\nu _\tau `$ (we omitted the index $`\nu _\tau `$) and:

$`F_{dec}(E_{\nu _\tau },E_\varphi ,E_{\nu _a})=f_{\nu _\tau }(E_{\nu _\tau })\left[1+f_\varphi (E_\varphi )\right]\left[1f_{\nu _a}(E_{\nu _a})\right]`$   
$`f_\varphi (E_\varphi )f_{\nu _a}(E_{\nu _a})\left[1f_{\nu _\tau }(E_{\nu _\tau })\right].`$ (10)   
The contribution of the decay term, $`I^{decay}`$, into the collision integral of eq. (5) is considerably simpler for numerical calculations than the contribution of scattering, $`I^{scat}`$, because the former is only one-dimensional while the scattering terms only can be reduced to two-dimensional integrals.

Four integro-differential kinetic equations for the distribution functions $`f_j(p,t)`$ $`(j=\nu _\tau ,\nu _\mu ,\nu _e,\varphi )`$ and eq. (4) are solved numerically by the method described in refs.  (see also the earlier papers ). Technical details of the calculations are presented in the following section.

It is convenient instead of time and momenta to use the following dimensionless variables:

$$x=ma(t),y_j=p_ja(t),$$
(11)   
where $`m`$ is an arbitrary parameter with dimension of mass, which we took as $`m=1`$ MeV, and the scale factor $`a(t)`$ is normalized so that $`a(t)=1/T_\nu =1/T_\gamma `$ at high temperatures or at early times. In terms of these variables the kinetic equations (5) can be rewritten as:

$$Hx_xf_j(x,y_j)=I_j^{scat}+I_j^{dec}.$$
(12)   
### 3.2 Initial conditions

For numerical integration of the kinetic equations we modified our program for stable massive tau neutrinos (see the details of the stable case calculation in our paper ) by including the decay terms eq. (7-9) in equations (12).

All the calculations of the two-body-reaction part of the collision integral remain the same as in our previous work. So we discuss below only the changes to the program related to the decay terms. We solve the system of kinetic equations (12) in the ’time’ interval $`x_{in}xx_{fin}`$ under several assumptions concerning decay. First, we assume that the decay was negligible before the initial time $`x_{in}`$ and take the initial majoron distribution function equal to zero, $`f_\varphi (x_{in},y)=0`$. We can always do this for sufficiently small $`x`$, because the collision rates in eq. (12) are very high due to the factor $`1/x^2`$ while the decay terms are proportional to $`x^2/\tau `$ and even for very small life-times, $`\tau 1`$ sec they are suppressed at small $`x`$. To determine the proper initial moment $`x_{in}`$ we approximately solved the kinetic equation for $`f_\varphi (x,y)`$ neglecting the contribution from $`\varphi `$ in the collision integral and assuming equilibrium distributions for all other participating particles. Then we can integrate this equation over $`\varphi `$ momentum and obtain the following equation for the energy density of majorons at early times:

$$_x(x^4\rho _\varphi )=2C_\varphi x^2\frac{m_{\nu _\tau }}{\sqrt{x^4\rho _{tot}}\tau }\frac{1}{2\pi ^2}_0^{\mathrm{}}\frac{ydy}{e^y1}\mathrm{log}\frac{1+\mathrm{exp}(\frac{m_{\nu _\tau }^2x^2}{4y})}{\mathrm{exp}(y)+\mathrm{exp}(\frac{m_{\nu _\tau }^2x^2}{4y})},$$
(13)   
where $`C_\varphi =1.2216.582\sqrt{3/(8\pi )}=2.777`$ is a constant and $`x^4\rho _{tot}=10.75\pi ^2/30`$ is the total energy density of photons, $`e^\pm `$, and the three neutrino species in comoving frame, and the lifetime, $`\tau `$, is given in seconds.

For a small value of the product $`\kappa =(m_{\nu _\tau }^2x^2/4)1`$ we can solve this equation and obtain:

$`\rho _\varphi `$ $`=`$ $`2C_\varphi {\displaystyle \frac{x^3}{3x^4}}{\displaystyle \frac{m_{\nu _\tau }}{\sqrt{\rho _{tot}}\tau }}({\displaystyle \frac{\zeta (3)}{2\pi ^2}}{\displaystyle \frac{\mathrm{log}(2)}{24}})(10.08\kappa +0.006\kappa ^2+O(\kappa ^3))`$ (14)   
$``$ $`0.31{\displaystyle \frac{m_{\nu _\tau }}{x\tau }}(10.08\kappa +0.006\kappa ^2+O(\kappa ^3)).`$   
From equation (14) we can find the energy density of scalars for a given initial time $`x_{in}`$ and compare it with the energy density of equilibrium massless neutrinos $`\rho _{eq}=0.57573/x^4`$. We will assume that scalars are absent if $`\rho _\varphi /\rho _{eq}<0.01`$. For $`\nu _\tau `$ lifetime in the range $`1<\tau /\text{sec}<\mathrm{}`$ we take the initial time $`x_{in}=0.1`$ for all values of the mass. For $`0.1<\tau /\text{sec}<1`$ we take $`x_{in}=0.1`$ for $`\nu _\tau `$ mass in the range $`0.1<m_{\nu _\tau }/\text{MeV}<2`$ and $`x_{in}=0.05`$ for $`2<m_{\nu _\tau }/\text{MeV}<20`$. For $`\tau /\text{sec}<0.1`$ we always take $`x_{in}=0.05`$.

### 3.3 Choice of the final time

The contributions to the collision integral in the r.h.s. of eq. (12) from two-body reactions are suppressed at large $`x`$ by at least the factor $`1/x^2`$. We find that at $`x50`$ all the quantities, which do not feel the decay, reach their asymptotic values so, to be on the safe side, we choose the final time for the reactions $`x_{scat}=100`$. Beyond this time all collision terms are taken to be zero.

The decay continues until all the tau-neutrinos have disappeared. In our runs we assume that the decay terminates at $`x_{decay}`$ at which $`\rho _{\nu _\tau }/\rho _{eq}<10^4`$, where $`\rho _{eq}=0.57573/x^4`$ is the energy density of one massless equilibrium neutrino. As the final time $`x_{fin}`$ we take the maximum of the two, $`x_{decay}`$ and $`x_{scat}`$.

The nucleosynthesis code of ref.  requires the final time $`x2000`$ for the total energy density, photon temperature and $`np`$ rates. We applied a separate program, which calculates these quantities, using the energy conservation law and our final values for the distribution functions $`f_{\nu _i}`$ and the temperature $`T_\gamma `$ at time $`x_{fin}`$.

### 3.4 Momentum grid

The particles produced by $`\nu _\tau `$ decays have momenta around $`y=xm_{\nu _\tau }/2`$. The decay becomes operative and starts to diminish the $`\nu _\tau `$ number density at the time $`t\tau _{\nu _\tau }`$ or at $`x_D\sqrt{\tau _{\nu _\tau }}`$, if the inverse decay is not important, i.e. if the mass and life-time are sufficiently large. Thus for the decay products we should take into account momenta up to $`y_D=C\sqrt{\tau _{\nu _\tau }}m_{\nu _\tau }`$, where the coefficient $`C`$ depends on the time when the decay is completed.

On the other hand the distribution functions of the tau neutrino and of the massless neutrino, which does not participate in the decay, are suppressed at high momenta. Already at $`y=20`$ we have $`\rho (y>20)/\rho 10^5`$ for these particles, and we therefore choose the momentum cut-off equal to $`y_L=20`$. The cases $`y_D>y_L`$ and $`y_D<y_L`$ are considered separately.

When $`y_D<y_L`$ (it corresponds to $`\tau _{\nu _\tau }<0.1`$ sec for $`m_{\nu _\tau }>1`$ MeV and $`\tau _{\nu _\tau }<1`$ sec for $`m_{\nu _\tau }<1`$ MeV) we take a logarithmic grid for the Majorons in the region $`0.01<y<20`$ and a linear grid for the neutrinos in the region $`0<y<20`$.

Since the collision integrals for two-body reactions are two-dimensional, their calculation takes most of the CPU time, and correspondingly we confined ourselves to a 100 point grid for the reaction terms. On the other hand, the decay terms are one-dimensional and we took a 1600 point grid in their calculation.

When $`y_D>y_L`$ (it corresponds to $`\tau _{\nu _\tau }>0.1`$ sec for $`m_{\nu _\tau }>1`$ MeV and $`\tau _{\nu _\tau }>1`$ sec for $`m_{\nu _\tau }<1`$ MeV) we took a linear grid for the decay products in the interval $`0<y<y_D`$ while for the neutrino, which does not participate in the decay, the linear grid was taken in the range $`0<y<20`$. In the interval $`0<y<20`$ we took 100 points in the grid, and expanded with an equally spaced grid with the same $`dy`$ for the decay products up to $`y=y_D`$. We calculated the reaction terms in the collision integrals only in the interval $`0<y<20`$, because the decay into modes with $`y>20`$ becomes important when the contribution of reactions into the collision integral is already suppressed by the factor $`1/x^4`$ with $`x1`$ for any mass.

At $`y=0`$ the analytical expressions for the collision integral for $`\nu _\tau `$ differ from those at $`y0`$ due to $`0/0`$-uncertainty arising from the factor $`1/E_1p_1`$ in front of the integrals and similar vanishing factors in the D-functions (see their definition in our previous papers ). Because of that, we consider the point $`y=0`$ separately. This permits to compare the collision integrals at $`y=0`$ with those in nearby points, $`y1`$, in order to check that numerical errors are small in the region of small momentum, $`y<1`$. This is important for a massive tau-neutrino, because its distribution function rapidly changes in the course of evolution.

An important check of our program is the comparison of the run with the large life-time, $`\tau _{\nu _\tau }=400`$ sec, with that with an infinite life-time, $`\tau _{\nu _\tau }=\mathrm{}`$, (when the decay terms in the collision integral are switched off) for masses in the interval $`0.1<m_{\nu _\tau }/\text{MeV}<20`$. The results of both these runs are very close for all masses.

### 3.5 Nucleosynthesis

In order to extract the implications of a massive decaying $`\nu _\tau `$ on light element abundances, we have modified the standard nucleosynthesis code (ref. ) in the following way. First we import the final photon temperature, which is different from $`1.401/x`$, in order to specify the final value of the baryon to photon ratio, $`\eta _{10}=10^{10}n_B/n_\gamma `$. We calculate various quantities at the correct photon temperature and import the values to the code. These imported quantities are the total energy density $`\rho _{tot}`$, the 6 weak interaction rates for $`(np)`$-reactions, and finally $`d(\mathrm{ln}a^3)/dT_\gamma `$ (with the account of neutrinos) which governs the evolution of the photon temperature. The excessive number of equivalent neutrinos, $`\mathrm{\Delta }N`$, which may be both positive and negative, is interpolated from a look-up table as a function of <sup>4</sup>He and $`\eta _{10}`$.

## 4 Results

The impact on nucleosynthesis of a heavy decaying tau neutrino strongly depends upon the decay channel. In the case of the decay into $`\nu _\mu +\varphi `$ the most important effect is the overall change in the total energy density and the corresponding change of the universe cooling rate. Nonequilibrium corrections to the spectra of $`\nu _e`$ are relatively weak for a small life-time, such that practically all $`\nu _\tau `$ have already decayed at the moment of neutron-proton freezing, $`T0.5`$ MeV. For a larger life-time some nonequilibrium $`\nu _e`$ would come from annihilation $`\nu _\tau +\nu _\tau \overline{\nu }_e+\nu _e`$ and, as is well known, would directly change the frozen $`n/p`$-ratio. The distortion of the $`\nu _e`$ spectrum is much stronger in the case of the decay $`\nu _\tau \nu _e+\varphi `$. Correspondingly the results for the light element abundances are considerably different for these two different decay channels.

In figs. 1 the total energy density of all three neutrino species and of the scalar $`\varphi `$ is presented as a function of the tau neutrino mass and life-time for asymptotically large values of ”time” $`x`$. It is given in terms of the effective number of massless equilibrium neutrinos, $`N_{eff}=(\rho _{\nu _\tau }+\rho _{\nu _\tau }+\rho _{\nu _\tau }+\rho _\varphi )/\rho _{eq}`$, where $`\rho _{eq}=(7\pi ^2/120)m_0^4/x^4`$ is the energy density of one massless neutrino. To avoid confusion we should stress, that this effective number of neutrinos, $`N_{eff}`$, has nothing to do with the excessive number of equivalent neutrinos, $`\mathrm{\Delta }N`$, which is introduced to describe the change in the primordial <sup>4</sup>He abundance (see below). The behavior of the ”iso-neutrino” curves is quite transparent: for a large life-time the relative contribution from nonrelativistic $`\nu _\tau `$ into $`\rho `$ increases as the first power of the scale factor $`ax`$. For larger masses the $`\nu _\tau `$ becomes nonrelativistic earlier, and its contribution into the energy density is bigger. On the other hand, very heavy tau-neutrinos annihilate more efficiently and thus their frozen number should be smaller than the number density of less massive ones. This explains the extrema in fig. 1a. For small masses and lifetimes a plateau is reached, where a fraction of percent uncertainty makes the iso-curves wiggle. This uncertainty arises from merging of different programs for different mass-lifetime regions. Somewhat less trivial is the fact, that for larger values of mass and relatively small life-times the energy density can be smaller than in the standard model. This part of the graph is enlarged and presented in fig. 1b. A decline in the energy density in this range is related to the early decay of $`\nu _\tau `$, when its number density is not yet frozen with respect to annihilation and is Boltzmann suppressed.

In figs. 2 and 3 the time evolution of the energy (a) and number (b) densities of different species are presented. The behavior of $`\nu _\mu `$ and $`\varphi `$ reflect the behavior of $`\nu _\tau `$. For a small mass and life-time (fig. 2) $`\nu _\tau `$ remains relativistic and its number and energy densities follow the equilibrium massless ones. For $`x>1`$, when the decay starts to operate, some decrease in $`n_{\nu _\tau }`$ is observed due to filling of the previously empty scalar states. However, the energy density, $`\rho _{\nu _\tau }`$, does not decrease starting from $`x5`$ due to the above mentioned effect of domination of nonrelativistic particles. Later, at $`x>20`$, the Boltzmann suppression becomes essential and both $`n_{\nu _\tau }`$ and $`\rho _{\nu _\tau }`$ quickly drop down. For large mass and life-time (fig. 3) the features in the picture are more pronounced. The decrease in $`n_{\nu _\tau }`$ and $`\rho _{\nu _\tau }`$ for $`0.1<x<0.5`$ is related to the Boltzmann suppression of almost equilibrium $`\nu _\tau `$. After that, the $`\nu _\tau `$-annihilation becomes frozen and the energy density goes up. Above $`x=3`$ the decay comes into play and $`\nu _\tau `$ disappear from the primeval plasma. To test the consistency of our calculations we checked the conservation of particle number densities and found it valid with an accuracy better than 1%.

Our main results for the decay channel $`\nu _\tau \nu _\mu +\varphi `$ are presented in figs. 4. We calculated the production of <sup>4</sup>He for different masses and life-times, using the nucleosynthesis code modified as described in section 3.5, and expressed the variation of the mass fraction of <sup>4</sup>He in terms of the extra number of neutrino species which gives rise to the same amount of <sup>4</sup>He. The results are presented in two graphical forms: as a three dimensional graph $`\mathrm{\Delta }N=F(m_{\nu _\tau },\tau _{\nu _\tau })`$ (fig. 4a) and as a series of curves $`\mathrm{\Delta }N=F(m_{\nu _\tau })`$ for different values of $`\tau _{\nu _\tau }`$ (fig. 4b). For large masses and low life-times $`\mathrm{\Delta }N`$ is negative. This is related to the decrease of the energy density discussed above. Thus if $`m_{\nu _\tau }=10`$ MeV and $`\tau _{\nu _\tau }=0.1`$ sec, the effective number of neutrino species at nucleosynthesis is only 2.5.

Comparison with the results of other groups shows a rather strong deviation. We ascribe this to the simplifying approximations made in the earlier papers, which might give rise to a significant difference with the exact calculations, and to a better accuracy of our numerical calculations, which is typically at the fraction of per cent level. For example in the case of $`m_{\nu _\tau }=14`$ MeV and $`\tau _{\nu _\tau }=0.1`$ sec we obtain for the energy density $`\rho /\rho _{\nu _0}^{eq}=2.9`$, while the group  obtained 2.5. In the limit of small life-times and masses our results are close to 3.57 (this is the energy density of three light neutrinos and one scalar), while the results of  are close to 3.9. The effective number of neutrino species found from <sup>4</sup>He in our case is $`3+\mathrm{\Delta }N=2.9`$ for $`m_{\nu _\tau }=10`$ MeV and $`\tau _{\nu _\tau }=1`$ sec, while that found in ref.  is 3.1. The difference is also large for $`m_{\nu _\tau }=10`$ MeV and $`\tau _{\nu _\tau }=0.01`$ sec: we find $`3+\mathrm{\Delta }N=2.66`$ and the authors of  obtained 2.86.

Consideration of the decay $`\nu _\tau \nu _e+\varphi `$ proceeds essentially along the same lines. The energy densities of different particle species in this case are almost the same as in the case of the decay $`\nu _\tau \nu _\mu +\varphi `$, so we do not present them. However, the influence on light element abundances is now very different from the previous case, because the nonequilibrium electronic neutrinos produced by the decay strongly shift the $`n/p`$ ratio. In particular an excess (with respect to the equilibrium distribution) of low energy electronic neutrinos would create a smaller $`n/p`$-ratio and thus a correspondingly smaller mass fraction of <sup>4</sup>He.

The light element abundances for different values of $`m_{\nu _\tau }`$, $`\tau _{\nu _\tau }`$ and the baryon-to-photon ratio $`\eta _{10}=10^{10}n_B/n_\gamma `$ are summarized in figs. 5-11. Massive unstable $`\nu _\tau `$ may act both ways, diminishing and increasing the mass fraction of <sup>4</sup>He. In terms of the effective number of neutrino species it may correspond to a large and positive $`\mathrm{\Delta }N`$, up to $`\mathrm{\Delta }N=2.5`$, or even to the negative $`\mathrm{\Delta }N=2`$ (see fig. 7). As seen on the $`\mathrm{\Delta }N`$ iso-curves (figs. 8), there is a minimum in $`\mathrm{\Delta }N`$ for lifetime $`0.05`$ sec and mass $`4`$ MeV.

Also for this channel we disagree with previous papers. E.g. for $`m=0.6`$ MeV and $`\tau =100`$ sec we find $`Y({}_{}{}^{4}\text{He})0.244`$, whereas ref.  obtains $`Y({}_{}{}^{4}\text{He})0.20`$.

It was shown in ref.  that late decaying $`\nu _\tau `$ with $`\tau _{\nu _\tau }=10^310^4`$ sec and $`m_{\nu _\tau }>3.6`$ MeV would distort the deuterium abundance strongly if the decay proceed into electronic neutrinos. These $`\nu _e`$ would create excessive neutrons through the reaction $`\nu _e+pn+e^+`$, which would form extra <sup>2</sup>H. This is seen explicitly on fig. 9b, where <sup>2</sup>H clearly increases as a function of lifetime. The deuterium production goes quadratically with the baryon density, and we see, that this effect is much less pronounced for low $`\eta _{10}`$.

In figs. 9-11a the mass fraction of <sup>4</sup>He is presented as a function of $`m_{\nu _\tau }`$ for several values of $`\tau _{\nu _\tau }`$ for different baryon-to-photon ratio, $`\eta _{10}=10^{10}n_B/n_\gamma =3,5,7`$. In figs. 9-11b,c the similar graphs for the primordial <sup>2</sup>H and <sup>7</sup>Li are presented.

More graphs showing various elements (<sup>2</sup>H, <sup>4</sup>He and <sup>7</sup>Li) and as functions of mass and lifetime, for both channels $`\nu _\tau \nu _\mu +\varphi `$ and $`\nu _\tau \nu _e+\varphi `$ and for $`\eta _{10}=1,3,5,7,9`$ can be found on the web-page: http://tac.dk/~sthansen/decay/ together with plots of the n-p reaction rates.

Finally we changed the initial condition, allowing the majorons to be in equilibrium at $`x_{in}`$, $`f_\varphi (x_{in})=1/(\mathrm{exp}(y)1)`$. This situation corresponds to the majorons being produced by some other mechanism prior to the $`\nu _\tau `$ decay as discussed in section 2. In this case the inverse decay is more efficient than for $`f_\varphi (x_{in})=0`$ and $`\rho _{\nu _\tau }`$ decreases slower as is seen from fig. 12. The change in $`\mathrm{\Delta }N`$, as compared to the case when $`f_\varphi (x_{in})=0`$, varies between 0.4 and 1.0, $`\mathrm{\Delta }N_{f_\varphi =eq}=\mathrm{\Delta }N_{f_\varphi =0}+(0.41.0)`$, depending on mass and lifetime. In particular for long lifetime this difference goes to 0.57 as expected for all masses.

## 5 The role of possible neutrino oscillations

The considerations of this paper neglect a possible influence of neutrino oscillations on primordial nucleosynthesis. If the Super Kamiokande result  is true, then $`\nu _\mu `$ is strongly mixed with another neutrino, either $`\nu _\tau `$ or a sterile one. If $`\nu _\mu `$ is mixed with $`\nu _\tau `$, then the mass of the latter is small, $`m_{\nu _\tau }<160`$ keV and our results obtained for larger masses would be irrelevant. In the other possible case, when oscillations are between $`\nu _\mu `$ and a sterile one, $`\nu _s`$, the $`\nu _\tau `$ mass is not restricted other than by the direct laboratory limit  and our calculations make sense and permit to improve the latter significantly. Still the oscillations into a sterile neutrino state would have an important impact on nucleosynthesis and would change our results obtained without oscillations. We postpone a detailed numerical calculations of the influence of neutrino oscillations on primordial nucleosynthesis for the future study and give only some simple estimates here (a brief discussion of the impact of the Super Kamiokande results on nucleosynthesis is also given in ref. ).

In accordance with the measurement  the mixing angle is surprisingly large, $`\mathrm{sin}^22\theta =0.81`$ and vacuum oscillations are not suppressed at all. Interactions with the thermal bath generally result in a suppression of oscillations, except for a possible MSW-resonance. This has been studied in a large number of works and a list of references can be found in the recent papers . It is sufficient for our purposes to use the simple results of the early papers , that the oscillations $`\nu _\mu \nu _s`$ would be fully developed if $`\mathrm{sin}^42\theta |\delta m^2|310^3\text{eV}^2`$  and $`\mathrm{sin}^22\theta \delta m^2>310^6\text{eV}^2`$ . Thus for $`|\delta m^2|10^3`$ eV<sup>2</sup> the oscillations $`\nu _\mu \nu _s`$, in the case that $`\nu _\mu `$ were produced by the usual weak interaction processes in the primeval plasma, would not change our results considerably. The oscillations of muonic neutrinos created in the $`\nu _\tau `$ decays would not influence nucleosynthesis if the $`\nu _\tau `$ life-time is large enough. In this limit the oscillations do not change the total energy density of the plasma in comparison with the no-oscillation case. The decays with a short life-time and the subsequent oscillations would somewhat amplify the production of extra neutrino species in the primeval plasma and might be noticeable in nucleosynthesis.

The oscillations might be very strongly suppressed if the cosmic lepton asymmetry (in the case considered it is muon asymmetry) is large in comparison with the usually accepted value $`10^910^{10}`$. The calculations made in this work should be correspondingly modified to take into account the change in the distribution function, $`f_{\nu _\mu }`$, due to the non-vanishing muonic chemical potential. However, even an increase of muon asymmetry by 6-7 orders of magnitude would have a negligible effect on the kinetics of nucleosynthesis. This is not so for the oscillations of $`\nu _e\nu _s`$, when a large electric asymmetry might be generated by the oscillations and could shift the equilibrium value of $`n/p`$-ratio (for a recent work and the list of references see ).

If one sterile neutrino species exists, then it is rather natural to expect three sterile ones. This possibility could be realized for example if the mass matrix of neutrinos contains both Majorana and Dirac mass terms . In this case, if $`\nu _\tau `$ indeed is heavy and unstable, the anomalies in neutrino physics should/could be explained by the oscillations in the sector $`\nu _e\nu _\mu \nu _{s1}\nu _{s2}`$. Due to these oscillations there might be effectively up to two extra neutrino species at nucleosynthesis, but a massive and unstable $`\nu _\tau `$ which could effectively create a negative $`\mathrm{\Delta }N`$ would permit to avoid a contradiction with the data.

For the case of neutrino interactions with majorons the contribution into the neutrino refraction index from majorons in the medium would strongly suppress oscillations , if the $`\nu _\tau \nu _\mu \varphi `$-coupling constant is larger than $`10^610^7`$. It corresponds to the very short life-time $`\tau _{\nu _\tau }10^6`$ sec $`(10^7/g)^2(\text{MeV}/m_{\nu _\tau })`$.

Thus to summarize, our calculations should be considerably modified if neutrino oscillations, as indicated by Super Kamiokande data, indeed exist. However, first, it seems premature to make a strong conclusion based on these data and, second, there exist regions in the parameter space for which our results with the neglect of oscillations remain relevant.

Acknowledgments

It is a pleasure to thank J.W.F. Valle for useful discussions. We also thank K. Kainulainen for helpful comments. The work of AD and SH was supported in part by the Danish National Science Research Council through grant 11-9640-1 and in part by Danmarks Grundforskningsfond through its support of the Theoretical Astrophysical Center. The work of DS was supported in part by the Russian Foundation for Fundamental Research through grants 97-02-17064A and 98-02-17493A. SP was supported by Conselleria d’Educació i Ciència of Generalitat Valenciana, by DGICYT grant PB95-1077 and by the EEC under the TMR contract ERBFMRX-CT96-0090. SP and DS thank TAC for hospitality when this work was done.

Figure Captions:

Fig. 1 Asymptotic values for the effective number of neutrinos $`N_{eff}=(\rho _{\nu _\tau }+\rho _{\nu _\tau }+\rho _{\nu _\tau }+\rho _\varphi )/\rho _{eq}`$, where $`\rho _{eq}=(7\pi ^2/120)m_0^4/x^4`$ is the energy density for one massless neutrino. In fig. 1a we present the region $`0.1<\tau /\text{sec}<400`$ and $`0.1<m_{\nu _\tau }/\text{MeV}<20`$. In fig. 1b we show the large mass region with decay times in the region $`0.01<\tau /\text{sec}<1`$. In this region the effective number of neutrinos can be less than 3.

Fig. 2 For $`m_{\nu _\tau }=0.1`$ MeV and $`\tau _{\nu _\tau }=1`$ sec we present the evolution in time $`x`$ of the energy densities (fig. 2a) and number of particles densities (fig. 2b) in units of massless equilibrium neutrinos ($`\rho _{eq}=(7\pi ^2/120)m_0^4/x^4`$ and $`n_{eq}=0.183m_0^3/x^3`$). In the time region $`1<x<6`$ the massive $`\nu _\tau `$ decays to $`\nu _\mu `$ and the majoron $`\varphi `$, but because of the small mass, the inverse decay-process delays the complete vanishing of $`\nu _\tau `$ until time $`x=25`$.

Fig. 3 For $`m_{\nu _\tau }=10`$ MeV and $`\tau _{\nu _\tau }=10`$ sec we present the evolution in time $`x`$ of energy densities (fig. 3a) and number of particles densities (fig. 3b) in units of massless equilibrium neutrinos ($`\rho _{eq}=(7\pi ^2/120)m_0^4/x^4`$ and $`n_{eq}=0.183m_0^3/x^3`$). In the time region $`0.1<x<3`$ the decay does not play any role, in the region $`x<1`$ the annihilation of $`\nu _\tau `$ into other particles is important. Massive $`\nu _\tau `$ decays to $`\nu _\mu `$ and majoron $`\varphi `$ around $`x=3`$.

Fig. 4 Relative number of equivalent massless neutrino species, $`\mathrm{\Delta }N=N_{eqv}3`$, as a function of $`\nu _\tau `$ mass and lifetime $`\tau `$, found from <sup>4</sup>He. In fig. 4a we present a 2-dimensional surface for $`\mathrm{\Delta }N(m_{\nu _\tau },\tau )`$. In fig. 4b $`\mathrm{\Delta }N`$ as a function of $`\nu _\tau `$ mass for different lifetimes including the asymptotic value for stable massive $`\nu _\tau `$ (decay time $`\tau =\mathrm{}`$) is presented.

Fig. 5 <sup>4</sup>He for different lifetimes as a function of the $`\nu _\tau `$ mass. $`\eta _{10}=3`$ (fig. 5a). <sup>2</sup>H for different lifetimes (fig. 5b). <sup>7</sup>Li for different lifetimes (fig. 5c).

Fig. 6 Fig. 6a shows $`\mathrm{\Delta }N`$ contours found from <sup>4</sup>He as a function of $`\nu _\tau `$ mass and lifetime. $`\eta _{10}=3`$. The contours shown correspond to $`\mathrm{\Delta }N=`$ -0.6, -0.4, -0.2, -0.1, 0.0, 0.1, 0.2, 0.4, 0.6. Fig. 6b shows the same as fig. 6a with small lifetimes $`0.001\tau /\text{sec}0.1`$.

$`\nu _\tau \nu _e+\varphi `$

Fig. 7a $`\mathrm{\Delta }N`$ surface found from <sup>4</sup>He as a function of $`\nu _\tau `$ mass and lifetime. $`\eta _{10}=3`$.

Fig. 7b $`\mathrm{\Delta }N`$ for different lifetimes as a function of $`\nu _\tau `$ mass. $`\eta _{10}=3`$.

Fig. 8a $`\mathrm{\Delta }N`$ contours for $`\eta _{10}=3`$ as a function of $`\nu _\tau `$ mass and lifetime. Contours follow the increase with the step of $`0.2`$.

Fig. 8b Same as fig. 8a with small lifetimes $`0.001\tau /\text{sec}0.1`$.

Fig. 8c <sup>4</sup>He contours for $`\eta _{10}=3`$ as a function of $`\nu _\tau `$ mass and lifetime. Contours follow the increase with the step of $`0.05`$.

Fig. 9 <sup>4</sup>He (9a), <sup>2</sup>H (9b) and <sup>7</sup>Li (9c) as functions of the $`\nu _\tau `$ mass for different lifetimes. $`\eta _{10}=3`$.

Fig. 10 <sup>4</sup>He (10a), <sup>2</sup>H (10b) and <sup>7</sup>Li (10c) as functions of the $`\nu _\tau `$ mass for different lifetimes. $`\eta _{10}=5`$.

Fig. 11 <sup>4</sup>He (11a), <sup>2</sup>H (11b) and <sup>7</sup>Li (11c) as functions of the $`\nu _\tau `$ mass for different lifetimes. $`\eta _{10}=7`$.

Fig. 12  Energy density of $`\nu _\tau `$ and $`\varphi `$ as functions of time $`x`$ for the cases with $`f_\varphi (x_{in})=1/(\mathrm{exp}(y)1)`$ and $`f_\varphi (x_{in})=0`$.

Figure 1a.

Figure 1b.

Figure 2a.

Figure 2b.

Figure 3a.

Figure 3b.

Figure 4a.

Figure 4b.

Figure 5a.

Figure 5b.

Figure 5c.

Figure 6a.

Figure 6b.

Figure 7a.

Figure 7b.

Figure 8a.

Figure 8b.

Figure 8c.

Figure 9a.

Figure 9b.

Figure 9c.

Figure 10a.

Figure 10b.

Figure 10c.

Figure 11a.

Figure 11b.

Figure 11c.

Figure 12.