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

TL;DR

This paper refines the calculation of dark matter relic abundance produced via freeze-in through a light vector mediator, highlighting the importance of plasma effects like plasmon decay, which significantly influence detection prospects.

Contribution

It introduces an additional production channel via plasmon decay and updates relic abundance calculations, emphasizing plasma effects in freeze-in dark matter models.

Findings

Plasmon decay is a dominant production channel for sub-MeV dark matter.

Including plasma effects reduces the expected signals in dark matter searches.

Dark matter exhibits a non-thermal phase space distribution affecting cosmology.

Abstract

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.

Equations

$$Y_{\chi }=\frac{n_{\chi }}{s}\sim \frac{n_{\mathrm{SM}}^2⟨\sigma v⟩}{sH}\sim 2\times 10^{-4}\frac{g_{\chi }^2{g}_{\mathrm{SM}}^2{M}_{\mathrm{Pl}}}{T},$$

$$\mathcal{L}\supset$$

$$+e{J}_{\mathrm{EM}}^{\mu }{A}_{\mu }+g_{\chi }\bar{\chi }{\gamma }^{\mu }\chi A_{\mu }^{\prime }+\bar{\chi }(i\partial -m_{\chi })\chi ,$$

$${\sigma }_{T,\chi \chi }$$

$$g_{\chi }\lesssim 4\times 10^{-5}\times \left(\frac{v}{10^{-3}}\right)\times {\left(\frac{m_{\chi }}{1\text{MeV}}\right)}^{3/4}\times {\left(\frac{10}{\ln \left(m_{\chi }^2{v}^2/m_{A^{\prime }}^2\right)}\right)}^{1/4}.$$

$${\mathcal{L}}_{\mathrm{IM}}\supset$$

$$+J_{\mathrm{EM}}^{\mu }\left(e{\tilde{A}}_{\mu }+\frac{e\kappa m_{A^{\prime }}^2}{m_{A^{\prime }}^2-m_A^2}{\widetilde{A^{\prime }}}_{\mu }\right)+g_{\chi }\bar{\chi }{\gamma }^{\mu }\chi \left({\widetilde{A^{\prime }}}_{\mu }-\frac{\kappa m_A^2}{m_{A^{\prime }}^2-m_A^2}{\tilde{A}}_{\mu }\right).$$

$${\mathcal{L}}_{\mathrm{IM}}\supset J_{\mathrm{EM}}^{\mu }\left(e{\tilde{A}}_{\mu }\right)+g_{\chi }\bar{\chi }{\gamma }^{\mu }\chi \left({\widetilde{A^{\prime }}}_{\mu }+\kappa {\tilde{A}}_{\mu }\right),$$

$$Y\equiv n_{\chi }/s=4.35\times 10^{-7}\times \left(\frac{1\mathrm{MeV}}{m_{\chi }}\right).$$

$$\frac{d{n}_{\text{DM}}}{da}+\frac{3n_{\text{DM}}}{a}=\frac{2}{aH}\left({⟨\sigma v⟩}_{e^{+}{e}^{-}\to \chi \bar{\chi }}{n}_e^2+{⟨\Gamma ⟩}_{{\gamma }^{\ast }\to \chi \bar{\chi }}{n}_{{\gamma }^{\ast }}\right).$$

$$\displaystyle\sum_{\mathrm{d.o.f.}}\left|\mathcal{M}\right|^{2}_{e^{+}e^{-}\leftrightarrow\chi\bar{\chi}}=\frac{32Q^{2}e^{4}}{(p_{e^{+}}+p_{e^{-}})^{4}}\Big{(}$$

$$\displaystyle+m_{e}^{2}(p_{\chi}\cdot p_{\bar{\chi}})+m_{\chi}^{2}(p_{e^{+}}\cdot p_{e^{-}})+2m_{e}^{2}m_{\chi}^{2}\Big{)},$$

$${⟨\sigma v⟩}_{e^{+}{e}^{-}\to \chi \bar{\chi }}{n}_e^2=\int \frac{\text{mathchar}22d^3{p}_{e^{+}}}{2E_{e^{+}}}\frac{\text{mathchar}22d^3{p}_{e^{-}}}{2E_{e^{-}}}$$

$$\times \sum _{d.o.f.}{\left|\mathcal{M}\right|}_{e^{+}{e}^{-}\to \chi \bar{\chi }}^2(2\pi {)}^4{\delta }^{(4)}(p_{e^{+}}+p_{e^{-}}-p_{\chi }-p_{\bar{\chi }})$$

$$\int \frac{d^3{q}_{12}d{s}_{12}}{2E_{12}}{\delta }^{(4)}(q_{12}-p_1-p_2)=1,$$

$${⟨\sigma v⟩}_{e^{+}{e}^{-}\to \chi \bar{\chi }}{n}_e^2$$

$$\times \sum _{d.o.f.}{\left|\mathcal{M}\right|}_{e^{+}{e}^{-}\to \chi \bar{\chi }}^2(2\pi {)}^4{\delta }^{(4)}(p_{e^{+}}+p_{e^{-}}-p_{\chi }-p_{\bar{\chi }}){\delta }^{(4)}(q_{\chi \bar{\chi }}-p_{\chi }-p_{\bar{\chi }}).$$

$${\Phi }_{\chi \bar{\chi }}(s_{\chi \bar{\chi }}){\left|\mathcal{M}\right|}_{\text{CM}}^2(s_{\chi \bar{\chi }})\equiv \int \frac{\text{mathchar}22d^3{p}_{\chi }}{2E_{\chi }}\int \frac{\text{mathchar}22d^3{p}_{\bar{\chi }}}{2E_{\bar{\chi }}}(2\pi {)}^4{\delta }^{(4)}(q_{\chi \bar{\chi }}-p_{\chi }-p_{\bar{\chi }})\sum _{d.o.f.}{\left|\mathcal{M}\right|}_{e^{+}{e}^{-}\to \chi \bar{\chi }}^2$$

