# 1 fig2

Active-sterile neutrino oscillations in the early Universe is a fascinating possibility with far-reaching consequences e.g. for nucleosynthesis \[1-12\] and CMB radiation . Nucleosynthesis considerations in particular have made it possible to place stringent constraints on model building aimed at understanding the observed neutrino anomalies in terrestrial observations. In the very first papers  it was observed that the mixing with an active species ($`SU(2)`$-doublet) endows the sterile ($`SU(2)`$-singlet) neutrino with effective interactions, which can be strong enough to bring the sterile species in equilibrium. The ensuing excess energy density would result in a failure of the nucleosynthesis explanation of the observed light element abundances . This line of reasoning was put to a solid computational foundation in refs. , and these results were later reproduced in ref. . Of particular interest was the observation that nucleosynthesis is in conflict with $`\nu _\mu \nu _s`$-oscillation solution to the atmospheric neutrino problem .

Already in  it was noted that nucleosynthesis constraints  depend on the reasonable assumption that the leptonic asymmetries are not many orders of magnitudes larger than the baryonic one. More specifically, e.g. for a mass squared difference $`|\delta m^2|=10^4\text{eV}^2`$ a large initial asymmetry, $`L_\nu ^{\mathrm{in}}>\mathrm{\hspace{0.33em}10}^5`$ (here $`L_\alpha =(N_\alpha N_{\overline{\alpha }})/N_\gamma `$) would suppress the effective mixing angle so much that the equilibration would never take place. This observation was later revived by Foot and Volkas , who basing on this and another previously observed effect, an exponential growth of leptonic asymmetry , suggested an interesting way to circumvent the nucleosynthesis constraints without invoking unnatural initial conditions . Their scenario assumes a novel mass-mixing scenario, where a $`\nu _\tau \nu _s`$-mixing, with carefully chosen parameters, produces a large leptonic asymmetry (but does not equilibrate $`\nu _s`$), which suppresses the subsequent $`\nu _\mu \nu _s`$-mixing angle and thereby prevents the $`\nu _s`$-equilibration from taking place. Some details concerning the growth of the asymmetry in this scenario are still under debate .

It was later observed by Shi  that the period of exponential growth exhibits chaotic features and therefore, while the amplitude of the final asymmetry is robust, its sign appeared to be essentially arbitrary. This raises some interesting questions: for example, is the sign of $`L_\nu `$ sensitive to the fluctuations in the initial conditions, like in the baryon asymmetry? If so, one should expect a large suppression in the effective asymmetry present at the important epoch for the $`\nu _\mu \nu _s`$-oscillations due to diffusion effects. For this purpose it is important to establish the extent of possible chaotic or regular regions in the parameter space. In this letter we have studied the dependence of the sign of $`L_\nu `$ on neutrino mixing parameters. We find a rather clearcut division of the parameter space into non chaotic and partly or completely chaotic<sup>1</sup><sup>1</sup>1We do not claim here that the system exhibits chaoticity in the mathematical sense of the definition of chaos; we merely mean that the system is sensitive to small variations of the parameters. regions. In chaotic regions the final sign of the asymmetry is indeed found to be highly sensitive also to fluctuations in the initial conditions.

Another, more direct consequence follows from the fact that sign($`L`$) affects the computed <sup>4</sup>He abundance, either directly in the case of $`\nu _e\nu _s`$ oscillations, or when induced by large $`L_\nu `$ created in $`\nu _\mu \nu _s`$ or $`\nu _\tau \nu _s`$ oscillations and later transferred to $`\nu _e`$-sector via active-active oscillations  . It then follows that possible chaotic behavior will constrain our chances to draw any definite conclusions about the effects of sterile neutrinos on Big Bang nucleosynthesis, as we will discuss below.   

In the early Universe neutrinos experience frequent scatterings, which tend to bring their distributions into thermal equilibrium. The requisite mathematical formalism is therefore very different from the one particle approach valid for description of accelerator physics (beams) and even solar neutrinos. Indeed, the objects of interest are the (reduced) density matrices for the neutrino and antineutrino ensembles

