# Solar Reflection of Dark Matter

Haipeng An  Department of Physics, Tsinghua University, Beijing 100084, China
Center for High Energy Physics, Tsinghua University, Beijing 100084, China
Center for High Energy Physics, Peking University, Beijing 100871, China
Haoming Nie  Department of Physics, Tsinghua University, Beijing 100084, China
Maxim Pospelov  William I. Fine Theoretical Physics Institute, School of
Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA
Josef Pradler  Institute of High Energy Physics, Austrian Academy of Sciences,
1050 Vienna, Austria  Adam Ritz  Department of Physics and Astronomy,
University of Victoria, Victoria, BC V8P 5C2, Canada

(August 2021)

###### Abstract

The scattering of light dark matter off thermal electrons inside the Sun
produces a “fast” sub-component of the dark matter flux that may be detectable
in underground experiments. We update and extend previous work by analyzing
the signatures of dark matter candidates which scatter via light mediators.
Using numerical simulations of the dark matter-electron interaction in the
solar interior, we determine the energy spectrum of the reflected flux, and
calculate the expected rates for direct detection experiments. We find that
large Xenon-based experiments (such as XENON1T) provide the strongest direct
limits for dark matter masses below a few MeV, reaching a sensitivity to the
effective dark matter charge of $\sim 10^{-9}e$.

##  I Introduction

The evidence for dark matter (DM), through its gravitational signatures on
multiple astrophysical and cosmological scales, continues to stand as one of
the primary motivations for physics beyond the Standard Model. As the
technology of direct detection experiments improves, and experimental
sensitivity to dark matter in the galactic halo is pushed down toward the
threshold for elastic scattering of neutrino background fluxes (i.e. the
“neutrino floor”), there is increasing attention focused on a variety of DM
models that go beyond the traditional WIMP paradigm. Coupled to the broadening
scope of DM searches is the identification of “blind spots” for direct
detection, and associated experimental efforts to address them.

Until recently, one such problematic topic in direct detection has been the
realm of light WIMPs, where the momentum transfer in scattering often falls
below the detection threshold for the majority of DM experiments searching for
nuclear recoil. Light dark matter may exist well below the low mass end of the
traditional thermal relic WIMP window Lee and Weinberg (1977), due to the
possibility of light mediators playing a role in freeze-out in the early
universe Boehm and Fayet (2004); Pospelov  _et al._ (2008); Knapen  _et al._
(2017). In such models, scattering of DM on electrons may provide better
prospects for detection. To address the direct detection of DM in the MeV and
sub-MeV mass range via its scattering on electrons, new detection technologies
with ever decreasing energy thresholds (reaching down to $O({\rm eV})$ scales)
are being developed Battaglieri  _et al._ (2017); Barak  _et al._ (2020);
Aguilar-Arevalo  _et al._ (2019); Agnese  _et al._ (2014). At the current
time, these devices are relatively small, and as a consequence operate with
exposures that are many orders of magnitude smaller than the exposures
achieved in the large Xenon-based experiments Aprile  _et al._ (2017); Akerib
_et al._ (2017a); Tan  _et al._ (2016). An alternative pathway for generating
extra sensitivity is to account for sub-components of DM in full halo
distribution that are considerably more energetic than the primary galactic
flux, and therefore are visible to conventional underground detectors that
have higher thresholds, but also large volumes and exposures Kouvaris (2015);
An  _et al._ (2017); Emken  _et al._ (2018); Bringmann and Pospelov (2019);
Ema  _et al._ (2019); Dent  _et al._ (2020); Cappiello and Beacom (2019);
Harnik  _et al._ (2021); Ge  _et al._ (2021); Chen  _et al._ (2020); Emken
(2021); Bell  _et al._ (2021). These energetic sub-components have been
identified after carefully considering the impact of stars, cosmic rays, and
other astrophysical processes on dark matter in the halo distribution via non-
gravitational scattering.

In An  _et al._ (2017), we pointed out that there is an important direct
detection sensitivity to DM in the 10 keV – 10 MeV mass range due to a
population of ‘reflected DM’ generated through scattering by more energetic
electrons in the Sun (or the Earth) prior to intersecting the detector.
Reflective (single or multiple) scattering by solar electrons allows the
kinetic energy of a sub-component of the ambient dark matter in the halo to be
lifted into the keV range, making such light dark matter candidates visible
above the detector energy thresholds. The scenario is illustrated in Fig. 1,
where MeV-scale DM scattering off free electrons in the Sun generates a more
energetic component of the flux. While the fraction of that flux impinging on
the Earth is necessarily subject to a geometric suppression factor, it was
shown in An  _et al._ (2017) that this more energetic component of the
spectrum provides new sensitivity of multiple existing experiments to MeV-
scale DM models.

Denoting the reduced DM-$e$ mass $\mu_{{\rm DM},e}$, a single scatter can
result in the energy of the reflected DM being raised subject to the
constraint,

$E_{\rm DM}^{\rm refl}<E_{\rm DM}^{\rm refl,max}=\frac{4E_{e}\mu_{{\rm
DM},e}}{m_{e}+m_{\rm DM}}=\frac{4E_{e}m_{{\rm DM}}m_{e}}{(m_{e}+m_{\rm
DM})^{2}},$ (1)

which is much higher than $E_{\rm DM}^{\rm halo}$ and comparable to the
typical solar electron kinetic energy $E_{e}\sim kT_{e}\sim O({\rm keV})$.
Thus $E_{\rm DM}^{\rm refl}$ can be above the detection threshold of several
existing direct detection experiments. In the limit of small scattering cross
sections (and a nearly transparent Sun for DM), the signal of the reflected
flux scales quadratically with the scattering cross section
$\propto\sigma^{2}$, as two scattering events, inside the Sun and inside the
detector, are involved.

Figure 1: A schematic illustration of the reflected dark matter flux generated
through solar scattering.

The reflection mechanism is most effective for dark matter with a mass of
${\cal O}({\rm MeV})$, e.g. comparable to the electron mass. It is important
to keep in mind that DM models with mass in the MeV or sub-MeV mass range are
subject to a number of cosmological constraints from BBN and the CMB if they
ever achieved thermal equilibrium with Standard Model degrees of freedom in
the universe. In this paper, we extend the analysis of An  _et al._ (2017) to
include benchmark sub-thermal MeV-scale models of dark matter interacting via
a light vector mediator and with extremely small effective coupling to
electrons. Very small couplings may actually be responsible for the DM
abundance via the so-called “freeze-in” mechanism. Calculation of the
reflected flux in this case requires a somewhat different analysis of
scattering inside the Sun, due to the longer-range interaction, which allows
us to model DM scattering in using the model of a millicharge particle
propagating through the solar electron-ion plasma. We assess the sensitivity
of the Xenon1T experiment to the reflected DM component in such scenarios, and
compare these results with experiments that have achieved the lowest energy
thresholds for detection (e.g. SENSEI Barak  _et al._ (2020)).

Calculation of the experimental sensitivity to Reflected DM involves two
primary steps that we will describe in this paper; (i) computation of the
spectrum of reflected DM from the Sun, for which we use a simulation described
in Sections II, III and IV to describe scattering mediated by both heavy and
light mediators; and (ii) the computation of the experimental direct detection
sensitivity to DM-electron scattering, which is somewhat novel as the DM can
be semi-relativistic, and our approach is described in Section V. We present
our final sensitivity results for DM models with both heavy and light
mediators in Section 5, and finish with some concluding remarks in Section VI.

##  II Overview of DM scattering inside the Sun

The thermal plasma inside the Sun is composed of electrons and ions, with
temperatures from about 1 eV at the surface to about 1 keV at the center, and
for the purpose of this study it will be safe to assume that it is known
exactly. To fix our notation, we write the general scattering rate for DM
inside the plasma as follows

$\displaystyle\Gamma_{\chi}=\frac{1}{2k_{1}^{0}}\int\frac{d^{3}k_{2}}{2k_{2}^{0}(2\pi)^{3}}\int\frac{d^{3}p_{1}}{2p_{1}^{0}(2\pi)^{3}}f(p_{1})\int\frac{d^{3}p_{2}}{2p_{2}^{0}(2\pi)^{3}}(2\pi)^{4}\delta^{4}(k_{1}+p_{1}-k_{2}-p_{2})\sum_{\rm
spin}|{\cal M}|^{2}\ ,$ (2)

where $k_{1}$ and $k_{2}$ are the four-momenta of the incoming and outgoing DM
particle, $p_{1}$ and $p_{2}$ are the four-momenta of the incoming and
outgoing electron or ion, $f$ is the momentum distribution of the electrons or
the ions, and ${\cal M}$ is the scattering matrix element. The sum is over the
final state spins, and the average is performed over the initial spins. We
will assume that the electron and ion number densities at each value of the
solar radius $r$ are distributed classically according to Maxwell-Boltzmann
statistics:

$\displaystyle n_{e,{\rm ion}}(r)=\int\frac{d^{3}p}{(2\pi)^{3}}f_{e,{\rm
ion}}(p,r)\ ,$ (3)

where

$\displaystyle f_{e,{\rm ion}}(p,r)=n_{e,{\rm
ion}}(r)\left(\frac{2\pi}{m_{e,{\rm
ion}}T(r)}\right)^{3/2}e^{-p^{2}/2m_{e,{\rm ion}}T(r)}\ .$ (4)

Here $n_{e,{\rm ion}}(r)$ is the radius-dependent number density that is given
by the standard solar model. We will suppress the explicit dependence on the
radius below.

