Lyman-break galaxies [433] at $z \gtrsim 3$; the excess of star-forming low-mass galaxies at $z \gtrsim 3$ reflects into an excess of red low-mass galaxies at $z \approx 0$ [427, 434].

(b) The number of galaxies with stellar mass $M_* \lesssim 10^{10} M_{\odot}$ in the Universe is systematically overpredicted (by a similar factor) by all theoretical CDM models [435–438] in the whole redshift range $1 \lesssim z \lesssim 3$ [439].

(c) The local distribution of galaxy rotation velocities (velocity function) predicted by CDM models is characterized by an excess of galaxies in the low-velocity end [440]. Since the rotation velocity directly probes the depth of the potential wells rather than the their gas or stellar content, the velocity distribution is less prone to the effect of energy feedback (which enters only in selection effects or in determining the shape of the velocity curve), this is considered as a challenging issue for solutions based on feedback effects.

As for satellite galaxies, long standing problems of CDM scenario (or, at least, points that require the improvement in the implementation of the different baryonic processes) include:

(d) The mass and luminosity distributions of satellite galaxies predicted by simulations and by semi-analytic models tend to exceed the observed abundances [46, 339, 441]. The improved implementation of processes involving baryons — such as heating or winds caused by Supernovae explosions or by the UV background — may be effective in suppressing star formation (SF) in the low-mass haloes of satellite galaxies at low redshifts $z \le 1$ [50, 417, 442–445] thus reducing the discrepancy. However, recent observational results for larger samples of central galaxies are still in tension with the predictions of CDM models [446]. In addition, reducing the number of detectable satellites through feedback processes results in hosts that are too massive to be accommodated within measured galactic rotation curves [447], and it still does not avoid over-producing low-mass galaxies at redshifts $z \gtrsim 1$.

(e) CDM models predict most (in fact, more than 80% [448]) of satellite galaxies to be quiescent, with low specific star formation rates, $\text{SSFR} \le 10^{-11} \text{ yr}^{-1}$, while observations indicate a much lower fraction ~30% [449–453] for satellite galaxies with stellar mass in the range $10^8 M_{\odot} \le M_* \le 10^{10} M_{\odot}$, as shown in figure 10 (right). The discrepancy is generally ascribed to the implementation of the various baryonic processes contributing to the quenching of star formation in satellites, like the stripping of hot gas [430, 454–457] due to the pressure of the gas in the host halo. However the persistence of such a problem in the most recent semi-analytic models [458–460] and N-body hydrodynamical simulations [460–462] indicates that it constitutes a major challenge for all CDM galaxy formation models, and a fundamental problem in understanding the evolution of low-mass galaxies [460]. A solution has been recently sought in the implementation of star formation [463], in reincorporation time-scales of the hot gas ejected by SNe or winds [464], or in a delayed and gradual stripping of hot gas from satellites [438, 465, 466], but none of the proposals proved to be completely successful. This is because a gentle mode of strangulation should be at work in order to support cooling and star formation in satellite galaxies for several Gyr. Phenomenological estimates of the quenching time scales required to comply with the observed quiescent fraction of satellites yield values in the range 3–6 Gyr [450, 467–469], and is not easy to obtain such a delayed and gentle quenching time scale from a physical modelling of gas stripping.
On the other hand, ref. [470] has shown that the long quenching time scales may be understood in terms of a starvation scenario. However, this paper mainly deals with the satellites in the local group and it is meant to explain the small values of the quiescent fraction in such a group, while the problem is that — when all the host haloes are considered — the CDM quiescent fraction is instead too high. Still they show that starvation can actually be relevant for satellites with masses larger than $10^8 M_\odot$, while for smaller galaxies the authors suggest ram-pressure stripping to be more efficient in getting — however — short quenching times (i.e., still below 2 Gyr).

While a solution in terms of refined treatment of baryonic processes is still possible, several authors have considered the possibility that the problem is rooted in the CDM power spectrum, and explored the effects of suppressed small-scale perturbations. Alternative DM models such as warm and mixed DM would naturally lead to a suppression of the power spectrum at small scales, resulting in a reduced dwarf galaxy abundances [471, 472]. This could be in better agreement with observations [335, 337, 473], as shown in figure 10 (left panels) for the case of a WDM spectrum corresponding to a thermal relic particle $m_X = 1$ keV. Assuming a WDM spectrum also yields a quiescent fraction of satellite dwarf galaxies in agreement with observations. In fact, the delay in the average collapse time of the first structures in WDM scenarios results into a prolonged star formation history in small mass objects, thus yielding a larger fraction of star-forming low-mass satellites at low redshifts compared to the CDM case [474], as shown in the right panel of figure 10.<fnref n="22" />

The improvement in the dwarf and satellite abundances and passive fraction compared to the CDM case illustrated in the figure requires a WDM suppression of the CDM power spectrum corresponding to thermal relics $1 \lesssim m_X \lesssim 3$ keV [334, 473, 474, 479]. Indeed, a thermal relic mass $m_X = 1.5$ keV, corresponds to that produced by (non-thermal) sterile neutrinos with mass $m_X \approx 6-7$ keV (depending on the production mechanism). These constitute the simplest candidates, see e.g. [480] for a Dark Matter interpretation of the origin of the recent unidentified X-ray line reported in stacked observations of X-ray clusters with the XMM-Newton X-ray Space telescope with both CCD instruments aboard the telescope, and the Perseus cluster with the Chandra X-ray Space Telescope [190, 191]. Such a value is consistent with lower limits ($m_X \ge 1$ keV for thermal relics at $2\sigma$ level [481]) derived from the density of high-redshift galaxies set by the two objects already detected at $z \approx 10$ by the Cluster Lensing And Supernova survey with Hubble (CLASH). However the thermal relic equivalent mass of sterile neutrinos inferred from X-ray observations is only marginally consistent with the limits set by different authors ($m_X \ge 2.3$ keV [482]; $m_X \ge 1.5$, [341]; $m_X \ge 1.8$, [313]) from the abundance of ultra-faint Milky Way satellites measured in the Sloan Digital Sky Survey, see e.g. [483], and definitely in tension with the constraint derived by comparing small scale structure in the Lyman-$\alpha$ forest of high-resolution ($z > 4$) quasar spectra with hydrodynamical N-body simulations which yield bounds ranging from $m_X \gtrsim 1.8$ to $m_X \gtrsim 3.3$ keV [327, 329, 331, 332]. Note however that various uncertainties may still affect the constraints (see discussions in [484-486]): the comparison of subhaloes to Milky Way dwarfs assumes a factor $\approx 4$ correction for the number of dwarfs being missed by current surveys, and lower correction factors would appreciably weaken the constraints; Lyman-$\alpha$ is also a challenging tool, as it requires disentangling the effects of pressure support and thermal

<footnote n="22">²²Note that the result on the right panel has not been published, but the quiescent fraction for this case was instead computed specifically for this White Paper. However, the model adopted is exactly that used to compute the results in ref. [474]. All details on the model can be found therein.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** *Top Left Panel:* the evolution of the K-band luminosity functions in the WDM cosmology (solid line) is compared with data from [432]. Dashed line refer to the standard CDM case. *Bottom Left Panel:* the evolution of the stellar mass function in the WDM cosmology (solid line) is compared with the standard CDM case (dashed); data are from [475] (squares), [435], and [439]. In both the upper and the lower left panels, the arrows show the range of magnitudes and stellar masses corresponding to the free streaming mass. *Right Panel:* satellite quiescent fraction at different stellar masses. Data are from [450], red triangles; [451], green squares; ref. [449] and [453], blue dots; [448], yellow crosses. These are compared with predictions from a WDM model [474] and from different CDM models in the literature: the MORGANA model [476], grey dot-dashed line]; the Millenium-based SAM [477], blue short-dashed line]; the [478] code (purple dashed line); the Millenium-based model [438], petrol solid line (taken from [451]); the Rome model [474], solid line.</figcaption></figure>

broadening from those caused by DM spectrum, as well as assumptions on the thermal history of the intergalactic medium and of the ionizing background.

## 4.3 Inner density profiles of small galaxies and the cusp-core problem (author: W. Evans)

Cosmological ΛCDM simulations provide galaxy-scale structures that are too compact. One manifestation of this is that the inner density profiles of simulated haloes follow a power-law cusp $\rho \sim r^{-\alpha}$ with $\alpha \approx 1$. This result was first obtained by [418, henceforth NFW], and it has survived the mammoth increases in computing power and resolution over the two succeeding decades, see e.g. [339, 487]. Cusps are a fundamental feature of dissipationless simulations, arising from the intrinsic coldness of the velocity distributions.

As this is an unambiguous prediction of ΛCDM, much effort has been expended on testing it. Large galaxies like the Milky Way are dominated by baryons in the center [488, 489]. The uncertain effects of the growth of baryonic disks, bars, and bulges may have altered and homogenized the Dark Matter distribution, possibly erasing any cusps. So, the best testing ground are dwarf and low surface brightness galaxies, as here any effects of baryonic feedback are smaller.

Dwarf spirals and low surface brightness galaxies are usually endowed with abundant neutral hydrogen. Here, the gas rotation curve can be directly linked to the matter distribution, dominated by Dark Matter. Possible problems may include non-circular motions, the ambiguous position of the center and the resolution of the data [490]. However, high
resolution Hα rotation curves, where available, show excellent consistency with gas rotation curves [491], suggesting beam smearing is not a significant obstacle. Recent surveys like THINGS (“The HI Nearby Galaxies Survey”) have acquired very high spectral and spatial resolution data from the VLA (Very Large Array) on 34 galaxies. The quality of this data has overcome most of the problems and led to robust conclusions. For example, [492] studied the dIrr galaxies IC 2574 and NGC 2366, correcting for non-circular motions, and concluded that both are best described by dark halo models with a kpc sized density core. Ref. [493] looked at seven dwarf galaxies in the THINGS sample and showed that the rotation curves are gently rising in the inner parts, consistent with core-like models with $\alpha = -0.29 \pm 0.07$, significantly different from the $\alpha = -1$ value from dissipationless simulations. This situation seems to persist in low surface brightness (LSB) galaxies as well. This is a particularly important regime to study as the large extent and low stellar density of LSB galaxies does not readily permit the Dark Matter to be re-distributed via feedback from star formation. Ref. [494] presented 30 high resolution HI rotation curves for LSB galaxies. They fitted the data with both cored pseudo-isothermal and cusped models, and found the former to be preferred. Even for the cusped fits, the distribution of concentrations has too low a mean to match the predictions of simulations.

Dwarf spheroidal (dSph) galaxies have also been studied intensively. These galaxies are highly DM-dominated, but have little or no HI gas, and so the mass profile must be deduced from stellar dynamical modelling. This procedure is susceptible to well-known difficulties, such as the mass-anisotropy degeneracy [495]. However, dSphs have the advantage that they are some of the closest galaxies to the Milky Way, and so samples of radial velocities of bright giant stars can be easily gathered. The size of such samples has increased from hundreds to thousands in recent years [496, 497]. It was soon realized that, for a single stellar population, only the mass within the half-light radius is well-constrained [498], and the data cannot discriminate between cores and cusps [499]. However, many dSphs are built from multiple populations of stars, typically at least metal-rich and metal-poor ones. If the surface brightness and line-of-sight kinematics of two or more populations can be measured, then the degeneracy in the data is broken. Ref. [500] used statistical techniques to separate the metal-rich and metal-poor populations in the Fornax and Sculptor dSphs. Using a simple estimator for the mass enclosed within the half-light radius of each population, they demonstrated that the central density slopes were shallower than NFW and consistent with cores of size several hundred kpc. Ref. [501] used full phase space distribution functions to model the multiple populations in Sculptor and demonstrated that cored halo models are very strongly preferred over cusped ones. Subsequently, ref. [502] found a way of recasting the arguments in terms of the virial theorem, which has the benefit of robustness and generality. Ref. [503] applied the virial method to Fornax and found a core size of $1_{-0.4}^{+0.8}$ kpc. At least for two of the largest dSphs around the Milky Way, therefore, the data point unambiguously towards cores. Although [504] did construct fits to Sculptor using NFW haloes, this work remains open to the objection that the data are not well matched in the center, where the velocity dispersions of the two populations are too high (a typical failing of all cusped models).

Given the unambiguous results from observations, a number of investigators have used simulations to explore whether baryonic feedback could erase a cusp and generate a core. Early calculations [505], assuming a single burst of star formation which generates a baryonic wind, showed that the energetic requirements are formidable. However, simulations of multiple bursts of star-formation demonstrated that cusps could be smoothed away [506]. For example, ref. [507] presented SPH and N-body simulations including star formation and
supernova driven outflows. They showed that repeated gas outflows and bursty star formation can cause the central density of the Dark Matter to become cored. The rapidly changing potential in simulations, caused by outflows of expanding gas from bursts of star formation, can transfer energy into the collisionless Dark Matter particles and generate a Dark Matter core [508]. It now seems clear that this process can operate in simulations, but it is not yet evident that the simulations do accurately reproduce the difficult physics of star formation and supernova driven winds. Feedback problems are most acute for LSB galaxies, as star formation is minimal and inefficient when gas density is low. LSBs are typically metal-poor and dust-free, both implying that they form comparatively few stars over a Hubble time. Additionally, supernova winds cannot expel gas from the high mass end of the LSB galaxy sequence. It is therefore a substantial challenge to use feedback to erase cusps in all LSB galaxies, even if the process can be made to work for dwarfs [509].

WDM (such as sterile neutrinos) does naturally produce cores. The Dark Matter phase space density is finite if the Dark Matter temperature is finite. If the Dark Matter is collisionless, then Liouville's theorem guarantees that the phase space density cannot increase. So, haloes are expected to have cores if the Dark Matter is not cold [75], though the cores do not appear to be large enough to satisfy the data on LSB galaxies [335, 342, 510]. WDM also delays the collapse of dwarf-sized haloes and their associated star formation to later epochs. It remains unclear whether WDM is more successful than CDM in reproducing the intrinsic properties of dwarf galaxies.

## 4.4 Too-big-to-fail (author: E. Papastergis)

### 4.4.1 Introduction and background

The too-big-to-fail (TBTF) problem is a small-scale cosmological challenge to the ΛCDM paradigm. It posits the issue that it may not be possible to simultaneously reproduce the observed number density and internal kinematics of dwarf galaxies within the ΛCDM model. The TBTF problem is closely related to the “missing satellites” and “cusp-core” problems, since it is formulated in terms of both the abundance of galaxies as well as their kinematic properties (see sections 4.2 and 4.3). At the same time, the TBTF challenge represents a stronger theoretical issue for ΛCDM than either of the two aforementioned problems on their own. This is why the TBTF problem has become one of the most active topics of research in extragalactic astronomy today.

The TBTF problem was first identified in the population of Milky Way (MW) satellites [53, 54]. In particular, it is natural to expect that the bright (or “classical”) MW satellites would be hosted by the most massive subhaloes of the MW. Boylan-Kolchin et al. used the Aquarius [511] and Via Lactea II [512] simulations to show that the largest subhaloes of a MW-sized halo are “too massive”<fnref n="23" /> to accommodate the observed kinematics of the bright MW satellites (figure 11). Following the initial formulation of the TBTF problem in the context of the MW, several non-cosmological solutions were quickly identified. For example, the TBTF problem would likely not occur if the virial mass of the MW was relatively low ($M_{vir,MW} \lesssim 1 \times 10^{12} M_{\odot}$), or if the MW was a statistical outlier [513–515].

However, the MW does not seem to be a unique case. Ref. [516] showed that the same TBTF problem is encountered in the satellite population of the Andromeda galaxy

<footnote n="23">²³To be more precise, the massive subhaloes in the Aquarius and VL II simulations were found to be “too dense” to accommodate the kinematics of the bright MW satellites. That is, the dynamical mass enclosed by the half-light radii of the MW satellites is too low to be compatible with the mass profiles of the largest subhaloes in a MW-sized system.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** The solid curves represent the rotation curves of massive subhaloes of a MW-sized halo in the Aquarius DM simulation. The datapoints correspond to measurements of the rotational velocity of bright MW satellites at the half-light radius, $V_{\text{rot}}(r_{1/2})$. Note how most of the expected hosts of the bright MW satellites are “too massive” to be accomodated within the measured kinematics of the satellites. (Reproduced with permission from [54].)</figcaption></figure>

(M31). In addition, [517] and [518] have argued based on SDSS data that the MW’s satellite population is typical for a galaxy of its size. Overall, it seems that a comprehensive solution of the TBTF problem should be applicable to the satellite population of any MW-like galaxy. Such a solution has recently been put forward by [55, 519, 520]. In particular, they showed that subhaloes in hydrodynamic simulations are more vulnerable to tidal stripping than their DM-only counterparts. Efficient stripping can lower the DM density in the inner parts of massive subhaloes compared to the DM-only case, resulting in kinematics that agree with the measurements for MW satellites, see figure 3 in [55].

This baryonic solution to the satellite TBTF problem has been well received by the community, because it is physically well-motivated and broadly applicable. At the same time, it relies on satellite-specific environmental processes that are not expected to operate on isolated galaxies. As a result, it is of key scientific importance to test whether the TBTF problem is encountered also for dwarf galaxies in the field. To this end, ref. [521] measured the stellar kinematics of a sample of seven non-satellite dwarfs in the Local Group (LG; within ~1.5 Mpc from the MW). They found that their kinematics do not differ from the kinematics of MW satellites of the same luminosity. This constitutes indirect evidence against the Zolotov & Brooks mechanism, which requires the massive subhaloes of the MW to be heavily stripped. Furthermore, [522] compared the stellar kinematics of the Kirby et al. [521] dwarfs with the kinematics of the largest haloes present in a simulated LG analog (ELVIS simulation [523]). They found once again that the expected hosts of the LG dwarfs were to massive to be compatible with the measured dwarf kinematics, see figure 7 in [522]. Their finding demonstrated that the TBTF problem is relevant for both satellites and non-satellite galaxies in the LG.

Even earlier, [524] obtained compelling evidence that the TBTF problem is also present for dwarf galaxies in the field. Ferrero et al. [524] assigned $\Lambda$CDM haloes to galaxies such that the number density of galaxies in SDSS as a function of their stellar mass is reproduced.
<figure><img src="image_13.png" /></figure>

**Figure 12.** The upper gray band represents the typical range for the rotation curves of haloes with $M_{vir} = 10^{10} M_{\odot}$. These haloes are the least massive haloes that are expected to host detectable galaxies in a $\Lambda$CDM universe, according to AM. The colored points represent the rotational velocity of a sample of dwarf galaxies, measured at the outermost point of their HI rotation curve ($V_{rot}(R_{out})$). Similarly to what seen in figure 11, the expected $\Lambda$CDM hosts of dwarf galaxies are “too massive” to be compatible with the galaxies’ observed kinematics. Unlike figure 11, most of the dwarfs plotted here are fairly isolated galaxies in the field. (Reproduced with permission from [524].)

A quantitative galaxy-halo connection under this constraint can be derived statistically, using the technique of abundance matching [AM; 525]. In a $\Lambda$CDM universe, AM predicts that all detectable field galaxies ($M_* \gtrsim 10^6 - 10^7 M_{\odot}$) should be hosted by relatively massive haloes ($V_{h,max} \gtrsim 45 \text{ km s}^{-1}$). This result is a straightforward consequence of the fact that the mass function of haloes in $\Lambda$CDM rises much faster than the stellar mass function of galaxies. In simpler terms, $\Lambda$CDM predicts that the abundance of haloes rises quickly with decreasing halo mass. Galaxy formation must therefore be restricted to relatively massive haloes with $V_{h,max} \gtrsim 45 \text{ km s}^{-1}$, because lower mass haloes are too numerous compared to the number of observed dwarf galaxies.

They then set out to test this theoretical expectation, by compiling a literature sample of dwarf galaxies with measured rotation curves in the 21cm emission line of atomic hydrogen (HI). Surprisingly, they found that the HI rotation curves of a large fraction of field dwarfs suggested that their host haloes have masses well below the expected “threshold” for galaxy formation (figure 12). The result of Ferrero et al. was subsequently confirmed and strengthened by the work of [447]. They combined an up-to-date literature sample of galaxies with HI rotation curves and a census of low-mass galaxies from the ALFALFA 21cm survey [526], to convincingly show that the TBTF problem is also present in the field [447, figures 6 & B1]. Unlike previous analyses, the work of Papastergis et al. [447] is based exclusively on HI data, and so it can circumvent many observational uncertainties related to optical measurements present in previous works.
### 4.4.2 Possible solutions of the TBTF problem within ΛCDM

There is now extensive evidence that the TBTF problem concerns not just satellites, but low-mass dwarf galaxies in general. Accordingly, a generic solution to the problem must be applicable to both the satellite and the field context. The most plausible solutions to the problem within ΛCDM are “baryonic”, in the sense that they invoke modifications to the abundance and internal structure of DM haloes due to baryonic processes:

(a) *Baryonic effects on the abundance of haloes:* there are two baryonic effects that can alter the number density of haloes compared to the DM-only case. First, haloes with $V_{h,max} \lesssim 70 \text{ km s}^{-1}$ can lose the majority of their baryonic content to galactic outflows powered by supernova feedback [527]. This effect is not captured by DM-only simulations, because the latter treat the total matter density of the Universe ($\Omega_{\text{mat}} = \Omega_{\text{DM}} + \Omega_{\text{B}} = 0.32$) as a dissipationless fluid. As a result, haloes in this mass range are expected to have about 20% lower mass than their DM-only counterparts. Second, haloes with $V_{h,max} \lesssim 25 \text{ km s}^{-1}$ are subject to cosmic reionization feedback. These lower mass haloes are unable to accrete gas after the Universe became reionized ($z \gtrsim 6-8$), and are therefore expected to be “dark” in most cases. The combined consequence of these two effects is to alter the abundance of detectable low-mass haloes compared to the result of a DM-only simulation, see figure 2 in [420]. The baryonic modification of the abundance of haloes can alleviate the TBTF problem, because the host haloes assigned through AM to dwarf galaxies are less massive than in the case of a DM-only simulation, see figure 3 in [420].

(b) *Baryonic effects on the rotation curves of haloes:* hydrodynamic simulations have shown that the repeated blow-out of gas from the central regions of dwarf galaxies can create shallow “cores” in the DM profile of their host haloes [508, 528]. This effect can alleviate the TBTF problem, because the rotation velocity in the inner parts of a cored halo is lower compared to its “cuspy” DM-only counterpart of the same total mass. Core creation through baryonic feedback has two properties that are important in the context of the TBTF problem: (i) The efficiency of core creation has a characteristic dependence on galaxy mass [529]. In particular, core formation is most efficient in dwarf galaxies with $V_{\text{rot}} \sim 50 \text{ km s}^{-1}$, but the efficiency drops for higher and lower mass galaxies [see 530, figure 6]. (ii) Core creation is only effective over the region where the bulk of star formation takes place, typically $R \lesssim 1 \text{ kpc}$ for dwarf galaxies in the field [507]. This means that the rotational velocity of a dwarf galaxy can be significantly lower than the DM-only value when measured at a small galactocentric radius, but the difference becomes small at radii $R \gtrsim 2 \text{ kpc}$.

Overall, it is still not entirely clear whether the baryonic effects described above can fully resolve the TBTF problem. For example, ref. [420] shows that the effect (a) can explain the observed kinematics of faint MW and M31 satellites, since most low-mass haloes with $V_{h,max} \lesssim 25 \text{ km s}^{-1}$ are expected to be dark (see their figure 4). However, this mechanism alone does not seem to be able to resolve the field TBTF problem, which is encountered at a larger mass scale ($V_{h,max} \gtrsim 35 - 40 \text{ km s}^{-1}$; see [447], figures 9 & B1). On the other hand, [531] have argued that effect (b) can be used to explain the kinematics of low-mass dwarfs, both in the case of satellites and isolated objects in the LG (their figure 2). At the same time, ref. [447] has used the results of a set of hydrodynamic simulations with efficient feedback [55, 507, 532] to argue that core creation cannot fully resolve the TBTF problem for field dwarfs (see their figure 10).
<figure><img src="image_14.png" /></figure>

