Making dark matter out of light: freeze-in from plasma effects

Dark matter (DM) could couple to particles in the Standard Model (SM) through a light vector mediator. In the limit of small coupling, this portal could be responsible for producing the observed DM abundance through a mechanism known as freeze-in. Furthermore, the requisite DM-SM couplings provide a concrete benchmark for direct and indirect searches for DM. In this paper, we present updated calculations of the relic abundance for DM produced by freeze-in through a light vector mediator. We identify an additional production channel: the decay of photons that acquire an in-medium plasma mass. These plasmon decays are a dominant channel for DM production for sub-MeV DM masses, and including this channel leads to a significant reduction in the predicted signal strength for DM searches. Accounting for production from both plasmon decays and annihilations of SM fermions, the DM acquires a highly non-thermal phase space distribution which impacts the cosmology at later times; these cosmological effects will be explored in a companion paper.


I. INTRODUCTION
One of the most well-studied mechanisms for setting the observed dark matter (DM) abundance is thermal freeze-out, where DM is in equilibrium with the Standard Model (SM) thermal bath at very early times. The DM abundance is then depleted through annihilations at later times until the DM drops out of chemical equilibrium. The appeal of this mechanism is that the final relic abundance is generally independent of the high-temperature initial conditions at reheating. Furthermore, producing the observed relic abundance requires a particular thermally averaged annihilation cross section in most thermal freeze-out scenarios, σv ∼ 10 −26 cm 3 /s. This weakscale cross section provides a target that can be probed by direct and indirect detection experiments. Assuming the relic abundance is set by annihilations to SM particles, then consistency with Big Bang Nucleosynthesis (BBN) generally requires that thermal freeze-out candidates have masses m χ 1 MeV [1][2][3]. The appealing simplicity of this scenario has led to an enormous number of DM searches targeting the thermal freeze-out mechanism, with a particular emphasis on weakly interacting massive particle (WIMP) candidates in the m χ ∼ GeV−TeV mass range. More recently, there has been a growing interest in m χ ∼ MeV−GeV thermal candidates where interactions with the SM or within a hidden sector deplete the DM density to the observed value [4][5][6][7][8][9][10][11][12][13][14][15][16].
The freeze-in mechanism for DM production is a compelling alternative to thermal freeze-out, where DM is instead produced by feeble, sub-Hubble interactions of SM particles [17][18][19][20][21][22]. If the dominant freeze-in process is annihilation of SM particles into DM via a light mediator, then many of the appealing features of thermal freeze-out are maintained. For annihilation through a mediator lighter than the DM, the thermal cross section typically scales as σv ∼ g 2 χ g 2 SM /(4πT ) 2 where g χ is the mediator-DM coupling, g SM is the mediator-SM coupling, and T is the SM temperature. With this scaling, DM freeze-in dominantly occurs at the lowest temperature where the process is kinematically accessible, and thus the mechanism is not sensitive to the reheat scale. 1 Freeze-in through a light vector mediator has emerged as a key benchmark for sub-GeV direct detection experiments. Producing the observed DM relic abundance implies a tiny value for the coupling constants, which is difficult to target with accelerator searches. However, sufficiently light mediators give rise to scattering cross sections that scale as σ ∝ 1/v 4 for relative velocity v, implying that the kinematics of the Milky Way (where v ∼ 10 −3 ) can enhance the detectability of DM coupling to a light mediator. If the mediator also couples to charged SM fermions, then the DM can scatter off of electrons or nuclei and may be detectable with the next generations of direct detection experiments [23][24][25][26][27][28][29][30][31][32][33][34] (see also Ref. [35] for a recent review). Indeed, recent experimental results by SENSEI [36,37], SuperCDMS [38], and DarkSide [39] are demonstrating significant progress towards achieving the sensitivity needed in the MeV-GeV mass range. It was also shown recently that Xenon1T [40] is for the first time constraining freeze-in in the GeV-TeV mass range [41].
In the keV−MeV DM mass range, freeze-in is the leading scenario that could be tested by proposed low-threshold direct detection experiments. Refs. [42,43] studied the possible direct detection cross sections in models of sub-MeV DM, finding that it would be difficult to observe thermal freeze-out scenarios (even purely within a dark sector) due to a combination of BBN, CMB, fifth force, and stellar emission constraints. Obtaining accurate predictions of freeze-in is thus an important step in the program to search for low-mass DM. While freeze-in from electronpositron annihilations via a light vector mediator has been studied in the past [23,44], in this work we thoroughly explore a previously overlooked production mechanism: freeze-in through plasma effects. The contribution of plasma effects to dark sector thermalization was estimated earlier in Refs. [45,46] and the effect on freeze-in via a heavy mediator was recently considered in Ref. [47] as we were in the late stages of completing this work, but it was not included in previous studies of freeze-in through a light vector mediator. We find that the plasma production of DM is a dominant channel for sub-MeV DM masses, and will therefore restrict our discussion to this mass range. The additional contribution to the relic abundance implies that the target cross section for direct detection is lower by roughly an order of magnitude for the lowest experimentally accessible DM masses.
The rest of this paper is organized as follows. We begin in Section II by reviewing the arguments for the simplest viable freeze-in models in the keV-MeV mass range: either pure millicharged DM arising from a DM hypercharge or effectively millicharged DM that is coupled to an ultralight dark photon mediator. These two scenarios are almost phenomenologically identical, with the key difference being that DM-DM scattering can be parametrically larger when dark photon interactions are present. These DM candidates have recently received considerable attention in the context of the anomalous 21 cm global signal [48][49][50][51][52]. In Section III we compute the DM relic abundance from freeze-in via a light mediator. We include the effects of plasmon decays for the first time, and show the impact for direct detection. We then present the calculation of the phase space distribution for freeze-in DM in Section IV. A summary of our results can be found in Section V. In a companion paper [53], we will apply the calculations of the phase-space distribution to cosmological observables, showing that the cosmic microwave background (CMB) and probes of large-scale structure (LSS) provide a strong complementary test of DM freeze-in for m χ ∼ keV−MeV. In particular, we find that existing cosmological constraints restrict m χ tens of keV for freeze-in via a light mediator, and it will be possible to probe even higher masses with planned experiments.