$$\rho _\nu \frac{1}{2}P_0(1+𝐏),\rho _{\overline{\nu }}\frac{1}{2}\overline{P}_0(1+\overline{𝐏}).$$
(1)   
Solving full momentum dependent kinetic equations for $`\rho _\nu (p)`$ and $`\rho _{\overline{\nu }}(p)`$  is obviously a very difficult task. Instead, we employ the momentum averaged equations for $`𝐏=𝐏(p)`$, with $`p3.15T`$, which should be expected to give a good approximation for the full system . (Our preliminary studies with full momentum dependent equations support this assumption). Moreover, for the parameters we are interested in, one can neglect the collision terms so that $`P_0`$ remains a constant and can be set to a unity. The coupled equations of motion then are (for definiteness we shall focus here on $`\nu _\tau \nu _s`$ oscillations; other cases are obtained from this by simple redefinitions<sup>2</sup><sup>2</sup>2Interested reader can find these redefinitions for example from .)

$`\dot{𝐏}`$ $`=`$ $`𝐕\times 𝐏D𝐏_T`$   
$`\dot{\overline{𝐏}}`$ $`=`$ $`\overline{𝐕}\times \overline{𝐏}D\overline{𝐏}_T`$ (2)   
where $`\dot{𝐏}d𝐏/dt`$ and we defined $`𝐏_TP_x\widehat{𝐱}+P_y\widehat{𝐲}`$. In the case of $`\nu _\tau \nu _s`$ oscillations the damping coefficient is $`D1.8G_F^2T^5`$ , and $`\overline{D}D`$ to a very high accuracy. It is convenient to decompose the rotation vector $`𝐕`$ as

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

$`V_x`$ $`=`$ $`\mathrm{\Delta }\mathrm{sin}2\theta `$   
$`V_0`$ $`=`$ $`\mathrm{\Delta }\mathrm{cos}2\theta +\delta V_\tau `$   
$`V_L`$ $`=`$ $`\sqrt{2}G_FN_\gamma L,`$ (4)   
where $`\theta `$ is the vacuum mixing angle, $`\mathrm{\Delta }\delta m^2/2p`$ and the photon number density $`N_\gamma 2\zeta (3)T^3/\pi ^2`$. The effective asymmetry $`L`$ appearing in the leading contribution $`V_L`$ to the neutrino effective potential is

$$L=\frac{1}{2}L_n+L_{\nu _e}+L_{\nu _\mu }+2L_{\nu _\tau }(P)$$
(5)   
where $`n`$ refers to neutrons and we have assumed electrical neutrality of the plasma. The remaining piece to the effective potential $`\delta V_\tau `$ is given by 

$$\delta V_\tau =\sqrt{2}G_FN_\gamma A_\tau \frac{pT}{2M_Z^2},$$
(6)   
with $`A_\tau =14\zeta (4)/\zeta (3)12.61`$. The rotation vector for antineutrinos is simply $`\overline{𝐕}(L)=𝐕(L)`$. The coupling of particle and antiparticle sectors occurs through the asymmetry term, where

$$L_{\nu _\tau }(P)=L_{\nu _\tau }^{\mathrm{in}}+\frac{3}{8}(P_z\overline{P}_z).$$
(7)   
with $`L_{\nu _\tau }^{\mathrm{in}}`$ being the initial $`\nu _\tau `$ asymmetry.

Even the simplified one state-quantum kinetic equations (2) are very difficult to handle numerically, because of the vast difference in the time scales involved (Hubble expansion rate, matter oscillation frequency and the width of the resonance, for example) on one hand and due to extremely strong coupling induced by the asymmetry term on the other. The so-called static approximation employed in  reduces the system to one first order differential equation. Unfortunately it is not really suitable for the treatment of oscillations at the resonance, since many of its basic assumptions – that the system is adiabatic, that the MSW-effect can be neglected, and that the rate of change of lepton number is dominated by the collisions – break down at the resonance.

These considerations emphasize the need for a very careful numerical approach. In practice, the accuracy is much improved if one makes separation between the large ($``$ number density) and small components ($``$ asymmetry) in equations (2). To this end we change the variables into

$$P_i^\pm P_i\pm \overline{P}_i,$$
(8)   
in terms of which (2) become

