format: text

text:
# August 2000 Comment on “Neutrino oscillations in the early universe: how can large lepton asymmetry be generated?”

## I Introduction and overview

Active–sterile neutrino oscillations in the early universe can generate large differences in the number densities of relic neutrinos and antineutrinos .***For the historical record, the reader should note that Ref. contains the claim that the neutrino asymmetry is always small (less than $`10^7`$) and hence unimportant for all parameters of interest. This contrasts starkly with the later studies . The resulting “lepton asymmetries” have important phenomenological ramifications: the suppression of sterile neutrino production because of large matter effects , a modification to primordial Helium synthesis if a sufficiently large asymmetry is created for the electron neutrino prior to or during the Big Bang Nucleosynthesis epoch , with also the possibility to describe inhomogeneous BBN scenarios .

We focus in this Comment on two flavour active–sterile oscillations. Furthermore we restrict our discussion to that region of oscillation parameter space which induces large lepton asymmetry creation well before neutrino decoupling at $`T1`$ MeV. An important issue is the “final value of the asymmetry”, that is, the steady state value attained at the end of the dynamical evolution. The asymmetry is conveniently defined to be

$$L_{\nu _\alpha }=\frac{n_{\nu _\alpha }n_{\overline{\nu }_\alpha }}{n_\gamma },$$
(1)
where $`n_\psi `$ is the number density of species $`\psi `$, and $`\alpha `$ is one of $`e`$, $`\mu `$ or $`\tau `$. In Ref.  it is argued that the final value for $`L_{\nu _\alpha }`$ is of order $`0.20.37`$ for typical parameter choices. This claim is recently disputed by A. D. Dolgov, S. H. Hansen, S. Pastor and D. V. Semikoz (DHPS) in Ref. , who conclude that although the asymmetry can rise by five orders of magnitude above the baryon asymmetry level of $`10^{10}`$, it cannot reach the $`0.20.37`$ range found in Ref. . The purpose of this Comment is to explain the flaws in the DHPS analysis.

The basic reason turns out to be very simple: DHPS neglect the pure Mikheyev–Smirnov–Wolfenstein (MSW) effect . As explained in Ref., pure MSW transitions dominate lepton asymmetry evolution at lower temperatures, after the period of initial explosive $`L_{\nu _\alpha }`$ growth. The approximate evolution equations derived by DHPS have sufficient validity to describe the onset of asymmetry growth, but they are not a good approximation after the explosive amplification phase ends. Their oversight of the role played by MSW transitions also leads them to exaggerate the importance of a “back-reaction term” \[$`B_1`$ in Eq. (51) of DHPS\]. As we will explain, when this term is large, the MSW effect is also dominant so that their Eq. (51) is no longer valid. The size of the back-reaction is thus a moot point.

We now clarify the structure of our paper, while simultaneously providing a summary of the main points to be made. These remarks will be of interest to those readers who want to know the gist of our argument but not the technical details.

* Section II contains a quick survey of the Quantum Kinetic Equations (QKEs) governing active–sterile neutrino evolution . These “almost exact” equations are not at issue; the debate is on their correct solutions. The QKEs simultaneously incorporate coherent matter-affected evolution, decoherence due to non-forward neutrino scattering, and repopulation of active neutrino distribution functions from the background plasma.
* Section III reviews how the essential physics encoded in the QKEs can be revealed. A fundamental issue is collision-dominated versus coherent matter-affected oscillation-driven amplification. The former pertains to conditions prior to and during the onset of asymmetry growth, in which case the evolution is well described by the static approximation of Refs. , and by the closely related adiabatic limit approximation developed in Refs. . Equation (51) of DHPS is identical to the leading order term of the static approximation if the $`B_1`$ term is neglected, and is thus able to correctly predict the existence of a critical temperature  at which asymmetry evolution enters a brief explosive growth phase. Subsequent evolution, however, ceases to be collision dominated. In addition, a significant asymmetry, typically $`10^5`$ when explosive growth stops, now exists in the plasma; the antisymmetric Wolfenstein term in the matter potential is sufficiently large to modify the neutrino and antineutrino oscillation patterns in noticeably dissimilar ways. A vitally important consequence of this is the separation of the neutrino and antineutrino MSW resonance momenta . Depending on the sign of $`L_{\nu _\alpha }`$, one resonance momentum remains in the body of the (almost) Fermi–Dirac distribution, while the other is rapidly relegated to the tail. Figures 1b and 2b depict this phenomenon for two different oscillation parameter choices. The separation means that the MSW effect, now dominant over the collision mechanism, has very asymmetric consequences for neutrinos and antineutrinos . For $`L_{\nu _\alpha }>0`$, the resonance in the body is able to efficiently convert $`\overline{\nu }_\alpha `$’s into $`\overline{\nu }_s`$’s, while the other is impotent in the scarcely populated tail. This important physics was not understood by DHPS. Indeed, when collisions are absent, neutrino density matrix evolution is driven only by matter-affected propagation (and the expansion of the universe of course). Equation (51) of DHPS does not reduce to the coherent matter-affected evolution equations, so it does not contain the pure MSW effect . The static approximation, which yields $`\frac{dL_{\nu _\alpha }}{dt}=0`$ for zero collision rate, also neglects the MSW effect. Thus these equations are valid only in the collision dominated regime. Figures 1a and 2a illustrate our arguments. The solid and dash-dotted lines represent respectively numerical solutions to the exact QKEsDHPS are unable to successfully integrate the QKEs, and do not accept these curves as accurate, although the numerical procedures used, and their stability, are thoroughly discussed in Ref. . and the static approximation equations. Clearly, the latter describes the explosive growth phase excellently, but underestimates the subsequent growth of $`L_{\nu _\alpha }`$. At steady state, the discrepancy is significant. The dashed curve reveals the physics behind the QKEs in a very amusing way. The generating code compares the decoherence or interaction length $`\mathrm{}_{\mathrm{int}}`$ with the matter-affected oscillation length $`\mathrm{}_m`$ at every time step, with both lengths evaluated at resonance. If $`r|_{\mathrm{res}}(\mathrm{}_{\mathrm{int}}/\mathrm{}_m)|_{\mathrm{res}}`$ is less than one, the static approximation equations are used to forward the system in time. For $`r|_{\mathrm{res}}>1`$, pure adiabatic MSW evolution is employed. The quantity $`r|_{\mathrm{res}}`$ therefore measures the relative importance of collisions and MSW transitions. In the event $`r|_{\mathrm{res}}1`$, the MSW effect is damped since collisions seriously disrupt the coherence. The opposite case has $`r|_{\mathrm{res}}1`$, implying that collisions are relatively negligible. The dashed curve in Fig. 1a agrees with the QKEs in both the high and low $`T`$ regimes extremely well. There is a slight discrepancy immediately after explosive growth, since neither effect dominates at this stage — the very visible bump in Fig. 2a is an artifact of this “either/or” approximation in this intermediate regime.
* The DHPS equations underestimate the asymmetry after explosive growth even more seriously than does the static approximation, because of their $`B_1`$ back-reaction term. Section IV compares Eq. (51) of DHPS with the static approximation. As already noted, and acknowledged by DHPS, these equations agree to leading order. We show, however, that the $`B_1`$ term is large only when the MSW effect has taken over as the dominant driver of asymmetry evolution. Since DHPS Eq. (51) and the static approximation both neglect this effect, they should not be used when MSW transitions are important. (We will also identify where the MSW effect is unwittingly removed in their derivation.) Furthermore, the expansion parameter adopted by DHPS in their approximation scheme is not in fact the correct choice at lower temperatures when there is a significant asymmetry in the plasma. This misunderstanding led DHPS to erroneously claim validity for their equations at later stages in the evolution, especially the role of the back-reaction term which is but another casualty following from a faulty expansion parameter. In fact, the dramatic termination of the asymmetry growth due to this bogus term has its origin in the expression $`(1+x)^11x`$. The left hand side corresponds to the static approximation (with a minor discrepancy), while the right hand side denotes the DHPS result, with $`x`$ identified as the $`B_1`$ term. The artificially strong cut-off arises because the right hand side is used even when $`x1`$.
* Section V contains comments on other less important errors and some misleading statements made by DHPS, while Sec. VI is a conclusion.
* The Appendix provides a translation of DHPS’s notation into a more common lingua franca used by other authors. The rather obscure nomenclature of DHPS proves to be the biggest obstacle to understanding their analysis.

## II Short review of the Quantum Kinetic Equations

We consider the case of active–sterile two state mixing where the weak eigenstates $`\nu _\alpha `$ ($`\alpha =e,\mu `$ or $`\tau `$) and $`\nu _s`$ are linear combinations of two mass eigenstates $`\nu _a`$ and $`\nu _b`$,

$$\nu _\alpha =\mathrm{cos}\theta _0\nu _a+\mathrm{sin}\theta _0\nu _b,\nu _s=\mathrm{sin}\theta _0\nu _a+\mathrm{cos}\theta _0\nu _b,$$
(2)
with $`\theta _0`$ as the vacuum mixing angle, and $`\mathrm{cos}2\theta _0>0`$ by definition. The compact notation

$$s\mathrm{sin}2\theta _0,c\mathrm{cos}2\theta _0,$$
(3)
and the convention $`\delta m^2m_b^2m_a^2`$ will be used from now on.

The one-body reduced density matrix  for a $`\nu _\alpha \nu _s`$ system of momentum $`p`$ in the early universe can be parameterised through

$$\rho (y)=\frac{1}{2}[P_0(y)+𝐏(y)\sigma ],$$
(4)
in which $`𝐏(y)=P_x(y)\widehat{𝐱}+P_y(y)\widehat{𝐲}+P_z(y)\widehat{𝐳}`$ is the “polarisation vector”, and $`\sigma =\sigma _x\widehat{𝐱}+\sigma _y\widehat{𝐲}+\sigma _z\widehat{𝐳}`$, where $`\sigma _i`$ are the Pauli matrices. The quantity $`y`$ is a dimensionless momentum variable that does not red-shift by definition:

$$y\frac{p}{T}p(t),$$
(5)
where $``$ is the Friedmann–Robertson–Walker scale factor and $`T`$ is the temperature. It is understood that the density matrices and the variables $`P_i(y)`$ depend also on time $`t`$ or, equivalently, temperature $`T`$, which, for $`m_eTm_\mu `$, are related through $`dt/dT=1/(T)M_P/5.44T^3`$, where $`=\dot{}/`$ is the Hubble parameter and $`M_P1.22\times 10^{22}`$ MeV is the Planck mass.

We normalise the density matrices to give momentum distribution functions