$$=\frac{Q^2{e}^4}{2\pi s_{\chi \bar{\chi }}^2}\sqrt{1-\frac{4m_{\chi }^2}{s_{\chi \bar{\chi }}}}\left(s_{\chi \bar{\chi }}^2+\frac{1}{3}(s_{\chi \bar{\chi }}-4m_e^2)(s_{\chi \bar{\chi }}-4m_{\chi }^2)+4s_{\chi \bar{\chi }}(m_{\chi }^2+m_e^2)\right),$$

$${⟨\sigma v⟩}_{e^{+}{e}^{-}\to \chi \bar{\chi }}{n}_e^2=\int \frac{d^3{q}_{\chi \bar{\chi }}d{s}_{\chi \bar{\chi }}}{2E_{\chi \bar{\chi }}}{e}^{-E_{\chi \bar{\chi }}/T}$$

$$\times$$

$${\Phi }_{e^{+}{e}^{-}}(s_{\chi \bar{\chi }})\equiv \frac{1}{8\pi }\sqrt{1-\frac{4m_e^2}{s_{\chi \bar{\chi }}}},$$

$${⟨\sigma v⟩}_{e^{+}{e}^{-}\to \chi \bar{\chi }}{n}_e^2$$

$${⟨\sigma v⟩}_{e^{+}{e}^{-}\to \chi \bar{\chi }}{n}_e^2=\frac{T}{(2\pi {)}^3}\int ds\sqrt{s}{\Phi }_{e^{+}{e}^{-}}(s){\Phi }_{\chi \bar{\chi }}(s){\left|\mathcal{M}\right|}_{\text{CM}}^2(s)K_1(\sqrt{s}/T)$$

$$i{\mathcal{M}}_{{\gamma }^{\ast }\to \chi \bar{\chi }}=iQe{\widetilde{ϵ}}_{\mu }(k)\bar{u}(p_{\chi }){\gamma }^{\mu }v(p_{\bar{\chi }}),$$