$`\dot{P}_x^+`$ $`=`$ $`V_0P_y^+V_LP_y^{}DP_x^+`$   
$`\dot{P}_y^+`$ $`=`$ $`V_0P_x^++V_LP_x^{}DP_y^+`$   
$`\dot{P}_z^+`$ $`=`$ $`V_xP_y^+`$   
$`\dot{P}_x^{}`$ $`=`$ $`V_0P_y^{}V_LP_y^+DP_x^{}`$   
$`\dot{P}_y^{}`$ $`=`$ $`V_0P_x^{}+V_LP_x^+DP_y^{}`$   
$`\dot{P}_z^{}`$ $`=`$ $`V_xP_y^{}.`$ (9)   
We have studied numerically the behaviour of the system described by (9) as a function of the oscillation parameters $`\delta m^2`$ and $`\mathrm{sin}^22\theta `$, and in particular the evolution of the asymmetry (7). Of crucial importance in this evolution is the occurrence of the resonance at $`V_0=0`$ if $`\delta m^2<0`$. Inserting the appropriate parameters, one finds that the resonance temperature is given by 

$$T_{\mathrm{res}}16.0(|\delta m^2|\mathrm{cos}2\theta )^{1/6}\mathrm{MeV},$$
(10)   
where $`\delta m^2`$ is given in units $`\mathrm{eV}^2`$. Far above the resonance the damping terms tend to suppress the off-diagonal elements $`𝐏_T`$ and moreover, the system is driven towards the initially stable fixed point $`L=0`$. As soon as the system passes the resonance however, $`L=0`$ becomes an unstable fixed point and two new locally stable and degenerate minima corresponding to the solutions of $`V_L+V_0=0`$ appear; these are given by the condition $`|L|11(|\delta m^2|/\mathrm{eV}^2)(\mathrm{MeV}/T)^4`$.

The system is roughly analogous to a ball rolling down a valley that branches to two, passing via a saddle configuration. Initially the branching into these two new valleys can be very shallow and it may stay that for a long time. Once on the side of the bifurcation, the ball still keeps on passing over the central barrier (the continuation of the old stable fixed point to an unstable extremum) until the barrier grows too high, or until friction (damping) reduces energy enough, and it gets trapped to one of the valleys. It is easy to picture in ones mind that in a case of a very shallow bifurcation a small change in the initial conditions ($`L^{\mathrm{in}}`$), or in the shape of the valleys (oscillation parameters), can very much affect which minimum the system finally chooses to settle in.   

In Fig. 1 we show the asymmetry $`L_{\nu _\tau }`$ (7) resulting from solving equations (9) for two sets of parameters. In Fig. 1a the resonance is rather narrow and only one oscillation occurs before the system is trapped into the minimum with a negative sign of $`L_{\nu _\tau }`$. The subsequent oscillations about this new local minimum are quickly washed away by the damping terms. In contrast, Fig. 1b shows an example of oscillation parameters with which the bifurcation into new local minimum is extremely slow and for a long time there is hardly any barrier between the two minima with opposite signs of $`L_{\nu _\tau }`$, and the system oscillates thousands of times before settling down to a minimum with positive $`L_{\nu _\tau }`$. After settling down, the further evolution of the asymmetry follows a power-like behaviour. These results agree well with those of ref. .

It is instructive to look a little more carefully into how the system approaches the resonance. Before the resonance the off-diagonal $`P^\pm `$ components are very near zero and $`P_z^{}`$ near the value $`L=0`$. Just before the resonance the $`P_{x,y}^+`$ components begin to increase, which triggers both the growth of $`P_{x,y}^{}`$ components and the decrease of $`P_z^+`$. As $`V_0`$ changes the sign at the resonance it creates an instability in the equation for $`P_y^{}`$, which eventually strongly pushes $`P_y^{}`$ to negative direction. The simple coupling of $`P_y^{}`$ to $`P_z^{}`$ in (9) then drags $`P_z^{}`$ along leading to a rapid growth of $`L_{\nu _\tau }`$. So far these phenomena have not much affected the evolution of $`P^+`$ variables (which have continued to grow). Eventually however, the exponential growth of $`P^{}`$ terms causes the $`V_L`$ term in the equations for $`P_x^+`$ and $`P_y^+`$, insofar neglible, become dominant. Large $`V_L`$ then forces $`P_x^+`$ and $`P_y^+`$ to change sign and grow to opposite direction until $`V_L`$ again changes the sign. Additionally, the ensuing oscillatory motion of the $`P_{x,y}^+`$ components induces the oscillation into other variables as well, leading to the exponentially large oscillation pattern observed in $`L_{\nu _\tau }`$.   

