Cosmology and prospects for sub-MeV dark matter in electron recoil experiments

Dark matter is poorly constrained by direct detection experiments at masses below \SI{1}{\MeV}. This is an important target for the next generation of experiments, and several methods have been proposed to probe this mass range. One class of such experiments will search for dark matter--electron recoils. However, simplified models with new light degrees of freedom coupled to electrons face significant pressure from cosmology, and the extent of these restrictions more generally is poorly understood. Here, we perform a systematic study of cosmological constraints on models with a heavy mediator in the context of an effective field theory. We include constraints from (i) disruption of primordial nucleosynthesis, (ii) overproduction of dark matter, and (iii) the effective number of neutrino species at recombination. We demonstrate the implications of our results for proposed electron recoil experiments, and highlight scenarios which may be amenable to direct detection.


I. INTRODUCTION
The identity of dark matter (DM) remains one of the most significant problems in cosmology and particle physics. Over the past few decades, experimental efforts to detect and characterize DM have been guided by the assumption that the dark species is a weakly interacting massive particle (WIMP). However, despite substantial improvements to experimental sensitivity, neither astrophysical nor terrestrially-produced DM has been definitively detected. Increasingly strong constraints have placed the WIMP paradigm under pressure [1][2][3][4][5], spurring the development of new models across the mass spectrum.
In the meantime, direct searches for DM have largely targeted the weak scale. Most extant direct detection experiments are designed to detect the scattering of DM with atomic nuclei, and due to kinematic limits, they have poor sensitivity to a DM particle with mass below 10 GeV [6][7][8]. Analyses of the phase space distribution of DM in dwarf spheroidal galaxies bound the mass of fermionic DM to m DM > ∼ 1 keV regardless of the production mechanism [9], and the Lyman-α forest imposes a comparable constraint on thermal relic DM of any kind [10]. But beyond these bounds, DM models with mass between 1 keV and 10 GeV are poorly constrained. Several well-motivated scenarios [e.g. asymmetric DM,11] naturally feature masses between 1 keV and 10 GeV, making this range an appealing target for future direct detection experiments [12].
This has driven much interest in novel detection methods suited to light DM particles, and several such experiments have been proposed in the last few years [13][14][15][16][17][18][19][20][21][22][23] (see sections IV-V of [24] for a review). These experiments are designed to be sensitive to the very small recoil energies characteristic of the scattering of light particles, and as such, many are designed to search for the scatter- * blehmann@ucsc.edu † profumo@ucsc.edu ing of DM with electrons instead of nuclei, a strategy first detailed in [13]. Several experiments now constrain DMelectron scattering at masses as low as ∼ 1 MeV [25][26][27][28][29][30].
The more recent proposal of [15], based on electrons in aluminum superconductors, is sensitive to deposited energies of order 1 meV, allowing for the detection of particles as light as 1 keV. However, although the most generic astrophysical constraints do not restrict DM at masses between 1 keV and 10 GeV, it is well known that particular models can be constrained by cosmological observables, especially for masses below 1 MeV [31]. In particular, light DM interacting with electrons risks running afoul of the following restrictions: • The DM must not significantly alter successful predictions of the ratios of light elemental abundances produced in big bang nucleosynthesis (BBN) [32][33][34]; • To accord with measurements of the effective number of neutrino species (N eff ), the thermal history of the DM species must not significantly alter the temperature ratio of photons and neutrinos at recombination [35]; • While a single species of DM particle may not account for the entirety of the present-day DM density, no species may be produced with an abundance exceeding that threshold.
In each case, such cosmological constraints bound the couplings between new species and Standard Model (SM) particles, which also determine the event rates in direct detection experiments. Thus, in a given model, the cosmological effects of light DM can be related to the direct detection cross section. Given an experimental proposal and a DM model, one can then determine the extent of the parameter space accessible to the experiment and consistent with cosmology. Such an approach has been applied to electron recoil experiments by [36] for a class of simplified models, and more recently in a variety of model-dependent instances [37][38][39][40][41][42][43][44][45].
In this work, we show that cosmological constraints on a new light (sub-MeV) species interacting with electrons can be greatly generalized with a small number of assumptions. Assuming a heavy mediator between DM and the SM, we study the cosmological implications of a light DM species in an effective field theory (EFT), and use the same EFT to evaluate direct detection prospects. We thus obtain model-independent cosmological limits on the scattering cross section of DM with electrons in an actual experiment. The model-independent methodology is similar in spirit to [46][47][48], but applied to directly connect cosmological constraints and detection prospects in the sub-MeV regime.
This paper is organized as follows. In section II, we describe our EFT framework for modeling light DM coupled to electrons. In section III, we derive model-independent cosmological constraints on the DM species. In section IV, we evaluate the DM-SM scattering cross section in our EFT, and compare cosmological bounds with prospects in a fiducial experiment. Finally, we discuss implications for direct detection experiments in section V. A complete set of constraints and tables of cross sections are placed after the end of the text.
Throughout this work, we denote a scalar DM field by φ and a fermionic DM field by ψ. When speaking about the DM species generally, without specifying its spin, we will denote it with χ.

II. EFFECTIVE INTERACTIONS OF SUB-MEV DARK MATTER
In this section, we build a theoretical framework to study the effective interactions of sub-MeV DM of spin 0 or 1 2 . We study DM candidates that are singlets under the SM gauge groups, and we consider both scalar and fermionic DM. We first specify the working assumptions of our EFT framework, and we thereafter develop the scalar and fermion cases separately.