A. The case for light vector mediators
The simplest observationally viable models for sub-MeV freeze-in through a light mediator can be divided into two classes, where (1) the DM only has interactions mediated by the SM photon or (2) the DM has interactions with an ultralight kinetically mixed dark photon. We note that models of millicharged DM [45,54] can fall under either category: they can arise as a limit of the dark-photon model where the dark photon is nearly massless, or they could be present as Dirac fermions with a tiny hypercharge. 2 For sub-MeV freeze-in to be relevant for direct detection, vector mediators are the only observationally viable option due to stringent constraints on other light mediators with the requisite couplings to the SM, as outlined below. For direct detection of freeze-in, the mediator masses must be sufficiently small compared to the typical momentum transfer for scattering processes. If the mediators are heavier, then they do not give rise to the v −4 enhancement that would render extremely feeble DM-SM interactions detectable on Earth. For nuclear recoils the relevant momentum scale is set by galactic kinematics q ∼ m χ v ∼ 10 −3 m χ , while for electron recoils the typical electron momentum in the target material is most relevant q ∼ αm e ≈ 4 keV, where m e is the electron mass and α is the electromagnetic fine structure constant. Thus for sub-MeV DM, the experimentally relevant mediators have masses below O(1) keV.
Assuming an annihilation cross section of SM fermions into DM with the form σv ∼ g 2 χ g 2 SM /(4πT ) 2 , the relic abundance can be estimated as where M Pl = 1/ √ 8πG is the reduced Planck mass and we assumed T ∼ MeV. Then for m χ ∼ MeV, we find that g χ g SM 10 −12 to saturate the relic abundance. This order-of-magnitude estimate is in agreement with more detailed calculations below. Since obtaining the relic abundance from freeze-in requires g χ g SM ∼ 10 −12 , g SM must be greater than 10 −12 if we require the dark sector to be perturbative (i.e. g χ 1). Weakly coupled, sub-keV mediators can be emitted in stars, affecting their luminosity and lifetime. The observed properties of stars lead to strong bounds on such mediators, which we summarize here (see also Refs. [42,43] where these bounds are collected and discussed in the context of sub-MeV DM models): • Scalars and pseudoscalars coupled to electrons − The strongest bound on a light scalar with interaction g φee φēe comes from helium ignition in red giants, with g φee 7 × 10 −16 for sub-keV masses [57]. For a sub-keV pseudoscalar, observations of white dwarfs give typical constraints of g aee 2 × 10 −13 [58][59][60]. A caveat for most stellar emission bounds is that when the coupling is increased, the new particle may be trapped within the star and would not lead to anomalous energy loss. However, this would still affect energy transport in the star, which can be constrained for the range of couplings relevant for freeze-in through this mediator [61,62].
• Scalars and pseudoscalars coupled to nucleons − Similar to the case of mediators coupling to electrons, red giants constrain g φnn 10 −12 for a scalar [57] and g ann few × 10 −10 for a pseudoscalar [58,63]. While the latter coupling appears at face value to be sufficiently large, freeze-in through baryons is largely suppressed after the QCD phase transition due to the low baryon number density. Therefore, in this case our estimate for the minimum g SM with T ∼ 1 MeV is much too low and freeze-in would have to occur with a larger value of g SM that is in tension with stellar bounds.
• Scalar mixing with the Higgs − The bounds here are similar to those in the two previous cases, and it has been shown in Ref. [64] that freeze-in through this portal is only a viable mechanism for producing all of the DM for DM masses above a few hundred MeV.
• Kinetically mixed dark photon − In this case, the stellar constraints on g SM decrease linearly with the mediator mass for masses below ∼ 100 eV [65,66] because of the in-medium plasma mass suppression of producing dark photons from SM interactions, as detailed in Eq. (5) and the surrounding discussion in Section II C. From the collected bounds on dark photons from Refs. [67], a dark photon can have g SM > 10 −12 when its mass is well below 1 eV. At even lower masses, the coupling could be ∼ 10 −3 for masses below 10 −14 eV.
• B − L vector − Stellar constraints on a B − L vector are similar to that for the dark photon. However, for eV-scale and lighter mediator masses, a B −L vector is also strongly constrained by fifth force searches (e.g. [68,69]), which limits the mediator-SM coupling to below 10 −12 .
Since we are focusing on the simplest benchmarks for direct detection, we do not consider more exotic possibilities with additional particles and interactions. From the bounds on new particles with the couplings described above, we conclude that freeze-in through a light mediator is viable either when the mediator is (1) the SM photon, and the DM has a tiny electric charge, or (2) when the mediator is an ultralight kinetically mixed dark photon. We discuss these two closely related scenarios in the rest of the section. In both cases, DM has an effective charge Qe (or millicharge Q) with respect to the SM photon. This parameter determines the relic abundance, irrespective of which of the two models is under consideration. Both models allow for heat and momentum transfer between SM particles and DM during epochs when the typical relative velocities are low (as discussed in Section IV C), which is relevant to observations of the CMB [70][71][72][73][74][75][76] and the cosmological 21 cm global signal [48][49][50][51][52]. The main phenomenological difference between these two possibilities is that DM-DM scattering via a dark photon can be parametrically larger than DM-DM scattering mediated by the SM photon, as discussed below. If present at a sufficient level, the DM self-scattering can play an important role in determining the DM phase space distribution at late times, well after freeze-in.

