# 1 The evolution of neutrino asymmetry at the resonance temperature for mixing parameters 𝛿⁢𝑚²=-10⁻²⁢"eV"² and sin²{2⁢𝜃}=10^{-7.5}. Solid line corresponds to initial value 𝜂=10⁻¹⁰ and dashed line to 𝜂=2×10⁻¹⁰.

The recent observation of strong zenith angle dependence of the atmospheric neutrino deficit by the Super-Kamiokande neutrino experiment has provided strong evidence for neutrino oscillations $`\nu _\mu \nu _X`$, where $`\nu _X`$ is either $`\nu _\tau `$ or a new, sterile neutrino $`\nu _s`$ . While the $`\nu _\mu \nu _s`$ solution presently is less favored by SK data , reconciling the existing data from all the neutrino experiments, including LSND, is not possible unless there exists at least one sterile neutrino mixing with the active neutrinos. Such mixing would have interesting consequences for primordial nucleosynthesis  and CMB radiation . For example, sterile neutrinos could be brought into equilibrium prior to nucleosynthesis, increasing the energy density of the universe and thereby neutron-to-proton freeze out temperature, leading to more helium-4 being produced. This scenario has been numerically studied and strong limits to neutrino mixing parameters have been obtained .

Under certain conditions active-sterile neutrino oscillations may also lead to exponential growth of neutrino asymmetry, as was first discovered by Barbieri and Dolgov . Later Foot and Volkas observed that by this mechanism very large asymmetries could be generated, which would have a significant effect on the primordial nucleosynthesis  by directly modifying directly the $`np`$ reactions. Moreover, they showed that an asymmetry generated by $`\nu _\tau \nu _s`$ mixing could suppress sterile neutrino production in $`\nu _\mu \nu _s`$-sector, loosening the bounds of , according to which the SK atmospheric deficit could not be explained by $`\nu _\mu \nu _s`$ mixing.

Later it was found that this asymmetry generation is chaotic in the sense that determining the sign of the final asymmetry sign($`L`$) does not simply follow from the initial conditions . This phenomenon was studied in , and it was shown that the indeterminacy is associated with a region of mixing parameters, where asymmetry is rapidly oscillating right after the resonance. As a consequence the the amount of Helium-4 produced cannot be precisely determined in such a scenario .

In a recent paper  it was claimed, however, that sign($`L`$) is completely determined by the initial asymmetry, and moreover, that there is only a slight growth of asymmetry after the resonance. In this article we clarify the origin of indeterminacy in sign($`L`$) and show that it is a real physical phenomenon, not disturbed by numerical inaccuracy. Instead, we will argue that the disagreement arises due to overly simplifying approximations used in . Finally, we will point out to a likely cause leading to observed suppression of the final magnitude of the asymmetry in .

In the early universe the coherent evolution of the neutrino states is interrupted by frequent decohering collisions. Therefore the evolution of the system needs to be studied using the density matrix formalism. We parameterize the reduced density matrices of the neutrino and anti-neutrino ensembles as

$$\rho _\nu \frac{1}{2}P_0(1+𝐏),\rho _{\overline{\nu }}\frac{1}{2}\overline{P}_0(1+\overline{𝐏}),$$
(1)   
where each matrix is assumed to be diagonal in momentum space, while each momentum state has $`2\times 2`$-mixing matrix structure in the flavour space. Solving the full momentum dependent kinetic equations for these density matrices numerically  is a very complicated task and all attempts published to date have used some approximations to simplify the problem. Here we use momentum averaged approximation i.e. we set $`𝐏(p)𝐏\left(p\right)`$, with $`p3.15T`$. This approach has been found to give a very good approximation of the $`\nu _s`$ equilibration , and it will be sufficient for the purposes of this letter. The coupled equations are then (in the case of $`\nu _\tau \nu _s`$ oscillations, other cases can be obtained easily by simple redefinitions which are found for example in ):