A. The EFT framework
We assume that DM is dominated by a single particle species with a mass below 1 MeV. The MeV scale is cosmologically significant as the scale of big-bang nucleosynthesis (BBN). The DM annihilation and scattering processes that we consider in this work always involve energy exchanges well below this scale, whether they take place in the early universe or in a laboratory today. Thus, this situation lends itself well to an effective low-energy description.
At energies well below the MeV scale, the only dynamical SM degrees of freedoms are electrons and positrons (e ± ), neutrinos (ν), and photons (γ). We assume further that there is no additional light degree of freedom besides the DM particle: all remaining new physics is presumed to lie well above the MeV scale, including any media-tors between DM and SM particles. Physics at sub-MeV scales is thus well described by an EFT in which only e ± , ν, γ, and the DM χ are dynamical degrees of freedom. This is the theoretical framework we employ for our analysis.
Before presenting the EFT in more detail, it is instructive to take a step back and discuss the conceptual starting point of our work: a renormalizable theory with DM as well as mediator fields in the spectrum. The EFT language powerfully encodes the many UVcomplete realizations which give the same low energy physics. We make three additional assumptions about the UV-complete theory, described below and graphically summarized in fig. 1: 1. The DM is stablilized by a Z 2 symmetry and is thus absolutely stable.
2. The couplings between mediators and SM fields respect electroweak gauge invariance, in the sense that the χ-e L coupling is equal to the χ-ν coupling. We make this assumption to clarify the impact of the DM species on N eff , as discussed in the next section. It does not influence the other constraints.
3. DM couples to the visible sector via mediator fields ζ i , with masses satisfying T BBN m ζi . When writing our EFT Lagrangian, we additionally assume that m ζi m weak 100 GeV. However, as we explain below, our results do not depend on this assumption-it simply guides the form of certain EFT operators we consider.
The last of these assumptions is made only to clarify how we should write the low-energy Lagrangian to accommodate lower mediator masses. Ultimately, our EFT will contain a mass scale Λ EFT which is related to the mediator masses, and we will deal only with Λ EFT rather than the m ζi . Each operator will appear suppressed by some power of Λ EFT , and with an additional coupling (Wilson coefficient) g. We ensure that we remain in the regime of validity of the EFT by enforcing Λ EFT T BBN , so it is convenient to assume that g ∼ O(1) and take Λ EFT to be the free parameter in our analysis. Small deviations of g from unity can then be absorbed by rescaling Λ EFT . But if g is not O(1) in a typical UV completion, and Λ EFT is not many orders of magnitude larger than T BBN , we have reason for caution: rescaling Λ EFT to absorb a very small g could violate the requirement that Λ EFT T BBN . This would make it difficult to interpret results computed with g = 1.
Thus, at smaller mediator masses, it is important to correctly estimate the normalization of the operators, and to separate g from any non-O(1) coupling typical of UV completions. When writing the operators in our EFT, we assume that m ζi m weak so that weak scale particles (i.e., the weak gauge bosons W and Z, the Higgs boson h, and the top quark t) can always be integrated out before the mediators. It is then possible to define an intermediate EFT with weak scale particles integrated out and mediators in the spectrum. This intermediate theory then guides our expectations for the size of the coupling in the effective theory after integrating out the mediators, a fact we will use in section II B.
Finally, note that we ignore any renormalizable couplings between the DM and SM fields, assuming that all interactions are encoded in the EFT. Notice that no such operators exist in the fermionic case under our assumptions, since we take the DM to be a SM singlet, and the Z 2 symmetry forbids the lepton portal operator φLH. In the scalar case, on the other hand, this is something we impose. However, as we will discuss shortly, this assumption has no consequences for the results of our analysis.
At energies at or below the scale of BBN, the effective Lagrangian schematically reads Here Λ EFT is the mass scale associated with the EFT, which reflects the scale of the heavy degrees of freedom in the theory; L SM is the SM Lagrangian with only the e ± , ν, and γ fields; and L DM is the DM free theory contribution. The form of L DM depends on whether the DM is a scalar φ or a fermion ψ. If the DM is a scalar, then and if the DM is a fermion, then The remaining (infinite) sum over the higher-dimensional operators in eq. (1) accounts for the effective interactions between DM and SM fields. In our analysis, we will retain terms up to dimension 6.
In the following subsection, we parametrize the interactions between DM and electrons. All operators consistent with a Z 2 symmetry have the schematic form where the function B I (χ) contains an even number of DM fields, and I denotes a set of Lorentz indices. We will eventually truncate all operators beyond dimension 6, so for our purposes, B I (χ) always contains two DM fields. This DM bilinear is multiplied by an electron bilinear, for which the independent Dirac structures can be fully enumerated: If the electron bilinear is not a Lorentz scalar, the contraction of its free Lorentz indices with the ones of the DM bilinear ensures that the full operator in eq. (4) is a Lorentz invariant. We now discuss the allowed operators for scalar and fermion DM.

B. EFT for scalar DM
To describe our EFT for scalar DM, we must enumerate all operators of the form up to some mass dimension. Note that φ carries no Lorentz indices or spinor indices. Thus, if the index set I carried by the electron bilinear is non-empty, the only option is to insert derivatives in the scalar bilinear so that all indices are contracted. A classification of all possible cases is provided in table I. Of the four resulting operators, two are dimension-5, while the other two include a derivative and are dimension-6. We use the notation Note that we omit the operator ∂ µ φ † φ + φ † ∂ µ φ ēγ µ e, since it vanishes under integration by parts and application of the equation of motion: Similarly, the operator ∂ µ φ † φ + φ † ∂ µ φ ēγ µ γ 5 e is redundant: integrating by parts again, we obtain The resulting integrand is proportional to O (φ) P , one of the other operators in our basis. Moreover, this contribution is dimension-6 while O (φ) P is dimension-5, so it is suppressed in the Lagrangian with an additional factor of Λ −1 EFT . In some cases, renormalizable operators are allowed, and might appear in addition to the effective operators discussed above. For instance, in the context of a Higgs portal model [see e.g. 49] the operator φ † φ H † H is allowed without affecting DM stability. After electroweak symmetry breaking (EWSB), this operator produces a cubic coupling φ † φ v h. Integrating out the SM Higgs boson generates an effective operator proportional to O (φ) S . Thus, adding renormalizable couplings does not introduce any new physical effects in our analysis. The only effect is to add a correction to the Wilson coefficient of a single operator, with a size typically smaller than the values we consider in our analysis.
Each of the operators in table I coupling φ to the SM is normalized with a factor of the electron Yukawa coupling y e . This factor generically appears in UV completions of the effective theory: if we are to couple the electron to φ without violating the Z 2 symmetry, then we must employ a scalar mediator ζ. In turn, coupling ζ to the lepton doublet L without breaking gauge invariance involves interaction terms of the form Thus, after EWSB, ζ mixes with the Higgs boson h. To construct an EFT from the Lagrangian in the broken phase, we must integrate out the mass eigenstates corresponding to (ζ, h), which will always produce a factor of y e in addition to the inverse of the mediator mass scale. This is important for power counting considerations. For instance, given that L EFT ⊃ O ∝ φ † φēe, one might naively account for the apparent parametric dependence on m e by writing e.g.
This would suggest that O ∝ Λ −2 EFT , but as we have argued, a typical UV completion would instead give rise to an effective operator O y e Λ −1 EFT φ † φēe. With the effective operators in the scalar case now enumerated, we can consider annihilation and scattering processes for each one. Matrix elements for 2 → 2 annihilation and scattering are given in table III. The corresponding cross sections are given in tables IV and V.