B. DM with photon-mediated interactions
If the DM is a Dirac fermion χ with a tiny hypercharge Q Y (the only gauge-invariant, renormalizable operator leading to a bare millicharge), then it can interact via the SM photon. After electroweak symmetry breaking, the DM obtains an electric charge given by eQ Y ≡ eQ (taking the convention where the Gell-Mann Nishijima formula reads Q = I 3 + Y ). Although there are also Z-mediated DM interactions, they are negligible for the relevant epochs where T m Z . This gives the simplest model of millicharged DM. It is difficult to incorporate such matter content into a Grand Unified Theory (GUT) [77]; however, this scenario is economical in that it requires that no additional particles be introduced to the SM aside from the DM itself.
The possibility that this DM candidate obtains its relic abundance by thermal freeze-out has been considered before in Ref. [78], where it was shown to be excluded by structure formation when all of the DM is produced this way. Thus, freeze-in is the simplest remaining possibility for producing this DM candidate, with g χ = eQ and g SM = e in the language of the previous subsections.
There are stellar emission bounds on this DM candidate because the DM can be pair produced by the decay of plasmons in stars, leading to additional energy loss. These bounds are shown as the shaded region in our summary plot, Fig. 8. Constraints on DM pair produced in SN1987a were derived in Refs. [45,79] and require Q 10 −9 for m χ up to a few MeV, which does not impact freeze-in. However, there are constraints for m χ below O(10) keV from emission in white dwarfs, horizontal branch stars, and red giants (see Appendix of Ref. [46]). Note that the range of m χ where stellar emission can constrain freeze-in is exponentially sensitive to assumptions about temperatures within the stars. In addition, the bounds derived are applicable in the weak coupling limit where the DM escapes cleanly from the star. For sufficiently large Q, DM emission could contribute to energy transport within the star and the effects have not been carefully studied in this regime. The couplings for freeze-in are large enough that they could be in this regime and stellar bounds on freeze-in should be regarded with care.
The relevant interactions for the relic abundance and phase space distribution in this model are SM annihilations and plasma decay into the DM. DM-SM scattering can become important at late times but, as we discuss in Section IV C, the effect must be small to be consistent with limits from the CMB. The DM self-scattering cross section is proportional to Q 4 , and we find it to be irrelevant for the phase space. Finally DM-photon scattering is also proportional to Q 4 and is not enhanced in the low-velocity limit, so it is also irrelevant.

C. DM with dark photon interactions
We next consider Dirac fermion DM coupled to a kinetically mixed dark photon A , with the vacuum Lagrangian given by where A is the SM photon, κ is the kinetic mixing parameter and χ is Dirac fermion DM. For the purposes of this discussion, we consider Abelian kinetic mixing, noting that non-Abelian kinetic mixing is also possible [80,81]. The mixing parameter κ could have any number of origins; for instance, it could be generated as a result of loop diagrams with heavy matter fields charged under both A and A [82] or from certain compactifications of type IIB strings [83,84]. Since the kinetic mixing term is a marginal operator, we take the point of view of a bottom-up effective field theory and we will treat it here as a small free parameter without specifying its origin. In this model, the combination of couplings relevant for the relic abundance is g χ g SM = g χ κe.
As discussed in Section II A, the dark photon mass must satisfy m A 1 eV in order to give a sufficient coupling for freeze-in while also evading existing bounds on stellar energy loss [67]. However, the requirements are even more stringent because unlike the model presented in Section II B there could be large A -mediated DM self interaction. For m A < eV, the mediator would be light enough to give rise to v −4 enhanced DM self-scattering in astrophysical environments, with a rate proportional to g 4 χ . Furthermore, as mentioned before, the freeze-in relic abundance is determined by the product g χ κe, meaning that large g χ can be compensated by reducing κ to give the same observed relic abundance. Thus a sizable DM self-interaction is possible, and could be relevant to astrophysical probes of self-interacting DM (SIDM). The effects of SIDM are typically parameterized by the momentum-transfer self-scattering cross section, which in the limit of a very light vector mediator is given by [85] where θ CM is the scattering angle in the center-of-mass (CM) frame, σ χχ is the self-interaction cross section, and α χ is the dark equivalent of the electromagnetic fine structure constant, α χ ≡ g 2 χ /4π. Typical bounds on SIDM require σ χχ /m χ < 1 − 10 cm 2 /g for systems ranging from dwarf galaxies where v ∼ 10 −4 to merging clusters where v ∼ 10 −2 (for a recent review, see Ref. [86]). While few simulation-based studies of self-interactions have been done in the ultralight mediator limit (see for instance Ref. [87]), we can estimate the expected bound. Taking the more restrictive limit of σ χχ /m χ ∼1 cm 2 /g, the bound is Since κeg χ 10 −12 is needed for sub-MeV freeze-in, the SIDM bounds imply that the kinetic mixing is κ 10 −7 for MeV-scale DM. For sub-eV dark photons, such large kinetic mixing is only possible when m A 10 −10 eV [67]. For even lighter DM, g χ is even more restricted so κ 10 −5 is required for freeze-in, which is possible when m A 10 −14 eV. Therefore, we are required to consider an "ultralight" dark photon [42]. Note that black hole superradiance constrains dark photons being present in the mass spectrum (in the small-coupling limit) between ∼ 10 −14 − 10 −11 eV and preliminarily between ∼ 10 −19 − 10 −17 eV [88].
Such a light dark photon is phenomenologically equivalent to the massless dark photon limit for all processes considered in this paper because the m A is much lower than the effective in-medium photon mass m A in the early universe. Then, following Appendix D of Ref. [42], the vacuum Lagrangian in Eq. (2) is modified with an additional term m 2 Rotating away the mixing term in the presence of m A and m A and rewriting in terms of the mass eigenstatesÃ andÃ , the in-medium Lagrangian is given by From this, we see that when m A m A , the interaction terms above reduce to meaning that DM has an effective millicharge parameter Q = κg χ /e, and the interactions are identical to those for a massless dark photon. Note that this suppression of the A -SM coupling in the m A m A limit is the source of the in-medium (plasma mass) suppression of the stellar constraints on dark photons [65,66] discussed in Section II A. Also note that this suppression 3 For simplicity we consider a constant m 2 A for the schematic purposes of this discussion, although the photon polarization tensor Π µν ( q, ω) (which gives rise to the in-medium effective mass) depends on the photon momentum q, energy ω, polarization, and thermal properties of the medium. For an on-shell mode with ω ∼ | q|, m 2 A would correspond to the plasma mass, as discussed in Section III B. For scattering processes with a highly off-shell mode, | q| ω, m 2 A is given by the Debye mass [89].
means that the dark photon is not abundantly produced by SM interactions in the early universe and does not contribute to the effective number of relativistic species, N eff . In the exactly massless A limit, we are free to perform a field redefinition on A → A + κA in the vacuum Lagrangian, Eq. (2), which eliminates the kinetic mixing term and generates a DM interaction term g χχ γ µ χ(A µ + κA µ ), which is again identical to having a millicharge Q = κg χ /e under U (1) EM .
The model considered here thus provides another realization of millicharged DM, and all of the stellar constraints discussed in the previous section apply. The only difference is the additional DM self-interaction via the A , which potentially leads to sizeable self-interactions.