**Figure 13.** Panel (a): the cumulative velocity function of galaxies measured by the ALFALFA survey (green dotted line; [447]), compared to the cumulative velocity function of haloes in four cosmological models: CDM and WDM with $m_{\text{WDM}}=1,2$, & 4 keV (blue, red, orange, and yellow lines respectively; ref. [334]). Panel (b)–(e): the solid line in each panel represents the AM relation between the measured HI rotational velocity of a galaxy, $V_{\text{rot,HI}}$, and the maximum rotational velocity of the host halo, $V_{\text{h,max}}$. The gray datapoints represent a sample of field dwarfs with resolved HI kinematics. Their y-axis position is determined by the width of their HI emission line, while their x-axis position represents the most massive halo that can be accommodated within their HI rotation curve. Note how the upper limits on $V_{\text{h,max}}$ set by the kinematics of dwarfs shift upwards as $m_{\text{WDM}}$ decreases, due to effect (b). Cosmological models where the AM relation (solid line) is inconsistent with the upper limits set by individual galaxies are expected to overestimate the number density of dwarf galaxies detected in the field. (Reproduced with permission from [447].)

### 4.4.3 Can Warm Dark Matter solve the TBTF problem?

Warm Dark Matter (WDM) has been regarded as a promising solution to the small-scale challenges of ΛCDM. The reason is that small-scale power in suppressed in a WDM universe, and as a consequence, structure formation at the mass scales of dwarf galaxies is curbed [334, 341, 482, 533, 534]. In the context of the TBTF problem, WDM has two advantageous features:

(1) **There are fewer low-mass haloes in a WDM universe.** The mass (or velocity) function of haloes in a WDM universe is suppressed with respect to a CDM universe. As far as the TBTF problem is concerned, this means that observed galaxies can be assigned to lower mass haloes than in CDM, without overestimating the galaxies’ number densities.

(2) **WDM haloes have low concentrations.** At fixed total mass, a WDM halo has a lower concentration than its CDM counterpart [e.g., 337, 535]. As a result, the rotation curve of a WDM halo rises more slowly than in CDM, lowering the rotational velocity measured at small galactocentric radii.

Figure 13 provides a graphic representation of how effects (1) and (2) can alleviate the field TBTF problem. Panels b-e show that both effects become more pronounced as the mass$^{24}$ of

<footnote>$^{24}$For all WDM particle masses quoted in this section, $m_{\text{WDM}}$, refer to thermal relic particles.</footnote>
the WDM particle decreases. In fact, for WDM models with relatively large particle masses, $m_{\text{WDM}} \ge 4$ keV, there is no difference between the CDM and WDM case as far as the TBTF problem is concerned. Even in a 2 keV WDM model the TBTF problem is still present, even though somewhat less severe. According to figure 13, it takes a *thermal relic* WDM model with a particle mass as low as $\approx 1$ keV (and certainly below 2 keV [536]) to fully resolve the issue. The resonantly produce sterile neutrino dark matter models, on the other hand, can achieve good agreement with the observations for large range of masses [537].

Figure 13 therefore illustrates the main challenge that WDM faces as a successful solution to the TBTF problem: *the required particle mass is quite low, $m_{\text{WDM}} \lesssim 2$ keV. It is unclear whether WDM models with such low particle masses are viable in view of the constraints posed by other astrophysical probes. For example, current constraints from measurements of the small-scale power spectrum of the Ly-$\alpha$ forest [538] put a lower bound at 2.1 keV if one does not impose any bias on the temperature of the intergalactic medium at redshifts $z = 5-6$. However, this bound becomes much stronger, $m_{\text{WDM}} > 3.3$ keV at 95% confidence, if one assumes powerlaw scaling of temperature with redshift, [331].$^{25}$ The suppression of small-scale power in a WDM model with $m_{\text{WDM}} \lesssim 2$ keV may be “too much of a good thing”. More specifically, for sufficiently low values of $m_{\text{WDM}}$ a MW-sized halo may not possess enough subhaloes to host the number of known ultrafaint satellites in the MW. Based on the consideration above, the authors of [482] derive a lower limit of $m_{\text{WDM}} > 2.3$ keV, while [341] obtains a value of $m_{\text{WDM}} > 1.5$ keV; see also ref. [313] for an alternative limit from satellite counts. The difference in results is partially due to different realizations of the MW-sized halo in simulations (as discussed e.g. in [535]). Furthermore, [334] and [539] have argued that the same is true in the field; WDM models with $m_{\text{WDM}} < 2$ keV severely underestimate the number density of the smallest field dwarfs with $V_{\text{rot}} \approx 10-20$ km s$^{-1}$ (their figures 2 and 6, respectively). On the other hand, non-thermal DM settings could provide a game-change, see e.g. [535, 537].

## 4.5 The kinematics and formation of subhaloes in Warm Dark Matter simulations (authors: M. Wang, L. Strigari)

There is debate as to whether the central densities of DM-dominated galaxies are less dense than the haloes that are predicted in DM-only simulations [500, 504]. This issue along with the lack of observed satellites around the Milky Way can be ameliorated by considering alternative Dark Matter models. One particularly intriguing scenario is Warm Dark Matter (WDM), which differs from CDM in that the Dark Matter has a non-negligible velocity dispersion. Generally, these WDM models come with two distinct observational signatures. First, there is a reduction in clustering and a delay in collapse times for structures on linear scales approximately smaller than the free-streaming scale of the WDM particle. Second, the velocity dispersion of WDM imposes an upper limit on the Dark Matter phase-space density.

Numerical simulations are now rapidly improving, allowing for more precise predictions for Dark Matter halo formation in a WDM universe. Understanding the predictions of these simulations, along with other alternative to $\Lambda$CDM models, might provide clues to discriminating between different particle Dark Matter models. For example, it is found that WDM models with thermal relic mass $\sim$ keV can generate a truncation on matter power spectrum at scale $\sim$ a few Mpc. This results in a suppression of the formation of small structure below the free-streaming scale of WDM and thus result in delayed formation of haloes. As

<footnote>$^{25}$Exclusion limits are quoted in this article as published in the original reference, and may be subject to limitations (observational or theoretical) of the specific method used to derive them.</footnote>
<figure><img src="image_15.png" /><figcaption>Figure 14. Histograms of subhalo properties for subhaloes that fit the kinematics of the Fornax dSph [542]. Here we show three properties of subhaloes at the present time: the maximum circular velocity ($V_{max}$, left panels), and the subhalo mass ($M_{sub}$, right panels). The navy histograms show the results from ΛCDM, the light blue histograms are for the WDM 2.3 keV simulation, the green histograms are for the WDM 1.6 keV simulation.</figcaption></figure>

a result these simulations generally indicate that WDM haloes are less concentrated on account of their typically later formation epochs [58, 341]. These models also imply a different reionization history [486] or star formation history at early times [540]. At present time at subgalactic scales, the stellar abundance and metallicity can thus be quite different from the CDM predictions [541].

WDM simulations also predict different properties for Galactic substructure than in ΛCDM simulations. It is thus pressing to understand how these predictions can be translated into tests of observations of Milky Way substructure in the form of dwarf satellite galaxies. Here we briefly address this question from the perspective of the kinematics and formation of subhaloes in WDM scenarios using N-body numerical simulations.

To address this question, in [542] we map subhaloes in WDM simulations to dwarf spiralodal (dSph) galaxies, using kinematic observations of the dSphs. A similar approach was taken in [543] for ΛCDM simulations. With the simple assumptions that these are spherical dynamically equilibrated systems, there is a strong degeneracy between the statistics of stellar orbits (i.e. whether velocity dispersions are isotropic, or are radially or tangentially biased) and the shape of the stellar and Dark Matter density profiles. Additional observational constraints like the measurement of stellar proper motion will be necessary to break these degeneracies. Nevertheless, we find that in WDM cosmology the subhaloes that host dSphs satellite galaxies exhibit different properties due to their shallower potential and late formation history.

In figure 14 we show the subhalo properties that provide good-fits to the Fornax galaxy kinematic data [542]. Here we show the maximum circular velocity $V_{max}$ at present time, and mass of the subhalo $M_{sub}$. We find that the distribution of these parameters from WDM simulations can be quite different than those from ΛCDM simulations.

We also find that the mass distributions of the WDM haloes that host the dSphs are different from those that host the dSphs in ΛCDM simulations. The results are that WDM models tend to reside in more massive subhaloes than in ΛCDM because its less concentrated profiles. In the case of Fornax, which is the brightest Milky Way dSph with high-quality stellar kinematic data, in ΛCDM the $V_{max}$ is very well-constrained between ~16–20 km/s with average about 17 km/s. However, for WDM models, the $V_{max}$ can be as large as ~ 28 km/s.
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Mass assembly history for the subhaloes that fit the kinematics of the Fornax dSph [542]. The gray lines show the subhaloes that are less luminosity than Fornax using abundance matching technique from [545], and the blue lines show those with luminosity that match Fornax. The vertical dotted lines show the infall times of the blue curves. In the left panel the dashed black lines include 68% of the simulation fits, and the black solid lines show the mean value.</figcaption></figure>

It is also interesting to determine the formation histories of subhaloes in WDM and ΛCDM. In figure 15 we show the mass assembly history for the Dark Matter subhaloes that match the kinematics of the Fornax dSph in both ΛCDM and WDM with 2.3 keV mass [542]. The blue lines show those subhaloes that match the luminosity of Fornax using abundance matching technique [544]. Interestingly we see that the WDM subhaloes are usually formed or gain significant amount of mass much later than the ΛCDM subhaloes.

These same results also can provide an indication of how rare a Fornax-like dSph is in both ΛCDM and WDM. In the ΛCDM simulations we found that only two out of 124 subhaloes that match the kinematics of Fornax also predict the luminosity that is observed. If the prediction of abundance matching holds, the bright luminosity of Fornax requires it to be hosted by a progenitor subhalo that is more massive than its present mass. Since figure 14 shows that the present-time mass of ΛCDM Fornax candidates are typically less massive than the WDM candidates, this suggests that Fornax subhaloes in the ΛCDM might suffer from severe tidal stripping and significant loss on Dark Matter halo mass. It is not clear whether this implies visible stellar tidal streams or not, because the stellar distribution is deeply embedded at the center of Dark Matter potential. However, for WDM the amount of Dark Matter stripping may be less significant because the present-time mass is larger than it is in ΛCDM. Currently, stellar tidal streams have not been observed in Fornax, and from the orbital properties indicate Fornax has not yet gone through many pericenter passages.

The kinematical properties of dSphs may thus provide unique clues to the nature of particle Dark Matter. In particular it may help distinguish between WDM and ΛCDM models. Future higher resolution hydrodynamical simulations of WDM and ΛCDM will be important to understand the impact of baryons on the formation of dSphs. Further, more detailed measurements of the kinematics of dSphs, including proper motions of stars within dSphs and proper motions of the dSphs within the Milky Way potential, hold a great amount of promise.
# 5 Observables related to keV neutrino Dark Matter

Section Editors:
Marco Drewes, George Fuller

## 5.1 Phase space analysis (author: D. Gorbunov)

Sterile neutrinos are well-motivated Dark Matter candidates. Since we need a substantial DM component in galaxies, any particles — DM candidates — must be capable of residing there at present. This implies two very general constraints on properties of the DM particles. First, the DM lifetime must exceed the age of the Universe. For the case of sterile neutrinos, this requirement imposes a constraint on the lifetime that is weaker than that from the searches for emission lines from the DM decay (see section 5.3). Second, the DM “particle” X of mass $M_X$ must be slow and “compact” enough to be confined inside a galaxy.

The velocity distribution of DM particles is supposed to be the same as that of galactic stars and cold gas clouds, which may be approximated by the Maxwell distribution

$$F_X(\mathbf{v}) = \frac{1}{(\sqrt{2\pi} M_X \sigma_X)^3} \exp\left(-\frac{\mathbf{v}^2}{2\sigma_X^2}\right),$$

with $\sigma_X = \sigma_X(\mathbf{x})$ being the position-dependent velocity dispersion of stars in the galaxy. The dispersion takes the smallest value in the galaxy center and grows outward approaching $10^{-4} - 10^{-3}$ for small (dwarf) and large (normal) galaxies, respectively. Dark matter particles of higher velocities escape from the galaxy, whose gravity force fails to keep them inside. As regards to the sterile neutrino, the dispersion determines the width of the monochromatic line expected in the galactic X-ray spectrum due to the radiative decay $N \to \nu\gamma$ as described in section 5.3.

The “compactness” refers to the quantum nature of the DM particle. Namely, to locate a boson inside the galaxy its de Broglie wave length $\lambda_B = 2\pi/(M_X v_X)$ must be smaller than the size of the smallest observed galaxies, the dwarf galaxies, that is about 1 kpc. This yields a rather eccentric lower limit on the mass of the boson DM particle at the level of $M_X \gtrsim 10^{-22}$ eV.

In case of the fermion DM, a similar limit applicable to sterile neutrinos, is much more constraining. Indeed, fermions are subject to the Pauli principle. Therefore, for the galactic DM, the maximum of the phase space density with respect to velocity (which happens at $\mathbf{v} = 0$),

$$f_X^{\max}(\mathbf{v}, \mathbf{x}) = \frac{\rho_X(\mathbf{x})}{M_X} F_X(0)$$

may not exceed the critical value

$$f_F^{\text{crit}} = \frac{g_X}{(2\pi)^3}, \qquad (5.1)$$

where $\rho_X(\mathbf{x})$ is the DM mass density and $g_X$ denotes the number of intrinsic degrees of freedom of the DM particle, that is $g_X = 2$ for the sterile neutrino (Majorana fermion). Thus the Pauli principle places the lower limit on the fermion DM mass,

$$\frac{(2\pi)^{3/8}}{g_X^{1/4}} \left(\frac{\rho_X}{\sigma_X^3}\right)^{1/4} \le M_X. \qquad (5.2)$$
This bound is independent of the production mechanism. It gives $M_X \gtrsim 25$ eV for the Milky Way galaxy, but becomes one order of magnitude stronger when applied to the most compact and dim structures, dwarf spheroidal galaxies.

Let us note that bound (5.2) refers to the present population of DM particles in galaxies, and hence is applicable to any fermion DM irrespectively of properties of the DM particles and the DM production mechanism operating in the early Universe. Actually, bound (5.2) can be strengthened by extending the consideration to the evolution of the phase space density of the DM particles during the structure formation as follows.

After decoupling from primordial plasma (or production in a non-thermal way) the DM spectrum is fixed, and later on its shape is frozen while the Universe expands: only the physical momenta get redshifted, $p \propto 1/a(t)$ with $a(t)$ being the scale factor, see e.g. [546]. Indeed, in the early Universe the plasma is homogeneous to a large extent (the spatial fluctuations are below about $10^{-4}$). Hence the distribution of the DM particles in the phase space is homogeneous too, and initially (e.g. at decoupling) is described in terms of the conformal momentum $\mathbf{k} = \mathbf{p}/a$ as $f_i(\mathbf{k})$. Inasmuch as the conformal $d^3\mathbf{x} d^3\mathbf{k}$ and physical $d^3\mathbf{X} d^3\mathbf{p}$ phase space volumes coincide,

$$d^3\mathbf{x} d^3\mathbf{k} = d^3(\mathbf{x} a) d^3\left(\frac{\mathbf{k}}{a}\right) = d^3\mathbf{X} d^3\mathbf{p},$$

the shape of the momentum distribution remains intact,

$$f(\mathbf{k}, t) = f_i(\mathbf{k}) = f_i\left(\frac{\mathbf{k}}{a(t)} a(t)\right) = f_i(\mathbf{p} a(t))$$

while the inhomogeneities are small.

However, later the distribution function becomes inhomogeneous, evolves in time, and looks different for DM particles in voids, filaments, galaxy clusters, dwarf galaxies, etc. Nevertheless, for the collisionless particles (as we believe the DM is and the sterile neutrino is for sure) the Liouville theorem guarantees that the phase space density remains constant in time along any given trajectory, though flows from one region to another. In practice one is interested not in this, the microscopic density, but rather in the density averaged over macroscopic regions of the phase space, the so-called coarse grained phase space density, $\tilde{f}(\mathbf{k}, \mathbf{x})$. Naturally, it grows in the underdense regions and gets diluted in the overdense regions, as illustrated in figure 16. Therefore, it always obeys the following inequality

$$\tilde{f}(\mathbf{k}, \mathbf{x}, t) \le \max_k f_i(\mathbf{k}), \qquad (5.3)$$

including the galactic DM component at present. Then l.h.s. of (5.3) can be extracted from observations, while r.h.s. is the prediction of a particular mechanism of the DM production operating in the early Universe.

The analysis of astronomical data allows to estimate for a particular galaxy the value of the following quantity

$$Q \equiv \frac{\rho_0}{\langle \mathbf{v}_\parallel^2 \rangle^{3/2}},$$

where $\rho$ is the mass density in the galaxy center region and $\langle \mathbf{v}_\parallel^2 \rangle$ is the averaged (over the objects in that region) squared velocity along the line of sight. For obtaining the strongest constraint on the DM particle mass from eq. (5.3) the relevant are the highest values of $Q$,
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Evolution of the phase space density in time: initially occupying the compact region (left panel) particles spread over the phase space (right panel). The volume remains intact, but the coarse grained phase space density decreases in the dense regions.</figcaption></figure>

which are provided by observations of the stellar dynamics in the smallest structures, dwarf spheroidal galaxies, [441]

$$Q = (0.005 - 0.02) \frac{M_{\odot}/\text{pc}^3}{(\text{km/s})}. \qquad (5.4)$$

The galaxies are compact and dim with mass strongly dominated even in the central part by the DM component. Hence, one has $\rho_0 = M_X n_X$ with $n_X$ standing for the number density of the DM particles in the galaxy center. For the spheroidal dwarfs one can substitute $\langle \mathbf{v}_{\parallel}^2 \rangle = \langle \mathbf{v}^2 \rangle / 3$, where $\langle \mathbf{v}^2 \rangle$ is the average squared velocity. For the DM particles with average squared momentum $\langle \mathbf{p}^2 \rangle = M_X^2 \langle \mathbf{v}^2 \rangle$.<fnref n="26" />

Collecting all terms together we arrive at

$$Q = 3^{3/2} M_X^4 \frac{n}{\langle \mathbf{p}^2 \rangle^{3/2}} \approx 3^{3/2} M_X^4 \tilde{f}(\mathbf{p}, \mathbf{X}, t_0),$$

where in the last equality we used the phase space distribution of DM particles in the galaxy center, $\mathbf{x} \simeq 0$ at present, $t = t_0$. Then from eq. (5.3) we obtain the lower limit on the DM particle mass known as the Tremaine-Gunn-type bound [75],

$$M_X \gtrsim \left( \frac{Q}{3^{3/2} \max \tilde{f}_i} \right)^{1/4}$$

This bound is valid for both bosons and fermions, but $\max \tilde{f}_i$ of course depends on the production mechanism. For the former it supersedes what we have above from the de Broglie waves, while for the latter it coincides with that from the Pauli principle if $\max \tilde{f}_i$ reaches the critical value (5.1), but generally is somewhat more restrictive given $\max \tilde{f}_i$ in the denominator. Applying this limit to the sterile neutrino nonresonantly produced in the early Universe with the corresponding spectrum (see section 6.1), one obtains:

$$M_X \gtrsim 6 \text{ keV} \left( \frac{0.2}{\Omega_{\text{DM}}} \right)^{1/3} \left( \frac{Q}{0.005 \frac{M_{\odot}/\text{pc}^3}{(\text{km/s})}} \right)^{1/3} \left( \frac{g_*(T_d)}{43/4} \right) \qquad (5.5)$$

<footnote n="26">²⁶In ref. [313] it has been suggested that a Maxwell Boltzmann phase density estimator is more accurate to be applied to dwarf galaxies.</footnote>
with all factors of order one. Thus, literally using the estimates (5.4), one concludes that the light sterile neutrino, $M_X \lesssim 6$ keV nonresonantly produced in the early Universe, is disfavored [312, 313] from analysis of the phase space density as the main component of galactic DM. The limit suffers from the systematic uncertainties in the parameter estimates (5.4) [311], but is still much stronger as compared to that from the Pauli blocking. Similar limits exist for other fermion warm DM candidates, and are usually relevant for the non-thermal production mechanisms, see ref. [547] for example of light gravitino.

## 5.2 Lyman-α forest constraints (author: M. Viel)

The Lyman-α forest, the main manifestation of the intergalactic medium (IGM), is produced by intervening filaments of neutral hydrogen along the line-of-sight to a distant source, typically a quasar (QSO). Since its discovery in the early 70s, it has been used to investigate the distribution of matter and the physical conditions of such Lyman-α “clouds” of gas. However, in the late 90s there has been a paradigm shift: while this observable was in the beginning used as a local probe of baryonic matter and of the physics the absorbers, it became very soon apparent that the Lyma-α forest retains cosmological information and could be used as a tracer of the large scale structure. This new cosmological role for the IGM has been mainly allowed by two factors: on one side the new high-resolution high signal-to-noise quasar spectra and on the other side the progress made by the first set of hydrodynamic simulations [548]. A unique features of the IGM is that it is very complementary to other cosmological observables: it mainly probes small-medium scales from sub-Mpc up to few hundreds of Mpc and redshift in the range $z = 2-6$.

The observable quantity is the transmitted flux along the line-of-sight to a distant bright source, i.e. $F = \exp(-\tau)$, that needs to be related to underlying matter. In order to get some insight, we start from the definition of optical depth in redshift-space at $u$ (in km s⁻¹):

$$ \tau(u) = \frac{\sigma_{0,\alpha} c}{H(z)} \int_{-\infty}^{\infty} dy n_{\text{HI}}(y) \mathcal{V} \left[ u - y - v_{\parallel}^{\text{IGM}}(y), b(y) \right] dy, \quad (5.6) $$

where $\sigma_{0,\alpha} = 4.45 \times 10^{-18} \text{ cm}^2$ is the hydrogen Lyα cross-section, $H(z)$ the Hubble constant at $z$, $y$ the real-space coordinate (in km s⁻¹), $\mathcal{V}$ is the standard Voigt profile normalized in real-space, $b = (2k_B T/mc^2)^{1/2}$ is the velocity dispersion in units of c. In the highly ionized case (which is of interest here) the local density of neutral hydrogen can be easily related to the local gas/IGM density as:

$$ n_{\text{HI}}(\mathbf{x}, z) \approx 10^{-5} \bar{n}_{\text{IGM}}(z) \left( \frac{\Omega_{0b} h^2}{0.019} \right) \left( \frac{\Gamma_{-12}}{0.5} \right)^{-1} \left( \frac{T(\mathbf{x}, z)}{10^4 \text{K}} \right)^{-0.7} \left( \frac{1+z}{4} \right)^3 (1 + \delta_{\text{IGM}}(\mathbf{x}, z))^2, \quad (5.7) $$

with $\Gamma_{-12}$, the hydrogen photoionization rate in units of $10^{-12} \text{ s}^{-1}$ and $T(\mathbf{x}, z) = T_0(z)(1 + \delta_{\text{IGM}}(\mathbf{x}, z))^{\gamma(z)-1}$, where both the temperature at mean density $T_0$ and the adiabatic index $\gamma$ depend on the whole IGM ionization history. From the eqs. above, we have understood how the observed flux is derived from the underlying gas distribution. The latter can be modeled either analytically (with a log-normal model for example) or by means of high resolution hydrodynamic simulations that incorporate baryonic pressure (i.e. Jeans smoothing). In semi-analytical models [549, 550] this is usually done in Fourier space where the IGM component reads:

