Sterile Neutrinos as Dark Matter: Alternative Production Mechanisms in the Early Universe

: We study various production mechanisms of sterile neutrinos in the early universe beyond and within the standard model. We obtain the quantum kinetic equations for production and the distribution function of sterile-like neutrinos at freeze-out, from which we obtain free streaming lengths, equations of state and coarse grained phase space densities. In a simple extension beyond the standard model, in which neutrinos are Yukawa coupled to a Higgs-like scalar, we derive and solve the quantum kinetic equation for sterile production and analyze the freeze-out conditions and clustering properties of this dark matter constituent. We argue that in the mass basis, standard model processes that produce active neutrinos also yield sterile-like neutrinos, leading to various possible production channels. Hence, the ﬁnal distribution function of sterile-like neutrinos is a result of the various kinematically allowed production processes in the early universe. As an explicit example, we consider production of light sterile neutrinos from pion decay after the QCD phase transition, obtaining the quantum kinetic equation and the distribution function at freeze-out. A sterile-like neutrino with a mass in the keV range produced by this process is a suitable warm dark matter candidate with a free-streaming length of the order of few kpc consistent with cores in dwarf galaxies.


Introduction
In the concordance ΛCDM standard cosmological model, dark matter (DM) is composed of primordial particles, which are cold and collisionless [1,2]. In this cold dark matter (CDM) scenario, structure formation proceeds in a hierarchical "bottom up" approach: small scales become non-linear and collapse first, and their merger and accretion leads to structure on larger scales. CDM particles feature negligible small velocity dispersion, leading to a power spectrum that favors small scales. In this hierarchical scenario, dense clumps that survive the merger process form satellite galaxies. A popular dark matter candidate is a weakly-interacting, cold thermal relic (WIMP) [3]; however, decades long experimental searches for this candidate have not yet provided compelling observational evidence in its favor. Alternatives to the WIMP paradigm are the focus of current studies [4][5][6][7] due to several observational shortcomings of the standard cold dark matter cosmology (ΛCDM) at small scales. N-body simulations of cold dark matter produce dark matter profiles that generically feature cusps yet observations suggest a smooth-core profile [8][9][10] (core-cusp problem). The same type of N-body simulations also predict a large number of dark matter dominated satellites surrounding a typical galaxy, which is inconsistent with current observations [11][12][13] (missing satellites problem). Both the missing satellites and core-cusp problem can be simultaneously resolved by allowing some fraction of the dark matter to be "warm" (WDM) [14][15][16][17][18][19] with a massive "sterile" neutrino being one popular candidate [20][21][22][23][24]-other examples include Kaluza-Klein excitations from string compactification or axions [4,25], although baryonic physics may also be part of the solution. The "hotness" or "coldness" of a dark matter candidate is discussed in terms of its free streaming length, λ f s , which is the cut-off scale in the linear power spectrum of density perturbations. Cold dark matter (CDM) features small ( pc) λ f s that brings about cuspy profiles, whereas WDM features λ f s ∼ few kpc, which would lead to cored profiles. Recent (WDM) simulations including velocity dispersion suggest the formation of cores, but do not yet reliably constrain the mass of the (WDM) candidate in a model independent manner since the distribution function of these candidates is also an important quantity which determines the velocity dispersion and thereby the free streaming length [26][27][28].