Symbol Operator
Real case No   TABLE I. Operators coupling the electron to a dark scalar φ. The third column indicates whether or not the operator survives when φ is taken to be a real scalar.

C. EFT for fermion DM
If the DM is a fermion ψ, the structure of the EFT is similar to the scalar case. We again have a set of operators which are products of an electron bilinear and a ψ bilinear. Using generalized Fierz identities, it can be shown that operators of the form (ψO 1 e)(ēO 2 ψ) are redundant, in that they can be written as linear combinations of operators of the form (ψO 1 ψ)(ēO 2 e) [50]. Thus, we can construct a complete basis of effective operators by enumerating the possible insertions O 1 and O 2 . All of the electron bilinears from the scalar case appear here as well, and most of the possible ψ bilinears are obtained from these by making the replacement e → ψ.
In addition to these bilinears, we can form a spin-2 current at dimension 6, e.g. of the formψσ µν ψ. Since σ µν is antisymmetric, the other bilinear must not be symmetric in its Lorentz indices, so it must contain another insertion of σ µν . Thus, such an operator has the general form W µναβψ σ µν ψēσ αβ e. At dimension 6, the indices of W µναβ can come only from two factors of the metric or one factor of the Levi-Civita symbol ε. In the former case, again due to antisymmetry of σ µν , the only nontrivial contraction is If W is instead formed from the Levi-Civita symbol, then the operator has the form ε ρ1ρ2ρ3ρ4ψ σ µν ψēσ αβ e, where (ρ 1 , ρ 2 , ρ 3 , ρ 4 ) is a permutation of (µ, ν, α, β). Up to an overall sign, the indices ρ i can be rearranged into the latter order, so all such operators are proportional tō But ε µναβ σ αβ = −2iσ µν γ 5 , so if we simply add iσ µν γ 5 to our list of insertions, we can assume that W µναβ is a product of metric tensors. (We retain the factor of i to preserve Hermiticity.) Further, the argument above demonstrates that it is sufficient to place this insertion in only one of the two bilinears: the operator formed by inserting iσ µν γ 5 in both bilinears is redundant. We choose to place this insertion in the electron bilinear. The complete list of operators for fermionic DM is shown in table II. Matrix elements for 2 → 2 annihilation and scattering are given in table VI. The corresponding cross sections are given in tables VII and VIII.