III. RELIC ABUNDANCE FROM FREEZE-IN
Here we compute the relic abundance of DM from freeze-in. We begin by reproducing the contribution from annihilation of SM fermions ff → χχ that was previously calculated in Refs. [23,44]. Because freeze-in is peaked at low temperatures and this paper concerns sub-MeV DM, electrons are the primary source of DM for this channel; in the rest of this section we explicitly refer to freeze-in off electrons, noting we have numerically checked that adding heavier fermions (for instance muons) to the calculation changes the results by less than 1%. In addition to freeze-in off electrons, there is a contribution from plasmon decays, γ * → χχ, which we calculate for the first time. Photon annihilation into DM γγ → χχ is suppressed by an additional factor of Q 2 and can be safely neglected.
In what follows, we take the observed present-day relic DM abundance to be ω c ≡ Ω c h 2 = 0.12 [90]. After freeze-in, the DM density should scale like a −3 and it is common practice to compare this to another quantity that has the same scaling irrespective of changes to the SM bath temperature. In this work we choose to compare the number density to the entropy density. Taking the present-day CMB temperature to be 2.73 K, the observed yield is then For m χ 1 keV, the DM yield is much lower than the order unity yield for relativistic species, such that DM contributes negligibly to N eff . This is in contrast to other DM models, such as thermal freeze-out, where sub-MeV DM would generically inject a considerable amount of entropy to the photon or neutrino sectors and would violate observational bounds on N eff . 4 The low DM occupation number also implies that it is possible to self-consistently ignore backreactions that would reduce the DM number density, namely DM annihilation to electrons and inverse decays to plasmons. For instance, if we ignore the back-reaction, the solution for the number density of DM is significantly lower than the electron number density during the entirety of freeze-in in spite of the fact that the latter is becoming Boltzmann suppressed. Depletion of the DM number density through annihilation to dark photons χχ → γ γ is negligible for the same reason. In what follows, we solve the 0 th moment of the Boltzmann equation ignoring backreactions, noting that we have numerically checked that they are negligible. The relevant equation is then Here we are using a as our time variable. The relationship between a and the SM temperature T (which determines the DM production rate) is not adiabatic during freeze-in because the electrons are leaving the thermal bath at this time; this is discussed further in Appendix A. Note that we are solving for the total DM density which includes both χ andχ in the matter budget; assuming zero DM chemical potential, n DM = 2n χ = 2nχ, which accounts for the factor of two in Eq. (8). 5

A. Annihilations
In computing the DM relic abundance from annhilations of electron-positron pairs, we treat the two scenarios discussed in Section II as indistinguishable in the limit that m A → 0. We also ignore the in-medium photon mass for this process, which we find to be a percent level effect for s-channel annihilations happening at the relevant range of temperatures. In this limit, the matrix element squared is where we sum over both initial and final spin degrees of freedom (d.o.f.) without averaging and where Q is the effective millicharge in the dark photon case, Q = κg χ /e. The thermally averaged cross section appearing in Eq. (8) for this process is given by 3 . We assume that from the onset of freeze-in, the electrons have entered the non-relativistic regime where their phase space is given by a Maxwell-Boltzmann distribution with temperature T and zero chemical potential. As we will show, sub-MeV DM freeze-in through the annihilation channel is most effective at temperatures T m e where the effects of Fermi-Dirac statistics can be neglected. We also ignore Pauli blocking of the DM due to its low occupation number.
To evaluate the thermal cross section, we note that the primordial plasma has a preferred rest frame (where bulk motions average to zero), which breaks Lorentz invariance. The phase space 5 This factor is related to the usual factor of 1/2 that appears in the Boltzmann equation for Dirac fermions [92,93]; however, unlike the ordinary case of thermal DM, the change in the comoving DM density for freeze-in is independent of the DM number density (i.e. there is no factor of n 2 DM appearing in Eq. (8)) which accounts for the factor of four difference. factors of Eq. (10) are evaluated in a frame that is comoving with the plasma. Practically, we can perform the integration by inserting factors of unity, where q 12 is the effective bulk 4-momentum of the particles labelled 1 and 2 and s 12 can be thought of as the effective (Lorentz invariant) mass-squared of a single particle with that bulk 3-momentum and energy (i.e. here E 12 = s 12 + q 2 12 ). Inserting such a factor into Eq. (10) gives The integral over p χ and pχ does not depend on the frame of q χχ , so the two-body phase space of p χ and pχ can be evaluated in the CM frame of q χχ . We define and insert this into the expression for the thermally averaged cross section Again, we can evaluate the integral over p e + and p e − in the center-of-mass frame. Defining the thermally averaged cross section becomes We can write this result in terms of the first order modified Bessel function of the second kind where we have dropped the subscript on the integration variable s. Note that s is restricted to s > 4 max m 2 e , m 2 χ . The procedure above provides an alternate derivation of the well-known results from Ref. [92], and we have validated this method here because we use similar techniques to derive the full collision term for annihilation in Section IV A.

B. Plasmon decay
The early Universe is an optically thick plasma where photons acquire an in-medium mass; this can be understood classically as arising from the electrons' oscillatory response to a propagating electric field and the dynamical shielding of that electric field. This effective mass is also manifest in the photon propagator and the polarization vectors of external photon legs in the medium; in other words, the photon mass and wavefunction are renormalized in the plasma. The effective masses and dressed polarization functions for the transverse and longitudinal "plasmon" modes are shown in Fig. 1 and explicit formulae are provided in Appendix B. The effective mass for plasmons is closely related to the plasma frequency. For a relativistic plasma at zero chemical potential, the plasma frequency is ω p = eT /3 ≈ 0.1T where e is electric charge.
Plasmons can undergo decay provided that it is kinematically allowed. For instance, plasmons can decay to neutrino pairs through mixing with the Z boson [94]. Plasmons cannot decay to charged particles in the SM because their effective mass is also renormalized in the medium and it is always kinematically forbidden. However, this is not the case for millicharged DM where corrections to the mass are suppressed by powers of Q.
The effective matrix element that captures plasmons decaying to DM is where˜ µ (k) is the dressed polarization vector for the longitudinal and transverse plasmon modes as detailed in Appendix B, where we work in Coulomb gauge. We express this process in terms of the DM effective millicharge Q and in Appendix C we show explicitly that decaying through a dark photon gives the same effective matrix element in the limit m A → 0. In squaring and summing over polarizations, only the diagonal terms (LL, ++, and −−) contribute, where the photon four-momentum is given by K µ = ω(k), k µ with appropriate dispersion relations for transverse and longitudinal modes ω t (k) and ω (k) (see Appendix B), the DM fourmomentum is given by (E χ , p χ ) µ , θ is the angle between k and p χ , and Z t (k) and Z (k) are wavefunction renormalization factors (shown in Fig. 1) that are related to the dressed polarization vectors for the transverse and longitudinal modes. The thermally averaged decay rate is and can be evaluated directly. Taking the plasmons to be Bose-Einstein distributed, the longitudinal and transverse contributions to this rate are where the effective plasmon masses are m (k) 2 = ω (k) 2 − k 2 for the longitudinal modes and m t (k) 2 = ω t (k) 2 −k 2 for the tranverse ones. The final integrals over k can be computed numerically and the total plasmon contribution to decay is dominated by the transverse modes (note that we are working in Coulomb gauge). This is because the longitudinal mode has a finite range of k over which it can propagate, meaning that it has less available phase space than the transverse mode which has no restriction in k. Furthermore, the longitudinal mass and renormalization factors fall steeply within the range of k where this mode can propagate.