Cosmological Consequences of Decoupled Particles
While the study of kinetics in the early universe is available in the literature [1,72,73], in this section we summarize the cosmological consequences of decoupled particles [74][75][76][77][78][79] highlighting several aspects relevant to the next sections.
We consider a spatially flat FRW cosmology with length element The local dispersion relation for a particle is given by E(t) = p 0 (t) = m 2 + P 2 f (t) (2) where P f (t) = p c a(t) , and p c is the comoving momentum. A frozen distribution describing a particle that has been decoupled from the plasma is constant along geodesics, therefore, taking the distribution to be a function of the physical momentum P f and time, it obeys the Liouville equation or collisionless Boltzmann equation Obviously, a solution of this equation is where p c is the time independent comoving momentum. The physical phase space volume element is invariant, d 3 X f d 3 xP f = d 3 x c d 3 p c , where f , c refer to physical and comoving volumes, respectively.
In what follows, we consider general distributions as in Equation (5) unless specifically stated.
The kinetic energy momentum tensor associated with this frozen distribution is given by where g is the number of internal degrees of freedom, typically 1 ≤ g ≤ 4. Taking the distribution function to be isotropic, it follows that (8) where ρ is the energy density and P is the pressure. In summary, The pressure can be written in a manner more familiar from kinetic theory as where v f = P f /E is the physical (group) velocity of the particles measured by an observer at rest in the expanding cosmology. Covariant energy conservation follows fromṗ 0 = −H(t) P 2 f /p 0 , leading tȯ The number of particles per unit physical volume is and obeys dn(t) dt namely, the number of particles per unit comoving volume n(t) a 3 (t) is conserved. These are generic results for the kinetic energy momentum tensor and the particle density for any distribution function that obeys the collisionless Boltzmann Equation (4).
It is convenient to introduce the dimensionless ratios and Changing the integration variable in Equations (9)- (13) to P f = y T d (t), we find leading to the equation of state: In the ultrarelativistic and non-relativistic limits, x → 0 and x → ∞, respectively, we find In the ultrarelativistic limit, the energy density and pressure become describing radiation behavior. In the non-relativistic limit and the equation of state becomes corresponding to cold matter behavior. In the non-relativistic limit, it is convenient to write where ζ(3) = 1.2020569 . . . , g d is the number of ultrarelativistic degrees of freedom at decoupling, and n γ (t) is the photon number. The average squared velocity of the particle is given in the non-relativistic limit by Therefore, the equation of state in thermal equilibrium is given by where σ is the one-dimensional velocity dispersion given at redshift z by where, in general, f d depends on the mass of the particle. For a particle that decouples while non-relativistic, this is recognized as the generalization of the Lee-Weinberg lower bound [73,[80][81][82][83], whereas if the particle decouples in or out of LTE when it is ultrarelativistic, f d,a (y) does not depend on the mass. Equation (30) provides an upper bound, which is a generalization of the Cowsik-McClelland [73,84] bound. For the particle to be a suitable dark matter candidate, the free streaming length, namely the cutoff in the dark matter power spectrum, must be much smaller than the Hubble radius, this implies that the velocity dispersion must be small, namely the particle must be non-relativistic From Equation (24), this constraint yields where T d (t) is given by Equation (15). From Equations (16), (27) and (32), we obtain the following condition for the particle to be non-relativistic at redshift z m 2.958 Taking the relevant value of the redshift for large scale structure to be the redshift at which reionization occurs z s ∼ 10 , we find the following generalized constraint on the mass of the particle of species (a), which is a dark matter component 2.958 The left side of the inequality corresponds to the requirement that the particle be nonrelativistic at reionization (taking z s ∼ 10), namely a small velocity dispersion V 2 /c 2 1, corresponding to a free streaming length λ f s ∼ V 2 /c 2 d H much smaller than the Hubble radius (d H ), while the right hand side is the constraint from the requirement that the decoupled particle is a dark matter component, namely Equation (30) is fulfilled.

Free Streaming Length
The free streaming length gives an indication of the cut-off in the matter power spectrum.
The comoving free streaming wavevector is defined in analogy with the Jeans wavevector in the fluid description of perturbations, namely where V 2 (η) is given by Equation (24), which we have written as The cutoff scale in the power spectrum is the comoving free streaming length, namely where d H = 1/H 0 = 3 Gpc/h is the Hubble distance. This definition differs from the usual definition of the comoving free streaming distance l f s during matter domination by factors of O(1): During the matter dominated era, it follows that hence the free streaming distance from matter-radiation equality until a 0 O(1) is given by

Coarse Grained Phase Space Densities
In their seminal article, Tremaine and Gunn [85] argued that the coarse grained phase space density is always smaller than or equal to the maximum of the (fine grained) microscopic phase space density, which is the distribution function. A similar argument was presented by Dalcanton and Hogan [9,10].
However, it is straightforward to show that the resulting coarse grained phase space density is not a Liouville invariant. Instead, we define the coarse grained (dimensionless) primordial phase space density D ≡ n(t) which is Liouville invariant and where P 2 f is defined in Equation (24). Since the distribution function is frozen and is a solution of the collisionless Boltzmann (Liouville) Equation (4), it is clear that D is a constant, namely a Liouville invariant in the absence of self-gravity. Explicitly, D is given by When the particle becomes non-relativistic ρ(t) = m n(t) and V 2 = where Q DH is the phase-space density introduced by Dalcanton and Hogan [9,10] and the one-dimensional velocity dispersion σ is defined by Equation (25).
In the non-relativistic regime, D is related to the coarse grained phase space density Q TG introduced by Tremaine and Gunn [85] The observationally accessible quantity is the phase space density ρ/σ 3 , therefore, using ρ = m n for a decoupled particle that is non-relativistic today and Equation (25), we define the primordial phase space density where we found that keV 4 km/s 3 = 1.2723 10 8 M kpc 3 . During collisionless gravitational dynamics, phase mixing increases the density and velocity dispersions in such a way that the coarse grained phase space density either remains constant or diminishes, namely where D is given by Equation (42) for an arbitrary distribution function. For a particle that decouples when it is ultrarelativistic, D does not depend on the mass, hence Equation (47) yields a lower bound on the mass of the particle directly from the observed phase space density and the knowledge of the distribution function.
The results obtained in this section are of paramount importance to assess the cosmological consequences of any dark matter candidate from its frozen distribution function.

Production from Scalar Decay
As an explicit example, in this section we study [47] the model presented in references [48,49,86,87] as a simple extension of the minimal standard model with only one sterile neutrino; however, including more species is straightforward. The Lagrangian density is given by where L SM is the standard model Lagrangian, L α ; α = 1, 2, 3 are the standard model SU(2) lepton doublets, ν is a singlet sterile neutrino with a (Majorana) mass m, a real Higgs-like scalar χ with Yukawa coupling Y to the sterile neutrino, which in turn is Yukawa coupled to the active neutrinos via the Higgs doublet H, thereby building a see-saw mass matrix in terms of the vacuum expectation of this doublet [48,86,87]. For a sterile neutrino in the mass range keV, we consider the vacuum expectation value and mass M of the singlet scalar χ in the range ∼100 GeV as discussed in references [48,49,86]. If the sterile neutrino mass m s ∼ keV results from the vacuum expectation value χ ∼ 100 GeV, then the Yukawa coupling Y ∼ 10 −8 . In this section, we study the production of sterile neutrinos via the decay of the gauge singlet scalar field χ with Y 1 and M m, assuming that the scalar field χ is strongly coupled to the plasma and is in (LTE) with a Bose-Einstein distribution function where k is a comoving wavevector and where T 0 = T d a(t d ) would be the temperature of the plasma today. During the radiation dominated era the Hubble expansion rate is given by [73] H(t) 1.66g where g(t) is the effective number of relativistic degrees of freedom. The first step is to obtain the quantum kinetic equation for production.

Quantum Kinetic Equation for Production
With the Yukawa interaction L I = Y χ ν ν and χ a scalar field with mass M and ν either a Dirac or Majorana fermion field of mass m the quantum kinetic equation is obtained just as in Minkowski space time by obtaining the total transition probabilities per unit time of the decay and inverse decay processes.
The quantum kinetic equation for the neutrino population n(p; t) (and similarly for the antineutrino) is of the usual form where the gain and loss terms are obtained from the corresponding transition probabili- The gain term is determined by the decay reaction χ → ν ν corresponding to an initial state with N k quanta of the scalar field χ and n p,s , n q,s quanta of neutrinos and antineutrinos, respectively, with N k − 1, n p,s + 1, n q,s + 1 quanta in the final state. The corresponding Fock states are The loss term is determined by the inverse decay reaction ν ν → χ corresponding to an initial state with N k quanta of the scalar field χ and n p,s , n q,s quanta of neutrinos and antineutrinos, respectively, with N k = 1, n p,s − 1, n q,s − 1 quanta in the final state. The corresponding Fock states are The calculation of the matrix elements is standard, the fields are quantized in a volume V in terms of Fock creation-annihilation operators, with the corresponding spinor solutions for the neutrino fields. To lowest order in Y we find similarly, for the loss term we obtain, the frequencies Summing the |M f i | 2 over k, q, s, s and taking the infinite volume limit, we find the total transition probability for the gain term ∑ k, q,s,s where T is the total reaction time. For the loss term we find the same expression, but with the replacement N p+ q → 1 + N p+ q , (1 − n p ) → n p , (1 − n q ) → n q . Carrying out the angular integral, we obtain the final expression for the neutrino production rate where q ± are the roots of the equations In the limit M m we find these roots to be The extension to the cosmological case replaces the momenta by the physical momenta and for Majorana neutrinos n = n.
For Y 2 1 we expect that neutrinos will decouple early and their distribution function will freeze-out with n p , n p 1, an assumption that will be confirmed selfconsistently below.
Neglecting the neutrino population buildup in (59), and neglecting terms of order m 2 /M 2 1, taking the scalar field to be in LTE and replacing the momenta in (59) by their physical values, yields where in terms of the corresponding physical momenta. These values are determined by the kinematic thresholds for decay. For Y 1, neutrinos decouple at temperatures T d m, namely when they are still relativistic, therefore we can safely neglect terms of order m 2 /T 2 (t) < m 2 /T 2 d

1.
Under this assumption (to be confirmed self-consistently below), using P f (t)/T(t) = p/T 0 and neglecting terms of order m 2 /M 2 1, the kinetic equation above simplifies. Let us introduce the variables and leading to the following form of the quantum kinetic equation, where and the effective number of relativistic degrees of freedom depends on time through the temperature. Under the assumption M m, for example taking M ∼ 100 GeV, m ∼ keV, the exponential in the numerator inside the logarithm can be neglected for all y m 2 /M 2 ∼ 10 −16 . Although we are interested in the small momentum region of the distribution function (small y), the phase space suppression for small momentum in the integrals of the distribution function entails that we can safely neglect the contributions of such small values, therefore neglecting the numerator inside the logarithm in (69). To integrate the rate Equation (69), we must provide g(t). In the standard model, this function is approximately constant in large intervals and features sharp variations in the regions of temperature when relativistic degrees of freedom either decay, annihilate or become non-relativistic [73].
We will assume that g(t) is approximately constant in the region of (large) temperature before and during decoupling, replacing the value of g by its average g over the range in which the rate (69) is appreciable. Therefore, we replace With these approximations, Equation (69) is equivalent to the kinetic equation of references [48,87] and can be integrated exactly. We obtain, and Γ[a, b] is the incomplete Gamma function. The rate (69) vanishes as τ → 0, reaches a maximum and falls-off exponentially as τ/4y → ∞. The maximum rate is larger for smaller values of y, and this feature implies an enhancement of the distribution function at small momenta, which is important for the cosmological consequences studied in ref. [47].
The asymptotic behavior of the distribution function (72) for τ/4y 1 is given by n(y; τ) For y 1, the distribution function is largest, reaching its asymptotic value for τ 1, whereas for large values y 1, the distribution function is strongly suppressed and reaches its asymptotic form much later.
Hence, for small y, which is the region of the distribution function most relevant for small scale structure formation, decoupling occurs fairly fast at a decoupling temperature ∼ M m s . The dark matter transfer function depends more sensitively on the small momentum (small y) region of the distribution function. The analysis above shows that for y < 1, the distribution function freezes out on a "time" scale τ 1, hence a temperature scale T(t f ) ∼ M, where t f is the "freeze-out" time. Hence, if the effective number of relativistic degrees of freedom g(t) varies smoothly within the temperature range in which the population freezes-out ∼ M the approximation (71) is justified and reliable. Therefore, we take the value g as the effective number of relativistic degrees of freedom at decoupling, for example, for a freeze-out temperature T(t f ) ∼ M ∼ 100 GeV in the standard model g ∼ 100.
The distribution function at freeze-out is obtained by taking the τ/4y → ∞ limit in Equation (72), hence we introduce the distribution function at decoupling This is the unperturbed distribution function that enters in the Boltzmann-Vlasov equation that determines the evolution of density and gravitational perturbations, and is displayed in Figure 1. It is remarkable that for y 1 in striking contrast to the thermal Fermi-Dirac distribution function. The enhancement for small y indicates that this is a cold dark matter species.

The Cutoff in the Power Spectrum: Free Streaming Length
The comoving free-streaming wavevector plays the same role as the comoving Jeans wavevector in a fluid, but with the speed of sound replaced by the velocity dispersion of the decoupled particle (see Section 2), it is given by and its value today is where is the three dimensional velocity dispersion of the non-relativistic particles today and we introduced Using the relation where T cmb = 2.348 × 10 −4 eV is the temperature of the cosmic microwave background today, with g d = g is the effective number of relativistic degrees of freedom at decoupling, and Ω M h 2 = 0.105 for non-baryonic dark matter, leads to For the non-equilibrium normalized distribution function (75), we find Taking the number of relativistic degrees of freedom at decoupling g d = g (see Equation (71) and preceding discussion), we find leading to free streaming wavevector and wavelength at matter-radiation equality For comparison, the value of k f s (0) for a fermion that decoupled with a thermal relativistic relativistic Fermi-Dirac distribution function is The smaller velocity dispersion in the case of the non-equilibrium distribution function yields a ∼30% increase in the free streaming wavevector, and consequently, a decrease in the free streaming length. However, just as importantly in enhancing k f s (0) is decoupling at high temperature for a larger value of g d . A sterile neutrino produced via the Dodelson-Widrow [21] mechanism at decoupling temperature T d ∼ 150 MeV; g d ∼ 30 [21] yields a free streaming wavevector at matter-radiation equality λ (dw) f s (t eq ) ∼ 900 kpc.

Production from Standard Model Processes
Sterile neutrinos are SU(2) weak singlets that do not interact via standard model weak interactions, and couple to the flavor active neutrinos via an off diagonal mass matrix. Different models propose different types of (see-saw like) mass matrices. For our purposes, the only relevant aspect is that the diagonalization of the mass matrix leads to mass eigenstates, and these are ultimately the relevant degrees of freedom that describe dark matter in its various possible realizations, cold, warm or hot. Cold and warm dark matter would be described by heavier species, whereas the usual (light) active-like neutrinos could be a hot dark matter component depending on the absolute value of their masses. In this section, we focus on the production of sterile-like neutrinos from standard model couplings: in the mass basis, neutrinos couple via standard model charged and neutral current interactions, albeit with very small mixing matrix elements. If kinematically allowed, the same processes that produce active neutrinos will produce sterile-like neutrinos with much smaller branching ratios determined by the effective mixing angles.
Upon diagonalization of the mass matrix, the active (massless) flavor neutrino fields are related to the neutrino fields that create/annihilate mass eigenstates as H αs ν s,L (x) (88) where α = e, µ, τ, m = 1, 2, 3 refer to the light (atmospheric, solar) mass eigenstates with massess M m and s = 4, 5 · · · refer to sterile-like mass eigenstates of mass M s . In the mass basis the full Lagrangian density is where L mSM is the standard model Lagrangian augmented with a diagonal neutrino mass matrix for the active like neutrinos ν m of masses M m (m = 1, 2, 3), in this Lagragian density the charged and neutral interaction vertices are written in terms of the neutrino mass eigenstates with the mixing matrix U αm , L 0s is the free field Lagrangian density for the sterile-like neutrinos ν s of masses M s (s = 4, 5 · · · ) and to linear order in H αs where with P L = (1 − γ 5 )/2. From now on, we refer to sterile-like neutrinos instead of "sterile" neutrinos, because unlike the concept of a sterile neutrino, the sterile-like mass eigenstates neutrinos do interact with standard model degrees of freedom via charged and neutral current vertices.
In terms of the interaction vertices (90), we can obtain the quantum kinetic equations for production of massive neutrinos to order (|H αs |) 2 from the master equation of the form gain-loss described in the previous section, namely where n s (q; t) is the distribution function of the sterile neutrino. This is calculated for each process: the gain term describes the increase in the population n s (q; t) by the creation of a sterile-like state and the loss by the annihilation. We first obtain the gain and loss terms in Minkowski space-time, and follow by replacing the respective momenta and temperatures by the physical quantities in the expanding cosmology.
There are important loop corrections at finite temperature, self-energy corrections to the incoming and outgoing external "legs" as well as vertex corrections. There are also finite temperature corrections to the mixing angles arising from self-energy loops, and these tend to suppress the mixing matrix elements [36][37][38][88][89][90]; therefore, in medium the matrix elements H sl → H sl,e f f (T) and typically in absence of Mikheyev-Smirnov-Wolfenstein (MSW) resonances H sl,e f f (T) < H sl . At high temperatures, there are hard thermal loop corrections to the self-energy of fermions and vector bosons [91][92][93][94][95][96][97][98][99][100][101][102][103] that lead to novel collective excitations with masses ∝ gT where g is the gauge coupling. Photons and fermion-antifermion pairs form plasmon collective excitations with mass ∝ eT. For T > T EW , the W, Z vector bosons do not acquire a mass via the Higgs mechanism because the ensemble average of the Higgs field vanishes, but they acquire thermal masses ∝ g w T akin to the plasmon collective excitations. Thermal masses for collective excitations of W, Z could open up kinematic windows for decay into ν s above the electroweak transition.
The quantum kinetic equations from standard model vertices had been obtained in refs. [66,67], and the finite temperature effects on mixing angles and quantum kinetics at the electroweak scale have been studied in refs. [90,104]. Their general form is where the gain and loss rates are obtained as in Minkowski space time replacing the momenta and energies by the local time dependent counterparts in the expanding cosmology and the (LTE) distribution functions for the various bosonic or fermionic species (with vanishing chemical potentials). As a consequence of the linearity of the quantum kinetic Equation (93), which, in turn is a consequence of keeping only terms of O(|H αs | 2 ), O(| H ms | 2 ), the full quantum kinetic equation is a simple sum over all possible channels with total gain and loss rates therefore the gain and loss rates in the quantum kinetic Equation (93) are the total rates (94). For a general discussion of the various standard model processes, see ref. [67], and in refs. [90,104] production from vector boson decay near the electroweak scale T EW 100 GeV was studied in detail, primarily focused on heavier sterile-like neutrinos. These studies revealed novel resonances in absence of lepton asymmetries and unexpected helicity dependence on the production rates. The reader is referred to these articles for further details.
The main point of this discussion is that the final distribution function at decoupling for sterile-like neutrinos is the result of all the kinematically allowed processes that produced these species, hence such distribution function is a mixture of the various production channels.

Sterile Neutrinos from Pion Decay
It is clear from the discussion above that standard model interaction vertices that lead to the production of active neutrinos will also lead to the production of the sterile-like neutrinos as long as the processes are kinematically allowed. This implies a far broader range of production mechanisms than those that had been the focus in the literature. A reliable calculation of the production rate requires the inclusion of finite temperature corrections to masses and interaction vertices. In this section, we consider the specific example of production of sterile-like neutrinos from pion decay shortly after the QCD phase transition (crossover) into the confined phase.
The QCD phase transition occurs at a temperature scale T QCD 155 MeV at a time scale t 10 µsecs, below T QCD quarks and gluons are confined into hadrons, the lightest of which are π-mesons (pions). Therefore, these are the degrees of freedom that feature the largest population at temperatures below T QCD . Charged pions decay via charged current weak interactions into charged leptons and neutrinos, therefore, as per the discussion in the previous section, this decay process will also contribute to the production of sterile-like neutrinos from the process π + → µ + ν s via the mixing matrix element. In fact, in terrestrial accelerator experiments, charged pion decay is one of the primary sources of neutrinos.
In ref. [40][41][42][43], the authors proposed to study hadronic contributions to the production of sterile neutrinos by considering the self-energy corrections to the active (massless) neutrinos in terms of correlators of vector and axial-vector currents. While this is correct in principle, it is an impractical program: the confined phase of QCD is strongly coupled, and the description in terms of nearly free quarks is at best an uncontrolled simplification, casting doubts on the reliability of the conclusions in this reference.
Instead, in this section, we study the production from pion decay by relying on the well understood effective field theory description of weak interactions of pions, enhanced by the results of a systematic program that studied finite temperature corrections to the pion decay constant and mass implementing non-perturbative techniques from chiral perturbation theory, linear and non-linear sigma models and lattice gauge theory [105][106][107][108][109][110][111][112][113][114][115][116][117][118][119][120]. There are at least three important reasons that motivate this study: (1) it is a clear and relevant example of the quantum kinetic equation approach to production and freeze out that includes consistently finite temperature corrections. (2) Pion decay surely contributed to the production of sterile-like neutrinos in the early universe if such species of neutrinos do exist: the QCD transition to its confined phase undoubtedly happened, and the consequent hadronization resulted in baryons and low lying mesons, pions being the lightest, are the most populated degrees of freedom near the QCD scale. (3) The two leptonic decay channels π → µν s ; π → eν s feature different kinematics and thresholds, in particular the (V-A) vector-axial coupling results in chiral suppression of the e channel for vanishingly small M s . Since the neutrino produced in the decay is kinematically entangled with the emitted lepton [121], we expect that the distribution functions from each channel will display differences as a manifestation of this kinematic entanglement and dependence on the particular production channel.
At temperatures larger than the QCD phase transition T QCD ∼ 155 MeV [117][118][119][120], quarks and gluons are asymptotically free. Below this temperature, QCD bound states form on strong interaction time scales, the pion being the lightest. Lattice gauge theory calculations [117][118][119][120] suggest that the confinement-deconfinement transition is not a sharp phase transition, but a smooth yet rapid crossover at a temperature T c ∼ 155 MeV within a temperature range ∆T ± 10 MeV. This occurs in the radiation dominated epoch at t 10 µs within a time range ∆t ∼ 2-3 µs. This is much larger than the typical strong and electromagnetic interaction time scales 10 −22 s implying that pions that form shortly after the confining cross-over are brought to LTE via strong, electromagnetic and weak interactions on time scales much shorter than ∆t. After pions are produced, they reach LTE via π − π scattering on strong interaction time scales, then decay into leptons and maintain LTE via detailed balance with the inverse process since charged leptons and active-like neutrinos are also in LTE. In conclusion, for T 155 MeV, pions are present in the plasma in thermal/chemical equilibrium due to pions interacting on strong interaction time scales (10 −22 s) while the crossover transition occurs on the order of 2-3 µs.
If the pions (slowly) decay into heavy neutrinos ν s with very small mixing angles, detailed balance for π lν s will not be maintained if the sterile-like neutrinos are not in LTE. As the pion is the lowest lying bound state of QCD, it is a reasonable assumption that during the QCD confinement-deconfinemet crossover pions will be the most dominantly produced mesonic bound state and, during this time, pions will remain in LTE with the light active-like neutrinos and charged leptons by detailed balance π l α ν α . We focus on sterile-like neutrino ν s production from π → lν s , which is suppressed by the mixing matrix element |H αs | 2 1 with respect to the active neutrinos, and this channel does not maintain detailed balance.
A full study of sterile neutrino production through π decay in the early universe requires various finite temperature corrections, and a preliminary study that focused on the production of neutrinos in the keV mass range [66] has provided a thorough study, the results of which suggest a suitable warm dark matter candidate with free streaming lengths on the order of several kpc. Furthermore, this production mechanism may yield a substantial (if not the dominant) contribution to DM below the QCD scale. We summarize several of these results in this section.

Quantum Kinetic Equation for π → l ν s
It is generally accepted that in the early universe, where temperatures and densities are larger than the QCD scale (∼155 MeV), quarks and gluons are asymptotically free forming a quark-gluon plasma. As the universe expands and cools, quarks and gluons undergo two phase transitions: deconfinement/confinement and chiral symmetry breaking. Confinement and hadronization result predominantly in the formation of baryons and mesons, the lightest of which-the pions-are dominant and are the pseudo Goldsone bosons associated with chiral symmetry breaking. A recent lattice QCD calculation [117][118][119][120] suggests that this phase transition is not first order but a rapid crossover near a critical temperature T QCD ≈ 155 MeV. Pions thermalize in the plasma via strong, electromagnetic and weak interactions and are in local thermodynamic equilibrium. Their decay into leptons and active neutrinos is balanced by the inverse process as the leptons and active neutrinos are also in LTE, therefore pions remain in LTE with active neutrinos by detailed balance. However if the pions (slowly) decay into sterile neutrinos, detailed balance will not be maintained, as the latter are not expected to be in LTE.
We focus on sterile neutrino ν s production from π → lν s , which is suppressed by |H ls | 2 1 with respect to the active neutrinos and does not maintain detailed balance. We also restrict the analysis to a scenario with no lepton asymmetry, which sets the chemical potential of pions and leptons to zero. The interaction Hamiltonian responsible for this decay is where J π σ = i∂ σ π(x, t) is the pseudoscalar pion current. The time evolution of the population of sterile-like neutrinos can be described via a quantum kinetic equation that takes the usual form of Gain-Loss where the gain and loss contributions are obtained from the appropriate transition matrix elements |M f i | 2 obtained from the interaction Hamiltonian (95) (for details see ref. [66]). The quantum kinetic equation is obtained in ref. [66], it is given by where V ud are the elements of the CKM matrix, G F is Fermi's constant, f π the pion decay constant and p ± are obtained from the constraint This gives the solutions Note that these bounds coalesce at when m 2 π − m 2 l + m 2 s = 2m 2 π m 2 s and the rate, Equation (96), vanishes simply because this corresponds to the kinematic threshold. As discussed in the previous section, these results are extended to an expanding cosmology by replacing the momentum with the physical momentum, q → Q f = q c /a(t) and the temperature by T 0 /a(t).

Non-Thermal Sterile Neutrino Distribution Function
A body of work [107][108][109][110][111][112][114][115][116] has established that, when π's are present in the medium in LTE, the π decay constant, f π , and π mass vary with temperature for T T QCD , where T QCD is the critical temperature for the QCD phase transition. We account for these effects and make several simplifications by implementing the following: • The finite-temperature pion decay constant has been obtained in both non-linear sigma models [108][109][110][111] and Chiral perturbation theory [107,112,[114][115][116] with the result given as This result is required in the quantum kinetic equation since production begins near T QCD ∼ 155 MeV and continues until the distribution function freezes out. We assume prior to T QCD that there are no pions and that hadronization happens instantaneously at T ∼ T c ∼ 155 MeV. • The mass of the pion varies with temperature as described in detail in refs. [107,115,116]. The finite temperature corrections to the pion mass is calculated with electromagnetic corrections in chiral perturbation theory and its variation with temperature is shown in Figure 2 of [107]. In these references, it can be seen that between 50 and 150 MeV the pion mass only varies between 140 and 144 MeV. Since this change is so small, we neglect the temperature variation in the pion mass and simply use its average value: m π (T) = 142 MeV (see fig in ref. [107]). • We assume that the lepton asymmetry in the early universe is very small so that we may neglect the chemical potential in the distribution function of the pions and charged leptons. This asymmetry is required for the Shi-Fuller mechanism, but will not be present in these calculations. We will show a similar enhancement at low moment to SF, but the enhancement found here is with zero lepton asymmetry. • With the assumption that there is no lepton asymmetry, the contribution to thermodynamic quantities from π − → lν will be equal to that of π + →lν. In which case, the degrees of freedom will be set at g ν = 2 accounting for both equal particle and antiparticle contributions in the case of Dirac fermions and the two different sources (π ± ) for Majorana fermions. The different helicities have already been accounted by summing over spins in the evaluation of the matrix elements of the previous section. • We assume that there had been no production of sterile neutrinos prior to the hadronization period from any other mechanisms (such as scalar decays or DW). This allows us to set the initial distribution of the sterile neutrinos to zero in the kinetic equation, which implies that our results for the distribution function will be a lower bound for the distribution function. Any other prior sources could only enhance the population of sterile neutrinos. By neglecting the initial population, we can neglect the Pauli blocking factor of the ν's in the production term, and we can also neglect the loss term (see discussion below).
After the QCD phase transition, there is an abundance of pions present in the plasma in thermal/chemical equilibrium. The pions will decay, predominantly via π ± → l ± ν s (ν s ), producing sterile neutrinos which, assuming that sterile neutrinos had not been produced up to this point, will have a negligible distribution function. The reverse reaction (lν s → π) will not occur in any significant quantities also due to the assumption of null initial population and |U ls | 2 1; under these assumptions, we may neglect the loss terms in the kinetic equation. With these assumptions, we use the following distributions for the production terms in the quantum kinetic equation where k c is a comoving momentum as discussed in Section 2.
With these replacements, neglecting the loss terms and setting E l (p, q) = E π (p) − E s (q), the quantum kinetic equation becomes e −E ν (q)/T e E π (p)/T (e E π (p)/T − 1)(e −E s (q)/T e E π (p)/T + 1) where the limits of integration are given by Equation (98) and we have suppressed the dependence of physical momentum on time. The integral can be done by a simple substitution with the final result given here: where p ± are given by Equation (98). Making the following change of variables where T 0 is the temperature of the plasma today since the normalization is set by a(t 0 ) = 1.
The QCD phase transition occurs during the radiation dominated era so that the Hubble factor is given by Equation (51), inserting it into Equation (102) suggests the definition During the period shortly after hadronization when m µ T m π the relativistic degrees of freedom are g(t) ∼ 14.25 while in the regime m e T m µ the degrees of freedom count is g(t) ∼ 10.75 [73]. We expect the sterile neutrino decoupling to happen well above the electron mass (this will be justified later) and since the variation of g(t) is small we replace it with its average value, g(t) ∼ḡ = 12.5, so that we can neglect the time dependence of Λ.
These substitutions and variable changes lead to a more tractable form of the kinetic equation The population build up is obtained by integrating where we have neglected any early population of ν s and the value of τ 0 is determined by when the pions are produced, assumed almost immediately after the hadronization transition. Our assumption is that this happens instantaneous at the QCD phase transition and the pions reach equilibrium instantaneously. This is justified by the results of [117][118][119][120], which suggest a continuous transition that allows for thermalization on strong interaction time scales. As shown in [117][118][119][120], the continuous phase transition occurs at T QCD ∼ 155 MeV so that τ 0 = m π /T QCD = 0.92 ≈ 1. As we set τ 0 below this value, we expect that the population would increase, as there will simply be more time for production to occur; this will be confirmed in a subsequent section.
The rate equation and the resulting population buildup as a function of τ for different values of y has been studied in detail in ref. [66]; however, an estimate of the sterile neutrino decoupling temperature can be made by considering the pion distribution. As the plasma temperature cools to well below the pion mass, the pion distribution will go as f π = e −m π /T(t) , leading to a large suppression of the production rate at T 10 MeV, which is when we expect the sterile neutrinos to freeze out. This is indeed what is found numerically in the population build up calculations of ref. [66]. The light mass limit relevant for sterile-like neutrinos with m s keV can be studied analytically. As discussed above, we expect freeze out to occur on the order of T f ∼ 10-15 MeV and we will consider here light sterile neutrinos with m s O (MeV). Restricting attention to this mass range sets m s /T f 1 for the duration of sterile neutrino buildup and simplifies the kinetic equation considerably. For this particular production mechanism, it follows that m 2 s m 2 π − m 2 l and we introduce the parameter so that, upon expanding in small parameters m s /T f and m 2 s /(m 2 π − m 2 l ) lead to the following simplifications which is relevant for a wide range of light steriles that freeze out at m s /T f 1. Note that we are suppressing the m s dependence of ∆ and will do so for the remainder of this work (for m s 1 MeV). In this limit, the kinetic equation simplifies to We must ensure that the rate remains small in order to ignore the sterile's population build up and consequent Pauli blocking. The population scales with Λ, and if one were to compute the next order correction by including the first order buildup in the rate equation, the higher order correction would scale as Λ 2 and so on. Provided Λ 1 (discussed below), the first order correction will be sufficient and higher order perturbations will be calculated in future work.
In order to evaluate the integral analytically, several mild simplifications are made. As previously mentioned, we use the fact that g(t) varies slowly during the production process and a reasonable estimate is to instead use its average value (12.5). Additionally, if we are restricting our attention to the study of sterile neutrinos with m s 1 MeV, then the first bracketed term inside of the logarithm (which is independent of τ) simplifies considerably.
The remaining τ dependence in the logarithm is a result of the Bose-Einstein suppression of the pions' thermal distribution and the 1/y 2 dependence is a result of the phase space factors (with m s 1 MeV).
With these simplifications and by expanding the logarithms in a power series the integral can be carried out analytically. The final result is given as where Γ(z, ν) is the upper incomplete gamma function.
To get the frozen distribution, we take the long time limit, τ → ∞, to arrive at n(τ, τ 0 , y) Equation (112) is the decoupled distribution function of sterile neutrinos with m s 1 MeV arising from pion decay near the QCD phase transition. This distribution function is valid for a wide range of sterile neutrino masses as we have only assumed m s /T(t) 1 for the period of production/freezeout (T f ∼ 10-15 MeV), which is valid as long as we consider m s 1 MeV.
The distribution function is shown for several values of τ 0 in Figure 2, where it can be seen that decreasing the lower limit increases the value of the distribution function. This is interpreted quite simply: production of steriles begins sooner and so the overall population is larger. If there are pions present in the plasma prior to the hadronization transition, then this could be extended back to temperatures until the finite temperature corrections to the pion decay constant are no longer reliable: τ ∼ m π / √ 6 f π ∼ 0.623.  Equation (112) is the decoupled distribution function of sterile neutrinos with m s 1 MeV arising from pion decay near the QCD phase transition. This distribution function is valid for a wide range of sterile neutrino masses as we have only assumed m s /T(t) 1 for the period of production/freezeout (T f ∼ 10-15 MeV), which is valid as long as we consider m s 1 MeV.

Bounds from Dark Matter and Dwarf Spheroidal Galaxies
The analysis above shows that freezeout occurs between τ ∼ 3-5 or T ∼ 10-15 MeV, so that sterile-like neutrinos are relativistic at the time of decoupling. For the light sterile neutrinos we consider here, we relate the number of relativistic species at the time of sterile decoupling to the photon temperature today by the usual relation between the plasma and photon temperatures: For m s 0.01 eV, we may neglect terms of O((y/x) 2 ), therefore the sterile neutrinos are non-relativistic today yielding γ,0 With this, the contribution to the density today is obtained using the distribution function, and we find where When m s 1 MeV the moments, I n (m s ) do not vary significantly and, for this mass range, they may be approximated by their value at m s = 0. We work under the assumption that m s 1 MeV so that we may use the limit I n (0) in subsequent calculations, and the limiting values are listed in Table 1. Using the results of Section 2, if we consider sterile masses with m ν s m l then we may neglect the sterile mass in both ∆, Λ so that we arrive at Considering light scalars simplifies the scales, Λ, so that the appropriate scales in the problem are Λ π→lν (m s = 0) ≡ Λ l ; Λ µ = 0.032 * |H µs | 2 10 −5 ; Λ e = 1.7 * 10 −6 * |H es | 2 10 −5 (118) so that Using the values of n γ h 2 /ρ c = 1/25.67 eV and Ω DM h 2 = 0.1199 while assuming g ν s = 2 and g d =ḡ = 12.5 leads to the bounds m s |H µs | 2 10 −5 ≤ 0.739 keV; m s |H es | 2 10 −5 ≤ 7242 keV .
As discussed in refs. [9,[74][75][76], the dark matter phase space density decreases over the course of galactic evolution and the primordial phase space density may be used as an upper bound to obtain limits on the mass of dark matter. Using observational values for dwarf spheroidal galaxies from ref. [76] a set of bounds complementary to those from CMB measurements can be obtained. As discussed, imposing the condition D p ≥ D nr gives us the constraint D p ≥ 1 Assuming, as before, that the sterile neutrino mass is much smaller than the charged lepton mass renders the phase space density independent of the sterile neutrino mass. This leads to a bound on the mass given as Using the results from Section (2), the phase space density is given as so that the bound becomes which can serve as a complementary bound to the limits set from Ω DM . Values of the phase space density today are summarized in ref. [76] and using the data from the most compact dark matter halos leads to bounds on sterile neutrino dark matter. The halo radius (r h ), velocity dispersion (σ), phase space density today and the calculated bounds are summarized in Table 2 where we chose several of the most compact dwarf spheroidal galaxies (a more thorough list is available in [76]).
The analysis of refs. [66,67] suggests that if sterile neutrinos are responsible for the x-ray signal, then production from π → µν s is a mechanism consistent with the data within a narrow region, while sterile neutrinos produced from π → eν are not.

Equation of State and Free Streaming
The equation of state for an arbitrary dark matter candidate is characterized by the parameter w(T) = P /ρ. A light sterile neutrino (m s 1 MeV) freezes out while it is still relativistic since m s /T 1 during production/freezeout; therefore, the results of the previous section hold. This distribution will then determine at what temperature this species becomes non relativistic, which is rewritten here explicitly in terms of m s /T: Many fermionic dark matter candidates that freeze out at temperature T f are treated as being in LTE in the early universe so that their distribution functions are given by the standard form To compare the new distribution to thermal results, assume that thermal particles with the same mass also freezeout while relativistic. The equation of state arising from thermal distributions and the non-thermal distribution we obtain are plotted as a function of m s /T in Figure 3. Note that the non-thermal distribution equation of state parameter is smaller for all times. This is a reflection of the enhancement of small momentum so that the non-thermal distribution results in a dark matter species, which is colder and becomes non-relativistic much earlier than the thermal result. In summary, the thermal distribution produces particles that become non-relativistic when m s /T 1, whereas the pion decay mechanism produces particles that become non-relativistic when m s /T ∼ 1. This non-thermal distribution function produces a dark matter candidate that is colder than those produced at LTE.  The free streaming wavevector enters when one considers a linearized collisionless Bolzmann-Vlasov equation describing the evolution of gravitational perturbations, which ultimately lead towards structure formation [74][75][76][77][78]. The free streaming wavevector k f s leads to a cutoff in the linear power spectrum of density perturbations and is given by Modes with k < k f s lead to gravitation collapse in a manner akin to the Jeans instability. This is shown explicitly and discussed at length in refs. [74][75][76]. Assuming that a light sterile neutrino is the only dark matter (so that ρ ν s = ρ DM ) and using the results of Section 2 (for a non-relativistic species), the free streaming wave vector is given by 2 dy y 2 f d (y) dy y 4 f d (y) .
Using the latest values from Planck [29] sets the free streaming wave vector as where we have used the notation λ l f s (0) ≡ λ f s (0) π→lν . These values of the free-streaming length, namely the cutoff in the (DM) power spectrum, are consistent with the size of cores in dwarf galaxies.
Our focus in this study is to bring to the fore alternative production mechanisms of sterile neutrinos and to obtain their frozen distribution function from such processes, and from these their clustering properties for galaxy formation. A further assessment on whether sterile neutrinos produced by a particular mechanism are, indeed, suitable dark matter candidates would require assessing whether their mass range is consistent with the various bounds [122,123]

Conclusions
Sterile neutrinos are SU(2) singlets that do not couple directly to other degrees of freedom, but mix with active neutrinos via a mass matrix. A sterile-like mass eigenstate with a mass scale few keV could be a suitable warm dark matter candidate that could explain various discrepancies between N-body simulations of structure formation with (CDM) and observations of dwarf spheroidal galaxies. In this brief review, we have described alternative production mechanisms of sterile-like neutrinos beyond and within the standard model of particle physics. The main aspects discussed in this review are the following: • We began with a description of the dynamical properties of dark matter particles described by frozen distribution functions, such as effective equation of state, phase space density and the free-streaming length, which determines the cutoff in the dark matter power spectrum. • We considered a simple extension beyond the standard model in which sterile-like neutrinos couple to a Higgs-like scalar via a Yukawa coupling, obtaining the quantum kinetic equations for production, the distribution function at freeze out and from it the clustering properties of this (DM) candidate such as the abundance and free streaming length. • Within the standard model, we have argued that in the basis of mass eigenstates, standard model processes that produce active neutrinos also produce sterile-like neutrinos; however, these are suppressed by the small mixing angles. We discussed several possible production channels, obtained the general form of the quantum kinetic equations to leading order in the mixing angles and analyzed the possibility of thermalization. This analysis leads us to conclude that the final distribution function of sterile-like neutrinos is a result of the various kinematically allowed production processes available during the cosmological history of the early universe. • As specific example, we focused on the production of sterile-like mass eigenstates from the decay of pions after the QCD phase transition, including finite temperature effects. Shortly after the QCD transition at T QCD 155 MeV, t 10 µs pions being the lightest degrees of freedom in QCD feature large populations and their main decay channels (for charged pions) produce neutrinos. In the basis of mass eigenstates, these processes produce sterile-like neutrinos. We obtained and solved the quantum kinetic equations for production and freeze-out, and from the solution we obtained the clustering properties of sterile-like neutrinos produced via this mechanism. We find that these could be a suitable warm dark matter candidate with the correct abundance if the mass is of the order of few keV and their free streaming length of a few kpc consistent with the scale of cores in dwarf spheroidal galaxies.
Although direct observational evidence for keV sterile-like neutrinos remains elusive, these still remain, at least theoretically, as suitable warm dark matter candidates with several possible production channels beyond and within the standard model in the early universe.