A Thermal Neutrino Portal to Sub-MeV Dark Matter

Thermal relics lighter than an MeV contribute to the energy density of the universe at the time of nucleosynthesis and recombination. Constraints on extra radiation degrees of freedom typically exclude even the simplest of such dark sectors. We explore the possibility that a sub-MeV dark sector entered equilibrium with the Standard Model after neutrino-photon decoupling, which significantly weakens these constraints and naturally arises in the context of neutrino mass generation through the spontaneous breaking of lepton number. Acquiring an adequate dark matter abundance independently motivates the MeV-scale in these models through the coincidence of gravitational, matter-radiation equality, and neutrino mass scales, $(m_\text{Pl} / T^\text{MRE})^{1/4} \, m_\nu \sim \text{MeV}$. This class of scenarios will be decisively tested by future measurements of the cosmic microwave background and matter structure of the universe. While the dark sector dominantly interacts with Standard Model neutrinos, large couplings to nucleons are possible in principle, leading to observable signals at proposed low-threshold direct detection experiments.

The mass of dark matter (DM) is relatively unconstrained. Demanding that its de Broglie wavelength is smaller than the typical size of Dwarf galaxies requires m DM 10 −22 eV, while microlensing searches for massive composite objects imply that m DM 10 58 GeV [1][2][3]. However, if DM acquired its abundance through thermal contact with the Standard Model (SM) bath, the viable mass range is significantly reduced and a much sharper picture emerges. For concreteness, we define "thermal dark matter" in this manner: thermal dark matter: dark matter that acquired its cosmological abundance after entering thermal equilibrium with the Standard Model bath at temperatures much higher than the freeze-out temperature of numberchanging interactions.
The canonical example of this scenario is embodied by the Weakly Interacting Massive Particle (WIMP) paradigm, in which DM is assumed to be in thermal contact with the SM bath while relativistic before chemically (and later kinetically) decoupling from the SM while non-relativistic. For m DM keV, thermal DM is sufficiently cold such that the free-streaming length in the early universe does not suppress the growth of matter perturbations on scales larger than the observed structures in intergalactic gas [4,5]. For larger masses, perturbative unitarity requires m DM 100 TeV under the assumption of a standard thermal cosmological history [6]. Thus, the thermal DM paradigm drastically restricts the possible mass range.
Although no theoretical inconsistencies arise for small masses, m DM MeV is often quoted as a robust lower bound on the mass of any thermal relic [7][8][9][10][11][12][13][14]. Such limits are usually derived from indirect measurements of the expansion rate of the universe in the radiation-dominated epoch, which can be parametrized in terms of the effective number of neutrino species, N eff . Sub-MeV thermal DM is relativistic at the time of nucleosynthesis and can modify N eff . However, the successful predictions of standard Big Bang nucleosynthesis (BBN) and observations of the cosmic microwave background (CMB) constrain N eff to lie near the SM expectation, N eff 3.046 [15,16].
As originally pointed out in Refs. [17][18][19] and recently studied in the context of light DM in Ref. [20], constraints on sub-MeV relics can be alleviated if equilibration between the DM and SM sectors occurs after neutrinos have already decoupled from the photon bath. As we will argue below, this process of delayed equilibration is characteristic of thermal DM that is much lighter than a GeV. In this work, we investigate a concrete and predictive model in which this scenario naturally arises for DM thermally coupled to SM neutrinos. There has been resurged interest in models of light thermal DM that interacts with neutrinos [21][22][23][24][25][26], which has largely been driven by the fact that such interactions constitute a simple mechanism to evade strong constraints from late-time distortions of the CMB [27].
Although our investigation is warranted solely as a proof of concept for sub-MeV thermal relics, the consideration of such models is timely. Various experimental technologies have recently been proposed for the direct detection of thermal DM down to the keV-scale [28][29][30][31]. However, below an MeV, the landscape of cosmologically viable models that will be tested by these experiments is rather unclear and under-explored (see Refs. [13,32] for detailed investigations of some simplified models). While the most minimal versions of the models examined in this work do not give rise to observable signals at these low-threshold detectors, variations upon these scenarios yield detectable rates. We will investigate this in more detail towards the end of this work. Furthermore, as we will discuss below, our setup will be definitively tested by upcoming cosmological observations, such as CMB-S3/S4 (and to some degree 21-cm) experiments.
The remainder of this paper is structured as follows. In Sec. II, we review the standard considerations of sub-MeV thermal relics as studied in previous literature. We then discuss in detail how the standard constraints can be alleviated in a model-independent manner in Sec. III. In Sec. IV, we introduce a simple concrete model motivated by the observed masses and mixing angles of the SM neutrinos. These models Secs. V and VI. We briefly summarize our results and conclusions in Sec. VII. A more detailed discussion on some aspects of the model is presented in Appendix A.

II. REVIEW OF SUB-MEV THERMAL RELICS
In this section, we discuss the physics of light relics and their effects on the measurements of primordial light element abundances and the CMB. For the models considered in this work, the main impact of the new degrees of freedom is through their contribution to the Hubble expansion rate, H 8π 3 where m Pl 1.22 × 10 19 GeV is the Planck mass and we have assumed that the energy content of the universe is dominated by the radiation component, ρ rad . The radiation energy density includes contributions from SM particles (γ, e ± , ν) and the dark sector. It is conveniently parametrized by the effective number of neutrino species, N eff , such that ρ rad ≡ ρ γ 1 + (7/8) (ξ SM ν ) 4 N eff (T ) , where ξ SM ν (T ) = T SM ν /T γ is the neutrino-to-photon temperature ratio in the standard cosmology (see Eq. (8) below). Thus, N eff is simply the neutrino and dark sector contribution to the total radiation energy density, normalized to the photon bath. In contrast to the common definition of N eff as a late-time quantity (only to be evaluated at the time of recombination), N eff (T ) in Eq. (2) parametrizes the expansion rate at temperatures below a few MeV. N eff can be modified either by changing the actual number of degrees of freedom in the radiation bath or by altering T ν /T γ . The notation for these and other relevant temperature scales is compiled in Table I for convenience.
Novel evolution of N eff (T ) can modify the predictions of primordial nucleosynthesis and recombination.
The outcomes of these cosmological epochs have been precisely measured and therefore constrain nonstandard behavior of N eff . Below, we summarize the effects of varying N eff on aspects related to BBN and the CMB and then review how light dark sectors can run afoul of the resulting constraints.
Observations of the CMB power spectrum are also sensitive to the total radiation energy density at the time of recombination. Detailed analyses of this effect are presented in Refs. [40,41]. We summarize their arguments below. CMB temperature anisotropies on scales smaller than the diffusion length of photons at recombination are exponentially damped, a mechanism known as Silk or diffusion damping [42]. On the microscopic level, this corresponds to the stochastic process of photons Thomson scattering with free electrons. Hence, the diffusion distance, r d , can be written parametrically as r d ∼ √ N λ mfp , where N is the number of scatters, λ mfp ∼ 1/n e σ T is the photon mean free path, n e is the free electron number density, and σ T is the Thomson cross section. The diffusion length scale is therefore r d ∼ 1/Hλ mfp λ mfp ∼ 1/Hn e σ T . is dictated by the corresponding angular scale, θ s = r s /D A . Hence, the ratio of angular scales θ d /θ s = r d /r s ∼ H/n e σ T is independent of D A . The position of the first peak has been measured to a precision of 5 × 10 −4 [27]. Thus, fixing θ s to the observed value, the scaling argument above implies that larger N eff (and hence H) leads to a larger θ d , thereby suppressing power in the damping tail of the CMB. Note that the degree of damping at small angular scales is increased for larger N eff , even though the underlying physical diffusion length is decreased. This behavior is seen explicitly in full Boltzmann simulations [40,41]. The argument above also makes explicit the degeneracy between N eff and Y p ; since n e ∝ 1−Y p , the effect on r d /r s from decreasing Y p can be compensated by increasing N eff . This degeneracy is broken by considerations of BBN.
Measurements by the Planck satellite constrain the effective number of neutrino species at the time of last scattering with unprecedented precision, N eff 3.15 ± 0.23 at 68% confidence [27]. Although the inclusion of different cosmological datasets modifies this result slightly, we will take this value as a representative benchmark in our analysis. A recent direct measurement of the local Hubble constant, H 0 , is in tension with the inferred value from Planck data at the level of ∼ 3.4σ [43]. The inclusion of additional relativistic species at the time of recombination significantly alleviates the tension, favoring ∆N eff 0.4 [43][44][45][46]. This is not the case when the "preliminary" Planck measurements of high-polarization are included, which favor a standard cosmology, but it is possible that this dataset is plagued by low-level systematics [27,47].