$$f_{\nu _\alpha }(y)=\frac{1}{2}[P_0(y)+P_z(y)]f_{\mathrm{eq}}^0(y),f_{\nu _s}(y)=\frac{1}{2}[P_0(y)P_z(y)]f_{\mathrm{eq}}^0(y),$$
(6)
for $`\nu _\alpha `$ and $`\nu _s`$ respectively. The reference distribution $`f_{\mathrm{eq}}^0(y)`$ is of the Fermi–Dirac form \[see Eq. (12)\] with zero chemical potential. Similar expressions pertain to antineutrinos, with $`P_i(y)\overline{P}_i(y)`$. The number densities are obtained from the distribution functions through

$$n_\psi =\frac{T^3}{2\pi ^2}_0^{\mathrm{}}f_\psi (y)y^2𝑑y,$$
(7)
where $`\psi `$ denotes the species.

The evolution equations for $`𝐏(y)`$ and $`P_0(y)`$ are

$`{\displaystyle \frac{𝐏(y)}{t}}`$ $`=`$ $`𝐕(y)\times 𝐏(y)D(y)[P_x(y)\widehat{𝐱}+P_y(y)\widehat{𝐲}]+{\displaystyle \frac{P_0(y)}{t}}\widehat{𝐳},`$ (8)
$`{\displaystyle \frac{P_0(y)}{t}}`$ $``$ $`\mathrm{\Gamma }(y)\left\{{\displaystyle \frac{f_{\mathrm{eq}}(y)}{f_{\mathrm{eq}}^0(y)}}{\displaystyle \frac{1}{2}}[P_0(y)+P_z(y)]\right\}.`$ (9)
The damping or decoherence function $`D`$ is related to the total collision rate for $`\nu _\alpha `$ of momentum $`y`$ with the background plasma, $`\mathrm{\Gamma }(y)`$, via

$$D(y)=\frac{\mathrm{\Gamma }(y)}{2}\frac{1}{2}k_\alpha G_F^2T^5y,$$
(10)
with $`k_e1.27`$ and $`k_\mu =k_\tau 0.92`$ (for $`m_eTm_\mu `$).These values neglect the Pauli blocking factors. Including Pauli blocking factors lead corrections of order $`10\%`$ for these quantities. It is useful to define the quantity

$$\mathrm{}_{\mathrm{int}}D^1,$$
(11)
to be identified as the decoherence or interaction length.

The function $`f_{\mathrm{eq}}(y)`$ is the Fermi–Dirac distribution,

$$f_{\mathrm{eq}}(y)\frac{1}{1+e^{y\stackrel{~}{\mu }_\alpha }},$$
(12)
where $`\stackrel{~}{\mu }_\alpha \mu _{\nu _\alpha }/T`$ (and $`\stackrel{~}{\mu }_{\overline{\alpha }}\mu _{\overline{\nu }_\alpha }/T`$ for antineutrinos) is a dimensionless chemical potential which does not red-shift but depends on the neutrino asymmetry.

The matter potential vector , $`𝐕(y)=V_x(y)\widehat{𝐱}+V_z(y)\widehat{𝐳}`$, has coefficients

$$V_x=\frac{\mathrm{sin}2\theta _0}{\mathrm{}_0},V_z=\frac{\mathrm{cos}2\theta _0}{\mathrm{}_0}+V_\alpha ,$$
(13)
where $`\mathrm{}_02yT/\delta m^2`$ is the vacuum oscillation length. The quantity $`V_\alpha (y)`$ is the effective potential which takes the form

$$V_\alpha (y)=\pm \sqrt{2}G_Fn_\gamma L^{(\alpha )}\sqrt{2}G_Fn_\gamma \frac{A_\alpha T^2y}{M_W^2},$$
(14)
where the $`+()`$ sign corresponds to neutrino (antineutrino) ensemble, $`G_F`$ is the Fermi constant, $`M_W`$ is the $`W`$-boson mass, $`A_e17`$ and $`A_{\mu ,\tau }4.9`$ (for $`m_eTm_\mu `$), and $`L^{(\alpha )}`$ is the effective total lepton number (for the $`\alpha `$-neutrino species). This last quantity is essentially a sum of the individual asymmetries:

$$L^{(\alpha )}L_{\nu _\alpha }+L_{\nu _e}+L_{\nu _\mu }+L_{\nu _\tau }+\eta ,$$
(15)
with $`\eta `$ being a small term due to the asymmetry of the electrons and nucleons — typically of order $`|\eta |5\times 10^{10}`$, calculated from

$$\eta =L_e\frac{1}{2}L_N(\mathrm{for}\alpha =e),\eta =\frac{1}{2}L_N(\mathrm{for}\alpha =\mu ,\tau ).$$
(16)
For simplicity we shall henceforth write $`L`$ instead of $`L^{(\alpha )}`$.

We now introduce $`v_\alpha `$ and the dimensionless variables $`a`$ and $`b`$ by

$`v_\alpha `$ $``$ $`\mathrm{}_0V_\alpha =a+b,`$ (17)
$`a`$ $``$ $`\sqrt{2}G_Fn_\gamma L\mathrm{}_08.033Ly\left({\displaystyle \frac{\mathrm{eV}^2}{\delta m^2}}\right)\left({\displaystyle \frac{T}{\mathrm{MeV}}}\right)^4,`$ (18)
$`b`$ $``$ $`\sqrt{2}G_Fn_\gamma \mathrm{}_0{\displaystyle \frac{A_\alpha T^2y}{M_W^2}}=y^2\left({\displaystyle \frac{T}{T^\alpha }}\right)^6\left({\displaystyle \frac{\mathrm{eV}^2}{\delta m^2}}\right),`$ (19)
where, in the definition of $`v_\alpha `$, $`(+)`$ pertains to neutrinos (antineutrinos), and

$$T^\alpha 18.9(23.4)\mathrm{MeV},\mathrm{for}\alpha =e(\mu ,\tau ),$$
(20)
is a convenient reference temperature. In this new notation, the matter-affected oscillation length is given by

$$\mathrm{}_m\frac{1}{\sqrt{V_x^2+V_z^2}}=\frac{\mathrm{}_0}{\sqrt{s^2+\left(cv_\alpha \right)^2}}.$$
(21)
We will also have cause to consider the ratio

$$d\frac{\mathrm{}_0}{\mathrm{}_{\mathrm{int}}}=k_\alpha G_F^2\frac{T^6y^2}{\delta m^2}=\delta b$$
(22)
where $`\delta 0.8(2.0)\times 10^2`$ for $`\alpha =e(\mu ,\tau )`$. In addition, a very important role is played by

$$r\frac{\mathrm{}_{\mathrm{int}}}{\mathrm{}_m}=\frac{\sqrt{V_x^2+V_z^2}}{D}=\frac{\sqrt{s^2+(cv_\alpha )^2}}{d},$$
(23)
which is the ratio of the interaction and matter-affected oscillation lengths. This quantity parameterises the relative importance of decohering collisions and coherent oscillations on the evolution of the density matrices.

Note that the MSW resonance condition is given by $`V_z=0c=v_\alpha =a+b`$, where $`(+)`$ is used for neutrinos (antineutrinos). Evaluating $`r`$ at resonance,

$$r|_{\mathrm{res}}=\frac{V_x}{D}|_{\mathrm{res}}=\frac{s}{d}|_{\mathrm{res}},$$
(24)
we find that the relative importance of collisions and coherent oscillations is strongly dependent on $`s\mathrm{sin}2\theta _0`$.

To make our argument easier to follow, it is useful to distinguish two different asymptotic regimes for the resonance condition.

1. Negligible total lepton number $`|a||b|,c`$. This is the situation prior to the explosive growth phase, in which the resonance condition

$$c=b,$$
(25)
is identical for neutrinos and antineutrinos. Equation (25) can be recast in the form

$$T_{\mathrm{res}}^0(y)T^\alpha \left(\frac{|\delta m^2|}{\mathrm{eV}^2}\right)^{\frac{1}{6}}y^{\frac{1}{3}}c^{\frac{1}{6}},$$
(26)
where $`T_{\mathrm{res}}^0`$ is the temperature at which neutrinos or antineutrinos of momentum $`y`$ are at resonance. For the body of the Fermi–Dirac distribution, i.e., $`0.1y10`$, one can then deduce from this expression the temperature range in which oscillations will be most significant. Note, however, that since Eq. (14) holds only for $`T150\mathrm{MeV}`$, we are restricted to the regime $`|\delta m^2|/\mathrm{eV}^210^4`$ for $`c1`$.

Using the definition of $`T_{\mathrm{res}}^0`$ one may rewrite the function $`b`$ via

$$|b|=c\left(\frac{T}{T_{\mathrm{res}}^0}\right)^6.$$
(27)
This will be useful later. Equation (25) can also be rewritten in the form

$$y_{\mathrm{res}}^0=\left(\frac{T^\alpha }{T}\right)^3\left(\frac{|\delta m^2|}{\mathrm{eV}^2}\right)^{\frac{1}{2}}c^{\frac{1}{2}},$$
(28)
which shows that the resonance momentum moves from low to high values as the temperature decreases. At $`y=y_{\mathrm{res}}^0`$, the condition $`|a||b|=c`$ is equivalent to

$$L1.0(0.4)\times 10^6\left(\frac{|\delta m^2|}{\mathrm{eV}^2}\right)^{\frac{1}{3}}(y_{\mathrm{res}}^0)^{\frac{1}{3}}c^{\frac{1}{3}}L_t,$$
(29)
for $`\alpha =e(\mu ,\tau )`$. The most relevant choice for $`y_{\mathrm{res}}^0`$ is about 2, because asymmetry growth starts approximately at this point as the resonance momentum moves form low to high values with decreasing $`T`$. This value is approximately correct if sterile neutrino production has been negligible until that moment (which is true for small enough mixing angles).

If $`LL_t`$, the Wolfenstein term $`a`$ in the effective potential becomes dominant and two different resonances exist, one on either side of $`y_{\mathrm{res}}^0`$ in the momentum distribution. If $`L`$ is positive, the resonance at the lower (higher) value occurs only for antineutrinos (neutrinos), and vice versa if $`L`$ is negative. This “separation of the resonances” is crucial for understanding the mechanism of neutrino asymmetry generation.

2. Negligible finite-temperature term $`|b||a|,c`$. This situation holds after a substantial $`L`$ has been created. The resonance condition is now approximately given by

$$c=a,$$
(30)
where $`(+)`$ refers to neutrinos (antineutrinos). In the case $`L>0`$, the resonance occurs for antineutrinos (neutrinos) if $`\delta m^2<0`$ ($`\delta m^2>0`$). The opposite holds if $`L<0`$. It is important to recall that the asymmetry generation mechanism actually requires $`\delta m^2<0`$.

Equation (30) implies that the resonance momentum takes the value

$$y_{\mathrm{res}}^{\mathrm{low}}=y_{\mathrm{res}}^0\left[\frac{L_t}{|L|}\right]^{\frac{3}{4}}\frac{c}{8.033}\frac{|\delta m^2|}{\mathrm{eV}^2}\frac{1}{L(T/\mathrm{MeV})^4}.$$
(31)
Notice that $`y_{\mathrm{res}}^{\mathrm{low}}<y_{\mathrm{res}}^0`$ when $`|L|>L_t`$.

