Derivation of the sterile neutrino Boltzmann equation from quantum kinetics

An extensive, growing body of work has been penned on cosmologies that include one or more sterile neutrinos. Early entries in the literature formulated a Boltzmann-like equation describing sterile-neutrino production in a way that bypasses the numerical tracking of high-frequency complex phases, and meticulous quantum-kinetic analyses shortly thereafter put the formula on firmer ground. A new and more direct derivation is given here, showing that the equation follows almost immediately from a quantum relaxation-time approximation and an expansion in the mixing angle. Besides reproducing the desired result, the relaxation ansatz captures to a high degree of accuracy the interlaced dynamics of oscillations, decoherence, and plasma repopulation. Successes and limitations of the semiclassical equation are illustrated numerically and are shown to reflect the accuracy of the approximations employed in the derivation. The inclusion of interactions among the sterile neutrinos is also briefly addressed.


I. INTRODUCTION
Sterile neutrinos continue to be actively studied as sources of oscillation anomalies, as reconcilers of cosmic tensions, and as candidates for dark matter. In all these cases the cosmological abundance must be calculated, and so the dynamics of active-sterile mixing must be contended with. The essential challenge is that the full quantum-kinetic problem involves disparate time scales and the interplay of coherent (∝ G F ) and incoherent (∝ G 2 F ) effects. A Boltzmann-like equation is often used to calculate the nonthermal abundance of sterile neutrinos produced from active ones [1][2][3]: where f a(s) is the active (sterile) distribution function, Γ a is the scattering rate of active neutrinos, and θ m and ω m are the in-medium mixing angle and oscillation frequency. (We suppress dependence on momentum here and throughout.) Eq. (1) is a semiclassical approximation of the quantum kinetic equation (QKE) [4][5][6][7][8][9][10][11][12] i dρ dt for the density matrix ρ. Its computational appeal lies in the fact that, by packaging together the effects of the Hamiltonian H and collision term C as a single effective production rate, one can overlook the quantum phases and evolve only the classical densities. The first derivations of Eq. (1) (or variations of it) were based on single-particle arguments that equated the ν s production rate with the product of the ν a scattering rate and the probability of a ν a oscillating into a ν s [1,2]. Later analyses hearteningly arrived at similar * ljohns@physics.ucsd.edu formulas working from quantum-kinetic descriptions and judiciously applying approximations for the evolution in flavor space [13][14][15][16][17][18][19][20]. Our purpose here is to add another entry to the list, one that is complementary to the references just cited and whose virtue is the insight it gives into the quantum dynamics underlying the semiclassical behavior. Given the ongoing interest in sterile neutrinos, having a robust simplification of the quantum dynamics may prove useful for future applications. The guiding idea, which we support numerically, is that the evolution of ρ at small mixing angle is well described by the exponential decay of its deviations from equilibrium. As we demonstrate below, this simple ansatz leads promptly to Eq. (1).
In Sec. II we go through the derivation and discuss it in the context of other treatments. In Sec. III we present numerical comparisons of the Boltzmann and QKE solutions, highlighting the accuracy not only of Eq. (1) but also of the ansatz on which it is based. In Sec. IV we conclude.