C. Standard Light Relics
Neutrinos decouple from the photon bath at a temperature of T ν dec ∼ 2 MeV [48]. A set of sub-MeV hidden sector (HS) particles (collectively denoted as X) that is equilibrated with the SM at temperatures below T ν dec can lead to significant deviations in the observed value of N eff . The lightest stable particle of this HS constitutes the DM of the universe. For simplicity, we assume that X couples to the SM neutrinos and that all such particles have a common mass given by m X . We first consider the standard case where X equilibrates with the SM neutrinos before the point of neutrino-photon decoupling, as has been investigated in Refs. [7][8][9][10][11][12][13][14]. The temperature evolution of the neutrino bath is then easily derived from the conservation of comoving entropy density.
The effective number of relativistic degrees of freedom, g i * , in each bath (i = ν, X, γ) determines the entropy density, s i ≡ (2π 2 /45) g i * T 3 i , and energy density, ρ i ≡ (π 2 /30) g i * T 4 i , where T ≡ T γ . For three generations of left-handed SM neutrinos, At temperatures below T ν dec , the comoving entropy densities in the ν − X and photon bath are separately conserved. Using that s ν+X ≡ s ν + s X and s γ separately scale as a −3 (a is the scale factor), one finds where is the temperature of species i normalized to the photon temperature [49]. Treating electron-photon decoupling as instantaneous, we can approximate the number of relativistic degrees of freedom coupled to the photon bath as g γ * (T m e ) = 2 + (7/8) × 4 = 11/2 and g γ * (T m e ) = 2. Equating Eq. (4) at temperatures above and below m e , and using that ξ ν (T m e ) = 1, one recovers the standard result ξ ν (T m e ) 4 11 density of the ν − X bath. Again using Eq. (4), but for T ν m X and T ν m X , we find that For later convenience, we define ξ SM ν as the value of ξ ν assuming a standard cosmology (g X * = 0) such that Using the above results, the defining expression for N eff in Eq. (2) can be rewritten as In Eq. (9), we have assumed that X decouples instantaneously once its temperature drops below its mass (T X m X ), which is encapsulated by the Heaviside step function, Θ [49,50]. Note that Eq. (9) reduces to N eff 3 when g X * = 0 and ξ ν = ξ SM ν . In the SM, neutrino decoupling is not instantaneous, and e ± annihilations partially heat the neutrino bath, resulting in N eff 3.046 [15,16]. In Eq. (9), we have approximated 3.046 3. Substituting Eqs. (7) and (8) into Eq. (9), we find that if X equilibrates with the SM neutrinos at temperatures above T ν dec . If eV m X MeV, then Eq. (10) gives N eff 3.57 (N eff 3.79) at the time of nucleosynthesis (recombination) for g X *

As discussed in
Secs. II A and II B, this is excluded from considerations of BBN and Planck measurements of the CMB by more than 2σ. Furthermore, realistic models of light thermal DM often require g X * few, leading to even larger deviations in N eff . It is this basic insight that has driven many studies to claim that sub-MeV thermal DM is not cosmologically viable [7][8][9][10][11][12][13][14].