$`\dot{𝐏}`$ $`=`$ $`𝐕\times 𝐏(D+{\displaystyle \frac{d}{dt}}\mathrm{log}P_0)𝐏_T+(1P_z){\displaystyle \frac{d}{dt}}\mathrm{log}P_0\widehat{𝐳},`$   
$`\dot{\overline{𝐏}}`$ $`=`$ $`\overline{𝐕}\times \overline{𝐏}(\overline{D}+{\displaystyle \frac{d}{dt}}\mathrm{log}\overline{P}_0)\overline{𝐏}_T+(1\overline{P}_z){\displaystyle \frac{d}{dt}}\mathrm{log}\overline{P}_0\widehat{𝐳},`$   
$`\dot{P_0}`$ $`=`$ $`\mathrm{\Gamma }^{}(\nu _\tau \overline{\nu }_\tau \alpha \overline{\alpha })\left(n_{eq}^2n_{\nu _\tau }n_{\overline{\nu }_\tau }\right),`$   
$`\dot{\overline{P}_0}`$ $`=`$ $`\overline{\mathrm{\Gamma }}(\nu _\tau \overline{\nu }_\tau \alpha \overline{\alpha })\left(n_{eq}^2n_{\nu _\tau }n_{\overline{\nu }_\tau }\right),`$ (2)   
where $`\dot{x}dx/dt`$ and $`𝐏_T=P_x\widehat{𝐱}+P_y\widehat{𝐲}`$. The damping coefficients for particles and anti-particles are $`D\overline{D}1.8G_F^2T^5`$ very accurately. The rotation vector $`𝐕`$ is

$$𝐕=V_x\widehat{𝐱}+\left(V_0+V_L\right)\widehat{𝐳},$$
(3)   
where

$`V_x`$ $`=`$ $`{\displaystyle \frac{\delta m^2}{2p}}\mathrm{sin}2\theta ,`$   
$`V_0`$ $`=`$ $`{\displaystyle \frac{\delta m^2}{2p}}\mathrm{cos}2\theta +\delta V_\tau ,`$   
$`V_L`$ $`=`$ $`\sqrt{2}G_FN_\gamma L,`$ (4)   
where $`\theta `$ is the vacuum mixing angle, $`\delta m^2=m_{\nu _s}^2m_{\nu _\tau }^2`$, $`N_\gamma `$ is the photon number density and the effective asymmetry $`L`$ in the potential $`V_L`$ is given by

$$L=\frac{1}{2}L_n+L_{\nu _e}+L_{\nu _\mu }+2L_{\nu _\tau }(P)=\eta +2L_{\nu _\tau }(P),$$
(5)   
in the case of an electrically neutral plasma. Asymmetry $`L_{\nu _\tau }`$ is obtained from

$$L_{\nu _\tau }=\frac{3}{8}\left(P_0(1+P_z)\overline{P}_0(1+\overline{P}_z)\right)$$
(6)   
and $`L_n`$ is the neutron asymmetry. The potential term $`\delta V_\tau `$ is approximately 

$$\delta V_\tau =17.8G_FN_\gamma \frac{pT}{2M_Z^2}.$$
(7)   
The rotation vector for anti-neutrinos is simply $`\overline{𝐕}(L)=𝐕(L)`$.

Neutrino and anti-neutrino ensembles are very strongly coupled in Eq. (2) through the effective potential term $`V_L(L)`$, which makes their numerical solution particularly difficult; as long as neutrino asymmetry remains small, there is a large cancellation in the Eq. (6), leading to a potential loss of accuracy. To overcome this problem we define new variables

$$P_\alpha ^\pm P_\alpha \pm \overline{P}_\alpha .$$
(8)   
In terms of these the Eq. (2) for $`𝐏^\pm `$ becomes

$`\dot{P}_x^\pm `$ $`=`$ $`V_0P_y^\pm V_LP_y^{}\stackrel{~}{D}P_x^\pm ,`$   
$`\dot{P}_y^\pm `$ $`=`$ $`V_0P_x^\pm +V_LP_x^{}\stackrel{~}{D}P_y^\pm V_xP_z^\pm ,`$   
$`\dot{P}_z^\pm `$ $`=`$ $`V_xP_y^\pm +A_\pm \left(2P_z^+\right)A_{}P_z^\pm ,`$ (9)   
where we have defined

$`A_+`$ $``$ $`{\displaystyle \frac{d}{dt}}\mathrm{log}P_0^+,`$   
$`A_{}`$ $`=`$ $`\mathrm{\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}0},`$   
$`\stackrel{~}{D}_{}`$ $`=`$ $`D+A_+.`$ (10)   
Finally, since for the averaged interaction rates $`\mathrm{\Gamma }=\overline{\mathrm{\Gamma }}`$, the difference $`P_0^{}`$ is not affected by collisions and we find

$`\dot{P}_0^+`$ $`=`$ $`2\mathrm{\Gamma }\left(n_{eq}^2n_{\nu _\tau }n_{\overline{\nu }_\tau }\right),`$   
$`\dot{P}_0^{}`$ $`=`$ $`0.`$ (11)   
Our objective is to study whether the sign of the final asymmetry, sign($`L`$) , follows deterministically from the initial conditions. In  it was found to be chaotic, whereas the authors of  claim that sign($`L`$) is fully deterministic and equal to the sign of the initial neutrino asymmetry. Similar results have been reported by other groups as well .