II. DERIVATION
We begin by applying the quantum relaxation-time approximation to the collision term [5,16,21,22]: where Γ = (1/2) diag(Γ a , Γ s ) and ρ eq C is the H = 0 equilibrium. If the states do not mix, then the classical relaxation-time approximation is recovered, and the densities f a and f s approach their equilibrium values with time scales Γ −1 a and Γ −1 s respectively. We posit that the same approximation applies to the entire right-hand side of Eq. (2), with a single effective relaxation parameter replacing the individual scattering rates and the flavor equilibrium ρ eq F replacing the classical equilibrium. That is,  with Γ eff = (γ m /4) diag (1, 1). (The extra factor of 1/2 is added as a matter of preference.) Hence and in particular If the mixing angle is small, f eq a can safely be replaced in this equation by f a .
Using ρ = P 0 (1 + P · σ) /2, it follows from Eq. (5) that the polarization vector obeys At the same time, using H = (ω m /2) B m and setting Γ s = 0, Eqs. (2) and (3) imply where D = Γ a /2 is the decoherence rate and B m = sin 2θ m x − cos 2θ m z. The ansatz tells us that Eqs. (7) and Eqs. (8) can be set equal to each other at any moment in the evolution. For the sake of extracting γ m , we choose to equate them prior to significant sterile production, during which time P nearly equals z and P 0 and f a nearly equal f eq a . To first order in the deviations, P satisfies the eigenvalue equation Nontrivial solutions of Eq. (9) correspond to roots of the cubic equation Applying perturbation theory to zeroth order in sin 2 2θ m uncovers two of the roots, with values D±iω m . The third root, which is the purely real one that we seek, appears at first order: Plugging this into Eq. (6), we arrive at Eq. (1) as desired.
The analysis applies just as well to antineutrinos (or the right-handed states, if Majorana) as it does to neutrinos (or the left-handed states). If chemical potentials are involved in the dynamics, they are simply incorporated into the equilibrium distribution functions.
Eq. (9) was also considered in Ref. [4] (albeit not in the context of sterile neutrinos), Ref. [20] (albeit in a somewhat different form), and Refs. [13,16,17]. The last two were part of a series, along with Refs. [14,15,18,19], that provided significant insights into the dynamics of active-sterile oscillations. Vital to the derivation developed in those works is the approximation dP 0 /dt = 0, which was carefully shown in Ref. [19] to be justified despite its inconsistency with f a remaining near equilibrium during sterile production. We similarly find that the correct value of γ m is obtained by dropping the repopulation terms from the equation of motion obeyed by P, even though repopulation is crucial for accurately describing the evolution of the system as a whole. Our findings, based on the quantum relaxation-time approximation, are consistent in this regard with those of Ref. [19], based (in the words of the authors) on the "brute-force" approach.
Besides reproducing Eq. (1) with minimal effort, the preceding analysis also has the advantage of pointing to a deeper physical picture: the small-θ m dynamics is dominated by a flavor-space trajectory in which nonequilibrium deviations ρ eq F − ρ decay with (instantaneous) lifetime 2/γ m . Moreover, as we show in the next section, this approach sacrifices little in the way of accuracy for what it gains in simplicity.
But before moving on to the numerical analysis, let us briefly comment on the generalization to scenarios in which sterile neutrinos remain inert under the Standard Model couplings but interact via new ones. Reinstating Γ s = 0 in Eq. (3) leads to where D a,s = Γ a,s /2 and f eq s denotes the equilibrium that f s tends toward if the mixing is turned off. Because repopulation of f s drives P toward −z, the net effect of repopulation on P no longer takes the second-order form (Ṗ 0 /P 0 )(z − P).
However, suppose that Γ s is much faster than the active-sterile conversion rate. In that case f s is always very close to f eq s on the conversion time scale, and the new repopulation terms can once again be dropped from the equations so long as f s is consistently interpreted as being at the sterile-sector equilibrium value. The same relaxation ansatz can then be used as before, again leading to Eq. (9). The only difference now is that the decoherence rate in the expression for γ m should be interpreted as the sum D a + D s , in agreement with Refs. [23][24][25]. Note that the distinction between ρ eq C and ρ eq F is generally important to make, but is rendered moot when Γ s = 0.

