A Thermal Neutrino Portal to Sub-MeV Dark Matter

# A Thermal Neutrino Portal to Sub-MeV Dark Matter

## Abstract

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, . 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.

SLAC-PUB-17278

## I Introduction

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 , while microlensing searches for massive composite objects imply that  Allsman et al. (2001); Wyrzykowski et al. (2011); Tisserand et al. (2007). 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 number-changing 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 , 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 Viel et al. (2013); Baur et al. (2016). For larger masses, perturbative unitarity requires under the assumption of a standard thermal cosmological history Griest and Kamionkowski (1990). Thus, the thermal DM paradigm drastically restricts the possible mass range.

Although no theoretical inconsistencies arise for small masses, is often quoted as a robust lower bound on the mass of any thermal relic Ho and Scherrer (2013); Steigman (2013); Boehm et al. (2013); Nollett and Steigman (2014, 2015); Steigman and Nollett (2015); Green and Rajendran (2017); Serpico and Raffelt (2004). 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, . Sub-MeV thermal DM is relativistic at the time of nucleosynthesis and can modify . However, the successful predictions of standard Big Bang nucleosynthesis (BBN) and observations of the cosmic microwave background (CMB) constrain to lie near the SM expectation,  Mangano et al. (2002, 2005).

As originally pointed out in Refs. Chacko et al. (2005, 2004); Bartlett and Hall (1991) and recently studied in the context of light DM in Ref. Berlin and Blinov (2018), 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 Bertoni et al. (2015); Gonzalez Macias and Wudka (2015); Gonzï¿½lez-Macï¿½as et al. (2016); Batell et al. (2018a, b); Schmaltz and Weiner (2017), 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 Ade et al. (2016).

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 Hochberg et al. (2016a, b); Schutz and Zurek (2016); Knapen et al. (2017a). 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. Green and Rajendran (2017); Knapen et al. (2017b) 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 predict DM-neutrino couplings of size and independently motivate thermal DM near the MeV-scale through the coincidence of gravitational, matter-radiation equality, and neutrino mass scales, i.e., . We then turn to the cosmology and possible modes of detection in 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)1/2 ρ1/2radmPl , (1)

where is the Planck mass and we have assumed that the energy content of the universe is dominated by the radiation component, . The radiation energy density includes contributions from SM particles (, , ) and the dark sector. It is conveniently parametrized by the effective number of neutrino species, , such that

 ρrad≡ργ[1+(7/8) (ξSMν)4 Neff(T)] , (2)

where is the neutrino-to-photon temperature ratio in the standard cosmology (see Eq. (8) below). Thus, 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 as a late-time quantity (only to be evaluated at the time of recombination), in Eq. (2) parametrizes the expansion rate at temperatures below a few MeV. can be modified either by changing the actual number of degrees of freedom in the radiation bath or by altering . The notation for these and other relevant temperature scales is compiled in Table 1 for convenience.

Novel evolution of can modify the predictions of primordial nucleosynthesis and recombination. The outcomes of these cosmological epochs have been precisely measured and therefore constrain non-standard behavior of . Below, we summarize the effects of varying on aspects related to BBN and the CMB and then review how light dark sectors can run afoul of the resulting constraints.

### ii.1 Big Bang Nucleosynthesis

is constrained by observations of light nuclei abundances, as reviewed in, e.g., Ref. Patrignani et al. (2016). The abundances of helium-4, , and deuterium, , are measured with a precision of a few percent and therefore provide the most sensitive probes of the expansion rate during the epoch of nucleosynthesis. We now discuss these elements in turn.

In the early universe, neutrons and protons interconvert through weak processes such as . Once the temperature of the photon bath drops below the neutron-proton mass difference, , the neutron-proton ratio is approximately fixed, , where is the freeze-out temperature. Most of these neutrons are eventually converted into due to its large binding energy per mass (the remainder decays or ends up in deuterium or heavier nuclei). Hence, the mass fraction can be estimated by a simple counting argument, . Helium-4 is also produced in stars, but its primordial abundance can be observationally inferred, for instance, from measurements of recombination emission lines of ionized gas in low-metallicity dwarf galaxies Izotov et al. (2014).