$$\sum _{d.o.f.}{\left|\mathcal{M}\right|}_{{\gamma }^{\ast }\to \chi \bar{\chi }}^2=4Q^2{e}^2\times \{\begin{array}{ll}2Z_t(k)(p_{\chi }^2{\sin }^2\theta +{\omega }_t(k)E_{\chi }-k{p}_{\chi }\cos \theta )& ++\&--\\ Z_{\ell }(k)\frac{{\omega }_{\ell }(k{)}^2}{k^2}({\omega }_{\ell }(k)E_{\chi }-2E_{\chi }^2+k{p}_{\chi }\cos \theta )& \text{LL},\end{array}$$

$${⟨\Gamma ⟩}_{{\gamma }^{\ast }\to \chi \bar{\chi }}{n}_{{\gamma }^{\ast }}=\int \frac{\text{mathchar}22d^{3}k}{2\omega (k)}\frac{\text{mathchar}22d^3{p}_{\chi }}{2E_{\chi }}\frac{\text{mathchar}22d^3{p}_{\bar{\chi }}}{2E_{\bar{\chi }}}f\left(\omega (k)\right)(2\pi {)}^4{\delta }^{(4)}\left(K-p_{\chi }-p_{\bar{\chi }}\right)\sum _{d.o.f.}{\left|\mathcal{M}\right|}_{{\gamma }^{\ast }\to \chi \bar{\chi }}^2,$$

$${⟨\Gamma ⟩}_{{\gamma }_{\ell }^{\ast }\to \chi \bar{\chi }}{n}_{{\gamma }_{\ell }^{\ast }}=\frac{Q^2{e}^2}{(2\pi {)}^3}\int k^{2}dk\frac{Z_{\ell }(k){\omega }_{\ell }(k)(m_{\ell }(k{)}^2+2m_{\chi }^2)\sqrt{m_{\ell }(k{)}^2(m_{\ell }(k{)}^2-4m_{\chi }^2)}}{3m_{\ell }(k{)}^4\left(e^{{\omega }_{\ell }(k)/T}-1\right)}$$

$${⟨\Gamma ⟩}_{{\gamma }_t^{\ast }\to \chi \bar{\chi }}{n}_{{\gamma }_t^{\ast }}=\frac{4Q^2{e}^2}{(2\pi {)}^3}\int k^{2}dk\frac{Z_t(k)(m_t(k{)}^2-m_{\chi }^2)\sqrt{m_t(k{)}^2(m_t(k{)}^2-4m_{\chi }^2)}}{3{\omega }_t(k)m_t(k{)}^2\left(e^{{\omega }_t(k)/T}-1\right)},$$

$${\sigma }_e=\frac{16\pi Q^2{\alpha }^2{\mu }_{\chi e}^2}{(\alpha m_e{)}^4}.$$

$$\frac{\partial f_{\chi }}{\partial t}-H\frac{p_{\chi }^2}{E_{\chi }}\frac{\partial f_{\chi }}{\partial E_{\chi }}=\frac{C(p_{\chi },t)}{E_{\chi }},$$

$$f_{\chi }(p_{\chi },t)={\int }_{t_i}^{t}d{t}^{\prime }\frac{C\left(\frac{a(t)}{a(t^{\prime })}{p}_{\chi },t^{\prime }\right)}{\sqrt{\frac{a(t{)}^2}{a(t^{\prime }{)}^2}{p}_{\chi }^2+m_{\chi }^2}}={\int }_{a_i}^{a(t)}\frac{d{a}^{\prime }}{a^{\prime }H(a^{\prime })}\frac{C\left(\frac{a(t)}{a^{\prime }}{p}_{\chi },a^{\prime }\right)}{\sqrt{\frac{a(t{)}^2}{a^{\prime 2}}{p}_{\chi }^2+m_{\chi }^2}}.$$

$$C(p_{\chi },t{)}_{e^{+}{e}^{-}\to \chi \bar{\chi }}=\frac{1}{2(2\pi {)}^3}\int \frac{d^3{q}_{e^{+}{e}^{-}}d{s}_{e^{+}{e}^{-}}}{2E_{\bar{\chi }}2E_{e^{+}{e}^{-}}}$$

$$\times {\Phi }_{e^{+}{e}^{-}}(s_{e^{+}{e}^{-}}){\left|\mathcal{M}\right|}_{\text{CM}}^2(s_{e^{+}{e}^{-}}),$$

Full text

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, $\langle\sigma v\rangle\sim 10^{-26}$ cm3/s. This weak-scale 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_{\chi}\gtrsim 1$ MeV Boehm et al. (2013); Nollett and Steigman (2014, 2015). 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_{\chi}\sim$ GeV$-$TeV mass range. More recently, there has been a growing interest in $m_{\chi}\sim$ MeV$-$GeV thermal candidates where interactions with the SM or within a hidden sector deplete the DM density to the observed value Boehm and Fayet (2004); Pospelov et al. (2008); Feng and Kumar (2008); Hochberg et al. (2014, 2015, 2018a); Choi et al. (2017); D’Agnolo et al. (2018, 2017); D’Agnolo and Ruderman (2015); Pappadopulo et al. (2016); Cline et al. (2017); Kopp et al. (2016).

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 Asaka et al. (2006, 2007); Gopalakrishna et al. (2006); Page (2007); Hall et al. (2010); Bernal et al. (2017). 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 $\langle\sigma v\rangle\sim g_{\chi}^{2}g_{\rm SM}^{2}/(4\pi T)^{2}$ where $g_{\chi}$ is the mediator-DM coupling, $g_{\rm 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.111We assume the minimal scenario where the dark sector is not populated in abundance at reheating.

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 $\sigma\propto 1/v^{4}$ for relative velocity $v$, implying that the kinematics of the Milky Way (where $v\sim 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 Essig et al. (2012a, b, 2016); Hochberg et al. (2017); Derenzo et al. (2017); Hochberg et al. (2018b); Knapen et al. (2017a); Griffin et al. (2018); Schutz and Zurek (2016); Knapen et al. (2017b); Hochberg et al. (2016a, b) (see also Ref. Battaglieri et al. (2017) for a recent review). Indeed, recent experimental results by SENSEI Crisler et al. (2018); Abramoff et al. (2019), SuperCDMS Agnese et al. (2018), and DarkSide Agnes et al. (2018) are demonstrating significant progress towards achieving the sensitivity needed in the MeV-GeV mass range. It was also shown recently that Xenon1T Aprile et al. (2018) is for the first time constraining freeze-in in the GeV-TeV mass range Hambye et al. (2018).

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. Knapen et al. (2017c); Green and Rajendran (2017) 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 electron-positron annihilations via a light vector mediator has been studied in the past Chu et al. (2012); Essig et al. (2012a), 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. Davidson et al. (2000a); Vogel and Redondo (2014) and the effect on freeze-in via a heavy mediator was recently considered in Ref. An et al. (2018) 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 Bowman et al. (2018); Barkana (2018); Barkana et al. (2018); Berlin et al. (2018); Muñoz and Loeb (2018). 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 Dvorkin et al. , 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_{\chi}\sim$ keV$-$MeV. In particular, we find that existing cosmological constraints restrict $m_{\chi}\gtrsim$ tens of keV for freeze-in via a light mediator, and it will be possible to probe even higher masses with planned experiments.

II Models for sub-MeV freeze-in

II.1 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 Davidson et al. (2000a); Dubovsky et al. (2004) 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.222Other models that have been considered in the past require giving neutrinos small charges as well Foot et al. (1993), which we do not consider further due to strong experimental bounds on neutrino charge Chen et al. (2014).

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\sim m_{\chi}v\sim 10^{-3}m_{\chi}$, while for electron recoils the typical electron momentum in the target material is most relevant $q\sim\alpha m_{e}\approx 4$ keV, where $m_{e}$ is the electron mass and $\alpha$ is the electromagnetic fine structure constant. Thus for sub-MeV DM, the experimentally relevant mediators have masses below $\mathcal{O}(1)$ keV.

Assuming an annihilation cross section of SM fermions into DM with the form $\langle\sigma v\rangle\sim g_{\chi}^{2}g_{\rm SM}^{2}/(4\pi T)^{2}$, the relic abundance can be estimated as

[TABLE]

where $M_{\rm Pl}=1/\sqrt{8\pi G}$ is the reduced Planck mass and we assumed $T\sim$ MeV. Then for $m_{\chi}\sim$ MeV, we find that $g_{\chi}g_{\rm SM}\simeq 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_{\chi}g_{\rm SM}\sim 10^{-12}$, $g_{\rm SM}$ must be greater than $10^{-12}$ if we require the dark sector to be perturbative (i.e. $g_{\chi}~{}\lesssim~{}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. Knapen et al. (2017c); Green and Rajendran (2017) 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_{\phi ee}\phi\bar{e}e$ comes from helium ignition in red giants, with $g_{\phi ee}\lesssim 7\times 10^{-16}$ for sub-keV masses Hardy and Lasenby (2017). For a sub-keV pseudoscalar, observations of white dwarfs give typical constraints of $g_{aee}\lesssim 2\times 10^{-13}$ Raffelt (2008); Viaux et al. (2013); Miller Bertolami et al. (2014). 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 Carlson and Salati (1989); Raffelt and Starkman (1989).

• •

Scalars and pseudoscalars coupled to nucleons $-$ Similar to the case of mediators coupling to electrons, red giants constrain $g_{\phi nn}\lesssim 10^{-12}$ for a scalar Hardy and Lasenby (2017) and $g_{ann}\lesssim{\rm few}\times 10^{-10}$ for a pseudoscalar Raffelt (2008); Giannotti et al. (2017). 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_{\rm SM}$ with $T\sim 1$ MeV is much too low and freeze-in would have to occur with a larger value of $g_{\rm 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. Krnjaic (2017) 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_{\rm SM}$ decrease linearly with the mediator mass for masses below $\sim 100$ eV An et al. (2013a, b) 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.3. From the collected bounds on dark photons from Refs. Jaeckel (2012), a dark photon can have $g_{\rm SM}>10^{-12}$ when its mass is well below $1$ eV. At even lower masses, the coupling could be $\sim 10^{-3}$ for masses below $\lesssim 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. Murata and Tanaka (2015); Adelberger et al. (2003)), 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.3), which is relevant to observations of the CMB Dvorkin et al. (2014); Xu et al. (2018); Slatyer and Wu (2018); Kovetz et al. (2018); Boddy et al. (2018); Dubovsky and Gorbunov (2001); Boddy and Gluscevic (2018) and the cosmological 21 cm global signal Bowman et al. (2018); Barkana (2018); Barkana et al. (2018); Berlin et al. (2018); Muñoz and Loeb (2018). 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.

II.2 DM with photon-mediated interactions

If the DM is a Dirac fermion $\chi$ 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}\equiv 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\ll m_{Z}$. This gives the simplest model of millicharged DM. It is difficult to incorporate such matter content into a Grand Unified Theory (GUT) Okun et al. (1984); 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. McDermott et al. (2011), 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_{\chi}=eQ$ and $g_{\rm 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. Davidson et al. (2000a); Chang et al. (2018) and require $Q\lesssim 10^{-9}$ for $m_{\chi}$ up to a few MeV, which does not impact freeze-in. However, there are constraints for $m_{\chi}$ below $\mathcal{O}(10)$ keV from emission in white dwarfs, horizontal branch stars, and red giants (see Appendix of Ref. Vogel and Redondo (2014)). Note that the range of $m_{\chi}$ 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.3, 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.

II.3 DM with dark photon interactions

We next consider Dirac fermion DM coupled to a kinetically mixed dark photon $A^{\prime}$, with the vacuum Lagrangian given by

[TABLE]

where $A$ is the SM photon, $\kappa$ is the kinetic mixing parameter and $\chi$ is Dirac fermion DM. For the purposes of this discussion, we consider Abelian kinetic mixing, noting that non-Abelian kinetic mixing is also possible Barello et al. (2016); Argüelles et al. (2017). The mixing parameter $\kappa$ 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^{\prime}$ Dienes et al. (1997) or from certain compactifications of type IIB strings Abel et al. (2008); Goodsell et al. (2009). 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_{\chi}g_{\text{SM}}=g_{\chi}\kappa e$.

As discussed in Section II.1, the dark photon mass must satisfy $m_{A^{\prime}}\lesssim 1$ eV in order to give a sufficient coupling for freeze-in while also evading existing bounds on stellar energy loss Jaeckel (2012). However, the requirements are even more stringent because unlike the model presented in Section II.2 there could be large $A^{\prime}$-mediated DM self interaction. For $m_{A^{\prime}}<$ 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_{\chi}^{4}$. Furthermore, as mentioned before, the freeze-in relic abundance is determined by the product $g_{\chi}\kappa e$, meaning that large $g_{\chi}$ can be compensated by reducing $\kappa$ 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 Feng et al. (2010)

[TABLE]

where $\theta_{\text{CM}}$ is the scattering angle in the center-of-mass (CM) frame, $\sigma_{\chi\chi}$ is the self-interaction cross section, and $\alpha_{\chi}$ is the dark equivalent of the electromagnetic fine structure constant, $\alpha_{\chi}\equiv g_{\chi}^{2}/4\pi$. Typical bounds on SIDM require $\sigma_{\chi\chi}/m_{\chi}<1-10$ cm2/g for systems ranging from dwarf galaxies where $v\sim 10^{-4}$ to merging clusters where $v\sim 10^{-2}$ (for a recent review, see Ref. Tulin and Yu (2018)). While few simulation-based studies of self-interactions have been done in the ultralight mediator limit (see for instance Ref. Kummer et al. (2019)), we can estimate the expected bound. Taking the more restrictive limit of $\sigma_{\chi\chi}/m_{\chi}\sim$1 cm${}^{2}/$g, the bound is

[TABLE]

Since $\kappa eg_{\chi}\gtrsim 10^{-12}$ is needed for sub-MeV freeze-in, the SIDM bounds imply that the kinetic mixing is $\kappa\gtrsim 10^{-7}$ for MeV-scale DM. For sub-eV dark photons, such large kinetic mixing is only possible when $m_{A^{\prime}}\lesssim 10^{-10}$ eV Jaeckel (2012). For even lighter DM, $g_{\chi}$ is even more restricted so $\kappa\gtrsim 10^{-5}$ is required for freeze-in, which is possible when $m_{A^{\prime}}\lesssim 10^{-14}$ eV. Therefore, we are required to consider an “ultralight” dark photon Knapen et al. (2017c). Note that black hole superradiance constrains dark photons being present in the mass spectrum (in the small-coupling limit) between $\sim 10^{-14}-10^{-11}$ eV and preliminarily between $\sim 10^{-19}-10^{-17}$ eV Baryakhtar et al. (2017).

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^{\prime}}$ is much lower than the effective in-medium photon mass $m_{A}$ in the early universe. Then, following Appendix D of Ref. Knapen et al. (2017c), the vacuum Lagrangian in Eq. (2) is modified with an additional term $m_{A}^{2}A^{\mu}A_{\mu}/2$.333For simplicity we consider a constant $m_{A}^{2}$ for the schematic purposes of this discussion, although the photon polarization tensor $\Pi^{\mu\nu}(\vec{q},\omega)$ (which gives rise to the in-medium effective mass) depends on the photon momentum $\vec{q}$, energy $\omega$, polarization, and thermal properties of the medium. For an on-shell mode with $\omega\sim|\vec{q}|$, $m_{A}^{2}$ would correspond to the plasma mass, as discussed in Section III.2. For scattering processes with a highly off-shell mode, $|\vec{q}|\gg\omega$, $m_{A}^{2}$ is given by the Debye mass Blaizot et al. (1995). Rotating away the mixing term in the presence of $m_{A}$ and $m_{A^{\prime}}$ and rewriting in terms of the mass eigenstates $\tilde{A}$ and $\tilde{A}^{\prime}$, the in-medium Lagrangian is given by

[TABLE]

From this, we see that when $m_{A}\gg m_{A^{\prime}}$, the interaction terms above reduce to

[TABLE]

meaning that DM has an effective millicharge parameter $Q=\kappa g_{\chi}/e$, and the interactions are identical to those for a massless dark photon. Note that this suppression of the $A^{\prime}$-SM coupling in the $m_{A^{\prime}}\ll m_{A}$ limit is the source of the in-medium (plasma mass) suppression of the stellar constraints on dark photons An et al. (2013a, b) discussed in Section II.1. Also note that this suppression 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_{\text{eff}}$.

In the exactly massless $A^{\prime}$ limit, we are free to perform a field redefinition on $A^{\prime}\to A^{\prime}+\kappa A$ in the vacuum Lagrangian, Eq. (2), which eliminates the kinetic mixing term and generates a DM interaction term $g_{\chi}\bar{\chi}\gamma^{\mu}\chi(A^{\prime}_{\mu}+\kappa A_{\mu})$, which is again identical to having a millicharge $Q=\kappa g_{\chi}/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^{\prime}$, 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 $f\bar{f}\to\chi\bar{\chi}$ that was previously calculated in Refs. Essig et al. (2012a); Chu et al. (2012). 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, $\gamma^{*}\to\chi\bar{\chi}$, which we calculate for the first time. Photon annihilation into DM $\gamma\gamma\to\chi\bar{\chi}$ 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 $\omega_{c}\equiv\Omega_{c}h^{2}=0.12$ Aghanim et al. (2018). 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

[TABLE]

For $m_{\chi}\gtrsim 1$ keV, the DM yield is much lower than the order unity yield for relativistic species, such that DM contributes negligibly to $N_{\text{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_{\text{eff}}$.444An exception for thermal, sub-MeV DM was pointed out in Ref. Berlin and Blinov (2018), where the DM thermalizes with the SM thermal bath after neutrino-photon decoupling, reducing the contribution to $N_{\text{eff}}$. Furthermore, in this model changes to $N_{\text{eff}}$ that occur after DM thermalization are compensated by decoupling at a later time.

The low DM occupation number also implies that it is possible to self-consistently ignore back-reactions 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 $\chi\bar{\chi}\to\gamma^{{}^{\prime}}\gamma^{{}^{\prime}}$ is negligible for the same reason. In what follows, we solve the $0^{\text{th}}$ moment of the Boltzmann equation ignoring back-reactions, noting that we have numerically checked that they are negligible. The relevant equation is then

[TABLE]

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 $\chi$ and $\bar{\chi}$ in the matter budget; assuming zero DM chemical potential, $n_{\text{DM}}=2n_{\chi}=2n_{\bar{\chi}}$, which accounts for the factor of two in Eq. (8).555This factor is related to the usual factor of $1/2$ that appears in the Boltzmann equation for Dirac fermions Gondolo and Gelmini (1991); Srednicki et al. (1988); 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_{\text{DM}}^{2}$ appearing in Eq. (8)) which accounts for the factor of four difference.

III.1 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^{\prime}}\rightarrow 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

[TABLE]

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=\kappa g_{\chi}/e$. The thermally averaged cross section appearing in Eq. (8) for this process is given by

[TABLE]

where $\mathchar 22\relax\mkern-12.0mud^{3}p\equiv d^{3}p/(2\pi)^{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\lesssim 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 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,

[TABLE]

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}=\sqrt{s_{12}+\vec{q}_{12}^{\,2}}$). Inserting such a factor into Eq. (10) gives

[TABLE]

The integral over $p_{\chi}$ and $p_{\bar{\chi}}$ does not depend on the frame of $q_{\chi\bar{\chi}}$, so the two-body phase space of $p_{\chi}$ and $p_{\bar{\chi}}$ can be evaluated in the CM frame of $q_{\chi\bar{\chi}}$. We define

[TABLE]

and insert this into the expression for the thermally averaged cross section

[TABLE]

Again, we can evaluate the integral over $p_{e^{+}}$ and $p_{e^{-}}$ in the center-of-mass frame. Defining

[TABLE]

the thermally averaged cross section becomes

[TABLE]

We can write this result in terms of the first order modified Bessel function of the second kind $K_{1}(z)=z\int_{1}^{\infty}du\,e^{-zu}\sqrt{u^{2}-1}$ with $u=\sqrt{1+q_{\chi\bar{\chi}}^{2}/s_{\chi\bar{\chi}}}\,$ :

[TABLE]

where we have dropped the subscript on the integration variable $s$. Note that $s$ is restricted to $s>4\max\left(m_{e}^{2},m_{\chi}^{2}\right)$. The procedure above provides an alternate derivation of the well-known results from Ref. Gondolo and Gelmini (1991), and we have validated this method here because we use similar techniques to derive the full collision term for annihilation in Section IV.1.

III.2 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 $\omega_{p}=eT/3\approx 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 Braaten and Segel (1993). 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

[TABLE]

where $\tilde{\epsilon}_{\mu}(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^{\prime}}\rightarrow 0$. In squaring and summing over polarizations, only the diagonal terms ($LL$, $++$, and $--$) contribute,

[TABLE]

where the photon four-momentum is given by K^{\mu}=\big{(}\omega(k),\vec{k}\big{)}^{\mu} with appropriate dispersion relations for transverse and longitudinal modes $\omega_{t}(k)$ and $\omega_{\ell}(k)$ (see Appendix B), the DM four-momentum is given by $\left(E_{\chi},\vec{p}_{\chi}\right)^{\mu}$, $\theta$ is the angle between $\vec{k}$ and $\vec{p}_{\chi}$, and $Z_{t}(k)$ and $Z_{\ell}(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

[TABLE]

and can be evaluated directly. Taking the plasmons to be Bose-Einstein distributed, the longitudinal and transverse contributions to this rate are

[TABLE]

[TABLE]

where the effective plasmon masses are $m_{\ell}(k)^{2}=\omega_{\ell}(k)^{2}-k^{2}$ for the longitudinal modes and $m_{t}(k)^{2}=\omega_{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.

III.3 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, freeze-in from $e^{+}e^{-}$ annihilation is always kinematically allowed and this process only ends when the electron number density becomes Boltzmann suppressed, namely $T\lesssim 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_{\gamma^{*}}\sim\omega_{p}=2m_{\chi}$. Since $\omega_{p}\approx 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_{\chi}$ = MeV. The change to the freeze-in benchmark for direct detection is shown in Fig. 3,

where the cross section for electron recoils is

[TABLE]

Here $\mu_{\chi e}$ is the electron-DM reduced mass, $\mu_{\chi e}=m_{e}m_{\chi}/(m_{e}+m_{\chi})$. 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 Chuzhoy and Kolb (2009); Hu et al. (2017); Dunsky et al. (2018) 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 Dvorkin et al. .

We must solve the full Boltzmann equation in an expanding background, given by

[TABLE]

where $C(p_{\chi},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_{\chi}$ to leading order and the Boltzmann equation can be solved by direct integration Bae et al. (2018),

[TABLE]

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

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

[TABLE]

where $E_{\bar{\chi}}=\sqrt{m_{\chi}^{2}+p_{{\chi}}^{2}+q_{e^{+}e^{-}}^{2}-2p_{\chi}q_{e^{+}e^{-}}\cos\theta}$, $E_{e^{+}e^{-}}=\sqrt{s_{e^{+}e^{-}}+q^{2}_{e^{+}e^{-}}}$ and $\theta$ is the angle that $\vec{q}_{e^{+}e^{-}}$ makes with the unconstrained, unintegrated $\vec{p}_{\chi}$. Defining $x\equiv\cos\theta$ and dropping the subscript on the bulk electron momentum, we find

[TABLE]

Requiring that $x\in[-1,1]$ and switching integration variables,

[TABLE]

Then, to solve for the final phase space from annihilation, we can combine Eqs. (25) and (28). Note that because $p_{\chi}$ is fixed (rather than an integration variable), $s$ in the above integral is restricted to $s>\max\left(4m_{e}^{2},2m_{\chi}(E_{\chi}+m_{\chi})\right)$ 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.

IV.2 Phase space from plasmon decay

The collision term from plasmon decay,

[TABLE]

proceeds through direct computation. We find

[TABLE]

where the limits of the $k$ integral are determined by the requirement that $x_{0}=(2E_{\chi}\omega_{\ell,t}(k)-m_{\ell,t}(k)^{2})/2kp_{\chi}$ 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 momentum-squared 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.

IV.3 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_{\chi}/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 Dvorkin et al. .

In the limit of relativistic DM scattering with non-relativistic SM particles (the case after freeze-in until $T_{\gamma}\sim m_{\chi}$), the differential cross section with respect to the center-of-mass scattering angle $\theta_{\text{CM}}$ is given by

[TABLE]

where $p_{\text{CM}}\equiv|\vec{p}_{\text{CM}}|$ is the momentum in the CM frame. Here we have taken $p_{\chi}\ll 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$ in the cross section). The dependence on the Debye mass $m_{D}$ comes from the photon propagator for electric scattering in a medium Blaizot et al. (1995). The usual $t$-channel divergence is thus regulated in the forward-scattering limit by the Debye angle, defined as $\theta_{D}\equiv m_{D}/p_{\text{CM}}$. Once the plasma has become non-relativistic with $T_{\gamma}\lesssim m_{e}$, the Debye mass is given by

[TABLE]

in natural units, assuming $\Omega_{b}h^{2}=0.022$ Aghanim et al. (2018) 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 $\theta_{D}\ll 1$

[TABLE]

Since $m_{b}\gg m_{\chi}$ 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_{\chi}$. 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

[TABLE]

where $n_{p}\approx 1.5\times 10^{-10}\,T_{\gamma}^{3}$ and $p_{\chi}\approx 0.4\,p_{\gamma}\approx T_{\gamma}$. For $T_{\gamma}$ 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

[TABLE]

where $\mu_{\chi b}$ is the reduced mass of the DM and baryon, $\mu_{\chi b}=m_{\chi}m_{b}/(m_{\chi}+m_{b})$, $v$ is the relative velocity between DM and SM particles, and $p_{\rm CM}=\mu_{\chi b}v$. The momentum transfer cross section is

[TABLE]

where again we take the $\theta_{D}\ll 1$ limit. Note that the Coulomb logarithm appearing here differs from the one that appears in the often-quoted Ref. McDermott et al. (2011); 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. McDermott et al. (2011), 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 Dvorkin et al. (2014); Xu et al. (2018); Slatyer and Wu (2018); Kovetz et al. (2018); Boddy et al. (2018), which we explore in more detail in our companion paper Dvorkin et al. .

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 Dubovsky and Gorbunov (2001); Dvorkin et al. (2014); Boddy and Gluscevic (2018); Boddy et al. (2018). The CMB bounds require that the momentum transfer rate is slow compared to the rate of Hubble expansion at $z\approx 1100$, thus limiting the possible effect on the DM phase space. We calculate the bounds in detail in the companion paper Dvorkin et al. , 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^{\prime}$, scattering processes such as $e^{-}+\gamma\to e^{-}+A^{\prime}$ are also possible and scale only as kinetic mixing squared $\kappa^{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.3. Processes like $\chi+\gamma\to\chi+A^{\prime}$ scale as $Q^{2}g_{\chi}^{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_{\chi}^{2}$, even at the largest values of $g_{\chi}$ that saturate SIDM bounds.

IV.4 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_{\chi}^{4}$, and thus may be important if $g_{\chi}$ is sufficiently large compared to $\kappa$. 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 non-relativistic limit, the quantity $\langle a(t)^{2}p_{\chi}^{2}\rangle$ 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_{\rm therm}$ when the momentum transfer rate per DM particle and per Hubble time is order unity:

[TABLE]

where $v$ is the relative velocity between DM particles and $\sigma_{T,\,\chi\chi}$ 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\approx p_{\chi}/m_{\chi}\approx T_{\gamma}(z)/m_{\chi}$. 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 $\tfrac{3}{2}T^{\rm eff}_{\chi}(z)\equiv\langle p_{\chi}^{2}\rangle/(2m_{\chi})$, where $p_{\chi}$ 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,

[TABLE]

where $n_{\rm DM}(z)$ is the DM number density.

Fig. 7 shows the redshift of thermalization for two representative choices of $\kappa$ (thus fixing $g_{\chi}$ 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_{\rm therm}$ with the redshift of recombination $z\approx 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_{\chi}$. 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_{\chi}$ 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 $\kappa=10^{-3}$ (which is consistent with bounds on ultralight dark photons), $g_{\chi}$ is small enough that DM self-interactions are not important at recombination and the phase space is described by the results of Sections IV.1-IV.2. The comparison of the free-streaming and thermalized phase space can be seen in Fig. 5.

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_{\chi,\,e^{+}e^{-}}(p_{\chi})\sim e^{-p_{\chi}/T_{\chi,\,e^{+}e^{-}}}$, where $T_{\chi,\,e^{+}e^{-}}$ is an effective DM temperature inherited from the electrons; at late times, this exponential distribution persists because the DM does not thermalize 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\rightarrow 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\sim\mathcal{O}\left(10^{-11}\right)$. 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 $\kappa$), it may be possible for DM self-scattering to thermalize the DM phase space distribution. In this case, $Q=\kappa g_{\chi}/e$ (where $\kappa$ can take on a wide range of values) and DM self-scattering via the dark photon scales as $g_{\chi}^{4}$, meaning that with the appropriate choice of $\kappa$ and $g_{\chi}$ it is possible to efficiently self-scatter while still producing the observed relic abundance. The coupling $g_{\chi}$ cannot be arbitrarily large due to observational limits on SIDM in astrophysical systems; however, there is a range of $g_{\chi}$ 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 $\left<p_{\chi}^{2}\right>$ 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_{\chi}(p_{\chi})$ to the equivalent quantities for the SM photon bath, we find that the typical DM momentum is similar to the typical photon momentum, $\langle p_{\chi}\rangle\approx(0.4-0.7)\times\langle p_{\gamma}\rangle$ 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 $\sim 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 freeze-in will lead to reduced drag force. Taking into account the updated Debye logarithm (which may weaken existing limits by a factor of $\sim 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 Dvorkin et al. , which will place restrictions on the range of masses where DM freeze-in via a light mediator is observationally viable.

Acknowledgments

We thank Masha Baryakhtar, Asher Berlin, Simon Knapen, Jung-Tsung Li, Adrian Liu, Aneesh Manohar, Sam McDermott, Julian Muñoz, and Tomer Volansky for helpful discussions. We acknowledge the importance of equity and inclusion in this work and are committed to advancing such principles in our scientific communities. CD was supported by NSF grant AST-1813694, Department of Energy (DOE) grant DE-SC0019018, and the Dean’s Competitive Fund for Promising Scholarship at Harvard University. TL is supported by an Alfred P. Sloan Research Fellowship and Department of Energy (DOE) grant DE-SC0019195. Parts of this paper were prepared while at the KITP, supported by the National Science Foundation under Grant No. NSF PHY-1748958. KS is supported by a National Science Foundation Graduate Research Fellowship and a Hertz Foundation Fellowship. KS is grateful for the hospitality of the Lisanti group at Princeton University and of the Center for Cosmology and Particle Physics at New York University, where part of this work was completed.

Appendix A Evolution of the SM bath

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

[TABLE]

Here we distinguish between the neutrino and SM bath temperatures $T$ and $T_{\nu}$; 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_{\nu}\sim 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,

[TABLE]

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_{\nu}(T_{\nu})a(T)^{3}$ is constant by definition. The remaining constraint equation on the temperature evolution is then

[TABLE]

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_{\nu}=(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,

[TABLE]

with $M_{\mathrm{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.

Appendix B In-medium plasma properties

In this Appendix, we follow the discussion of Ref. Braaten and Segel (1993), where the case of plasmons decaying to neutrinos was considered. The key approximation developed in that work was to evaluate thermal quantities at typical velocities, where thermal integrals have the most support. Specifically, the typical electron velocity is given by $v_{*}=\omega_{1}/\omega_{p}$, defined in terms of the first mode frequency and plasma frequency,

[TABLE]

[TABLE]

where $f_{e}$ is the phase space density of electron-positron pairs. Protons can also be included but their contribution is negligible because protons are heavy and thus slow to respond to electric fields, and also because their number density is much lower than that of the electrons at the relevant epochs.

The electromagnetic polarization tensor can be written as a thermal integral and expressed in terms of the longitudinal and transverse polarization functions, $\Pi_{\ell}$ and $\Pi_{t}$, as

[TABLE]

where $\omega$ and $\vec{k}$ are the plasmon energy and wavevector, and where the vacuum transverse polarization vectors $\vec{\epsilon}_{\pm}$ are chosen to be orthogonal to the direction of propagation and normalized to unity. In terms of the quantities above, the polarization functions can be approximated as

[TABLE]

[TABLE]

These approximations are accurate up to $\mathcal{O}(\alpha)$ and up to $\mathcal{O}(k^{2})$ at small $k$ for all electron temperatures and densities.

The effective propagator can then be constructed; in Coulomb gauge, its nonzero components are

[TABLE]

The poles in the propagator yield the renormalized longitudinal and transverse dispersion relations for on-shell plasmons,

[TABLE]

while the residues of the poles are identified as a combination of dressed polarization four-vectors, $\tilde{\epsilon}^{\mu}(k)\tilde{\epsilon}^{\nu}(k)^{*}$, for the appropriate polarization. The dressed polarization vectors are given by

[TABLE]

Given the approximations for $\Pi_{\ell}$ and $\Pi_{t}$ and the dispersion relations, the residue functions can be written as

[TABLE]

[TABLE]

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\kappa J_{\rm EM}^{\mu}A^{\prime}_{\mu}$. In a thermal plasma, this coupling generates an in-medium mixing term in the Lagrangian given by $\kappa A_{\mu}\Pi^{\mu\nu}A^{\prime}_{\nu}$ where $\Pi^{\mu\nu}$ is the electromagnetic polarization tensor. The matrix element in the dark photon basis is then given by

[TABLE]

where $D_{A^{\prime}}^{\nu\alpha}$ is the dark photon propagator. Taking the $m_{A^{\prime}}=0$ limit and working in Coulomb gauge, the propagator is given by

[TABLE]

Here we are ignoring in-medium corrections on the dark photon propagator, which are suppressed by factors of $\kappa^{2}$. Contracting $D^{\nu\alpha}_{A^{\prime}}$ with $\Pi_{\mu\nu}$ yields a vertex

[TABLE]

In the second line, we have assumed on-shell transverse and longitudinal modes for the respective pieces of the vertex function and used the dispersion relations of Eq. (51). Contracting this with a dressed polarization vector for the external photon yields

[TABLE]

which gives the same result as the vertex obtained in the millicharge basis.

Appendix D Regulating forward scattering

The differential DM-baryon scattering cross section can be written with respect to the CM angle $\theta_{\text{CM}}$ as

[TABLE]

In the limit where all of the particles are non-relativistic and where $m_{A^{\prime}}\ll m_{D}$ (if a dark photon is even present in the theory), the matrix element squared for DM-baryon Coulomb scattering is given by666Note that if the dark photon mass becomes large enough that it poses a relevant scale in the problem, then an additional factor of $q^{4}/(q^{2}-m_{A^{\prime}}^{2})^{2}$ appears to account for the in-medium couplings in Eq. (5).

[TABLE]

where we averaged over initial spins and summed over final spins. Here $q$ is the momentum transfer four-vector which satisfies $q^{2}=-2p_{\text{CM}}^{2}(1-\cos\theta_{\text{CM}})$ in the CM frame, $p_{\text{CM}}=|\vec{p}_{\text{CM}}|$ is the magnitude of the 3-momentum in this frame, and $m_{D}$ is the Debye mass. This effective mass arises from considering the longitudinal polarization tensor of the plasma $\Pi^{00}$ with the appropriate photon kinematics ($\omega\ll|\vec{q}|$) Blaizot et al. (1995), which corresponds to screened Coulomb scattering. It can also be understood as the effective mass appearing in the screened electric potential, which takes the form of a Yukawa potential Blaizot et al. (1995); Braaten and Yuan (1991); Raffelt (1986) or as a scale appearing in the electric form factor for a thermal Gibbs ensemble of charged particles in the plasma Raffelt (1986). Note that the transverse polarization tensor $\Pi^{ij}$, which corresponds to the magnetic scattering mode, vanishes in the static $\omega\ll|\vec{q}|$ limit Blaizot et al. (1995); however, this mode of scattering is negligible for a non-relativistic plasma where its contribution is suppressed by factors of $v$ Baym et al. (1990); Braaten and Yuan (1991).

The Debye mass automatically regulates the forward scattering divergence in the transfer cross section

[TABLE]

where in the second equality we have taken the approximation $s=(m_{b}+m_{\chi})^{2}$ for non-relativistic particles and have also taken the approximation $p_{\text{CM}}\gg m_{D}$. In the CM frame $p_{\text{CM}}=\mu_{\chi b}v$ where $\mu_{\chi b}=m_{b}m_{\chi}/(m_{b}+m_{\chi})$ is the DM-baryon reduced mass and $v$ is the relative velocity. If we had cut the integral by hand at some angle $\theta_{D}$ (rather than including the Debye mass in the propagator) we would have obtained a logarithm $\ln 2/\theta_{D}$ so we identify the correct Debye angle as $m_{D}/p_{\text{CM}}$. This Coulomb logarithm also agrees with other DM-baryon scattering rates found in the literature, for instance in Refs. Davidson et al. (2000b); Dubovsky et al. (2004); Raffelt (1986); Vogel and Redondo (2014).

This procedure yields a different logarithm than Ref. McDermott et al. (2011), which has been used for recent CMB constraints on millicharged DM. In that work, the angular integral was cut by using the relation between impact parameter and scattering angle for (electric) Coulomb scattering, and requiring that the impact parameter for pairwise DM-baryon scattering not exceed the Debye length $\lambda_{D}=1/m_{D}$. This translated to a minimum scattering angle that depended on the DM millicharge, with $\theta_{\rm min}=2Q\alpha/(3T\lambda_{D})$. The corresponding minimum momentum transfer in that case would be $|\vec{q}|^{2}=4Q^{2}\alpha^{2}p_{\rm CM}^{2}m_{D}^{2}/(9T^{2})$. For freeze-in where $p_{\rm CM}\approx T$ and $Q<10^{-10}$, we see that $|\vec{q}|^{2}\ll m_{D}^{2}$ 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. McDermott et al. (2011) 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 $\sim 2-3$ larger than the one obtained with the procedure of Eq. (62). As a result, CMB limits on millicharged DM that use this result may be too strong.