$$ \delta_{0}^{\text{IGM}}(\mathbf{k}, z) = \frac{\delta_{0}^{\text{DM}}(\mathbf{k}, z)}{1 + k^2/k_j^2(z)} \equiv W_{\text{IGM}}(k, z) D_{+}(z) \delta_{0}^{\text{DM}}(\mathbf{k}) \quad (5.8) $$
with $D_+(z)$ the linear growing mode of DM density fluctuations and $\delta_0^{\text{DM}}(\mathbf{k})$ is the Fourier transformed DM linear overdensity at $z=0$. The low-pass filter $W_{\text{IGM}}(k, z) = (1+k^2/k_J^2)^{-1}$ depends on the comoving Jeans length in the following way:

$$k_J^{-1}(z) = H_0^{-1} \left[ \frac{2\gamma k_B T_0(z)}{3\mu m_p \Omega_{0m}(1+z)} \right]^{1/2}, \quad (5.9)$$

with $k_B$ the Boltzmann constant, $T_0$ the temperature at mean density, $\mu$ the molecular weight of the IGM, $\Omega_{0m}$ the present-day matter density parameter and $\gamma$ the ratio of specific heats.

It is more practical to consider flux fluctuations $\delta F = F/\langle F \rangle - 1$ to those of matter $\delta$ and eqs. above can be simplified to obtain the so-called Fluctuating Gunn-Peterson Approximation:

$$F(\mathbf{x}, z) = \exp \left[ -A \left( 1 + \delta^{\text{IGM}}(\mathbf{x}, z) \right)^{\beta} \right], \quad (5.10)$$

where $\beta \approx 2-0.7(\gamma-1)$, while $A$ is a normalization constant of order unity which determines the mean flux in the considered redshift interval. It is easy to derive, by Taylor expanding the equation above, that $\delta F(\mathbf{x}, z) \approx -A\beta\delta_0^{\text{IGM}}(\mathbf{x}, z)$. The flux power is usually measured in 1D by using high-resolution or low-resolution spectra, while it has been recently exploited by the SDSS-III/BOSS collaboration in full 3D to allow the discovery of Baryonic Acoustic Oscillations in Lyman-$\alpha$ flux [551, 552]. The relation between 1D and 3D power is:

$$P_{3D}(k) = -\frac{2\pi}{k} \frac{d}{dk} p_{1D}(k), \quad (5.11)$$

which is valid for both flux and matter power spectra and shows that the 1D power is an integral of the 3D and effects like a Warm Dark Matter (WDM) cutoff could impact at relatively larger scale and thereby could be constrained by IGM data.

From the considerations made above, it is pretty clear that the Lyman-$\alpha$ forest can be used to place tight constraints on cold DM coldness [553], however it must also be noticed that the temperature could affect flux power in a similar way: a hot IGM has less power than a cold one since lines will be broadened more in the first case. This signature could be degenerate with the one induced by WDM. Also, it should be kept in mind that structure formation is not directly sensitive to the properties of DM particles. It is mainly sensitive to their free streaming length. A bound on the free streaming horizon $\lambda_{\text{FS}}$ can be converted into a bound on the DM particles' mass if their momentum distribution is known. Heavy CDM particles are in good approximation at rest at all relevant times. Most of the constraints on lighter particles presented so far are obtained by assuming a thermal relic WDM model for which the power is suppressed below a given scale. That is, it has been assumed that the momentum distribution of the DM particles is proportional to a thermal spectrum. In this case, one can estimate [554]:

$$\lambda_{\text{FS}} \sim 1\text{Mpc} \frac{\text{keV}}{m_{\text{DM}}} \frac{\langle p_{\text{DM}} \rangle}{\langle p_{\nu} \rangle}, \quad (5.12)$$

where $\langle p_{\text{DM}} \rangle$ is the average momentum of the DM particles and $\langle p_{\nu} \rangle \sim 1$ keV is the comoving momentum of active neutrinos around the time when thermally produced sterile neutrino DM particles become non-relativistic. For highly non-thermal distribution functions, which can appear if the sterile neutrinos are produced resonantly (see section 6.1.4) or in a particle decay (see section 6.3), there is no simple relation like (5.12).

The first of such limits (notice that strong absorption systems were used before to place constraints on DM nature in [557]), is presented in [555] where a thermal relic of masses > 0.75
<table><caption>**Table 4.** Summary of limits obtained on the mass of a thermal relic (WDM candidate) by using a comprehensive set of Lyman-α forest data both at low and high resolution. All the numbers apart from [555] quote $2\sigma$ C.L. lower limits. In [59] also a frequentist analysis of the SDSS data is presented. In [328] the method used to derive the constraints is based on the effective bias method proposed by [556], in [329] there is an approximate hydrodynamical scheme but the likelihood analysis is performed in a wide parameter space, while in the other works the simulations are full hydro and the method is based on a Taylor expansion of the flux power.</caption><thead><tr><th>reference</th><th>limit (keV)</th><th>method</th><th>data</th></tr></thead><tbody><tr><td>[555]</td><td>0.75</td><td>N-body only, no marginalization</td><td>8 high res.</td></tr><tr><td>[328]</td><td>0.55</td><td>full hydro</td><td>30 high res.</td></tr><tr><td>[327]</td><td>2.5</td><td>approximate hydro</td><td>3000 low res.</td></tr><tr><td>[329]</td><td>2</td><td>full hydro, approximate method</td><td>3000 low res.</td></tr><tr><td>[332]</td><td>4</td><td>full hydro, approximate method</td><td>3000 low res., 55 high res.</td></tr><tr><td>[59]</td><td>2.2</td><td>full hydro, frequentist analysis</td><td>3000 low res.</td></tr><tr><td>[331]</td><td>3.3</td><td>full hydro, improved analysis</td><td>3000 low res., 25 high res., high z</td></tr></tbody></table>

keV is shown to be consistent with the data, however this number was obtained with a limited set of N-body only simulations and does not rely on a full marginalization over the many nuisance and astrophysical parameters. In [328], by using the data points (i.e. linear matter power) derived in [558], it was demonstrated that with full hydro and proper marginalization in the parameter space (although obtained with the effective bias model of [556], on a larger sample of QSO spectra) this number was relaxed to 0.55 keV ($2\sigma$ C.L.), in this paper also a lower limit on a non resonantly produced sterile neutrino was quantified to be 2 keV.

The main new results soon after these analyses were performed, was the advent of SDSS (Sloan Digital Sky Survey), which allowed an unprecedented measurement of linear matter power spectrum at $z=3$ with very good accuracy. This was derived from a set of about 3000 low resolution low signal-to-noise spectra [559] and by using a new method. In particular, the likelihood space was explored fully by accurately modeling also astrophysical and instrumental effects, while on the other side the calibration of the flux power was determined by running only a few hydro simulations, while the parameter space was explored by using approximate hydro schemes (that were computationally less expensive). The main point to stress is however that for the first time, the relatively wide range explored by the data $z=2.2-4.2$ allowed to break the degeneracies present between cosmology and thermal history, resulting in significantly tighter constraints. In [329], SDSS data returned the following constraints: 2.5 keV and 14 keV lower limits on a thermal relic and non-resonantly produced sterile neutrinos at a $2\sigma$ C.L., respectively. The same data set was then analyzed by [327] using a completely different method (albeit still in a Bayesian framework): a second-order Taylor expansion of the flux power in a relatively better explored (compared to [329]) hydrodynamical parameter space gave a ~20% weaker constraints of a mass of a thermal relic > 2 keV and of a sterile neutrino > 10 keV at $2\sigma$ C.L., respectively. In [560] some scenarios in which there is some entropy released is also explored along with non-zero lepton number models and confronted with IGM constraints. The investigation of how the flux power relates to the underlying matter power spectrum emphasized the need for QSO spectra at relatively high redshift, where the non-linear evolution is smaller and possibly closer to the linear cut-off. In [332] the SDSS data were complemented by a smaller set of 55 high resolution high redshift QSO spectra that resulted in a lower limit of 4 keV, however the higher redshift
bin was explored with hydro simulations that were not fully converged (although there was
a resolution correction implemented in the pipeline). This is the tightest limit on WDM
coldness ever obtained from IGM data.

After these results, in [59] the SDSS data set (without high resolution samples) was used
to get constraints in a frequentist framework for a non resonantly produced sterile neutrinos
obtaining $m > 8$ keV at 99.7% confidence limit in a pure ΛCDM model. Most notably,
this work also investigated mixed (cold and warm) models and found that any mass for a
sterile neutrinos could be allowed if the WDM fraction is below 60%. An important break-
through was made in [330] where it was realized that mixed models linear power spectra were
remarkably similar to those of a resonantly produced sterile neutrinos in the context of min-
imal extension of the standard model: in such cases the limits above showed that any sterile
neutrino mass > 2 keV would have been in agreement with the Lyman-α forest SDSS data.

The most recent and comprehensive analysis of Lyman-α forest data for WDM models
has been made in [331]. Compared to the previous works there are several improvements: the
hydrodynamic simulations have been run at very high resolution and do not need any reso-
lution correction nor they suffer of numerical fragmentation; the non-linearities both in the
matter power and flux power are modeled to unprecedented accuracy [561]; the data sets of
high-resolution spectra extends the previous ones (Keck HIRES and Magellan spectrograph);
the modeling carefully consider instrumental resolution, signal-to-noise and astrophysical
nuisance parameters (mainly thermal histories and ultraviolet background inhomogeneities).
The conclusion is that the lower limit is 3.3 keV at $2\sigma$ C.L. for a WDM thermal relic; this
has to be considered as a robust limit, obtained by marginalizing over a physically motivated
set of thermal histories and conservatively considering systematic errors both on the data
and modeling side. It is worth noticing that if one allows for a sudden jump in the tempera-
ture evolution the constraint above becomes weaker by about 1 keV. These (and almost all
other known) bounds rely on the assumption that the DM momentum distribution is not
too different from a thermal spectrum (or a superposition of two thermal spectra). Then a
relation of the type (5.12) can be used. Translating this number in a constraint on resonantly
produced sterile neutrinos with a spectrum as shown in figure 25 is not easy and could be
done by comparing the linear transfer function proper of a thermal relic of such mass with
the one of the wanted model. Consequences in terms of high redshift structure formation for
a thermal relic mass particle compatible with all the constraints have been recently explored
in [540], while forecasts in terms of 21 centimeter intensity mapping are presented in [562]
also for lower values of the mass.

In figure 17 we show the difference between a ΛCDM best fit model and a WDM model
of 2 keV in terms of agreement with flux power spectra at low and high resolution. It is clear
the power of the high redshift data that allows to severely constrain WDM coldness.

## 5.3 X-ray observations (authors: O. Ruchayskiy, T. Jeltema, A. Neronov, D. Iakubovskyi)

### 5.3.1 X-ray signals — overview

The main decay channel of sterile neutrinos with the mass below $2m_e$ is $N \rightarrow \nu\bar{\nu}$ (see figure 18, left) with different combinations of flavours. It determines the lifetime of the particle. The sterile neutrino DM has a sub-dominant radiative decay channel $N \rightarrow \nu\gamma$ (figure 18, right panel). The decay width of this process is about 128 times smaller than the main into active
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Best fit model for the data sets used in the analysis (SDSS+HIRES+MIKE) shown as green curves. We also show a WDM model that has the best fit values of the green model except for the WDM (thermal relic) mass of 2 keV (red dashed curves). These data span about two orders of magnitude in scale and the period 1.1-3.1 Gyrs after the Big Bang. From this plot is apparent how the WDM model does not fit the data at small scales and high redshift.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 18.** Decay channels of the sterile neutrino N with the mass below twice the electron mass. Left panel: dominant decay channel to three (anti)neutrinos. Right panel shows radiative decay channel that allows to look for the signal of sterile neutrino DM in the spectra of DM dominated objects.</figcaption></figure>

neutrinos $\nu_a$ and photon with energy $E = m_s/2$, with the width [321–323, 359, 563–565]

$$ \Gamma_{N \to \gamma \nu_a} = \frac{9 \alpha G_F^2}{256 \cdot 4\pi^4} \sin^2 2\theta m_s^5 = 5.5 \times 10^{-22} \theta^2 \left[ \frac{m_s}{1 \text{ keV}} \right]^5 \text{ s}^{-1}. \quad (5.13) $$
If the sterile neutrino is a main ingredient of the DM, it is potentially detectable in various X-ray observations [324, 566, 567].

There are several types of the X-ray signal produced by the sterile neutrino DM. First, decays of DM particles throughout the history of the Universe should produce a contribution to the diffuse X-ray background (XRB) [324, 566, 568, 569]. The DM decay contribution to the XRB is (see [569] and references therein)

$$F_{\text{XRB}} \approx \frac{\Gamma_{N \to \gamma \nu_a} \rho_{\text{DM}}^0}{2\pi H_0} \approx 8 \times 10^{-11} \left[ \frac{\theta^2}{10^{-11}} \right] \left[ \frac{m_s}{7 \text{ keV}} \right]^5 \frac{\text{erg}}{\text{cm}^2 \cdot \text{s} \cdot \text{sr}}, \quad (5.14)$$

where $\rho_{DM}^0$, $H_0$ are the DM density in the universe and the Hubble constant today. Neutrinos decaying at different red shifts produce a broad X-ray line with extended “red” tail. Such a feature in the XRB spectrum is, in principle, readily detectable (and distinguishable) from the broad-band continuum of observed XRB [570, 571]. The non-detection of the DM decay feature in the XRB signal enables to put a bound on $\theta$ and $m_s$ roughly at the level of [569] (see also [572])

$$\text{XRB bound:} \quad \Omega_s \sin^2(2\theta) \lesssim 3 \times 10^{-5} \left[ \frac{1 \text{ keV}}{m_s} \right]^5, \quad (5.15)$$

where $\Omega_s \le \Omega_{DM}$ is the present day density of sterile neutrino DM.

As first pointed out refs. [315, 566], clustering of the DM at small red shifts results in the enhancement of the DM decay signal in the direction of large mass concentrations, such as galaxy clusters, see also [358, 573–575]. Typical overdensity, $\mathcal{R}$, ($\mathcal{R} \equiv \rho/\rho_{DM}^0$) of a galaxy cluster is at the level of $\mathcal{R} \sim 10^3$, while typical cluster size is about $R_{cl} \sim \text{Mpc}$. An estimate for the DM decay flux from a galaxy cluster is

$$F_{cl} = \frac{\Gamma_{N \to \gamma \nu_a} M_{DM, FoV}}{4\pi D_L^2 m_s} \approx 1.4 \times 10^{-7} \left[ \frac{\theta^2}{10^{-11}} \right] \left[ \frac{m_s}{7 \text{ keV}} \right]^4 \left[ \frac{D_L}{100 \text{ Mpc}} \right]^{-2} \left[ \frac{M_{DM, FoV}}{10^{13} M_{\odot}} \right] \frac{\text{ph}}{\text{cm}^2 \text{s}}$$