Our goal is to obtain the differential scattering rate
$d\Gamma_{\chi}/dk_{2}d\cos\theta$ in the plasma frame as a function of the
incoming DM momentum $k_{1}$, where $\theta$ is the scattering angle, at each
value of the solar radius. With this information, we can simulate the
propagation of DM particles inside the Sun with the following strategy : we
first assume that the plasma distribution inside the Sun is isotropic. We then
discretize the radial direction into approximately ${N}_{\rm shell}=2000$
shells, and take the temperature and density in each shell as constant. We
further assume that the shells are thin enough that the probability for DM to
have one collision in each shell is much smaller than unity; this assumption
is justified a posteriori for the cross sections of interest. The initial
velocity of DM is a random variable that samples a truncated Gaussian
distribution. The width of the Gaussian $\exp(-v^{2}/(2v_{0}^{2})$ is given
approximately by $v_{0}\sim 10^{-3}\times c$, and the high-velocity truncation
occurs at the nominal escape velocity of $\sim 2\times 10^{-3}c$. The
reflected spectrum is not sensitive to this assumption. The impact parameter
is chosen randomly from 0 to $4R_{\odot}$ to incorporate gravitational
focusing.

Outside the Sun, the effect of gravitation is accounted for through angular
momentum and total energy conservation laws. Once inside the Sun, within each
shell, we assume the DM particle trajectory is a straight line. Gravity is
accounted for by calculating the change of velocity due to the force exerted
by the inner shells. Within each shell, with the temperature and electron and
ion densities taken from the solar model, we calculate the total and
differential scattering rates. The total probability determines if scattering
occurs in that shell. If scattering does not occur, the DM particle propagates
to the next shell on a straight trajectory. Gravitational bending of the
trajectory is accounted for at the boundary with the next shell, where the DM
velocity $\vec{v}_{\chi}$ is changed according to the gravitational influence
of the inner shells. Once scattering occurs in a shell, according to the total
probability, the differential scattering rate is used to randomly generate the
the new angle and momentum. This new after-scattering velocity is used to
propagate DM particle to the boundary of the next shell, and so forth. The
simulation continues until the DM particle leaves the Sun with a positive
total kinetic energy. If the total kinetic energy is negative, it will reenter
the Sun and repeat the above procedure. Finally, after shooting a sufficient
number of DM particles, we can statistically determine the energy distribution
of the reflected DM flux.

The final step in the simulation is to form the normalized distribution
$F_{A_{\rho}}$ of the DM flux traversing the “impact disc” of area
$A_{\rho}=\pi(4R_{\odot})^{2}$, $\int dE_{\chi}F_{A_{\rho}}=1$. The reflected
flux observed at Earth’s location is then determined An  _et al._ (2017) as
follows,

$\frac{d\Phi_{\rm reflected}}{dE_{\chi}}=\Phi_{\rm
halo}\times\frac{F_{A_{\rho}}A_{\rho}}{4\pi({\rm A.U.})^{2}},$ (5)

where the $A_{\rho}/(4\pi({\rm A.U.})^{2})$ ratio contains the main geometric
suppression of the reflected flux at the Earth. In addition, the kinetic
energy of the reflected particle at the Earth’s location is reduced according
to the change in the potential energy.

Figure 2: Left panel: temperature and electron density distributions inside
the Sun taken from Ref. Bahcall  _et al._ (2005). Right panel: temperature and
electron density distributions in the Sun’s corona reproduced from Ref. De La
Luz  _et al._ (2008).

##  III Reflected DM flux: scattering via contact interactions

In this section we describe in detail how to simulate the reflected energy
spectrum in the case of contact ($s$-wave) scattering. This treatment follows
previous work An  _et al._ (2017); Emken (2021). Taking DM-electron scattering
as the example, we can write the momentum-independent scattering matrix
element in terms of the total scattering cross section $\sigma_{\rm tot}$,

$\displaystyle\overline{|{\cal M}|^{2}}=16\pi(m_{\chi}+m_{e})^{2}\sigma_{\rm
tot}\ .$ (6)

Defining the momentum transfer as $q=k_{2}-k_{1}$, and replacing $d^{3}k_{2}$
by $d^{3}q$ we obtain,

$\displaystyle\Gamma_{\chi-e}=\frac{\pi(m_{\chi}+m_{e})^{2}\pi\sigma_{\rm
tot}}{m_{e}^{2}m_{\chi}^{2}}\int\frac{d^{3}p_{1}}{(2\pi)^{3}}f_{e}(p_{1})\frac{d^{3}q}{(2\pi)^{3}}\delta\left(q^{2}\left(\frac{1}{2m_{\chi}}+\frac{1}{2m_{e}}\right)+\frac{{\bf
k}_{1}\cdot{\bf q}}{m_{\chi}}-\frac{{\bf p}_{1}\cdot{\bf q}}{m_{e}}\right)\ .$
(7)

Integrating out the angular part of $\vec{p}_{1}$ first, we obtain

$\displaystyle\Gamma_{\chi}=\frac{(m_{e}+m_{\chi})^{2}\pi\sigma_{\rm
tot}}{m_{e}^{2}m_{\chi}^{2}}\int\frac{d^{3}q}{(2\pi)^{3}}\int_{p_{1{\min}}}\frac{p_{1}^{2}dp_{1}}{(2\pi)^{2}}f_{e}(p_{1})\frac{m_{e}}{p_{1}q}\
,$ (8)

where

$\displaystyle
p_{1,\min}=\left|q\left(\frac{m_{e}}{2m_{\chi}}+\frac{1}{2}\right)+\frac{m_{e}}{m_{\chi}}\frac{{\bf
k}_{1}\cdot{\bf q}}{q}\right|\ .$ (9)

Given the Maxwell-Boltzmann form of the distribution $f_{e}$, the integral
over $p_{1}$ can be carried out analytically so that the scattering rate takes
the form:

$\displaystyle\Gamma_{\chi}=\frac{(m_{e}+m_{\chi})^{2}\pi\sigma_{\rm
tot}n_{e}m_{e}}{m_{e}^{2}m_{\chi}^{2}}\left(\frac{2\pi}{m_{e}T}\right)^{1/2}\int\frac{d^{3}q}{(2\pi)^{3}}\frac{1}{q}\exp\left[-\frac{1}{2m_{e}T}\left(\frac{m_{e}}{2m_{\chi}}(2k_{1}\cos\theta_{q}+q)+\frac{q}{2}\right)^{2}\right]\
.$

It follows that the differential rate is given by

$\displaystyle\frac{d\Gamma_{\chi}}{dqd\cos\theta_{q}}\sim{\cal
N}q\exp\left[-\frac{1}{2m_{e}T}\left(\frac{m_{e}}{2m_{\chi}}(2k_{1}\cos\theta_{q}+q)+\frac{q}{2}\right)^{2}\right]\
,$ (11)

where ${\cal N}$ is a factor that does not dependent on the kinematic
variables of the process. The total scattering rate can then be calculated
directly from Eq. (2),

$\displaystyle\Gamma_{\chi}=\sigma_{\rm
tot}n_{e}\sqrt{\frac{T}{m_{e}}}\times\sqrt{\frac{1}{2\pi}}\left[2e^{-x_{\chi}^{2}/2}+(2\pi)^{1/2}(1+x_{\chi}^{2}){\rm
erf}\left(\frac{x_{\chi}}{\sqrt{2}}\right)\right]\ ,$ (12)

where $x_{\chi}=v_{\chi}(m_{e}/T)^{1/2}$. The probability for a DM particle to
scatter off an electron in an infinitesimal distance $\Delta l$ can then be
written as

$\displaystyle P=\Gamma_{\chi}\times\frac{\Delta l}{v_{\chi}}=\sigma_{\rm
tot}n_{e}\Delta
l\left[\sqrt{\frac{\pi}{2}}\frac{1}{x_{\chi}}e^{-x_{\chi}^{2}/2}+\frac{1+x_{\chi}^{2}}{x_{\chi}^{2}}{\rm
erf}\left(\frac{x_{\chi}}{\sqrt{2}}\right)\right]\ .$ (13)

With the differential scattering rate (11) and total scattering probability
(13) in hand, the reflected flux can be obtained following the strategy
presented in Sec. II. In Fig. 3 we show the resulting normalized flux for
$m_{\chi}=0.5$ MeV. The DM-electron cross section is varied from $10^{-39}$ to
$10^{-34}$ cm2. As the cross section increases, the DM particles experience
scattering with the electrons in the outer, more dilute layers of the Sun.
Consequently, the spectrum becomes softer. In Fig. 4, we compare the reflected
spectra for different values of $m_{\chi}$. In the non-relativistic limit, the
energy transfer-efficiency reaches its maximum when the masses of the
elastically scattering particles are equal. As seen in Fig. 4, the spectra for
$m_{\chi}=0.5$ MeV are harder than for $m_{\chi}=0.1$ and $2.5$ MeV. This
conclusion applies to both large and small values of the cross section.

Figure 3: Reflected DM energy spectrum for $m_{\chi}=0.5$ MeV. The
distributions result from cross sections ranging from $\sigma_{\rm
tot}=10^{-39}$ – $10^{-34}$cm2.  Figure 4: Reflected DM energy spectrum when
both DM mass and cross section are varied. The first number in paranthesis
corresponds to $m_{\chi}$ in MeV, and the second number to $\sigma_{\rm tot}$
in cm2.

###  III.1 The single scattering limit and the validation of the simulation

It is clear that when the value of the cross section is sufficiently small,
the probability of multiple scattering diminishes, and the single scattering
limit adequately describes the reflected spectrum. In this subsection, we
compare our simulations with a semi-analytic single-scattering formula in the
appropriate limit.

We can estimate the relevant cross section for the single-scattering
approximation by considering the scattering probability in the solar core. The
relevant mean free path is $l_{\rm fp}=[n_{e}\langle\sigma_{\rm
tot}v_{r}\rangle]^{-1}{\bar{v}_{\chi}}$ where ${\bar{v}_{r}}\sim{\bar{v}_{e}}$
and ${\bar{v}_{\chi}}$ are respectively the characteristic relative velocity
(approximately the electron velocity) and the DM velocity in the core. We can
then estimate the scattering probability as

$\displaystyle P_{s}\sim R_{\rm traj}/l_{\rm fp}\sim\rho_{\rm
core}/m_{p}\times R_{\rm traj}\times\sigma_{\rm
tot}\times\frac{\bar{v}_{e}}{\bar{v}_{\chi}}\sim\frac{\sigma_{\rm
tot}}{10^{-38}~{}{\rm cm}^{2}}\ ,$ (14)

where $\rho_{\rm core}$ is approximately $150$ gram$/$cm3. Inside the core,
the typical ratio of electron to DM velocities
${\bar{v}_{e}}/{\bar{v}_{\chi}}\sim 10$. $R_{\rm traj}$ refers to the typical
distance that a DM particle travels inside the core of the Sun, where the core
size is about one quarter of the solar radius.

Therefore, in the region where $\sigma_{\rm tot}\ll 10^{-38}$ cm2, the
probability of scattering for each DM particle traversing the Sun is much
smaller than one, and the reflected flux can be estimated by considering just
single scattering events. For such events, the DM trajectory prior scattering
is determined by gravity alone. As a result, the differential scattering
probability simplifies to

$\displaystyle\left.\frac{dP}{dE_{2}}\right|_{\rm trajectory}=\int_{\rm
trajectory}\frac{dl}{v_{\chi}}\frac{d\Gamma_{\chi}}{dE_{2}}(k_{1}(r(l)),T(r(l)),n_{e}(r(l)))\
,$ (15)

where $v_{\chi}$ is the speed of the DM particle at the moment of scattering.
The trajectory of the integral is solely determined by the initial velocity
and impact parameter; $E_{2}$ is the kinetic energy of the outgoing DM
particle. In the region of interest, we have $k_{2}\gg k_{1}$, and as a result
$q\approx k_{2}$. In this limit, we may integrate $d\cos\theta$ in Eq. (11)
and obtain

$\displaystyle\frac{d\Gamma}{dE_{2}}=\frac{\sigma_{\rm tot}n_{e}}{(2\pi
m_{e}T)^{1/2}}\frac{(m_{e}+m_{\chi})^{2}}{m_{e}m_{\chi}}\exp\left[-\frac{E_{2}}{T}\frac{(m_{e}+m_{\chi})^{2}}{4m_{e}m_{\chi}}\right]\
.$ (16)

One observes that in this limit, the differential rate is insensitive to
$k_{1}$. The normalized kinetic energy distribution of the reflected flux can
then be written as

$\displaystyle\frac{2}{R_{\rm max}^{2}}\int dR\,R\int
dv_{0}f_{\chi}(v_{0})\left.\frac{dP}{dE_{2}}\right|_{\rm trajectory}\ .$ (17)

In this expression, $R_{\rm max}$ is the maximum value of the impact
parameter, which we choose to be $4R_{\rm Sun}$ so as to include the
gravitationally focused component, and $f_{\chi}$ is the initial DM velocity
distribution. If the trajectory does not intersect with the Sun, the total
kinetic energy does not change. In Fig. 5 we show the normalized distributions
for $m_{D}=0.5$ MeV with $\sigma_{\rm tot}$ varied from $10^{-38.5}$ cm2 to
$10^{-37}$ cm2. The colored curves show the result using the simulation
described in Sec. II. The big dots show the result obtained using the single
scattering calculation. We observe that for $\sigma_{\rm tot}\ll 10^{-37}$
cm2, the single scattering results agree with the simulation. For larger
values of $\sigma_{\rm tot}$ one leaves the single-scattering regime and the
analytical approximation begins to over-estimate the reflected flux.

Figure 5: Comparison of the single-scattering approximation (heavy dots) with
numerical simulations (scatter points) for $m_{\chi}=0.5$ MeV and different
values of the cross section in units of cm2.

###  III.2 Including the effect of ions

The DM particles may also interact with ions, and the strength of the
corresponding cross section relative to the cross section on electrons is
model-dependent. While it is fixed for scattering via a dark photon type
mediator, it will be different for other mediators, in general. Since we
consider $m_{\chi}\ll m_{\rm nucleus}$, collisions with ions will only change
the direction of DM. If the collision happens when the DM particle goes inward
to core, the ions may shield DM from going into the core. However, if the
collision happens when the DM particle is on an outward trajectory without
having been accelerated non-gravitationally, the collision may turn DM back
into the core and enlarge its probability for up-scattering. Without a
detailed simulation it is therefore hard to discern whether ions enhance or
diminish the reflected flux.

The average velocity of ions inside the core of the Sun is about $10^{-3}$
times of the speed of light, and about one order of magnitude smaller than the
typical velocity of DM inside the core. Therefore, the scattering rate is
determined by the velocity of DM. As a result, relative to $\chi-$e
scattering, $\chi-$ion scattering suffers from an approximate suppression
factor of $v_{\chi}/v_{e}\sim 0.1$. Therefore, when the total scattering
probability on electrons is well below unity (single-scattering limit), ions
cannot change the reflected flux. In Fig. 6, we show the energy distributions
of the reflected flux for different choices of $\sigma_{\rm tot}^{\chi-e}$ and
$\sigma_{\rm tot}^{\chi-{\rm ion}}$. For each curve, the first number in the
legend shows the ratio $\sigma_{\rm tot}^{\chi-p}/\sigma_{\rm tot}^{\chi-{e}}$
and the second number shows the value of $\sigma_{\rm tot}^{\chi-{e}}$;
$m_{\chi}$ is chosen to be $0.5$ MeV. We observe that when $\sigma_{\rm
tot}^{\chi-e}$ is below roughly $10^{-38}$ cm2, the ions do not have any
noticeable impact. On increasing $\sigma_{\rm tot}^{\chi-e}$, ions suppress
the high-energy tail of the reflected spectrum. Therefore, we conclude that
the net effect of ions is to shield DM from the core.

Figure 6: Reflected spectrum upon inclusion of ions for $m_{\chi}$=0.5 MeV.
The first number in the legend is $\sigma_{\rm tot}^{\chi-{\rm
ion}}/\sigma_{\rm tot}^{\chi-{e}}$ and the second number is $\sigma_{\rm
tot}^{\chi-e}$ in units of cm2.

##  IV Reflected DM flux: scattering via light mediators

For reasons that we will discuss at length, DM scattering via a very light
mediator requires a more involved calculation of the reflected flux compared
to the case of contact interactions. Here we assume that the dark mediator has
a negligible mass compared to the momentum transfers considered, which in
practice amounts to a sub-eV range.111In some models it might be important
that we keep the mediator mass above a specific low threshold; for example,
the propagation of millicharged particles is affected by magnetic fields,
solar winds etc Dunsky  _et al._ (2019). A small dark photon mass above
$10^{-13}$eV, on the other hand, would suppress the effect of the magnetic
field beyond distances of 1000 km regardless of the effective charge. There
are two principal options to realize an almost massless dark mediator: one is
through kinetic mixing with the photon and the other is through a new force,
e.g. a vector boson coupled to the $U(1)_{e-\mu}$ or $U(1)_{e-\tau}$ currents
with a tiny coupling. In this section, we will discuss these two cases
separately.

###  IV.1 Light mediator through kinetic mixing

The Lagrangian in this case takes the form,

$\displaystyle{\cal
L}=-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-\frac{\kappa}{2}F_{\mu\nu}V^{\mu\nu}-\frac{1}{4}V_{\mu\nu}V^{\mu\nu}-\frac{1}{2}m_{V}^{2}V_{\mu}V^{\mu}+V_{\mu}J_{\chi}^{\mu}+A_{\mu}J^{\mu}\
,$ (18)

where $J_{\chi}^{\mu}$ is the dark sector current, and the corresponding
charges are subsumed into the definitions of currents. Since in this paper the
symbol $\varepsilon$ is reserved for the electric permittivity, we prefer to
denote the dark photon and the kinetic mixing parameter as $(V,\kappa)$ as
opposed to the more common notation of $(A^{\prime},\varepsilon)$. We will
refer to the charge of DM under the $V_{\mu}$ mediator as $e_{D}$.

In a general Feynman gauge, the propagator of the photon can be written as

$\displaystyle\Pi^{\mu\nu}=\frac{1}{\varepsilon_{L}({q^{0}}^{2}-q^{2})}\epsilon_{L}^{\mu}\epsilon_{L}^{\nu}+\frac{1}{\varepsilon_{T}{q^{0}}^{2}-q^{2}}\sum_{r}\epsilon_{Tr}^{\mu}\epsilon_{Tr}^{\nu}+\frac{q^{\mu}q^{\nu}}{\xi({q^{0}}^{2}-q^{2})^{2}}\
.$ (19)

In this expression, $\varepsilon_{L(T)}$ denotes the longitudinal (transverse)
electromagnetic permittivity, $q=|{\bf q}|$, and the summation is over
transverse indices. The last term, the gauge-dependent part of the propagator,
does not give a contribution to the scattering matrix element. In
predominantly Coulomb-type scattering, it is the 0-component of $J^{\mu}$ that
contributes the most. For non-relativistic (NR) scattering, the matrix element
can be written as

$\displaystyle{\cal M}$ $\displaystyle=$ $\displaystyle\langle
J_{\chi}^{0}\rangle\frac{1}{q^{2}-m_{V}^{2}}\kappa
q^{2}\frac{{\epsilon_{L}^{0}}^{2}}{\varepsilon_{L}({q^{0}}^{2}-q^{2})}\langle
J^{0}\rangle$ (20) $\displaystyle=$ $\displaystyle\frac{\kappa
q^{2}}{\varepsilon_{L}q^{2}(q^{2}-m_{V}^{2})}\langle
J_{\chi}^{0}\rangle\langle J^{0}\rangle\ .$

Given that in NR scattering the energy transfer is subdominant to the momentum
transfer, $q^{0}\ll|\vec{q}|$, we have

$\displaystyle{\cal M}\approx\frac{\kappa\langle J_{\chi}^{0}\rangle\langle
J^{0}\rangle}{\varepsilon_{L}(q^{2}+m_{V}^{2})}\ .$ (21)

It is therefore important to calculate the longitudinal permittivity
$\varepsilon_{L}$. We will do so using linear response theory, the details of
which are presented in Appendix A.

To present the results, let us introduce a dimensionless quantity $A$,

$A=\frac{q^{0}}{q}\left(\frac{m_{e}}{T}\right)^{1/2},$ (22)

in terms of which, the longitudinal permittivity is given by

$\displaystyle\varepsilon_{L}$ $\displaystyle=$ $\displaystyle
1+\frac{e^{2}n_{e}}{Tq^{2}}\left[1-2A\int\frac{xdx}{(2\pi)^{1/2}}\tanh^{-1}(x/A)e^{-x^{2}/2}\right]\
.$ (23) $\displaystyle\equiv$ $\displaystyle
1+\frac{e^{2}n_{e}}{Tq^{2}}F_{1}(A)$

where function $F_{1}$ is defined as

$\displaystyle
F_{1}(A)=1-\left(\frac{\pi}{2}\right)^{1/2}Ae^{-A^{2}/2}\left[{\rm
erfi}(A/\sqrt{2})-i\right]\ .$ (24)

In the system of units used in this paper $e^{2}=4\pi\alpha$. Not
surprisingly, the main dimensionful parameter controlling the permittivity is
the Debye length $\ell_{\rm Debye}$ defined here as $\mu_{\rm
Debye}^{2}=\ell_{\rm Debye}^{-2}=4\pi\alpha n_{e}/T$.

The corresponding result for the transverse permittivity is given by

$\displaystyle\varepsilon_{T}=1-\frac{e^{2}n_{e}}{Tq^{2}}\frac{1}{A}\left(\frac{\pi}{2}\right)^{1/2}e^{-A^{2}/2}\left[{\rm
erfi}\left(\frac{A}{\sqrt{2}}\right)-i\right]=1+\frac{e^{2}n_{e}}{Tq^{2}}\frac{1}{A^{2}}(F_{1}(A)-1)\
.$ (25)

As an explicit check of this result, one can take the limit where the photon
is nearly on-shell ($A\gg 1$) to confirm that both permittivities have a
common limit, $\varepsilon_{L(T)}\to 1-\omega^{2}_{p}/{q^{0}}^{2}$, with the
usual non-relativistic expression for plasma frequency,
$\omega_{p}^{2}=4\pi\alpha m_{e}^{-1}n_{e}$.

The above result can easily be generalized to include ions, leading to the
general form,

$\displaystyle\varepsilon_{L}=1+\frac{\mu_{\rm
Debye}^{2}}{q^{2}}\left[F_{1}(A_{e})+\sum_{i}\frac{Z_{i}^{2}n_{i}}{n_{e}}F_{1}(A_{i})\right]\
,$ (26)

where

$\displaystyle A_{i}=\frac{q^{0}}{q}\left(\frac{m_{i}}{T}\right)^{1/2}\ ,$
(27)

and $Z_{i}$, $m_{i}$ and $n_{i}$ are the charges, masses and number densities
of the ions. As a result, the effect of plasma screening differs if the target
in the collision is different, which is a signature of model dependence.
Furthermore, charged particles with different masses also have different
screening effects.

####  IV.1.1 Screening effects in $\chi$-e and $\chi$-ion collisions

In the case of a $\chi$-e collision, and for the hard collisions that we are
interested in,

$\displaystyle q^{0}\sim T,\;\;q\sim(Tm_{e})^{1/2}.$ (28)

As a result, $A_{e}$ is of order one. For the contribution from ions, since
$m_{i}\sim$ GeV, we have $A_{i}\gg 1$, and as discussed in in the previous
subsection,

$\displaystyle F_{1}(A_{i})\rightarrow-\frac{1}{A_{i}^{2}}\ll 1\ ,$ (29)

Therefore, for $\chi-e$ collisions, the contribution from ions to
$\varepsilon_{L}$ is suppressed by a factor of ${\cal O}(10^{3})$ relative to
the contribution from electrons, and is therefore negligible.

In the case of $\chi-$ion collisions we have $q^{0}/q\sim v_{i}$, where
$v_{i}$ is the typical speed of the ions. Therefore, for the contribution from
the electrons, we have $A_{e}\sim v_{i}/v_{e}\ll 1$. As a result, and as
discussed above, $F_{1}(A_{e})\approx 1$, and the screening contribution from
electrons can be described by regular Debye screening. For the contributions
from the ions, it is easy to see that $A_{i}$ can be of order unity. Therefore
we need to rely on the full result for $F_{1}$ as shown in Eq. (65).

###  IV.2 Scattering rates

####  IV.2.1 $\chi-$e scattering

From Eqs. (2) and (21) the scattering rate can written as

$\displaystyle\Gamma_{\chi-e}$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}k_{2}}{(2\pi)^{3}}\int\frac{d^{3}p_{1}}{(2\pi)^{3}}f_{e}(p_{1})\frac{\kappa^{2}e_{D}^{2}e^{2}q^{4}}{|q^{2}+\mu_{\rm
Debye}^{2}F(q^{0}/q)|^{2}(q^{2}+m_{V}^{2})^{2}}$ (30) $\displaystyle\times
2\pi\delta\left(\frac{(k_{1}^{2}-k_{2}^{2})}{2m_{\chi}}+\frac{p_{1}^{2}-(\mathbf{k}_{1}+\mathbf{p}_{1}-\mathbf{k}_{2})^{2}}{2m_{e}}\right)\
,$