## III Understanding the lepton asymmetry evolution curves

The Quantum Kinetic Equations can be numerically integrated and lepton asymmetry evolution thereby extracted. In terms of the density matrices, the asymmetry is given by

$$L_{\nu _\alpha }=\frac{1}{2n_\gamma }\frac{T^3}{2\pi ^2}_0^{\mathrm{}}\left[P_0(y)+P_z(y)\overline{P}_0(y)\overline{P}_z(y)\right]f_{\mathrm{eq}}^0(y)y^2𝑑y.$$
(32)
Taking the time derivative of Eq. (32), and imposing $`\alpha +s`$ lepton charge conservation, i.e.,

$$_0^{\mathrm{}}\left[\frac{P_0(y)}{t}\frac{\overline{P}_0(y)}{t}\right]f_{\mathrm{eq}}^0(y)y^2𝑑y\frac{}{t}\left(\frac{n_{\nu _\alpha }+n_{\nu _s}n_{\overline{\nu }_\alpha }n_{\overline{\nu }_s}}{n_\gamma }\right)=0,$$
(33)
one easily deduces that

$$\frac{dL_{\nu _\alpha }}{dt}=\frac{1}{2n_\gamma }\frac{T^3}{2\pi ^2}_0^{\mathrm{}}V_x(P_y\overline{P}_y)f_{\mathrm{eq}}^0y^2𝑑y.$$
(34)
In principle, one could use Eq. (32) to directly calculate the asymmetry; Equation (34) is formally redundant. However, for very small asymmetries, the right hand side of the former equation is the difference of two numbers which are large relative to $`L_{\nu _\alpha }`$. Equation (34) should therefore be used to circumvent the numerical stability problems potentially posed by the direct procedure.

DHPS were apparently unable to achieve successful numerical integration of the QKEs for $`|\delta m^2|>10^7eV^2`$. We would guess that at least part of the reason for this was that they did not use Eq. (34).

The solid lines in Figs. 1a and 2a show our numerically stable asymmetry evolution curves from integrating the QKEs. There are four phases in the evolution: Above the critical temperature $`T_c`$, the asymmetry $`L_{\nu _\alpha }`$ evolves from its small initial value to the $`\eta 10^{10}`$ level such that $`L\eta `$. At $`T=T_c`$, explosive growth occurs for a short period, during which the asymmetry rises to roughly $`10^5`$ for the typical choice $`\delta m^21`$ eV2. After this phase, asymmetry growth continues at a more leisurely power law rate. One can measure from the graphs that $`L_{\nu _\alpha }T^4`$. Further evolution brings the asymmetry to a steady state value in the range $`0.20.37`$ depending on the oscillation parameters chosen.

The post explosive growth behaviour is disputed by DHPS. In this section, we review the physical explanation for the $`L_{\nu _\alpha }`$ curves obtained by solving the QKEs, and demonstrate concordance between these curves and those generated by certain approximate evolution equations which highlight the essential physical processes occurring in the various regimes.

### A Collision dominated evolution and the static approximation

The neutrino collision rate $`\mathrm{\Gamma }`$ increases as $`T^5`$. We expect, therefore, that at high enough temperatures, collisions are sufficiently frequent to seriously disrupt the coherent matter-affected oscillatory evolution of the ensemble. The finite temperature $`b`$ term in the effective potential also scales as $`T^5`$, and predominates at high $`T`$ if the initial neutrino asymmetry is small (say of the order of the baryon asymmetry). We focus first on this regime.

A useful and intuitively reasonable measure of the importance of collisional decoherence is afforded by the quantity $`r|_{\mathrm{res}}(\mathrm{}_{\mathrm{int}}/\mathrm{}_m)|_{\mathrm{res}}`$ defined in Eq. (24). The collision dominated regime then corresponds to $`r|_{\mathrm{res}}1`$. In this regime, one can develop a physical picture for distribution function evolution that yields approximate time-development equations for $`f_{\nu _\alpha }`$ and $`f_{\nu _s}`$ that reproduce the QKE behaviour exceedingly well. This technique has been called the “static approximation” . From another perspective, the resulting equations have been shown to arise in the adiabatic limit approximation to the QKEs in the large damping or collision dominated epoch . The approximate equations feature only the distribution functions, not the coherences $`P_{x,y}`$, so we are working in a Boltzmann limit.

We will not repeat the derivations in Refs. , but present just the results. Define:

$$z_{\alpha ,s}\frac{f_{\nu _{\alpha ,s}}}{f_{\mathrm{eq}}^0},\overline{z}_{\alpha ,s}\frac{f_{\overline{\nu }_{\alpha ,s}}}{f_{\mathrm{eq}}^0}.$$
(35)
The $`z^{}s`$ are the actual $`\nu _\alpha `$ and $`\nu _s`$ distribution functions normalised by the Fermi–Dirac distribution with zero chemical potential. In the static approximation, the sterile neutrino $`z`$’s obey the rate equations

$$\frac{z_s}{T}=\frac{\mathrm{\Gamma }_{\alpha s}}{T}\left[z_\alpha z_s\right],\frac{\overline{z}_s}{T}=\frac{\overline{\mathrm{\Gamma }}_{\alpha s}}{T}\left[\overline{z}_\alpha \overline{z}_s\right],$$
(36)
where

$$\mathrm{\Gamma }_{\alpha s}=\frac{1}{4}\frac{\mathrm{\Gamma }s^2}{s^2+d^2+(cb+a)^2},\overline{\mathrm{\Gamma }}_{\alpha s}=\frac{1}{4}\frac{\mathrm{\Gamma }s^2}{s^2+d^2+(cba)^2},$$
(37)
are the transition rates.

Using lepton number conservation and the following definitions:

$$z_s^+\frac{z_s+\overline{z}_s}{2},z_s^{}\frac{z_s\overline{z}_s}{2},$$
(38)
we see that the neutrino asymmetry evolution equation is given by

$$\frac{dL_{\nu _\alpha }}{dt}=\frac{1}{n_\gamma }\frac{T^3}{2\pi ^2}_0^{\mathrm{}}\frac{z_s^{}}{t}f_{\mathrm{eq}}^0y^2𝑑y.$$
(39)
From Eq. (36) the rate of change of $`z_s^{}`$ is obtained by solving the coupled equations

$`{\displaystyle \frac{z_s^+}{t}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }_{\alpha s}+\overline{\mathrm{\Gamma }}_{\alpha s}}{2}}[z_\alpha ^+z_s^+]+{\displaystyle \frac{\mathrm{\Gamma }_{\alpha s}\overline{\mathrm{\Gamma }}_{\alpha s}}{2}}\left[z_\alpha ^{}z_s^{}\right],`$ (40)
$`{\displaystyle \frac{z_s^{}}{t}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }_{\alpha s}\overline{\mathrm{\Gamma }}_{\alpha s}}{2}}[z_\alpha ^+z_s^+]+{\displaystyle \frac{\mathrm{\Gamma }_{\alpha s}+\overline{\mathrm{\Gamma }}_{\alpha s}}{2}}\left[z_\alpha ^{}z_s^{}\right].`$ (41)
Adopting the instantaneous repopulation approximation

$$z_\alpha =\frac{1+e^y}{1+e^{y\stackrel{~}{\mu }_\alpha }},\overline{z}_\alpha =\frac{1+e^y}{1+e^{y+\stackrel{~}{\mu }_\alpha }},$$
(42)
and expanding these expressions to first order in the chemical potential, we finally obtain

$$\frac{dL}{d(T/\mathrm{MeV})}=_0^{\mathrm{}}\frac{AL+BL_{\nu _\alpha }+B^{}}{\mathrm{\Delta }(s^2)}𝑑y,$$
(43)
with

$`A`$ $``$ $`1.27s^2k_\alpha \left({\displaystyle \frac{\mathrm{eV}^2}{\delta m^2}}\right)\left({\displaystyle \frac{T}{\mathrm{MeV}}}\right)^6(bc)y^4(1+e^y)^1(1z_s^+),`$ (44)
$`B`$ $``$ $`0.09s^2k_\alpha \left({\displaystyle \frac{T}{\mathrm{MeV}}}\right)^2[s^2+d^2+(cb)^2+a^2]y^3e^y(1+e^y)^2,`$ (45)
$`B^{}`$ $``$ $`0.03s^2k_\alpha \left({\displaystyle \frac{T}{\mathrm{MeV}}}\right)^2[s^2+d^2+(cb)^2+a^2]y^3(1+e^y)^1z_s^{},`$ (46)
where we have used $`n_\gamma =2\zeta (3)T^3/\pi ^2T^3/4.1`$, and

$$\mathrm{\Delta }(s^2)=[s^2+d^2+(cb+a)^2][s^2+d^2+(cba)^2].$$
(47)
Equations (40) and (43) form a coupled system that must be simultaneously solved to track asymmetry evolution and sterile neutrino production.

### B Repopulation, the damped MSW effect, and the adiabatic limit in the collision dominated epoch

References  show how the static approximation follows from the QKEs in the adiabatic limit. These results are briefly reviewed, with the focus on repopulation handling. We will see that repopulation is not a critical issue for asymmetry evolution in the collision dominated regime, although it does play an interesting role in the otherwise mundane $`T>T_c`$ epoch. Importantly, we explain how the static approximation excludes the MSW effect.

One may extract the collision dominated adiabatic limit beginning with the four coupled QKEs \[Eq. (8)\] . These calculations show that the proper incorporation of repopulation affects the evolution of $`L_{\nu _\alpha }`$ at sub-leading order only; an excellent description of the asymmetry evolution in the said limit may still be obtained from less rigorous repopulation treatments. We will not repeat here the somewhat involved algebra of Ref., but instead explain the gist of the result through a simpler argument. The key is to consider repopulation in two extreme limits: a full decoupling limit and a thermal equilibrium limit.

The first limit corresponds to completely switching off the refilling interactions, i.e., $`\frac{P_0}{t}=0`$ , while the other assumes $`f_{\nu _\alpha }=f_{\mathrm{eq}}`$, where the right hand side includes the appropriate chemical potential . In comparison, the latter limit is much closer to the truth since weak interaction induced scattering occurs frequently in the epoch under consideration.

In the full decoupling limit, the QKEs \[Eq. (8)\] reduce to the homogeneous system

$$\frac{𝐏}{t}=𝒦𝐏,$$
(48)
where the $`3\times 3`$ matrix $`𝒦`$ is given by