III. COSMOLOGICAL CONSTRAINTS
Cosmological constraints on DM are typically modeldependent. However, the broad class of models which we consider admits only a very restricted set of thermal histories for the DM species, which allows us to derive general cosmological constraints in the context of our EFT.
We divide the thermal histories into two cases: either the DM is in thermal equilibrium with the SM at high temperatures, and freezes out below some temperature; or it never attains thermal equilibrium, and the abundance is instead set non-thermally. It is possible that the dark species only enters equilibrium at late times, but this scenario mirrors the thermal freeze-out case in almost every respect.
In the freeze-out scenario, two constraints are particularly robust: first, if the DM remains in thermal equilibrium during the epoch of big bang nucleosynthesis (BBN), its effect on the Hubble parameter is generally sufficient to perturb light elemental abundances [34]. Second, if at some temperature the DM is in thermal equilibrium with electrons and not with neutrinos, or vice versa, then entropy can be transferred from the DM to neutrinos alone or to electrons and photons alone. This changes the temperature ratio of the two thermal baths, which modifies the effective number of neutrino species, N eff , as determined from the cosmic microwave background (CMB) [35].
Finally, in the case of out-of-equilibrium (non-thermal) production, the DM never attains thermal equilibrium, and so may evade these two constraints. However, if the coupling to electrons is too large, DM will be overproduced even under the most generous assumptions.
Note that new light species are also subject to constraints from energy loss in stars and supernovae [51]. However, these constraints rely on complicated microphysical inputs that must be computed in detail for each model. Moreover, supernova temperatures lie up to an order of magnitude above the scale of BBN, requiring our effective theory to be valid at higher energies. Thus, we do not evaluate these constraints explicitly, but simplistic estimates suggest that they are at best comparable in strength to our cosmological constraints over the mass range of interest.
We now examine each of our constraints in more detail.
A. Freeze-out and primordial nucleosynthesis Light element abundances today are a sensitive probe of cosmology at scales near 1 MeV. If an additional light species is assumed to be in thermal equilibrium at these scales, the standard predictions of big bang nucleosynthesis (BBN) are modified, with observable consequences. Since thermal equilibrium in turn depends on DM interactions, light element abundances translate to stringent constraints on the interaction rates.
In a broad class of models, the DM species is in thermal equilibrium with the SM bath at high temperatures, and eventually drops out of equilibrium below some freezeout temperature, T FO . In our framework, freeze-out is a generic requirement of any scenario in which DM is in thermal equilibrium with electrons at temperatures T < ∼ 1 MeV, since the EFT is valid in this regime. If the DM species freezes out during or after BBN, and the DM species is in equilibrium at higher temperatures, then the predictions of light element abundances are generally perturbed to a degree incompatible with their measured values [32][33][34]52]. The ratios of these abundances are set by the temperatures at which interconversion processes freeze out, which depend in turn on the Hubble parameter H. Since H is sensitive to the energy density, adding a new species that stays in equilibrium and remains relativistic for much of the epoch of BBN has a significant impact on the produced light element abundances. Note that in a small range of our parameter space, equilibrium during BBN is consistent with observables if the dark species enters equilibrium at a specific time during BBN [40]. This is a very narrow exception to our framework, so we neglect it for the remainder of this work.
The temperature at which freeze-out occurs is fixed by the DM mass and the couplings. The prospect of experimental detection by any particular apparatus places a lower bound on the scattering cross section χe − → χe − . However, for a given interaction, the scattering cross section is directly related to the annihilation cross section χχ → e + e − which regulates the thermodynamics of the DM species in the early universe. A lower bound on the scattering cross section thus corresponds to a lower bound on the annihilation cross section, which translates to an upper bound on the freeze-out temperature. Thus, if DM is in thermal equilibrium with the SM at high temperatures, BBN constrains the scattering cross section with electrons. For our purposes, we will simply consider a model to be ruled out by light element abundances if it predicts that DM is in equilibrium at T = 1 MeV. A more detailed treatment would involve numerical integration of a system of Boltzmann equations to compute the yield of each light element, and would include the effects of energy injection due to ongoing DM annihilation [see e.g. 34]. However, ultimately, such a calculation would only slightly shift the minimum allowed freeze-out temperature, corresponding to a very small correction to our constraints.
The freeze-out temperature and relic density for a given model are also found by solving the Boltzmann equation, but in a relatively simple incarnation. In our framework, we have only a single DM species χ which interacts with electrons exclusively through 2 → 2 processes. For this case, using Maxwell-Boltzmann statistics, the Boltzmann equation takes the form where x ≡ m χ /T parametrizes cosmic time; σ is the cross section forχχ → e + e − ; Y ≡ n/s is the abundance of χ, where n is the number density and s the entropy density of χ; and Y eq and n eq are the equilibrium abundance and number density of χ, respectively. We identify Γ A ≡ n eq σ |v| as the annihilation rate of χ when in equilibrium. The thermally-averaged cross section can be obtained as [53] σ It is clear from eq. (15) that the abundance will stabilize once Γ A /H < ∼ 1. This condition gives an estimate of the temperature T FO at which χ departs from equilibrium, and thus allows us to test whether a set of parameter values is consistent with BBN observables.
In general, when studying the decoupling of χ, it is important to consider the coupling to neutrinos as well as electrons. If χ has a non-negligible coupling to neutrinos, it is conceivable that the DM could be kept in equilibrium at later times via thermal contact with the neutrino bath, which would tend to strengthen our constraints. However, the coupling to neutrinos can always be set to zero independent of the coupling to electrons: we assume χ couples to the neutrino only via the SU(2) L doublet, and χ can couple independently to e R and to e L . Thus, when evaluating BBN constraints, we ignore thermal contact with neutrinos in order to obtain the most conservative limits.