where $\mu_{\rm Debye}^{2}=e^{2}n_{e}/T$ is the Debye screening parameter, and

$\displaystyle F(q^{0}/q)=F_{1}(A_{e})+b_{p}F_{1}(A_{p})+4b_{\rm
He}F_{1}(A_{\rm He})\approx F_{1}(A_{e})\ ,$ (31)

where $b_{p,{\rm He}}=n_{p,{\rm He}}/n_{e}$, and $F_{1}(A)$ is defined as
$\Delta\varepsilon_{L}/\Delta\varepsilon_{L}^{(0)}$; its real and imaginary
parts are shown in Fig. (12).

Following the derivation in the contact interaction case, the $e-\chi$
scattering rate can be written as

$\displaystyle\Gamma_{\chi-e}=\frac{\kappa^{2}e^{2}e_{D}^{2}n_{e}}{(2\pi)^{3/2}}\left(\frac{m_{e}}{T}\right)^{1/2}\\!\\!\int
dqqd\cos\theta_{qk_{1}}\frac{1}{|q^{2}+\mu_{\rm
Debye}^{2}F_{1}(A_{e})|^{2}}\frac{q^{4}}{|q^{2}+m_{V}^{2}|^{2}}\exp\left({-\frac{p_{1\rm
min}^{2}}{2m_{e}T}}\right),\qquad$ (32)