A. Temperature Evolution and Effective Number of Neutrino Species
In Sec. II, we noted that a single sub-MeV degree of freedom that is equilibrated with the SM below the temperature of neutrino-photon decoupling, T ν dec ∼ 2 MeV, can lead to deviations in N eff that are in conflict with considerations of BBN and the CMB. In this section, we illustrate that if light relics enter equilibrium with the SM at temperatures below T ν dec , then such constraints are significantly relaxed [17][18][19][20].
Let us assume that a similar collection of sub-MeV particles (X) equilibrates with the SM neutrino bath while relativistic but after neutrino-photon decoupling. The assumption of relativistic equilibration is not strictly necessary, but simplifies the estimates below (see Sec. III C). As summarized in Table I,  In the standard WIMP framework (red), dark matter is assumed to be in equilibrium with the Standard Model bath long before freeze-out. Dark matter produced through freeze-in (yellow) is assumed to have a negligible abundance at early times and never fully equilibrates with the Standard Model. We propose a scenario (blue) that alleviates strong constraints from measurements of the effective number of neutrino species and is much more akin to the WIMP paradigm, in which an initially cold (compared to the photon bath) population of sub-MeV particles relativistically equilibrates with the Standard Model bath after neutrino-photon decoupling and before freeze-out.
Similar behavior is also expected for standard WIMPs, although the temperature at equilibration (T X eq ) is typically much larger.
define T X eq m X and T X dec ∼ m X as the temperature of the photon bath at which X enters and exits equilibrium with neutrinos, respectively, and T BBN ∼ (10 − 50) keV as the temperature at which nucleosynthesis has effectively concluded. We will be interested in the case where the HS is initially colder than the SM bath. A schematic representation of the cosmological evolution of the HS comoving number density is shown in Fig. 1. Contrary to DM that is produced via freeze-in [51], we assume that the HS is fully relativistic while equilibrating with the SM, analogous to the thermal history of a standard WIMP.
For concreteness, we assume that X equilibrates with the SM neutrinos after neutrino-photon and electronphoton decoupling, i.e., T X eq T ν dec , m e ∼ MeV. An example of the temperature evolution of the neutrino and HS baths is shown in Fig. 2. These results were obtained by numerically solving the Boltzmann equations for the X and ν energy densities. Analytic approximations will be derived below. If HS-SM equilibration occurs through decays and inverse-decays of a HS species into neutrinos (X ↔ νν), then the relevant Boltzmann equations arė where Γ dec is the decay rate for X → νν, the superscript "eq'" denotes an equilibrium distribution, and (blue) sectors for an initial temperature ratio of ξ 0 X = 0.3. Compared to standard cosmology, neutrino-dark matter equilibration and decoupling cools and heats the neutrino population relative to its expected value in the Standard Model, respectively. The horizontal gray dashed lines correspond to the approximate analytic estimates of Eqs. (16) and (19). (Right) Evolution of the effective number of neutrino species in the case that dark matter equilibrates with neutrinos after (solid blue) or before (dotted blue) neutrino-photon decoupling. The horizontal gray dashed lines correspond to the approximate analytic estimates given in Eqs. (17) and (20). For concreteness, we have taken the hidden sector to be made up of a 10 keV Majorana fermion and a 5 keV real scalar.
we have been explicit at which temperature the equilibrium number/energy densities should be evaluated.
In writing the above equations, we have neglected Bose-enhancement and Pauli-blocking factors. Including these effects modifies the collision term by O(1) factors, but does not significantly change our results. Eq. (11) can be solved numerically for the evolution of T X,ν as a function of the photon temperature, T . The time variable can be traded for the photon temperature through the relation [52] where ρ tot ≡ ρ γ + ρ ν + ρ X and similarly for the pressure density, P tot . In Eq. (11), we have neglected chemical potentials, assuming that the interactions between the HS and neutrino baths enable each species to rapidly track equilibrium distributions dictated by T X,ν . This is a good approximation for the model described below in Sec. IV, since chemical potentials are suppressed by number-changing reactions involving a light spin-0 mediator. In particular, for O(1) couplings and keV-scale masses in the HS scalar potential of Sec. IV C, 4 → 2 self-interactions involving the spin-0 mediator decouple well after DM freeze-out.
In the left panel of Fig. 2, we show the cosmological evolution of the neutrino and HS temperatures normalized to that of the photon bath as solid red and blue lines, respectively, assuming that the HS consists of a 10 keV Majorana fermion and a 5 keV real scalar. For comparison, we also display the temperature evolution of the neutrino bath in the SM (dotted red), assuming that no new light thermal relics are present (g X * = 0). The initial HS-SM temperature ratio is fixed to ξ X = 0.3, such that the HS is initially much colder than the SM neutrino and photon populations. Energy conservation then implies that ν − X equilibration cools (heats) the neutrino (X) bath at T ∼ T X eq . If this occurs after neutrino-photon decoupling, this leaves the photon bath unaffected. Later, when the temperature drops below m X and the HS decouples, X dumps its entropy back into the neutrinos, reheating them to a temperature slightly above the SM expectation.
These two processes, equilibration (neutrino cooling) and decoupling (neutrino heating), have counteracting effects on the neutrino temperature, which lead to a partial cancellation and a significant reduction in modifications to N eff , whose evolution is shown as the solid blue line in the right panel of Fig. 2. If the HS is initially colder than the SM bath, this cancellation is a direct consequence of thermodynamics and does not constitute a tuning of the model. For comparison, we also show the temperature evolution of N eff , taking the standard assumption that equilibration occurs before neutrino-photon decoupling (dotted blue), as in Sec. II C and Refs. [7][8][9][10][11][12][13][14]. If equilibration occurs after neutrino-photon decoupling, deviations in N eff are significantly reduced. We now derive analytic approximations for the asymptotic behavior of ξ ν,X and N eff , which are shown as the horizontal gray dashed lines in Fig. 2.
As we will soon see, N eff is sensitive to the initial value of ξ X ≡ T X /T before X − ν equilibration or electron-photon decoupling, but, similar to DM production via freeze-in, it is insensitive to the particular value of ξ X as long as ξ X 1 [51]. We define ξ 0 X ≡ ξ X (T T X eq , m e ) as this initial temperature ratio.
As mentioned above, for simplicity, we assume that electron decoupling occurs before DM equilibration.
Comoving entropy is conserved as electrons decouple from the photon plasma. Electron annihilations heat photons relative to the neutrino and X baths. Hence, as in Sec. II C, for T X eq T m e , we have Along with Eq. (9), this implies that N eff is given by This is the standard result for an uncoupled population of dark radiation.
If the HS and neutrino baths equilibrate while X and ν are relativistic, the sum of their comoving energy densities, ρ ν+X a 4 , is approximately conserved. This can be seen from Eq. (11), which implies that d ρ ν+X a 4 /dt = ρ ν+X a 4 H (1 − 3 w), where w ≡ P ν+X /ρ ν+X . When T ν m ν and T X m X , we have w 1/3 and d(ρ ν+X a 4 )/dt 0. Therefore, before and immediately after X − ν equilibration, where we have used s γ ∝ a −3 . Equating this expression at temperatures above and below T X eq , we find ξ νX (T X dec T T X eq ) 4 11 where ξ νX ≡ ξ ν = ξ X is the temperature ratio when X is equilibrated with the SM neutrino bath. Comparing the above expression to the standard result of Eq. (8), we see that for ξ 0 X 1, ν−X equilibration significantly lowers the temperature of the neutrino bath, i.e., ξ νX ξ SM ν . Eqs. (9) and (16) then imply that during X − ν equilibration and before X becomes non-relativistic. Note that Eq. (17) is identical to the expression of Eq. (14). This is consistent with the fact that d(ρ ν+X a 4 )/dt 0 and that N eff is defined in terms of the total radiation energy density.
We use conservation of entropy when X becomes non-relativistic and decouples, since this process occurs in equilibrium. Hence, just before and after X becomes non-relativistic. Equating this expression above and below T X dec ∼ m X and using Eqs. (16) and (9), we find ξ ν T T X dec 4 11 and Note that in the ξ 0 X 1 limit and taking T X dec ∼ m X , Eqs. (14), (17), and (20) reduce to N eff (T m X ) 3 (21) and where in the inequality we have imposed g X * 1 for any light HS.
Compared to the standard result of Eq. (10), the deviation in N eff away from its SM expectation is significantly reduced in Eq. (20) for ξ 0 X 1. As mentioned previously, if T X eq T ν dec , then ν − X equilibration drains the neutrino bath of energy, lowering its temperature compared to that of photons.
Later, when X becomes non-relativistic and decouples, it reheats the neutrinos to a temperature close to the SM expectation. These processes have counteracting effects on ξ ν , such that the neutrino bath is reheated to a smaller degree than if T X eq T ν dec . However, as seen from Eq. (19), even for ξ 0 X 0, there is an irreducible heating of the neutrino bath since equilibration of two initially decoupled gases leads to an overall increase in the comoving entropy of the ν − X system. In the left (right) panels of Fig. 2, the horizontal gray dashed lines correspond to the approximate values given by Eqs. (16) and (19) (Eqs. (17) and (20)). The numerical solutions are in good agreement with these approximate expressions, which warrants their use in the remainder of this work. We also note that a similar cancellation arises when a sub-MeV relic equilibrates directly with the photon bath after neutrino-photon decoupling, but we will not explore such models in this work.
Values of g X * (the effective number of sub-MeV dark sector states that equilibrate with neutrinos) and ξ 0 X (the initial dark sector-to-photon temperature ratio) compatible with the effective number of neutrino species at the time of nucleosynthesis (green) and recombination (blue). Regions compatible with BBN are shown for scenarios in which dark matter decouples from neutrinos before (T X dec T BBN ) and after (T X dec T BBN ) the end of nucleosynthesis.
We also highlight parameter space that alleviates the tension between Planck and local measurements of the Hubble parameter, H0. The representative model space (red) corresponds to a dark sector with a dark matter scalar or Majorana fermion and a scalar mediator. The vertical dashed gray line corresponds to the standard assumption that X equilibrates with neutrinos before neutrino-photon decoupling (ξ 0 X 1).
Equations (14), (17), and (20) imply that constraints from nucleosynthesis and the CMB can be alleviated if T X eq T ν dec and ξ 0 X 1. In Fig. 3, we highlight regions of parameter space in the g X * − ξ 0 X plane that are compatible with measurements of N eff . If T X eq T ν dec , then ξ 0 X = 1 in general and its value encapsulates the sensitivity of our setup to physics in the ultraviolet. For instance, if X was initially in thermal equilibrium with the SM but decoupled at T Λ QCD before reentering equilibrium at T T ν dec , then ξ 0 X ∼ (10/100) 1/3 ∼ 0.5. More generally, ξ 0 X = 1 arises in theories of asymmetric reheating of the DM and SM sectors [53]. Throughout this work, we take ξ 0 X to be a free parameter of the low-energy theory. Note that physics at low-energies is insensitive to this temperature ratio as long as ξ 0 X 1. This is analogous to the level of ultraviolet-sensitivity for DM produced from freeze-in processes, where one typically assumes a negligible initial DM abundance at early times [51].
For T X eq T ν dec , N eff transitions from Eq. (17) to Eq. (20) near the decoupling temperature, T X dec ∼ m X . As a result, limits from nucleosynthesis depend on the ordering of T X dec ∼ m X and T BBN ∼ (10 − 50) keV. Regions compatible with BBN are shown in Fig. 3 for both of the temperature orderings T X dec T BBN and T X dec T BBN . For T X dec T BBN , N eff is static during BBN and is given only by the expression in Eqs. (14) and (17). However, for T X dec T BBN , N eff evolves from the form given in Eqs. (14) and (17) to that of Eq. (20) during nucleosynthesis. Detailed studies of BBN, which demand N eff 2.85 ± 0.28 within 1σ, often assume a single fixed value of N eff throughout the entire formation of light nuclei [33]. However, as we have seen, this is not generally the case for a light HS that equilibrates and decouples from the SM during nucleosynthesis [54]. In deriving a constraint, we demand that N eff never deviates from the best-fit constant value by more than 2σ, i.e., |N eff (T ) − 2.85| ≤ 0.56 for T T BBN . We note that this is most likely overly conservative, since for ξ 0 X 1 and values of T X dec only slightly greater than T BBN , significant deviations in the expansion rate will only occur at the end of nucleosynthesis. For instance, this could potentially lead to slight changes in the deuterium or 7 Li abundance without affecting the production of 4 He. It would be interesting to consider the bounds from detailed investigations of BBN, while assuming time-variations of N eff in this manner. We leave such considerations to future work [55].
Cold DM is necessarily non-relativistic at the time of recombination, i.e., eV T X dec ∼ m X . To remain consistent with Planck measurements of the CMB within 2σ, we demand that |N eff (T ) − 3.15| 0.46 for T m X , where we take the form for N eff given in Eq. (20) [27]. Note that this CMB bound on N eff assumes standard nucleosynthesis, which is modified in the delayed equilibration scenario, as described above. A more realistic approach would be to fit both Y p and N eff to the CMB power spectrum. This can significantly expand the allowed parameter space due to the Y p -N eff degeneracy described in Sec. II B. Also shown in Fig. 3 are regions of parameter space that alleviate the tension between Planck and local measurements of the Hubble parameter, H 0 . As a representative favored range, we take N eff 3.4 ± 0.05 [43][44][45][46]. Models of light thermal DM require a stable species and a light mediator. We highlight regions of parameter space corresponding to the presence of two real scalars in the HS (g X * = 2), or a light Majorana fermion and a real scalar (g X * = 2.75). The standard case of T X eq T ν dec corresponds to the limit ξ 0 X 1, which is in strong tension with measurements of both the CMB and primordial nuclei abundances for g X * 1.