C. Couplings for freeze-in
In solving the zeroth moment of the Boltzmann equation for the DM relic abundance, we find that the relative contributions from e + e − annihilation and plasmon decays are starkly different in different mass ranges, as illustrated in Fig. 2. This can be understood by considering the fact that freeze-in is dominant at low temperatures, provided that it is kinematically allowed and that the population the DM is freezing in from has a sufficient abundance. For sub-MeV DM, freezein from e + e − annihilation is always kinematically allowed and this process only ends when the electron number density becomes Boltzmann suppressed, namely T m e . Meanwhile, the plasmon abundance is not Boltzmann suppressed but the mass runs with temperature, so freeze-in through plasmon decay ends when it is no longer kinematically allowed, namely when m γ * ∼ ω p = 2m χ . Since ω p ≈ 0.1T in the relativistic limit, plasmon decay to millicharged DM shuts off at an earlier time compared to annihiliation. These two criteria are shown in Fig. 2 and indeed we see that plasmon decays are more dominant in determining the relic abundance for lower mass DM because the decays are active for a longer period of time.
In terms of the effective millicharge needed to produce the observed DM relic abundance, we find that including plasmon decays leads to a significant reduction in coupling for keV-mass DM while the effect is small once m χ = MeV. The change to the freeze-in benchmark for direct detection is shown in Fig. 3, where the cross section for electron recoils is Here µ χe is the electron-DM reduced mass, µ χe = m e m χ /(m e + m χ ). At the lowest mass where proposed low-threshold direct detection experiments are sensitive, the plasmon decay channel for DM production lowers the expected signal strength by roughly an order of magnitude. It has been noted in the literature [95][96][97] that millicharged DM could be efficiently accelerated in supernova remnants, which would lead to an accelerated component of dark cosmic rays and eject DM from the disk. Both of these effects can lead to substantial changes to the predicted direct detection rates and sensitivities of proposed experiments shown above. However, the conclusions are highly sensitive to aspects of cosmic ray physics which are not fully understood, such as the injection of particles into the diffusive shock acceleration process. The predictions would also be sensitive to whether the DM obtains its effective millicharge through a kinetic mixing portal; in this case, the dark photon mass and couplings can affect the acceleration, and an exploration of these effects is beyond the scope of this work.

IV. DARK MATTER PHASE SPACE DISTRIBUTION
Since freeze-in DM is so weakly coupled to the SM, it does not thermalize with the SM during freeze-in and the phase space distribution can deviate substantially from a thermal distribution. While this has no clear impact on direct detection, since galaxy assembly is expected to significantly alter the DM velocity distribution, it does affect DM free-streaming and DM-SM scattering in the early universe. Here we compute the full phase space distributions needed to determine the cosmological observables; the signatures, constraints, and detection prospects will be presented in a companion paper [53]. We must solve the full Boltzmann equation in an expanding background, given by where C(p χ , t) is the collision term, which encapsulates all interactions that affect the phase space. At early times, the interactions that determine the phase space evolution are e + e − annihilation and plasmon decay. We have checked numerically that heavier fermion annihilation processes (for instance the annihilation of muon-antimuon pairs) affect the phase space by a negligible amount because they occur only at early times when freeze-in is less efficient. Scattering has a negligible impact on the phase space during freeze-in since the DM occupation number is much smaller than that of electrons or plasmons. Neglecting the small effect of scattering during freeze-in, the collision term is independent of f χ to leading order and the Boltzmann equation can be solved by direct integration [98], Here the factors of a in the integrand keep track of redshifting of momentum due to expansion. We use the scale factor a as our time variable rather than the common choice of using the SM temperature because it is not evolving adiabatically as the electron-positron pairs leave the bath during freeze-in. The temperature evolution and the evolution of the Hubble parameter are detailed in Appendix A.
After freeze-in ends, the DM momenta redshift and the phase space distribution is constant in comoving momentum. However, at late times DM-SM and DM-DM scattering eventually can become important since the scattering cross sections are peaked at low relative velocities. The effects of DM-SM scattering on the phase space are generally negligible for the allowed parameter space, but DM self-scattering can lead to thermalization of the DM phase space distribution. Whether this occurs is model-dependent, and we discuss the conditions for this to occur in Section IV D.

A. Phase space from annihilation
The computation of the full collision term from annihilation proceeds similarly to the computation of its zeroth moment. Once again, inserting a factor of unity as defined in Eq. (11), we find where Eχ = m 2 χ + p 2 χ + q 2 e + e − − 2p χ q e + e − cos θ, E e + e − = s e + e − + q 2 e + e − and θ is the angle that q e + e − makes with the unconstrained, unintegrated p χ . Defining x ≡ cos θ and dropping the subscript on the bulk electron momentum, we find Requiring that x ∈ [−1, 1] and switching integration variables, Then, to solve for the final phase space from annihilation, we can combine Eqs. (25) and (28). Note that because p χ is fixed (rather than an integration variable), s in the above integral is restricted to s > max 4m 2 e , 2m χ (E χ + m χ ) unlike in the integral for determining the thermally averaged cross section. The resulting evolution of the phase space distribution is shown in the left panel of Fig. 4.