where

$\displaystyle p_{1\rm min}$ $\displaystyle=$
$\displaystyle\left|q\left(\frac{1}{2}+\frac{1}{2}\frac{m_{e}}{m_{\chi}}\right)+k_{1}\left(\frac{m_{e}}{m_{\chi}}\right)\cos\theta_{qk_{1}}\right|,$
$\displaystyle A_{e}$ $\displaystyle=$
$\displaystyle\frac{m_{e}}{m_{\chi}}\left(\frac{k_{1}}{\sqrt{m_{e}T}}\cos\theta_{qk_{1}}+\frac{1}{2}\frac{q}{\sqrt{m_{e}T}}\right).$
(33)

The rate $\Gamma_{\chi-e}$ can be simplified as follows,

$\displaystyle\Gamma_{\chi-e}=\frac{n_{e}\kappa^{2}e^{2}e_{D}^{2}}{(2\pi
T)^{3/2}m_{e}^{1/2}}\int
dx_{q}dc_{q}\exp\left[-\frac{1}{2}\left(\frac{1}{2}x_{q}+A_{e}\right)^{2}\right]\frac{1}{|x_{q}^{2}+\frac{\mu_{\rm
Debye}^{2}}{m_{e}T}F_{1}(A_{e})|^{2}}\frac{x_{q}^{4}}{(x_{q}^{2}+\frac{m_{V}^{2}}{m_{e}T})^{2}},$
(34)

where

$\displaystyle
A_{e}=v_{1}\left(\frac{m_{e}}{T}\right)^{1/2}c_{q}+\frac{1}{2}\frac{m_{e}}{m_{\chi}}x_{q}\
,$ (35)

with $c_{q}=\cos\theta_{qk_{1}}$.

####  IV.2.2 $\chi-$ion scattering

In the case of the $\chi-$ion scattering, we have

$\displaystyle\Gamma_{\chi-{\rm ion}}$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}k_{2}}{(2\pi)^{3}}\int\frac{d^{3}p_{1}}{(2\pi)^{3}}f_{\rm
ion}(p_{1})\frac{\kappa^{2}Z^{2}e_{D}^{2}e^{2}q^{4}}{|q^{2}+\mu_{\rm
Debye}^{2}F(q^{0}/q)|^{2}(q^{2}+m_{V}^{2})^{2}}$ (36) $\displaystyle\times
2\pi\delta\left(\frac{(k_{1}^{2}-k_{2}^{2})}{2m_{\chi}}+\frac{p_{1}^{2}-(\mathbf{k}_{1}+\mathbf{p}_{1}-\mathbf{k}_{2})^{2}}{2m_{\rm
ion}}\right)\ ,$

where in this case

$\displaystyle F(q^{0}/q)=F_{1}(A_{e})+b_{p}F_{1}(A_{p})+4b_{\rm
He}F_{1}(A_{\rm He})\approx 1+b_{p}F_{1}(A_{p})+4b_{\rm He}F_{1}(A_{\rm He})\
.$ (37)

Therefore, one finds

$\displaystyle\Gamma_{\chi-{\rm
ion}}=\frac{n_{i}\kappa^{2}e^{2}e_{D}^{2}Z_{i}^{2}}{(2\pi
T)^{3/2}m_{i}^{1/2}}\int
dx_{q}dc_{q}\exp\left[-\frac{1}{2}\left(\frac{1}{2}x_{q}+A_{i}\right)^{2}\right]\frac{1}{|x_{q}^{2}+\frac{\mu_{\rm
Debye}^{2}}{m_{i}T}F(q^{0}/q)|^{2}}\frac{x_{q}^{4}}{(x_{q}^{2}+\frac{m_{V}^{2}}{m_{i}T})^{2}},$
(38)

where

$\displaystyle F(q^{0}/q)=1+\sum_{i}F_{1}(A_{i})\ .$ (39)

Both rates, $\Gamma_{\chi-{\rm ion}}$ and $\Gamma_{\chi-e}$, exhibit a strong
dependence on scattering angle/momentum transfer.

####  IV.2.3 Numerical results in the small $m_{V}$ limit

The $m_{V}\to 0$ case is the most complex since Coulomb scattering is strongly
enhanced in the forward region. Even with screening, the total scattering
rate,

$\displaystyle\Gamma_{\chi-e}=\frac{n_{e}\kappa^{2}e^{2}e_{D}^{2}}{(2\pi
T)^{3/2}m_{e}^{1/2}}\int
dx_{q}dc_{q}\exp\left[-\frac{1}{2}\left(\frac{1}{2}x_{q}+A_{e}\right)^{2}\right]\frac{1}{|x_{q}^{2}+\frac{\mu_{\rm
Debye}^{2}}{m_{e}T}F_{1}(A_{e})|^{2}}\ ,$ (40)

is still too large. However, since we want the DM particle to be accelerated
inside the Sun, we are more interested in large angle scattering. Therefore we
can add an IR cutoff to $x_{q}$ that sets

$\displaystyle x_{q}=\frac{q}{\sqrt{m_{e}T}}=\frac{\Delta
v_{1}m_{\chi}}{\sqrt{m_{e}T}}>\frac{\zeta
v_{1}m_{\chi}}{\sqrt{m_{e}T}}=\zeta\frac{m_{\chi}}{m_{e}}v_{1}\sqrt{\frac{m_{e}}{T}}\
.$ (41)

One observes that the physical meaning of the cut-off is that for each
scattering event we neglect the case that the magnitude of the change of ${\bf
v}_{1}$ is smaller than $\zeta v_{1}$.

Then at each step of the simulation, just as in the case of contact
interactions, we use the total rate (with the IR cut-off in place) to
calculate the scattering probability. If scattering occurs, just as before we
use the differential rate $d\Gamma_{\chi-e}/dc_{q}dx_{q}$ to generate the
direction and the final momentum of the outgoing dark matter.

With the IR cutoff in place, scattering events can be separated into two
categories; namely, hard scattering with $x_{q}>x_{q}^{\rm cut}$, and soft
scattering with $x_{q}<x_{q}^{\rm cut}$. Hard scattering events can accelerate
the DM particles and generate a hard reflected spectrum whereas soft
scatterings will provide effective friction and soften the reflected spectrum.
The smaller the value of $x_{q}$, the more the viscous effects are included in
the simulation. In Fig. 7, we show the reflected spectrum with $m_{V}=0$, and
$e_{D}\kappa=10^{-10}$ for different choices of $\zeta$, from $\zeta=5$ to
100. The total number of DM particles we shoot into Sun is $10^{7}$, and the
initial conditions are chosen as described in Sec. II. For $\zeta=100$ (as
shown by the purple curve), only very hard scattering events remain, and the
results are qualitatively similar to the contact interaction case as discussed
in the previous sections. One observes clearly that decreasing $x_{q}$ softens
the spectrum and also enhances the soft spectrum overall. Furthermore, we can
see that the curves for $\zeta=5$ and $\zeta=10$ completely coincide with each
other. Therefore, we can conclude that a $x_{q}$ cut with $\zeta<10$ is
sufficiently accurate for us to capture the relevant soft scattering effects.
In practice, our computational power allows us to use $\zeta$ as small as 0.5,
which should be able to capture all the soft scattering effects. Comparing the
blue dotted curve with the reflected curves in contact scattering case, we
observe that in the light mediator case, the spectrum is much softer, as
expected. This suggests that the use of detectors with low thresholds would be
advantageous in the search for the reflected flux in light mediator DM models.

Figure 7: The dependence of simulated reflected spectrum on the IR cutoff
(labeled by values of $\zeta$). Small values for the cutoff lead to a
“migration” of the hard tail of the reflected spectrum to softer values due to
the effective viscosity. The similarity of the blue and red curves indicates
that values of $\zeta<10$ are adequate for obtaining realistic results.