B. Effective number of neutrinos in CMB
Another powerful constraint applicable to a new light species is the effective number of neutrino species, N eff , as measured from CMB. N eff characterizes the contributions to the radiation energy density at recombination from relativistic species apart from photons, and is defined by In the absence of any other relativistic species, N eff 3. The SM actually predicts N eff = 3.046, accounting for the three neutrino species and for small effects due to non-idealities in the decoupling process [54,55]. This is consistent with analyses of Planck data, which find N eff 3.1 ± 0.2 [56]. Additional species are strongly disfavored. A single additional relativistic degree of freedom, i.e., a real scalar, is weakly consistent with current limits. However, CMB stage 4 experiments are expected to measure ∆N eff ≡ N eff −3.046 to within ±0.03, which is just sensitive enough to probe the minimum contribution from a real scalar at 1σ [57].
But a new species need not be relativistic at recombination to alter N eff . The introduction of a light DM species can change N eff by modifying the ratio of the photon and neutrino temperatures [31,35,58,59], and hence the ratio of energy densities in eq. (17). In the absence of additional species, the chemical decoupling of electrons and neutrinos takes place at T 0 D ≈ 2.3 MeV [60]. Any entropy transferred from DM to electrons after this decoupling leads to heating of the photon bath, and any entropy transferred to neutrinos heats the neutrino bath. If the new species transfers entropy differentially to the photon and neutrino baths at any time after the two baths decouple, the temperature ratio of the baths is modified. Note that ∆N eff thus depends on the relative size of the couplings to electrons and neutrinos, as pointed out in [58] and detailed extensively in [61].
Typically, the DM will transfer its entropy to one or both baths as a consequence of the conservation of comoving entropy density: when the DM becomes nonrelativistic while still in thermal equilibrium, the associated entropy must be transferred to any relativistic species to which it is still coupled. Thus, these species are heated when the DM becomes non-relativistic. Now, suppose that a sub-MeV DM species is coupled to electrons and neutrinos when T < T 0 D . If the DM species decouples from one and only one of these two relativistic species before it becomes non-relativistic itself, then the DM will reheat only one of the two baths, changing the temperature ratio. An exception to this rule occurs when the DM enters equilibrium with one bath below T 0 D , so that the DM accepts entropy of the same order that it loses upon decoupling later on [37]. We will discuss this scenario further in section V.
We now examine the calculation of N eff in detail. We will write T XY to denote the temperature at which species X and Y lose direct thermal contact, i.e., the temperature below which Γ(X ↔ Y )/H < 1 in our effective theory. The species X and Y might be kept in thermal equilibrium by a third species Z in our framework, i.e., through processes X ↔ Z and Z ↔ Y that remain active. We define T D to be the actual temperature at which electrons and neutrinos drop out of thermal equilibrium with one another once all inter-conversion processes have frozen out, including multi-step processes involving the dark species. Thus, in the standard scenario, MeV, but the introduction of a new species can keep electrons and neutrinos in thermal equilibrium at lower temperatures.
In particular, suppose that DM decouples from electrons instantaneously at a temperature T χe , and from neutrinos at a temperature T χν . If T 0 D < min{T χe , T χν }, then any entropy transferred to either photons or neutrinos can be shared between the two, so DM reheats these species equally, and the standard calculation is un-changed. However, if T χe < T 0 D < T χν , then χ remains in thermal contact with photons while relativistic, reheating the photon bath but not the neutrino bath. This increases the photon temperature, reducing N eff . Similarly, if T χν < T 0 D < T χe , then the reverse is true: DM reheats the neutrino bath, and N eff increases.
The only other possibility is max{T χe , T χν } < T 0 D , in which case χ acts as a thermodynamic mediator between electrons and neutrinos below T 0 D . In this situation, electrons and neutrinos remain in thermal equilibrium until the temperature falls below T D = max {T χe , T χν }. If the electron is still relativistic throughout this process, then the impact on N eff is determined by the ordering of T χe and T χν . But if T D < ∼ m e , the impact on N eff is quite different: photons and neutrinos are still in thermal contact while electrons become non-relativistic, so the electron also transfers some of its entropy to the neutrino bath. As we will see shortly, this can have a dramatic impact on N eff .
To calculate N eff , we follow the procedure described in [62]. In our scenario, the DM species is non-relativistic at recombination, so we assume that N eff is not modified by any additional degrees of freedom at recombination. Then, given the temperature ratio of the neutrino and photon baths at recombination, N eff is given by where N ν is the number of SM neutrinos (3). In turn, we can determine the temperature ratio from conservation of comoving entropy density.
Recall that the entropy density of a relativistic bosonic species i with g i internal degrees of freedom is given by 2π 2 g i T 3 /45. Away from the relativistic limit, denoting the true entropy density by s i , we say that this species has g s ≡ s i /(2π 2 T 3 /45) entropic degrees of freedom. Now, let g (γ) s and g (ν) s denote the entropic degrees of freedom in equilibrium with photons and neutrinos, respectively. Then g (α) s is given explicitly by where x i = m i /T α , and I indexes all species in equilibrium with species α (γ or ν). The sign in the denominator is determined by the statistics of species i. It can be shown [62] that if no entropy leaves the photon or neutrino baths after they decouple, then However, in our scenario, it is possible for entropy to leave one of the two baths below T D : suppose the DM decouples from one of the two baths above T D , and decouples from the other below T D , but while still relativistic. At this second decoupling, the DM's remaining ∆N eff as a function of the two decoupling temperatures Tχe and Tχν , assuming that χ is a Dirac fermion with mass 100 keV. Side and top panels show entropic degrees of freedom as a function of temperature. Gray shaded area indicates the region consistent with current data at 2σ. Labeled regions can be understood qualitatively as follows. Region A: Tχe, Tχν > T 0 D . Thus any entropy transferred by χ is shared between the γ and ν baths before they decouple. The standard calculation of N eff is unaltered. Region B: Tχe < T 0 D < Tχν . However, χ and e ± are relativistic at both decoupling events, so little entropy is transferred to either the γ or the ν bath. Region C: Now e ± becomes non-relativistic while still in thermal contact with the relativistic χ. The entropy ordinarily transferred by e ± to γ is now shared with χ, so γ is reheated less efficiently, and N eff increases. Region D: Here χ is relativistic below both T 0 D and Tχν , but becomes non-relativistic before Tχe is reached. Thus, χ reheats the γ bath exclusively upon becoming non-relativistic, decreasing N eff . Region E: χ becomes non-relativistic above both Tχν and Tχe, so it reheats both baths. The impact on N eff in this region comes from the delayed e ± -ν decoupling (see text). Region F: Tχe > Tχν , and χ is relativistic at Tχe. Thus, in addition to the delayed e ± -ν decoupling, χ reheats the ν bath. Region G: The electron and χ are relativistic at Tχe, so here the impact on N eff is due to χ reheating the ν bath. entropy leaves the bath to which it was last coupled. This only happens if T χe < T D ≤ T χν or T χν < T D ≤ T χe .
To account for this possibility, we modify the calculation of the temperature ratio as follows. Let us assume for the moment that T χe < T D ≤ T χν . Conservation of comoving entropy density in a thermal bath α amounts to the assertion that g (α) s | T T 3 a 3 is constant, where a is the scale factor. For T < T D , comoving entropy density is conserved in each bath except when T γ = T χe , so the temperatures of the two baths satisfy where the k i are constants. Generally, T rec < T χe , so Thus, to determine the temperature ratio, it is sufficient to identify the ratio k 1 /k 3 , which can be done in two stages. First, since T ν and T γ are equal at T D , we must have Similarly, at T χe , g s changes discontinuously while T γ is continuous in a. Thus, k 3 must satisfy where T ± χe denotes a temperature just above or below T χe . Now we have A similar calculation applies if T χν < T D ≤ T χe . Note that eq. (25) still assumes that χ does not enter equilibrium below T D , an exception we discuss further in section V. From eq. (25), it is easy to see why low DM decoupling temperatures can have a large impact on N eff . In the standard scenario, g (γ) s | T D includes photons (2) and relativistic electrons ( 7 8 × 4), which gives But if neutrinos and photons remain in thermal contact after electrons become non-relativistic, then g (γ) s | T D includes only photons, and the above ratio is increased to 1. This increases N eff by a factor of (11/4) 4/3 ≈ 3.9, already leading to N eff ≈ 12. If T χe < T χν , then χ reheats the photon bath when it becomes non-relativistic, reducing N eff . But if T χν < T χe , then χ reheats the neutrino bath, increasing N eff even further. The impact of relative decoupling temperatures on N eff is shown in fig. 2.
This approach assumes that the decouplings take place instantaneously, which is generally a good approximation. However, the approximation is poor when the decoupling process overlaps the range of temperatures during which a species becomes non-relativistic. In this case, the entropy of the species is changing rapidly, so it is difficult to estimate the amount of entropy transferred to other relativistic species before decoupling is complete. The temperature ratio can be determined precisely by numerical methods [see e.g. 61,63], and while that lies outside the scope of the present work, we note that instantaneous decoupling should be an effective approximation away from a narrow range of temperatures T χe and T χν , corresponding to a very small span of Λ EFT values in our parameter space.
To translate these results into constraints on the coupling between χ and electrons, we must make an assumption about the coupling between χ and neutrinos. If the coupling to neutrinos is very small, then χ may maintain thermal contact with electrons after decoupling from neutrinos. On the other hand, if the coupling to neutrinos is very large, then χ may remain in thermal contact with neutrinos after decoupling from electrons. In our case, we will assume that χ couples to ν exclusively by coupling to the lepton doublet (e L , ν e ) T . That is, we will assume that the χ-ν coupling is the same as the χ-e L coupling.
Even in this framework, the impact on N eff depends on the relative strengths of the χ-e L and χ-e R couplings. A non-zero coupling to e R tends to keep χ in equilibrium with electrons to lower temperatures, meaning that χ typically reheats the photon bath. This reduces the temperature ratio of eq. (20), producing ∆N eff < 0. However, if χ stays in equilibrium long enough to modify T D , then we can obtain ∆N eff > 0, as discussed above. Either way, increasing the coupling to e R only strengthens the effect, so we neglect this coupling to obtain conservative constraints. Note that this is different from our assump-tion in evaluating BBN constraints, where conservative constraints are obtained by neglecting the coupling to e L .