$(M_{DM,FoV} \approx \mathcal{R}\rho_{DM}^0 R_{cl} D_{\theta}^2 \Omega_{FoV}$ is the mass of DM within telescope's field of view (FoV), $D_L$, $D_{\theta}$ are the luminosity and angular diameter distances to the cluster). The DM decay flux from the direction of individual clusters is typically of the same order as the DM flux from the XRB contribution:

$$\frac{m_s F_{cl} / (2\Omega_{FoV})}{F_{\text{XRB}}} \sim \mathcal{R} R_{cl} H_0 \sim 1. \quad (5.16)$$

because numerically $\mathcal{R} \sim (R_{cl}H_0)^{-1}$. The results in refs. [566, 572] imply that foreground objects dominate at the level of prefactors.

Although the total DM decay flux from the XRB contribution and from the individual clusters are of the same order, the spectral shape of an expected signal is quite different. This makes stacked observations of galaxy clusters and observations of field galaxies more sensitive [566]. The flux from the cluster is in the form of a narrow line with the width determined by the spectral resolution of an X-ray detector or by the velocity dispersion of the DM halo of the cluster (if the detector energy resolution is better than $E/\Delta E \sim 10^2$). At the same time, the DM decay contribution into XRB is produced by the decays at red shifts $z \sim 0 \div 1$ and, as a result the DM decay line is broadened to $\Delta E \sim m_s/2$. Thus, even if the compact DM sources at $z \sim 0$ give just moderate enhancement of the DM decay flux, the enhancement of the signal in the narrow energy band centered on the line energy $E = m_s/2$
could be large [345, 569]. The results presented in refs. [566, 572] suggest that there is an even larger enhancement due to the dominant contributino from foreground objects.

This enhancement is however, tempered by the fact that most of the galaxy clusters show strong continuum and line emission from the hot intracluster gas. The temperatures of the intracluster medium are in the range $T_{gas} \sim G R \rho_{DM}^{0} R_{cl}^{2} m_{p} \sim 1-10$ keV. The diffuse continuum emission flux from the nearby cluster cores is orders of magnitude stronger than the XRB emission from behind the core [566, 576–582]. Thus, an increase of the flux in the DM decay line in the narrow energy band of the width about the energy resolution of the X-ray telescope is accompanied by an equally strong increase of the X-ray continuum contribution to the flux in the same energy band. An additional problem is the presence of emission lines from the hot intracluster medium in the spectra of clusters [578, 583–585]. These lines could be confused with the DM decay line, especially while the spectral resolution of the X-ray instruments is much larger than the line broadening.

Nearby DM mass concentrations, like the dwarf spheroidal (dSph) satellites of the Milky Way, have much smaller DM masses and their luminosity in the DM decay line is orders of magnitude lower than that of the galaxy clusters [586]. However, their proximity to the Milky Way makes increases the flux of the line, so that at the end the flux from the direction of individual dSph galaxies is expected to be comparable to the flux from individual galaxy clusters [345, 573, 574, 587–597]. Re-using the equation (5.3.1) for the flux estimate for the case of the dSph galaxies one finds

$$F_{dSph} \approx 1.4 \times 10^{-7} \left[ \frac{\theta^2}{10^{-11}} \right] \left[ \frac{m_s}{7 \text{ keV}} \right]^4 \left[ \frac{D_L}{100 \text{ kpc}} \right]^{-2} \left[ \frac{M_{DM,FoV}}{10^7 M_{\odot}} \right] \frac{\text{ph}}{\text{cm}^2\text{s}} \quad (5.17)$$

In addition the flux from the DM decays in the largest isolated mass concentrations, galaxy clusters, should produce a potentially detectable signal from individual clusters [566, 581, 598, 599] and from the entire galaxy cluster population [190]. Finally, the nearest DM halos, which are the Milky Way and other nearby galaxies, are also expected to produce DM decay signal with the strength comparable to the signal from the galaxy clusters, in spite of the much lower mass of the DM in these structures [347, 566, 572, 581, 587–589, 591, 594–597, 600–602].

The DM decay flux from the Milky Way halo is also comparable to the flux from isolated distant sources, like dSph galaxies or galaxy clusters [345]. The flux from DM decay in the Milky Way within the telescope field-of-view $\Omega_{FoV}$,

$$F = \frac{\Gamma \Omega_{FoV} S}{4 \pi m_{DM}} \approx 1.1 \times 10^{-2} \left[ \frac{\theta^2}{10^{-11}} \right] \left[ \frac{S}{10^{22} \text{ GeV/cm}^2} \right] \left[ \frac{m_s}{7 \text{ keV}} \right]^4 \frac{\text{ph}}{\text{cm}^2\text{s sr}} \quad (5.18)$$

is determined by the column density of the DM

$$S = \int_{0}^{\infty} \rho_{DM} \left( \sqrt{r_{\odot}^{2} - 2zr_{\odot} \cos \phi + z^{2}} \right) dz \quad (5.19)$$

where $r_{\odot} = 8.5$ kpc is the distance from the Sun to the center of our Galaxy and the “off-Galactic Center” angle $\phi$ is expressed through the galactic coordinates $(l, b)$ as $\cos \phi = \cos b \cos l$.

Search of the DM decay signal in the keV–MeV mass range was conducted using all the four types of the signal (the XRB contribution, the signal from galaxy clusters, from the
<table><caption>**Table 5.** Summary of existing X-ray observations of DM decays in different objects performed by different groups.</caption><thead><tr><th>ref.</th><th>Object</th><th>Instrument</th><th>Cleaned exp, ks</th></tr></thead><tbody><tr><td>[569]</td><td>Diffuse X-ray background</td><td>HEAO-1, XMM-Newton/EPIC</td><td>224, 1450</td></tr><tr><td>[581]</td><td>Coma, Virgo</td><td>XMM-Newton/EPIC</td><td>20, 40</td></tr><tr><td>[345]</td><td>Large Magellanic Cloud</td><td>XMM-Newton/EPIC</td><td>20</td></tr><tr><td>[598]</td><td>Milky Way</td><td>Chandra/ACIS-S3</td><td>Not specified</td></tr><tr><td>[600]</td><td>M31 (central 5')</td><td>XMM-Newton/EPIC</td><td>35</td></tr><tr><td>[598]</td><td>Abell 520</td><td>Chandra/ACIS-S3</td><td>67</td></tr><tr><td>[587]</td><td>Milky Way, Ursa Minor</td><td>XMM-Newton/EPIC</td><td>547, 7</td></tr><tr><td>[572]</td><td>Milky Way</td><td>Chandra/ACIS, X-ray microcalorimeter</td><td>1500, 0.1</td></tr><tr><td>[599]</td><td>1E 0657-56 ("Bullet cluster")</td><td>Chandra/ACIS-I</td><td>450</td></tr><tr><td>[602]</td><td>Milky Way</td><td>X-ray microcalorimeter</td><td>0.1</td></tr><tr><td>[601]</td><td>Milky Way</td><td>INTEGRAL/SPI</td><td>5500</td></tr><tr><td>[346]</td><td>M31 (central 5 – 13')</td><td>XMM-Newton/EPIC</td><td>130</td></tr><tr><td>[347]</td><td>Milky Way</td><td>INTEGRAL/SPI</td><td>12200</td></tr><tr><td>[588]</td><td>Ursa Minor</td><td>Suzaku/XIS</td><td>70</td></tr><tr><td>[589]</td><td>Draco</td><td>Chandra/ACIS-S</td><td>32</td></tr><tr><td>[590]</td><td>Willman 1</td><td>Chandra/ACIS-I</td><td>100</td></tr><tr><td>[604, 605]</td><td>Galactic center</td><td>Suzaku/XIS</td><td>190</td></tr><tr><td>[591]</td><td>M31, Fornax, Sculptor</td><td>XMM-Newton/EPIC, Chandra/ACIS</td><td>400, 50, 162</td></tr><tr><td>[592]</td><td>Willman 1</td><td>Chandra/ACIS-I</td><td>100</td></tr><tr><td>[593]</td><td>Segue 1</td><td>Swift/XRT</td><td>5</td></tr></tbody></table>

dSph galaxies and from the Milky Way halo) and a wide range of X-ray telescopes: XMM-Newton [345, 346, 569, 581, 587, 594, 600], Chandra [589–591, 598, 598, 599], Suzaku [588, 595], Swift [593], INTEGRAL [347, 601], HEAO-1 [569] and Fermi/GBM [603], as well as a rocket-borne X-ray microcalorimeter [572, 602]. A summary of the reported observations is given in tables 5–9.

Most of the searches of the DM decay line in X-rays did not result in positive detections of the line. The non-detection results are conventionally expressed in the form of the upper limit on the DM sterile neutrino mixing angle as a function of the mass. These limits are shown in figure 19.

### 5.3.2 3.5 keV line

The detection of an unidentified line was reported recently in the stacked spectrum of galaxy clusters [190], in the individual spectra of nearby galaxy clusters [190, 191, 622] (see also [612]), in the Andromeda galaxy [191], and in the Galactic Center region [608–
<table><caption>**Table 6.** Summary of existing X-ray observations of DM decays in different objects performed by different groups (extended).</caption><thead><tr><th>ref.</th><th>Object</th><th>Instrument</th><th>Cleaned exp, ks</th></tr></thead><tbody><tr><td>[606]</td><td>M33</td><td>XMM-Newton/EPIC</td><td>20-30</td></tr><tr><td>[485]</td><td>M31 (12 – 28' off-center)</td><td>Chandra/ACIS-I</td><td>53</td></tr><tr><td>[594]</td><td>Willman 1</td><td>XMM-Newton/EPIC</td><td>60</td></tr><tr><td>[595]</td><td>Ursa Minor, Draco</td><td>Suzaku/XIS</td><td>200, 200</td></tr><tr><td>[607]</td><td>Stacked galaxies</td><td>XMM-Newton/EPIC</td><td>8500</td></tr><tr><td>[313]</td><td>M31</td><td>Chandra/ACIS-I</td><td>404</td></tr><tr><td>[190]</td><td>Stacked clusters</td><td>XMM-Newton/EPIC, Chandra/ACIS</td><td>6000, 1370</td></tr><tr><td>[608]</td><td>M31, Perseus, Milky Way</td><td>XMM-Newton/EPIC</td><td>1226, 264, 7850</td></tr><tr><td>[609]</td><td>Galactic center</td><td>Chandra/ACIS</td><td>751</td></tr><tr><td>[610]</td><td>Galactic center, M31, Tycho</td><td>XMM-Newton/EPIC</td><td>690, 490, 175</td></tr><tr><td>[596]</td><td>Stacked dSphs</td><td>XMM-Newton/EPIC</td><td>410</td></tr><tr><td>[611]</td><td>Stacked galaxies</td><td>XMM-Newton/EPIC, Chandra/ACIS-I</td><td>14600, 15000</td></tr><tr><td>[612]</td><td>Perseus, Coma, Ophiuchus, Virgo</td><td>Suzaku/XIS</td><td>740, 164, 83, 90</td></tr><tr><td>[613]</td><td>Galactic center, Perseus</td><td>XMM-Newton/EPIC</td><td>Not specified</td></tr><tr><td>[614]</td><td>Perseus</td><td>Suzaku/XIS</td><td>520</td></tr><tr><td>[615]</td><td>Galactic Bulge</td><td>Suzaku/XIS</td><td>Not specified</td></tr><tr><td>[603, 616]</td><td>Milky Way</td><td>Fermi/GBM</td><td>4600</td></tr><tr><td>[617]</td><td>Milky Way</td><td>Suzaku/XIS</td><td>31500</td></tr><tr><td>[597]</td><td>Draco</td><td>XMM-Newton/EPIC</td><td>87</td></tr></tbody></table>

610], see also [623] for a recent review. The position of the line is $E = 3.55$ keV with an uncertainty in position ~ 0.05 keV. If the line is interpreted as originating from a two-body decay of a DM particle, then the latter has its mass at about $m_s \simeq 7.1$ keV and the lifetime $\tau_{DM} \sim 10^{27.8 \pm 0.3}$ sec [191] (and more narrow range based on the results of [190]). Using eq. (5.13) and converting the lifetime to the sterile neutrino mixing angle, one finds $\sin^2(2\theta) \simeq (2-20) \times 10^{-11}$ [191] or $\sin^2(2\theta) \simeq (4-11) \times 10^{-11}$ (based on results of [190]). The estimates of the errors on the mixing angle $\theta^2$ derived from different types of data show a visible scatter, which is, however, consistent with the quoted errors. For example, ref. [190, figure 14] shows the scatter in the possible mixing angles derived from the statistics of the data and not including the systematic error due to the uncertainties of the cluster mass estimates. The 3.5 keV line was also detected in the Galactic Center [608–610] with a flux that is consistent with the DM interpretation of the line in clusters and in the Andromeda galaxy [608]. Namely, it is sufficiently strong (to be consistent with the cluster detection) and sufficiently weak to be consistent with the non-detection from the outskirts of the Milky Way [191]. However, in the Galactic Center the flux of the line is also consistent with the expectations for plasma emission lines with no need for a dark matter component [608, 610].
<table><caption>**Table 7.** Summary of existing X-ray observations of DM decays in dwarf galaxies performed by different groups.</caption><thead><tr><th>ref.</th><th>Object</th><th>Instrument</th><th>Cleaned exp, ks</th></tr></thead><tbody><tr><td>[345]</td><td>Large Magellanic Cloud</td><td>XMM-Newton/EPIC</td><td>20</td></tr><tr><td>[587]</td><td>Ursa Minor</td><td>XMM-Newton/EPIC</td><td>7</td></tr><tr><td>[591]</td><td>Fornax</td><td>XMM-Newton/EPIC</td><td>50</td></tr><tr><td>[596]</td><td>Stacked dSphs</td><td>XMM-Newton/EPIC</td><td>410</td></tr><tr><td>[594]</td><td>Willman 1</td><td>XMM-Newton/EPIC</td><td>60</td></tr><tr><td>[597]</td><td>Draco</td><td>XMM-Newton/EPIC</td><td>87</td></tr><tr><td>[588]</td><td>Ursa Minor</td><td>Suzaku/XIS</td><td>70</td></tr><tr><td>[595]</td><td>Ursa Minor, Draco</td><td>Suzaku/XIS</td><td>200, 200</td></tr><tr><td>[589]</td><td>Draco</td><td>Chandra/ACIS-S</td><td>32</td></tr><tr><td>[590]</td><td>Willman 1</td><td>Chandra/ACIS-I</td><td>100</td></tr><tr><td>[591]</td><td>Sculptor</td><td>Chandra/ACIS</td><td>162</td></tr><tr><td>[592]</td><td>Willman 1</td><td>Chandra/ACIS-I</td><td>100</td></tr><tr><td>[593]</td><td>Segue 1</td><td>Swift/XRT</td><td>5</td></tr><tr><td>[618]</td><td>Draco</td><td>XMM-Newton/EPIC</td><td>2100 (MOS), 580 (PN)</td></tr><tr><td>[619]</td><td>Draco</td><td>XMM-Newton/EPIC</td><td>1990 (MOS), 650 (PN)</td></tr></tbody></table>

We note that the line is also recently detected in the *NuStar* blank-sky [620] and Galactic center [621] datasets. While it may have an instrumental origin, related to Solar activity [620], it was shown in [620] that in the ‘no sun’ dataset the line is present with the same intensity, and it appears several times stronger near the Galactic center, as figure 4 of [621] demonstrates.

A difficulty in interpreting the origin of a weak emission line is inherent uncertainty in the astrophysical backgrounds, in particular in the flux of plasma emission lines. The strongest uncertainty comes from two potassium lines, K XVIII at 3.47 and 3.51 keV. Given the spectral resolution of the XMM-Newton, within the systematic uncertainty the flux could be attributed to emission from these K XVIII plasma lines. Ref. [610] argued that considering a larger range of plasma temperatures reduces the tension between the observed 3.5 keV line flux in clusters and the expectations from known plasma lines (see however the subsequent discussion in [624–626]). The interpretation of the 3.5 keV line as a plasma line would imply that its surface brightness profile must trace the density of the plasma (more precisely: the distribution of potassium). If it is a DM line it should, on the other hand, trace the overall distribution of DM, which dominates the gravitating matter in Perseus and other clusters. This point is disputed: while the analysis in [608] suggests that the line traces the overall matter distribution (pointing towards DM), the authors of [613] conclude that the morphology is incompatible with the DM interpretation. Explicitly, the authors of ref. [613] argue that the morphology of the 3.5 keV signal is incompatible with DM in the Galactic Center with a spatial distribution that places strong upper limits on the possible dark matter flux. The authors of [613] also show that the 3.5 keV emission in Perseus is concentrated on the cool-core and less extended than expected for dark matter decay. Possible flaws in the analysis in [613] are discussed in [627]. Specifically, ref. [627] discusses some challenges with the template analysis of [613], which the authors of [613] however do not agree upon. The
<table><caption>**Table 8.** Summary of existing X-ray observations of DM decays in spiral galaxies and Milky Way halo performed by different groups.</caption><thead><tr><th>ref.</th><th>Object</th><th>Instrument</th><th>Cleaned exp, ks</th></tr></thead><tbody><tr><td>[600]</td><td>M31 (central 5')</td><td>XMM-Newton/EPIC</td><td>35</td></tr><tr><td>[587]</td><td>Milky Way</td><td>XMM-Newton/EPIC</td><td>547</td></tr><tr><td>[346]</td><td>M31 (central 5 – 13')</td><td>XMM-Newton/EPIC</td><td>130</td></tr><tr><td>[591]</td><td>M31</td><td>XMM-Newton/EPIC</td><td>400</td></tr><tr><td>[606]</td><td>M33</td><td>XMM-Newton/EPIC</td><td>20-30</td></tr><tr><td>[607]</td><td>Stacked galaxies</td><td>XMM-Newton/EPIC</td><td>8500</td></tr><tr><td>[313]</td><td>M31</td><td>Chandra/ACIS-I</td><td>404</td></tr><tr><td>[608]</td><td>M31, Milky Way</td><td>XMM-Newton/EPIC</td><td>1226, 7850</td></tr><tr><td>[610]</td><td>Galactic center, M31</td><td>XMM-Newton/EPIC</td><td>690, 490</td></tr><tr><td>[611]</td><td>Stacked galaxies</td><td>XMM-Newton/EPIC</td><td>14600</td></tr><tr><td>[613]</td><td>Galactic center</td><td>XMM-Newton/EPIC</td><td>Not specified</td></tr><tr><td>[572]</td><td>Milky Way</td><td>Chandra/ACIS</td><td>1500</td></tr><tr><td>[598]</td><td>Milky Way</td><td>Chandra/ACIS-S3</td><td>Not specified</td></tr><tr><td>[485]</td><td>M31 (12 – 28' off-center)</td><td>Chandra/ACIS-I</td><td>53</td></tr><tr><td>[609]</td><td>Galactic center</td><td>Chandra/ACIS</td><td>751</td></tr><tr><td>[611]</td><td>Stacked galaxies</td><td>Chandra/ACIS-I</td><td>15000</td></tr><tr><td>[604, 605]</td><td>Galactic center</td><td>Suzaku/XIS</td><td>190</td></tr><tr><td>[615]</td><td>Galactic Buldge</td><td>Suzaku/XIS</td><td>Not specified</td></tr><tr><td>[617]</td><td>Milky Way</td><td>Suzaku/XIS</td><td>31500</td></tr><tr><td>[602]</td><td>Milky Way</td><td>X-ray microcalorimeter</td><td>0.1</td></tr><tr><td>[601]</td><td>Milky Way</td><td>INTEGRAL/SPI</td><td>5500</td></tr><tr><td>[347]</td><td>Milky Way</td><td>INTEGRAL/SPI</td><td>12200</td></tr><tr><td>[603, 616]</td><td>Milky Way</td><td>Fermi/GBM</td><td>4600</td></tr><tr><td>[620]</td><td>Milky Way</td><td>NuSTAR</td><td>7500</td></tr><tr><td>[621]</td><td>Galactic center</td><td>NuSTAR</td><td>400</td></tr></tbody></table>

Potassium interpretation also cannot explain the origin of the line in the Andromeda galaxy reported in [191], the significance of which is, however, disputed [610, 624].

Systematic errors in instrumental calibration and/or systematics induced by the analysis procedure may impact the significance of weak lines. The calibration systematics was explored in [191] who demonstrated that no line is detected in an extremely long exposure combination of the off-center observations of the Milky Way ("blank sky" dataset). A 3.5 keV line is not detected in the stacked dwarf galaxies [596] arguing against a global calibration problem. It has also not been seen in a long exposure towards M31 reported in ref. [313]. In order to assess possible systematic effects, the ref. [610] considered the XMM-Newton spectra of the Tycho supernova remnant, assuming that this system displays similar plasma emission lines to those detected in clusters and the GC but would have no expected DM signal. A line at
<table><caption>**Table 9.** Summary of existing X-ray observations of DM decays in galaxy clusters and diffuse X-ray background performed by different groups.</caption><thead><tr><th>ref.</th><th>Object</th><th>Instrument</th><th>Cleaned exp, ks</th></tr></thead><tbody><tr><td>[569]</td><td>Diffuse X-ray background</td><td>HEAO-1</td><td>224</td></tr><tr><td>[569]</td><td>Diffuse X-ray background</td><td>XMM-Newton/EPIC</td><td>1450</td></tr><tr><td>[581]</td><td>Coma, Virgo</td><td>XMM-Newton/EPIC</td><td>20, 40</td></tr><tr><td>[190]</td><td>Stacked clusters</td><td>XMM-Newton/EPIC</td><td>6000</td></tr><tr><td>[608]</td><td>Perseus</td><td>XMM-Newton/EPIC</td><td>264</td></tr><tr><td>[613]</td><td>Perseus</td><td>XMM-Newton/EPIC</td><td>Not specified</td></tr><tr><td>[598]</td><td>Abell 520</td><td>Chandra/ACIS-S3</td><td>67</td></tr><tr><td>[599]</td><td>1E 0657-56 ("Bullet cluster")</td><td>Chandra/ACIS-I</td><td>450</td></tr><tr><td>[190]</td><td>Stacked clusters</td><td>Chandra/ACIS</td><td>1370</td></tr><tr><td>[612]</td><td>Perseus, Coma, Ophiuchus, Virgo</td><td>Suzaku/XIS</td><td>740, 164, 83, 90</td></tr><tr><td>[614]</td><td>Perseus</td><td>Suzaku/XIS</td><td>520</td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>**Figure 19.** Limits on the mixing angle $\sin^2(2\theta)$ as a function of sterile neutrino DM mass. The bounds are based on the works [313, 345–347, 485, 569, 572, 581, 587–590, 592–594, 596, 598, 598–602, 606, 607, 616]. All bounds are smoothed and additionally divided by a factor of 2 to take into account possible uncertainties in the DM content of a given object. A lower bound on $\sin^2\theta$ for a given DM mass can be imposed if the DM is produced via active-sterile mixing, see section 6.1.3, but is model dependent.</figcaption></figure>

~ 3.55 keV is detected in Tycho (with a significance of $\Delta\chi^2 = 19$).$^{27}$ Ref. [610] argued for either systematic errors in the estimated significance and flux of the 3.5 keV line or a systematic

<footnote>$^{27}$The significance of the signal found in [613] was communicated in private by T. Jeltema.</footnote>
underprediction of the flux expected from weak plasma lines like K XVIII, claiming that the K XVIII explanation is at odds with theoretical predictions of elemental abundances from Type Ia supernovae. While Tycho does show emission from similar elemental lines as seen in clusters and the GC, supernova remnants have a combination of thermal and non-thermal X-ray emission and the ionization conditions are not necessarily the same as in clusters, making the comparison less than clear cut. However, a 3.5 keV line is not necessarily expected in supernova remnants and the detection of one may point to instrumental or modeling systematics.

Ref. [614] explored the systematic for the Suzaku XIS instrument (somewhat similar to EPIC instruments of the XMM-Newton). They find that systematic errors in the instrumental response of XIS are at the level of half the detected line flux in clusters; they also suggest that the modeling of a large number of plasma emission lines which overlap given the instrumental energy resolution as in the cluster and GC analyses may artificially suppress the continuum in line free regions creating spurious weak signals (this does not explain, however, the origin of the line in the Andromeda galaxy spectrum, where the 3.5 keV line is the strongest one in the range 3-4 keV, see discussion in [610, 624]).

Ref. [611] has developed a new method of data analysis, and applied it to the stacked spectra of galaxies. It uses a non-physical (spline) model of the background. As a narrow line cannot be modeled by a small number of splines, the line gets simulated using the response matrices of the instrument, subtracted from the spectrum and the remaining signal is fitted by the splines. The sensitivity of this method has only been tested on simulated spectra. The simulations show in particular that the sensitivity is the lowest in the range around 3.5 keV (with the exact value of systematic uncertainty not reported by the authors). It can be seen that the simulated line with $sin^2(2\theta) = 7 \times 10^{-12}$ is not recovered at all, while the higher values of the mixing angle are recovered with about $\Delta sin^2(2\theta) \approx 1.6 \times 10^{-11}$ lower value (the negative residual in the spectrum propagated linearly). Thus, the lowest limit from [191] cannot be probed with this method. In addition, the method of [611] constructs the “signal estimator” by dividing an x-ray signal in each pixel by the expected DM signal. The amount of DM in each pixel is roughly estimated by combination of two scaling relations (one is between stellar and halo masses, [628] and another concentration-mass relation [629]) and it is not clear how strong is the systematic uncertainty from wrong DM estimates. Finally, the recovered continuum residuals shown in [611] indicate a large systematic uncertainty in the continuum that fluctuates up to levels of the signal, which should affect the claimed sensitivity.

Ref. [612] has searched for a line in four individual clusters, Perseus, Coma, Virgo, and Ophiuchus is deep Suzaku spectra. While they detect a line near 3.5 keV in Perseus, they do not detect a line in the other three clusters. Ref. [612] claims that this result is inconsistent with a DM interpretation, however, they do not provide sufficient analysis of uncertainties in DM densities of these clusters and therefore it is not clear, how robust this conclusion is. Moreover, the results of ref. [612] appear to be inconsistent with those reported in refs. [614, 627].

The authors of ref. [622] studied central parts of 19 galaxy clusters observed by XMM-Newton having the largest known DM column density in XMM-Newton Field-of-View. In 8 of them (confirming previous detections in Perseus [190] and Coma [612] clusters), the new line has been detected at $> 2\sigma$ level. The average position of the new line in these objects is 3.52±0.08 keV. The decaying DM lifetime imposed from the new line observations and non-observations in these objects is $(3.5 - 6) \times 10^{27}$ s consistent with previous detections of [190, 608, 608]. Recently, ref. [630] did not find the line in a stack of 33 clusters observed with Chandra and Hitomi did not have a significant detection in Perseus [631].
More recently, two papers have argued for potential astrophysical origins in the 3.5 keV line [632, 633]. Ref. [632] analyze high-resolution solar flare spectra and find that the coronal potassium abundance is a factor of 9-11 higher than the photospheric abundances assumed in the above analyzes, offering an explanation as to why the flux of these lines may be brighter than expected. It is, however, not clear why the solar corona should exhibit similar isotopic abundances to the IGM, as the metallicity is generally expected to be higher near stars. Ref. [633] argue that charge exchange, which has thus far not been included in the spectral modeling, can lead to enhanced S XVI lines near 3.46 keV. Though this is in principle possible, the precise energy of the S XVI lines does not agree with the unidentified line, and charge exchange has so far not been detected in x-ray cluster plasma.

As discussed above, an observation of the line in the direction of a dSph satellite of the Milky Way would be a ‘smoking gun’ for the DM interpretation of the signal, due to the absence of X-ray emitting gas in these systems. Ref. [596] has reported a search of the DM decay line in the stacked sample of dSph galaxies. Non-detection of the signal at 3.5 keV in the stacked spectrum of dSph imposes a $2\sigma$ upper limit on the mixing angle at the level $\theta^2 < 1.1 \times 10^{-11}$ and rules out the central value of $\theta^2$ reported in the ref. [190]. The significance of this conclusion is claimed to be $3-4\sigma$, depending on the assumption about the contribution of the Milky Way to the line signal. The analysis in ref. [313], which is claimed to be more constraining, only leads to a $1-2\sigma$ discrepancy.

Most recently, a data from a very long ($\sim 1.4$ Msec) XMM-Newton observation of Draco dSph has become available. A statistically significant emission line from Draco dSph has not been detected [618, 619]. This limits the lifetime of a decaying DM to $\tau > (7-9) \times 10^{27}$ sec at 95% CL the interval corresponds to the uncertainty of the DM column density in the direction of Draco). However, the analysis in [619] found a positive spectral residual (above the carefully modeled continuum) at $E = 3.54 \pm 0.06$ keV with a $2.3\sigma$ significance in the spectrum of the EPIC PN camera. The two MOS cameras show less-significant or no positive deviations, consistently within $1\sigma$ with PN [619]. An analysis of the same data by other authors in [618], found no excess near 3.5 keV in any of the three EPIC instruments. As a result, the significance of the exclusion is currently under debate.

### 5.3.3 Other line candidates in keV range

In addition to the new line at 3.5 keV described in section 5.3.2 above, other potential signs of the radiatively decaying DM have been discussed:

1. According to [604], the intensity of Fe XXVI Ly-γ line at 8.7 keV observed in Suzaku/XIS spectrum of the Milky Way center [605] cannot be explained by standard ionization and recombination processes, and the DM decay may be a possible explanation of this apparent excess.

2. According to section 1.4 of [615], two faint extra line-like excesses at 9.4 and 10.1 keV have been detected in the combined Suzaku/XIS spectrum of Galactic Bulge region. The positions of these excesses do not coincide with any bright astrophysical or instrumental line and their intensities can be explained in frames of decaying DM origin (see right figure 8 of [615]).

### 5.4 Laboratory constraints (author: O. Dragoun)

The β-ray spectroscopy is known to provide model independent upper limits on the effective mass of the electron neutrino, $m^2_{(\nu_e)} = \sum_i j U_{ei} j^2 \cdot m_i^2$, where $U_{ei}$ are the elements of the
<figure><img src="image_21.png" /><figcaption>**Figure 20.** Comparison of sensitivities of current and planned high energy missions towards the searches for weak lines.</figcaption></figure>

neutrino mixing matrix and $m_i$ are the neutrino mass states (see e.g. [634]). However, the measurements turned out to be an uneasy task. The massive neutrino changes the $\beta$-spectrum only it its uppermost part having an extremely low relative intensity. Therefore, the electron spectrometers intended for the neutrino mass determination should exhibit simultaneously the large luminosity $L$ (the product of the relative input solid angle $\Omega$ and the radioactive source area $S$), high energy resolution $\Delta E_{\text{instr}}$ and sufficiently low background. The large source area is needed to keep the sources as thin as possible. Otherwise the electron energy losses within the source deteriorate the measured $\beta$-spectrum shape.

The first limit $m_{\nu_e} < 5$ keV was derived from $\beta$-spectrum of $^{35}$S (with a $\beta$-decay energy of $Q_\beta = 167$ keV) taken with a magnetic spectrometer in 1948 [635]. The measurement of $^{3}$H $\beta$-spectrum ($Q_\beta = 18.6$ keV) with a setup of two orders of magnitude higher luminosity and the magnetic spectrometer set to $\Delta E_{\text{instr}} = 40$ eV yielded in 1972 the limit $m_{\nu_e} < 55$ eV at 90% CL [636]. The present limit is $m_{\nu_e} < 2$ eV [85], while the KATRIN experiment aims at 0.2 eV sensitivity [637].

In 1980, a long time before the undoubted evidence for neutrino oscillations, Shrock examined a possibility of searching for neutrino mass states $m_i$ also in $\beta$-spectra [127]. An admixture of each of such states should produce a specific discontinuity (a kink) in the $\beta$-
spectrum at energy $E_0 - m_i$, where $E_0$ is the endpoint energy of a particular $\beta$-emitter (i.e. the maximum energy of $\beta$-particles in the case when all $m_i = 0$). The relative intensity of the kink observed at $E_0 - m_i$ would determine the value of $jU_{e4}j^2$. The idea stimulated several investigators to search for an admixture of heavier neutrinos in their $\beta$-spectra. The decay of $^{64}$Cu that proceeds via both $\beta^+$- and $\beta^-$-branches with the endpoints energies of 653 and 579 keV, respectively, was explored with a magnetic spectrometer. The overall resolution including contribution from the energy losses within the radioactive source was determined by means of conversion electron spectroscopy (see e.g. [638]) and amounted 1 keV at 300 keV. The result was $jU_{e4}j^2 < 8 \cdot 10^{-3}$ at 90% CL for $110 \text{ keV} \le m_4 \le 450 \text{ keV}$ [639]. Also the investigation of the $^{20}$F $\beta$-decay ($Q_\beta = 7.0 \text{ MeV}$) with a magnetic spectrometer adjusted for the resolution of 3 keV at 4 MeV did not reveal any admixture of sterile neutrinos. The upper limit of $jU_{e4}j^2$ ranged from 0.59 to 0.18 percent for $400 \text{ keV} \le m_4 \le 2.8 \text{ MeV}$ [640].

Simpson investigated the $\beta$-spectrum of tritium implanted into a Si(Li) detector and observed a distortion in the spectrum part below 1.5 keV. He interpreted this distortion as the evidence of a heavy neutrino emission with the mass of about 17.1 keV and mixing probability of 3% [641]. This surprising finding prompted numerous studies of $^3$H and other $\beta$-emitters like $^{14}$C, $^{35}$S, and $^{241}$Pu. For some time, the effect was seen in several $\beta$-spectra taken with semiconductor spectrometers but it has never been recorded by magnetic and electrostatic instruments. Finally, the study of the $\beta$-spectrum of $^{63}$Ni ($Q_\beta = 66.9 \text{ keV}$) with a magnetic spectrometer put the admixture of the 17 keV neutrino below $5 \cdot 10^{-4}$ at 95% CL and found $jU_{e4}j^2 < 1 \cdot 10^{-3}$ for all $m_4$ between 4 and 30 keV [642]. At present, there are no $\beta$-spectroscopic indications for sterile neutrinos [643–647]. The relevant data are summarized in [85]. The best upper limits of $jU_{e4}j^2$ are shown in figure 21. Current investigators of sterile neutrinos may benefit from a detailed discussion of experimental issues concerning the 17 keV neutrino [648]. As an example, the scattering of $\beta$-particles from $^{35}$S decay on a thin diaphragm in front of the Si(Li) detector can produce a false 0.3% admixture of the 17 keV neutrino [649].

The reliable analysis of measured $\beta$-spectra requires precision knowledge of the spectrometer resolution function (the response to a mono-energetic signal). Although this function can be calculated it is advisable to check it for the actual instrument by means of an electron gun. Even better is to apply a suitable conversion electron line, where the natural width and shake up/off satellites were properly taken into account. Most of the systematic errors in $\beta$-spectroscopic searches for active as well as sterile neutrinos originated from incorrect treatment of the electron energy losses within radioactive sources and spectrometers. Corrections of measured $\beta$-spectra having no clear physical meaning not only decrease sensitivity of the search but can lead to erroneous conclusions.

In the study of $^{241}$Pu $\beta$-spectrum ($Q_\beta = 20.8 \text{ keV}$) with two electrostatic spectrometers the electron energy losses within a vacuum evaporated source were incorporated by means of the Monte Carlo simulation of individual elastic and inelastic scattering events. The quality of the approximation was tested by measuring the shape of the 3.7 keV conversion electron line of the $^{169}$Yb source covered subsequently with a thin plutonium layer (see figure 22). This approach allowed describing the measured part of the $^{241}$Pu $\beta$-spectrum (down to 2 keV) without any artificial fitting parameter [657].

Although the kink caused by the keV sterile neutrinos would change the $\beta$-spectrum part of several keV width, it is advisable to apply sufficient spectrometer resolution in order not to miss some unexpected spectrum features. For example the lines at energies of about 260 and 510 eV were observed on $\beta$-spectra of (ground state-to-ground state) decays of $^{241}$Pu
<figure><img src="image_22.png" /><figcaption>**Figure 21.** The best upper limits on the admixture $jU_{e4}j^2$ of sterile neutrinos derived from measured $\beta$-ray spectra. The limits are at the 90% CL except those for $^{64}$Cu [639] and $^{3}$H [650] that are at the 90% CL. The numbers refer to the following $\beta$-emitters and their investigators: 1 → $^{3}$H [651], 2 → $^{3}$H [650], 3 → $^{187}$Re [652], 4 → $^{3}$H [653], 5 → $^{3}$H [654], 6 → $^{63}$Ni [642], 7 → $^{63}$Ni [655], 8 → $^{35}$S [656], 9 → $^{64}$Cu [639], 10 → $^{20}$F [640]. The figure is reproduced from [634].</figcaption></figure>

and $^{63}$Ni that corresponded to the KLL Auger lines of carbon and oxygen [658]. The K-shell vacancies were created by impact of $\beta$-particles emerging from the sources on atoms in the contamination overlayer (see also figure 23). In $\beta$-spectra taken with a worse resolution the Auger lines would smear and could be misinterpreted as an anomaly.
<figure><img src="image_23.png" /><figcaption>**Figure 22.** The K-shell internal conversion electron line of the 63.1 keV transition in <sup>169</sup>Tm measured with an original and modified <sup>169</sup>Yb source [657]. The original source was prepared by vacuum evaporation on Al backing and was contaminated by a usual hydrocarbon overlayer. Subsequently, the source was covered by a thin plutonium layer in which the 3.7 keV electrons suffered additional elastic and inelastic scattering. The agreement of measured (points) and calculated shapes (smooth curves) demonstrates correctness of the electron loss function in plutonium applied in calculations of the <sup>241</sup>Pu β-spectrum shape. The natural width of this conversion line amounting 32 eV was taken into account.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 23.** The KLL-Auger lines of carbon and oxygen superimposed on continuous β-spectra of <sup>241</sup>Pu and <sup>63</sup>Ni [658]. The radioactive sources were prepared by evaporation in vacuum and measured with an electrostatic spectrometer.</figcaption></figure>
# 6 Constraining keV neutrino production mechanisms

Section Editors:
Marco Drewes, Fedor Bezrukov, George Fuller

Neutrino states that are considered “sterile” can (by definition) not feel any of the known forces of nature except gravity. That is, they carry no charges and are singlets under all gauge groups of the Standard Model. In order to compose the observed DM, they must, however, have some interactions with other particles in order to be produced in the early Universe. One can roughly distinguish three different ways how this can be realized.

* The sterile states can be produced via the weak interaction if they mix with ordinary neutrinos. This mixing occurs very generically because there is no a priori reason why the neutrino mass and interaction eigenstates should be aligned in flavor space. If the observed neutrino masses are generated via the (type-I) seesaw mechanism (see section 2.5 for details), then there is definitely mixing between active and sterile neutrinos. Production mechanisms that are based on this mixing are described in sections 6.1 and 6.2.

* Many theories beyond the SM predict the existence of additional scalar fields in Nature. If sterile neutrinos have Yukawa couplings to a scalar field, then they can be produced during its decay in the early Universe. This possibility is explored in section 6.3.

* If the SM is embedded into a more general theory of Nature (such as a Grand Unified Theory), then there are usually additional gauge symmetries that are “broken” at some energy scale above the reach of the LHC. This gives masses of the order of the energy scale where the symmetry is broken to the gauge bosons, effectively switching off the corresponding gauge interaction at lower energies. Thus, the neutrino states that appear to be sterile at the energies which can be probed in experiments may interact via new gauge interactions at higher energies. This allows to produce them in the high temperature phase of the early Universe. This idea is discussed in section 6.4.

## 6.1 Thermal production: overview (authors: M. Drewes, G. Fuller, A.V. Patwardhan)

### 6.1.1 Motivation

All that we know about DM is consistent with the hypothesis that it is composed of long-lived, massive particles without long-range interactions to each other or to SM particles. Neutrinos are the only particles in the SM that fulfil these requirements. They are, however, too light: thermally produced neutrinos $\nu_\alpha$ are relativistic when structures in the Universe begin to form, and contribute to the total DM density as Hot Dark Matter. This is clearly disfavoured by structure formation, e.g., see section 4.1. Sterile neutrinos, on the other hand, can be heavier and “colder”. This makes them natural DM candidates. Their interaction strength with all known particles is “weaker than weak”, i.e., suppressed with respect to the weak interaction by the active-sterile mixing angles $\theta_{I\alpha}$. Though $N_I$ are unstable particles and decay via the weak interaction of their $\nu_L$-component, a sufficiently small mixing angle $\theta_{I\alpha}$ can guarantee the required longevity. The emission line that is expected from the decay $N_I \to \nu_\alpha\gamma$ is one of the main testable predictions of sterile neutrino scenarios, and is discussed in section 5.3.

Whether sterile neutrinos are hot, warm, or cold Dark Matter depends on the way they were produced in the early Universe. The sterile neutrino production scenarios all break
down into two broad categories: (1) those where sterile neutrinos are produced copiously early on, perhaps even attaining thermal and chemical equilibrium like the active neutrinos, but with early decoupling and subsequent dilution of their relic densities down to acceptable DM values; and (2) those where there are no sterile neutrinos present to begin with. In the minimal scenario described by the seesaw Lagrangian, eq. (2.21), $\nu_R$ are gauge singlets, and their only interactions with other particles at low energies are mediated by their mixing $\theta$ with ordinary neutrinos. Then thermal production via this mixing is the only way to produce the mass eigenstates $N_I$. For values of $M_I$ and $\theta_{I\alpha}$ that are consistent with phase space considerations (see section 5.1) and the longevity constraint, the $N_I$ never reach thermal equilibrium in the early Universe. Hence, the minimal scenario falls into category (2). The production by mixing is absolutely unavoidable for any non-zero value of $\theta_{I\alpha}$, therefore the $N_I$ population produced in this way is the minimal Dark Matter content of the Universe for a given parameter choice. If the $N_I$ in addition have other interactions at high energy, then these can produce additional $N_I$ populations. The most studied examples of these are described in detail in section 6.3 and section 6.4.

In this section we concentrate on the minimalist (and unavoidable) model for production of a relic sterile neutrino number density and energy spectrum that could make these particles suitable DM candidates, namely, the thermal production via mixing with ordinary neutrinos [177].$^{28}$ This model is “minimalist” in that it requires only a minimal extension of the Standard Model, to massive active and sterile neutrinos which mix in vacuum, but which otherwise have only Standard Model interactions (i.e., the weak interaction in the case of active flavors). This guarantees that this model falls into the second category, i.e., no sterile neutrinos in the early Universe to begin with.

### 6.1.2 Active-sterile neutrino oscillations

Before going into more details of the thermal production mechanism, we explain the underlying physics in a simple quantum mechanical picture. In the following, we assume that only one species of heavy neutrinos contributes to the Dark Matter density and refer to this mass eigenstate as $\nu_s \equiv N_1$.

The earliest theoretical exploration of neutrino flavor mixing was done by Pontecorvo in the 1960s [78], and it already involved speculation about mixing between active and sterile neutrinos. Let us, for simplicity, assume that there is only one active neutrino flavor $\nu_\alpha$ and a sterile state $\nu_s$, with their mixing in vacuum described by

$$|\nu_\alpha\rangle = \cos\theta |\nu_1\rangle + \sin\theta |\nu_2\rangle, \qquad (6.1)$$

$$|\nu_s\rangle = -\sin\theta |\nu_1\rangle + \cos\theta |\nu_2\rangle, \qquad (6.2)$$

where $|\nu_1\rangle$ and $|\nu_2\rangle$ are the neutrino energy/mass eigenstates (to be identified with $\nu_i$ and $N_I$ defined in 2.5) with corresponding vacuum mass eigenvalues $m_1$ and $m_2$ (to be identified with the light $m_i$ and heavy $M_I$ defined in 2.5). Here $\theta$ is the vacuum mixing angle. Of course, the unitary transformation between the three active neutrinos and any sterile states could be quite complicated, with many parameters, but this simple 2 × 2 example will serve to illustrate the essence of how sterile neutrinos can be produced in the early Universe.

Quantum mechanical systems can be thought of as evolving with time in two ways. A state can evolve coherently in a smooth, continuous fashion according to a Schrödinger-like

<footnote>$^{28}$Note that a minimum amount of thermally produced right handed neutrinos even exists if neutrinos are Dirac particles [659].</footnote>
equation. Or, a state could change abruptly when a “measurement” is made, and the state is “reduced” into one of the eigenstates of the observable in question. The latter process is sometimes associated with de-coherence, because quantum mechanical information in the initial state is lost in the reduction/measurement process. Now we understand, of course, that de-coherence is really a manifestation of mixing and transferring quantum information over many states in the environment, a fundamentally many body process.

Correspondingly, there are two ways that a neutrino can change its flavor. The first way is through coherent neutrino flavor oscillations. As the above transformation equation shows, the neutrino propagation states, i.e., the energy/mass states, are not coincident with the weak interaction eigenstates, i.e., the flavor states. As a neutrino propagates through the plasma of the early Universe it can coherently forward scatter on particles that carry weak charge, for example, electrons and positrons, other charged leptons and their anti-particles, quarks, nucleons, and other neutrinos. In so doing, the neutrino will acquire an effective mass, akin to the way a photon acquires an index of refraction and effective mass as it propagates through glass. Also because of this, the neutrino will typically have an in-medium mixing angle $\theta_m$ that will generally differ from the vacuum mixing angle. At any point along its path a coherently propagating neutrino will have an amplitude to be either active or sterile, which will depend on this mixing angle $\theta_m$, and the oscillation phase.

The second way neutrinos can change their flavor is through direction- or energy-changing scattering processes. These will act like “measurements.” Active neutrinos with energy $E_ν$ will have scattering cross sections $\sigma \sim G_F^2 E_ν^2$, and overall scattering rates $\Gamma_{\nu_\alpha} \sim \sigma \cdot \text{flux} \sim G_F^2 T^5$, where $G_F$ is the Fermi constant that sets the strength of the weak interaction, $E_ν$ is the typical energy of active neutrinos which is of order the temperature $T$, and the number densities (or fluxes) of thermally equilibrated light particles are $\sim T^3$. A neutrino will scatter only through its active component.

In a single particle sense it works like this. A neutrino is set to a weak interaction (flavor) eigenstate when it suffers a scattering event. Subsequently, it propagates coherently toward its next scattering target, building up quantum mechanical flavor phase in a way governed by the in-medium unitary transformation between the flavor states and the instantaneous energy eigenstates:

$$|\nu_\alpha\rangle = \cos \theta_m(t) |\nu_1(t)\rangle + \sin \theta_m(t) |\nu_2(t)\rangle, \qquad (6.3)$$

$$|\nu_s\rangle = -\sin \theta_m(t) |\nu_1(t)\rangle + \cos \theta_m(t) |\nu_2(t)\rangle, \qquad (6.4)$$

where $t$ is the age of the Universe. When this neutrino comes upon the next scattering target, after one mean free path, the target “asks:” What flavor are you? That collapses the neutrino state into one of the flavor eigenstates, either sterile $|\nu_s\rangle$ or active $|\nu_\alpha\rangle$. In the latter case, i.e., wave-function collapse into an active flavor state, the neutrino actually suffers a scattering event. If we start out with all active neutrinos, then as the Universe expands, active neutrino scatterings will necessarily build up a reservoir of sterile neutrinos.

Though there are many models for producing a DM relic density of sterile neutrinos [66, 171, 177, 188, 189, 239, 313–316, 320, 324, 344, 397, 398, 480, 660–685], if active-sterile neutrino flavor mixing exists, then scattering-induced decoherence will necessarily contribute to the relic density [686]. Some of these models run afoul of X-ray observations or large-scale structure considerations or both, as described in refs. [186, 315, 326, 566, 687, 688] and also in section 5.3. However, many of them still remain viable.
### 6.1.3 De-cohering scatterings

At each real or potential scattering event, the probability that the active neutrino collapses into a sterile state is $\propto \sin^2 2\theta_m(t)$. Of course, as Pontecorvo pointed out, a sterile neutrino really is not “sterile,” since this neutrino will also build up flavor phase as it propagates and it can have a non-zero amplitude to be active. Its interaction rate will be $\Gamma_{\nu_s} \sim G_F^2 T^5 \sin^2[2\theta_m]$, and sterile neutrinos can be converted back to active states via scattering-induced decoherence.

As mentioned before, active neutrinos in the early Universe can acquire an effective mass through coherent forward scattering off of particles in the background plasma. This leads to a modification in the effective active-sterile mixing angle. For a neutrino state with momentum $p$, the effective in-medium mixing angle at a plasma temperature $T$ is given by

$$ \sin^2(2\theta_m) = \frac{\Delta^2(p) \sin^2(2\theta)}{\Delta^2(p) \sin^2(2\theta) + [\Delta(p) \cos(2\theta) - V_D - V_T]^2}, \quad (6.5) $$

where $\Delta(p) = \Delta m^2/(2p)$, with $\Delta m^2$ being the appropriate mass-squared splitting in vacuum. $V_D$ and $V_T$ are the finite density and finite temperature potentials, respectively, felt by the active neutrino. The finite density potential, $V_D$, arises as a result of asymmetries (particle vs. antiparticle) between weakly interacting particles, i.e., baryonic and leptonic asymmetries. The finite temperature potential, $V_T$, on the other hand, results from higher-order corrections and neutrino forward scattering off of thermally created particle-antiparticle pairs. It carries a negative sign, and for $T \lesssim M_W$ (the W-boson mass), is given by $V_T = -G_{\text{eff}}^2 p T^4$, where $G_{\text{eff}} \sim \mathcal{O}(G_F)$ can be taken to represent some overall neutrino coupling, summed over the various species of particles in the background plasma.

The negative sign and the strong temperature dependence mean that the thermal potential can cause heavy suppression of the effective in-medium mixing angle at high temperatures. For a keV–MeV mass sterile neutrino, $|V_T| \sim \Delta(p)$ at $T \sim 0.1–1$ GeV. In the limit of negligible lepton number, i.e., $V_D \sim 0$, the ratio of sterile neutrino interaction rate to the Hubble expansion rate at higher temperatures is $\Gamma_{\nu_s}/H \propto T^{-9}$, whereas at lower temperatures, where the vacuum oscillation term dominates, $\Gamma_{\nu_s}/H \propto T^3$. Owing to this high-temperature suppression, sterile neutrinos with vacuum mixing angles smaller than $\sin^2(2\theta) \sim 10^{-6}$ (10 keV/$M_1$) can never attain thermal equilibrium with the plasma, in the absence of additional non-standard-model interactions. Therefore, in such scenarios, their relic density cannot be set via a freeze-out process, unlike ordinary active neutrinos.

Such mixing with active neutrinos can, however, lead to sterile neutrinos being produced in the early Universe via scattering-induced decoherence, as was first pointed out in refs. [689, 690], and later on applied by Dodelson and Widrow [177] to the DM problem. Active neutrino scattering-induced decoherence into sterile states is driven by the active neutrino scattering rate $\Gamma_{\nu_\alpha} \sim G_F^2 T^5$. However, this scattering also gives rise to the active neutrino matter potentials described above, as well as quantum damping, both of which serve to inhibit active-sterile neutrino conversion. The sterile neutrino production rate via this mechanism can also be shown to peak at temperatures of $\sim 0.1–1$ GeV. Details of this mechanism are given in subsequent sections.

### 6.1.4 MSW-effect and resonant conversion

If there are considerable lepton numbers $\mathcal{L}_\alpha$ in the plasma, then $V_D \neq 0$, and ordinary neutrinos can have mass level crossings with sterile neutrinos. For a sterile neutrino mass in
the range of interest, i.e., keV–MeV, this occurs in the same temperature range ~ 0.1–1 GeV where the scattering-induced production peaks. These mass level crossings are essentially Mikheyev-Smirnov-Wolfenstein (MSW) resonances [318, 319]. These resonances, where $sin^2(2\theta_m) \rightarrow 1$, occur where the effective in-medium mass-squared of an active neutrino matches the vacuum mass-squared of the sterile species, $M_1^2$. From eq. (6.5), the MSW resonance condition can be obtained as $\Delta m^2 \cos(2\theta) = 2xT (V_D + V_T)$, or,

$$M_1^2 \approx \frac{4\sqrt{2}\zeta(3)}{\pi^2} G_F \mathcal{L}_\alpha x_{\text{res}} T^4 - 2r_\alpha G_F^2 x_{\text{res}}^2 T^6, \qquad (6.6)$$

where we have used $\Delta m^2 \cos(2\theta) \approx M_1^2$, and where the first and second terms on the right hand side in eq. (6.6) are $2xT V_D$, and $2xT V_T$, respectively. Here the scaled neutrino energy/momentum is $x \equiv p/T$ (the symbol $\epsilon$ for this quantity is also prevalent in the literature). This quantity becomes a co-moving invariant at late times when the neutrinos decouple and freely stream along geodesics — in that regime, the active neutrino temperature parameter will scale like the inverse of the scale factor, as will the momentum magnitude of the neutrino. In eq. (6.6), omitting the contribution from charged leptons and hadrons, cf. eq. (6.24), $\mathcal{L}_\alpha = 2L_\alpha + \sum_{\beta \neq \alpha} L_\beta$, with $\alpha, \beta = e, \mu, \tau$, corresponding to $\nu_e, \nu_\mu, \nu_\tau$, and where we produce the sterile neutrinos through the $\nu_\alpha \rightleftharpoons \nu_s$ mixing channel, and where the net lepton numbers in the individual neutrino flavors are, $L_\alpha \equiv (n_{\nu_\alpha} - n_{\bar{\nu}_\alpha})/n_\gamma$. Here, the photon number density is $n_\gamma = 2\zeta(3)T^3/\pi^2$, and the thermal neutrino number densities $n_{\nu_\alpha} = T^3 F_2(\eta_\alpha)/(2\pi^2)$, with $F_2(\eta_\alpha)$ the relativistic Fermi integral of order 2 and argument $\eta_\alpha$, the $\nu_\alpha$ degeneracy parameter which, in the dilute limit, is related to the corresponding lepton number by $1.46L_\alpha \approx \eta_\alpha$. The thermal potential arises from neutrino scattering on seas of thermally-created particle-antiparticle pairs in the early Universe, and its magnitude therefore depends on the extent to which the corresponding charged lepton flavor is populated, which in turn depends on the temperature. For example, at $T \approx 10$ MeV, the second term in eq. (6.6), corresponding to the thermal potential, has a parameter $r_\alpha$, which is roughly 79 for $\alpha = e$, and roughly 22 for $\alpha = \mu, \tau$ ($r_\mu$ increases to $r_e$ between $T \approx 10$ MeV and $T \approx 200$ MeV, as the themal $\mu^\pm$ seas become increasingly populated). Figure 24 gives examples of possible level crossing scenarios. In this figure we show the temperature dependence of the effective in-medium mass for $\nu_e$ for constant values of potential lepton number $\mathcal{L}_e$. If $\nu_e$'s are converted to $\nu_s$ at these level crossings, then the potential lepton number $\mathcal{L}_e$ will decrease, and the shape of the in-medium effective mass curve will change.

As this figure shows, an active neutrino can experience two MSW resonances as the Universe expands. This is because the first, positive term in eq. (6.6) dominates at lower temperature and scales like $T^4$, whereas the second, negative thermal term scales like $T^6$. Of course, at high enough temperature $T$, this negative term dominates and the potential will decrease with increasing temperature. We can solve eq. (6.6) for the resonant scaled energies of the active neutrinos,

$$x_{\text{res}}(t) = \frac{\sqrt{2}\zeta(3)}{\pi^2 r_{\alpha} G_F} \frac{\mathcal{L}(t)}{T^2(t)} \left( 1 \pm \sqrt{1 - \frac{\pi^4 r_{\alpha} M_1^2}{4\zeta^2(3) \mathcal{L}^2(t) T^2(t)}} \right), \qquad (6.7)$$

where $t$ is the age of the Universe. Clearly, resonances can occur only when the argument of the square root in this equation is greater or equal to zero, implying that the following condition be met:

$$|\mathcal{L}(t)| \cdot T(t) \geq \frac{\pi^2}{2\zeta(3)} \sqrt{r_{\alpha} M_1^2}. \qquad (6.8)$$
<figure><img src="image_25.png" /><figcaption>**Figure 24.** The effective in-medium mass, $\sqrt{2xT(V_D + V_T)}$, for an electron flavor neutrino, $\nu_e$, is shown as a function of temperature in the early Universe. Two values of initial potential lepton number for the $\nu_e \rightleftharpoons \nu_s$ channel are considered, $\mathcal{L}_e = 7 \times 10^{-4}$ (solid, blue curves) and $\mathcal{L}_e = 3 \times 10^{-4}$ (dot-dashed, black curve). Effective mass curves are shown for the labeled values of scaled neutrino energy $x = E_\nu/T$ (also refered to as $\epsilon$ in literature). As an example, the temperature track (dashed, red) of a sterile neutrino with mass 7.1 keV is shown – because this is sterile the track is independent of temperature. MSW resonances occur at mass level crossings, where the effective in-medium mass of the $\nu_e$ matches the sterile neutrino mass.</figcaption></figure>

Whether active neutrinos $\nu_\alpha$ are converted efficiently to sterile neutrinos $\nu_s$ at these resonances depends on several factors. First, scattering-induced de-coherent active-sterile transformation will be proportional to the product of the active neutrino scattering rate and the effective in-medium $\sin^2(2\theta_m)$, the latter being maximal (unity) at resonance. The number of active neutrinos converted to sterile species at resonance will then depend on a comparison of the scattering rate and the “resonance width,” roughly the length of time the effective in-medium mixing angle is maximal, i.e., $\delta t \sim H^{-1} \tan(2\theta)$, where the causal horizon size is $H^{-1} \sim m_{\text{pl}}/T^2$, with $m_{\text{pl}}$ the Planck mass. The point is that gravitation is weak ($m_{\text{pl}}$ large), the expansion rate of the Universe is slow, and the causal horizon is therefore large. In turn, this means that active-sterile conversion can begin efficiently. In coherent MSW language, we would say that the conversion process at resonance is *adiabatic* initially. However, efficient, adiabatic conversion $\nu_\alpha \to \nu_s$ means that the net lepton number $\mathcal{L}$ and, hence, the overall potential is reduced. Rapid reduction of $\mathcal{L}$ renders neutrino evolution through resonances non-adiabatic and, hence, inefficient. As the temperature $T$ drops and lepton number $\mathcal{L}$ is reduced the condition for the existence of resonances in eq. (6.8) will eventually be violated. At low enough temperature, a combination of non-adiabaticity or cessation of resonances, plus a reduced active neutrino scattering rate, imply that production of sterile neutrino relic
<figure><img src="image_26.png" /><figcaption>**Figure 25.** Example for a resonantly produced DM spectrum. The dashed-dotted line represents a momentum distribution that one would expect for $L_\alpha = 0$ (Dodelson-Widrow case), the dashed line is the colder contribution from the resonance. The solid line is the combined spectrum. This illustrative figure is similar to [186]; state of the art spectra can be found in refs. [317, 691].</figcaption></figure>

density ceases. In the end, the relic density of sterile neutrinos, and the relic energy spectrum of these, will then be the result of a sometimes complicated nonlinear interplay of the histories of active neutrino scattering rates and in-medium active-sterile mixing angle $\theta_m(t)$ [668].

Due to the dependence on $x$, different momentum modes go through the resonance at different temperatures. Therefore, resonant production of a sterile neutrino relic density in the right range to be DM leaves the sterile neutrinos with a distorted energy spectrum, with a significant fraction of the sterile neutrino population skewed toward lower $x$ than would be the case for a Fermi-Dirac black body-shaped distribution function. (In the examples in figure 24 note that for a given $L_c$ only lower values of $x$ can be resonant!) Therefore, resonantly produced sterile neutrino DM is generally “colder” than sterile neutrino DM particles with the same mass produced non-resonantly via a straight Dodelson-Widrow scenario, e.g., see figure 25.

To give some scale for the parameters involved, accounting for all of the DM with a 7.1 keV mass sterile neutrino, with a vacuum mixing angle $\sin^2\theta \sim 10^{-10}$ or $10^{-11}$ as suggested by X-ray observations [190, 191], via resonant, Shi-Fuller mechanism production, would require an initial net lepton number $\mathcal{L} \sim 10^{-5}-10^{-4}$ [315, 316, 480]. A non-resonant, zero-lepton number Dodelson-Widrow production scenario would need a much larger vacuum mixing angle at a given mass to produce all of the DM, and this would mean a radiative decay rate in conflict with X-ray bounds. Such a scenario could work, however, if the 7.1 keV sterile
neutrino were only ~15% of the total DM, thereby allowing for larger active-sterile mixing angles to not be excluded.

Finally, note that a lepton number $\mathcal{L} \sim 10^{-4}$, though well below what can be probed with light element abundances, i.e., $\mathcal{L} \sim 0.1$, is nevertheless very large compared with the measured baryon number $\eta_B \approx 6.11 \times 10^{-10}$. This brings up a very thorny issue. Where would such a “large” lepton number come from? We have problems in explaining the origin of the baryon number as well. However, the existence of several sterile neutrino generations, the non-equilibrium nature of active-sterile neutrino mixing, and CP-violation in the neutrino sector, can produce the Sakharov conditions required to produce lepton numbers and baryon numbers. Moreover, our ignorance of the sterile neutrino mass scale can be turned on its head, and (large) masses of two of these particles chosen so that the baryon and lepton numbers can be produced. The lepton number so produced could then drive resonant production of a lighter, DM candidate sterile neutrino. One such scenario is the $\nu$MSM, where the required $\mathcal{L}$ can be produced in the decay of heavier sterile neutrinos, which previously generated the baryon asymmetry of the Universe [188, 397, 398].

## 6.2 Thermal production: state of the art (authors: M. Drewes, M. Laine)

### 6.2.1 Examples of complete frameworks

The physical picture outlined in section 6.1 is sufficient for understanding most of the effects playing a role, but does as such not yet offer for a quantitative method for determining the DM abundance and spectrum in the early Universe. Indeed the quantities introduced, the potentials $V_T$ and $V_D$ and the active neutrino interaction rate $\Gamma_{\nu_\alpha}$, need to be put together in a way which accounts for contributions both from active neutrinos and antineutrinos; correctly incorporates Pauli blocking factors; tracks spin dependent effects if present; permits for the inclusion of hadronic contributions that are very important at the temperatures $T \gtrsim 100$ MeV where the production peaks; and systematically accounts for the dilution of the various lepton asymmetries during the resonant production of DM particles.

Theoretically speaking, the computation needed in order to complete this task is a problem in non-equilibrium statistical physics. In a typical approach that can be applied if most degrees of freedom are in local thermal equilibrium, one attempts to factorize such problems into two parts. The first part is the postulation of equations of motion that should capture the time evolution of the non-equilibrium degrees of freedom. The second part is the determination of the coefficients appearing in these equations (such as $V_T$, $V_D$, $\Gamma_{\nu_\alpha}$), which is normally carried out in thermal equilibrium. With both ingredients as well as initial conditions in place, the non-equilibrium equations of motion can be integrated in order to obtain the desired abundance.

One example of a non-equilibrium framework, often used for freeze-out computations in cosmology, is that of Boltzmann equations. In this case the dynamical degrees of freedom are classical phase space distribution functions, and the coefficients characterizing the equations are vacuum scattering amplitudes. However, since quantum coherence between states of different flavors is crucial to neutrino oscillations, and since the presence of a medium can significantly change the dispersion relations of various particles even in the absence of “real” scatterings, the classical distribution functions in the Boltzmann equations are not necessarily suitable for the present problem.

Flavor oscillations, thermal modifications of dispersion relations, and real scatterings responsible for de-coherence and the sterile neutrino production, can be taken into account
through matrix-valued generalizations of the Boltzmann equations called density matrix equa-
tions [295]. In this approach, active and sterile neutrinos are thought of as on-shell degrees
of freedom, whose non-equilibrium time evolution is governed by a coupled set of linear
equations of motion. The coefficients parametrizing these equations characterize three ba-
sic properties of neutrinos: the mixing angles; the potentials V<sub>T</sub> and V<sub>D</sub>; and a damping
coefficient, reflecting thermal scatterings of left-handed neutrinos with the Standard Model
particles constituting the heat bath.

The basic quantities in the density matrix framework are expectation values of the
ladder operators a<sub>I</sub>, a<sub>I</sub><sup>†</sup> that appear in the plane wave expansion

$$
N_I = \sum_h \int \frac{d^3 \mathbf{p}}{(2\pi)^3} \frac{1}{2\sqrt{\mathbf{p}^2 + M_I^2}} \left( u_{I,\mathbf{p}}^h e^{-i\mathbf{p}\mathbf{x}} a_{I,h}(\mathbf{p},\mathbf{t}) + v_{I,\mathbf{p}}^h e^{i\mathbf{p}\mathbf{x}} a_{I,h}^{\dagger}(\mathbf{p},\mathbf{t}) \right),
$$

$$
\nu_i = \sum_h \int \frac{d^3 \mathbf{p}}{(2\pi)^3} \frac{1}{2\sqrt{\mathbf{p}^2 + m_i^2}} \left( u_{i,\mathbf{p}}^h e^{-i\mathbf{p}\mathbf{x}} b_{i,h}(\mathbf{p},\mathbf{t}) + v_{i,\mathbf{p}}^h e^{i\mathbf{p}\mathbf{x}} b_{i,h}^{\dagger}(\mathbf{p},\mathbf{t}) \right).
$$

Here $h, h' = \pm$ label the helicities of the particles. The elements of the density matrix can
be defined as

$$
(\rho_{NN})_{IJ}^{hh'} \propto \frac{\langle a_I^h(p)^\dagger a_J^{h'}(p) \rangle}{V_s}, \quad (\rho_{\nu\nu})_{ij}^{hh'} \propto \frac{\langle b_i^h(p)^\dagger b_j^{h'}(p) \rangle}{V_s}, \quad (\rho_{\nu N})_{iJ}^{hh'} \propto \frac{\langle b_i^h(p)^\dagger a_J^{h'}(p) \rangle}{V_s}. \quad (6.9)
$$