###  IV.3 Scattering from the Solar Corona

While the main focus of our work is a derivation of the sensitivity to the
weakest couplings, it is also important to appreciate the limitations of our
method when the couplings and the cross sections become large, and the mean
free path for $\chi$ becomes so small that DM scatters too quickly to pass
into the interior of the Sun. In such a situation, one can entertain the
possibility of DM reflecting from the solar corona, that has a small column
density, but rather large temperatures (that are indeed $\sim$ two order of
magnitude larger than the temperature of the photosphere). The temperature and
electron distributions in the solar atmosphere are shown in the lower panel of
Fig. 2). We have also simulated the reflected flux from the corona, finding it
to be very small, but more energetic than the flux reflected from the
photosphere.

We will consider the possibility of a signal in the SENSEI surface run, which
will be discussed in detail in Sec. V.2. In this case, the kinetic energy of
the reflected DM needs to be above several eV. To obtain a qualitative
estimate of the reflected spectrum, we can use the approximation described in
Sec. III.1. Then, similar to Eq. (16), one obtains the differential scattering
rate

$\displaystyle\frac{d\Gamma_{\chi-e}}{dE_{2}}\approx\frac{\kappa^{2}e^{2}e_{D}^{2}n_{e}m_{e}}{4\pi(2\pi
m_{e}T)^{1/2}m_{\chi}E_{2}^{2}}\exp\left[-\frac{E_{2}}{T}\frac{(m_{\chi}+m_{e})^{2}}{4m_{\chi}m_{e}}\right]\
.$ (42)

The Sun’s corona extends to radii $10^{6}-10^{7}$ km with an electron density
smaller than $10^{9}$ cm-3 De La Luz  _et al._ (2008). For example, from Eq.
(42) one concludes that for $\kappa e_{D}\sim 10^{-4}$, the probability for DM
to be reflected with, say, $E_{2}>15$ eV (leading to a multi-electron signal
in a silicon detector) within the corona is smaller than one percent.
Therefore, we may use the single scattering approximation to estimate the
reflected spectrum with the differential scattering probability given by Eq.
(15).

Due to the Sun’s gravity, the trajectories of DM particles are hyperbolic. If
the initial impact parameter is small enough the DM particle may hit the
surface of the Sun. Once this happens, for a milli-charge $\kappa e_{D}$ as
large as $10^{-5}$ the DM particles can quickly achieve thermal equilibrium
with the electrons at the surface of the Sun, and then escape with eV-scale
kinetic energy in a random direction. The corona-accelerated DM particles may
also hit the surface of the Sun and lose their acquired kinetic energy. Such
circumstances complicate the simulation of the accelerated spectrum. As we are
primarily interested in a first order estimation of the corona-reflected flux,
we will neglect these complications. Concretely, we simply approximate the
trajectories as straight lines, and neglect the change of the velocity due to
gravity. The differential scattering probability then becomes,

$\displaystyle\frac{dP}{dE_{2}}\approx\left\langle\frac{1}{v_{\chi}}\right\rangle\int_{0}^{R_{\rm
c}}\frac{4\pi r^{2}dr}{\pi R_{\rm c}^{2}}\frac{d\Gamma_{\chi-e}}{dE_{2}}(r)\
,$ (43)

where $R_{c}$ is the radius enclosed by the corona; we choose
$R_{c}=1.69\times 10^{6}$ km (the solar radius is $6.9\times 10^{5}$ km). The
differential scattering probability for $m_{\chi}=0.2$ and 0.5 MeV with
$\kappa e_{D}=10^{-5}$ is shown in Fig. 8. As shown in the lower panel of Fig.
2, the coronal temperature remains around $10^{6}$ kelvin at a height of
$10^{6}$ km. The the index for the power law decay of the density $n_{e}$ at
this height is slower than $-3$. Therefore, we do not expect the contributions
to $dP/dE_{2}$ from the region beyond $r=1.69\times 10^{6}$ km to be entirely
negligible. However, given limited data on the coronal temperature and density
profiles at even larger radii, we limit $R_{c}=1.69\times 10^{6}$ km as our
largest radius. In this sense, our sensitivity result is conservative.

Figure 8: Differential scattering probability for $m_{\chi}=0.5$ MeV (red) and
0.2 MeV (blue) with $\kappa e_{D}=10^{-5}$.

##  V Direct detection of the solar reflected flux

In this section we shall use the solar reflected DM fluxes and spectrum
obtained above to set constraints from direct detection experiments. We focus
on the XENON1T experiment, and begin by introducing the standard components
that go into the computation of these limits.

The DM-electron cross section is conventionally written in terms of a
reference value, where the DM-electron three-momentum transfer is evaluated at
the typical atomic scale $q_{0}=\alpha m_{e}$,
$\bar{\sigma}_{e}\equiv\mu_{e}^{2}\overline{|\mathcal{M}(q=q_{0})|^{2}}/(16\pi
m_{\chi}^{2}m_{e}^{2})$, where $\overline{|\mathcal{M}(q)|^{2}}$ is the
squared matrix element for the scattering on a free electron and $\mu_{e}$ is
the DM-electron reduced mass. When the interactions are mediated by a
kinetically mixed vector (dark-photon), the cross section reads

$\displaystyle\bar{\sigma}_{e}=\frac{16\pi\kappa^{2}\alpha\alpha_{D}\mu_{e}^{2}}{(q_{0}^{2}+m_{V}^{2})^{2}},$
(44)

where $m_{V}$ is the dark photon mass and $\alpha_{D}=e_{D}^{2}/4\pi$ with
$e_{D}$ the gauge coupling in the dark sector. The actual dependence of the
cross section on momentum transfer is shifted into a “DM form factor”,

$\displaystyle|F_{\rm
DM}(q)|^{2}\equiv\frac{\overline{|\mathcal{M}(q)|^{2}}}{\overline{|\mathcal{M}(q=q_{0})|^{2}}}=\begin{cases}1&\text{contact}\\\
\frac{m_{V}^{2}+q_{0}^{2}}{m_{V}^{2}+q^{2}}&\text{light med.}\end{cases}\,.$
(45)

Notably, in the limit $m_{V}\to 0$ one recovers the millicharged DM model
where $Q_{\rm eff}=e_{D}\kappa/e$ is the fractional charge and $F_{\rm
DM}=(\alpha m_{e}/q)^{2}$.

###  V.1 Limits from XENON1T

Figure 9: Direct detection limits derived in this work from XENON1T S2 Aprile
_et al._ (2019a) and S1+S2 Aprile  _et al._ (2020a) data from solar reflected
DM in the limit of point-like interactions $F_{\rm DM}=1$ (orange shaded
regions as labeled.) The blue shaded regions are previous direct detection
limits using the galactic flux component from XENON1T Aprile  _et al._
(2019a), XENON10 Essig  _et al._ (2017), and SENSEI Barak  _et al._ (2020).
The gray shaded regions to the left are stellar cooling constraints from red
giant (RG) stars Chang  _et al._ (2021) and from the proto-neutron star of
SN1987a Chang  _et al._ (2018) for a dark photon-mediated model with
$\alpha_{D}=0.5$ and $m_{V}=3m_{\rm DM}$. The freeze-out curve is for a
complex scalar DM model An  _et al._ (2017). Left panel: only scattering with
electrons is assumed; the solid orange line shows the data-driven limit on the
charge yield where the lowest energy deposition is 0.19 keV Akerib  _et al._
(2017b). Right panel: various values for the scattering cross section on ions,
in relation to electrons are assumed as labeled. The additional scattering
channel has the net effect of weakening the limits, because particles are
reflected at larger, hence colder, solar radii.

Figure 10: Direct detection limits from solar reflected DM as in Fig. 9 but
for massless mediator interactions or millicharged DM, $F_{\rm DM}=(\alpha
m_{e}/q)^{2}$. The standard galactic direct detection limits are from SENSEI
Barak  _et al._ (2020), XENON10 Essig  _et al._ (2017), and DAMIC Aguilar-
Arevalo  _et al._ (2019). The astrophysical SN1987a Chang  _et al._ (2018) and
RG Vogel and Redondo (2014) cooling constraints are shown as before. The
freeze-in line is taken from Chang  _et al._ (2021); see also Dvorkin  _et
al._ (2019). Left panel: constraints in terms of the DM-electron reference
cross section evaluated at a momentum transfer $q_{0}=\alpha m_{e}$. Right
panel: constraints in terms of the millicharge measured in units of the
elementary charge $Q_{\rm eff}$; in the dark photon model with negligible
vector mass, $Q_{\rm eff}=\kappa e_{D}/e$holds, where $e_{D}$ is the gauge
coupling in the dark sector.

The velocity averaged electron recoil cross section for ionization from an
initial atomic state with principal and angular quantum numbers $n$ and $l$
can be written in a standard form Essig  _et al._ (2012),

$\displaystyle\frac{d\langle\sigma_{n,l}v\rangle}{d\ln E_{e}}$
$\displaystyle=\frac{\bar{\sigma}_{e}}{8\mu_{e}^{2}}\int dq\,\left[q|F_{\rm
DM}(q)|^{2}|f^{\rm ion}_{nl}(p_{e},q)|^{2}\eta(v_{\rm min}(q,\Delta
E_{n,l}))\right].$ (46)

Here, $|f^{\rm ion}_{nl}(p_{e},q)|^{2}$ is an ionization form factor which
depends on momentum transfer $q$ and final state electron momentrum through
$E_{e}=p_{e}^{2}/2m_{e}$. The flux-average of the inverse speed is denoted by
$\eta(v_{\rm min})$. In the non-relativistic limit it reads,

$\displaystyle\eta(E_{\rm min})=\int_{E_{\rm min}}dE\sqrt{\frac{m_{\rm
DM}}{2E}}\,\frac{1}{\Phi_{\rm halo}}\frac{d\Phi_{\rm refl}}{dE}\ ,$ (47)

where $E_{\rm min}$ is the minimum incident DM energy required to produce an
electron of recoil energy $E_{e}$ after lifting it into the continuum by an
investment of the binding energy $|E_{n,l}|$. In (47) we normalize the
integral to the halo flux $\Phi_{\rm halo}$; the total ionization rate from
the $(n,l)$ orbital is then given by $dR_{nl}/d\ln E_{e}=N_{T}\Phi_{\rm
halo}d\langle\sigma_{nl}v\rangle/d\ln E_{e}$, where $N_{T}$ is the target
density.

The energy deposition $\Delta E=E_{e}+|E_{n,l}|$ induces an ionization signal
in xenon (S2), and, for $E\gtrsim 1$ keVee, the scintillation signal (S1). In
this work we use XENON1T S1+S2 data from Aprile  _et al._ (2020b) and the
S2-only data from Aprile  _et al._ (2019b). We follow our previous work in
modeling the formation of the S1 and S2 signals and limit setting procedure An
_et al._ (2017). For the ionization form factor $|f^{\rm
ion}_{nl}(p_{e},q)|^{2}$ we use the calculation described in Essig  _et al._
(2020). It is based on numerical solutions of the Schroedinger equation for
bound and continuum electrons, and hence yields appropriate results accross
the entire momentum-transfer regime.