C. Non-thermal production
A viable model of DM must (partially) account for, but not exceed, the observed DM density of Ω DM h 2 0.12 [64]. If the DM is produced by thermal freeze-out, then a larger annihilation cross section reduces the relic density, so larger couplings conducive to direct detection are less likely to overproduce DM. But in the alternative scenario, if DM is produced out of equilibrium, the relic density increases with the annihilation cross section. In this case, overproduction is an important consideration.
If the DM species never attains thermal equilibrium with the SM, the abundance of DM will evolve toward its equilibrium value, but once Γ A /H < ∼ 1, the abundance will stay fixed. For renormalizable interactions, this outof-equilibrium production process is the standard freezein mechanism [65]. Out-of-equilibrium production has also been studied for non-renormalizable operators in the context of so-called ultraviolet freeze-in [66]. For temperatures below ∼ 10 MeV, within the constraints of our framework, such non-thermal production represents the only alternative to the freeze-out scenario.
The relic density of non-thermal DM is determined using the Boltzmann equation, much like the freeze-out case. The only difference is that the DM species χ is not in thermal equilibrium with e ± , and thus we cannot assume that χ has an equilibrium phase space density. Instead, we assume that the density of χ is negligible, such that the f 2 χ term drops out of the Boltzmann equation. In other words, starting from eq. (15), we approximate Y /Y eq 0, which gives Y (x) n eq (x) σ|v| (x)/H(x). It follows that the out-of-equilibrium yield can be estimated as As with freeze-out, the relic density in the non-thermal case is determined by the DM mass and couplings with SM particles. However, there is also a dependence on initial conditions in the form of x min and Y (x min ). In the freeze-out scenario, the abundance of DM in the early universe is simply the equilibrium abundance: equilibrium effectively erases the initial condition. But in the non-thermal scenario, equilibrium is never attained, so the dependence on the initial abundance is retained. Typically, when DM is produced by SM annihilations out of equilibrium, one calculates the relic density by fixing the DM density to zero at very early times and evolving non-thermally. This procedure requires that the interactions considered are renormalizable, in order for the production process to be modeled consistently at very high temperatures. Our effective operators are nonrenormalizable, so we cannot determine the relic density precisely in the non-thermal case: the result depends on the choice of UV completion. However, we can still place a lower bound on the relic density. We require that our effective theory is valid at scales below ∼ 10 MeV, so if we fix the abundance to some value at 10 MeV, we can determine the resulting relic abundance. In particular, by fixing the initial abundance to zero, we necessarily underestimate the relic density. This corresponds to a choice of x min and the condition that Y (x min ) = 0. With this initial condition, we can exclude models on the basis of their relic densities even when they never attain thermal equilibrium with the SM. Further, these constraints are determined entirely by conditions below T BBN , and are thus completely independent of the UV completion.
Note that if Λ EFT is sufficiently small, then even with this initial condition, the DM species will thermalize with the SM between T BBN and the present day. In this case, the relic density is set by the standard freeze-out paradigm, and eq. (27) is not valid. Even if the DM species does not quite enter thermal equilibrium, as long as it attains a non-negligible abundance, eq. (27) can significantly overpredict the relic density. Thus, while eq. (27) is useful to understand the qualitative features of the non-thermal relic density, we evaluate the constraint by numerically solving eq. (15).
As in the previous cases, we need to specify the coupling to neutrinos to perform these calculations consistently. Since the neutrino bath has a temperature comparable to the electron bath, a light χ can be effectively produced by neutrinos as well as electrons. Thus, a coupling between ν and χ can significantly affect the relic abundance. However, as with the coupling to electrons, the relic density is not monotonic in the coupling to neutrinos. If the DM never enters thermal equilibrium with any SM species, then a coupling to neutrinos tends to enhance the relic abundance by providing another production channel. On the other hand, if DM does enter equilibrium with neutrinos, then a larger coupling to neutrinos keeps it in equilibrium longer, reducing the relic abundance. However, at most of the points of interest in our parameter space, the constraint is driven by outof-equilibrium production, so we neglect the coupling to neutrinos when evaluating the relic density.

IV. CONSTRAINTS AND DETECTION RATES
The constraints we place on sub-MeV DM are relevant for direct detection experiments based on elastic electron-DM recoils. In principle, there are many such experiments, but they share several important features. Generically, electron recoil experiments prepare a lowtemperature collection of electrons for scattering with galactic halo DM, and by whatever mechanism, the experiment is sensitive to deposited recoil energies between some E min and E max . We calculate the detector sensitivity following [15], but the results are typical of electron recoil experiments with very low thresholds.