III. NUMERICAL ANALYSIS
In this section we numerically study the validity of the quantum relaxation-time approximation, including the sterile neutrino Boltzmann equation implied by it. For simplicity we begin by assuming time-independent mixing and scattering parameters ω, sin 2 2θ, and D. We then go on to introduce a time-dependent potential and follow the system through resonance.  8)] for two choices of the ratio D/ω. In both cases the mixing angle is taken to be θ = π/100, and in the Boltzmann equation f a is always set equal to f eq a . The insets show that the solutions are discrepant at very early times before the QKE solution settles into the decay mode [Eq. (6)] on which the Boltzmann equation is based. Once it does so, both solutions grow linearly in time, as expected when f s f eq a . As production proceeds, the early-time discrepancy becomes less important as a fraction of the sterile abundance, and on times longer than γ −1 the evolution as a function of t = γt is virtually independent of the chosen  parameters: For cosmological applications it is typically undesirable to have the density of sterile neutrinos approach the thermal value, and the insets in Fig. 1 are therefore the relevant comparison. Achieving agreement on these shorter time scales requires θ to be small. Fig. 2 illustrates this point: the same quantities are plotted here as in Fig. 1, but now with θ = π/5. Early-time discrepancies are greatly exacerbated.
The relaxation ansatz asserts that f a should decay exponentially toward f eq a . Since a small mixing angle inhibits f a from ever deviating greatly from the equilibrium value, the ansatz also implies a delicate near-cancellation between the growth of P 0 and the decay of P = |P|. Fig. 3 verifies that both of these expectations are indeed borne out in the case D/ω = 10 −2 . The result is similar for stronger damping.
Interestingly, it was shown in Ref. [26] that the Boltzmann equation can be derived from the assumption that ρ as and ρ sa are both constant. Despite its expedience, that approximation is not an accurate description of how the active-sterile coherence develops, particularly on a γ −1 time scale. Fig. 4 shows that in fact the real part declines throughout production and is well fit by the relaxation ansatz. The imaginary part is similar.
Cosmological production of sterile neutrinos involves time-dependent parameters of course, and many scenarios of interest involve resonant mixing in particular [25,[27][28][29][30][31][32][33]. It is well-known that the Boltzmann equation is inadequate when the system passes adiabatically and coherently through a resonance [28,30,34]. To illustrate this claim, and to make a connection with the foregoing analysis, we add a potential V (t)z to ωB, with The adiabaticity parameter [34] is defined to be α = ωH sin 2θ tan 2θ, where H is the potential scale height, We also define κ = D −1 /H, or in terms of the adiabaticity, While α-proportional to the number of oscillation lengths that fit within a resonance width-sets the probability of a neutrino coherently transitioning between energy eigenstates at resonance, the coherence parameter κ-the number of mean free paths that fit within a resonance width-indicates the degree to which scattering affects the evolution.
In Fig. 5 we compare the Boltzmann and QKE solutions for several choices of D/ω and α. Each panel on the right shows f a (t) calculated with the same parameters used in the panel to its left. The extent to which f a dips away from the equilibrium value at resonance indicates the inapplicability of the relaxation ansatz during this period of production. The top three rows show nonadiabatic resonances with increasing rates of decoherence. The bottom row shows an adiabatic resonance with a very large decoherence rate. Regardless of the adiabaticity, production through the resonance is flattened out as D/ω increases, and the accuracy of the Boltzmann solution improves. The key to the Boltzmann equation being a valid approximation is always that κ be small, to ensure that the relaxation ansatz for ρ is not too badly violated at resonance.

IV. SUMMARY
A few different derivations of the sterile neutrino Boltzmann equation can be found in the literature. Aside from its simplicity, the relaxation-time approach is notable in that it is based on an accurate description of the full quantum dynamics of active-sterile mixing. Replicating Eq. (1) is one implication, but the approximation describes equally well the repopulation of the active species and the decay of coherence.
We have numerically shown that the limitations of Eq. (1) reflect the assumptions used here to derive it: the Boltzmann solution works best when θ m is small, and it comes up short when resonance coherently steers the system away from relaxation. Cosmological sterile-neutrino production is just one scenario of interest that is covered, at least partially, by this range of validity. A similar analysis can be applied to the mixing of active states in the presence of unequal chemical potentials-a circumstance that is realized in supernovae, binary neutron-star mergers, and possibly the early universe. In light of the computational obstacles faced in these environments, there is a clear need for good quantum-kinetic approximations.