The key ingredient in the physics leading to the growth of the asymmetry is the appearance of the resonance. Indeed, if the squared mass difference $`\delta m^2<0`$, the effective potential $`V_0\pm V_L`$ goes through zero at the resonance temperature

$$T_c16.0(|\delta m^2|\mathrm{cos}2\theta )^{1/6}\text{MeV},$$
(12)   
where $`V_L0`$ is assumed, as effective asymmetry $`L`$ is driven to zero well before the resonance. After the resonance the balance of the system abruptly changes which leads to a rapidly growing $`L`$. We have numerically solved Eqs. (9) and (11), and examples of the results are shown in Fig. 1.

The behaviour of the system can be understood by a simple analogy of a ball rolling down a valley. After the resonance temperature $`T_c`$ the originally stable valley at $`L=0`$ becomes a ridge line separating two new, degenerate valleys corresponding to solutions of $`V_0\pm V_L=0`$. The system may first oscillate from one valley to another passing over the ridge line. However, because of friction (represented by damping terms) it will eventually settle into one or the other of the new valleys. It is easy to imagine that when there are many oscillations, even a very small difference in initial conditions may grow to a large phase difference at the time of settling down. In Fig. 1 these effects are demonstrated for two initial values: $`\eta =10^{10}`$ (solid line) and $`\eta =2\times 10^{10}`$ (dashed line) and oscillation parameters $`\delta m^2=10^2\text{eV}^2`$ and $`\mathrm{sin}^22\theta =10^{7.5}`$. For these parameters the resonance occurs at temperature $`T_c=7.41\text{MeV}`$.

Let us point out that change of sign($`L`$) can be effected either by variations in the oscillation parameters $`\delta m^2`$ and $`\mathrm{sin}^22\theta `$ or variations in the initial asymmetry $`\eta `$ as was shown in . These two cases are physically quite different of course. Varying oscillation parameters changes the shape of valleys forming after the resonance. This is an important issue because there will always be some experimental uncertainty in the measurements of masses and mixings, leading to unavoidable uncertainty of SBBN predictions in this scenario . Variations in $`\eta `$ are physical, for example due to local inhomogenieties in the baryon asymmetry created by early phase transitions, and correspond to deviations in the initial conditions (speed and direction of the ball) upon entering the resonance region.

In a recent article by Dolgov *et.al.*, results were presented which are in sharp disagreement with ours. In particular, it was claimed that the sign of neutrino asymmetry is completely deterministic. Moreover, they suggested that the chaotic behavior seen in  would be due to accumulated errors in the numerical codes. This suggestion is not correct. The numerical errors are completely under control in our computations. We have checked that allowing for the local error tolerance magnitudes larger than actually employed to get our results, have no effects on them. Rather we will now argue that the disagreement is due to ill-justified analytic approximations used in  to simplify the problem, which lead to artificial sign determinacy.

First approximation made in  was to neglect the term $`V_xP_z^{}`$ in the equation for $`P_y^{}`$ under the assumption that it is small in comparision to $`V_LP_y^+`$-term. However, just before the resonance effective asymmetry $`L`$ is driven to zero, so that effective potential $`V_L`$ is in fact very small, and $`V_xP_z^{}`$ should not be expected to be subdominant. We have indeed found that this term is of crucial importance for the initial asymmetry growth, as it prevents $`L`$ from getting arbitrarily small value before the resonance, as will be discussed below.

Second approximation made in  has even more dramatic effects. There it was argued that because Eq. (9) for $`P_x^{}`$ and $`P_y^{}`$ can be scaled to a form $`\dot{P}_{x,y}^{}=Q(a+b+c)`$, where $`Q5.6\times 10^4\sqrt{|\mathrm{cos}2\theta \delta m^2|}`$ is a large parameter, it should be safe to set the derivatives $`\dot{P}_{x,y}^{}`$ to zero. The rapid oscillations seen in our solutions indicate this approximation breaks down at the resonance temperature.

To study the effect of these approximations quantitatively, we introduce them into our equations. Dropping the term $`V_xP_z^{}`$ and setting $`\dot{P}_{x,y}^{}`$ to zero in the equations (9) we then find the constraints