B. General Model-Building
We have demonstrated that constraints on sub-MeV thermal relics are weakened when the HS equilibrates with the SM after neutrino-photon decoupling. We would like to understand if this naturally occurs in models of light thermal DM. It has long been appreciated that thermal DM which couples to the SM solely through the electroweak force must be heavier than the GeV-scale. The so-called Lee-Weinberg bound relates the mass of thermal DM to the weak scale (m W ), the temperature at matter-radiation equality (T MRE ∼ 0.8 eV), [56]. Equivalently, thermal DM that is lighter than a GeV often requires the presence of new light mediators [57]. It is therefore natural to expect that sub-MeV thermal DM, denoted by χ, is accompanied by additional HS mediators, ϕ, that are nearby in mass. In this case, there are two processes that can equilibrate the two sectors: scattering between HS and SM states, and decays of ϕ into the SM. As we will show, the temperature dependence of either of these processes generically predicts that a light HS enters thermal equilibrium with the SM while relativistic. This is illustrated in Fig. 4. The equilibration point is independent of HS mass scales for scattering, but for decays, it occurs later as HS masses are lowered. If this proceeds at temperatures below a few MeV, the mechanism described in Sec. III A is realized and modifications to N eff during nucleosynthesis and recombination are reduced. At temperatures much greater than m χ or m ϕ , we parametrize the rate for scattering and decays/inversedecays as where α eq is the effective coupling governing equilibration and the factor of m ϕ /T in the decay rate is a time-dilation factor. Comparing either process to the Hubble parameter, H ∼ T 2 /m Pl , demonstrates that the rate for equilibration overcomes the expansion rate at temperatures below T X eq ∼ α eq m Pl (scattering) for scattering and decays, respectively, where T X eq denotes the temperature at which the DM and SM sectors equilibrate. If we parametrize the rate for DM annihilation during freeze-out as σv ∼ α 2 FO /m 2 χ , then χ acquires an abundance in agreement with the observed DM energy density for where α FO is the effective coupling governing freeze-out. Using this relation in Eq. (24) allows us to write m χ in terms of T X eq , Equation (26) implies that χ and ϕ equilibrate with the SM after neutrino-photon decoupling (T X eq Bounds on warm DM typically exclude m χ few × keV [4,5]. Therefore, m χ keV along with Eq. (27) motivates α FO α eq . This can be accomplished if the processes governing freeze-out are enhanced compared to those governing equilibration. This is a natural hierarchy, for instance, in models of secluded DM [58], those involving freeze-out through resonant annihilations [59], or strongly interacting hidden sectors [60].
Once χ and/or ϕ become non-relativistic, Γ scatt and Γ dec are either suppressed by Boltzmann or T /m χ,ϕ factors. At this point, the equilibration rate quickly drops below Hubble expansion and the HS decouples from the SM. This behavior can be contrasted with equilibration through the exchange of a heavy mediator, in which case the rate governing equilibration always falls faster in temperature than H ∼ T 2 /m Pl . This is typical of the weak processes that maintain ν-e equilibrium where Γ scatt ∼ G 2 F T 5 . Schematic examples of these scenarios are shown in Fig. 4.
The presence of light mediators is strongly motivated for sub-GeV thermal DM. Thermalization through these light mediators generically predicts that DM enters equilibrium with the SM while relativistic and before DM freeze-out, as highlighted in Fig. 4. If DM is sufficiently light and there exists a hierarchy between the couplings governing freeze-out and those governing scattering/decays, then the HS equilibrates with the SM after neutrino-photon decoupling, alleviating constraints from measurements of N eff . In Sec. IV, we turn our attention to a concrete model that explicitly realizes this mechanism. However, as an aside, we first briefly comment on scenarios in which the HS instead does not equilibrate with the SM bath until it is semi-or non-relativistic.

C. Non-Relativistic Equilibration
In the previous sections, we focused on a scenario that is closely related to the standard WIMP paradigm: the HS and SM baths are in equilibrium at temperatures much greater than the DM mass, with chemical decoupling from the SM occurring at temperatures much lower than the DM mass. This is to be contrasted with freeze-in production, in which case DM never fully equilibrates with the SM [51]. Although it is not the central focus of this work, an interesting situation may arise between these two extremes, where the HS fully equilibrates with the SM while the DM is semi-or non-relativistic, but before freeze-out of number-changing interactions. We briefly comment on this possibility here.
A few of these cosmological scenarios are shown in Fig. 5. The blue lines correspond to models in which DM fully equilibrates with the SM neutrino bath after neutrino-photon decoupling but well before thermal freeze-out. The cosmology denoted by the solid blue line was already discussed in detail in Sec. III A, in which the DM is relativistic during HS-SM equilibration. This case is most analogous to the WIMP paradigm, and simple analytic approximations for the evolution of the HS/neutrino temperatures and N eff were derived Sec. III A, we noted that the increase in N eff at late times is due to an irreducible heating of the neutrino bath since the equilibration of two initially decoupled gases leads to an overall increase in the comoving entropy of the ν − X system, i.e., where Q is the heat exchanged between the two sectors. If the HS is equilibrated to semi-or non-relativistic temperatures, instead of relativistic ones, the overall heat transfer and entropy increase are reduced, leading to a corresponding decrease in the overall heating of the neutrino bath once the HS becomes non-relativistic.
As a result, modifications to N eff at late times are suppressed compared to relativistic equilibration, as shown explicitly in Fig. 6. Although it is beyond the scope of this study, such models constitute an interesting possibility for light, predictive, thermal-like DM. In the next section and the remainder of this work, we will instead focus on an explicit realization of the cosmological scenarios involving relativistic equilibration, as discussed in Sec. III A.

IV. SUB-MEV DARK MATTER WITH A MAJORON MEDIATOR
The measurement of neutrino oscillations has firmly established the presence of neutrino masses and mixing amongst the different flavor eigenstates. Along with the gravitational observations of DM, the discovery of neutrino masses strongly motivates the existence of physics beyond the SM. We now outline a minimal model that realizes the mechanism described in the previous sections. This model generates the neutrino mass splittings and mixing angles, along with the parameters of the DM sector, through the spontaneous breaking of lepton number. In Sec. IV A, we discuss the basic framework that is needed to generate the appropriate parameters in the neutrino sector. In Sec. IV B, we extend the model to include a stable neutral lepton, which will play the role of DM. We briefly discuss the details of the Higgs sector in Sec. IV C. A more detailed discussion concerning the explicit forms for the masses and interactions of the HS particles is given in Appendix A.

A. Neutrino Sector
The SM lacks the necessary ingredients to explain the observed neutrino masses and mixing angles. A simple solution is to include the dimension-five Weinberg operator, (LH) 2 /Λ uv [61]. Below the scale of electroweak symmetry breaking, this operator generates neutrino masses parametrically of the form m ν ∼ v 2 /Λ uv , where v 246 GeV is the SM Higgs vacuum expectation value (vev) and Λ uv is the effective scale of new physics. A natural microscopic realization of this operator is the so-called seesaw mechanism, which introduces right-handed neutrinos that are uncharged under the SM gauge group [62][63][64][65][66]. If neutrinos are Majorana, then m ν = 0 breaks lepton number, U (1) L . The global U (1) L symmetry can be broken explicitly, as in minimal seesaw models with an explicit Majorana mass for the right-handed neutrinos, or spontaneously when a U (1) L -charged scalar acquires an expectation value. In the latter case, right-handed neutrino masses are generated dynamically, and the seesaw mechanism can be implemented. Such models involve majorons, the pseudo-Nambu-Goldstone bosons (pNGBs) of U (1) L [67][68][69]. This light pseudoscalar will play the role of the mediator between the visible and dark sectors.
In writing down the model, we follow the notation and conventions of Refs. [70] and [71]. We introduce a complex scalar, σ, of lepton number L = 2, where we have assumed that σ acquires a non-zero vev, σ = f / √ 2. S and J are the real and imaginary excitations of σ, where J (often dubbed the majoron) is the Goldstone boson of spontaneous U (1) L -breaking.
In the presence of suppressed terms that softly break lepton number, J is a pseudo-Goldstone and acquires a small mass. Soft U (1) L -breaking terms can arise in the scalar potential, which is examined in Sec. IV C and Appendix A 2. While we naturally expect m J f , we will not specify the exact form of U (1) L -breaking and treat the majoron mass, m J , as a free parameter of the low-energy theory. A discussion of how such masses may arise from gravitational effects in a more complete theory is provided in Appendix A 3.
We introduce three generations of right-handed neutrinos, N , with lepton number L = −1. The most general renormalizable and U (1) L -symmetric Lagrangian coupling σ and N to the SM lepton sector is then given by where m D ≡ y ν v/ √ 2 and M N ≡ y N f / √ 2 are 3 × 3 mass matrices. Diagonalizing M νN gives rise to the neutrino mass basis, n i (i = 1, 2, . . . , 6), with masses m i . We define the unitary matrix V that diagonalizes the full active-sterile neutrino mass matrix by where V relates the gauge and mass eigenstates. 1. This is made explicit by the Casas-Ibarra parametrization as discussed in Appendix A 1 [72]. The interactions of the neutrino mass eigenstates (n i ) with the scalar degrees of freedom take the parametric form where h is the SM Higgs field. The explicit forms of these couplings, along with ones involving SM gauge bosons, are given in Appendix A 1. The most important feature of the above interactions is their proportionality to the neutrino masses, which is characteristic of the Higgs mechanism. In general, there may be other contributions to the masses of the sterile neutrinos, for instance originating from Dirac masses with additional L = +1 sterile neutrinos. In this case, the mass parameters m 4,5,6 written in these interactions are implicitly assumed to be the piece given by the scale f , i.e., ∼ f × ∂m 4,5,6 /∂f . However, it is important to keep in mind that M N f is still possible in extended models. We will return to this point later in Mass-mixing in the neutrino sector also induces interactions of the sterile states with electroweak currents and generates couplings of S and J to charged leptons and quarks via neutrino loops. These interactions are typically too small to be phenomenologically relevant, but we discuss them briefly in Secs. IV C and VI C as well as in Appendix A.

B. Dark Matter Sector
The model described in the previous section involves a viable mechanism for neutrino mass generation.
The new particles include a naturally light pseudo-Nambu-Goldstone boson, J, that couples to neutrinos. This is precisely the setup required to realize a viable cosmology for sub-MeV DM as described in Secs. II and III. To complete the model, we introduce an additional Weyl fermion, χ, of lepton number L = −1 and charged under an additional Z 2 . The Z 2 prevents χ from mass-mixing with the active or sterile neutrinos and stabilizes χ, which will serve as our DM candidate. The only renormalizable term consistent with the above symmetries is The phase of χ can be chosen such that the Yukawa coupling, λ χ , is purely real. Below the scale of U (1) Lbreaking, χ acquires a mass, In four-component notation, the interactions of the Majorana fermion, χ, with J and S are given by where the couplings are defined as C. Scalar Sector The U (1) L -preserving renormalizable scalar potential is given by This potential does not generate a mass for the majoron, J. However, soft U (1) L -breaking terms such as can give rise to a radiatively-stable mass for J. The full potential is then given by We fix the phase of σ such that its vev, f , is real, leaving a single physical phase in the couplings µ σ and a σ . This phase leads to CP-violating mixing of J with S and h. The details of mass-diagonalization and constraints on the scalar potential parameters are discussed in Appendix A 2.
As we will illustrate in Sec. V, delayed equilibration of the majoron sector is achieved for m J m S m h .
We will assume that the mixing angles in the scalar sector are small, such that they do not significantly impact physics in the DM sector. Indeed, we will show in Sec. VI C and in Appendix A that the Higgs mixing with light states is strongly constrained by stellar cooling, rare meson decays, and Higgs decays, implying that the scalar mixing angles are suppressed. In this hierarchical limit, the scalar mass eigenstates (ϕ 1,2,3 ) are nearly aligned with the gauge basis (J, S, h), with masses This assumption will be relaxed in Sec. VI D when we consider possible signals in futuristic low-threshold direct detection experiments. For λ σ ∼ O(1), the mass of the CP-even scalar, S, is near the scale of U (1) L -breaking, m S ∼ f . For simplicity, we will fix m S = f in estimates and numerical results below.
We also note that tree-level mixing between J, S, and the SM Higgs, h, is not solely responsible for interactions between the HS and the electrically charged SM fermions. Additional contributions arise from diagrams involving loops of active/sterile neutrinos and electroweak gauge bosons. We will not discuss these contributions in detail and instead refer the interested reader to the relevant sections of Refs. [70] and [71]. For instance, the radiatively induced Yukawa couplings between J, S and the SM quarks and charged leptons are naturally of size, where f is a charged SM fermion. The effect of this coupling is analogous to S − h and J − h mass-mixing with an effective angle, θ eff , given by where we have taken m ν ∼ 0.1 eV. As a result, tree-level contributions to J, S − h mixing are only phenomenologically relevant for sin θ 10 −15 . As we will discuss below, the suppressed size of these radiative interactions makes them irrelevant for the physics governing early universe cosmology and the signals discussed in Secs. V and VI. We will come back to these couplings in Sec. VI C, where we discuss effects of J and S on the physics of stellar cooling.

A. Equilibration
In this section, we will discuss aspects related to the equilibration of DM with the SM. DM, χ, is assumed to equilibrate with the SM neutrinos, ν, while both sectors are relativistic. Since the majoron, J, is a pseudo-Goldstone of U (1) L , we naturally take m χ ∝ f m J . In this case, χ freezes out through annihilations into pairs of on-shell majorons, χχ → JJ, followed by J → νν, as shown in Fig. 7. From the interactions given in Sec. IV B, the non-relativistic cross section for this process is where v (not to be confused with the SM Higgs vev) is the relative DM velocity, and we have defined the mass ratio r J ≡ m J /m χ < 1. In Eq. (44), we have also taken the limit that m S f m χ , m J . This form suggests that χ acquires an abundance in agreement with the observed DM energy density for Hence, f ∼ 10 MeV − 1 GeV for m χ ∼ keV − MeV. In the minimal model described in Sec. IV, the masses of the sterile neutrinos and HS scalars are also governed by the U (1) L -scale, f , and therefore, we parametrically These parametric estimates for the relevant mass scales suggest that processes involving N , J, and S are all potentially relevant when considering equilibration between the DM and SM sectors. We will assume that rates for scattering processes in the HS, such as χχ ↔ JJ, are large compared to reactions involving both HS and SM species. Therefore, equilibration between the SM and a single species in the HS rapidly equilibrates all of the lightest particles in the HS, namely χ and J. As noted in Sec. III, sub-MeV thermal relics are viable provided that the HS equilibrates with the SM at temperatures below T ν dec ∼ 2 MeV. Therefore, it is imperative that processes involving SM neutrinos and N , J, and S do not equilibrate before this point.
We now proceed to discuss these various processes in detail.
In the limit that m J,S eV, the decay rates of J and S into SM neutrinos are where the sum is over the three active neutrino flavors. From examining the Boltzmann equations in Eq. (11), the effective energy transfer rates from decays and inverse-decays that can be compared to Hubble expansion are Γ X eq (J ↔ ν ν) m J n eq J (T ν ) ρ eq ν (T ν ) Γ(J → ν ν) where n eq J,S is the equilibrium number density of J, S, respectively [53,73]. These processes are able to maintain kinetic equilibrium between the HS and SM if Γ X eq (J, S ↔ ν ν) H. The ratio, peaks at temperatures comparable to the mass of the decaying particle, T ∼ m J,S . For concreteness, let us assume that m S m χ m J . We find that equilibration occurs at temperatures T ν m χ through J ↔ νν or through S ↔ νν decays if In the above estimates, we have set m ν ∼ 0.1 eV, m S ∼ f , and have fixed f to the thermally-favored value in Eq. (45). Eqs. (49) and (50) imply that for m χ ∼ keV − MeV and m J 10 −2 m χ , equilibration through J ↔ νν dominates over S ↔ νν.
Decays of the sterile neutrinos (N ↔ Jν) are also potentially able to equilibrate the two sectors. For the simplest choices of mixing parameters (R = 1 in Appendix A 1), each generation of N couples to a single generation of ν. For M N m J , the corresponding decay rate is Γ X eq (N ↔ J ν) is given by the analogous form of Eq. (47). We find that these decays efficiently equilibrate the DM and SM sectors if where, once again, we have fixed f to the thermally-favored value in Eq.
After fixing f to the cosmologically-favored value in Eq. (45), this implies that Jν ↔ Jν never maintains equilibrium between DM and the SM for T RH TeV × (m χ /100 keV) 2 .
Other scattering processes include χν ↔ χν through J and S exchange, Jν ↔ Zν, and St ↔ ht, where t is the SM top quark. We find that the rates of equilibration for these reactions are subdominant compared to the ones considered above since they are suppressed by additional small couplings. i.e., χχ → JJ, followed by J → νν (see Fig. 7). For convenience, we repeat the form for the non-relativistic cross section from Eq. (44), In calculating the relic abundance of χ, we follow the semi-analytic approach as detailed in Refs. [50] and [74], where ξ X is evaluated at freeze-out, b ≡ σv/v 2 as in Eq. (54), and As before, X collectively denotes the light species in the HS (χ and J). x f is the value of x ≡ m χ /T at freeze-out, and can be solved numerically through the relation where c ∼ O(1) is a constant chosen by matching to numerical solutions of the Boltzmann equation.
In Fig. 8, we show the value of f as a function of the DM mass, m χ , that is needed for an adequate freezeout abundance of χ, assuming that the HS equilibrates with the SM neutrinos at relativistic temperatures, i.e., T m χ , m J . We have taken m χ > m J , and the thickness of the contour in Fig. 8    After fixing m 1 , the masses of the other SM neutrinos are given by the observed mass splittings [75,76].
The regions in Fig. 9 were obtained by solving the Boltzmann equations in Eq. (11) to find T X eq . The qualitative behavior can also be obtained by comparing the rate of J ↔ νν with the Hubble expansion rate.
We conclude this section with a brief derivation of the scaling in Fig. 9. If we enforce that T X eq m χ m J , then Eq. (59) reduces to Eq. (60) implies that the sub-MeV scale for thermal DM is a natural consequence of the smallness of the observed neutrino masses. This numerical coincidence is surprising, since the MeV-scale has been motivated here in a completely independent manner, compared to the discussion in the beginning of this work. Hence, the framework and model described in the previous sections self-consistently motivate thermal DM below the MeV-scale.

VI. SIGNALS AND CONSTRAINTS
We now discuss signals and constraints for the model outlined in Secs. IV and V. These include cosmological and astrophysical considerations of the CMB, the small-and large-scale structure of matter, neutrino scattering in the early universe, DM self-interactions, and stellar cooling. We also briefly explore the possibility of observing more direct signals in terrestrial searches for light DM, sterile neutrinos, or majorons.
While many of the models are already tightly constrained by existing measurements, there remain viable regions of parameter space that will be decisively tested in the near future. This is illustrated explicitly in For every value of the dark matter mass, mχ, and dark matter-majoron mass ratio, mχ/mJ , the lepton number breaking scale, f , is fixed to reproduce the correct relic abundance, as in Fig. 8. Requiring that the hidden sector equilibrates with the neutrino bath at a given temperature sets a lower bound on the neutrino masses; in the blue shaded regions, this lower bound exceeds the upper limit on mν set by CMB measurements for ξX T X eq /mχ = 1, 3.
In the red shaded regions, dark matter free-streaming or acoustic oscillations in the hidden sector result in a cutoff in the matter power spectrum that is inconsistent with the smallest observed dark matter substructures. Since the smallest halo mass is subject to uncertainty, we show the resulting constraint for M cutoff = 10 9 M (solid red) and we fix ξ X T X eq = (1 − 3) × m χ , so that the DM sector equilibrates with the SM before χ is non-relativistic (well before freeze-out), analogous to the standard picture for thermal WIMPs. We also fix the lightest SM neutrino mass, m 1 , and the scale of U (1) L -breaking, f , as in Figs. 8 and 9 so that χ makes up the entire DM abundance at late times.

A. CMB
The general framework discussed in Sec. III will be decisively tested by observations of the CMB in various The CMB also constrains these models through indirect measurements of SM neutrino masses. Because the majoron is the pseudo-Goldstone of lepton number, its interactions with neutrinos are set by m ν /f , which in turn determines the equilibration temperature, T X eq , as described in Sec. V (see Fig. 9). For  [27,78]. However, it has been noted that uncertainties in the CMB lensing amplitude can significantly weaken these cosmological limits [78]. Hence, for simplicity, we show only the Planck TT constraint in Fig. 10, for various choices of the equilibration temperature.
B. Structure Formation

Dark Matter Free-Streaming and Acoustic Oscillations
The models considered throughout this work can lead to observable deviations in the observed matter power spectrum. Light DM that remains coupled to HS or SM radiation until late times can suppress power at small scales via two distinct mechanisms: free-streaming and acoustic oscillations. These processes wash out structure below a characteristic comoving length scale, λ cutoff , which sets a lower bound on the present day mass of the smallest gravitationally collapsed DM structures, where ρ DM = 1.26 × 10 −6 GeV cm −3 is the present cosmological DM energy density [33]. The cutoff scale is determined by solving Boltzmann equations describing the coupled DM-radiation system during the epoch of DM decoupling and free-streaming, which modifies the initial primordial matter power spectrum [79][80][81].
Here, we merely estimate the cutoff scales for the two effects following Refs. [21,82,83]. The scale that enters Eq. (61) is then given by the larger of the two lengths associated with free-streaming (λ FS ) and acoustic oscillations (λ AO ), We now discuss each of these in turn.
Once χ kinetically decouples from the radiation bath (either from HS majorons or SM neutrinos), it begins to freely diffuse across the universe, suppressing matter perturbations smaller than the free-streaming scale, λ FS . This length scale is defined as the comoving distance traversed by DM from the time of decoupling (assumed to occur during radiation domination) until matter-radiation equality, where a is the scale factor, v χ = p χ /E χ is the physical velocity of χ, t KD and t MRE are the cosmological times associated with DM kinetic decoupling and matter-radiation equality, respectively, and c FS is an O(1) number. There is some ambiguity in c FS due to different conventions and O(1) factors that appear in the Boltzmann equation treatment of free-streaming [81]. For example, in Ref. [83], c FS = 1/2, while Ref. [81] finds c FS = π/(2 √ 6) 0.64. In evaluating λ FS , we take c FS = 1/2. To simplify the evaluation of Eq. (63), let us assume that χ kinetically decouples while non-relativistic at a photon temperature of T KD O(MeV).
In this case, Eq. (63) can be simplified to where g eff * is defined as in Eq. (56), T 0 2.3 × 10 −4 eV is the present day photon temperature, T KD is the temperature of the photon bath at DM kinetic decoupling, and ξ X and g eff * are evaluated at T KD . Density fluctuations of the DM fluid that enter the horizon while DM is kinetically coupled to SM neutrinos and/or relativistic majorons oscillate with the radiation bath, similar to the baryonic acoustic oscillations in the baryon-photon plasma. The amplitude of these modes is damped due to their coupling to radiation. As a result, they do not undergo the usual logarithmic growth during radiation domination [81]. This results in suppressed power on scales smaller than the comoving horizon at decoupling, where a KD and H KD are the scale factor and Hubble parameter at DM kinetic decoupling. Once again assuming that DM kinetic decoupling occurs at temperatures T KD m χ , O(MeV), Eq. (65) is approximately where g eff * is evaluated at T KD , as in Eq. (64). In order to evaluate Eqs. (64) and (66), we need to determine the photon temperature at kinetic decoupling, T KD . The DM, χ, chemically decouples when χχ ↔ JJ freezes out (see Sec. V B), but remains in kinetic equilibrium with the SM bath directly through χν ↔ χν or indirectly through χJ ↔ χJ (+ J ↔ νν). Since χJ ↔ χJ is governed by the same couplings as χχ ↔ JJ, the fact that χχ ↔ JJ freezes out at T X ∼ m χ /10 implies that χJ ↔ χJ decouples at T X ∼ m J /10. For T X ∼ m J /10 and m J ∼ m χ , the rate for χJ ↔ χJ is enhanced over that of χν ↔ χν by approximately Hence, χν ↔ χν decouples well before χJ ↔ χJ, and we expect χJ ↔ χJ to dictate T KD . In the limit that m χ T X , m J , the differential rate for this scattering process is approximately where p J is the momentum of J in the center of mass frame and t is the usual Mandelstam variable.
We follow Refs. [21,83] in calculating the temperature at kinetic decoupling, T KD . We estimate T KD by equating the momentum relaxation rate for χJ ↔ χJ (denoted by γ) to the Hubble expansion rate, where γ(χJ ↔ χJ) is defined as Above, f J is the phase-space density of J, and dσ/dt is as given in Eq. (68). In the non-relativistic limit and taking f m χ , m J , this becomes Eqs. (69) and (71) allow us to estimate the kinetic decoupling temperature, T KD , through the relation where g eff * and ξ X are evaluated at T KD , and in the second equality we have fixed f to the thermally-favored value, as shown in Eq. (45) and Fig. 8.
The minimum halo mass, M cutoff , can be calculated using Eqs. The minimum halo mass constraint sets a lower limit on the DM mass of m χ (10 − 50) keV, for the thermal relic parameter space shown in Fig. 10. This is a stronger bound compared to the often-quoted limit on warm DM [4], which is usually assumed to have decoupled from the SM while relativistic at large temperatures. In the present model, the momentum of χ redshifts less between chemical decoupling and matter-radiation equality because χ remains coupled to the radiation bath of J and ν until late times.
As seen in Fig. 10, the bound becomes more severe for larger values of m χ /m J since χ decouples later (see Eq. (72)) as m J → 0. Furthermore, for m χ /m J few, the cutoff in the power spectrum (λ cutoff ) is controlled by free-streaming, while for larger values of m χ /m J , acoustic oscillations in the HS dominate. This can be understood by taking the ratio of Eqs. (64) and (66). For m χ ∼ O(10) keV, we find where in the second equality we have used Eq. (72). As a result, acoustic oscillations dominate over freestreaming in controlling the matter power spectrum cutoff for m χ /m J few. These limits will be improved in the near-future with, e.g., observations of the 21-cm hydrogen line in the cosmic dark ages [91][92][93]. For instance, an order of magnitude improvement in the sensitivity to λ cutoff would probe most of the remaining parameter space in Fig. 10.
Various studies have examined the effect of DM-neutrino scattering (χν → χν) on the matter power spectrum [82,[94][95][96][97][98][99][100]. We previously showed in Eq. (67) that this process decouples well before χJ ↔ χJ and therefore is not relevant for structure formation. However, for completeness we will compare the upper limits derived in the works listed above to the scattering rate for χν ↔ χν in our model. Majoron exchange dominates this process, since m J m S ; the low-energy cross section takes the parametric form where the T 4 temperature dependence arises from the CP-odd nature of the interaction between the majoron and the non-relativistic χ. For sufficiently large scattering rates, DM and neutrinos are tightly coupled in the early universe, altering the observed matter power spectrum, for instance, in large galaxy surveys. These effects constrain the size of the DM-neutrino opacity, Q ≡ σv(χν → χν) /m χ , where the temperature scaling of Q is parametrized as either constant, Q ∝ T 0 , or falling as the temperature squared, Q ∝ T 2 . In the case of constant scaling, the strongest bounds lead to the constraint Q 10 −33 cm 2 /GeV [98]. Since the predicted rate in Eq. (74) falls as T 4 , we conservatively compare the upper bound from Ref. [98] to the value predicted in our model at temperatures near matter-radiation equality, T ∼ eV, which gives the strongest possible constraint. We find that the predicted rate in our model is many orders of magnitude below this observational limit throughout the relevant parameter space shown in Fig. 10.

Dark Matter and Neutrino Self-Interactions
Non-standard neutrino interactions mediated by new forces (such as the majoron) can also alter the behavior of fluctuations in the photon and baryon fluids during the early universe. In the standard cosmology, neutrinos diffuse freely after decoupling from the photon plasma at temperatures of a few MeV until they become non-relativistic well after recombination. Such free-streaming radiation creates anisotropic shear stress, which, through gravity, suppresses the amplitude and shifts the phase of acoustic modes in the CMB that enter the horizon during this epoch [40,41,101]. However, if self-interactions (or interactions with another species) allow neutrinos to form a tightly coupled fluid before matter-radiation equality, the point at which they begin free-streaming is delayed. As a result, the strength of anisotropic stress is reduced compared to the SM expectation, and the power in subhorizon fluctuations is correspondingly increased and shifted in phase towards smaller angular scales.
Recent studies have investigated the effects of neutrino self-interactions (νν → νν) on the CMB, where the strength of the neutrino opacity is parametrized in terms of the dimensionful coefficient of a four-fermion operator, G eff [46,102,103]. These analyses have found that G eff 1/(50 MeV) 2 is consistent with data from Planck, the Sloan Digital Sky Survey, and local measurements of the Hubble parameter. In particular, for the models considered in Sec. IV, elastic neutrino scattering proceeds through the exchange of the light spin-0 mediators, J and S. In the limit that m ν eV m J,S , the relevant cross section is parametrically where the effective coupling is given by Since m J m S , elastic neutrino scattering is dominantly governed by majoron exchange, so that G eff ∼ This is orders of magnitude below the upper bound derived in from Refs. [46,[102][103][104]. We note that the ν − J coupling in the early universe also delays neutrino free-streaming until J becomes non-relativistic. The bound on delayed free-streaming in Refs. [102,103] can be stated in terms of a lower limit on the redshift at neutrino decoupling: z ν dec > 1.3 × 10 5 . For the masses m J keV, as considered in this work, ν decouples from J well before this epoch.
J and S exchange also gives rise to DM self-scattering (χχ → χχ). The self-scattering cross section per DM mass is bounded from observations of the dynamics and structures of galaxy clusters to be σ/m χ cm 2 /g, where the characteristic value of the relative DM velocity is v 2 ∼ 10 −5 [105][106][107]. We follow the discussion in Refs. [32,108,109] to calculate the viscosity cross section for the self-scattering of identical DM particles.
For m S ∼ f m χ and in the limit that v m J /m χ 1, DM self-scattering is dominated by majoron exchange, For m χ keV, this rate is maximized for m χ ∼ keV and f 30 MeV, where f has been fixed to the thermally-favored value in Fig. 8. This gives σ(χχ → χχ)/m χ 10 −6 cm 2 /g, which is orders of magnitude below the inferred upper bound.

C. Stellar Cooling
New particles coupled to the SM can lead to additional energy loss mechanisms in stellar systems, such as supernovae, red giants, and horizontal branch stars. One of the most powerful constraints on new light degrees of freedom comes from the observed cooling rate of SN1987A [110]. For m J 10 MeV, annihilations of SM neutrinos into a light majoron (νν → J) can lead to qualitative changes in the measured neutrino burst duration. Supernova bounds on majorons have been studied in detail in Refs. [111][112][113][114]. Here we estimate an upper bound on the J − ν coupling as follows. The energy loss rate per unit volume scales as [114], where Γ J ∼ m 2 ν m J /f 2 is the zero-temperature majoron decay rate (see Eq. (46)) and n ν is the neutrino number density for a given ν flavor. It is important to distinguish between electron and the heavy flavor neutrinos in the core. The former have a large chemical potential, µ νe 200 MeV, with n νe ∼ µ 3 νe , while the latter have a thermal population, such that n νµ,τ ∼ T 3 SN , where T SN ∼ 30 MeV is the core temperature. The larger electron neutrino density leads to a stronger constraint on model parameters (unless the electron-neutrino-like mass eigenstate is massless). A conservative bound on the anomalous cooling rate is obtained by requiring that the instantaneous majoron-luminosity, L J , does not exceed the total neutrino-luminosity of L ν = 3 × 10 52 erg/s [110]: where R c 10 km is the core radius and we have taken m ν = 0.1 eV to maximize the energy loss. Our estimate is in good agreement with the dedicated analyses performed in Refs. [111][112][113][114]. The lower bound on f in Eq. (79) is orders of magnitude below the thermally-favored values in Fig. 8. Other relevant processes involving neutrinos include neutrino annihilation into pairs of majorons, i.e., νν → JJ. However, compared to single majoron production, this rate is suppressed by an additional factor of (m ν /f ) 2 1. Finally, we note that right-handed neutrinos with a mass of M N ∼ 200 MeV can help restart stalled shock-fronts and facilitate supernovae explosions [115]. This is precisely in the cosmologically motivated region in Fig. 8 for branch stars in Ref. [116].

D. Direct Searches
Another avenue in exploring these models consists of direct searches for the light HS mediators (J, S, N ) and/or DM (χ). As discussed in detail in Ref. [71], limits on majoron-SM couplings are obtained from searches for flavor-violating processes, such as neutrinoless double beta decay, K → πJ, and µ → eJ, which constrain See Ref. [121] for a comprehensive review of such searches. While these limits are not sensitive to the natural parameter space of these models, they exclude non-trivial forms of the active-sterile mixing matrix, R (see Eq. (A8)), that lead to enhanced mixing in the neutrino sector.
Recent years have seen an increased focus on new experimental technologies to explore the sub-GeV DM frontier [122]. Of particular interest in this work are futuristic detectors proposed to detect elastic recoils of nucleons or electrons from DM as light as ∼ O(keV), corresponding to ∼ O(meV) energy depositions [28][29][30][31].
In this section, we investigate the potential sensitivity of these experiments to the classes of models discussed throughout this work.
The strength of χ − SM elastic scattering is controlled by the size of the S − h and J − h mixing angles, α and β, respectively (defined in Appendix A 2). For the cosmologically-favored parameter space in Fig. 8, Eq. (A27) suggests that for m χ ∼ (1 − 100) keV, β 10 −16 − 10 −12 is needed to avoid tachyonic states in the HS scalar spectrum. The prospects for such couplings to yield detectable rates is minuscule, and hence, J-mediated interactions with charged SM fermions are negligible within the context of direct detection experiments. In contrast, the S − h mixing angle, α, is not as constrained, so we focus on S-mediated interactions. The Yukawa coupling of the SM fermions to S is given by This can be matched onto a low-energy theory involving nucleons (n) and pions (π ± ) [123], angles larger than α ∼ 10 −6 . If α is set to its maximally allowed value and f is fixed to the thermal line in Fig. 8, we find that the DM-nucleon elastic scattering rate is well below the irreducible neutrino background, σ p 10 −50 cm 2 , while the electron scattering rate is many orders of magnitude below the sensitivities of futuristic proposed technologies [122].
We now consider variations upon these minimal models. We will first propose a modification in which the scalar mediator S is lighter than χ and the scale f and possesses additional couplings to the SM. As in Ref. [32], we assume that S also couples directly to SM QCD, through an interaction of the form where Λ is the cutoff of the effective theory. This interaction could be generated, for instance, from direct couplings to a vector-like generation of heavy quarks. As before, this can be mapped onto a theory involving nucleons and pions at low energies. Parametrically, this is of the form L ∼ y n Snn + y n m n S ∂ µ π † ∂ µ π .
As shown explicitly in Ref. [32], these couplings can lead to detectable rates in proposed low-threshold detectors for m χ ∼ m S ∼ 100 keV, without conflicting with cosmological, astrophysical, or terrestrial constraints.
In order to enlarge the viable parameter space, we propose a slight modification of the model in Ref. [32], which we now outline.
Compared to canonical WIMPs, physics at temperatures much greater than ∼ MeV is not directly important for models of sub-MeV thermal relics. In light of this, we will consider a low reheat temperature of the universe following inflation, T RH . The requirement of radiation domination during BBN implies that T RH few MeV [126,127]. We will take for concreteness. This is also motivated in models involving gravitinos and/or moduli [128][129][130]. We now ask: what are the maximum allowed values of the nucleon coupling, y n , such that the DM and visible sectors do not equilibrate before neutrino-photon decoupling? The decays and inverse-decays, J ↔ νν, are still assumed to equilibrate the two sectors below a few MeV. We find that processes involving protons, p, and pions, π, such as Sp ↔ γp and Sπ ↔ γπ do not equilibrate the two sectors before neutrino-photon decoupling provided that y n 10 −5 −10 −3 , where the lower (upper) part of the range corresponds to T RH ∼ 10 (5) MeV, respectively. By closing a loop of charged nucleons or pions, these couplings also generate an interaction with photons, which (modulo tuning) is naturally of size We demand that the processes S ↔ γγ also does not prematurely equilibrate the DM and visible sectors.
This leads to the additional upper bound y n 10 −4 (m S /100 keV) −1/2 . Hence, in order for equilibration to occur below the temperature of neutrino-photon decoupling, we will conservatively require that y n 10 −5 .
An exhaustive study of the constraints on DM-nucleon couplings in the context of MeV-scale particles has recently been presented in Ref. [32]. Here, we summarize the most relevant bounds. Considerations of cooling of horizontal branch stars constrain y n 10 −10 . However, this limit rapidly diminishes for m S 100 keV. For masses above ∼ 200 keV, the dominant constraints are from measurements of the meson decays, K → πS, leading to y n 10 −5 . For y n 10 −7 , S is produced but trapped in supernova, and bounds from anomalous cooling are evaded. Therefore, limits from meson decays and stellar/supernovae cooling restrict the nucleon coupling to be in the range 10 −7 y n 10 −5 (viable range) , for m S 100 keV. As argued above, for couplings of this size, DM-SM equilibration in the early universe is still driven by the neutrino-majoron coupling, as in our minimal scenario of Sec. V A. The DM-proton elastic scattering cross section is roughly For m S 100 keV, and taking y n ∼ 10 −6 , we have where we have fixed f to the thermally-favored value in Fig. 8. Proposed experiments, such as superfluid helium targets, are projected to be sensitive to cross sections as small as σ p DD ∼ 10 −42 cm 2 in this mass range [122].
In recent years, there has been growing interest in exploring new cosmological paradigms and modes of detection for particle dark matter in the keV − GeV mass range. For such light masses, dark matter that is of a thermal origin is strongly constrained from a plethora of cosmological and astrophysical considerations, including nucleosynthesis, the cosmic microwave background, structure formation, and stellar cooling. In particular, sub-MeV thermal relics that were in equilibrium with the Standard Model bath at temperatures below an MeV necessarily contribute to deviations in the expansion rate of the universe at the time of nucleosynthesis and/or recombination relative to the standard cosmology. As a result, models of sub-MeV thermal dark matter are usually thought to be either excluded or require involved model-building to evade these constraints.
We have focused on a class of models that naturally evade such claims. For instance, if a cold hidden sector equilibrates with the Standard Model after neutrino-photon decoupling, deviations in the expansion rate of the universe are strongly suppressed, alleviating the corresponding bounds from measurements of the effective number of neutrino species. Although this statement applies to dark matter that equilibrates either with neutrinos or photons, we have focused on interactions with the Standard Model neutrino sector. This is motivated, in part, by the fact that constraints derived from stellar cooling are much stronger for new light forces that couple directly to electromagnetism.
We studied concrete realizations of the above scenario where the dark sector masses and interactions, as well as the observed neutrino masses and mixing angles, are generated at a single scale corresponding to the spontaneous breaking of lepton number in the Standard Model. The pseudo-Goldstone boson associated with this breaking is the majoron, which is the mediator responsible for equilibrating the dark matter and Standard Model sectors in the early universe. These models independently motivate the sub-MeV scale; demanding that thermal dark matter freezes out with an adequate abundance implies that its mass is parametrically related to the Planck mass, the temperature at matter-radiation equality, and the measured neutrino masses by m DM ∼ (m Pl /T MRE ) 1/4 m ν ∼ MeV. Along with considerations of structure formation, this restricts the viable mass range to m DM ∼ 10 keV − MeV and the majoron-neutrino interaction strength to be at the 10 −10 to 10 −9 level.
Despite the suppressed size of such interactions, this class of models will be decisively tested in the near future. For instance, thermal relics that relativistically equilibrate with any Standard Model species after neutrino-photon decoupling lead to an irreducible deviation in the effective number of neutrino species above the projected sensitivity of future CMB-S4 experiments. Improved measurements of the small-and large-scale structure of the universe will also probe these models, potentially testing most of the remaining parameter space. Furthermore, it is possible to introduce a large coupling of the majoron to nucleons which preserves the viability of the cosmology provided that the reheat temperature of the universe is small (∼ 10 MeV). In this case, dark matter detection is possible at recently proposed low-threshold direct detection experiments aimed at exploring the sub-GeV dark matter frontier. which implies that where R is any complex matrix such that R R T = 1. Solving for m D gives Eq. (A4). As stated in Ref. [72], continuous forms of R (not including reflections) can be parametrized in terms of three complex angles. In the seesaw limit, V takes the form It is straightforward to check that Eqs. (32) and (A8) hold to leading order in d /d h . The off-diagonal entries in Eq. (A8) parametrize the active-sterile neutrino mixing.
Electroweak-and U (1) L -breaking leads to mixing amongst the neutrino states. We now switch to fourcomponent notation and denote the Majorana neutrino mass eigenstates as n i , i = 1, 2, . . . , 6, with mass m i , such that n 1,2,3 are SM-like, and n 4,5,6 are sterile-like. We parametrize the couplings of these states to the scalar sector as where h is the physical SM Higgs field. As shown in Ref. [70], the effective couplings are where following Ref. [70], we define In general, there may be other contributions to the masses of the sterile neutrinos. In this case, the mass parameters m 4,5,6 written in Eq. (A10) are interpreted as the piece given by the scale f , i.e., ∼ f × ∂m/∂f .
The interactions of the electroweak gauge bosons with the neutrinos are given by Za γ 5 n j + g (ij) where the couplings are defined as where the diagonal entries correspond to the masses of the unmixed fields: The mass matrix M 2 ϕ is diagonalized in the mass eigenstate basis, given by ϕ 1,2,3 . In the limit of small mixing the flavor eigenstates are related to ϕ 1,2,3 via where the small angles α, β, and γ are defined by .

(A25)
Large mixing in the scalar sector can lead to tachyonic masses. The most stringent constraint is obtained in the S − J sector (since J is the lightest state and mixing with h is suppressed by the large Higgs mass).
Requiring that the S − J eigenstates have positive masses bounds the mixing as This constraint can also be translated into a bound on β which limits the size of tree-level interactions of the majoron with charged SM fermions (see Appendix A 4).

Scale of Lepton-Number Breaking and Planckian Effects
In this Appendix, we briefly comment on two theoretical aspects of the majoron construction described above. We have introduced two new energy scales associated with spontaneous and explicit U (1) L -breaking, f and m J , respectively, with f m J . As we saw in Sec. V, considerations of DM-SM equilibration require f to be much smaller than the electroweak scale, i.e., f v 246 GeV.
The first issue associated with these new energy scales is the radiative stability of f . Quantum corrections will generically shift the mass term (and the resulting vev) of the U (1) L -breaking field, σ, to the UV cutoff of the theory, i.e., Λ v. As with the SM Higgs hierarchy problem, supersymmetry can be used to regulate the sensitivity to UV physics. If the HS (including σ) is sequestered from the supersymmetry-breaking sector, a naturally small f can be radiatively induced through interactions with the SM via the right-handed (s)neutrino [18]. However, a small supersymmetry-breaking scale in the HS also implies the presence of new light degrees of freedom (e.g. the superpartners of χ and ϕ) that can play an important role in cosmology.
A detailed investigation of this scenario is beyond the scope of this work.
The second puzzling feature of the majoron construction is the origin of the scale m J . If U (1) L was an exact symmetry (at least classically), the majoron would be massless, so m J > 0 requires an explicit breaking of U (1) L . While the hierarchy m J f is protected by the fact that J is a pNGB, it is interesting to ask why m J < f in the first place if they are completely unrelated. Global symmetries are expected to be absent in theories of quantum gravity. A simplified argument is that a scattering process with a global charge in the initial state can destroy the charge in an intermediate black hole state. The black hole cannot carry global charge, so it decays democratically via Hawking radiation [131,132]. This means that the low-energy effective field theories should have Planck-scale violations of global symmetries. This is a well-known problem in axion models with a global Peccei-Quinn U (1) [133][134][135]. Thus U (1) L -breaking effects should also appear in the low energy description [136,137].
If the Planck-scale effects are unsuppressed, then one expects mass terms ∼ m 2 Pl σ 2 to appear, which would remove any pNGB from the spectrum. Thus, if we want a light majoron, Planck effects should enter through marginal or irrelevant operators ∼ 1/m n Pl , n ≥ 0. The standard way to ensure this is to engineer U (1) L to be an accidental symmetry, i.e., one that is a consequence of gauge charge assignments as in Ref. [137]. This can be accomplished, e.g., using a gauged U (1) B−L with an additional scalar field ϕ, such that the leading U (1) L -breaking term is When ϕ gets a vev, these operators can be mapped onto the L-breaking terms in the σ potential in Eq. (40).
Note that for a given charge assignment with this minimal field content, only one of the potential terms is generated at dimension five. This means that it is a reasonable approximation to turn them on one at a time in this minimal framework.
Is there a mass-scale that is singled out by the Planck-suppressed operators? The answer depends on what the natural scale for spontaneous B − L breaking is. At the very least, one needs to account for existing bounds on the B − L gauge boson, a type of Z which has been extensively studied, see, e.g. Refs [138][139][140]. The LHC constrains m Z /g Z > 6 − 100 TeV through dilepton resonance searches and bounds on four-fermion contact interactions [141,142]. Letting ϕ = v B−L / √ 2 and m Z = q ϕ g Z v B−L , the above experimental bound implies where the range depends on the mass of the Z . In the minimal scenario with q ϕ = Q B−L [ϕ] = 4/3, the µ 2 σ σ 2 term is generated from a dimension-five operator The experimental bound then suggests a very rough lower limit on the majoron mass This bound can be much weaker if the mass is generated by an operator with a higher dimension or if its Wilson coefficient is not O(1). It can be larger if U (1) L is explicitly broken at a scale Λ < m Pl , e.g., the GUT scale. Thus, the natural size for the majoron mass (if B − L is broken near the weak-scale and the scale of explicit breaking is m Pl ) is near the keV-scale under the above assumptions. This was also noted in Ref. [136]. While this link is tenuous at best, it is reassuring that an internally consistent picture for the scales f and m J seems attainable.

Interactions with Charged Fermions
The mixing of the dark sector scalars with the SM Higgs gives rise to S and J coupling to SM fermions.
These interactions can be summarized by where we approximated ϕ 1 J, ϕ 2 S and ϕ 3 h. Note that the interactions of the 125 GeV Higgs-like state are suppressed relative to the SM expectation by an effective mixing cos θ eff 1 − α 2 /2 − β 2 /2 .
The strongest constraints on the scalar potential parameters come from rare meson and invisible Higgs decays. These were recently analyzed in Ref. [125] in the context of Higgs-portal coupled dark sectors. A detailed discussion of flavor physics constraints is presented in Ref. [143]. An invisibly-decaying light scalar, ϕ, that mixes with the SM Higgs contributes to the invisible decay modes B ± → K ± ϕ and K ± → π ± ϕ, whenever kinematically allowed. We are interested in J and S that are much lighter than m B − m K and m K − m π , so that the observed limits on these rare decay modes constrain the effective mixing sin 2 θ eff < 9 × 10 −6 (B ± → K ± + inv.) (A35) sin 2 θ eff < 3 × 10 −8 (K ± → π ± + inv.).
Measurements of the Higgs properties at the LHC also constrain the parameters of the scalar potential.
For example, since χ, N , S, and J are much lighter than h, there are new invisible decay modes. The invisible branching fraction of the Higgs is constrained to be less than 0.23 at 95% confidence level [144,145], leading to the bound where the terms correspond to h → χχ, h → SS, JJ, and h → n i n j , respectively.
Interactions of S and J with the charged SM fermions are also generated by loops of neutrinos via couplings in Eqs. (A9) and (A12) [70,71]. Their characteristic size (see Eq. (42)) corresponds to a tiny effective mixing of ∼ 10 −15 . Thus, even with the stringent constraints on the mixing angles, the tree-level interactions of S and J with charged SM fermions can be much larger than those induced by loops.