Intimate Relationship Between Sterile Neutrino Dark Matter and $\Delta N_{\rm eff}$

The self-interacting neutrino hypothesis is well motivated for addressing the tension between the origin of sterile neutrino dark matter and indirect detection constraints. It can also result in a number of testable signals from the laboratories to the cosmos. We explore a model of neutrino self-interaction mediated by a Majoron-like scalar with sub-MeV mass, and show that explaining the relic density of sterile neutrino dark matter implies a lower bound on the amount of extra radiation in early universe, in particular $\Delta N_{\rm eff}>0.12$ at the CMB epoch. This lower bound will be further strengthened with an improved $X$-ray search at the Athena observatory. Such an intimate relationship will be unambiguously tested by the upcoming CMB-S4 project.

Introduction -The origin of dark matter (DM) poses a fascinating puzzle. A sterile neutrino with a keV-scale mass and small mixing with active neutrinos is a simple, longstanding DM candidate. Along with other traditional candidates (e.g., WIMPs/axions), sterile neutrinos are highly testable. In the minimal setup, the activesterile neutrino mixing that can account for the observed DM relic density [1] is already in tension with searches for DM decaying into monochromatic X-rays [2]. An elegant solution for alleviating this tension is to introduce a new self interaction among active neutrinos. This keeps neutrinos in thermal equilibrium with themselves longer in the early Universe and facilitates efficient sterile neutrino DM (SνDM) production [3,4]. This mechanism is attractive because it predicts a number of novel signatures for laboratory tests [5][6][7]. Interestingly, the self-interacting neutrino scenario emerges from a broad class of wellmotivated beyond-Standard-Model frameworks [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24].
We will soon enter a new precision era of cosmology. The Cosmic Microwave Background -Stage 4 (CMB-S4) project [25] is expected to measure cosmological parameters to an unprecedented level. In particular, it could restrict the amount of extra radiation (in units of extra neutrinos, ∆N eff ) during the CMB epoch to a fewpercent uncertainty, while it is still allowed (perhaps even favored) to be order one by the current data [26,27]. This could shed light on existing tensions resulting from various ways of measuring the Hubble parameter [28][29][30][31], as well as fundamental theories that can accommodate deviations from the standard cosmology. Moreover, ∆N eff measurement during the big-bang nucleosynthesis (BBN) epoch puts further restrictions on room for new physics.
In this letter, we point out an intriguing relationship between the two important ingredients of our Universethe origin of SνDM and ∆N eff for both CMB and BBN. The connection is established by a light mediator with a feeble coupling to active neutrinos such that it enters the thermal bath after BBN but decays away before CMB. SνDM then is dominantly produced from active neutrino self scattering through the exchange of on-shell mediator. The post-BBN thermalization of the mediator inevitably modifies the standard cosmology predictions by making new contributions to ∆N cmb eff and ∆N bbn eff . As a result, the viable parameter space for SνDM relic density is highly correlated with the predictions of ∆N eff . We further point out a new implication of the upcoming X-ray search by the Athena experiment for the next-generation precision measurements in cosmology-this scenario connects different frontiers of cosmological and astrophysical probes.
Models -SνDM is a linear combination of a gauge singlet fermion and the active neutrino states from the Standard Model [1], ν 4 = ν s cos θ + ν a sin θ, where ν 4 denotes a mass eigenstate and θ is the active-sterile mixing angle in vacuum. The sterile neutrino is a decaying DM candidate whose longevity is attributed to its small mass/mixing. We introduce a new force among active neutrinos, which can arise from dimension-six gauge invariant operators at energy scales well below the weak scale. The scalar field φ could be either real or complex. If real, φ could be identified as the Majoron [8,9], where λ = iM ν /f and f is the spontaneous lepton number breaking scale. If complex, φ also contains the radial mode and restores B − L symmetry at high scales [5]. The impact of active neutrino self-interactions on the relic density of SνDM was first explored in [3]. The new force enables more frequent active neutrino scattering than normal weak interactions, thereby, SνDM can be produced with a smaller mixing angle than required by Dodelson-Widrow [1]. Ref. [3] focuses on m φ between MeV-GeV, where the new interaction is sufficiently strong to fully thermalize φ before BBN. The present work broadens the focus and opens up the possibility of a sub-MeV mediator φ. The constraint from