$$𝒦=\left(\begin{array}{ccc}D& V_z& 0\\ V_z& D& V_x\\ 0& V_x& 0\end{array}\right)=\frac{1}{\mathrm{}_0}\left(\begin{array}{ccc}d& (c\pm ab)& 0\\ (c\pm ab)& d& s\\ 0& s& 0\end{array}\right).$$
(49)
Rewriting Eq. (48) in the instantaneous diagonal basis, and taking the adiabatic limit by ignoring the time derivatives of the diagonalisation matrix , the evolution is now driven by the eigenvalues $`\lambda _{1,2,3}`$ of $`𝒦`$. In most cases, $`\lambda _{1,2}`$ form a complex conjugate pair. Since we are interested in the collision dominated regime which includes requiring that $`V_x/D1`$, these are well approximated by

$$\lambda _{1,2}D\pm i\sqrt{V_x^2+V_z^2}=\mathrm{}_{\mathrm{int}}^1\pm i\mathrm{}_m^1.$$
(50)
The remaining eigenvalue $`\lambda _3`$ is small, real and negative, and obeys

$$\lambda _3=\frac{DV_x^2}{(D\lambda _3)^2+V_x^2+V_z^2}\frac{DV_x^2}{D^2+V_x^2+V_z^2}=D\frac{\left(\frac{V_x}{D}\right)^2}{1+\frac{V_x^2+V_z^2}{D^2}}.$$
(51)
Since the real parts of $`\lambda _{1,2}`$ are negative and much larger in magnitude than $`\lambda _3`$, the solution lies in the direction of the eigenvector associated with $`\lambda _3`$, while the transverse components are strongly damped (which is another way to define collision dominance ).

It is very important to observe that the oscillatory aspect of the evolution, driven by the imaginary parts of $`\lambda _{1,2}`$, is severely damped in the regime concerned. This is precisely where the pure MSW effect is removed from the dynamics. The DHPS procedure also contains a step which nullifies MSW transitions, a point they failed to notice. Note that the ratio of the imaginary to the real parts of $`\lambda _{1,2}`$ is approximately $`\mathrm{}_{\mathrm{int}}/\mathrm{}_m`$. When this ratio is not small, many oscillations happen within one interaction length, thus allowing coherent MSW transitions to take place, particularly at low temperatures where $`\mathrm{}_{\mathrm{int}}`$ is large. Using the diagonalisation matrix, one can also see that the oscillation amplitude is large. See Refs.  for further discussions on extracting the MSW effect from the QKEs.

The sterile neutrino production equation obtained in the fully decoupled limit is simply

$$\frac{z_s}{t}\frac{1}{2}\lambda _3[z_\alpha z_s],$$
(52)
which is in agreement with the static approximation. The related asymmetry evolution equation is also of static approximation form.

In the opposite thermal equilibrium limit, the QKEs become also a homogeneous system of three differential equations if chemical potentials are ignored in the active neutrino distribution functions, so that $`\frac{z_\alpha }{t}=0`$ and thus $`\frac{P_0}{t}=\frac{P_z}{t}`$, with the matrix $`𝒦`$ given by

$$𝒦=\left(\begin{array}{ccc}D& V_z& 0\\ V_z& D& V_x\\ 0& \frac{V_x}{2}& 0\end{array}\right)=\frac{1}{\mathrm{}_0}\left(\begin{array}{ccc}d& (c\pm ab)& 0\\ (c\pm ab)& d& s\\ 0& \frac{s}{2}& 0\end{array}\right).$$
(53)
Repeating the same procedure as before yields only one difference: the denominator of $`\lambda _3`$ now requires the substitution $`V_x^2V_x^2/2`$. This, however, has no practical importance in the collision dominated regime, since the condition $`r|_{\mathrm{res}}(\mathrm{}_{\mathrm{int}}/\mathrm{}_m)|_{\mathrm{res}}=(V_x/D)|_{\mathrm{res}}1`$ is required to validate the approximate equations at resonance. As $`r|_{\mathrm{res}}1`$, collision dominance collapses locally; MSW transitions become the chief amplification mode .

To conclude, we find that the rate of sterile neutrino production and asymmetry generation is the same at leading order in $`V_x^2`$, independent of the description of the active neutrino distributions: full decoupling , thermal equilibrium , or the more general equation . Deviations from the initial thermal equilibrium distribution have no noticeable effects on the asymmetry generation mechanism. This statement is qualitatively understood, because strongly damped oscillations do not produce significant spectral distortions.

### C Comparison of the QKEs and the static approximation

The dash-dotted curves in Figs. 1a and 2a show the integration of the static approximation equations (40) and (43), which are in excellent agreement with the QKE results (solid lines) prior to and during explosive growth. The essential physics in the collision dominated regime is very well captured by the static approximation: the MSW effect is damped, and rate equations for the neutrino distribution functions govern the evolution. In terms of the QKEs’ adiabatic limit, the transition rates in these equations are simply the eigenvalue $`\lambda _3`$.

Moreover, one can easily infer from Eq. (43) the existence of a critical temperature $`T_c`$ and its approximate value by examining the dominant term $`A`$ . This term has a fixed point at $`L=0`$ that is stable above $`T_c`$, and unstable below it where deviations from $`L=0`$ will cause runaway positive feedback. The transition is controlled by the function $`(bc)`$ in Eq. (44). Since $`bT^6`$, the $`A`$ term is generally positive for $`\delta m^2<0`$ at high temperatures, and $`L`$ is destroyed. At $`T_c`$, the “effective $`A`$” \[the $`y`$ integration in Eq. (43) must be performed to determine the net effect of $`A`$\] becomes negative, leading to a spurt of exponential growth.

The subdominant $`B`$ and $`B^{}`$ terms in Eq. (43) control the extent to which the condition $`L=0`$ is approximated prior to $`T_c`$; realistically $`L=0`$ does not exactly define a fixed point since $`B`$ and $`B^{}`$ are not proportional to $`L`$. Numerically, we find that $`L10^{15}`$ is a typical value. Interestingly, DHPS report values as low as $`10^{100}`$. The reason for this is clear: the $`B`$ and $`B^{}`$ terms in Eq. (43) arise from the nonzero chemical potentials in $`z_\alpha `$ and $`\overline{z}_\alpha `$. These terms are absent in the DHPS equations since chemical potentials are neglected in the repopulation term, even when they attempt to numerically solve the QKEs \[see Eq. (9) of DHPS and following discussion\]. In their case $`L=0`$ is much closer to being a true fixed point. This is a relatively minor issue, and not germane to the controversy over the final asymmetry value. Nevertheless, it is of some numerical interest since tracking values of order $`10^{15}`$, rather than $`10^{100}`$, is a much easier task. Also, if $`10^{100}`$ was the correct number, one must consider the presence of statistical fluctuations in the early universe as they would be able to change randomly the sign of the asymmetry at different spatial points . In the presence of chemical potentials in the active neutrino distributions, however, this effect can be excluded  for a large range of mixing parameters (roughly for $`\mathrm{sin}^22\theta _010^6`$).

On the practical front, the neglect of chemical potentials when numerically integrating the QKEs leads to an artificially higher sensitivity to small effects on the solution (numerical error, oscillatory terms). DHPS observed this sensitivity, but did not realise that the inclusion of chemical potentials was remedial from this perspective. This may be one reason why DHPS found that their results were chaotic in sign for some mixing parameters, while an analysis performed with chemical potentials did not . DHPS took this synthetic chaoticity as grounds to not rely on their own numerical solutions to the QKEs, or those of others .§§§Curiously, DHPS relied on the results of Ref. , which also reported rapid oscillations in the sign of the neutrino asymmetry for $`\delta m^2<10^7`$ eV2, although the equations used in Ref.  are no more accurate than those in Refs. . DHPS may not have been aware of this incongruity since Ref.  does not contain any figure of the asymmetry behaviour, although rapid sign changes are clearly reported.

After explosive growth, Figs. 1a and 2a show that the static approximation equations significantly underestimate the subsequent augmentation of $`L`$. This is the epoch in which the adiabatic MSW effect takes over as the principal propellant.

### D MSW dominated regime

Collision dominance and thus the static approximation collapse when $`\mathrm{}_{\mathrm{int}}/\mathrm{}_m`$ approaches and becomes larger than one . This is particularly evident in the adiabatic limit approach . Consider the eigenvalues of the matrix $`𝒦`$ at resonance:

$$\lambda |_{\mathrm{res}}=D,\frac{D}{2}\pm i\frac{\sqrt{4V_x^2D^2}}{2}.$$
(54)
When collision dominance breaks down, i.e., $`(V_x/D)|_{\mathrm{res}}1`$, the oscillatory imaginary parts are not preferentially damped. The real eigenvalue $`\lambda _3`$ and the real parts of the complex conjugate pair are now comparable in size. The imaginary components drive the MSW effect , inducing many large amplitude oscillation cycles within the decoherence time scale.

The simplest way to explore the implications of the MSW effect is through the approach of Ref.. We begin with a very important observation: after the initial burst of $`L`$ creation, the neutrino and antineutrino resonances become separated in momentum space. This is illustrated in Figs. 1b and 2b, and is easy to see analytically from Eqs. (25) and (30). At $`T_c`$, Eq. (25) holds, with neutrinos and antineutrinos having the same resonance momentum. After a substantial $`L`$ has been created, the resonance condition evolves to being given by Eq. (30), and the resonance momenta are separated. If $`L>0`$, then the neutrino resonance quickly moves to the tail of the distribution, while the antineutrino resonance stays in the body. The MSW effect therefore has very asymmetric consequences for neutrinos and antineutrinos: the former are unaffected, while the latter continue to be processed into sterile states, thereby adding to the asymmetry. We shall henceforth assume $`L>0`$ for definiteness.

Let us focus now on the large region of parameter space for which the MSW transitions are adiabatic. Here, complete conversion of antineutrinos at the resonance momentum into sterile states means that the asymmetry is augmented, while the position of the resonance evolves with $`L`$. Indeed, the rate of asymmetry growth depends on how quickly the resonance moves through the distribution. In the case of a small resonance width, it is easy to make this connection heuristically, which in turn gives rise to a simple evolution equation

$$\frac{dL_{\nu _\alpha }}{dT}=\frac{1}{n_\gamma }\frac{T^3}{2\pi ^2}[f_{\overline{\nu }_\alpha }(y_{\mathrm{res}})f_{\overline{\nu }_s}(y_{\mathrm{res}})]y_{\mathrm{res}}^2\frac{dy_{\mathrm{res}}}{dT}X\frac{dy_{\mathrm{res}}}{dT}.$$
(55)
Using the resonance condition of Eq. (30), one obtains from Eq. (55) the expression

$$\frac{dL_{\nu _\alpha }}{dT}=\frac{4}{T}\frac{Xy_{\mathrm{res}}}{1+\frac{Xy_{\mathrm{res}}}{L_{\nu _\alpha }}},$$
(56)
as a non-linear asymmetry evolution equation, for the case where $`\frac{dy_{\mathrm{res}}}{dT}<0`$. Repopulation is incorporated as per the procedure described in Ref..