FIG. 4.
A comparison of the phase space evolution of DM being produced by e + e − annihilation (left) and γ * decay (right) at m χ = 40 keV. The momenta shown here are comoving, P χ ≡ ap χ where a = 1 corresponds to T = 1 MeV. The phase space is normalized arbitrarily for the purposes of comparing the P χ -dependence side by side. Over time, the comoving phase space converges to its final frozen-in shape. The phase space from annihilation is similar to that of the thermal electrons from which they inherit their kinematics. Meanwhile, the phase space from plasmon decay is highly peaked at low P χ because freeze-in through this channel occurs predominantly at threshold when ω p ∼ 2m χ and the decay is peaked when the plasmon is "at rest," k → 0.

B. Phase space from plasmon decay
The collision term from plasmon decay, proceeds through direct computation. We find where the limits of the k integral are determined by the requirement that x 0 = (2E χ ω ,t (k) − m ,t (k) 2 )/2kp χ lies in the range [−1, 1]. The limits of integration cannot be solved for in closed form because of the nontrivial dispersion relations, so the phase space must be determined numerically.
The evolution of the phase space from plasmon decays is shown in the right panel of Fig. 4, and our results for the combined phase space can be found in Fig. 5. The distributions are noticeably nonthermal due to plasmon decays. Fig. 6 compares the average momentum and momentumsquared of the DM to the SM photons, which serves as a useful metric to determine the DM free-streaming and suppression of the growth of structure. The phase space is normalized to the comoving DM relic abundance for each mass depicted. The plasmon contribution dominates more at low masses than at high masses because freeze-in through this channel persists for longer at lower masses, ending when the plasmon mass is at threshold, ω p ∼ 2m χ . Also shown (dashed lines) are the phase space distributions that would arise if the DM could thermalize within its own sector, conserving P 2 χ for non-relativistic DM.