The direct detection limits from solar reflection for contact interactions
($F_{\rm DM}=1$) are shown in Fig. 9 with the orange shaded exclusion regions
for XENON1T S2 Aprile  _et al._ (2019a) and S1+S2 Aprile  _et al._ (2020a) as
labeled. As can be seen, they attain their maximum strength for $m_{\rm
DM}\simeq m_{e}$ corresponding to the point of most efficient energy transfer
in the sun. A variety of complementary constraints are described in the figure
caption. The limiting constraint (solid orange line) shows the data-driven
limit on the charge yield where we demand an energy deposition of at least
0.19 keV congruent with available data Akerib  _et al._ (2017b). As shown in
the right panel, the additional scattering on ions has the net effect of
weakening the limits, because particles are reflected at larger, hence colder,
solar radii.

It is instructive to compare the peak sensitivities of XENON1T to galactic
dark matter (blue region) with the sensitivity to reflected dark matter
(orange region). The strongest constraints on galactic dark matter scattering
on electrons, $\sigma\sim 10^{-40}\,{\rm cm}^{2}$, are derived for
$m_{\chi}\sim 200-300$ MeV. For the reflected dark matter flux, the best
sensitivity is for $\sigma\sim 5\times 10^{-38}\,{\rm cm}^{2}$ for $S1+S2$ and
$\sigma\sim 10^{-37}\,{\rm cm}^{2}$ for S2 only at masses that are three
orders of magnitude smaller, $m_{\chi}\sim 100-300$ keV. Notice that the
lowering of sensitivity by a factor of $10^{3}$ is composed from the two main
factors, large and small. The geometric loss, proportional to the inverse
solid angle subtended by the Sun, $1/\Delta\Omega_{\odot}\sim 10^{6}$, is
mitigated by the much larger number densities of dark matter at $m_{\chi}\sim
100-300$ keV, leading to an overall loss in sensitivity by only a factor of
$10^{3}$. (This rough comparison is possible because the same experiment is
used in both cases, XENON1T.)

Results for the millicharged case ($F_{\rm DM}=(\alpha m_{e}/q)^{2}$) are
shown in Fig. 10. Complementary limits are again described in the figure
caption. In the left panel, constraints in terms of the DM-electron reference
cross section is shown. In the right panel, the same constraints are shown in
terms of the fractional elementary charge $Q_{\rm eff}=\kappa e_{D}/e$. We
observe that the limits of $Q_{\rm eff}$ are strong, and relatively mass-
independent, surpassing the SN bounds by a slim margin. It is easy to see that
the best limits are imposed by the S2 only signal, and this is a direct
consequence of the softness of the reflected spectrum in this case.

###  V.2 Detecting the reflected flux from the solar corona

Figure 11:  An extension of the right panel of Fig. 10 to larger effective
couplings $Q_{\rm eff}$. The currently allowed white region for $Q_{\rm
eff}\gtrsim 5\times 10^{-6}$ opens up from below due to trapping of DM inside
the SN, while it is bounded from above by the SLAC millicharge experiment
Prinz  _et al._ (1998). A surface run by the SENSEI experiment with a small
exposure of 0.019 g-days Crisler  _et al._ (2018) and zero events in the
$5e^{-}$ bin allows us to place a superseded constraint using DM reflected
from the solar corona (solid orange line.) A projection to 1 kg-yr exposures
assuming a background rate of $100$/kg/day (zero) in the $5e^{-}$ bin gives
the upper (lower) orange dashed curves. The lower curve can be regarded as the
most optimal reach for a SENSEI-type surface detector; the actual sensitivity
is likely to be limited significantly by backgrounds. The galactic DM
constraints obtained with underground searches are the same as in Fig. 10, but
disappear for values in excess of $Q_{\rm eff}\sim 10^{-5}$ Emken  _et al._
(2019) as indicated by the hatching.

We close the exploration of direct detection sensitivity by commenting on the
possibility of probing larger cross sections. In particular, for millicharged
DM there exists a currently unconstrained region above the SN1987a bound
$Q_{\rm eff}\simeq 5\times 10^{-6}$ Chang  _et al._ (2018). This is due to the
trapping of pair-produced $\chi$-particles inside the supernova. The region
has a ceiling at $Q_{\rm eff}\lesssim\text{few}\times 10^{-5}$ due to the
direct search at the SLAC millicharge experiment Prinz  _et al._ (1998). Both
bounds are independent of the DM-nature of $\chi$. In this parameter region,
reflected DM experiences re-scattering in the Earth’s atmosphere and in any
overburden of a direct detection experiment. In particular, we do not expect
limits from the XENON1T experiment, as the diffusive process softens and
diminishes the spectrum too much. Prospects for probing this region may
nevertheless exist through the operation of various low-threshold detectors
above ground.

For example, the SENSEI semiconductor direct detection experiment had a
surface commissioning run with a very small exposure of $0.019$ gram-days
Crisler  _et al._ (2018). The threshold for ionizing an electron is
$E_{\mathrm{gap}}=1.11~{}\mathrm{eV}$ and it takes an average of
$\varepsilon=3.6~{}\mathrm{eV}$ for the creation of an electron-hole pair (at
least, as measured in high-energy recoils). For example, the corona-reflected
flux is capable of creating a charge multiplicity of at least $5e^{-}$ in a
scattering with $E_{e}\geq 15.5$ eV electron recoil energy, and for which the
SENSEI run did not observe any events; (the integer charge multiplicity may be
estimated as
$Q\left(E_{e}\right)=1+\left\lfloor\left(E_{e}-E_{\mathrm{gap}}\right)/\varepsilon\right\rfloor$
Essig  _et al._ (2016)).

In Fig. 11 we show the resulting limit for the SENSEI surface run using the
$5e^{-}$ bin which had zero counts (orange solid line). As can be seen, it is
superseded by the SLAC experiment. However, we also show the results of two
projections for an effective exposure of $1$ kg-day for a SENSEI-type
experiment. For the upper dashed curve, we assume a background-rate in the
$5e^{-}$ bin of $100$ /kg/day (dru). This appears to be a reasonable ballpark
figure based on a simple extrapolation of the $2,3,$ and 4 electron rates for
the SENSEI surface-run. The lower dashed curve assumes (an unrealistic) zero
observed events in the $5e^{-}$ bin. This should therefore be regarded as the
optimal reach for a SENSEI-type surface detector. We obtain these lines
assuming no re-scattering in the Earth’s atmosphere. We caution that this
approximation may be subject to corrections, especially for the line of
current sensitivity. The evaluation of semiconductor ionization rates, in a
generalization of (46) to crystals, follows Essig  _et al._ (2016) from where
we also take the ionization form factor. Complementary constraints from other
sources are again described in the figure caption.

##  VI Concluding Remarks

The fast flux of MeV-scale dark matter reflected from the Sun is an inevitable
consequence of scattering on solar electrons. This sub-dominant, but highly
energetic component of the total dark matter flux on Earth provides a novel
pathway to study light dark matter scattering on electrons. While the
viability of low mass DM motivates R&D involving new materials and new ways of
registering ever lower recoil energies in direct detection relevant for lower
values of $m_{\chi}$ Knapen  _et al._ (2017), the solar reflection mechanism
provide a nontrivial level of sensitivity using data from the largest and
cleanest dark matter detectors based on Xenon dual phase technology. As a
consequence, dark matter electron scattering limits can be imposed now,
without relying on the success of future experimental efforts. A galactic DM
particle with mass in the 100 keV to 1 MeV range has sub-eV kinetic energy and
falls below the experimental threshold of any existing experiment. The
reflected flux provides a concrete means of probing this part of the parameter
space, and current data from XENON1T imposes limits down to cross sections of
$10^{-37}$cm2 in the case of contact scattering on electrons.

While the results for contact DM-electron scattering are an update to our
earlier work An  _et al._ (2017) (see also Emken (2021)), the derivation of
the reflected flux for DM interactions induced by an ultra-light mediator, and
the resulting limits from direct detection presented in this paper, represent
a nontrivial extension. The method can now be applied to the most viable sub-
MeV DM models, which achieve their sub-thermal relic abundance via freeze-in.
The technical challenge to overcome was the strong enhancement of small angle
scattering, which leads to the predominant occurrence of soft collisions. The
forward singularity is of course regularized by Debye screening. In this
paper, we have described a consistent derivation of this screening for the
relevant range of energy and momentum transfers, $q^{0},q$, and implemented
this approach in simulations to generate the reflected spectrum. The numerical
treatment was made viable by introducing an artificial IR cutoff, which we
varied to confirm that the final spectrum was independent of this parameter.

The reflected spectrum for a nearly massless mediator is considerably softer
than in the case of contact interaction. This is expected, given that frequent
soft scatterings downgrade the energetic component of the spectrum. As a
result, most of the reflected spectrum leads to ionization in XENON1T that
falls below the threshold for registering the S1 signal, and the best
constraints are imposed by examining the S2 only part of the signal. The
resulting limits in the mass interval for $m_{\chi}$ from 0.1-to-1 MeV vary
from $10^{-35}$cm2 to $10^{-33}$cm2, when the cross sections are evaluated at
$q=\alpha m_{e}$. The corresponding limits on the effective charge are around
the $Q_{\rm eff}=10^{-9}$ benchmark, and are slightly stronger than the
conservative limits from analysis of SN cooling. These are the most sensitive
constraints by some margin on DM in the sub-MeV mass range obtained using
terrestrial probes. Despite the strength of these bounds, we find that the
freeze-in prediction for $Q_{\rm eff}(m_{\chi})$ is beyond reach using the
reflected flux with existing technology.222Note however, that this need not
imply an overproduction of dark matter as the $\chi\bar{\chi}\to 2V$
annihilation channel remains open, and with the right choice of $e_{d}$ can
ensure correct DM relic abundance for the parameters at the current level of
sensitivity.

Note Added: As this paper was being finalized, we became aware of recent
complementary work exploring cosmological constraints on millicharge DM Buen-
Abad  _et al._ (2021); Nguyen  _et al._ (2021). The resulting sensitivity to
the electron scattering cross section of sub-MeV dark matter is quite similar.
However, the results are complementary as the direct detection limits on
reflected DM degrade rather weakly in the case of fractional DM abundance,
i.e. even a sub-1% fraction of the total DM will have its scattering on
electrons constrained via solar reflection, while cosmological bounds are
fully relaxed.

###### Acknowledgements.

HA is supported by NSFC under Grant No. 11975134, the National Key Research
and Development Program of China under Grant No. 2017YFA0402204 and the
Tsinghua University Initiative Scientific Research Program. MP is supported in
part by U.S. Department of Energy (Grant No. DE-SC0011842). JP is supported by
the New Frontiers Program of the Austrian Academy of Sciences and by the
Austrian Science Fund (FWF) Grant No. FG 1. AR is supported in part by NSERC,
Canada.

## Appendices

##  Appendix A Momentum transfer-dependent screening from electrons

The Boltzmann equation for the distribution of non-relativistic electrons can
be written as