$`0`$ $`=`$ $`V_0P_y^{}V_LP_y^+\stackrel{~}{D}P_x^{},`$   
$`0`$ $`=`$ $`V_0P_x^{}+V_LP_x^+\stackrel{~}{D}P_y^{}.`$ (13)   
From these equations, one can solve the evolution of $`P_x^{}`$ and $`P_y^{}`$ algebraically with the result

$`P_x^{}`$ $`=`$ $`{\displaystyle \frac{V_L}{V_0^2+\stackrel{~}{D}^2}}\left(V_0P_x^++\stackrel{~}{D}P_y^+\right)`$   
$`P_y^{}`$ $`=`$ $`{\displaystyle \frac{V_L}{V_0^2+\stackrel{~}{D}^2}}\left(V_0P_y^++\stackrel{~}{D}P_x^+\right).`$ (14)   
The remaining variables in Eqs. (9) and (11) are then solved numerically. In Fig. 2 we plot the results of a computation with (solid line) and without (dashed line) the implementation of the constraints (14) for same parameters as in Fig. 1. The constrained solutions which fall on top of each others in the figure are indeed fully deterministic and display no oscillation. This is to be expected because the solutions (14) cannot produce sign changing oscillations in $`L`$, since when $`L`$ goes to zero so does $`P_y^{}`$ and hence $`\dot{P}_z^{}`$ is then strongly suppressed in the constrained case. This prohibits $`L`$ from changing sign, so that sign($`L`$) is decided by the initial value of $`L`$. When the evolution of the off-diagonal components is neglected, the mechanism of the asymmetry growth is different, and instead of oscillating between the two new valleys, the system slowly rolls down to one of them after the resonance. This is in fact common characteristic to so called static approximations , where the differences between off-diagonals in the neutrino and antineutrino density matrices are ignored.

Let us now consider in detail how the complicated interplay between the various terms in the Eq. (9) leads to the initial exponential growth and oscillations of the asymmetry at the resonance temperature. Before the resonance $`𝐏^+`$ components are practically independent of $`𝐏^{}`$ components as effective asymmetry $`L`$ is very small. However, off-diagonal components $`P_x^+`$ and $`P_y^+`$ do grow to fairly large values. When $`V_0`$ changes sign at the resonance $`P_x^+`$ is rapidly driven to change sign, while the evolution of $`P_y^+`$ remains unaffected due to the additional term $`V_xP_z^+`$. $`P_z^+`$ stays near its initial value until the off-diagonals begin to grow, after which it begins to diminish, signalling the sterile neutrino production. Overall the evolution of $`𝐏^+`$ components is smooth (see Fig. 3), and the noticeable direct effect of the resonance is the changing of $`P_x^+`$.

The evolution of $`𝐏^{}`$ components is more complicated. Before the resonance neutrino and anti-neutrino ensembles follow each other closely in the sense that off-diagonals $`P_{x,y}^{}`$ are small. The main force which is keeping the difference of neutrino and antineutrino off-diagonals small is the potential $`V_L`$ and not the effect of damping terms. This can be confirmed by explicitly plotting individual terms appearing in the derivatives $`\dot{P}_{x,y}^{}`$ (see Fig. 4). The underlying mechanism driving $`\dot{P}_{x,y}^{}`$ close to zero before the resonance is then the cancellations between the effective potential terms. The $`V_xP_z^{}`$ term in Eq. (9) prevents $`L`$ from stabilizing to an arbitrarily small value, since when $`L`$ and off-diagonal components $`P_{x,y}^{}`$ are driven towards zero before the resonance, the small difference $`P_z^{}`$ becomes important. The value of $`V_L`$ together with differences $`P_{x,y}^{}`$ will cancel the effect of a non-zero $`P_z^{}`$. In this way the terms $`V_0P_y^{}`$ and $`V_LP_y^+`$ cancel each other in the equation for $`P_x^{}`$, and $`V_0P_x^{}`$ together with $`V_LP_x^+`$ cancel $`V_xP_z^{}`$ in the equation for $`P_y^{}`$.

After $`V_0`$ changes sign the cancellation of the potentials in $`\dot{P}_x^{}`$ equation will not work, which turns around the effect of $`V_L`$ in the $`P_x^{}`$ equation. Moreover, as $`P_x^+`$ changes sign due to the resonance, the effect of $`V_L`$ changes also in the $`P_y^{}`$ equation. The difference between both off-diagonal components is now growing due to a non-zero value of $`L`$ leading to a self supporting exponential growth of $`L`$. We wish to stress that the magnitudes of all the effects discussed here are well above the numerical accuracy. It is true that before the resonance $`L`$ goes very close to zero, but this is not relevant. What is relevant instead is that $`P_z^{}`$ remains at a value of order $`\eta `$.