Here, V is some arbitrarily chosen unit volume and s the entropy density of the Universe. For I = J or i = j these objects simply give the occupation number of particles of flavor I in the quantum state in which the expectation value is calculated, the off-diagonal elements characterize correlations between different flavors. Together, they form a “density matrix” in flavor space for each momentum mode and helicity,

$$
\rho = \begin{pmatrix} \rho_{\nu\nu} & \rho_{\nu N} \\ \rho_{N\nu} & \rho_{NN} \end{pmatrix}. \qquad (6.10)
$$

Equations of this type have been widely used to describe the propagation of active neutrinos
in a medium, particularly at low temperatures $T \lesssim 10$ MeV.

Let us for a moment assume that the neutrinos reside in a thermal bath with adiabatically changing temperature. It is usually also assumed that one can neglect all correlations with $h \neq h'$ because the Hamiltonian commutes with the helicity operator. Then the density matrix follows an equation of motion of the form

$$
\frac{1}{\mathcal{H}X} \frac{d}{dX} \rho = -i[H_{\text{eff}}, \rho] - \frac{1}{2}\{\Gamma_{\text{eff}}, \rho - \rho^{\text{eq}}\}. \quad (6.11)
$$

Here, $\rho^{\text{eq}}$ is the equilibrium density matrix and $X = M/T$ a dimensionless time variable, where $M$ is an arbitrarily chosen mass scale. It is convenient to identify $M$ with the DM particle's mass. The function