The result of solving Eq. (56) was first presented in Fig. 1 of Ref. , and compared with the QKE solution for temperatures $`TT_c/2`$. Excellent agreement was obtained, and the final asymmetry was found to be in the $`0.20.37`$ range depending on the oscillation parameters. Numerical solutions of the QKEs and of Eq. (56) yielded completely consistent results for the temperature range where the latter equation was expected to be a good approximation.

Notice in particular that when $`L_{\nu _\alpha }1`$, Eq. (56) simplifies to

$$\frac{dL_{\nu _\alpha }}{dT}\frac{4L_{\nu _\alpha }}{T},$$
(57)
which immediately implies that $`L_{\nu _\alpha }T^4`$ in that regime. The $`T^4`$ behaviour that can be empirically measured from the QKE curves is fully explained by the adiabatic MSW effect.

As the asymmetry continues to grow, the antineutrino resonance evolves out of the body into the high momentum end of the distribution. The asymmetry becomes frozen at some final steady state value, whose approximate magnitude can be easily understood by integrating the Fermi–Dirac distribution from $`y=y_{\mathrm{res}}^{\mathrm{low}}0`$ to $`y=\mathrm{}`$:

$$L_{\nu _\alpha }^{\mathrm{final}}\frac{1}{4\zeta (3)}\left(\frac{T_{\nu _\alpha }}{T_\gamma }\right)^3_0^{\mathrm{}}\frac{y^2dy}{1+e^y}=\frac{3}{8}\left(\frac{T_{\nu _\alpha }}{T_\gamma }\right)^3.$$
(58)
The temperature ratio takes care of reheating due to $`e^+e^{}`$ annihilations at $`Tm_e0.5`$ MeV. Equation (58) gives the approximate magnitude only, because the distribution changes with time as the asymmetry is created. Numerically, the final values found by incorporating proper thermalisation effects are quite close to this estimate.

### E Combining the collision dominated and MSW dominated regimes.

Figures 1a and 2a summarise neatly the overall picture of asymmetry evolution, as discussed in Sec. I. The solid and dash-dotted lines represent solutions to the exact QKEs and the static approximation equations respectively. Clearly, the latter offers an excellent description for the explosive growth phase, but underrates the subsequent growth of $`L_{\nu _\alpha }`$. The dashed line is generated by a code that computes the ratio between the resonance interaction and matter-affected oscillation lengths, $`r|_{\mathrm{res}}(\mathrm{}_{\mathrm{int}}/\mathrm{}_m)|_{\mathrm{res}}`$, at every time step. If $`r|_{\mathrm{res}}<1`$, the static approximation equations are used to advance the system. For $`r|_{\mathrm{res}}>1`$, pure adiabatic MSW evolution is employed via Eq. (56). In both the high and low $`T`$ regimes, the dashed line and the QKEs are in superb agreement. The small disparity directly after explosive growth arises from the fact that neither effect is dominant in this intermediate stage. We conclude that the physics of asymmetry growth has been thoroughly understood.

## IV Comparison with the DHPS results

In Sec. IV A we “deconstruct” the DHPS equations, while we devote Sec. IV B to pinpointing where in their derivation the MSW effect is unknowingly neglected.

### A An analysis of the DHPS equations

For comparison with the static approximation, we rewrite DHPS’s evolution equations in the notation of this Comment using the “conversion relations” given in the Appendix.

Equation (42) of DHPS describes sterile neutrino production. Its translation is

$$\frac{z_s^+}{(T/\mathrm{MeV})}=\frac{s^2\mathrm{\Gamma }}{4T}\frac{d^2+(bc)^2+a^2}{\mathrm{\Delta }(0)}[1z_s^+].$$
(59)
Comparing this with Eq. (40), one sees that it agrees with the latter’s leading term, when (i) chemical potentials are neglected and thus $`z_\alpha ^+1`$, and (ii) the replacement

$$\mathrm{\Delta }(s^2)\mathrm{\Delta }(0),$$
(60)
is made.

The neutrino asymmetry evolution equation \[Eq. (51) of DHPS\] is equivalently

$$\frac{1}{L}\frac{dL}{d(T/\mathrm{MeV})}=_0^{\mathrm{}}\frac{A}{\mathrm{\Delta }(0)}\left[1\frac{s^2}{4}\frac{d^2+(cb)^2+a^2}{\mathrm{\Delta }(0)}\right]𝑑y.$$
(61)
The first term on the right hand side is again connected to the leading $`A`$ term of the static approximation \[Eq. (43)\], with the substitution Eq. (60). The $`B`$ and $`B^{}`$ terms are absent as expected, since DHPS neglect chemical potentials in the $`\nu _\alpha `$ and $`\overline{\nu }_\alpha `$ distribution functions.

The second term in the square brackets is DHPS’s back-reaction term $`B_1`$. This supposedly new term has, in fact, an interesting interpretation. Consider the function $`\mathrm{\Delta }(s^2/8)`$. One can perform the expansion

$`{\displaystyle \frac{A}{\mathrm{\Delta }(s^2/8)}}`$ $`=`$ $`{\displaystyle \frac{A}{\mathrm{\Delta }(0)}}\left[1{\displaystyle \frac{s^2}{4}}{\displaystyle \frac{d^2+(cb)^2+a^2}{\mathrm{\Delta }(0)}}\right]`$ (63)
$`+𝒪({\displaystyle \frac{s^2}{d^2+(cb+a)^2}},{\displaystyle \frac{s^2}{d^2+(cba)^2}})^2.`$
Thus the back-reaction term is but the first order correction in an expansion of $`\mathrm{\Delta }(s^2/8)`$. But we have already discussed in Sec. III B this sort of expansion, in the context of the $`V_x^2`$ term in the denominator of the eigenvalue $`\lambda _3`$, where the full decoupling and the instantaneous repopulation limits yielded $`\mathrm{\Delta }(s^2)`$ and $`\mathrm{\Delta }(s^2/2)`$ respectively. We argued that this distinction was of no practical importance, because the numerical difference between these two cases was only significant when the collision dominance assumption used in their derivation broke down. In a nutshell: (i) The first order terms are only large when the MSW effect is important. (ii) When the MSW effect is important, equations such as obtained from the static approximation, and DHPS (51), are no longer a good approximation to the QKEs.

Let us now evaluate Eq. (63) for the antineutrino ($`cba=0`$) resonance:

$$\frac{A}{\mathrm{\Delta }(s^2/8)}|_{\mathrm{res}}\frac{A}{\mathrm{\Delta }(0)}|_{\mathrm{res}}\left[1\frac{s^2}{4}\frac{d_{\mathrm{res}}^2+2a^2}{d_{\mathrm{res}}^2(d_{\mathrm{res}}^2+4a^2)}\right].$$
(64)
When the asymmetry grows at the critical temperature $`T_c`$, the resonance condition passes from a regime in which $`bc`$ when $`L`$ is negligible, to $`ac`$. In this second regime the term $`d_{\mathrm{res}}^2a^2`$ and the last expression becomes simply,

$$\frac{A}{\mathrm{\Delta }(s^2/8)}|_{\mathrm{res}}\frac{A}{\mathrm{\Delta }(0)}|_{\mathrm{res}}\left[1\frac{s^2}{8d_{\mathrm{res}}^2}\right].$$
(65)
So we see that the expansion parameter on resonance is $`s^2/d_{\mathrm{res}}^2`$, i.e., the square of the ratio between the local interaction and matter-affected oscillation lengths $`(r|_{\mathrm{res}})^2`$. As explained earlier, a small $`r|_{\mathrm{res}}`$ denotes that collisions are frequent enough to interrupt coherent evolution and thus to prevent the MSW effect, while a large $`r|_{\mathrm{res}}`$ means MSW conversion can occur as the collisional mechanism switches off. Of course this transition is gradual.

Moreover, using the equality $`d=\delta b`$, one can also see, using Eqs. (27) and (28), that

$$\frac{s^2}{d_{\mathrm{res}}^2}=\left(\frac{s}{\delta }\right)^2\left(\frac{y_{\mathrm{res}}^0}{y_{\mathrm{res}}}\right)^4=\left(\frac{s}{\delta }\right)^2\left(\frac{|L|}{L_t}\right)^3,$$
(66)
where $`y_{\mathrm{res}}^0`$ is the resonant momentum when the asymmetry is negligible. Recall that $`y_{\mathrm{res}}^0`$ is very small at high $`T`$, but grows as the temperature decreases, reaching $`2`$ at $`T_c`$. The quantity $`y_{\mathrm{res}}`$ is the general resonance momentum, which attains a minimum value depending on $`\delta m^2`$ (smaller for low $`\delta m^2`$) when the asymmetry grows (see Figs. 1b and 2b). It is thus clear from Eq. (66) that if one continues to use $`B_1`$ in the evolution equation when the asymmetry is large and $`s^2/d_{\mathrm{res}}^2𝒪(1)`$, the rise of the asymmetry will be artificially made slower. From that moment its evolution will be approximately $`|L|50L_t`$ and from the equations Eq.(29) and Eq.(28) this implies that $`|L|1/T`$, a behaviour observed by DHPS in their results. In this way the solution can grow just one order of magnitude since the onset of the instability.

In summary, the back-reaction term $`B_1`$ is generated by an expansion in $`s/d_{\mathrm{res}}`$, and not simply $`s/\delta `$ as DHPS believe. Their evolution equation is therefore invalid for temperatures at which $`s^2/d_{\mathrm{res}}^2>1`$. As emphasised before , both the static approximation and the DHPS equations belong to a family of collision dominated idealisations which require $`s/d1`$ (plus other conditions) for their validity. When $`s/d1`$, the MSW effect is missed in both cases. Clearly DHPS did not properly appreciate the constraints on the real expansion parameter adopted in their derivation, and hence trusted their equations’ validity for all temperatures.

For very large $`\delta m^2`$, DHPS do actually observe a large final value. This is because $`|L|_{\mathrm{fin}}10^5(|\delta m^2|/\mathrm{eV}^2)^{1/3}(T_c/T_f)10^4(|\delta m^2|/\mathrm{eV}^2)^{1/3}`$ where $`T_f`$ is the freezing temperature corresponding to a resonant momentum (of antineutrinos if $`L`$ is positive, of neutrinos if $`L`$ is negative) well in the tail of the distribution ($`y_{\mathrm{res}}^f10`$) and is very weakly dependent on the exact choice of $`y_{\mathrm{res}}^f`$. Thus for $`|\delta m^2|10^9\mathrm{eV}^2`$ they find $`|L|_{\mathrm{fin}}0.1`$, and are therefore still able to state that active–sterile neutrino oscillations can yield large asymmetries. However, this high $`\delta m^2`$ asymmetry rise is spurious. For $`\delta m^210^4\mathrm{eV}^2`$, the effective potential as given by Eq. (14) is meaningless, since the critical temperature is now much larger than $`100`$ MeV. The conditions of the plasma are different, especially at temperatures above the quark–hadron phase transition. The dynamics of neutrino oscillations in this regime is actually an open problem.