This Work
FIG. 1. The Big Picture: full parameter space that can account for the relic density of SνDM (thick black curve). We choose m4 = 4 keV and sin 2 2θ = 5 × 10 −9 . The colored regions show the experimental constraints on neutrino selfinteraction mediated by a real scalar φ, including leptonic meson decays [5,6,17,32] (purple), invisible Z width [5,7] (orange), double beta decay [33][34][35] (green), and ∆N eff for successful BBN [36] (blue). We will explore the SνDM production mechanism in the the sub-MeV m φ region marked by "This Work" which could be fully tested by the upcoming CMB-S4 and Athena experiments.
BBN is lifted if φ does not thermalize before/during the BBN epoch and has no primordial population from other high-scale physics [37][38][39][40][41][42][43][44][45]. This is achieved with a feeble coupling to neutrinos and a phase space suppression for νν → φ to occur at temperatures much higher than m φ . The condition for not thermalizing φ before BBN sets an upper bound on λ that scales as 1/m φ , whereas correct DM relic density requires λ to scale as √ m φ in the small coupling regime. As a result, viable parameter space opens up for m φ between keV and MeV scales. Remarkably, our result, in tandem with that from Ref. [3], nicely forms the "Big Picture" of Fig. 1. The thick black curve depicts the entire range of the theory space where the relic density of SνDM is explained, aided by neutrino self-interaction. Viable mediator mass range divides into two windows, above and below MeV scale. They correspond to different histories of the neutrino and SνDM sector in the early universe. In the heavier mediator region, the strongest probes are from terrestrial experiments. In contrast, the sub-MeV mediator regime will be best explored by cosmological observations of ∆N eff , to be detailed below.
Cosmology -SνDM production is governed by where and f 4 /f a are phase space distribution functions of sterile/active neutrinos, respectively. We introduce z ≡ m φ /T sc ν to label time, where T sc ν is defined as the neutrino temperature in the standard cosmology without φ; z expands linearly with the Universe's radius. The actual neutrino temperature T ν is affected by the interaction in Eq. (1). After weak-interaction decoupling, T ν deviates from T sc ν as the φ population builds. Γ is the thermally-averaged scattering rate for an active neutrino, and H is the Hubble parameter. Since we are interested in small m φ and λ, the effective active-sterile mixing angle is approximately equal to the vacuum one throughout the dominant SνDM production window.
The distribution function f a , after active neutrinos thermalize with φ, is characterized by T ν and a chemical potential µ ν , To derive their z dependence, we introduce three time scales, • z A : the onset of BBN -neutrinos have just decoupled from weak interaction. We assume φ has not yet equilibrated with neutrinos; • z B : the ν-φ system has just equilibrated. In general, φ remains ultra-relativistic; • z C : all φ have decayed away (T m φ ), mostly into neutrinos, with a small fraction decaying into SνDM.
Before z A , neutrinos and photons obey T ν = T sc ν = T γ . The ν-φ interaction freezes in a sub-thermal population of φ via the Boltzmann equation where Y φ = n φ /s is the yield, n φ is the number density of φ, and s is the entropy density of the Universe. The thermally-averaged rate is defined as γ φ→νν ≡ n eq ν (m φ /E)Γ φ→νν where Γ φ = 3λ 2 m φ /(32π) is the decay rate at rest. Here we assume φ couples universally to all neutrinos, λ ee,µµ,τ τ = λ. In the small m φ limit, Assuming no φ population at early times, integrating Eq. (4) up to z A determines the population of φ and its contribution to ∆N bbn eff . From z A -z B , the ν-φ interaction induces chemical equilibrium. At z B , T ν deviates from T sc ν ; a chemical potential is also generated, µ φ = 2µ ν = 2µν for real scalar φ, or µ φ = 2µ ν and µ φ * = 2µν for complex φ. Without CP violation, µ φ = µ * φ . At z B , T ν and µ ν are dictated by energy/lepton-number conservation [48][49][50][51]. Entropy is 1 10 ≈ 0.57 not preserved here; the φ thermalization process is irreversible. In our supplemental material, we provide details of this matching. The most important epoch for SνDM production is from z B -z C , by the νν → νν 4 process with on-shell, schannel φ exchange. This is terminated after φ becomes non-relativistic and freezes out. Entropy/lepton-number conservation allow us to derive T ν (z) and µ ν (z). The total energy of the Universe increases because some φ decay while turning non-relativistic. The neutrino selfinteraction rate also depends on T ν (z) and µ ν (z), where ω = m 2 φ /(4ET ν ). With the above, we integrate Eq. (2) up to z C to obtain the relic abundance of SνDM, where s 0 = 2891.2 cm −3 and ρ 0 = 1.05 × 10 −5 h −2 GeV/cm −3 are entropy and critical energy densities today [52]. The total entropy density at time z C (we take z C = 30) is where T γ 1.39T sc ν = 1.39m φ /z C , and The ν-φ dynamics also allow one to calculate ∆N bbn eff and ∆N cmb eff . A positive ∆N bbn eff is generated by the nonthermal production of φ at temperatures above the MeV scale, where ρ ν (z A ) = 9ζ(3)T ν (z A ) 4 /(4π 2 ) and ρ φ (z A ) can be calculated by solving Eq. (4). The thermal reaction rate γ φ↔νν is proportional to (λm φ ) 2 , implying weaker BBN constraints for light φ. By CMB times, φ particles have decayed away. The contribution to ∆N cmb eff is due to the non-standard distribution function of active neutrinos, 2 where In short, neutrino self interactions in the early universe via a light mediator φ imply a tight correlation between the origin of SνDM and the amount of extra radiation. The post-BBN thermalization of φ provides a novel mechanism that accounts for the correct relic density of SνDM. Meanwhile, φ thermalization distorts the active neutrino phase-space distribution, causing ∆N eff > 0 for both CMB and BBN. It is then of great interest to quantify this interplay, especially the prediction for ∆N cmb eff , which will serve as a well-motivated target for the upcoming CMB-S4.
Results and Discussions -We select three benchmark SνDM masses, m 4 = {2.5, 4, 7.1} keV, with the correspondingly largest active-sterile mixing parameter, sin 2 2θ 5 × 10 −9 , 10 −9 , 7 × 10 −11 , consistent with current X-ray search limits [2,[53][54][55]. In Fig. 2, we scan over (m φ , λ), focusing on m φ between keV and MeV, and derive the black solid curves that can account for the observed DM relic abundance, Ω ν4 h 2 = 0.1186. For m φ m 4 , all the curves follow λ ∼ m 1/2 φ , matching the result for m φ 1 MeV [3]. As m φ → m 4 , the curves bend up, where larger coupling λ is required, compensating for the phase-space suppression of the φ → νν 4 decay. We also present contours of constant ∆N bbn eff and ∆N cmb eff generated by the φ-ν interaction. Their intersections with the black contours dictate the correlation between Ω 4 and ∆N eff . For real (complex) scalar φ, we find the following predictions 0 < ∆N bbn eff < 0.57 (1.14) and 0.12 (0.19) < ∆N cmb eff < 0.9 (1.95). The maximal values correspond to strong couplings with φ fully thermalized by BBN, whereas the minimal values correspond to sufficiently small λ and/or m φ and negligible pre-BBN population of φ. In the lower-left half of each λ-m φ plane, both ∆N bbn eff and ∆N cmb eff follow the same parametrical dependence as the φ production rate γ φ↔νν in Eq. (4) and depend only on the product (λm φ ).
Thanks to the late thermalization of φ, the predicted ∆N cmb eff values are always higher than ∆N bbn eff . However, we find ∆N bbn eff is not always negligible and it serves as an independent test of non-standard neutrino cosmology. The prior-to-BBN production of φ particles can increase the expansion rate of the universe, leading to a higher neutron-to-proton ratio when BBN begins, and thus higher helium and deuterium abundances [17]. The state-of-the-art limit, ∆N bbn eff 0.5 [36], marginally excludes the case of fully-thermalized real scalar φ prior to BBN and sets a stronger constraint for complex φ. Note that the scalar can only mediate neutrino selfinteractions, ensuring that weak interactions are unaffected, and hence heat transfer from the plasma to the neutrino sea is only indirectly impacted. The changes produced in the primoridial abundances of helium and deuterium are well within observational uncertainties, and hence the state-of-the-art limit applies [56]. Precision calculations of neutrino decoupling taking into account of entropy transfer from the plasma to ν (and φ) could improve the accuracy of the ∆N bbn eff prediction to O(1%) level [57,58].
From Fig. 2 precision at the few-percent level. The two classes of probes are highly complementary and will further restrict or discover SνDM.
The range of SνDM mass shown in Fig. 3 is 2 keV < m 4 < 10 keV. Lower values of m 4 are inconsistent with the phase space density constraint from dwarf galaxies [2,60,61]. A stronger lower bound on m 4 may be set using Lyman-α forest [62] and Milky Way satellite dwarf galaxy counts [63]. In our scenario, the DM freestreaming length can be substantially shorter than regular warm DM. In numerical calculation, we note that SνDM is dominantly produced around z ∼ 3 − 4 where φ already starts to turn non-relativistic (see Fig. 2 in supplemental material A). The initial DM velocity is thus suppressed if m φ ∼ m 4 . Interestingly, such a neardegeneracy exactly occurs where ∆N eff is minimized (see Fig. 2). As a result, small scale constraints are relaxed.
So far we have assumed that φ couples equally to all active neutrinos. Our conclusions qualitatively hold if this assumption is relaxed to generic flavor non-universal and/or off-diagonal couplings. By the time of φ thermalization, there is sufficient time for neutrino flavor conversions to occur, which necessarily equilibrates φ with all the active neutrinos.
In summary, we explored the impact of active neutrino self-interaction, mediated by an ultralight scalar φ, on the relic density of SνDM and the amount of extra radiation (∆N eff ) in the early Universe. Strong neutrino self-interactions can facilitate efficient production of DM and generate a net ∆N eff at the same time. Our study reveals an intimate relationship between these two impor-   tant quantities and provides a well-motivated target for the future cosmological and astrophysical experiments, including CMB-S4 and Athena. These results, together with earlier work, depicts the complete parameter space for SνDM, with neutrino self-interactions mediated by a scalar whose mass spans over more than six orders of magnitude. They serve as well a motivated target for future experimental tests, from laboratories to the cosmos. Observations made in this work could extend to other neutrino self-interaction models mediated by vector bosons, such as gauged U (1) Lµ−Lτ [4,67,68] and other anomaly free gauge theories of baryon and lepton number.

A. Matching Conditions
In this section, we provide the details on the matching conditions among the three time scales defined in the main text, z A , z B , z C . This allows us to derive the time (z) dependence of the active neutrino temperature T ν (z) and chemical potential µ ν (z).
At time z A , the active neutrinos have just decoupled from the thermal plasma of other SM particles, thus T ν (z A ) = T sc ν (z A ) 1 MeV and µ ν (z A ) = 0. The corresponding energy and number densities of active neutrinos are Here we count neutrinos and antineutrinos separately, i.e., ρν = ρ ν , nν = n ν . Meanwhile, there also exist a sub-thermal population of φ that was frozen in via its interaction with active neutrinos. The corresponding values, denoted by ρ φ (z A ) and n φ (z A ), can be obtained by numerically solving Eq. (5) of the main text. The time z B is defined when φ has just entered chemical equilibrium with active neutrinos. At this moment, the φ-ν system shares the same temperature T ν (z B ), and also develops chemical potential where µ φ (z B ) = 2µ ν (z B ). The corresponding energy, number, and en-tropy densities are where we treat φ as an ultra-relativistic species (m φ T ν ), which is a good approximation throughout the parameter space we explore in this work. Again, in the case where φ is a complex scalar, the densities of the φ * degree of freedom are counted separately. In the absence of CP violation, The matching conditions between z A and z B are dictated by energy and lepton number conservation [48]. If φ is a real scalar, it comes to chemical equilibrium with both ν andν, thus Alternatively, if φ is a complex scalar, it only comes into equilibrium with neutrinos (φ * equilibrates withν), thus By definition z A /z B = T sc ν (z B )/T sc ν (z A ), and noting that T ν (z A ) = T sc ν (z A ), the above set of equations (14) or (15) allows us to derive the ratios In particular, in the case of real φ and with a negligible population of φ at time z A (corresponds to a tiny λ), we can derive consistent with the findings of Ref. [50].
As the next step, we consider the evolution of T ν and µ ν between z B and z C , where z C is defined as a temperature much lower than m φ . During this epoch, the ν-φ system remains in chemical equilibrium. As a result, we can describe their thermal distributions at any time z where z B < z < z C . Because neutrino masses are still negligible, their energy, number and entropy densities are described by the same functions as those in Eq. (13), by simply replacing z B → z. On the other hand, the mass of φ becomes important and the corresponding distribution functions are evaluated using can derive the ratio of neutrino temperatures in this model to that in standard cosmology, T ν (z)/T sc ν (z). For example, in the limit where z = z C 1 such that ρ φ , n φ , s φ → 0 are all Boltzmann suppressed, for real scalar φ and negligible population of φ at z A , we obtain also consistent with Ref. [50]. In Fig. 4, we show the evolution of T ν and µ ν as functions of z for three sets of model parameters. For very small λ, the asymptotic values are consistent with the results in Eq. (17) and (22). For sufficiently large λ, the freeze in population of φ (via Eq. (5) of the main text) before BBN is non-negligible which affects the matching results from z A to z C . The trend is also consistent with the limit where φ already thermalizes before BBN. In that case, we would have z A = z B and µ ν (z B ) = 0.
With the above results, we calculate the SνDM relic density using Eq. (7) in the main text. In Fig. 5, we plot the time dependence of dΩ ν4 /d ln z for a few sets of model parameters. In all cases, we find that SνDM is dominantly produced around z 3 − 4 where φ already starts to become non-relativistic. If φ and ν 4 is nearly degenerate (m φ m 4 ), SνDM will be produced with a velocity v c. Furthermore, the above results also allow us to calculate ∆N eff at CMB time. After φ decays away at z C 1, the ratios T ν (z C )/T sc ν (z C ) and µ ν (z C )/T ν (z C ) approach to constant values throughout the later evolution of the Universe, as shown in Fig. 4. We plug these functions into Eq. (11) of the main text. In the case of negligible φ population by the time of BBN, the value of ∆N cmb eff is 0.12 and 0.19 for real and complex φ case, respectively. The result for the real scalar case is consistent with the findings in Refs. [48].
We observe that ∆N bbn eff and ∆N cmb eff depend on model parameter space λ and m φ only through the combination λm φ as long as the product is small (corresponding to the low-left half of Fig. 2 of the main text). This dependence and the parametrical correlation are shown in Fig. 6 for real and complex scalars. We see that for sufficiently small λm φ , there is no impact on ∆N bbn eff (blue curves) -no abundance of φ has developed by the time of BBN. On the contrary, even for λm φ ≈ 10 −9 keV, ∆N cmb eff has a lowest value, 0.12 or 0.19, depending on whether φ is real or complex due to distortions of the active neutrino distributions. For both BBN and CMB, large enough λm φ results in a complete saturation of ∆N eff and these curves level off to their maxima. This approximation is only valid for m φ 1 MeV, where φ cannot decay away before BBN.

B. Minimal ∆N eff
Here we provide some further context for the procedure used in generating Fig. 3 of the main text. For a given m 4 , we determine the largest allowed sin 2 (2θ) assuming current limits from X-ray searches [2,[53][54][55], as well as the projected upper limit by Athena [59]. For a concrete choice, let us consider m 4 = 2.5 keV. The current limit on this mass is sin 2 (2θ) < 5 × 10 −9 , whereas Athena will have sensitivity at the level of sin 2 (2θ) = 10 −10 . We determine, for those two values of sin 2 (2θ), the preferred region of λ vs. m φ parameter space for which ν 4 can account for all of the DM observed today -this is shown as the two solid black lines in Fig. 7. Each of those curves, near the minimum in λ realized, attains some minimum value of both ∆N cmb eff and ∆N bbn eff . For the case we have considered here, m 4 = 2.5 keV, we find that with current limits on sin 2 (2θ), ∆N bbn eff (∆N cmb eff ) must be at least 0.02 (0.14), and with future constraints (assuming Athena does not observe a signal compatible with m 4 = 2.5 keV) will require ∆N bbn eff > 0.4 (∆N cmb eff > 0.6). In Fig. 7 we have assumed φ to be a real scalar, but the procedure is the same for complex φ. These minimum values are what is shown for that choice of m 4 in the left (current constraint on sin 2 (2θ)) and right (post-Athena constraint) panels of Fig. 3.