$$
\mathcal{H} \equiv - \frac{\partial}{\partial X} \sqrt{\frac{45}{4\pi^3 g} \frac{m_P}{2M^2} X} \qquad (6.12)
$$

can be identified with the Hubble parameter if the number of degrees of freedom *g* is constant
during the evolution. The commutator contains an effective Hamiltonian *H*<sub>eff</sub> that describes
coherent flavor oscillations. In vacuum, this term is simply given by the neutrino mass matrix. Above we have defined $\rho$ in the “mass basis” in flavor space; in this basis $H_{\text{eff}}$ is diagonal in vacuum. The equation of motion (6.11) is invariant under rotations in flavor space and can be transformed into any other basis for convenience.

The anticommutator comes from the interactions and gives rise to “dissipation”, i.e. processes that change the number of particles of a given species in a given mode. It involves decoherence and leads to scattering-induced sterile neutrino production. The flavor basis in which $\Gamma_{\text{eff}}$ is diagonal is called “interaction basis”, in general it is not identical to the one where $H_{\text{eff}}$ is diagonal. It is precisely this misalignment that leads to flavor oscillations: particles are produced in a state that corresponds to a $\rho$ which is diagonal in the interaction basis, then they start oscillating due to the presence of $H_{\text{eff}}$. At finite temperature and density, there are corrections to $H_{\text{eff}}$ due to $V_D$ and $V_T$. Let us split $H_{\text{eff}} = H_{\text{vac}} + H_T$ into a vacuum term and a finite temperature/density correction. $H_{\text{vac}}$ is determined by the vacuum mass matrix and diagonal in the mass basis, $H_T$ comes from the interactions and is diagonal in the interaction basis. The total $H_{\text{eff}}$ defines the effective (thermal) mass basis at a given temperature and density, it rotates when the temperature changes because the relative size of $H_{\text{vac}}$ and $H_T$ changes. In particular, in the limit of very high temperatures $T \gg M_I$, $H_{\text{vac}} \ll H_T$ can be neglected, and the effective mass basis becomes identical to the flavor basis. Then there are no flavor oscillations in spite of the fact that neutrinos have large “thermal masses”, simply because these thermal mass terms are diagonal in the flavor basis.

In practice, one usually makes a number of simplifications. Since active neutrino masses are very small, oscillations among active neutrinos can usually be neglected at the temperatures that are relevant for DM production. This allows to set $m_\nu = 0$, and $H_{\text{eff}}$ takes the form

$$H_{\text{eff}} = \begin{pmatrix} 0 & 0 \\ 0 & \frac{M_I^2}{2|\mathbf{p}|} \end{pmatrix} + H_T. \qquad (6.13)$$

Since a rotation among the active flavors does not change $H_{\text{vac}}$, it is common to use the flavor basis $(\rho_{\nu\nu})_{\alpha\beta}^{hh}$, where $H_T$ is diagonal, in the active sector. This is done by a block-diagonal rotation matrix diag$(U_\nu, 1)$. We can neglect all elements with $\alpha \neq \beta$ in this approximation. The active neutrinos are highly relativistic. This allows to interpret the helicity states of the Majorana field $\nu_i$ as “particles” and “antiparticles”. Their sum is the total particle number and their difference is a measure for the lepton asymmetry in the plasma (neglecting the asymmetry carried by charged leptons). If the active neutrinos are in thermal equilibrium at $T > 1.1$ MeV, the deviation of their phase space distribution function from equilibrium can be neglected, and the occupation numbers for all modes are simply given by a single number $T$. In kinetic equilibrium, one in addition needs to specify the asymmetries $L_\alpha$. They are crucial for sterile neutrino DM production because they lead to a non-zero $V_D$ and can trigger the MSW resonance. Since the active neutrino distribution functions remain close to equilibrium, one need not track the occupation numbers of individual momentum modes and can define the integrated quantities:

$$n_{\nu_{\alpha}} = \int \frac{d^3 \mathbf{p}}{(2\pi)^3} [(\rho_{\nu\nu})_{\alpha\alpha}^{++} + (\rho_{\nu\nu})_{\alpha\alpha}^{--}], \qquad (6.14)$$

$$L_{\alpha} = \int \frac{d^3 \mathbf{p}}{(2\pi)^3} [(\rho_{\nu\nu})_{\alpha\alpha}^{++} - (\rho_{\nu\nu})_{\alpha\alpha}^{--}]. \qquad (6.15)$$

We are interested in the production of sterile neutrinos with masses $M_I$ in the keV range and mixing angles $U_{I\alpha}^2 < 10^{-8}$ that are consistent with searches for emission from decaying
DM. Comparing a keV mass to the observed $\Omega_{DM}$ shows that the occupation numbers $\rho_{NN}$ at all times must remain significantly below their equilibrium values. The same can be said about $\rho_{\nu N}$ and $\rho_{N\nu} = \rho_{\nu N}^\dagger$. On the other hand, the equilibrium density matrix $\rho^{eq}$ can in good approximation be replaced by a unit matrix as long as the temperature is much larger than all particle masses, which is the case in the regime $T \sim 100$ MeV where thermal sterile neutrino DM production peaks. If $L_\alpha = 0$, then the off-diagonal elements $\rho_{\nu N}$ oscillate very rapidly due to the very different effective masses of active and sterile neutrinos. Then their effect on the evolution of the abundances $\rho_{NN}$ averages out. For $L_\alpha \neq 0$ the nonzero $V_D$ can lead to a level crossing between the active and sterile neutrino quasiparticle dispersion relations in the plasma, so this argument does not hold. However, in this case the element of $H_{\text{eff}}$ that multiplies an element of $\rho_{\nu N}$ on the r.h.s. of (6.11) becomes very small (because it is given by the effective mass splitting). That is, the elements of $\rho_{\nu N}$ (which are small numbers) enter the r.h.s. of (6.11) only as in a product with another small number, and this product can be neglected.<fnref n="29" />

Even though the density matrix formalism can describe much of the physics relevant to the problem, it still omits uncertainties of order unity. In particular, it assumes that active lepton asymmetries are only carried by neutrinos, omitting the fact that weak interactions are in equilibrium and lepton asymmetries can also be carried by charged leptons (with charge neutrality taken care of by leptons of other generations and/or by light hadrons such as pions).

A framework which includes all effects in Standard Model interactions, with the price of restricting to order $\theta_{I\alpha}^2$ in the active-sterile mixing angles, has been presented in ref. [692] and generalized to systems close to equilibrium and with three generations of lepton asymmetries in ref. [691].<fnref n="30" /> In this approximation, we can track the full density matrix of the system, denoted by $\hat{\rho}$. Its time evolution is determined by the quantum kinetic equation

$$i \frac{d\hat{\rho}_I(t)}{dt} = [\hat{H}_I(t), \hat{\rho}_I(t)], \qquad (6.16)$$

which is equivalent to the Schrödinger equation. Here $\hat{\rho}_I \equiv \exp(i\hat{H}_0t)\hat{\rho}\exp(-i\hat{H}_0t)$ is the density matrix in the interaction picture and $\hat{H}_I \equiv \exp(i\hat{H}_0t)\hat{H}_{\text{int}}\exp(-i\hat{H}_0t)$, where $\hat{H}_{\text{int}}$ is the interaction Hamiltonian proportional to the neutrino Yukawa couplings (i.e. $\theta_{I\alpha}$). The “free” part $\hat{H}_0$ includes all other interactions. Equation (6.16) is exact and contains all interactions.<fnref n="31" /> All other kinetic equations that are being used in the literature correspond to different approximations to (6.16). For instance, equation (6.11) can be interpreted as the limit where the sterile neutrinos are on-shell and the r.h.s. has been perturbatively expanded in $\theta$. This approximation allows to include coherent quantum oscillations and corrections to the (quasi)particle dispersion relations from interactions with the plasma. Equation (6.29) can then be viewed as a limit of (6.11) in which coherent oscillations between different flavors have been neglected.

Equation (6.16) can be solved perturbatively,