To find out the extent of chaotic and/or regular behaviour of sign($`L`$), we have scanned through the parameter space depicted in Fig. 2, which shows the sign of the final asymmetry $`L`$ with the initial value $`L^{\mathrm{in}}=10^{10}`$. As can be seen, the structure of sign($`L`$) is highly complex. In the upper left hand corner, extending downwards to large $`\theta `$, there is a regular region with no change in the asymmetry. Its existence is relatively easy to understand: this is the region where only the very first oscillation is carried out before the sign of the asymmetry is fixed. Since the direction of the first oscillation is determined by the sign of the initial asymmetry $`L^{\mathrm{in}}`$ (not necessarily the initial $`\nu _\tau `$-asymmetry), the sign of final neutrino asymmetry in this part of the parameter space should indeed be regular and fully determined.

The bands seen in the left hand side of Fig. 2 are formed as the system goes through two or more oscillations. In this region the number of oscillations is slowly increased as $`\theta `$ grows leading to less determined sign($`L`$) but it still can hardly be described as chaotic yet.

In addition to the two more or less regular regions there are regions where sign($`L`$) appears to be chaotic. The interval $`10^2<|\delta m^2|<\mathrm{\hspace{0.33em}1}`$ contains a very complicated structure. For $`|\delta m^2|>\mathrm{\hspace{0.33em}1}`$ one may discern some tendency for positive $`L`$ to prevail, while the region with $`|\delta m^2|<\mathrm{\hspace{0.33em}10}^2`$ appears to be pure white noise. In the lower right-hand corner of Fig. 2 (below the gray line), with $`|\delta m^2|/\mathrm{sin}^22\theta <\mathrm{\hspace{0.33em}10}`$, oscillations in $`L`$ will continue past the neutrino freeze-out and will not settle into any definite value. The boundary of the regular region above which $`L`$ is positive is given approximately by

$$\left(\mathrm{log}\frac{\delta m^2}{10^{5.8}\text{eV}^2}\right)^20.44\left(\mathrm{log}\frac{\mathrm{sin}^22\theta }{10^{4.5}}\right)^21.$$
(11)   
We conjecture that the chaotic behaviour occurs only when the oscillating period is long. We have also explored in finer detail a restricted region of parameters in the chaotic regime, without finding any structure. It is possible however, that the final sign of $`L`$ in this region is affected by the accumulated numerical error originating from the extremely high number of oscillation periods. In this sense proving a true chaoticity is of course not possible. Nevertheless, if the system is sensitive to numerical error, it should be expected to be sensitive to parameter fluctuations as well, so the general pattern of rapid sign($`L`$) fluctuations is expected to be robust. Finally, the region in the upper right-hand corner does not correspond to a large asymmetry, but it is merely the region where $`\nu _s`$ is fully equilibrated  and the absolute value of final $`L`$ is very small.

Changing the initial value for $`L^{\mathrm{in}}`$ does not change the picture qualitatively, although the structures evident in Fig. 2 shift slightly to the left when $`L^{\mathrm{in}}`$ is increased. Moreover, the changes saturate at $`L^{\mathrm{in}}>\mathrm{\hspace{0.33em}10}^9`$. These changes, or their absence, are very interesting however: In Fig. 3 we have plotted the value of $`L`$ at the temperature $`T=T_{\mathrm{res}}2.5`$, as a function of $`L^{\mathrm{in}}`$ for three representative choices of parameters. The first set, with $`\mathrm{sin}^22\theta =10^7`$ and $`\delta m^2=10^4`$ eV<sup>2</sup> corresponds to the stable region with positive $`L`$ in Fig. 2. As expected, $`L`$ remains positive independently of the initial value $`L^{\mathrm{in}}`$. In fact the dependence turned out to be smooth and linear showing that not only is the sign robust, but also that the numerical solution is well under control.