$\displaystyle Df=\frac{\partial f}{\partial t}+\mathbf{v}\cdot\frac{\partial
f}{\partial\mathbf{x}}+\dot{\mathbf{k}}\cdot\frac{\partial
f}{\partial\mathbf{k}}=0\ ,$ (48)

where we neglect collisions. The solution takes the form $f=f^{(0)}\propto
e^{-k^{2}/2m_{e}T}$. We now place the system in a weak and slowly varying
electromagnetic field, so that

$\displaystyle\dot{\mathbf{k}}=e(\mathbf{E}+\mathbf{v}\times\mathbf{B})\ .$
(49)

To study the response, $f$ can be expanded as follows

$\displaystyle f=f^{(0)}+f^{(1)}\ ,$ (50)

so the first order expansion of (48) becomes

$\displaystyle\frac{\partial f^{(1)}}{\partial
t}+\mathbf{v}\cdot\frac{\partial
f^{(1)}}{\partial\mathbf{x}}=-e(\mathbf{E}+\mathbf{v}\times\mathbf{B})\cdot\frac{\partial
f^{(0)}}{\partial\mathbf{k}}=-e\mathbf{E}\cdot\hat{\mathbf{k}}\frac{\partial
f^{(0)}}{\partial k}\ ,$ (51)

where we use the fact that $f^{(0)}$ is isotropic in $\mathbf{k}$. We will
also neglect the contribution from the magnetic field due to the small
velocity. Then the solution of $f^{(1)}$ can be formally written as

$\displaystyle f^{(1)}({\bf x,k},t)=-e\frac{\partial f^{(0)}}{\partial
k}\hat{\bf k}\cdot\int_{-\infty}^{t}dt^{\prime}\,{\bf E}(t^{\prime},{\bf
x}-{\bf v}(t-t^{\prime}))e^{\eta(t^{\prime}-t)}\ ,$ (52)

where $\eta$ is a positive infinitesimal number. The induced electric current
is given by

$\displaystyle{\bf j}_{\rm ind}(t,{\bf x})$ $\displaystyle=$ $\displaystyle
e\int\frac{d^{3}k}{(2\pi)^{3}}{\bf v}f^{(1)}({\bf x},{\bf k},t)$ (53)
$\displaystyle=$
$\displaystyle-e^{2}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{df^{(0)}(k)}{dk}{\bf
v}\int_{-\infty}^{t}dt^{\prime}\,\hat{\bf k}\cdot{\bf E}(t^{\prime},{\bf
x}-{\bf v}(t-t^{\prime}))e^{\eta(t^{\prime}-t)}$ $\displaystyle=$
$\displaystyle-e^{2}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{df^{(0)}(k)}{dk}{\bf
v}\int_{0}^{\infty}d\tau\,\hat{\bf k}\cdot{\bf E}(t-\tau,{\bf x}-{\bf
v}\tau)e^{-\eta\tau}\ .$

We will need the Fourier transform of ${\bf j}_{\rm ind}$ which can be written
as follows,

$\displaystyle{\bf j}_{\rm ind}(q^{0},{\bf q})$ $\displaystyle=$
$\displaystyle\int dt\,d{\bf x}\,e^{iq^{0}t-i{\bf q\cdot x}}{\bf j}_{\rm
ind}(t,{\bf x})$ $\displaystyle=$
$\displaystyle-e^{2}\int_{0}^{\infty}d\tau\int dt~{}e^{iq^{0}t-i{\bf q\cdot
v}\tau}\int\frac{d{\bf k}}{(2\pi)^{3}}\frac{df^{(0)}}{dk}{\bf v}\left[\hat{\bf
k}\cdot{\bf E}(t-\tau,{\bf q})\right]e^{-\eta\tau}$ $\displaystyle=$
$\displaystyle-e^{2}\int d\tau e^{(iq^{0}-i{\bf q\cdot
v}-\eta)\tau}\int\frac{d{\bf k}}{(2\pi)^{3}}\frac{df^{(0)}}{dk}{\bf
v}\left[\hat{\bf k}\cdot{\bf E}(q^{0},{\bf q})\right]$ $\displaystyle=$
$\displaystyle\int\frac{d{\bf k}}{(2\pi)^{3}}\frac{e^{2}}{iq^{0}-i{\bf q\cdot
v}-\eta}\frac{df^{(0)}}{dk}\frac{\bf k}{(k^{2}+m_{e}^{2})^{1/2}}\left[\hat{\bf
k}\cdot{\bf E}(q^{0},{\bf q})\right]$ $\displaystyle\approx$
$\displaystyle\int\frac{d{\bf
k}}{(2\pi)^{3}}\frac{df^{(0)}}{dk}\frac{(-i)e^{2}}{q^{0}-{\bf q}\cdot{\bf
k}/m_{e}+i\eta}\frac{\bf k}{m_{e}}\left[\frac{\bf k}{k}\cdot{\bf E}(q^{0},{\bf
q})\right],$

where we use the non-relativistic approximation in the last step.

The number density of electrons and their phase space distribution function in
the zeroth order are given by

$\displaystyle n_{e}=\int\frac{d{\bf k}}{(2\pi)^{3}}f^{(0)}(k)\
,~{}~{}f^{(0)}(k)=n_{e}\left(\frac{2\pi}{m_{e}T}\right)^{3/2}e^{-k^{2}/2m_{e}T}\
,$ (55)

They determine the polarization,

$\displaystyle{\bf P}=\frac{i}{q^{0}}{\bf j}_{\rm
ind}=-e^{2}n_{e}\frac{(2\pi)^{3/2}}{(m_{e}T)^{5/2}}\int\frac{d{\bf
k}}{(2\pi)^{3}}e^{-k^{2}/2m_{e}T}\frac{{\bf k}[{\bf k}\cdot{\bf E}(q^{0},{\bf
q})]}{m_{e}q^{0}(q^{0}-{\bf q\cdot k}/m_{e}+i\eta)}\ .$ (56)

Therefore, the permittivity tensor can be written as follows,

$\displaystyle\varepsilon^{ij}=\delta^{ij}-\frac{e^{2}n_{e}}{q^{0}}\frac{(2\pi)^{3/2}}{m_{e}^{7/2}T^{5/2}}\int\frac{d{\bf
k}}{(2\pi)^{3}}e^{-k^{2}/2m_{e}T}\frac{k^{i}k^{j}}{(q^{0}-{\bf q\cdot
k}/m_{e}+i\eta)}\ .$ (57)

To determine the Debye screening scale we extract the longitudinal part of
$\varepsilon^{ij}$

$\displaystyle\varepsilon_{L}=\varepsilon^{ij}\frac{q^{i}q^{j}}{q^{2}}=1-\frac{e^{2}n_{e}}{T}\frac{1}{m_{e}^{2}}\left(\frac{2\pi}{m_{e}T}\right)^{3/2}\int\frac{e^{-k^{2}/m_{e}T}k^{2}dk\,d\cos\theta}{(2\pi)^{2}}\frac{k^{2}\cos^{2}\theta}{q^{0}(q^{0}-kq\cos\theta/m_{e}+i\eta)}\
,$

and perform the integration over $\cos\theta$:

$\displaystyle\int
d\cos\theta\frac{\cos^{2}\theta}{q^{0}-kq\cos\theta/m_{e}+i\eta}=\frac{1}{q^{0}+i\eta}\left[\frac{2}{b^{3}}{\rm
tanh}^{-1}[b]-\frac{2}{b^{2}}\right]\ ,$ (59)

where

$\displaystyle b=\frac{kq}{m_{e}(q^{0}+i\eta)}\ .$ (60)

Therefore, we obtain

$\displaystyle\varepsilon_{L}=1+e^{2}n_{e}\left(\frac{2\pi}{m_{e}T}\right)^{3/2}\frac{1}{m_{e}^{2}T}\int\frac{k^{4}dk}{(2\pi)^{2}}\frac{1}{q^{0}(q^{0}+i\eta)}\left[\frac{2}{b^{2}}-\frac{2}{b^{3}}\tanh^{-1}b\right]e^{-k^{2}/2m_{e}T}\
.$ (61)

In the limit of vanishing energy transfer, $q^{0}=0$, we should recover the
standard literature result. It turns out, that the correct procedure is to
take the $\eta\rightarrow 0$ limit first. Then, as $q^{0}\rightarrow 0$, all
$q^{0}$ factors cancel, and we are left with

$\displaystyle\varepsilon_{L}$ $\displaystyle\rightarrow$ $\displaystyle
1+e^{2}n_{e}\left(\frac{2\pi}{m_{e}T}\right)^{3/2}\frac{1}{m_{e}^{2}T}\int\frac{k^{4}dk}{(2\pi)^{2}}e^{-k^{2}/2m_{e}T}2\left(\frac{m_{e}}{kq}\right)^{2}$
(62) $\displaystyle=$ $\displaystyle 1+\frac{e^{2}n_{e}}{Tq^{2}}\ ,$

which agrees with the standard electrostatic formula.

For $q^{0}>0$, we first neglect all the $\eta$’s outside the logarithm, and
obtain

$\displaystyle\varepsilon_{L}$ $\displaystyle=$ $\displaystyle
1+\frac{e^{2}n_{e}}{Tq^{2}}\left[1-2A\int\frac{xdx}{(2\pi)^{1/2}}\tanh^{-1}(x/A)e^{-x^{2}/2}\right]$
(63) $\displaystyle\equiv$ $\displaystyle
1+\frac{e^{2}n_{e}}{Tq^{2}}F_{1}(A),$

where

$\displaystyle A=\frac{q^{0}}{q}\left(\frac{m_{e}}{T}\right)^{1/2}$ (64)

and

$\displaystyle
F_{1}(A)=1-\left(\frac{\pi}{2}\right)^{1/2}Ae^{-A^{2}/2}\left[{\rm
erfi}(A/\sqrt{2})-i\right]\ .$ (65)

The real and imaginary parts of $F_{1}$ are shown in Fig. 12, which indicate
that the correction from non-zero $q^{0}$ can be substantial and of order
unity.

Figure 12: Dispersive (red) and absorptive (blue) parts of the longitudinal
permittivity $\varepsilon_{L}$ as a function of the parameter $A$.

Similarly, we may derive the correction to the transverse mode. We project
$\varepsilon^{ij}$ onto the transverse polarization vector, and exploiting the
fact that the plasma is isotropic and homogeneous, take the average of the two
polarizations. This way we obtain

$\displaystyle\varepsilon_{T}=1-\frac{e^{2}n_{e}}{T}\frac{1}{m_{e}^{2}}\left(\frac{2\pi}{m_{e}T}\right)^{3/2}\int\frac{e^{-k^{2}/m_{e}T}k^{2}dkd\cos\theta}{(2\pi)^{2}}\frac{\frac{1}{2}k^{2}\sin^{2}\theta}{q^{0}(q^{0}-kq\cos\theta/m_{e}+i\eta)}\
.$

On computing the $\theta$-integral, we find