### B How DHPS neglect the MSW effect

In extracting their simplified evolution equations from the QKEs, DHPS make use of two small expansion parameters: (i) $`1/Q`$, and (ii) $`\mathrm{sin}2\theta `$. We now examine their merits.

1. Small $`1/Q`$ expansion. The derivation begins with an attempt to solve analytically two out of a system of eight coupled differential equations, namely, Eqs. (27) and (28) in DHPS:

$`h_{}^{}/Q`$ $`=`$ $`Ul_{}\gamma h_{}VZl_+,`$ (67)
$`l_{}^{}/Q`$ $`=`$ $`{\displaystyle \frac{F}{2}}(a_{}s_{})Uh_{}\gamma l_{}+VZh_+,`$ (68)
where prime denotes differentiation with respect to $`\tau `$. These expressions are equivalent to

$`{\displaystyle \frac{2p}{\delta m^2}}{\displaystyle \frac{}{t}}(P_x\overline{P}_x)`$ $`=`$ $`(1b)(P_y\overline{P}_y)d(P_x\overline{P}_x)+a(P_y+\overline{P}_y),`$ (69)
$`{\displaystyle \frac{2p}{\delta m^2}}{\displaystyle \frac{}{t}}(P_y\overline{P}_y)`$ $`=`$ $`s(P_z\overline{P}_z)(1b)(P_x\overline{P}_x)d(P_y\overline{P}_y)a(P_x+\overline{P}_x),`$ (70)
in our language, with $`c1`$. Solutions to Eq. (67) are generically oscillatory. To this end, DHPS devise an oscillation averaging procedure that consists of setting all terms containing $`1/Q`$ (i.e., the $`\tau `$-derivatives) to zero, from which they obtain what are actually “steady state” solutions to the variables $`h_{}`$ and $`l_{}`$. We now explain why this method is flawed.

Firstly, the variable $`P_y`$ enters the big picture as the “driving force” for $`\nu _s`$ production $`\frac{z_s}{t}`$ \[see Eqs. (34) and (39)\]. In the context of solving differential equations, one may ignore contributions from oscillatory components in the driving force $`P_y`$ only if the oscillation amplitude is much smaller than the size of its steady “main term”. Otherwise, oscillations propel the evolution, producing a marked change in the integrated variable $`z_s`$ (i.e., the MSW effect) when the oscillation frequency is a minimum at resonance, in the same way that $`\frac{d}{dt}\left(\frac{1}{\omega }\mathrm{sin}\omega t\right)=\mathrm{cos}\omega t`$.

In their derivation, DHPS do not consider at all the magnitudes of these oscillatory terms, and discard them on the basis that they impair computational efficiency, when in fact even the adiabatic solution to $`P_y`$ (and $`\overline{P}_y`$) is purely oscillatory at resonance in the post explosive growth environment where a relatively large $`r|_{\mathrm{res}}(\mathrm{}_{\mathrm{int}}/\mathrm{}_m)|_{\mathrm{res}}`$ ratio prevails . Their equations cannot describe the sharp change in $`z_s`$ at the low temperature resonance because the dominant and strongly frequency-dependent term responsible for the MSW effect has been removed. Indeed, one cannot know the sizes of the steady and oscillatory components in $`P_y`$ or $`l_{}`$ without some input of the initial conditions; to eliminate dynamical variables in the equation’s differential form is premature and may be dangerous.We also mention a point noted in Ref. , that when the asymmetry has grown sufficiently to separate the neutrino and antineutrino resonances, the differences in the distributions, and their derivatives, can be as large as their sums, simply because one distribution is on resonance while the other is far from it.

Moreover, the process of artificially factoring out a large quantity by redefining the integration variable and then exploiting its inverse as an expansion parameter is not in fact a well-justified manœuvre to obtaining steady state solutions. To see this, one needs only to consider the case of exponential growth/decay, $`\frac{dx}{dt}=Ax`$, where the approximation $`\frac{1}{A}\frac{dx}{dt}0`$ and thus the conclusion $`x0`$ crumble for a positive $`A`$.

Fortunately, steady state solutions to $`h_{}`$ and $`l_{}`$ are indeed valid for some oscillation parameters $`(\delta m^2,\mathrm{sin}2\theta )`$ when collisions are frequent enough to quickly damp the oscillations, even at resonance, provided that the adiabatic condition is also satisfied . Thus DHPS have managed to reproduce results (albeit from a somewhat dubious procedure) that are consistent with the static approximation, in the manner that the latter may also be derived in the $`ds`$ limit by setting $`\frac{dP_x}{dt}\frac{dP_y}{dt}0`$. But, like the static approximation, their equations do not apply to the low temperature region after explosive growth where steady state solutions do not hold at resonance.

2. Small $`\mathrm{sin}2\theta `$ expansion. This approximation is used on two occasions, first as a pretext for disregarding the variable $`P_z\overline{P}_z`$ (or $`a_{}s_{}`$) in Eq. (67) \[Eq. (28) in DHPS\], and then in evaluating the eigenvalues and other properties of the matrix $``$ in their Eq. (33).

The first exploit is precarious from the perspective of solving differential equations, since $`P_{x,y,z,0}`$ (and therefore $`a_\pm ,s_\pm ,h_\pm ,l_\pm `$) are interdependent dynamical variables and one has no prior knowledge on their sizes. Indeed, it is far more seemly to assume a priori the $`P_z\overline{P}_z`$ term to be of some substance at resonance when the MSW effect kicks in at low temperatures, given that only antineutrinos (for $`L>0`$) can be converted into sterile states in this epoch by the said mechanism. The neglect of this term is also related to DHPS’s observation of the asymmetry reaching arbitrarily small values prior to explosive growth, to be discussed later in Sec. V.

In the second instance, the expansion procedure employed to diagonalise the matrix

$$=\left(\begin{array}{cccc}0& 0& 0& F\\ 0& 2\gamma & 0& F\\ 0& 0& \stackrel{~}{\gamma }& \stackrel{~}{U}\\ \frac{F}{2}& \frac{F}{2}& \stackrel{~}{U}& \stackrel{~}{\gamma }\end{array}\right),\stackrel{~}{U}=U\left(1\frac{D^2}{\gamma ^2+U^2}\right),\stackrel{~}{\gamma }=\gamma \left(1+\frac{D^2}{\gamma ^2+U^2}\right),$$
(71)
where $`Fs`$, $`\gamma =\delta /\tau ^2d`$, $`DaL`$ and $`U=1/\tau ^211b`$, is not itself at flaw. However, it does have its own limitations which DHPS did not notice.

Equation (71) suggests that the term $`F`$ (or $`s`$) should be compared with $`\gamma `$, $`U`$, $`\stackrel{~}{\gamma }`$ and $`\stackrel{~}{U}`$ if it is to be treated as a perturbation. DHPS consider only $`s`$ and $`\delta `$, presumably because $`\gamma =\delta `$ at resonance, which they define as occurring at $`\tau =1`$ such that $`U=0`$ (or $`c=b`$). However, this resonance condition assumes a small asymmetry, i.e., $`|b||a|`$. Thus it is able to describe the high temperature regime where the term $`b`$ dominates, but not after a sizeable asymmetry has been created, where the now separated neutrino and antineutrino resonances occur at any two $`\tau `$’s other than $`\tau =1`$. The ratio $`s/\delta `$ is meaningless in this second regime.

Lastly, it must be stressed that MSW flavour conversion is not a perturbative effect. Small $`\mathrm{sin}2\theta `$ expansion equals no MSW effect, in the spirit that treating the off-diagonal $`\mathrm{sin}2\theta `$ in the Schrödinger Equation for a two-flavour collisionless system as a perturbation will not yield the correct behaviour at resonance.

## V Other comments

### A The mystery of the oscillation in sign

DHPS allege that results presented by other independent groups are in conflict on the issue of oscillations in the sign of the asymmetry at the critical temperature. They use this to cast doubt on the credibility of some of these works. The true status of this subject is actually very clear, as we will explain below.

To begin with, one must realise that two different systems are discussed in the literature. The first is a simplified scenario in which all neutrinos have the same momentum, usually taken to be the zero chemical potential thermal average, $`p=p3.15T`$. The second is the realistic case where the neutrinos are described by a Fermi–Dirac distribution. The mean momentum simplification should be viewed as a “toy model” which is not always a good approximation to the realistic case. However, it does give some useful qualitative insight.

Another critical issue is whether the exact QKEs or an idealisation such as the static approximation are used in the calculations.

Using the mean momentum approximation, Ref. observed oscillations of sign for some parameter points and apparently drew the incorrect conclusion that oscillations were a generic feature of lepton number generation. The simplified system was revisited in Ref., where oscillations were found to occur for only part of the parameter space (see Fig. 2 of Ref. ). In both papers, the mean momentum QKEs were solved Recently a detailed analytic study of the mean momentum model was presented in Ref...

In the realistic Fermi–Dirac case, the issue of the sign of $`L_{\nu _\alpha }`$ was first discussed in Ref. in the context of the static approximation; in the region of parameter space where this approximation is valid, the sign does not oscillate. However, as also pointed out in that paper, the static approximation breaks down at large $`\mathrm{sin}^22\theta `$ due to the rapid rise in lepton number during the explosive growth phase. Similar results were obtained in Ref..

The oscillation issue was addressed in detail in Ref. , where comparisons of the static approximation and the QKEs showed that they gave concordant results for $`\mathrm{sin}^22\theta 10^6`$. For larger values, sign oscillations were found using the exact QKEs, but not the static approximation simply because the latter is invalid for too large a $`\mathrm{sin}^22\theta `$. Also, Refs.  and  are not in mutual disagreement, since the former examined the realistic case with a Fermi–Dirac distribution, while the latter considered the mean momentum approximation.

### B Abrupt cut-off at large mixing angles?

In an attempt to locate the region of parameter space for which neutrino asymmetry generation is possible, DHPS found the following approximate upper and lower bounds on the so-called “region of instability”:

$`\mathrm{sin}^22\theta _0\sqrt{|\delta m^2|/\mathrm{eV}^2}`$ $``$ $`10^6,`$ (72)
$`\mathrm{sin}^22\theta _0(|\delta m^2|/\mathrm{eV}^2)^{\frac{1}{6}}`$ $``$ $`4\times 10^{14}.`$ (73)
While the second constraint agrees with that found in Ref. , the first is too severe compared with previous calculations .\**\**\**The lower bound in Ref.  is superficially more stringent than Eq. (73) ($`5\times 10^{10}`$ instead of $`4\times 10^{14}`$) simply because a more restrictive condition — that the asymmetry must grow to the $`10^5`$ level — was imposed in Ref. , while DHPS demand only amplification by as little as one order of magnitude. However, both generate the same slope on a $`(\mathrm{log}\delta m^2,\mathrm{log}\mathrm{sin}2\theta `$) graph. In addition, DHPS observe an abrupt transition from the “unstable” to the “stable” (i.e., no growth) regions. This is also at odds with earlier works, and is most likely a byproduct of numerical errors for the following reason.