$$\hat{\rho}_I(t) = \hat{\rho}_0 - i \int_0^t dt' [\hat{H}_I(t'), \hat{\rho}_0] - \int_0^t dt' \int_0^{t'} dt'' [\hat{H}_I(t'), [\hat{H}_I(t''), \hat{\rho}_0]] + \dots, \qquad (6.17)$$

<footnote n="29">²⁹Here “smallness” can be quantified in powers of the mixing angle $\theta$: the production rate $\Gamma_N$ for sterile neutrinos, which is the lower right block of the matrix $\Gamma_{\text{eff}}$ in (6.11), is of order $\propto \theta^2$. This sets the macroscopic time scale $1/\Gamma_N$ in the system. Oscillations can be considered “fast” if they are faster than this. If the oscillation time scale is comparable to $\Gamma_N$, then the oscillations in principle have to be tracked. In this case, however, the element of $H_{\text{eff}}$ that multiplies an element of $\rho_{\nu N}$ is also $\sim \theta^2$.</footnote>

<footnote n="30">³⁰A different approach based on similar physical assumptions was proposed and used to compute DM spectra in ref. [317].</footnote>

<footnote n="31">³¹Note that the derivative on the l.h.s. should be a covariant derivative in curved space.</footnote>
where $\hat{\rho}_0 \equiv \hat{\rho}(0) = \hat{\rho}_1(0)$. This is equivalent to standard time-dependent perturbation theory. The strictly perturbative approach breaks down as soon as $\rho_{NN}$ approaches equilibrium, but by an appropriate choice of $\hat{\rho}_0$ these so-called secular terms can be avoided, so that the system can be described both near and far from equilibrium, with the former case properly including Pauli blocking. Subsequently the sterile neutrino production rate can be traced out from the solution. Simultaneously, the lepton number operator evolves as [693]:

$$ \hat{L}_{\alpha}(t) = \int_{\mathbf{x}} i \hat{N}_I F_{I\alpha}^{\dagger} \hat{\phi}^{\dagger} \hat{\ell}_{\alpha} + \text{h.c.}, \qquad (6.18) $$

where $\tilde{\phi} \equiv i\sigma_2\phi^*$ is the conjugate Higgs doublet, $\hat{\ell}_{\alpha}$ is the left-handed lepton doublet, and $F_{I\alpha} \propto \theta_{I\alpha}$ are Yukawa couplings. Tracing this equation with the density matrix yields an equation for the depletion of lepton asymmetries. Taken together the equations results in a coupled system, parametrized by $V_T$, $V_D$, and $\Gamma_{\nu_{\alpha}}$, which is valid both near and far from equilibrium to $O(\theta_{I\alpha}^2)$ but in principle to any order in SM couplings such as $\alpha_s$ [691].

In principle there are also frameworks which can be extended to higher orders in $\theta_{I\alpha}^2$. This could be relevant, for instance, for crosschecking that resonant production is treated fully consistently (cf. footnote 29). One possible starting point for this is the Schwinger-Keldysh formalism of non-equilibrium quantum field theory [694–697], in which all properties of the system can be expressed in terms of correlation functions of the quantum fields, without reference to particles or asymptotic states. This framework allows for a first principles description of sterile neutrinos in a dense plasma, and significant progress has been made towards this goal [401, 685, 698–714], in particular in the context of leptogenesis. However, even though exact equations can be written down, it remains a challenge to develop a systematic approximation scheme for their practical solution.

### 6.2.2 Matter potentials and active neutrino interaction rate

As has been discussed in the previous sections, a key ingredient in the DM production computation is the determination of the potentials $V_T$ and $V_D$ and the active neutrino damping coefficient, or interaction rate, denoted by $\Gamma_{\nu_{\alpha}}$. We start with a discussion of how these are normally defined in a way which permits for their practical determination, and then discuss the best available values that can be extracted from literature.

Let us denote by $P \equiv (E, \mathbf{p})$ the four-momentum of a left-handed (active) neutrino and by $\Sigma$ its self-energy. The matter potentials $V_T$ and $V_D$ parametrize the dispersion relation of the active neutrino. This can be obtained by searching for non-trivial solutions of [715–718]:

$$ \{\not P - \text{Re}\,\not\!D P_L\} \psi(\mathbf{p}) = 0, \qquad (6.19) $$

where $\psi$ is a Dirac spinor, $P_L \equiv (1-\gamma_5)/2$ is a chiral projector, and Re denotes the continuous (non-cut) part of the self-energy. Following ref. [719] and for future reference, we have factorized the chiral projector $P_L$ explicitly in eq. (6.19).

Because only left-handed neutrinos experience weak interactions within the Standard Model (the vertices that they attach to are proportional to $\gamma^\mu P_L$), the self-energy must respect chirality. Therefore it is of the form $\not D = a \not P + b \psi$ [719], guaranteeing that the left and right chiral components decouple from each other. The part $a$ can in general be omitted in comparison with the tree-level term in eq. (6.19), whereas the function $b$ defines the quantities that we are interested in.
More precisely, we can write [715]

$$ \mathcal{Z}P_L \equiv \left[ a \not P + \left( b_r - i \frac{\Gamma_{\nu_\alpha}}{2} \right) \psi \right] P_L, \qquad (6.20) $$

where $u \equiv (1, 0)$ is the plasma four-velocity; $a = a_r + ia_i$ is a complex function; and $b_r$ is a real function. The real part can be interpreted as

$$ b_r = -V_T - V_D, \qquad (6.21) $$

where $V_D$ is defined to be the part which is odd in chemical potentials. Moreover, assuming the structure in eq. (6.20), the scattering rate can be projected out as

$$ \Gamma_{\nu_{\alpha}} = -\mathrm{Tr} \left\{ \frac{E \not P - P^2 \psi}{p^2} \mathrm{Im} \mathcal{Z} P_L \right\}. \qquad (6.22) $$

According to ref. [720], this formula corresponds to the classic results given in refs. [720–722].

The scattering rate has many interpretations. On one hand it is often interpreted as the relaxation rate that drives the active neutrino population to thermal equilibrium [723]; more importantly for us, as discussed above, it directly determines the DM production rate,

$$ \Gamma_{\nu_s} \approx \frac{\sin^2(2\theta)}{4} \frac{M_I^4 \Gamma_{\nu_\alpha}}{(M_I^2 + 2E b_r + b_r^2)^2 + E^2 \Gamma_{\nu_\alpha}^2} + \mathcal{O}(\theta^4). \qquad (6.23) $$

This result can be compared to the expression (6.5) for the production in the quantum mechanical limit. General discussions concerning other interpretations of $\Gamma_{\nu_\alpha}$ can be found in refs. [706, 724].

Turning to numerical evaluations, the finite temperature part $V_T$ originates at 1-loop level and has a relatively simple expression, available both for $T \ll M_W$ [715] and more generally [716–718]. In contrast the determination of the finite density potential $V_D$, also addressed in ref. [715], turns out to be subtle. Because of $Z^0$-boson exchange, which is in thermal equilibrium in the temperature range of interest, it includes a contribution from hadrons; in a notation similar to eq. (6.7),

$$ \mathcal{L}_{\alpha} = 2L_{\nu_{\alpha}} + \sum_{\beta \neq \alpha} L_{\nu_{\beta}} + \left(\frac{1}{2} + 2\sin^2\theta_w\right) L_{e_{\alpha}} - \left(\frac{1}{2} - 2\sin^2\theta_w\right) \sum_{\beta \neq \alpha} L_{e_{\beta}} \\ + \left(\frac{1}{2} - \frac{4}{3}\sin^2\theta_w\right) \sum_{i=u,c} L_i - \left(\frac{1}{2} - \frac{2}{3}\sin^2\theta_w\right) \sum_{i=d,s,b} L_i, \qquad (6.24) $$

where $e_\alpha$ denotes charged leptons; $u, c, d, s, b$ stand for quark flavors; and $\theta_w$ is the weak mixing angle. Because the hadronic part couples differently to up and down-type quarks, the hadronic contribution contains a part which is not proportional to the baryon number density. Electric charge neutrality of the Standard Model plasma requires that the hadronic part is present even if the baryon asymmetry of the plasma were zero, because charged hadrons are needed to neutralize the contribution from charged leptons. The way to systematically include these effects, relating in particular the hadronic contributions to so-called quark number susceptibilities which can be measured with methods of lattice QCD, has been worked out in ref. [691], and lattice results have been included in the estimate of $V_D$ in ref. [317].
Finally, a precise determination of $\Gamma_{\nu_{\alpha}}$ is considerably more complicated than that of $V_T$ and $V_D$, including a large number of 2 $\leftrightarrow$ 2 scatterings at $T \ll M_W$. In order to present numerical results, let us define the projections

$$I_P \equiv -\mathrm{Tr} \left\{ \not p \operatorname{Im} \not P_L \right\}_{E=\sqrt{p^2+M_I^2}}, \quad I_u \equiv -\mathrm{Tr} \left\{ \not \psi \operatorname{Im} \not P_L \right\}_{E=\sqrt{p^2+M_I^2}}. \qquad (6.25)$$

For dimensionless combinations, it is natural to factor out the combination $G_F^2 T^4$ (here $G_F$ denotes the Fermi constant) as well as appropriate powers of $E$:

$$\hat{I}_P \equiv \frac{I_P}{G_F^2 T^4 E^2}, \qquad \hat{I}_u \equiv \frac{I_u}{G_F^2 T^4 E}. \qquad (6.26)$$

Then

$$\Gamma_{\nu_{\alpha}} = G_F^2 T^4 E \left\{ \hat{I}_P + \frac{M_I^2}{p^2} (\hat{I}_P - \hat{I}_u) \right\}_{E=\sqrt{p^2+M_I^2}}. \qquad (6.27)$$

Values for $\hat{I}_P$ and $\hat{I}_u$, as a function of $M_I$, $p$ and the active neutrino flavor, can be found on the web site http://www.laine.itp.unibe.ch/neutrino-rate/, together with an illustration of the magnitude of hadronic effects according to the prescription in ref. [320]. It should be noted that for purely leptonic effects at low temperatures $T \lesssim 20$ MeV, these results are $\sim 5\%$ smaller than previous expressions in the literature [715], possibly due to the fact that Pauli blocking factors were omitted in ref. [715].

### 6.2.3 Open questions

Even though the available expressions for $V_T$, $V_D$, and $\Gamma_{\nu_{\alpha}}$ as explained in section 6.2.2 and one of the frameworks discussed in section 6.2.1 permit to compute DM abundances as a function of $M_I$ and $\sin^2(2\theta)$, it may be questioned what the theoretical accuracy of such computations is. Let us list a number of issues to be kept in mind on this point:

(i) First and foremost, there is the question of how precisely hadronic contributions have been included. Thermal production peaks at temperatures $T \gtrsim 100$ MeV, where the Standard Model plasma undergoes a smooth crossover from a confined to a deconfined medium. The effective number of degrees of freedom changes by a factor $\sim 3$, which illustrates that hadronic effects dominate the plasma behavior. Roughly speaking, hadronic effects enter the DM production computation at three points: in determining the thermal history of the Universe through the temperature dependence of the equation of state (i.e. the temperature dependence of the Hubble rate); through their effect on the finite density potential $V_D$; and through the thermal scatterings that active neutrinos undergo with hadronic plasma constituents. In the most advanced current computations [317, 691], hadronic effects on the equation of state have been modelled according to ref. [725]. Lattice results concerning quark number susceptibilities have been included in the estimate of $V_D$ in ref. [317]. The most difficult issue is including hadronic thermal scatterings; even though they can be related to well-defined mesonic spectral functions [692], it is very hard to extract such spectral functions from first principles lattice measurements. Chiral effective theory inspired model computations have been attempted in refs. [317, 726], whereas in ref. [691] a smooth interpolation between confined and deconfined phases, tracking the hadronic contribution to entropy density, was adopted. In total, the hadronic uncertainties from all these ingredients probably remain on the $\sim 20\%$ level or so.
(ii) At temperatures above 10 MeV, where sterile neutrino DM production takes place, the lepton asymmetries carried by the different active lepton generations do not equilibrate with each other. That is, three different lepton asymmetries need to be tracked. Depending on the flavor structure of the neutrino Yukawa couplings, all or only some of these may contribute to resonant DM production. A general framework for describing this rather complicated system, including a proper account of the non-linear backreaction whereby the lepton asymmetries get depleted during resonant DM production, has been developed in ref. [691]. However the results depend significantly (~100%) on the flavor structures of the lepton asymmetries and Yukawa couplings (a collection of spectra for different choices can be found in figure 6 of ref. [691]). In order to get a unique answer, these parameters would need to be fixed from other considerations, such as leptogenesis and phenomenological constraints.

(iii) Although possible in principle [188, 395, 398], a theoretically reliable computation of the lepton asymmetries generated through the non-DM sterile neutrino oscillations is far from being established in practice. We recall that, since lepton asymmetries much larger than the baryon asymmetry are required, the lepton asymmetries need to be generated below a temperature of 160 GeV, so that the Higgs mechanism is active and sphaleron processes have switched off. Carrying out the computation in this setting is complicated, because the system has multiple mass scales and because two almost degenerate sterile neutrino generations participate in the dynamics.

(iv) There are also a number more refined “theoretical issues” that may be worth considering. For instance, one may worry about the gauge dependence of the quantities introduced ($V_T, V_D, \Gamma_{\nu_\alpha}$), which is not automatically guaranteed, given that the self-energy is not needed at the active neutrino on-shell point but rather the sterile neutrino one. Or, one may worry about an additional particle branch in the active neutrino dispersion relation, called the plasmino or abnormal or hole branch [727]. That it appears for active neutrinos in the Higgs phase has been demonstrated in ref. [718], so one may wonder whether the corresponding quasiparticles should be included as states in the density matrix formalism. Finally, in the context of leptogenesis, where sterile neutrino production is considered at high temperatures, it has been realized that to obtain correct leading-order results requires an all-orders resummation of the loop expansion, in order to account for the so-called Landau-Pomeranchuk-Migdal (LPM) effect [728, 729], but so far this has not yet been included in any Higgs phase computations.

We remark that solving the theoretical issues of point (iv) essentially requires giving up the simple picture of sterile neutrino production only through active-sterile oscillations, and with it also frameworks such as the density matrix one. Indeed the way that gauge independence is restored and the need for LPM resummation makes its appearance is that sterile neutrinos can also be produced via direct $1 \leftrightarrow 2$ scatterings involving Higgs bosons or longitudinal polarizations of gauge bosons (Goldstone modes). Such processes become important once $T \gtrsim$ a few GeV, and in general their inclusion requires a fully quantum field theoretic framework.

## 6.3 Production by particle decays (authors: F. Bezrukov, A. Merle, M. Totzauer)

Sterile neutrino appearance by neutrino flavor transitions, as considered in the previous subsection, is not the only mechanism by which sterile neutrino DM could have been produced
in the early Universe. A cosmological population of these particles can be generated in decays of some heavy particles as well. A simple example is given by a singlet (with respect to the Standard Model) scalar S that can decay into sterile neutrinos. A generic interaction Lagrangian is given by

$$ \mathcal{L}_{\text{int}} = \frac{y}{2} \overline{(\nu_R)^c} \nu_R S + \text{h.c.}, \qquad (6.28) $$

where sterile flavor indices are not shown explicitly, but can be added, if necessary. The general idea is fairly simple: S-particles that already exist or are created in some process could decay into sterile neutrinos at some point during the cosmological evolution. Now, the amount of neutrino produced does in general not depend on the active-sterile mixing angle $\theta$, but is instead controlled by the coupling $y$. Furthermore, for a given mass, the sterile DM can be colder or warmer, depending upon details of the particular model used. Therefore, the parameter range which is unrestricted by various observational constraints widens up in the mass-mixing angle parameter space.

There are several variations of the described mechanism, which will be described in more detail below. In particular, the “parent” S-particle can itself be in or out of thermal equilibrium at the time the sterile neutrinos are produced and, depending on the model, the decay into the sterile neutrinos may be the only possible decay channel [730] or one of many possible decay channels into light particles [664, 670, 731]. Also, the parent particle can be different from a scalar, e.g. a vector [732, 733] or fermion [679]. In addition, if the sterile-active neutrino mixing is present, then one has to take into account the mechanisms described in section 6.1 as a competing production mode [686].

### 6.3.1 Decay in thermal equilibrium

Let us first analyze what happens if the parent particle S is in thermal equilibrium for all relevant temperatures. That can happen if the S-particle is strongly mixed with the neutral Higgs boson, so that all decay channels in SM particles are present and maintain it in thermal equilibrium until temperatures below its mass.

The scalar S can be either real or complex. If S acquires non-zero vacuum expectation value (VEV), then the interaction eq. (6.28) leads not only to the neutrino production, but generates also the Majorana masses for them, $M_1 = y \langle S \rangle$. In this case, and if S is complex, the model contains (pseudo) Nambu-Goldstone boson, called the Majoron [203, 734]. As a result, the sterile neutrinos may be not stable enough to form the DM because they decay into a Majoron, if it is lighter,<fnref n="32" /> and an active neutrino. While this case is not interesting from the point of view of sterile neutrinos being the DM, the case of decaying sterile neutrinos cannot a priori be excluded and by this motivates their laboratory searches in the region of comparatively large mixing angles where they would be considered overabundant otherwise. For charged scalars, there can be significant thermal corrections to the DM production rate [685]. These affect both, the equation of motion for the scalar condensate [737] as well as the properties and kinematics of quasiparticles in the plasma [738, 739]. Moreover, the gauge interactions open up additional channels for DM production in $2 \to 2$ scatterings in addition to the decay [685].

If the S-field is real, in turn, the Majoron does not exist. A natural example of such a field S in eq. (6.28) would be the inflaton $\phi$. While nothing forbids the inflaton to have a VEV to generate Majorana masses for neutrinos (or even to be complex), in popular “mainstream”

<footnote n="32">³²The Majoron is not necessarily massless, and may in turn form DM, for the recent discussion see e.g. refs. [735, 736].</footnote>
models the inflaton is a real scalar field with zero VEV. The inflaton field $\phi$ is the driving force behind the Big Bang, and in its decay all matter and entropy in the Universe is produced. Therefore, it is natural to assume that sterile neutrinos can be produced as well [664]. While in a model considered in ref. [664] the inflaton field did have a VEV and by that provided the single source for the scale of electroweak symmetry breaking and for Majorana masses of sterile neutrinos, in general this is not necessary. Sterile neutrino production in the inflaton decay is always possible, subject to the value of the coupling constant $y$.

The remainder of this subsection is devoted to sterile neutrino production from a thermalized population of $\phi$-particles. Even thermalized quanta of the inflaton field are possible, as in the model of ref. [664]. The population of sterile neutrinos which originate from decays of $\phi$-particles is governed by the standard Boltzmann kinetic description. The result is definite as long as $\phi$-particles are in thermal equilibrium, and in a rough approximation the distribution function of the scalars translates into that of sterile neutrinos.

The distribution function of sterile neutrinos $f(p, t)$, where $p$ is the momentum and $t$ is the time, can be found from the solution of the kinetic equation [663, 664, 670]:

$$ \frac{\partial f}{\partial t} - H p \frac{\partial f}{\partial p} = \frac{2 m_{\phi} \Gamma}{p^2} \int_{p+m_{\phi}^2/(4p)}^{\infty} f_{\phi}(E) dE, \qquad (6.29) $$

where it was assumed that the inverse decays $N_1N_1 \to \phi$ can be neglected (this is true for small Yukawa couplings $y < 10^{-7}$). In eq. (6.29), $H$ is the Hubble constant, $f_{\phi}(E)$ is the distribution of $\phi$-particles, and $\Gamma = m_{\phi}f^2/(16\pi)$ is their partial width for the $\phi \to N_1N_1$ decay channel. For the case when the effective number of degrees of freedom is time-independent and $f_{\phi}(E)$ is thermal, a semi-analytic solution to eq. (6.29) can be found easily. At $t \to \infty$ it takes the form

$$ f(x) = \frac{16\Gamma M_0}{3m_{\phi}^2} x^2 \int_1^{\infty} \frac{(y-1)^{3/2} dy}{e^{xy}-1}, \qquad (6.30) $$

where $x = p/T$ and $M_0 \approx M_{\text{Pl}}/(1.66\sqrt{g_*})$, leading to a number density of

$$ n = \int \frac{d^3 p}{(2\pi)^3} f(p) = \frac{3\Gamma M_0 \zeta(5)}{2\pi m_{\phi}^2} T^3 \qquad (6.31) $$

and to an average momentum of created sterile neutrinos immediately after their production equal to $\langle p \rangle = \pi^6/(378\zeta(5))T = 2.45T$, which is about 20% smaller than that for an equilibrium thermal distribution, $p_T = 3.15T$.

The resulting abundance of sterile neutrinos is given by [664]:

$$ \Omega_{N_1} \sim \frac{y^2}{S_*} \frac{M_{\text{Pl}}}{m_{\phi}} \frac{M_1}{\text{keV}}. \qquad (6.32) $$

It is important to note that a major contribution to $\Omega_{N_1}$ is coming from the epoch when $T \sim m_{\phi}$, since later on the number density of $\phi$-particles is exponentially suppressed if they remain in thermal equilibrium. The dilution factor

$$ S_* \equiv S \frac{g_*(T_{\text{prod}})}{\tilde{g}_{*0}} = S \frac{g_*(T_{\text{prod}})}{3.9} \qquad (6.33) $$

in eq. (6.32) is the ratio of the number of effective degrees of freedom at production time, $g_*(T_{\text{prod}})$, till now, $\tilde{g}_{*0} = 3.9$ (this number takes into account the change of the temperature
of neutrinos at freezeout). $S_*$ also contains the entropy production factor $S$ if this process happens during the QCD transition, see section 6.4 and ref. [662]. A more refined analysis is required if the number of d.o.f. changes significantly at $T \sim T_{\text{prod}} \sim m_{\phi}/3$, see [664].

Today, the averaged momentum of sterile neutrinos produced in decays of particles at $T_{\text{prod}}$ is given by [664]:

$$ \frac{\langle p \rangle}{T_{\gamma}} = \frac{\pi^6 S_*^{-1/3}}{378 \zeta(5)} \approx 2.45 S_*^{-1/3}. \qquad (6.34) $$

Comparing this to non-resonant production with $\langle p_{\text{NRP}} \rangle/T_{\gamma} = 3.15(4/11)^{1/3}$,<fnref n="33" /> we can see, that the mass bounds from structure formation<fnref n="34" /> are weaker than for the NRP case [677]

$$ m_{\text{Decay}} = \frac{\langle p \rangle}{\langle p_{\text{NRP}} \rangle} m_{\text{NRP}} = \frac{2.45}{3.15} \left( \frac{1}{S} \frac{10.75}{g_{*}(T_{\text{prod}})} \right)^{1/3} m_{\text{NRP}}. \qquad (6.35) $$

This is consistent with resonant production being generally threatened by Lyman-α bounds [740]. However, it should be noted that, taking into account the dilution during the QCD transition, $\langle p_{\text{NRP}} \rangle/T_{\gamma}$ can in fact be smaller than $3.15(4/11)^{1/3}$ by about 20% [344].

**Light inflaton model [664].** Let us now turn out attention to one specific example model that makes use of the results described above. Let us assume, that the field $\phi$ in (6.28) is the field that leads to the inflationary expansion of the Universe.

First, let us briefly discuss the processes right after inflation and before reheating. They correspond to a particle creation by a coherently oscillating classical scalar field. Using the results of refs. [741, 742], we find that the ratio of the number density of light fermions produced to the entropy scales as $n_{\nu}/n_{\text{rad}} \sim y^{3/2}$. Fermions will be cold at the production, $\langle p \rangle \sim y^{1/2}T$. However, very light fermions in a keV mass range will contribute to the DM abundance appreciably for rather large values of the Yukawa coupling, $y \sim 10^{-2}$. We will not consider such types of models further in what follows.

The mass of the field and its VEV in this model are related as $m_{\phi} = \sqrt{2\lambda_{\phi}} \langle \phi \rangle$, where $\lambda_{\phi}$ is the quartic scalar self-coupling. For a sufficiently light inflaton, with a mass even smaller than the Higgs mass, the inflaton does not decay completely but instead thermalizes and is then in thermal equilibrium down to rather low temperatures, $T \ll m_{\phi}$. In that case, eq. (6.32) yields $\Omega_{N_1} \sim 10^{25} \lambda_{\phi} (M_1/m_{\phi})^3$. Sterile neutrinos can constitute all the DM e.g. for $m_{\phi} \approx 0.3 \text{ GeV}$ and $M_1 \approx 20 \text{ keV}$, if $\lambda_{\phi} \sim 10^{-13}$. Light inflatons could be specifically searched for at LHC, see refs. [731, 743]. Small values for the inflaton self-coupling, quoted above, are required by the CMB constraints in models with minimal coupling to gravity, and such constraints date back to the epoch when $\phi^4$ inflationary models were not yet ruled out by observations. Otherwise, with a non-minimal coupling to gravity, $\lambda_{\phi}$ could be larger for the inflaton field [743] and may even approach unity [744].

For an inflaton mass of $m_{\phi} < \mathcal{O}(500) \text{ MeV}$, taking $g_* = \text{const.}$ is not a good approximation, since exactly in this region $g_* \sim 60$ at $T \sim 1 \text{ GeV}$ to $g_* \sim 10$ at

$$ m_{\text{NRP}} = m_{\text{TR}} \left( \frac{1}{\Omega_{\text{DM}} h^2} \frac{m_{\text{TR}}}{94 \text{ eV}} \right)^{1/3} . $$

<footnote n="33">³³It has been pointed out in ref. [344] that the approximation $\langle p_{\text{NRP}} \rangle/T_{\gamma} = 3.15(4/11)^{1/3}$ does not hold exactly and $\langle p_{\text{NRP}} \rangle$ is in fact slightly smaller.</footnote>

<footnote n="34">³⁴There are two reference numbers that are usually present in the analysis of the structure formation bounds described in sections 5.1, 5.2, bound for the NRP produced sterile neutrino and bound for the thermal relic, which are related as</footnote>
$T \sim 1$ MeV, because of the disappearance of all the quark and gluon degrees of freedom. In this case a numerical solution of eq. (6.29) is required, with the input of the hadronic equation of state, which is however not known exactly. In [664], numerical calculations have been presented with the use of $g_*(T)$ constructed in [320] on the basis of the hadron gas model at low temperatures, and from the available information on lattice simulations and perturbative computations. It was found that the abundance, eq. (6.31), would change by a factor of $\tilde{S}_*(m_\phi) \approx 0.9$ at $m_\phi = 70$ MeV to 0.4 at 500 MeV. The average momentum stays almost unchanged in this interval of $m_\phi$. For higher inflaton masses, a good approximation to $\tilde{S}_*(m_\phi)$ is $\tilde{S}_*(m_\phi) \simeq [10.75/g_*(m_\phi/3)]^{3/2}$, leading to an average momentum of $\langle p \rangle \simeq 2.45T_\nu(10.75/g_*(m_\phi/3))^{1/3}$. The function $g_*(T)$ can be read off from figure 2 of ref. [320], and $\langle p \rangle$ could be as low as $1.5T$ taking into account entropy dilution by decays of heavier sterile neutrinos.

### 6.3.2 Production from generic scalar singlet decays

Another possibility is to produce the sterile neutrinos from the decay of a general singlet scalar $S$. This idea was first proposed in refs. [663, 670], where the scalar was assumed to enter thermal equilibrium and either decay immediately or freezing-out before decaying to sterile neutrinos. In this case the distribution of the sterile neutrino can be found approximately as a sum of two components [670]:

$$f(x, r) = f_1(x, r_f) + f_2(x, r), \qquad (6.36)$$

where $x = p/T$ is the neutrino momentum and $r = m_S/T$ is a time-parameter. At earlier times (lower momenta) the sterile neutrino is produced from in-equilibrium decays of the parent particle. At larger momenta (later times), sterile neutrinos are produced from the decay of the parent particle which is in thermal equilibrium. In the most generic situation the process is governed by three dimensionless parameters, $r_f = m_\phi/T_f$ giving the moment of the freeze-out of the parent particle, $\Lambda = \Gamma_{tot}M_0/m_\phi^2$ characterizing the total decay width of the parent particle, and by branching of the decay into the sterile neutrino $B = \Gamma(\phi \to N_1N_1)/\Gamma_{tot}$. Then the in-equilibrium decay contributes

$$f_1(x, r_f) = 2B\Lambda \left[ \frac{r_f^3}{3x^2} \ln \frac{1}{1 - e^{-x-(r_f^2/4x)}} + \frac{8x^2}{3} \int_1^{1+(r_f^2/4x^2)} \frac{(z-1)^{3/2} dz}{e^{xz}-1} \right]. \qquad (6.37)$$

The out-of-equilibrium component is$^{35}$

$$f_2(x, r) = \frac{B}{x^2} \left[ \int_{|(r_f^2/4x) - x|}^{\infty} x_{\phi} f_{\phi}(x_{\phi}, r_f) dx_{\phi} - \int_{r_f}^{r} \frac{r'}{2x} \left( \frac{r'^2}{4x} - x \right) f_{\phi} \left( \left| \frac{r'^2}{4x} - x \right|, r' \right) dr' \right]. \qquad (6.38)$$

with the concentration of the $\phi$ particle decaying as

$$f_{\phi}(x_{\phi}, r) = \frac{1}{e^{\sqrt{x_{\phi}^2 + r_f^2}} - 1} \left( \frac{r + \sqrt{x_{\phi}^2 + r_f^2}}{r_f + \sqrt{x_{\phi}^2 + r_f^2}} \right)^{\Lambda x_{\phi}^2} \times e^{-\Lambda(r\sqrt{x_{\phi}^2 + r_f^2} - r_f\sqrt{x_{\phi}^2 + r_f^2})}. \qquad (6.39)$$

Note, that in analysis of the present day momentum distributions, the overall momenta are divided by the factor $S_*^{1/3}$ from (6.33). An important observation [730] is that these

<footnote>$^{35}$One extra term present in [670] vanishes at late times $r \to \infty$.</footnote>
<table><caption>**Table 10.** Analytically accessible limiting cases of the production of sterile neutrinos from singlet scalar decay. The quantity $M_0$ is the rescaled Planck mass, given by $M_0 = (\frac{45M_{\text{Pl}}^2}{4\pi^3 g_* (T_{\text{prod}})})^{1/2} = 7.35 g_*^{-1/2} (T_{\text{prod}}) \times 10^{18}$ GeV, where $T_{\text{prod}}$ is the temperature at which the production of the scalar and its subsequent decay takes place. Furthermore, $T_{\text{FO}}$ is the freeze-out temperature that depends on $m_S$ and $\lambda$, while $K_2$ denotes the second Bessel function of second kind.</caption><thead><tr><th>Regime</th><th>Necessary condition for λ</th><th>Necessary approximation</th><th>Relic yield Y<sub>νS</sub> = n/s (t → ∞)</th></tr></thead><tbody><tr><td>Scalar freeze-in</td><td>small enough not to equilibrate the scalar</td><td rowspan="3">small enough to ensure<br>g<sub>*</sub> = const.<br>is a good approximation</td><td><sup>135</sup>&frasl;<sub>1024π<sup>5</sup></sub> <sup>M<sub>0</sub>(T<sub>prod</sub>)</sup>&frasl;<sub>m<sub>S</sub></sub> <sup>λ<sup>2</sup></sup>&frasl;<sub>g<sub>*</sub>(T<sub>prod</sub>)</sub></td></tr><tr><td>Scalar freeze-out + decay in equilibrium</td><td>large enough to equilibrate the scalar</td><td><sup>135</sup>&frasl;<sub>64π<sup>4</sup></sub> <sup>M<sub>0</sub>(T<sub>prod</sub>)</sup>&frasl;<sub>m<sub>S</sub></sub> <sup>y<sup>2</sup></sup>&frasl;<sub>g<sub>*</sub>(T<sub>prod</sub>)</sub></td></tr><tr><td>Scalar freeze-out + decay out of equilibrium</td><td>large enough to equilibrate the scalar</td><td><sup>45m<sub>S</sub><sup>2</sup>K<sub>2</sub>(m<sub>S</sub>/T<sub>FO</sub>)</sup>&frasl;<sub>4π<sup>4</sup>g<sub>*</sub>(T<sub>prod</sub>)T<sub>FO</sub><sup>2</sup></sub></td></tr></tbody></table>

two components can have contributions to the sterile neutrino population with significantly different momenta (see figure 27). In general, if the amount of thermal degrees of freedom $g_*$ changed during the production process, an explicit solution of the Boltzman equations is required [730]. Another subtlety appears at early times, when the temperature of the plasma exceeds mass of the decaying particle [685]. This case is mostly relevant if the decay width of the scalar is very large, so that the scalar should be considered as an intermediate state.

An alternative proposal was put forward in ref. [676], which pointed out that the scalar could also freeze-in and produce a suitable abundance. In ref. [683] this mechanism has been extended to the case of frozen-in scalars lighter than the Higgs mass, motivated by the recent hints for an X-ray line at 3.55 keV. The most general numerical study up to now, covering all the cases mentioned, was presented in ref. [730], which in particular confirmed earlier estimates indicating the sterile neutrinos produced from scalar decays can have a comparatively cold spectrum [669, 688, 732], see also ref. [745] for a recent review. While we here focus on the case of a total singlet $S$, note that electrically charged scalars could resemble the case of an equilibrated singlet scalar [746], since for the DM production it does not matter by which interaction the scalar is kept in equilibrium.

**The model.** From the aforementioned variety of models for scalar decay we pick one which is quite illustrative and allows for a semi-analytic discussion of some relevant aspects. This model is a minimal extension of the SM, featuring one real scalar (which we denote $S$), singlet under the SM and one right handed neutrino (denoted $N_1$). Again, the sterile neutrino is coupled to the scalar via a Yukawa-like interaction, while the scalar is coupled to the Higgs sector via the most general potential $2\lambda(H^\dagger H)S^2$ allowed after imposing a $\mathbb{Z}_4$-symmetry, thus forbidding the decay to the SM particles and setting the parameter $B = 1$ in eqs. (6.37), (6.38).$^{36}$

We stick to the case where the mass of the scalar is above the EW scale. As in the preceding part of section 6.3, the production of sterile neutrinos is exclusively governed by

<footnote>$^{36}$For detailed comments on this assumption, see [676, 730] and references therein.</footnote>
the reactions

$$SS \leftrightarrow hh \quad \text{and} \qquad (6.40)$$

$$S \rightarrow N_1 N_1, \qquad (6.41)$$

to lowest order in perturbation theory.<fnref n="37" /> Again we neglect the inverse reaction of two sterile neutrinos producing one scalar. This is justified by kinematical arguments and by the low number density of the sterile neutrinos, which implies that we can use (6.29) to convert the distribution function of the scalar, once known, into the distribution of the sterile neutrino.

The exact interplay of the reactions (6.40) and (6.41) is governed by two couplings, $\lambda_S$ and $y$. While $\lambda_S$ determines whether or not the scalar ever enters equilibrium and — if so — when it decouples from the plasma, the Yukawa coupling $y$ controls the decay width of the scalar and by that virtue has major impact on the production time of the sterile neutrinos.

For different regimes in the parameter space spanned by these couplings, we will calculate the distribution function $f_S$ of the scalar to subsequently obtain the distribution function of the sterile neutrino. Note, that the distribution function contains all relevant information: apart from the number density, the spectral form will also allow to derive consequences for structure formation, BBN, CMB, and so on.

**Cosmological dynamics.** Let us now turn to the description of the dynamics of the scalar in the different regimes. The abstract Boltzmann equation is:

$$\frac{\partial f_S(p,t)}{\partial t} = Q(p,t, m_H, m_S, \lambda) - \mathcal{P}(p,t, m_S, y) f_S(p,t) \\ - \mathcal{R}(p,t, m_H, m_S, \lambda) f_S(p,t) \int d^3p' \mathcal{S}(p',t, m_H, m_S) f_S(p',t), \qquad (6.42)$$

in the limiting case where $g_*$ can be approximated as being constant during production.<fnref n="38" /> The functions $P, Q, R, S$ encode all the kinematics of the processes. The term $Q$ describes the production of scalars from Higgs bosons and hence does not depend on $f_S$, while the integral term (containing $R$ and $S$) describes the backreaction and therefore scales as $f_S^2$. The depletion of scalars by their decay into two sterile neutrinos is given by the term containing $P$ and scaling linearly with $f_S$. For three limiting cases, the combination of (6.29) and (6.42) can be solved analytically for the relic abundance of sterile neutrinos. One of these cases is basically equivalent to the one discussed in eqs. (6.30) to (6.32). For the sake of completeness, we will summarize all three of them in table 10.

**Numerical results.** We have solved eq. (6.42) numerically for different regimes of the couplings $\lambda$ and $y$. Figure 26 shows lines in the space spanned by $\lambda$ and $y$ reproducing the relic abundance as observed by Planck [1] for different masses $m_S$ of the sterile neutrino. It is visible that there is a gap in $\lambda$, around which it is not possible to obtain the required value of the relic density without violating the Tremaine-Gunn (TG) bound [75]. Accordingly, the possible parameter space is clearly divided into the part to the left ($\lambda \lesssim 10^{-6.5}$), where the scalar freezes in and the part with $\lambda \gtrsim 10^{-5.5}$, where the scalar freezes out. In the latter regime, there exists a rather global upper limit on the Yukawa coupling $y$. If the scalar equilibrates and the decay width of the scalar, $\propto y^2$, is too large, the particle number density of sterile neutrinos will be too high to be consistent with TG bound. As one of the major results of the

<footnote n="37">³⁷For the case of $m_S < v_{EW}$ studies on the level of particle number densities can be found in [683].</footnote>

<footnote n="38">³⁸This can directly be mapped to [730, eq. (10)], where the authors performed a change of variables on (6.42).</footnote>
<figure><img src="image_27.png" /><figcaption>**Figure 26.** Isoabundance lines in the plane spanned by the Higgs portal $\lambda$ and the Yukawa coupling $y$. The abundance as observed by Planck can be reproduced with masses in the keV-range.</figcaption></figure>

numerical study, we present an intermediary freeze-out case in figure 27. In this case, both the decay *in and out* of equilibrium contribute in roughly equal parts to the final abundance. However, since the production times of these two components are significantly different, the spectrum features *two* peaks. This shows that the simple free-streaming estimates might *not* be sufficient and more elaborate analyses have to be done when assessing the observations from structure formation. The most advanced discussion of this mechanism, featuring a detailed account of structure formation bounds, has recently been provided in [747].

## 6.4 Dilution of thermally produced DM (author: F. Bezrukov)

In this section we discuss how to obtain the proper amount of the light keV scale DM if the DM candidate reached thermal equilibrium at some high temperature in the early Universe. This does not happen in νMSM, where the only source of interactions for $N_1$ is its minuscule Yukawa coupling constant, but is readily achieved in most high energy extensions of SM. Specifically, this happens if $N_1$ participates in gauge interactions of an extension of the SM gauge group. A recent detailed discussion of this possibility without reference to specific particle physics models is given in [748]. Gauge interactions If the “sterile” neutrinos are charged under some gauge group, their equilibration in the early universe is almost unavoidable unless their freeze-out temperature exceeds the maximal temperature of the primordial plasma. The freeze-out temperature of the sterile neutrino $T_{\text{f,c}} \sim g_*^{1/6} (M_{W_R}/M_W)^{4/3} \times 1$ MeV goes up with the mass $M_{W_R}$ of the new gauge interactions, so for significantly heavy $M_{W_R}$ it may exceed the preheating temperature after inflation.<fnref n="39" />

According to the standard lore if a light species reaches thermal equilibrium and freezes out while still being relativistic, then its abundance (or number density in comoving volume) is fixed in a unique way at the moment of freeze-out, and its contribution to the energy

<footnote n="39">³⁹This temperature is unknown; it may be constrained by CMB observations in models where reheating is mainly driven by perturbative processes [749, 750].</footnote>
<figure><img src="image_28.png" /><figcaption>**Figure 27.** Evolution of the distribution function over the timelike parameter $r = m_S/T$ (left panel) and evolution of the particle number densities (right panel). The plot of the distribution function in the left panel shows $x^2f_N(x)$ to give a intuitive picture of the number of particles with the corresponding momentum. In the right panel, $\tilde{n}_S$ ($y=0$) shows how the scalar would evolve if it was stable while $\tilde{n}_S^{\text{therm}}$ depicts how a thermalized particle would evolve. Comparing both panels, it is clearly visible how the decay in equilibrium shapes the peak around lower momenta early on while the second peaks comes from converting the frozen-out relic density of scalars at late times. The two scales result in a rather high average momentum of $\langle x \rangle = \langle p/T \rangle \approx 39.4$ which, however, is not a very good measure to describe the sum of two peaks lying at higher / lower momenta.</figcaption></figure>

balance at present, when it becomes non-relativistic, is simply proportional to the mass of the particle. This simplified analysis gives a bound of $\lesssim 10$ eV for the DM mass. However, a more detailed analysis [671] should take into account the subtleties of the expansion history of the Universe after the freeze-out moment, giving the contribution of the sterile neutrino $N_1$ to the present energy density of the Universe

$$ \Omega_{N_1} \approx 0.265 \frac{1}{S} \left( \frac{10.75}{g_*(T_{f,c})} \right) \left( \frac{M_1}{1 \text{ keV}} \right) \times 90. \qquad (6.43) $$

Here, $S$ is the amount of entropy that was released in non-equilibrium processes after the freeze-out$^{40}$ (i.e. $s_0/a_0^3 = Ss_f/a_f^3$), and $g_*(T_{f,c})$ is the effective number of degrees of freedom in the plasma at the moment of the DM neutrino freeze-out. We will not discuss adding 1000 d.o.f. at freeze-out, but will consider the possibility of entropy release after the DM freeze-out.

The entropy can be produced in out-of-equilibrium decays of some heavy and relatively long-lived particle [19, 751]. In most models there are several sterile neutrinos in addition to the DM $N_1$, and the heavier $N_{2,3}$ can conveniently used for entropy generation (though, one can imagine models with other source of entropy). Assuming that they decay shortly before nucleosynthesis at MeV temperature, their decay leads to the DM abundance [673]

$$ \Omega_{N_1} \approx 0.265 \left( \frac{M_1}{1 \text{ keV}} \right) \left( \frac{1.6 \text{ GeV}}{M_{N_{2,3}}} \right) \left( \frac{1 \text{ sec}}{\tau_{N_{2,3}}} \right)^{1/2} \left( \frac{g_*(T_{f,2/3})}{g_*(T_{f,1})} \right). \qquad (6.44) $$

<footnote>$^{40}$In fact, a factor of $S = 1-2$ is present in the plain $\nu$MSM as well, and it helps to widen the allowed area of parameters [662].</footnote>
Thus, for given masses of DM sterile neutrino $N_1$ and entropy generating one $N_{2,3}$ it is possible to define the lifetime $\tau_{N_{2,3}}$, leading to the proper DM abundance today. There are two main constraints on the allowed values of the parameters:

* Successful nucleosynthesis requires at least $\tau_{N_{2,3}} \lesssim 1$ sec, with (6.44) leading to the lower bound on the heavy sterile neutrino mass.

* The diluter should be relativistic at its freeze-out, $T_{f,h} > M_{2,3}$ (otherwise the entropy release is diminished). This can be converted into the lower bound on the scale of the gauge interactions, that lead to the thermalization of the diluter.

Simple estimates for these bounds give the heavy sterile neutrinos heavier than $M_{2,3} \gtrsim (M_1/1 \text{ keV}) \times 1.6 \text{ GeV}$ and the mass of the additional gauge boson, leading to thermalization of $N_{2,3}$, larger than $M_{W_R} \gtrsim (M_1/1 \text{ GeV})^{3/4} \times 10 \text{ TeV}$. In [674] it was shown that the decay of the heavy left neutrinos leading to the entropy dilution can also lead to leptogenesis, which proceeds via usual thermal high temperature scenario.

More careful analysis for a minimal left-right symmetric theory was made in [673]. Taking into account two diluting particles $N_2$ and $N_3$, and tuning the flavour structure to separate the freeze-out temperatures $T_{f,h}$ and $T_{f,c}$ of the dilutors and DM sterile neutrinos, the constraints can be relaxed in a very specific corner of the parameter space.

The momentum distribution of the produced sterile neutrinos is also cooled by the entropy dilution, leading to the weaker bound from the Lyman-$\alpha$. The spectrum of the DM corresponds to thermal spectrum with the temperature, reduced by a factor of $S^{-1/3}$. This corresponds to the thermal relic case in the terminology of the majority of the Lyman-$\alpha$ paper, and allows to use the bounds from the recent papers without modifications. The latest bound from [331] gives

$$M_1 > 3.3 \text{ keV}(2\sigma) \text{ and } 2.5 \text{ keV}(3\sigma),$$

while the older very conservative bound from [59] gives for thermal relic, $M_1 > 1.5 \text{ keV}$ at 99.7% C.L.

Last but not least, in the models with entropy dilution one should check that the lifetime of the DM sterile neutrino is still compatible with the X-ray observations, as far as additional interactions often lead to increase in the $N_1 \to \nu\gamma$ decay rate. For example in the LR symmetric models a bound on the $W_R - W_L$ mixing appears in addition to the bound on the sterile-active neutrino mixing.

To conclude the section we should note variations of the scenario. In [671, 673] entropy dilution was achieved in the context of Left-Right symmetric extensions of the SM. A model with radiative generation of active neutrino masses was analyzed in [752, 753], model with composite neutrinos in [754]. Note, however, that quite generally a large amount of entropy dilution is not easy to get in agreement with big bang nucleosynthesis [755].

Of course, similar setup with the DM abundance given by (6.43) can be realized with DM candidates that are not sterile neutrinos, but some other light stable particle, see [756] as a recent example.
# 7 keV neutrino theory and model building (particle physics)

Section Editors:
Alexander Merle, Viviana Niro

In this section, we will address some of the most generic mechanisms used in the literature to motivate light sterile neutrinos in general and sterile neutrinos with keV masses in particular. While there are some similarities to eV-sterile neutrinos, we will in particular consider the small active-sterile mixing, which is characteristic for sterile neutrinos with keV-scale masses acting as Dark Matter (DM).

## 7.1 General principles of keV neutrino model building (authors: A. Merle, V. Niro)

In the SM of particle physics, sterile neutrinos are viewed as mass eigenstates which have predominantly right-handed chirality, i.e. they do not couple to SM gauge bosons.<fnref n="41" /> These fields are typically called in the literature “right-handed” neutrinos if they are heavy (GeV or higher) and “sterile neutrinos” if they are light, with mass around eV-keV or MeV. However, this distinction is to some extend artificial and we will use the two terms in this section in an interchangable way.

These types of fermions can obtain mass from a Majorana mass term, cf. section 1, which does a priori not have any relation to the electroweak symmetry breaking scale. Thus, it will be natural to expect the mass of sterile neutrinos to be very heavy. For this reason, it is necessary to develop mechanisms which can yield light sterile neutrinos and which can also protect a potentially small sterile neutrino mass from higher order corrections.

In the literature one can find a number of example mechanisms that achieve this goal. To illustrate the generic principles, we can divide the known mechanisms into two categories: mechanisms which suppress a naturally large sterile neutrino mass and mechanisms which give a small but non-zero mass correction to a sterile neutrino mass that is exactly zero at leading order. In ref. [66], these two categories have been named *top-down* and *bottom-up* scheme, respectively, see figure 28. This division may not be entirely exhaustive, as we will discuss in section 7.4.1 and section 7.4.2, but it nevertheless illustrates a generic principle: in the framework of a quantum field theory it is impossible to obtain absolute mass scales, but we can nevertheless explain strong *mass patterns*, i.e., hierarchies or (quasi-)degeneracies between masses.

The mechanisms presented in section 7.2 all rely on the principle of suppressing one sterile neutrino mass eigenvalue to a value much smaller than its “natural” scale (and thus smaller than the masses of other sterile neutrinos present in the model). Such suppressions typically come from new physics affecting the different generations of sterile neutrinos in different ways, which may break mass degeneracies or strongly enhance a mass splitting that is already present. On the other hand, mechanisms based on discrete or continuous flavor symmetries, see section 7.3, typically predict one strictly massless sterile neutrino in the limit of exact symmetry, along with one or several massive ones. Once (part of) the symmetry is broken, e.g. by higher order corrections, this zero mass eigenvalue receives corrections proportional to the symmetry breaking parameter which must be smaller than the mass scale of the heavier sterile neutrinos, as otherwise we would not have the symmetry in the first place. Other models where the suppression is obtained with a more complicated mechanisms

<footnote n="41">⁴¹Note that this might be different in settings with an extended gauge groups, such as left-right models.</footnote>
<figure><img src="image_29.png" /><figcaption>**Figure 28.** The two generic mass shifting schemes used in models for light sterile neutrinos. Most of the mechanisms proposed ultimately apply some variant of one of those schemes. In both cases, for $M_1$, the “natural” mass scales are indicated as well as how the shift to the keV scale. (Figure similar to figure 4 in ref. [66].)</figcaption></figure>

are presented in section 7.4. In general, all these models generate a hierarchy between the light (nearly massless) sterile neutrino and the heavier sterile neutrinos involved.

One final point that is important for all mechanisms is the stability of the mass generation mechanism. Whichever mechanism we consider, it will generically receive corrections, e.g. from quantum loop contributions or from higher order terms, neglected at the order under consideration. While such corrections are often small in practice, there may be cases in which they can spoil the conditions required for the mechanism to work. This should be cross-checked in order to verify the validity of a certain mechanism. In most of the examples presented, this check has been performed in the respective references, see for example the models presented in sections 7.3.2 and 7.3.3, but we will not always explicitly comment on it in this manuscript. However, we invite the reader to consult the original references which provide many more details than we could possible illustrate in the framework of this section.

## 7.2 Models based on suppression mechanisms

In this section, we will introduce the first class of mechanisms which can motivate a light sterile neutrino mass scale, namely those based on the principle of suppressing a naturally large sterile neutrino mass to a value at the keV scale.

### 7.2.1 The split seesaw mechanism and its extensions (author: R. Takahashi)

We start our discussion with two models based on a suppression mechanism for the neutrino Yukawa couplings, both of which rely on the existence of extra spatial dimensions. The models are the split [672] and the separate [757] seesaw mechanisms, which are considered in a flat five-dimensional space-time. The five-dimensional space is assumed to be compactified on $S^1/Z_2$ and the five-dimensional coordinate is denoted as $y \equiv x^5$. The fundamental region is given by $y \in [0, \ell]$. There are two fixed points at $y = 0, \ell$. We consider branes (=four-dimensional subspaces) at the fixed points, i.e., one brane at $y = 0$ is a standard model (SM) brane and the other at $y = \ell$ is a hidden brane. The SM particles reside at the SM brane.