The second set, with $`\mathrm{sin}^22\theta =10^6`$ and $`\delta m^2=1`$ eV<sup>2</sup>, lies in the intermediate region where positive $`L`$ predominates, and the same dominance is seen as a function of the inital value $`L^{\mathrm{in}}`$. The last set with $`\mathrm{sin}^22\theta =10^6`$ and $`\delta m^2=10^3`$ eV<sup>2</sup> corresponds to the chaotic region. It is evident that sign($`L`$) is very sensitive to initial conditions, displaying clear randomness as a function of $`L^{\mathrm{in}}`$.   

The final value of sign($`L`$) has consequences for both atmospheric neutrinos and primordial nucleosynthesis. It has been proposed that Super-Kamiokande results for atmospheric neutrinos, which lie in the forbidden zone , might still allow a active-sterile mixing solution if the asymmetry growth is taken into account . Although the oscillation parameters in the case of atmospheric neutrinos are in the region where asymmetry growth is not expected, it has been argued that other neutrino oscillations could induce a large asymmetry in the active-sterile sector which the oscillations cannot damp.

If the outcome of neutrino oscillations is highly chaotic, the validity of such a scenario might be suspect. However, we found a large region in the parameter space where sign($`L`$) is very robust with respect to small variations of the mixing parameters. No chaoticity should be expected there with respect to other small perturbations, such as local perturbations in $`L^{\mathrm{in}}`$, either. It is in these stable domains where one would expect that the mechanism of ref.  can be successful.

In the region where sign($`L`$) is chaotic in the oscillation parameter space, it was also found to be sensitive to fluctuations in $`L^{\mathrm{in}}`$; these are predicted to be generated for example during the QCD phase transition, or in scenarios of electroweak baryogenesis . In such case causally disconnected regions would be expected to develop large asymmetries with a random sign distribution. It has been argued that then the nucleosynthesis constraint on active-sterile mixing would be even more stringent, because of additional MSW conversion taking place in the boundaries of domains with different sign($`L`$. However, our results indicate that the new constraints obtained in  may be overly optimistic, because for a large part of their excluded region we have found sign($`L`$) to be stable against small fluctuations; hence in no domain formation should be expected to occur in the first place.

Determining the sign of $`L`$ is important also for considering the effect of the electron (anti)neutrino spectrum distortions on the light element abundances . When the momentum spectrum gets distorted from its thermal equilibrium value the neutron to proton freezing ratio will change. Direct $`\nu _e\nu _s`$ oscillations obviously can induce such distortions, but also scenarios where large asymmetry is first generated in $`\nu _{\mu ,\tau }\nu _s`$ and then transferred to electron neutrino via $`\nu _e\nu _{\mu ,\tau }`$ oscillation, could have considerable effects on the electron neutrino spectrum. It turns out that positive sign($`L`$) has the effect of decreasing and negative sign($`L`$) of increasing <sup>4</sup>He abundance , so that the difference is $`\mathrm{\Delta }Y10^2`$, with some dependence on the oscillation parameters. Because the oscillation parameters cannot be measured with an arbitrary accuracy, it follows that in the region where the sign($`L`$) is chaotic, the role of resonant active-sterile neutrino mixing in Big-Bang Nucleosynthesis can not be reliably estimated. Rather, in this region, depicted in Fig. 2, there always remains an arbitrariness in the <sup>4</sup>He abundance given by $`\mathrm{\Delta }Y10^2`$, which should be considered as a source of systematic error.

In the region where the sign is stable, more concrete conclusions can be drawn. However, in this region $`L`$ is positive and only a rather small negative shift in the helium abundance $`0.005<\mathrm{\Delta }Y<0`$ was found  for these parameters. Interestingly enough, such a shift could ameliorate the apparent conflict of the nucleosynthesis theory viz-a-viz observations .

Our results in this paper are based on an averaged momentum description of the neutrino ensemble. Some effects, like the diffusion of the asymmetry between different momentum states, would seem to indicate the need for using full momentum-dependent kinetic equations. This is rather hard, since one has to deal with exponential growth in every momentum state and the width of the resonance is, for most of the parameter space, so small that one needs a very large number of bins to complete the task. Our preliminary results with momentum dependent kinetic equations support the results presented here.   

This work has been supported by the Academy of Finland under the contract 101-35224.