Primordial nucleosynthesis is the dominant source of deuterium, since it is destroyed in stellar processes. Its abundance provides an additional handle on constraining the expansion rate at temperatures below an MeV. Deuterium also plays a crucial role in the production of through such reactions as followed by . Due to the small values of the deuterium binding energy () and baryon-to-photon ratio (), the production of light nuclei is delayed until , a phenomenon known as the “deuterium bottleneck.” However, unlike , once produced, deuterium is easily destroyed. Deuterium burning proceeds through the same reactions as mentioned above until . Its primordial abundance can be determined, e.g., through observations of absorption spectra of distant quasars Cooke et al. (2014).

Modifications to correspond to changes in the Hubble expansion rate. For , the expansion rate is enhanced, so that weak processes that convert freeze out earlier (at a larger temperature, ). As a result, the neutron-proton ratio, , is increased, leading to a larger primordial abundance with  Bernstein et al. (1989); Cyburt et al. (2016). Deviations in also modify the predicted abundance of deuterium. An increased cosmological expansion rate corresponds to a shorter time-scale for efficient deuterium burning during . Hence, for , the predicted deuterium abundance is increased.

If the baryon density is fixed by the observed nuclear abundances, recent detailed studies have determined  Cyburt et al. (2016) and  Foley et al. (2017) within during nucleosynthesis. The spread in the inferred value of is largely determined by the uncertainty in the primordial value of . This can be seen using  Aver et al. (2015) and the parametric relation  Bernstein et al. (1989). The best-fit central value of additionally depends on the inferred baryon-to-photon ratio, which is largely driven by the observed deuterium abundance.

### ii.2 Cosmic Microwave Background

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. Bashinsky and Seljak (2004); Hou et al. (2013). 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 Silk (1968). On the microscopic level, this corresponds to the stochastic process of photons Thomson scattering with free electrons. Hence, the diffusion distance, , can be written parametrically as , where is the number of scatters, is the photon mean free path, is the free electron number density, and is the Thomson cross section. The diffusion length scale is therefore . A larger (and correspondingly larger ) decreases the diffusion damping distance scale. As a result, photons travel a shorter average distance out of overdensities. However, observations of the CMB measure the angular scale of diffusion, , where is the angular distance to the surface of last scattering. is not independently determined, since it depends on the evolution of dark energy from recombination to present. The dependence on can be eliminated by considering the length scale of the sound horizon, , at the time of recombination. The position of the first acoustic peak in the CMB power spectrum is dictated by the corresponding angular scale, . Hence, the ratio of angular scales is independent of . The position of the first peak has been measured to a precision of  Ade et al. (2016). Thus, fixing to the observed value, the scaling argument above implies that larger (and hence ) leads to a larger , thereby suppressing power in the damping tail of the CMB. Note that the degree of damping at small angular scales is increased for larger , even though the underlying physical diffusion length is decreased. This behavior is seen explicitly in full Boltzmann simulations Bashinsky and Seljak (2004); Hou et al. (2013). The argument above also makes explicit the degeneracy between and ; since , the effect on from decreasing can be compensated by increasing . 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, at 68% confidence Ade et al. (2016). 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, , is in tension with the inferred value from Planck data at the level of  Riess et al. (2016). The inclusion of additional relativistic species at the time of recombination significantly alleviates the tension, favoring  Riess et al. (2016); Bernal et al. (2016); Brust et al. (2017); Oldengott et al. (2017). 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 Ade et al. (2016); Aghanim et al. (2016).

### ii.3 Standard Light Relics

Neutrinos decouple from the photon bath at a temperature of  Enqvist et al. (1992). A set of sub-MeV hidden sector (HS) particles (collectively denoted as ) that is equilibrated with the SM at temperatures below can lead to significant deviations in the observed value of . The lightest stable particle of this HS constitutes the DM of the universe. For simplicity, we assume that couples to the SM neutrinos and that all such particles have a common mass given by . We first consider the standard case where equilibrates with the SM neutrinos before the point of neutrino-photon decoupling, as has been investigated in Refs. Ho and Scherrer (2013); Steigman (2013); Boehm et al. (2013); Nollett and Steigman (2014, 2015); Steigman and Nollett (2015); Green and Rajendran (2017); Serpico and Raffelt (2004). 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, , in each bath () determines the entropy density, , and energy density, , where . For three generations of left-handed SM neutrinos,

 gν∗=(7/8)×3×2=21/4 . (3)

At temperatures below , the comoving entropy densities in the and photon bath are separately conserved. Using that and separately scale as ( is the scale factor), one finds

 gν∗+gX∗gγ∗ ξ3ν=constant, (4)

where

 ξi≡Ti/T (5)