C. Effect of DM-SM scattering
We argue here that the effect of DM-SM scattering on the DM phase-space distribution is small from freeze-in until the onset of recombination. The relevant quantity is the momentum-transfer rate, which we estimate in the limits where the DM is relativistic and non-relativistic. We do not consider scattering by relativistic, charged SM particles because this is only relevant for electrons during freeze-in; during freeze-in, the number density of DM is many orders of magnitude smaller than the number density of electrons and the effect of electron-DM scattering is suppressed by n χ /n e relative to the dominant effect of electron-positron annihilations on the phase space. As outlined below, DM-SM scattering becomes more important at low velocities, corresponding to later cosmological times. This can affect CMB anisotropies and the cosmological 21 cm signal, and we provide more detailed calculations in that context in our companion paper [53].
In the limit of relativistic DM scattering with non-relativistic SM particles (the case after freezein until T γ ∼ m χ ), the differential cross section with respect to the center-of-mass scattering angle θ CM is given by where p CM ≡ | p CM | is the momentum in the CM frame. Here we have taken p χ m e , which is a good approximation after freeze-in has ended. In this approximation, the dependence on the SM particle mass drops out, making scattering with electrons and protons equally important (we refer to them collectively as "baryons," in the remainder of this discussion, hence the subscript b  6. A comparison between moments of the DM phase space and the SM photon phase space as a function of DM mass. For reference, the moments for the SM photon are p γ = 2.7 T γ and p 2 γ = 10.35 T 2 γ . While the DM phase space is not thermal, these moments can be thought of as relating to the DM effective temperature, which will have ramifications for the subsequent cosmology. As the DM mass rises, the effective temperature increases because e + e − annihilations become more important than plasmon decays and have a comparatively fatter high-p χ tail. At even larger masses where m χ is comparable to m e , that high-p χ tail is suppressed because the DM mass becomes relevant to the kinematics of annihilation, causing the effective temperature to drop. in the cross section). The dependence on the Debye mass m D comes from the photon propagator for electric scattering in a medium [89]. The usual t-channel divergence is thus regulated in the forward-scattering limit by the Debye angle, defined as θ D ≡ m D /p CM . Once the plasma has become non-relativistic with T γ m e , the Debye mass is given by in natural units, assuming Ω b h 2 = 0.022 [90] and that the ionization fraction is unity. The momentum transfer cross section is defined for DM self-scattering in Eq. (3) and the analogous definition applies for scattering between DM and SM particles. For relativistic DM, we find that in the limit of the Debye angle θ D 1 Since m b m χ and the baryons are non-relativistic, the DM momentum in the CM frame can be approximated by the DM momentum in the comoving frame, p χ . As illustrated in Fig. 6, the typical DM momentum is comparable to the SM photon temperature, with both quantities redshifting after freeze-in. Therefore, we can estimate the momentum transfer rate per DM particle and per Hubble time as where n p ≈ 1.5 × 10 −10 T 3 γ and p χ ≈ 0.4 p γ ≈ T γ . For T γ in the keV-MeV range and Q < 10 −10 for freeze-in, this rate is tiny and thus scattering in this regime has a negligible effect on the DM phase space.
For scattering of non-relativistic DM and charged SM particles, the differential cross section is instead given by where µ χb is the reduced mass of the DM and baryon, µ χb = m χ m b /(m χ + m b ), v is the relative velocity between DM and SM particles, and p CM = µ χb v. The momentum transfer cross section is where again we take the θ D 1 limit. Note that the Coulomb logarithm appearing here differs from the one that appears in the often-quoted Ref. [78]; however, that reference did not include the Debye mass in the photon propagator, as discussed in Appendix D. Compared to the Coulomb logarithm in Ref. [78], our treatment of the Debye mass results in a factor of 2.5 − 3 smaller momentum transfer rate at recombination; this will translate to a weaker CMB bound on generic millicharged DM than has been reported previously [70][71][72][73][74], which we explore in more detail in our companion paper [53].
Given the velocity scaling in Eq. (37), momentum transfer is most important at late times. For freeze-in couplings, there may be a substantial effect at the recombination epoch. In particular, momentum transfer during this epoch leads to a drag force between the DM and baryon fluids, which can affect CMB anisotropies [70,[74][75][76]. The CMB bounds require that the momentum transfer rate is slow compared to the rate of Hubble expansion at z ≈ 1100, thus limiting the possible effect on the DM phase space. We calculate the bounds in detail in the companion paper [53], properly accounting for the velocity distribution for freeze-in DM with the updated Coulomb logarithm.
In addition to DM-baryon scattering as discussed above, DM-photon scattering is possible. However, these processes do not have the low-velocity v −4 enhancement in the rate and the cross section scales as Q 4 , so the effects are negligible. In the model with a dark photon A , scattering processes such as e − + γ → e − + A are also possible and scale only as kinetic mixing squared κ 2 . However, these processes are still negligible compared to DM-baryon scattering since they lack the low-v enhancement and have an additional large suppression due to the in-medium kinetic mixing effects, as discussed in Section II C. Processes like χ + γ → χ + A scale as Q 2 g 2 χ ; these also lack the v −4 enhancement and any enhancement (relative to DM-baryon scattering) from the large photon-to-baryon ratio is more than compensated by the factor of g 2 χ , even at the largest values of g χ that saturate SIDM bounds.

D. Effect of DM-DM scattering
In the absence of a dark photon, DM self scattering is proportional to Q 4 , rendering it entirely negligible. However, self-interactions of the DM can effectively thermalize the phase space distribution in the model with a dark photon. The rate for dark photon mediated DM scattering is proportional to g 4 χ , and thus may be important if g χ is sufficiently large compared to κ. Similar to DM-baryon scattering, the cross section scales as 1/v 4 and so these effects are most important at later times when the DM is cooler. Sufficient levels of self-scattering will convert a free-streaming phase space distribution into a Maxwell-Boltzmann or Gaussian velocity distribution. In the nonrelativistic limit, the quantity a(t) 2 p 2 χ will remain the same after this process (by conservation of comoving energy), although other moments of the phase space differ.
To determine when self-scattering becomes important, we estimate the redshift z therm when the momentum transfer rate per DM particle and per Hubble time is order unity: n χ σ T, χχ v H(z therm ) = 1 (38) where v is the relative velocity between DM particles and σ T, χχ is the self-scattering momentum transfer cross section given in Eq. (3), with the dark photon mass regulating the forward scattering instead of the Debye mass that is present for DM-baryon scattering. Using the ratio of the average DM momentum to the photon momentum in Fig. 6, we approximate the relative velocity as v ≈ p χ /m χ ≈ T γ (z)/m χ . In this estimate, we have assumed that DM is non-relativistic at the time self-interactions become important. The self scattering randomizes the DM velocities while preserving the average kinetic energy where p χ is physical momentum and the average momentum-squared is given in Fig. 6. After self-scattering becomes significant, the DM phase space is described by a thermal Maxwell-Boltzmann distribution, where n DM (z) is the DM number density. Fig. 7 shows the redshift of thermalization for two representative choices of κ (thus fixing g χ to yield the observed relic abundance), where we see the assumption of non-relativistic DM is a reasonably good approximation in our estimates. Since the phase space calculations here will be an input to determining CMB constraints on freeze-in DM, we compare z therm with the redshift of recombination z ≈ 1100. For constraints from structure formation, a range of redshifts will be relevant. We also show some fiducial limits from SIDM, which give upper bounds on g χ . Fig. 7 illustrates that the DM phase space at the time of recombination depends sensitively on the model parameters and on the robustness of SIDM limits in different astrophysical systems. For the largest values of g χ consistent with the weaker assumed SIDM bounds, the DM phase space is described by a Maxwell-Boltzmann distribution at the time of recombination for all the DM masses we consider. However, for κ = 10 −3 (which is consistent with bounds on ultralight dark photons), g χ is small enough that DM self-interactions are not important at recombination and the phase space is described by the results of Sections IV A-IV B. The comparison of the free-streaming and thermalized phase space can be seen in Fig. 5. The freeze-in relic abundance is determined by Q = g χ κ/e and we show z therm assuming two values of κ (where g χ is fixed to obtain the DM relic abundance). The epoch when DM self-thermalization becomes relevant is highly sensitive to the choice of couplings, which can yield different results for CMB observables depending on whether thermalization occurs before recombination. Note that DM halo formation is neglected in this estimate. Also shown are bounds on DM self-thermalization which come from the SIDM limits on g χ in Eq. (4). For illustration, we assume σ T, χχ 1 cm 2 /g for scattering via an ultralight mediator and show both v ∼ 10 −3 and v ∼ 10 −4 , speeds relevant to a halo the size of the Milky Way and to a dwarf galaxy. In this figure we have taken m A = 10 −14 eV, which is sufficiently light that the constraints on the kinetic mixing parameter κ are rather weak.

V. RESULTS AND DISCUSSION
In this paper, we have shown that DM freeze-in through a light vector mediator is substantially affected by plasmon decay, which constitutes a new production channel. This is an efficient way of producing sub-MeV DM and is dominant over SM fermion annihilation for masses below a few hundred keV. To account for this extra production channel, the couplings between the DM and the SM must be reduced in order to obtain the observed relic abundance of DM. For the lightest DM masses that are accessible to low-threshold direct detection experiments, the predicted cross section is lowered by roughly an order of magnitude. Updated predictions for freeze-in through a light vector mediator are shown in Fig. 8.
The presence of this channel also affects the DM phase space. In the absence of plasmon decays, the DM is never technically thermal but it acquires a distribution that appears thermal by inheriting the electron phase space distribution at the time of production. At early times f χ, e + e − (p χ ) ∼ e −pχ/T χ, e + e − , where T χ, e + e − is an effective DM temperature inherited from the electrons; at late times, this exponential distribution persists because the DM does not thermalize

Super-CDMS G2
Fre eze -in FIG. 8. Summary plot including early-universe plasma effects for the parameter space of sub-MeV freeze-in DM. The correct DM relic abundance is obtained for couplings on the freeze-in line. We show constraints coming from emission of DM pairs in white dwarf, horizontal branch and red giant stars [46], while bounds from emission of DM pairs in supernovae apply for Q 10 −9 [79]. Dotted lines are projected sensitivities of proposed direct detection experiments as in Fig. 3.
to give the Maxwell-Boltzmann distribution that would be expected for non-relativistic matter in equilibrium. On the other hand, the plasmon decay channel yields a DM phase space distribution that never appears thermal, which can be attributed to the running of the plasmon mass with temperature and the fact that plasmon decays occur dominantly as the plasmon wavenumber k → 0. For DM masses where plasmon decays are the dominant production mode, the phase space is peaked at low momentum and has a long tail; for DM masses where contributions from both channels are important, the phase space distribution is bimodal.
Though the DM is born with a highly non-thermal distribution, it may be possible for the DM to thermalize with itself under the right circumstances. For DM that is only charged under the SM U (1) EM with millicharge Q, the thermalization rate is suppressed by a factor of Q 4 where the requisite Q to produce the DM relic abundance is Q ∼ O 10 −11 . If the DM is also charged under a dark U (1) gauge group that kinetically mixes with the SM U (1) EM (with mixing parameter κ), it may be possible for DM self-scattering to thermalize the DM phase space distribution. In this case, Q = κg χ /e (where κ can take on a wide range of values) and DM self-scattering via the dark photon scales as g 4 χ , meaning that with the appropriate choice of κ and g χ it is possible to efficiently selfscatter while still producing the observed relic abundance. The coupling g χ cannot be arbitrarily large due to observational limits on SIDM in astrophysical systems; however, there is a range of g χ where self-scattering thermalizes the DM before recombination and where the SIDM bounds are simultaneously satisfied. Energy is conserved within the DM fluid, so for non-relativistic DM p 2 χ will be conserved and the resulting distribution has a well-defined notion of temperature.
Although the freeze-in DM phase space distribution may not be thermal, it is still informative to take moments of the distribution. When comparing the first and second moments of f χ (p χ ) to the equivalent quantities for the SM photon bath, we find that the typical DM momentum is similar to the typical photon momentum, p χ ≈ (0.4 − 0.7) × p γ depending on the DM mass. In other words, the DM is born considerably warmer than what is typically assumed for cold DM initial conditions. This will have ramifications for cosmology in two key ways: • Freeze-in DM will behave like warm DM, leading to suppression of the matter power spectrum below some physical scale roughly corresponding to the free-streaming length. This effect is not already captured by existing limits on warm DM, where different DM phase space distributions are assumed. To understand this suppression quantitatively, a Boltzmann code is necessary that accounts for the potentially nonthermal phase space from freeze-in. Having understood this, it will be possible to constrain DM freeze-in via a light vector mediator using probes of the matter power spectrum and the halo mass function.
• Existing CMB limits on DM with an effective millicharge do not straightforwardly apply to the case of freeze-in. These limits stem from a DM-baryon drag; because the drag is highly sensitive to the relative DM-baryon velocity (the cross section scales like ∼ v −4 ), modifications to the DM phase space can substantially alter the size of the effect. Existing limits have made the assumption of cold dark matter, and the larger DM velocities for freezein will lead to reduced drag force. Taking into account the updated Debye logarithm (which may weaken existing limits by a factor of ∼ 2 − 3), the limit on freeze-in will be further reduced compared to previously reported results.
Both of these effects will be thoroughly explored in our companion paper [53], which will place restrictions on the range of masses where DM freeze-in via a light mediator is observationally viable. Throughout this work, we take the properties of the SM thermal bath to be given by their equilibrium values at zero chemical potential. The photons and neutrinos are relativistic gases with energy and entropy densities Here we distinguish between the neutrino and SM bath temperatures T and T ν ; in this work we assume that the neutrinos kinetically decouple at a temperature that is higher than relevant for sub-MeV freeze-in and that their temperature evolves adiabatically T ν ∼ 1/a during this epoch, which is a good approximation at the percent level. We also ignore the negligible neutrino masses. Meanwhile, the electrons are transitioning from being relativistic to being non-relativistic, so we use the unapproximated expressions for the energy and entropy density, Throughout the evolution of the SM bath, we require conservation of entropy. Since we are assuming adiabatic evolution of the neutrino temperature, its entropy s ν (T ν )a(T ) 3 is constant by definition. The remaining constraint equation on the temperature evolution is then (s γ (T ) + s e (T )) a 3 = const., which yields a smooth temperature evolution T (a), as shown in Fig. 9. After the electrons have fully left the bath, we recover the usual result T ν = (4/11) 1/3 T . We can then use this temperature evolution to evolve the Hubble parameter smoothly through the transition as the electrons leave the thermal bath, H 2 (a) = ρ e (T (a)) + ρ γ (T (a)) + ρ ν (T ν (a)) 3M 2 Pl (A4) with M Pl the reduced Planck mass. Both the temperature and Hubble evolution feed into the calculations of the DM relic abundance and phase space in the main body of the text.
The poles in the propagator yield the renormalized longitudinal and transverse dispersion relations for on-shell plasmons, ω (k) 2 = ω (k) 2 k 2 Π (ω (k), k) ω t (k) 2 = k 2 + Π t (ω t (k), k), while the residues of the poles are identified as a combination of dressed polarization four-vectors, µ (k)˜ ν (k) * , for the appropriate polarization. The dressed polarization vectors are given bỹ Given the approximations for Π and Π t and the dispersion relations, the residue functions can be written as .

(B11)
Appendix C: Plasmon decays through a dark photon In this Appendix, we show that plasmon decays in the millicharge basis (Eq. (6)) are identical to decays in the basis where the dark photon has a coupling eκJ µ EM A µ . In a thermal plasma, this coupling generates an in-medium mixing term in the Lagrangian given by κA µ Π µν A ν where Π µν is the electromagnetic polarization tensor. The matrix element in the dark photon basis is then given by iM = iκg χ˜ µ (k)Π µν ω, k D να A ω, k ū(p χ )γ α v(pχ) ≡ iκg χ˜ µ (k)ū(p χ )γ α v(pχ) Γ α µ , where D να A is the dark photon propagator. Taking the m A = 0 limit and working in Coulomb gauge, the propagator is given by which gives the same result as the vertex obtained in the millicharge basis. millicharge, with θ min = 2Qα/(3T λ D ). The corresponding minimum momentum transfer in that case would be | q| 2 = 4Q 2 α 2 p 2 CM m 2 D /(9T 2 ). For freeze-in where p CM ≈ T and Q < 10 −10 , we see that | q| 2 m 2 D and so we expect that the Yukawa-like form of the effective potential leads to a strong screening effect for modes of such large spatial size. In other words, the requirement of Ref. [78] may not be restrictive enough because DM-baryon scattering is suppressed by factors of Q relative to the strong collective effects in the plasma that give rise to the Debye mass. Because forward scattering is so peaked, the resulting transfer cross section is highly sensitive to the limits of integration and their procedure yields a transfer cross section that is a factor of ∼ 2 − 3 larger than the one obtained with the procedure of Eq. (D3). As a result, CMB limits on millicharged DM that use this result may be too strong.