Increasing $`\mathrm{sin}2\theta `$ for a fixed $`\delta m^2`$ intensifies sterile neutrino production, which, in turn, has the effect of shifting the critical momentum to values larger than $`2`$ . By Eq. (26), the critical temperature decreases correspondingly, and explosive asymmetry generation is thereby delayed. (The growth curve also becomes smoother but a substantial asymmetry is still eventually attained.) It follows that for very large mixing angles, no explosive growth can occur since the critical momentum is then banished to the tail of the distribution.

DHPS note correctly the connection between suppression of asymmetry amplification for a large $`\mathrm{sin}2\theta `$ and a substantial sterile neutrino distribution manufactured prior to the onset of an otherwise explosive growth. The latter may be approximated by a theta function in momentum space, such as

$$z_s^+(y)=\left\{1\mathrm{exp}\left[4.3(2.25)\times 10^4s^2\left(\frac{\delta m^2}{\mathrm{eV}^2}\right)^{\frac{1}{2}}\right]\right\}\theta (y_{\mathrm{res}}^0y),$$
(74)
derived from Eq. (36) for $`\alpha =e(\mu ,\tau )`$. Substituting Eq.(72) into the above, one finds that DHPS’s upper bound corresponds to having a final sterile neutrino distribution of $`z_s^+=0.04(0.02)`$. DHPS found an expression similar to Eq. (74) in the limit $`s1`$ for $`\alpha =e`$ \[see their Eq. (54)\]. They argue that such accumulation of sterile neutrinos is able to destroy completely the asymmetry generation since the factor $`1z_s^+`$ appearing in the term $`A`$ in the growth rate \[Eq. 43\] now favors stability.

As explained before, in the regime prior to the onset of asymmetry generation, the static approximation and DHPS evolution equations are in complete agreement at leading order; this suppression effect is to be distinguished from DHPS’s supposedly new but as yet negligible “back-reaction”. On the other hand, any interested person can numerically integrate the static approximation equations at leading order and find no abrupt cut-off even when the sterile neutrino distribution increases to the level of $`z_s^+0.04`$. Therefore the severity of DHPS’s upper bound has to be imputed to numerical errors.††††††The calculation of $`z_s^+`$ should be done numerically, coupled with the asymmetry evolution equation. The analytical expression in Eq. (74) approximates the production with a theta function which is actually a sharp jump around the resonance, valid only before the onset of neutrino asymmetry generation (see Ref.  for more details).

Physically, the negative contribution from a non-negligible sterile neutrino distribution associated with a large $`\mathrm{sin}2\theta `$ is compensated for by an increase in the critical momentum. In more descriptive words, the asymmetry generation “waits” for the onset of instability when this compensation occurs at lower temperatures (see Ref.  for further discussions).

### C Chemical potentials and final value of neutrino asymmetry

As noted in the Sec. III C, the role of chemical potentials prior to the onset of explosive growth is to prevent the asymmetry from being diminished to arbitrarily small values. This function is described by the term $`B`$ in the static approximation \[Eq. (43)\], or, equivalently, contained in the variable $`s(z_\alpha \overline{z}_\alpha )`$ (or $`Fa_{}`$ in DHPS’s notation) in the QKEs \[Eq. (67)\], which DHPS ignored. Another mechanism that can inhibit drastic lepton number destruction is the presence of a sterile neutrino asymmetry, created by oscillations from an initial $`\alpha `$-neutrino asymmetry. This effect is manifested in a nonzero $`s(z_s\overline{z}_s)`$ (or $`Fs_{}`$) in the QKEs, and corresponds to the term $`B^{}`$ in the static approximation. In the case of a zero initial asymmetry, this term is negligible. On the other hand, the term $`B`$ due to finite chemical potentials must never be neglected. The fact that DHPS omit both $`B`$ and $`B^{}`$ causes their solution to collapse to values of order $`10^{100}`$ prior to the onset of asymmetry growth.

Furthermore, observe that the condition of lepton number conservation \[Eq. (33)\] and the expressions for $`P_0`$ and $`\overline{P}_0`$ in Eq. (8) together lead to the following corollary:

$$_0^{\mathrm{}}\mathrm{\Gamma }(y)\left[\left(z_\alpha ^{\mathrm{eq}}(\mu )\overline{z}_\alpha ^{\mathrm{eq}}(\overline{\mu })\right)\left(z_\alpha \overline{z}_\alpha \right)\right]f_{\mathrm{eq}}^o(y)y^2𝑑y0.$$
(75)
Clearly, if one neglects the chemical potential $`\mu _\alpha `$, the remaining integral is no longer vanishing, but becomes a fictitious friction force proportional to $`\mathrm{\Gamma }L_{\nu _\alpha }`$ in Eq. (34) through the otherwise absent $`\frac{P_0}{t}\frac{\overline{P}_0}{t}`$ term. This point was noted in Ref.  and speculated to be culpable for DHPS’s distorted results. However (and we are arguing in DHPS’s favour), this phantom force does not in fact appear in DHPS’s final evolution equation (51), as they specifically define $`_0^{\mathrm{}}\mathrm{\Gamma }(y)\left(z_\alpha \overline{z}_\alpha \right)f_{\mathrm{eq}}^o(y)y^2𝑑y0`$ to be their lepton charge conservation condition \[their Eq. (30)\] which, in the absence of chemical potentials, is actually mathematically consistent with Eq. (33). In effect, DHPS have taken, though fortuitously, the appropriate measures to negate the artificial friction force inherent in their formulation.

## VI Conclusions

We review how the numerical solution to the QKEs supplies a lepton asymmetry growth curve whose features can be physically understood and reproduced by approximate evolution equations. In the higher temperature collision dominated epoch, the static approximation mimics the explosive amplification extremely accurately, and allows one to understand why a critical temperature exists for $`\delta m^2<0`$ and $`\mathrm{cos}2\theta _01`$. As the collision rate decreases with temperature, the MSW effect eventually takes over as the dominant amplifier, since the transitions are adiabatic for a large range of mixing parameters. Also, the separation of the neutrino and antineutrino resonance momenta after a significant $`L`$ has been created allows efficient conversion of antineutrinos (for $`L>0`$) into sterile states, while oscillations are matter-suppressed for neutrinos (and vice versa for $`L<0`$). The roles of the static approximation and the MSW effect are succinctly depicted in Figs. 1a and 2a.

The major flaw in the DHPS analysis is their neglect of the pure MSW effect in the post explosive growth epoch. Their Eqs. (42) and (51) are approximations to the QKEs valid in the collision dominated regime and for a sufficiently small $`L`$. Like the static approximation, they cannot describe the MSW effect, and therefore underrate the asymmetry growth at lower temperatures. However, the DHPS underestimation is much more severe than the static case, because their Eq. (51) contains a bogus back-reaction term $`B_1`$ that becomes large and artificially terminates the asymmetry growth under precisely the same conditions that imply MSW dominance. This can be understood from the symbolic relation between the two approximation schemes: $`\frac{1}{(1+x)}1x`$, where the left hand side represents the static situation and the right hand side correlates with DHPS. The back reaction term $`x`$ leads to a severe, but fictitious, cut-off when the MSW effect dominates at $`x1`$.

We conclude that origin of the $`0.20.37`$ range for the final steady state value of the asymmetry is well understood analytically and physically, in a way which is completely consistent with brute force numerical solutions of the Quantum Kinetic Equations.

###### Acknowledgements.

This work was supported by the Australian Research Council, the Commonwealth of Australia and Istituto Nazionale di Fisica Nucleare (INFN). We thank R. Buras, S. Hansen, K. Kainulainen and D. Semikoz for interesting correspondence.

## A Translation of DHPS notation

We give conversion relations between the DHPS nomenclature and the notation in this Comment, denoted by the sign “$``$”, with DHPS on the left and our equivalent on the right.

The effective potential is

$`V_{eff}^a`$ $`=`$ $`\pm C_1\eta G_FT^3+C_2^a{\displaystyle \frac{G_F^2T^4E}{\alpha }}V_\alpha ,`$ (A1)
$`C_1`$ $`=`$ $`{\displaystyle \frac{2\sqrt{2}\zeta (3)}{\pi ^2}}0.345,C_2^\alpha ={\displaystyle \frac{2\sqrt{2}\zeta (3)}{\pi ^2}}{\displaystyle \frac{A_\alpha }{G_FM_W^2}}0.5846(0.1623),`$ (A2)
for $`\alpha =e(\mu ,\tau )`$, and $`E`$ is the neutrino energy. The quantity

$$\eta L^{(\alpha )}$$
(A3)
is the so-called effective asymmetry (abbreviated to $`L`$ in the body of this Comment).

The neutrino matrix and its expansion in terms of the Pauli matrices are

$$\rho (p)=\left(\begin{array}{cc}\rho _{aa}& \rho _{as}\\ \rho _{sa}& \rho _{ss}\end{array}\right)\frac{f_{\mathrm{eq}}^0}{2}\left(\begin{array}{cc}P_0+P_z& P_xiP_y\\ P_x+iP_y& P_0P_z\end{array}\right).$$
(A4)
For the four matrix elements, DHPS write down the following kinetic equations:

$`i\dot{\rho }_{aa}`$ $`=`$ $`F_0(\rho _{sa}\rho _{as})/2i\mathrm{\Gamma }_0(\rho _{aa}f_{eq}^0),`$ (A5)
$`i\dot{\rho }_{ss}`$ $`=`$ $`F_0(\rho _{sa}\rho _{as})/2,`$ (A6)
$`i\dot{\rho }_{as}`$ $`=`$ $`W_0\rho _{as}+F_0(\rho _{ss}\rho _{aa})/2i\mathrm{\Gamma }_1\rho _{as},`$ (A7)
$`i\dot{\rho }_{sa}`$ $`=`$ $`W_0\rho _{sa}F_0(\rho _{ss}\rho _{aa})/2i\mathrm{\Gamma }_1\rho _{sa},`$ (A8)
with the quantities $`W_0`$ and $`F_0`$ given by the expressions

$$F_0=\frac{\delta m^2}{2p}\mathrm{sin}2\theta _0V_x,W_0=\frac{\delta m^2}{2p}\mathrm{cos}2\theta _0+V_{eff}^aV_z,$$
(A9)
and we have used $`Ep`$. The rates $`\mathrm{\Gamma }_0`$ and $`\mathrm{\Gamma }_1`$ for $`\alpha =e(\mu ,\tau )`$, and their translations are