is the temperature of species normalized to the photon temperature Feng et al. (2008). Treating electron-photon decoupling as instantaneous, we can approximate the number of relativistic degrees of freedom coupled to the photon bath as and . Equating Eq. (4) at temperatures above and below , and using that , one recovers the standard result

 ξν(T≲me)≃(411)1/3≃0.7 . (6)

When becomes non-relativistic, it heats up the SM neutrinos and negligibly contributes to the entropy density of the bath. Again using Eq. (4), but for and , we find that

 ξν(Tν≲mX)≃(411)1/3(1+gX∗gν∗)1/3 . (7)

For later convenience, we define as the value of assuming a standard cosmology () such that

 ξSMν≡{1 , T≳me(4/11)1/3 , T≲me . (8)

Using the above results, the defining expression for in Eq. (2) can be rewritten as

 (9)

In Eq. (9), we have assumed that decouples instantaneously once its temperature drops below its mass (), which is encapsulated by the Heaviside step function,  Feng et al. (2008); Berlin et al. (2016). Note that Eq. (9) reduces to when and . In the SM, neutrino decoupling is not instantaneous, and annihilations partially heat the neutrino bath, resulting in  Mangano et al. (2002, 2005). In Eq. (9), we have approximated . Substituting Eqs. (7) and (8) into Eq. (9), we find that

 Neff≃{3(1+gX∗/gν∗) , Tν≳mX3(1+gX∗/gν∗)4/3 , Tν≲mX , (10)

if equilibrates with the SM neutrinos at temperatures above . If , then Eq. (10) gives () at the time of nucleosynthesis (recombination) for . As discussed in Secs. II.1 and II.2, this is excluded from considerations of BBN and Planck measurements of the CMB by more than . Furthermore, realistic models of light thermal DM often require , leading to even larger deviations in . It is this basic insight that has driven many studies to claim that sub-MeV thermal DM is not cosmologically viable Ho and Scherrer (2013); Steigman (2013); Boehm et al. (2013); Nollett and Steigman (2014, 2015); Steigman and Nollett (2015); Green and Rajendran (2017); Serpico and Raffelt (2004).

## Iii Delayed Equilibration

### iii.1 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, , can lead to deviations in 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 , then such constraints are significantly relaxed Chacko et al. (2005, 2004); Bartlett and Hall (1991); Berlin and Blinov (2018).

Let us assume that a similar collection of sub-MeV particles () 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.3). As summarized in Table 1, we define and as the temperature of the photon bath at which enters and exits equilibrium with neutrinos, respectively, and 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 Hall et al. (2010), 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 equilibrates with the SM neutrinos after neutrino-photon and electron-photon decoupling, i.e., . 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 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 (), then the relevant Boltzmann equations are

 ˙ρeqX(TX)+3H(ρeq% X(TX)+PeqX(TX)) ≃−ΓdecmX(neqX(TX)−neqX(Tν)) ˙ρeqν(Tν)+4Hρeqν(Tν) ≃+ΓdecmX(neqX(TX)−neqX(Tν)) , (11)

where is the decay rate for , the superscript “eq’” denotes an equilibrium distribution, and 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 factors, but does not significantly change our results. Eq. (III.1) can be solved numerically for the evolution of as a function of the photon temperature, . The time variable can be traded for the photon temperature through the relation Venumadhav et al. (2016)

 ˙T=−3H(dρtotdT)−1(ρtot+Ptot) , (12)

where and similarly for the pressure density, . In Eq. (III.1), 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 . 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 couplings and keV-scale masses in the HS scalar potential of Sec. IV.3, 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 (). The initial HS-SM temperature ratio is fixed to , such that the HS is initially much colder than the SM neutrino and photon populations. Energy conservation then implies that equilibration cools (heats) the neutrino () bath at . If this occurs after neutrino-photon decoupling, this leaves the photon bath unaffected. Later, when the temperature drops below and the HS decouples, 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 , 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 , taking the standard assumption that equilibration occurs before neutrino-photon decoupling (dotted blue), as in Sec. II.3 and Refs. Ho and Scherrer (2013); Steigman (2013); Boehm et al. (2013); Nollett and Steigman (2014, 2015); Steigman and Nollett (2015); Green and Rajendran (2017); Serpico and Raffelt (2004). If equilibration occurs after neutrino-photon decoupling, deviations in are significantly reduced. We now derive analytic approximations for the asymptotic behavior of and , which are shown as the horizontal gray dashed lines in Fig. 2.

As we will soon see, is sensitive to the initial value of before equilibration or electron-photon decoupling, but, similar to DM production via freeze-in, it is insensitive to the particular value of as long as  Hall et al. (2010). We define 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 baths. Hence, as in Sec. II.3, for , we have

 ξν(TXeq≲T≲me)≃(411)1/3,ξX(TXeq≲T≲me)≃(411)1/3ξ0X . (13)

Along with Eq. (9), this implies that is given by

 Neff(T≳TXeq)≃3(1+gX∗gν∗ ξ04X). (14)

This is the standard result for an uncoupled population of dark radiation.

If the HS and neutrino baths equilibrate while and are relativistic, the sum of their comoving energy densities, , is approximately conserved. This can be seen from Eq. (III.1), which implies that , where . When and , we have and . Therefore,

 gν∗ξ4ν+gX∗ξ4X(gγ∗)4/3=constant, (15)

before and immediately after equilibration, where we have used . Equating this expression at temperatures above and below , we find

 ξνX(TXdec≲T≲TXeq)≃(411)1/3(gν∗+gX∗ξ04Xgν∗+gX∗)1/4 , (16)

where is the temperature ratio when is equilibrated with the SM neutrino bath. Comparing the above expression to the standard result of Eq. (8), we see that for , equilibration significantly lowers the temperature of the neutrino bath, i.e., . Eqs. (9) and (16) then imply that

 Neff(TXdec≲T≲TXeq)≃3(1+gX∗gν∗ ξ04X), (17)

during equilibration and before becomes non-relativistic. Note that Eq. (17) is identical to the expression of Eq. (14). This is consistent with the fact that and that is defined in terms of the total radiation energy density.

We use conservation of entropy when becomes non-relativistic and decouples, since this process occurs in equilibrium. Hence,

 gν∗ξ3ν+gX∗ξ3Xgγ∗=%constant, (18)

just before and after becomes non-relativistic. Equating this expression above and below and using Eqs. (16) and (9), we find

 ξν(T≲TXdec)≃(411)1/3(1+gX∗gν∗)1/12(1+gX∗gν∗ξ04X)1/4 , (19)

and

 Neff(T≲TXdec)≃3(1+gX∗gν∗)1/3(1+gX∗gν∗ξ04X) . (20)

Note that in the limit and taking , Eqs. (14), (17), and (20) reduce to

 Neff(T≳mX)≃3 (21)

and

 Neff(T≲mX)≃3(1+gX∗/gν∗)1/3≳3.18 , (22)

where in the inequality we have imposed for any light HS.

Compared to the standard result of Eq. (10), the deviation in away from its SM expectation is significantly reduced in Eq. (20) for . As mentioned previously, if , then equilibration drains the neutrino bath of energy, lowering its temperature compared to that of photons. Later, when 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 . However, as seen from Eq. (19), even for , 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 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.

Equations (14), (17), and (20) imply that constraints from nucleosynthesis and the CMB can be alleviated if and . In Fig. 3, we highlight regions of parameter space in the plane that are compatible with measurements of . If , then in general and its value encapsulates the sensitivity of our setup to physics in the ultraviolet. For instance, if was initially in thermal equilibrium with the SM but decoupled at before reentering equilibrium at , then . More generally, arises in theories of asymmetric reheating of the DM and SM sectors Adshead et al. (2016). Throughout this work, we take 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 . 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 Hall et al. (2010).

For , transitions from Eq. (17) to Eq. (20) near the decoupling temperature, . As a result, limits from nucleosynthesis depend on the ordering of and . Regions compatible with BBN are shown in Fig. 3 for both of the temperature orderings and . For , is static during BBN and is given only by the expression in Eqs. (14) and (17). However, for , evolves from the form given in Eqs. (14) and (17) to that of Eq. (20) during nucleosynthesis. Detailed studies of BBN, which demand within , often assume a single fixed value of throughout the entire formation of light nuclei Patrignani et al. (2016). However, as we have seen, this is not generally the case for a light HS that equilibrates and decouples from the SM during nucleosynthesis Hufnagel et al. (2018). In deriving a constraint, we demand that never deviates from the best-fit constant value by more than , i.e., for . We note that this is most likely overly conservative, since for and values of only slightly greater than , 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 abundance without affecting the production of . It would be interesting to consider the bounds from detailed investigations of BBN, while assuming time-variations of in this manner. We leave such considerations to future work Berlin et al. ().

Cold DM is necessarily non-relativistic at the time of recombination, i.e., . To remain consistent with Planck measurements of the CMB within , we demand that for , where we take the form for given in Eq. (20Ade et al. (2016). Note that this CMB bound on assumes standard nucleosynthesis, which is modified in the delayed equilibration scenario, as described above. A more realistic approach would be to fit both and to the CMB power spectrum. This can significantly expand the allowed parameter space due to the - degeneracy described in Sec. II.2. Also shown in Fig. 3 are regions of parameter space that alleviate the tension between Planck and local measurements of the Hubble parameter, . As a representative favored range, we take  Riess et al. (2016); Bernal et al. (2016); Brust et al. (2017); Oldengott et al. (2017). 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 (), or a light Majorana fermion and a real scalar (). The standard case of corresponds to the limit , which is in strong tension with measurements of both the CMB and primordial nuclei abundances for .

### iii.2 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 (), the temperature at matter-radiation equality (), and the Planck mass (), such that  Lee and Weinberg (1977). Equivalently, thermal DM that is lighter than a GeV often requires the presence of new light mediators Boehm and Fayet (2004). 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.1 is realized and modifications to during nucleosynthesis and recombination are reduced.

At temperatures much greater than or , we parametrize the rate for scattering and decays/inverse-decays as

 Γscatt ∼α2eqT  (scattering) Γdec ∼αeqm2φ/T  (decays) , (23)

where is the effective coupling governing equilibration and the factor of in the decay rate is a time-dilation factor. Comparing either process to the Hubble parameter, , demonstrates that the rate for equilibration overcomes the expansion rate at temperatures below

 TXeq ∼αeqmPl  (scattering) TXeq ∼(αeqm2φmPl)1/3  (decays) (24)

for scattering and decays, respectively, where denotes the temperature at which the DM and SM sectors equilibrate. If we parametrize the rate for DM annihilation during freeze-out as , then acquires an abundance in agreement with the observed DM energy density for

 mχ∼αFO(TMREmPl)1/2 , (25)

where is the effective coupling governing freeze-out. Using this relation in Eq. (III.2) allows us to write in terms of ,

 mχ∼{(αFO/αeq) (T% MRE TXeq)1/2  (scattering)(αFO/αeq)1/3 (mχ/mφ)2/3 (TXeq/mPl)1/6 TXeq  (% decays) . (26)

Equation (26) implies that and equilibrate with the SM after neutrino-photon decoupling () if

 mχ≲keV×{(αFO/α%eq)  (scattering)(αFO/105αeq) (mχ/mφ)2/3  (decays) . (27)

Bounds on warm DM typically exclude  Viel et al. (2013); Baur et al. (2016). Therefore, along with Eq. (27) motivates . 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 Pospelov et al. (2008), those involving freeze-out through resonant annihilations Griest and Seckel (1991), or strongly interacting hidden sectors Hochberg et al. (2014). Once and/or become non-relativistic, and are either suppressed by Boltzmann or 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 . This is typical of the weak processes that maintain - equilibrium where . 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 . 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.

### iii.3 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 Hall et al. (2010). 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.1, 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 were derived in Sec. III.1. If the HS and neutrino baths equilibrate while DM is semi- or non-relativistic, is no longer conserved. Instead, the system of Boltzmann equations in Eq. (III.1) must be solved numerically. Such models are shown as the dashed and dotted blue contours in Fig. 5.

We show the temperature evolution of for these generalized scenarios in Fig. 6, analogous to the right panel of Fig. 2. The various contours correspond to the examples shown in Fig. 5. For each of these lines in Fig. 6, HS-SM equilibration occurs after neutrino-photon decoupling. The solid blue contour corresponds to HS-SM equilibration while the DM is relativistic, as studied in Sec. III.1. For the dashed and dotted blue contours, equilibration occurs instead when the DM is semi- or non-relativistic, as illustrated in Fig. 5. In Sec. III.1, we noted that the increase in 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 system, i.e.,

 dSν+X=dQ(1TX−1Tν)>0, (28)

where 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 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.1.

## 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.1, we discuss the basic framework that is needed to generate the appropriate parameters in the neutrino sector. In Sec. IV.2, 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.3. A more detailed discussion concerning the explicit forms for the masses and interactions of the HS particles is given in Appendix A.

### iv.1 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,  Weinberg (1979). Below the scale of electroweak symmetry breaking, this operator generates neutrino masses parametrically of the form , where is the SM Higgs vacuum expectation value (vev) and 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 Schechter and Valle (1980); Gell-Mann et al. (1979); Mohapatra and Senjanovic (1980); Yanagida (1979); Minkowski (1977). If neutrinos are Majorana, then breaks lepton number, . The global symmetry can be broken explicitly, as in minimal seesaw models with an explicit Majorana mass for the right-handed neutrinos, or spontaneously when a -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  Chikashige et al. (1981); Gelmini and Roncadelli (1981); Schechter and Valle (1982). 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. Pilaftsis (1994) and Garcia-Cely and Heeck (2017). We introduce a complex scalar, , of lepton number ,

 σ=1√2 (f+S+iJ) , (29)

where we have assumed that acquires a non-zero vev, . and are the real and imaginary excitations of , where (often dubbed the majoron) is the Goldstone boson of spontaneous -breaking. In the presence of suppressed terms that softly break lepton number, is a pseudo-Goldstone and acquires a small mass. Soft -breaking terms can arise in the scalar potential, which is examined in Sec. IV.3 and Appendix A.2. While we naturally expect , we will not specify the exact form of -breaking and treat the majoron mass, , 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, , with lepton number . The most general renormalizable and -symmetric Lagrangian coupling and to the SM lepton sector is then given by

 −L⊃yνLNH+12yNσN2+h% .c. , (30)

where two-component spinor and flavor indices are implied. Above, and are the SM lepton and Higgs doublets, respectively. Below the scale of electroweak and -breaking, the interactions in Eq. (30) give rise to the neutrino mass matrix in the basis,

 MνN=(0mDmTDMN) , (31)

where and are mass matrices. Diagonalizing gives rise to the neutrino mass basis, (), with masses . We define the unitary matrix that diagonalizes the full active-sterile neutrino mass matrix by

 V†MνNV∗=diag(m1,…,m6) , (32)

where relates the gauge and mass eigenstates.1

In the seesaw limit (), and are SM-like and sterile-like neutrino species, respectively, with masses schematically of the form and . The off-diagonal entries of correspond to active-sterile mixing and are suppressed by . This is made explicit by the Casas-Ibarra parametrization as discussed in Appendix A.1 Casas and Ibarra (2001). The interactions of the neutrino mass eigenstates () with the scalar degrees of freedom take the parametric form

 −L∼(mimj)1/2(S+iJf+hv)ninj+h.c. , (33)

where 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 sterile neutrinos. In this case, the mass parameters written in these interactions are implicitly assumed to be the piece given by the scale , i.e., . However, it is important to keep in mind that is still possible in extended models. We will return to this point later in Sec. V.1.

Mass-mixing in the neutrino sector also induces interactions of the sterile states with electroweak currents and generates couplings of and 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.3 and VI.3 as well as in Appendix A.

### iv.2 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, , 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 and charged under an additional . The 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

 −L⊃12λχσχ2+h.c. (34)

The phase of can be chosen such that the Yukawa coupling, , is purely real. Below the scale of -breaking, acquires a mass,

 mχ=λχf√2 . (35)

In four-component notation, the interactions of the Majorana fermion, , with and are given by

 L⊃ λχS S ¯¯¯¯χχ+λχJ J ¯¯¯¯χiγ5χ , (36)

where the couplings are defined as

 λχS ≡−λχ2√2=−mχ2f, λχJ ≡λχ2√2=mχ2f . (37)

### iv.3 Scalar Sector

The -preserving renormalizable scalar potential is given by

 VL(H,σ)=−μ2H|H|2+λH|H|4−μ2σ|σ|2+λσ|σ|4+λσH|σ|2|H|2 . (38)

This potential does not generate a mass for the majoron, . However, soft -breaking terms such as

 V⧸L=−(μ′2σσ2+aσσ|H|2+h.c.) (39)

can give rise to a radiatively-stable mass for . The full potential is then given by

 V=VL+V⧸L . (40)

We fix the phase of such that its vev, , is real, leaving a single physical phase in the couplings and . This phase leads to CP-violating mixing of with and . 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 . 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.3 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 () are nearly aligned with the gauge basis (), with masses

 m21 ≃m2J≃4Reμ′2σ+Reaσv2/√2f m22 ≃m2S≃2λσf2 m23 ≃m2h≃2λHv2 . (41)

This assumption will be relaxed in Sec. VI.4 when we consider possible signals in futuristic low-threshold direct detection experiments. For , the mass of the CP-even scalar, , is near t