$\displaystyle\varepsilon_{T}$ $\displaystyle=$ $\displaystyle
1-\frac{e^{2}n_{e}}{T}\frac{1}{m_{e}^{2}}\left(\frac{2\pi}{m_{e}T}\right)^{3/2}\int\frac{e^{-k^{2}/m_{e}T}k^{4}dk}{(2\pi)^{2}q^{0}(q^{0}+i\eta)}\frac{(b^{2}-1)\tanh^{-1}[b]+b}{b^{3}}$
(67) $\displaystyle=$ $\displaystyle
1-\frac{e^{2}n_{e}}{Tq^{2}}A\int\frac{xdx}{(2\pi)^{1/2}}e^{-x^{2}/2}\left[\left(\frac{x^{2}}{A^{2}}-1\right)\tanh^{-1}\left(\frac{x}{A}\right)+\frac{x}{A}\right]$
$\displaystyle=$ $\displaystyle
1-\frac{e^{2}n_{e}}{Tq^{2}}\frac{1}{A}\left(\frac{\pi}{2}\right)^{1/2}e^{-A^{2}/2}\left[{\rm
erfi}\left(\frac{A}{\sqrt{2}}\right)-i\right]\ .$

In the non-relativistic limit, for small $q$, $\varepsilon_{T}$ should have a
zero at $(q^{0})^{2}=\omega_{p}^{2}$, where $\omega_{p}^{2}=e^{2}n_{e}/m_{e}$
is the plasma frequency. It is easy to see that this is indeed the case by
taking the limit $A\gg 1$:

$\displaystyle\left(\frac{\pi}{2}\right)^{1/2}\frac{1}{A}e^{-A^{2}/2}{\rm
erfi}(A/\sqrt{2})\rightarrow\frac{1}{A^{2}}\ .$ (68)

Therefore, one finds the expected pole in this limit

$\displaystyle\varepsilon_{T}\rightarrow
1-\frac{e^{2}n_{e}}{m_{e}{q^{0}}^{2}}\ .$ (69)

Similarly, one can find the large $A$ limit of $\varepsilon_{L}$, which
coincides with the expression above.

Finally, we would like to comment on a numerical subtlety in implementing
$F_{1}(A)$ in simulations. The real part of $F_{1}$ contains the error
function with an imaginary argument, which is not defined in ROOT. Therefore,
in order to have a feasible numerical scheme, we use the following empirical
formula,

$\displaystyle{\rm
Re}F_{1}^{\prime}=\left(1-\frac{x^{3}}{5}\right)\exp\left(-\frac{x^{2}}{2}\right)-\frac{1}{x^{2}}\exp\left(-\frac{1}{x^{2}}\right)\
.$ (70)

A comparison between ${\rm Re}F_{1}$ and ${\rm Re}F_{1}^{\prime}$ indicating
their similarity is shown in Fig. 13.

Figure 13: Comparison between ${\rm Re}F_{1}$ and ${\rm Re}F_{1}^{\prime}$.

## References

  * Lee and Weinberg (1977) B. W. Lee and S. Weinberg, Phys. Rev. Lett. 39, 165 (1977).
  * Boehm and Fayet (2004) C. Boehm and P. Fayet, Nucl.Phys. B683, 219 (2004), arXiv:hep-ph/0305261 [hep-ph] .
  * Pospelov  _et al._ (2008) M. Pospelov, A. Ritz, and M. B. Voloshin, Phys.Lett. B662, 53 (2008), arXiv:0711.4866 [hep-ph] .
  * Knapen  _et al._ (2017) S. Knapen, T. Lin, and K. M. Zurek, Phys. Rev. D 96, 115021 (2017), arXiv:1709.07882 [hep-ph] .
  * Battaglieri  _et al._ (2017) M. Battaglieri _et al._ , (2017), arXiv:1707.04591 [hep-ph] .
  * Barak  _et al._ (2020) L. Barak _et al._ (SENSEI), Phys. Rev. Lett. 125, 171802 (2020), arXiv:2004.11378 [astro-ph.CO] .
  * Aguilar-Arevalo  _et al._ (2019) A. Aguilar-Arevalo _et al._ (DAMIC), Phys. Rev. Lett. 123, 181802 (2019), arXiv:1907.12628 [astro-ph.CO] .
  * Agnese  _et al._ (2014) R. Agnese _et al._ (SuperCDMSSoudan Collaboration), Phys.Rev.Lett. 112, 041302 (2014), arXiv:1309.3259 [physics.ins-det] .
  * Aprile  _et al._ (2017) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 119, 181301 (2017), arXiv:1705.06655 [astro-ph.CO] .
  * Akerib  _et al._ (2017a) D. S. Akerib _et al._ (LUX), Phys. Rev. Lett. 118, 021303 (2017a), arXiv:1608.07648 [astro-ph.CO] .
  * Tan  _et al._ (2016) A. Tan _et al._ (PandaX), Phys. Rev. D93, 122009 (2016), arXiv:1602.06563 [hep-ex] .
  * Kouvaris (2015) C. Kouvaris, Phys. Rev. D 92, 075001 (2015), arXiv:1506.04316 [hep-ph] .
  * An  _et al._ (2017) H. An, M. Pospelov, J. Pradler, and A. Ritz, (2017), arXiv:1708.03642 [hep-ph] .
  * Emken  _et al._ (2018) T. Emken, C. Kouvaris, and N. G. Nielsen, Phys. Rev. D 97, 063007 (2018), arXiv:1709.06573 [hep-ph] .
  * Bringmann and Pospelov (2019) T. Bringmann and M. Pospelov, Phys. Rev. Lett. 122, 171801 (2019), arXiv:1810.10543 [hep-ph] .
  * Ema  _et al._ (2019) Y. Ema, F. Sala, and R. Sato, Phys. Rev. Lett. 122, 181802 (2019), arXiv:1811.00520 [hep-ph] .
  * Dent  _et al._ (2020) J. B. Dent, B. Dutta, J. L. Newstead, and I. M. Shoemaker, Phys. Rev. D 101, 116007 (2020), arXiv:1907.03782 [hep-ph] .
  * Cappiello and Beacom (2019) C. Cappiello and J. F. Beacom, Phys. Rev. D 100, 103011 (2019), arXiv:1906.11283 [hep-ph] .
  * Harnik  _et al._ (2021) R. Harnik, R. Plestid, M. Pospelov, and H. Ramani, Phys. Rev. D 103, 075029 (2021), arXiv:2010.11190 [hep-ph] .
  * Ge  _et al._ (2021) S.-F. Ge, J. Liu, Q. Yuan, and N. Zhou, Phys. Rev. Lett. 126, 091804 (2021), arXiv:2005.09480 [hep-ph] .
  * Chen  _et al._ (2020) Y. Chen, J. Shu, X. Xue, G. Yuan, and Q. Yuan, (2020), arXiv:2006.12447 [hep-ph] .
  * Emken (2021) T. Emken, (2021), arXiv:2102.12483 [hep-ph] .
  * Bell  _et al._ (2021) N. F. Bell, J. B. Dent, B. Dutta, S. Ghosh, J. Kumar, J. L. Newstead, and I. M. Shoemaker, (2021), arXiv:2108.00583 [hep-ph] .
  * Bahcall  _et al._ (2005) J. N. Bahcall, A. M. Serenelli, and S. Basu, Astrophys. J. 621, L85 (2005), arXiv:astro-ph/0412440 [astro-ph] .
  * De La Luz  _et al._ (2008) V. De La Luz, A. Lara, E. Mendoza, and M. Shimojo, Geofisica Internacional 47, 197 (2008).
  * Dunsky  _et al._ (2019) D. Dunsky, L. J. Hall, and K. Harigaya, JCAP 07, 015 (2019), arXiv:1812.11116 [astro-ph.HE] .
  * Aprile  _et al._ (2019a) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 123, 251801 (2019a), arXiv:1907.11485 [hep-ex] .
  * Aprile  _et al._ (2020a) E. Aprile _et al._ (XENON), Phys. Rev. D 102, 072004 (2020a), arXiv:2006.09721 [hep-ex] .
  * Essig  _et al._ (2017) R. Essig, T. Volansky, and T.-T. Yu, Phys. Rev. D96, 043017 (2017), arXiv:1703.00910 [hep-ph] .
  * Chang  _et al._ (2021) J. H. Chang, R. Essig, and A. Reinert, JHEP 03, 141 (2021), arXiv:1911.03389 [hep-ph] .
  * Chang  _et al._ (2018) J. H. Chang, R. Essig, and S. D. McDermott, JHEP 09, 051 (2018), arXiv:1803.00993 [hep-ph] .
  * Akerib  _et al._ (2017b) D. S. Akerib _et al._ (LUX), Phys. Rev. D 96, 112011 (2017b), arXiv:1709.00800 [physics.ins-det] .
  * Vogel and Redondo (2014) H. Vogel and J. Redondo, JCAP 1402, 029 (2014), arXiv:1311.2600 [hep-ph] .
  * Dvorkin  _et al._ (2019) C. Dvorkin, T. Lin, and K. Schutz, Phys. Rev. D 99, 115009 (2019), arXiv:1902.08623 [hep-ph] .
  * Essig  _et al._ (2012) R. Essig, J. Mardon, and T. Volansky, Phys. Rev. D85, 076007 (2012), arXiv:1108.5383 [hep-ph] .
  * Aprile  _et al._ (2020b) E. Aprile _et al._ (XENON), (2020b), arXiv:2006.09721 [hep-ex] .
  * Aprile  _et al._ (2019b) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 123, 251801 (2019b), arXiv:1907.11485 [hep-ex] .
  * Essig  _et al._ (2020) R. Essig, J. Pradler, M. Sholapurkar, and T.-T. Yu, Phys. Rev. Lett. 124, 021801 (2020), arXiv:1908.10881 [hep-ph] .
  * Prinz  _et al._ (1998) A. A. Prinz _et al._ , Phys. Rev. Lett. 81, 1175 (1998), arXiv:hep-ex/9804008 [hep-ex] .
  * Crisler  _et al._ (2018) M. Crisler, R. Essig, J. Estrada, G. Fernandez, J. Tiffenberg, M. Sofo haro, T. Volansky, and T.-T. Yu (SENSEI), (2018), arXiv:1804.00088 [hep-ex] .
  * Emken  _et al._ (2019) T. Emken, R. Essig, C. Kouvaris, and M. Sholapurkar, JCAP 09, 070 (2019), arXiv:1905.06348 [hep-ph] .
  * Essig  _et al._ (2016) R. Essig, M. Fernandez-Serra, J. Mardon, A. Soto, T. Volansky, and T.-T. Yu, JHEP 05, 046 (2016), arXiv:1509.01598 [hep-ph] .
  * Buen-Abad  _et al._ (2021) M. A. Buen-Abad, R. Essig, D. McKeen, and Y.-M. Zhong, (2021), arXiv:2107.12377 [astro-ph.CO] .
  * Nguyen  _et al._ (2021) D. Nguyen, D. Sarnaaik, K. K. Boddy, E. O. Nadler, and V. Gluscevic, (2021), arXiv:2107.12380 [astro-ph.CO] .