A. Estimation of the event rate
In the proposal of [15], the detector is constructed from an aluminum superconductor. At low temperatures, electrons move through the detector with velocities of order the Fermi velocity v F , and with the appropriate instrumentation, recoil energies as low as 1 meV may be detectable. We now review the calculation of the detection rate, following [15] and [67].
To compute the detection rate, we will consider scattering events at fixed recoil energy E R . We label the initial and final DM momenta by p 1 and p 3 , and the initial and final electron momenta by p 2 and p 4 . We do the same for the energies, so that We define the 3-momentum transfer by q = p 1 − p 3 . We denote 4-momenta by P i , and we write q = |q| and p i = |p i |. We denote the local DM number density by n χ , and the scattering rate by Γ = n e σv rel . The event rate per unit detector mass is where f χ (v χ ) is the local DM velocity distribution in the lab frame. We take the velocity distribution to be a Maxwell-Boltzmann distribution in the galactic frame with rms velocity 220 km/s and a cutoff at the halo escape velocity v esc 500 km/s. We then determine f χ (v χ ) by taking the Earth velocity to be 244 km/s in the galactic frame [68]. Now we turn to the evaluation of the scattering rate Γ(v χ , E R ). Observe that Γ not only contains the scattering cross section, but also accounts for the effects of Pauli blocking, effectively controlling the available phase space for scattering events. Following [67], we estimate Γ by Here, δ 4 P is a Dirac delta enforcing conservation of 4momentum; δ E fixes the recoil energy, setting where |M| 2 is the matrix element for the scattering process.
In many cases of interest, W is independent of the initial and final momenta of the target (p 2 and p 4 ), in which case the rate factorizes as where S accounts for Pauli blocking, and is given explicitly by (32) In our EFT, W is not generally independent of the target momenta. However, we can treat scattering in the non-relativistic limit, where such independence is guaranteed: the denominator in eq. (30) is independent of the momenta to first order, and can be replaced with 16m 2 χ m 2 e . The squared matrix element depends on the momenta only through the Mandelstam variables s and t, which have non-relativistic limits so |M| 2 is also independent of p 2 and p 4 to first order. Thus, for the remainder of this work, we will consider W to be a function of p 1 and p 3 only, and factorize the rate as in eq. (31).
We work in the low-temperature limit, where f FD reduces to a Heaviside step function, where E F ≈ 11.7 eV is the Fermi energy of aluminum. In this case, S(E R , q) can be evaluated explicitly. We perform the p 4 integral using the 3-momentumconservation delta function, and we use the remaining energy-conservation delta function to integrate over cos θ 2 . This leaves a 1-dimensional integral, This integral can be evaluated directly by comparing the arguments of the Heaviside functions. The result is where E 2 M = 2m e E R − q 2 2 /(4q 2 ) and E 2 S is given by To actually evaluate the rate in eq. (29), we change coordinates to (E R , q). Since there is no dependence on the azimuthal angle, we obtain and the limits of integration are q − < q < q + , where Under this change of coordinates, in the non-relativistic limit, t 2p 2 1 − 2m χ E R − q 2 . In particular, this means that W depends on p 1 and p 3 only through q, p 1 , and E R . EFTψ ψēe (g = 1). Background contours show scattering cross section, labeled as log 10 (σscat/cm 2 ). Black, DD: direct detection sensitivity (95% CL) with 1 kg yr exposure. Green, CMB: constraint from N eff . Orange, BBN: constraint from light element abundances. Blue, RD: constraint from relic density.
Then the differential scattering rate dΓ/dE R in eq. (28) is given by The limits of the E R integral in eq. (28) are set by the lower and upper thresholds of the detector, which we take to be 1 meV and 1 eV, respectively. Note that there are kinematical constraints on the minimum DM velocity (E 1 ) required to deliver a given recoil energy E R . Thus, the cutoff in the velocity distribution effectively imposes a maximum E R at fixed m χ .

B. Detection prospects and constraints by operator
We now examine our cosmological constraints in relation to the projected experimental reach for each of the operators in tables I and II. Figures 4 to 7 show cosmological constraints alongside projected 95% CL direct detection constraints with a 1 kg yr exposure. For purposes of discussion, we show constraints for O (ψ) SS in fig. 3.
All of the interactions considered for ψ a Dirac fermion can also be evaluated for ψ a Majorana fermion, and we do not consider matrix elements for Majorana fermions separately. Rather, we can directly relate our cosmological constraints on a Dirac fermion to the Majorana case. Whereas the relic density is controlled by n ψ Γ A = n 2 ψ σ |v| for a Dirac fermion, this expression doublecounts the phase space for a Majorana fermion. Since the relic density is inversely proportional to the annihilation rate, it follows that the relic density of a Majorana fermion is simply twice that of a Dirac fermion with the same mass and interactions [53,69].
The annihilation rate also sets the freeze-out temperature for a species in equilibrium with the SM, via the condition Γ A H. In general, Γ A ∼ Λ 4−k EFT for a dimensionk operator. All of our operators with DM a fermion are dimension-6, so to go from the Dirac case to the Majorana case, it is sufficient to make the replacement In principle, the value of N eff is also different in the Majorana case, but in nearly the entire excluded parameter space, ∆N eff is large compared with experimental uncertainty, sufficient to rule out a Majorana fermion as well as a Dirac fermion. Thus, in sum, the cosmological constraint curves in fig. 3 are shifted up slightly by a factor of √ 2 in the Majorana case, while the direct detection projections are unchanged.
Note that in the figure, the vertical axis is the suppression scale Λ EFT , effectively corresponding to inverse coupling. Thus, a stronger constraint line appears higher on the plot, and excludes the parameter space below. The vertical axis in each plot gives the value of Λ EFT alone, and the coupling g is taken to be 1. This is distinct from fixing g/Λ EFT or g/Λ 2 EFT , since we must have Λ EFT T BBN at all points regardless of the value of the coupling. Otherwise, the EFT would be applied outside its regime of validity.
The projected direct detection reach (DD, black) is generally the lowest line in each figure, i.e., the weakest constraint. The next line, stronger at low masses by greater than an order of magnitude in Λ EFT , is the constraint from light element ratios (BBN, orange). A comparable constraint is obtained from N eff (CMB, green). The final constraint is from overproduction of DM (RD, blue). Note that for some operators, there are narrow islands of parameter space where the N eff constraint is weakened. In these regions, the impact on N eff is transitioning between ∆N eff < 0 and ∆N eff > 0, as in fig. 2. Similarly, some regions with small Λ EFT are not ruled out by overproduction, since the DM thermalizes and freezes out at a lower abundance. However, we find no region of parameter space in which the projected direct detection constraints exceed all three cosmological probes.
Simplistically, this suggests that any model with a heavy mediator detectable by such an experiment is ruled out by cosmology. However, there remain possible exceptions to these constraints, as we discuss in the following section.