$`\mathrm{\Gamma }_0`$ $`=`$ $`2\mathrm{\Gamma }_1=g_a{\displaystyle \frac{180\zeta (3)}{7\pi ^4}}G_F^2T^5{\displaystyle \frac{y}{3.15}}=1.13(0.794)G_F^2T^5y,`$ (A10)
$`\mathrm{\Gamma }_0`$ $``$ $`\mathrm{\Gamma },\mathrm{\Gamma }_1D.`$ (A11)
However, the coefficient $`g_a/3.15`$ is smaller than $`k_\alpha `$ \[see Eq. (37)\] found in Ref. . This discrepancy simply leads to a different correspondence of the final asymmetry with a given choice of parameters $`(\mathrm{sin}^22\theta _0,\delta m^2)`$, and is thus irrelevant for the solution of the controversy.

Are the two sets of kinetic equations (8) and (A5) equivalent? Rewriting the former in terms of $`P_{x,y,z,0}`$ by way of Eq. (A4), one finds that they are indeed concordant, save for the replacement $`V_zV_z`$ which is just a matter of convention. Thus we conclude that approximation schemes in DHPS and in Refs.  begin with identical QKEs. However, as discussed in the main text, DHPS use the Fermi–Dirac distribution with zero chemical potential as the active neutrino equilibrium distribution function.

Moving on, DHPS replace the variables $`\rho _{ij}`$ with ones closely connected to $`P_0`$ and $`𝐏`$:

$`\rho _{\alpha \alpha }f_{\mathrm{eq}}^0[1+a]`$ $``$ $`a{\displaystyle \frac{1}{2}}(P_0+P_z)1=z_\alpha 1,`$ (A12)
$`\rho _{ss}f_{\mathrm{eq}}^0[1+s]`$ $``$ $`s{\displaystyle \frac{1}{2}}(P_0P_z)1=z_s1,`$ (A13)
$`\rho _{\alpha s}f_{\mathrm{eq}}^0[h+il]`$ $``$ $`h{\displaystyle \frac{P_x}{2}},l{\displaystyle \frac{P_y}{2}}.`$ (A14)
They then take sums and differences of corresponding functions for neutrinos and antineutrinos, i.e., for any generic variable $`X`$ and its antineutrino counterpart $`\overline{X}`$, they define

$$X^\pm \frac{X\pm \overline{X}}{2}.$$
(A15)
The time integration variable is changed to the inverse of the dimensionless temperature

$$x\frac{\mathrm{MeV}}{T},$$
(A16)
which generates a factor $`x`$. Therefore DHPS introduce the following new quantities:

$`\gamma `$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }_1}{x}}{\displaystyle \frac{D}{x}},`$ (A17)
$`F`$ $`=`$ $`{\displaystyle \frac{F_0}{x}}={\displaystyle \frac{\delta m^2}{2px}}\mathrm{sin}2\theta _0{\displaystyle \frac{V_x}{x}},`$ (A18)
$`W`$ $`=`$ $`{\displaystyle \frac{W_0}{x}}={\displaystyle \frac{\delta m^2}{2px}}\mathrm{cos}2\theta _0+{\displaystyle \frac{V_{eff}^a}{x}}{\displaystyle \frac{V_z}{x}}.`$ (A19)
They also find it useful to separate $`W`$ in two different contributions, $`W=U\pm VZ`$, with

$`U`$ $`=`$ $`{\displaystyle \frac{\delta m^2}{2px}}\left(\mathrm{cos}2\theta _0+C_2^a{\displaystyle \frac{G_F^2T^4E}{\alpha }}\right){\displaystyle \frac{\delta m^2}{2px}}(cb),`$ (A20)
$`VZ`$ $`=`$ $`{\displaystyle \frac{C_1\eta G_FT^3}{x}}{\displaystyle \frac{\delta m^2}{2px}}a,`$ (A21)
where

$$V=\frac{29.6}{x^2},Z=10^{10}\left(\eta _0_0^{\mathrm{}}\frac{dy}{4\pi ^2}y^2f_{\mathrm{eq}}^0a_{}\right).$$
(A22)
The quantity $`Z`$, as defined above, is clearly connected to the lepton asymmetry $`\eta `$, although their exact relationship is not specified by DHPS. This we provide here for completeness:

$$\eta =\frac{4\pi ^2}{\zeta (3)}Z.$$
(A23)
With these definitions, DHPS rewrite the QKEs in the new variables $`s_\pm `$, $`a_\pm `$, $`h_\pm `$, $`l_\pm `$:

$`s_\pm ^{}`$ $`=`$ $`Fl_\pm ,`$ (A24)
$`a_\pm ^{}`$ $`=`$ $`Fl_\pm 2\gamma _+a_\pm 2\gamma _{}a_{},`$ (A25)
$`h_\pm ^{}`$ $`=`$ $`Ul_\pm VZl_{}\gamma _+h_\pm \gamma _{}h_{},`$ (A26)
$`l_\pm ^{}`$ $`=`$ $`{\displaystyle \frac{F}{2}}(a_\pm s_\pm )Uh_\pm +VZh_{}\gamma _+l_\pm \gamma _{}l_{},`$ (A27)
in which terms containing $`\gamma _{}=(\mathrm{\Gamma }_1\overline{\mathrm{\Gamma }}_1)/2`$ are later discarded. This is also the case in Eq. (8). Thus we stress again that the DHPS kinetic equations are equivalent to Eq. (8).

Now DHPS begin to extract their evolution equations by seeking approximate solutions to Eq. (A24). More new notations and symbols are introduced during the derivation, in particular, the integration variable $`x`$ is replaced with $`\tau `$ by way of the function $`U`$:

$$U=\mathrm{cos}2\theta _0\frac{\delta m^2}{2px}\left(\frac{1}{\tau ^2}1\right),$$
(A28)
such that $`U=0`$ at $`\tau =1`$. Then it follows that

$$\tau =\xi \frac{x^3}{y}\sqrt{\frac{c}{b}}=\left(\frac{T_{\mathrm{res}}^0}{T}\right)^3,$$
(A29)
where $`\xi 6.5(12.8)\times 10^3\sqrt{\mathrm{cos}2\theta _0\frac{|\delta m^2|}{\mathrm{eV}^2}}`$, for $`\alpha =e(\mu ,\tau )`$.

At this point DHPS divide all the equations by the quantity $`M`$, where

$$M\frac{\delta m^2c}{2px}1.12\times 10^9\frac{x^2}{y}\frac{|\delta m^2|}{\mathrm{eV}^2},$$
(A30)
and make the replacements $`\frac{1}{M}(\gamma ,F,U,VZ)(\gamma ,F,U,VZ)`$, leading to

$`\gamma ={\displaystyle \frac{2p}{c\delta m^2}}\mathrm{\Gamma }_1={\displaystyle \frac{\delta }{\tau ^2}}cdd,`$ (A31)
$`F=\mathrm{tan}2\theta s,`$ (A32)
$`U={\displaystyle \frac{1}{\tau ^2}}1{\displaystyle \frac{b}{c}}1=\left({\displaystyle \frac{T}{T_{\mathrm{res}}^0}}\right)^61,`$ (A33)
$`VZ{\displaystyle \frac{a}{c}}a.`$ (A34)
These new definitions combine to produce the change

$$\frac{d}{dx}\frac{1}{Q}\frac{d}{d\tau },$$
(A35)
in the left hand side of the kinetic equations, where

$$Q=M\left(\frac{d\tau }{dx}\right)^15.7\times 10^4\sqrt{\frac{|\delta m^2|}{\mathrm{eV}^2}\mathrm{cos}2\theta _0}.$$
(A36)
More new definitions follow from here, whose translations are listed below (with $`c1`$):

$`D`$ $`=`$ $`ZVa,`$ (A37)
$`\sigma ^2`$ $`=`$ $`\gamma ^2+U^2d^2+(cb)^2,`$ (A38)
$`\stackrel{~}{U}`$ $`=`$ $`U\left(1{\displaystyle \frac{D^2}{\sigma ^2}}\right)(b1)\left[1{\displaystyle \frac{a^2}{d^2+(1b)^2}}\right],`$ (A39)
$`\stackrel{~}{\gamma }`$ $`=`$ $`\gamma \left(1+{\displaystyle \frac{D^2}{\sigma ^2}}\right)d\left[1+{\displaystyle \frac{a^2}{d^2+(1b)^2}}\right],`$ (A40)
$`\stackrel{~}{\sigma }`$ $`=`$ $`\stackrel{~}{\gamma }^2+\stackrel{~}{U}^2,`$ (A41)
$`t`$ $`=`$ $`{\displaystyle \frac{1}{\tau }}\left({\displaystyle \frac{T}{T_{\mathrm{res}}^0}}\right)^3=\sqrt{{\displaystyle \frac{b}{c}}}\sqrt{b}.`$ (A42)
The quantity $`\tau `$ is inversely proportional to $`y`$. DHPS thus define the new variable

$$q=y\tau \left(\frac{T_{\mathrm{res}}^0|_{y=1}}{T}\right)^3.$$
(A43)
Moreover they introduce a new function $`b_0`$ that expresses the difference between sum distribution of active neutrinos and sum distribution of sterile neutrinos:

$$b_0s^+1z_s^+,$$
(A44)
such that their final sterile neutrino production equation \[Eq. (42) in DHPS\] is

$$\frac{db_0}{d\tau }=\frac{QF^2\stackrel{~}{\gamma }}{2\stackrel{~}{\sigma }^2}b_0,$$
(A45)
while the asymmetry evolution equation \[their Eq. (51)\] is given by

$$\frac{1}{Z}\frac{dZ}{dq}=\delta Bq^{5/3}_0^{\mathrm{}}𝑑t\frac{t^4(t^21)f_{eq}(tq)b_0(1/t)}{\sigma ^2\stackrel{~}{\sigma }^2}\left[1\frac{F^2(\sigma ^2+D^2)}{4\sigma ^2\stackrel{~}{\sigma }^2}\right].$$
(A46)
With some trivial algebra one can show that

$`\sigma ^2\stackrel{~}{\sigma }^2`$ $`=`$ $`[\gamma ^2+(U+D)^2][\gamma ^2+(UD)^2],`$ (A47)
$`\sigma ^2+D^2`$ $`=`$ $`\gamma ^2+U^2+D^2.`$ (A48)
Using the conversion relations, the DHPS sterile neutrino production equation translates to

$$\frac{z_s^+}{(T/\mathrm{MeV})}=\frac{s^2\mathrm{\Gamma }}{4T}\frac{d^2+(bc)^2+a^2}{\mathrm{\Delta }(0)}[1z_s^+],$$
(A49)
while their neutrino asymmetry evolution equation is equivalent to

$$\frac{1}{L}\frac{dL}{d(T/\mathrm{MeV})}=_0^{\mathrm{}}\frac{A}{\mathrm{\Delta }(0)}\left[1\frac{s^2}{4}\frac{d^2+(cb)^2+a^2}{\mathrm{\Delta }(0)}\right]𝑑y,$$
(A50)
in our language.

Note Added

Shortly after completion of this paper, R. Buras and D. V. Semikoz informed us of their work (hep-ph/0008263), where they also concluded that the main result of DHPS was incorrect.