It is interesting to note that the dynamics explained above will always lead to an initial growth of $`L`$ into the direction given by the sign of $`\eta `$. This results agrees with the analytical considerations in . Nevertheless, as we have seen, this does not guarantee that the final sign is that of $`\eta `$. One should bear in mind that the beginning of the asymmetry growth is a very complicated phenomenon, where all the variables and almost all the terms in the evolution equations are important to the outcome of the resonance and a complete account of the dynamics of the off-diagonal elements is of crucial importance. This makes it of course very difficult, if not impossible, to find any sensible analytical approximation to Eq. (9).

There still remains one contradiction with respect to the results of  and the constrained solution obtained here; namely, as seen from Fig. 2, even the constrained solutions yield $`L`$ which grows to a large value. The quantitative study of this effect is beyond the validity of the method employed here. We may however point out a further approximation made in , which appears to be the reason why they see a much weaker asymmetry growth. The collision term in  for the active-active component, $`\rho _{aa}`$ of the density matrix (1) is of the form

$$\mathrm{\Gamma }(p)(\rho _{aa}(p)f_{eq}),$$
(15)   
where in  $`f_{eq}`$ was taken to be the free Fermi distribution function $`f_{eq}=[1+\mathrm{exp}(p/T)]^1`$. The authors in  note themselves that $`f_{eq}`$ should actually be the distribution function which includes a chemical potential, $`f_{eq}(\mu )=[1+\mathrm{exp}(p/T\mu /T)]^1`$, but they assume the difference to be minor because the $`\mu /T`$ is small. However, one can show that $`\mu /T0.7L_{\nu _\tau }`$ and expand $`f_{eq}(L_{\nu _\tau })`$ to give

$$f_{eq}(L_{\nu _\tau })f_{eq}(0)\pm 0.7L_{\nu _\tau }\frac{f_{eq}(0)}{1+\mathrm{exp}(p/T)}.$$
(16)   
Approximating the second term of the expansion to zero in (15) then means that the system is seeking the free Fermi distribution instead of the correct equilibrium distribution which includes an asymmetry. This gives rise to an artificial force proportional to $`\mathrm{\Gamma }(p)L_{\nu _\tau }`$ resisting the growth of the asymmetry. Taking this into account, together with their previous assumptions leading to constraints (14), according to which the rate of asymmetry generation itself is weak

$$\dot{P}_z^{}P_y^{}V_LL,$$
(17)   
it appears likely that this force is able stop the asymmetry growth before the system has reached the true bottom of the valley. We chose not to pursue this issue further, since even the previous approximations has rendered the system unphysical by denying the possibility of chaoticity.

In this letter we have considered the indeterminacy or chaoticity in the sign of the neutrino asymmetry $`L_{\nu _\tau }`$ arising from $`\nu _\tau \nu _s`$ mixing in the early universe. We carefully discussed the dynamics giving rise to the growth of asymmetry, and unravelled the mechanism leading to the uncertaintity in sign($`L`$) . We confirmed that in the region of the parameter space identified in  the system is very sensitive to small variations in the initial conditions and mixing parameters. We have carefully checked that our numerical methods are higly accurate, so that the effects of the variations leading to a sign indeterminacy in the final asymmetry are physical.

As we pointed out, our results are in contradiction with the recent claims , according to which the sign of asymmetry is completely deterministic and the asymmetry growth is small compared to the results of . We have showed that the results of  are in fact unphysical and that they arise because of oversimplyfying approximations which artificially stabilize the dynamics responsible for the chaoticity.

Indeed, we found that all the terms in the evolution equations and all the components of the density matrix are important for the dynamics of the system. Hence it appears unlikely that any simplifying analytic approximation can be found that would describe the system adequately. We have employed the momentum averaged equations, as they are sufficient to study the chaoticity of the asymmetry growth and to resolve the validity of approximations imposed in . Our unpublished results with a momentum dependent code support the results of this letter, as well as do the results obtained by other groups working with momentum dependent equations , although the chaotic region appears to be somewhat smaller, when momentum dependence is included.

Acknowledgements

The author wishes to thank Kari Enqvist and Kimmo Kainulainen for enjoyable collaboration on issues related to this letter, and in particular KK for suggesting this research and his patient help in preparation of this manuscript. It is a pleasure to acknowledge the hospitality of NORDITA, Copenhagen, during the time this work was completed.