V. DISCUSSION AND CONCLUSIONS
In this work, we have derived cosmological constraints on a broad class of sub-MeV DM models that can be compared directly with detection prospects in electron recoil detectors. We now revisit the generality of our constraints, point out possible exceptions, and discuss the outlook for sub-MeV DM at electron recoil experiments.
Effectively, our goal has been to derive cosmological constraints on the scattering cross section between electrons and sub-MeV DM. Cosmology is mainly sensitive to the DM annihilation cross section, and in order to connect the two cross sections, we have produced these constraints in the context of an EFT. We have enumerated the possible thermal histories for a single DM species in this framework. If the DM is in thermal equilibrium with electrons at high temperatures, then light element abundances and N eff constrain the freeze-out temperature, and thereby constrain the interactions between χ and the SM. In the alternative scenario, if the DM is out of equilibrium at early times, a lower bound can be placed on the relic density, providing an independent constraint on the interactions. In both cases, a constraint is placed on the coupling between DM and electrons, assuming a specific form for the interaction.
In general, the form of the operator coupling electrons to DM affects the relationship between the annihilation cross section at early times and the scattering cross section today. Typically, then, constraints obtained by these methods are model-dependent. However, if the DM-SM mediator has a mass above ∼ 10 MeV, then our approach is quite general: our results are only sensitive to physical processes at lower temperatures, where the EFT is valid and cosmological history is well-established. Still, beyond the mediator mass, there are a few possible exceptions to the constraints derived here.
First, some of these constraints can be evaded with an extended dark sector. The overproduction constraint can be weakened, since such models provide a mechanism to deplete the DM relic density (e.g. through cannibalism, [70]). However, even in this case, the existence of a light DM species is enough for the BBN and N eff bounds to remain effective-adding additional dark degrees of freedom does nothing to improve the situation. One could still escape these constraints by assuming that a phase transition takes place in an extended dark sector between T BBN and the present day, such that the EFT is not valid in both epochs.
Another class of exceptions consists of models in which the dark species enters thermal equilibrium with the SM below T BBN , and thus below T D , the temperature of neutrino-photon decoupling. In this case, the entropy transferred to the SM bath upon freeze-out can be comparable to the entropy accepted upon equilibration, so the constraint from N eff can be circumvented [37]. This scenario is possible only in a very limited segment of the heavy-mediator parameter space, which we estimate as follows. We set the abundance of DM to zero at 1 MeV, and then determine the minimum value of Λ EFT below which DM thermalizes before the temperature drops to 0.5 MeV, thus still influencing BBN. Above this value of Λ EFT , it is possible to evade bounds from BBN and N eff , depending on initial conditions. Typically, this minimal value of Λ EFT is about one decade weaker than the BBN limit, and still out of reach of direct detection projections across most of our mass range.
Note that the overproduction bound already assumes an initial condition with zero DM abundance, so it cannot be evaded in this way. This is an example of the utility of the several overlapping constraints: the most conservative assumptions are different for each constraint, and correspondingly, exceptions apply differently as well. It is thus important to consider all of our constraints simultaneously, even in cases where the overproduction constraint appears to dominate. Finally, it might be possible to evade our constraints by taking some arbitrary linear combination of the effective operators in tables I and II. In principle, in this high-dimensional parameter space, there might be points for which interference of the matrix elements in tables III and VI conspires to reduce the DM annihilation or production cross section while preserving the scattering cross section. Then each of our cosmological constraints would be weakened, while the projected direct detection constraints would be maintained. However, in order for this to work, the Wilson coefficients would have to be engineered to produce such a cancellation.
In light of these constraints, the outlook for electron recoil experiments is brightest for DM masses 1 MeV < ∼ m χ < ∼ 1 GeV or for mediator masses m ζ 10 MeV. A light mediator certainly remains a possibility, but is subject to additional constraints [see e.g. 43]. The case of a light mediator is thus best studied in the context of simplified models, as in the analysis of [36]. Inelastic scattering may also improve direct detection prospects relative to cosmological constraints, and, of course, DM masses above ∼ 1 MeV remain an interesting target. However, if DM is dominantly composed of a single light species, and interacts dominantly with electrons via a heavy mediator, then cosmological constraints compromise the prospects of proposed experiments. ACKNOWLEDGMENTS BVL and SP are partially supported by the U.S. Department of Energy grant number DE-SC0010107. This work made use of the FeynCalc package [71,72] for evaluation of cross sections and the ColorBrewer tool [73] for selection of colorblind-friendly palettes. We thank Hiren Patel for cross-checking matrix elements and cross sections, and we thank Yonit Hochberg for valuable input regarding superconducting detectors. We thank Miguel Escudero, Rouven Essig, Juri Fiaschi, and Robert Scherrer for comments on an earlier version of this manuscript. We thank the Galileo Galilei Institute for hospitality while parts of this work were completed. Finally, we are especially grateful to Francesco D'Eramo, who originally conceived of this project and was deeply involved in most of the execution. Although circumstances compelled Francesco to leave the project later on, we remain deeply appreciative of his input and encouragement.     A vanish if φ is taken to be a real scalar. The matrix elements for scattering, φe − → φe − , are obtained from these by the substitution s ↔ t.      TT vanish if ψ is taken to be a Majorana fermion. For brevity, we define m 2 ± ≡ m 2 e ± m 2 ψ . The matrix elements for scattering, ψe − → ψe − , are obtained from these by the substitution s ↔ t.    TT vanish if ψ is taken to be a Majorana fermion.For brevity, we define m 2 ± ≡ m 2 e ± m 2 ψ and s ± i ≡ s ± m 2 i .