Directly Detecting MeV-scale Dark Matter via Solar Reflection

If dark matter (DM) particles are lighter than a few MeV/$c^2$ and can scatter off electrons, their interaction within the solar interior results in a considerable hardening of the spectrum of galactic dark matter received on Earth. For a large range of the mass vs. cross section parameter space, $\{m_e, \sigma_e\}$, the"reflected"component of the DM flux is far more energetic than the endpoint of the ambient galactic DM energy distribution, making it detectable with existing DM detectors sensitive to an energy deposition of $10-10^3$ eV. After numerically simulating the small reflected component of the DM flux, we calculate its subsequent signal due to scattering on detector electrons, deriving new constraints on $\sigma_e$ in the MeV and sub-MeV range using existing data from the XENON10/100, LUX, PandaX-II, and XENON1T experiments, as well as making projections for future low threshold direct detection experiments.

Introduction.Astrophysics and cosmology provide one of the strongest arguments for an extension to the Standard Model (SM) of particle physics, through the need for dark matter (DM).The 'theory-space' for dark matter remains vast, motivating a range of experimental approaches.A well-motivated class of models achieve the required relic abundance through thermal freeze-out during the early radiation-dominated epoch, which points to particles with weak-scale interactions -weakly interacting massive particles (WIMPs) -with the required annihilation rate σ ann v ∼ 10 −36 cm 2 (c = 1 from now on).A range of direct detection experiments, searching for the elastic scattering of such DM particles in the galactic halo on nuclei, have now pushed the limit down to the scale of σ n ∼ 10 −46 cm 2 for weak-scale masses [1].
In this Letter, we point out that further direct detection sensitivity to DM in the 10 keV -10 MeV mass range is possible through consideration of 'reflected DM' initially scattered by more energetic electrons in the Sun (or the Earth) prior to scattering in the detector.This being much higher than E halo DM and indeed comparable to the typical solar electron kinetic energy E e ∼ kT e ∼ O(keV).Thus E refl DM can be above the detection threshold for a number of existing experiments, including XENON10, XENON100, LUX, PandaX-II and XENON1T.
The basic scenario is summarized in Fig. 1.DM scattering off free electrons in the Sun generates a new (more energetic) component of the flux impinging on the Earth.While there is necessarily a geometric suppression factor, associated with re-scattering in the direction of the Earth, we find that this is still sufficient to produce new levels of sensitivity to MeV and sub-MeV dark matter, where no direct detection constraints previously existed.The limits and projected sensitivity from electron scattering at a number of experiments are summarized in  [20,23].Filled contours reflect current limits, while dashed contours denote future projections.The thick gray relic density contour is for the DM model in Eq. ( 5).A vertical line at 100 keV indicates a schematic lower limit on m DM from stellar energy loss while the more model-dependent cosmological N eff constraint is not shown (see text).

Fig. 2.
Solar Reflection of Light DM.DM scattering on particles inside the Sun has been extensively studied as an ingredient for the indirect signature of DM annihilation to high energy neutrinos.The evolution of DM that intercepts the Sun depends crucially on its mass.Given a large enough elastic cross section on nuclei, WIMP dark matter with mass above a few GeV can be efficiently captured and thermalized.However, for light DM, the capture process is less efficient, and DM tends to re-scatter at larger radii and evaporate.The 'evaporated' component of the DM flux impinging on the Earth may help improve sensitivity to σ n [24], and, as we are going to show, the effect mediated by σ e is even more pronounced for MeV and sub-MeV mass reflected DM; for a detailed comparison between DM scattering on electrons vs. nucleons inside the sun see [25].
Depending on the scattering cross section σ e , and thus the mean free path, reflection may occur after just one or two interactions, or after partial thermalization through multiple scatters within the Sun.The reflected DM flux will be determined via a simulation which tracks the kinematics after initial entry into the Sun.We will assume a velocity-independent s-wave cross section, but it is notable that the relative importance of the reflected flux would be enhanced for models with a power-like dependence of the cross section on the relative electron-DM velocity, σ e ∝ (v rel ) n , such as would occur e.g. for scattering via higher multipoles.We note in passing that energy-loss or transfer from or inside the sun due to the scattering is negligible for the considered parameter region.
To determine the reflected contribution to the DM flux, the incoming velocity is assumed to follow a Maxwell-Boltzmann distribution with an expectation value of 10 −3 , and an escape velocity cut-off at 2 × 10 −3 .This velocity is negligible compared to solar electrons, and thus DM that scatters in the Sun acquires E recoil DM ∼ T .To gain some intuition, we note first that the probability of scattering off electrons in the solar core is approximately σ e × R core × n core e ∼ σ e /pb, and thus the Sun scatters efficiently if σ e 10 −36 cm 2 .In this optically thick regime, scattering occurs in the convective zone at a characteristic radius R scatt given implicitly by σ e R Rscatt n e (R)dR ∼ O (1).It follows that the electron temperature, and thus the recoil energy, will depend on σ e which in turn determines R scatt , through the radiustemperature relation [26].As the cross section is reduced, R scatt also decreases and E refl,max DM increases as scattering occurs in hotter regions of the core.Further decreasing the cross section ultimately increases the mean free path ∼ (σ e n e ) −1 beyond the solar radius, and the strength of the reflected flux is suppressed.The scattering probability and the background DM flux in the halo, defined through the number density and average velocity as Φ halo ≡ n DM v halo DM , may be combined into a simple estimate for the reflected DM flux incident on the Earth, 2 σ e n core e R core , σ e 1 pb, S g Rscatt 1 A.U. 2 , σ e 1 pb.
(2) In the estimate (2), the overall coefficient of 1/4 has a geometric origin from πR 2 /(4π(1 A.U.) 2 ).S g denotes the gravitational focussing effect that enhances the area at spatial infinity subtended by the effective solar scattering disk πR 2 scatt .For example, at R scatt ∼ R , we have (10), given the value of the solar escape velocity v esc .We note that the overall energy extracted from the Sun by reflected DM does not exceed ∼ 10T × πR 2 Φ halo , and therefore is not constrained by solar energetics being many orders of magnitude below solar luminosity.
Taking a representative choice of m DM ∼ 3 MeV, one can estimate the maximum value of the recoil energy distribution to be ∼ 0.5 T (R scatt ) at σ e 10 −36 cm 2 .For example, a single scatter would accelerate a 3 MeV DM particle up to ∼ 100 eV energy for σ e ∼ 10 −33 cm 2 (R scatt =0.8R ).The reflected flux (2) in this optically thick regime is 10 5 cm −2 s −1 , leading to O(20) ionizations/day in 1kg of Xe.This constitutes a detectable signal, and motivates a more detailed analysis.
Our preliminary estimates (2) need to be augmented to include the possibility of multiple scattering, which can significantly impact the energy of the reflected particles.Since this is difficult to treat analytically, we will make use of a simulation to determine the energy spec- trum and intensity of the reflected DM flux.The simulation scans the initial velocity and impact parameter to determine the initial trajectory into the Sun.The step size was chosen as 0.01R , and the Standard Solar Model [26] was used to determine the temperature, density and elemental abundance at each given radius.For a given cross section σ e , the scattering rate was then determined probabilistically.If DM does not scatter, it propagates to the next step with velocity shifted according to the gravitational potential.If DM scatters, the electron momentum was generated according to the temperature distribution, and the new trajectory determined by first boosting to the DM-electron rest frame, and assuming an s-wave cross section.The gravitational effect on the trajectory was included after each nontrivial scattering.This process was repeated until the DM particle exits the Sun.
We find that it is sufficient to limit our simulations by a maximal impact parameter ρ max = 4R .Outside that range, only the slowest DM particles will enter the Sun, giving a highly subdominant contribution to the reflected flux.Thus, we simulate the energy distribution F Aρ (E) of particles interacting with (or missing) the Sun initially collected from the A ρ = 16πR 2 impact area.After accounting for the gravitational redshift, and the resulting reflected DM flux at Earth determined via As there is some arbitrariness in A ρ , the simulated reflected flux contains an admixture of the initial unscattered distribution.This does not affect subsequent calculations because this component stays below detec-tion thresholds.Fig. 3 shows the final kinetic energy distribution at Earth for 3 MeV DM particles.For σ e ∼ 1 nb, the distribution turns over close to 100 eV, consistent with naive estimates.Moreover, tracking the trajectories indicates that DM does indeed have a higher probability to enter the core region if the cross section is below about 10 −34 cm 2 .Despite the lower cross-section, the enhanced core temperature can in turn lead to less scatters for DM to exit the Sun, resulting in the observed enhancement in the tail of the distribution as the cross-section decreases.However, the effect eventually turns off once the cross section drops well below a pb, as the mean free path and thus the collision rate becomes too low.
Direct detection via electron scattering.With the reflected DM flux and velocity distribution in hand, the scattering signatures can be determined along the lines of the DM-electron scattering analysis of [19,20], with the modifications outlined below.We consider DM scattering off bound electrons in the detector, having fixed energy E e = m e − |E B |, with binding energy E B and a range of momenta.The process of interest corresponds to atomic ionization DM + A → DM + A + + e − with DM threemomentum transfer q.To match the literature, we write the differential scattering rate as a function of electron recoil energy in terms of a reference cross-section σ e [20], where the DM form factor F DM can be taken to 1 if the interaction is short range.We only consider cases where the angular dependence is trivial, q = | q|.The dimensionless atomic form factor describing the strength of the ionization process from atomic state n, l is given by We evaluate the latter using radial Hartree-Fock atomic wavefunctions R nl (r) [27] in ψ nlm ( r) = R nl (r)Y lm (r) and the plane wave approximation | p e = e i p e • r , including a Sommerfeld factor with effective charge Z eff = 1 [19]; p e = 2m e E R,e .When m DM 0.1 MeV, q • r 1 is possible.In order to avoid spurious contributions to f nl from potential numerical non-orthogonality in p e |1|nlm , we subtract the identity operator, and evaluate p e |e i q• r − 1|nlm in these cases instead.The event rate from level (n, l) is then determined by evaluating the average over the incoming energy spectrum of the reflected DM component, that in the nonrelativistic limit is η(E min ) = Emin dE(m DM /(2E)) ) as follows: the number of produced quanta at the interaction point is N Q = E dep./W with W = 13.7 eV [28,29], partitioned into n e ionized electrons escaping the interaction point and n γ scintillation photons.The latter follow a binomial distribution with N Q trials and single event probability f e,γ = n e,γ /N Q .For the purpose of setting limits we only use data above E dep.= 0.19 keV for computing n e , corresponding to the lowest measured charge yield [30] (together with [31]; see also [32,33]), and determine the light output selfconsistently by demanding conservation of N Q .
The detected signals are related by N Q = S1/g 1 + S2/g 2 where g 1 is the light collection efficiency and g 2 is the electron scintillation response times the electron extraction efficiency at the gas-liquid interface.For computing S1 we use the respective values g 1 = 0.12, 0.1134, 0.144, 0.1 PE/γ for XENON100 [34], PandaX-II (run 10) [35], XENON1T [1], and LZ [36].For computing S2 we use the respective values g 2 = 20, 12.1 PE/e − for XENON100 [37] and LUX [38]; for XENON10, the data has already been converted from S2 to the number of electrons [39].S1 is sampled from a binomial distribution with n γ trials and detection probability g 1 ; a Gaussian PMT resolution of σ PMT / ñγ = 0.4 PE in detected photons ñγ is included.For S2 we assume an average 80% electron drift survival probability and apply a representative Gaussian width of σ S2 / √ ñe = 7 PE [40] in the conversion of successfully drifted electrons, ñe , to S2.After accounting for detection efficiencies, and respecting the nominal thresholds in the various experiments, the generated signals are compared to data as reported in [1,3,34,38] and [37,39,41] for S1 and S2-only, respectively.Exemplary spectra for S2 in XENON100 and for S1 in XENON1T are shown in Fig. 4. In the final step, we use the 'p max method' [39,42] to arrive at the limits in the plane of σ e and m DM .
To complete this analysis, we highlight the principal reach of future direct detection experiments (making optimistic assumptions.)For LZ, the next generation liquid xenon experiment [36], we assume, for simplicity, 100% detection efficiency in the acceptance region S1 ≥ 3 PE and include the solar neutrino generated background in the electron recoil band [6].For future semiconductor experiments, we employ the ionization form factor computed in [21] and apply it to a straightforward generalization of (4); we then follow the recommendations of [21] to obtain the projections for SENSEI [43] (superCDMS [44]) with 100 g-yr (10 kg-yr) background-free exposure and 2e − (1e − ) ionization threshold.The results are summarized in Fig. 2. Further details are found in a supplement.
Constraints on Light DM Models.To demonstrate the application of our analysis, we consider a complex scalar dark matter candidate interacting with the electron vector current, ( This model has been analyzed thoroughly, in particular when the interaction is rendered UV-complete via introduction of a kinetically mixed 'dark photon' [11,45,46].The p-wave annihilation channel allows this model to escape stringent CMB constraints [47].Carrying out the standard freeze-out calculation, and adjusting the coupling in σ ann v rel to reproduce the correct relic abundance as a function of m χ , we arrive at the scattering cross section given by , where v 2 e = 1 − m 2 e /m 2 χ .When m χ is close to or below m e , a more accurate thermal average is required, which we implement numerically following Refs.[48,49].The resulting contour is plotted in Fig. 2, and one observes that the reflected DM scattering analysis excludes m χ < 2 MeV region, while higher masses are currently allowed. Going further afield in 'model space', there is now an increased focus on variants of the thermal relic (or WIMP) paradigm, that can ensure the correct relic abundance over the MeV mass range, e.g., SIMPs [50], EL-DERS [51], and models utilizing freeze-in production with very light mediators (so that F DM (q) = (αm e /q) 2 ).The latter case is of interest, as the target parameters correspond to σe ∼ 10 −37 − 10 −38 cm 2 , for m χ ∼ 100 − 1000 keV [21], which provide a challenging goal for future experiments.
Discussion.We have analyzed the direct detection sensitivity to DM-electron scattering, via an energetic 'reflected DM' flux produced through re-scattering in the Sun.This leads to new sensitivity at the sub-pb level for light dark matter in the sub-MeV mass range.Similar rescattering can also occur within the Earth, which would be of particular interest in producing daily modulation.However, the up-scattering effect would be less significant due to the lower electron temperature.
The limits shown in Fig. 2 apply to all DM models with significant scattering cross sections on electrons.However, models in this mass range are subject to a number of powerful indirect constraints.Besides the CMBanisotropy-derived limits on annhilation of DM, there are constraints from stellar energy loss, and the measured radiation energy density, N eff , as well as from primordial nucleosynthesis (BBN) [52][53][54].A universally safe way of escaping the BBN and N eff bounds is to consider m χ > few MeV.Internally thermalized DM models with a lower mass can avoid the constraint on N eff (which in these models is generally shifted below 3), by annihilating into a mixture of SM states (e.g.photons) and neutrinolike dark radiation, as there are compensating effects on the number of equivalent neutrinos [52,53].We emphasize that the new constraints derived on σ e in this paper are direct, and largely independent of additional particle content in the early universe.
We conclude by emphasizing that 'reflected DM' is an intrinsic contribution to the DM flux, and can be probed by all upcoming experiments with sensitivity to electron scattering, e.g.SENSEI, CRESST-III [55], Su-perCDMS, LZ, and CDEX-1T [56].We leave a study of other DM models as well as an investigation of potential signal/background discrimination for future work [57].

MONTE CARLO SIMULATION FOR DM-REFLECTION INSIDE THE SUN
In this section, we provide details on the Monte Carlo program that simulates the reflection of the DM particles from solar electrons.Our program derives both the spectrum of reflected DM and the magnitude of the reflected flux.The simulated flux is then used to derive the expected signal in dark matter detectors along the procedure described in the main text with additional details provided in the subsequent sections.
Our starting point is the standard Maxwell-Boltzmann velocity distribution for galactic DM as it standardly used in direct detection analyses.Irrespective of the details of the assumed galactic velocity distribution, once DM reaches the sun, its velocity will become dominated by the sun's gravity.One of the variables in the simulation is the impact parameter ρ of an incoming DM particle.In the absence of gravitational focussing the relevant range of ρ is limited from above by the solar radius R .The simulation scans the range of impact parameters 0 ≤ ρ ≤ 4R ; larger impact parameters do not change the result as the majority of DM particles in that case misses the Sun.We have numerically checked that the resulting magnitude of the reflected flux and the energy spectrum of the DM coming leaving the Sun do not change if we change the range of the impact parameter from 4R to 3R .The initial conditions of each DM particle are generated by randomly sampling the velocity from the Maxwell-Boltzmann distribution and by choosing the impact parameter evenly from a disc with the radius of 4R .
When a DM particle is outside the Sun, we calculate its trajectory analytically following Newton's law.Using the classical trajectory we determine the incident angle and velocity at the surface of the Sun.Once the DM particle is inside the Sun, for a given cross section, we calculate the mean free path l fp (r) of the DM particle at each radial location r inside the Sun, The input used for this calculation (the electron density inside the Sun n e ) is obtained from the standard solar model of Bahcall [26], with competing solar models such as [59] yielding identical results (see Fig. 5).We then compare 0.1 × l fp with l 0 ≡ 0.01 × R , and choose the smaller to be the step size (l step ) for tracing the trajectory inside the Sun.The probability of DM particle to scatter with an electron within one step is Then we generate a number ξ with a flat random distribution from 0 to 1.If ξ < P c , there is no scattering and we move the DM particle to the end point of this step while changing its velocity according to the gravitational potential.
If ξ > P c the DM particle scatters with an electron within this step.We randomly generate the initial energy and momentum of the electron according to the local temperature of the Sun.Then we boost the DM particle and the electron to their center-of-mass frame.Then following the differential cross section we randomly generate the directions of the outgoing DM particle and the electron.Then we boost them back to the solar frame.(The interactions considered in this paper are contact-type, as the mass of the mediator particle is assumed to be larger than maximum momentum transfer.This simplifies the distribution of the final state momenta in the collision.) At each step we monitor if the DM particle leaves the Sun.Once it is out of the Sun, we calculate its kinetic energy plus the potential energy from the solar gravity and then put it into a histogram.Since the initial impact parameter is from 0 to 4R , there is a chance that the DM particle never passes through the Sun.If this happens we put the initial kinetic energy of this particle into the histogram.Then we normalize the histogram to get a normalized distribution of the energy spectrum of solar reflected DM, which is F Aρ (E) above Eq.( 3) in the paper.
The normalized histograms for m DM = 3 MeV for different cross sections are shown in Fig. 3 in the paper.For larger values of σ e , the DM particle prefers to collide with the electrons in the outer layers of the Sun, and therefore the energy it acquires from the Sun is relatively smaller due to lower temperatures.Whereas in the case of smaller cross sections, the DM particle can penetrate deeper and acquire larger energy through collisions with hotter electrons.This explains why the red curve in the Fig. 3 drops earlier compared to the rest.
A brief summary of the main features of the reflected flux of the DM particles is shown in Tab.I.For each value of the dark matter mass and scattering cross section the average energy of the reflected DM, the endpoint of the reflected DM spectrum (defined as the upper limit of the energy interval containing 95% of the reflected flux), and the total value of the reflected DM flux at the earth position are shown.The endpoint energy in the reflected spectrum should be compared with the galactic endpoint, m DM v 2 esc /2, where v esc is the escape velocity, quantifying the hardening of the reflected spectrum.This ratio is found to be in the range 100-7000 for the parameters listed in the table.

EVALUATION OF THE DM ELECTRON SCATTERING CROSS SECTION
Here we provide the details of our evaluation of the DM(χ)-electron(e) scattering cross section.We treat the electron recoil nonrelativistically, E e = m e + E R,e with E R,e m e , but allow for the general case when the incoming DM particle may be relativistic.From the definition of momentum transfer q = p χ − p χ , energy conservation gives the total amount of energy lost by DM in the collision and hence deposited in the detector, The scattering cross section on bound electrons is conveniently normalized to the non-relativistic cross section σ e on a free electron (see, e.g.[19]), , where the square of scattering amplitude is summed (averaged) over intial (final) state spins.The latter is evaluated at a momentum transfer characteristic for an atomic process, q = | q| = αm e , with any remaining q-dependence and/or χ-energy dependence absorbed by a DM form factor For the purpose of this paper, we consider the simple case of a contact interaction for which F DM = 1 while more general cases are obtained in a straightforward manner.The differential electron recoil rate resulting from ionization of atomic state with principal and angular quantum numbers n and l can be brought into the form, where ∆E e = E R,e + |E B |.The δ-function can be used to perform the angular part of the integral in A factor 1/v in the usual expression of galactic DM-electron scattering is being replaced by the more general factor m 2 χ /(p χ E χ ).The ionization form factor is given by, where p ≡ p e − q.For very light DM, | q • r| 1 may be attained in the evaluation of (11).If the final state wave function is approximated by a plane or Coulomb wave, and not the exact atomic electron wave function, the bound and scattering states are not guaranteed to be orthogonal, nlm| p e = 0. To avoid spurious contributions in the approximations employed we modify the transition matrix element by subtracting the unity operator, where χ nl is the fourier transform of the radial Hartree-Fock wave function R nl , χ nl (q) = 4π(−1) l i l dr r 2 R nl (r)j l (qr).Explicitly, one finds, The first term on the right hand side is the DM-electron scattering form factor previously obtained in the literature [19] while the subsequent terms originate from the subtraction.The remaining angular integral is evaluated as, d cos θ χ nl p 2 e + q 2 − 2p e q cos θ P l p e − q cos θ p 2 e + q 2 − 2p e q cos θ .
Since two unit vectors are dotted as argument of the Legendre polynomial P l , the latter integrand is always wellbehaved.One can verify that the subtraction works by noting that, (q → 0), (15) and that P l (1) = 1 so that Eq. ( 13) indeed vanishes identically for q → 0. We have also verified this numerically in our computer code.It turns out that the subtraction has only a relatively minor influence on the derived bounds, see Fig. 5.We note in passing that the case considered here is quite different from the case when DM is slow and | q • r| 1 for which the ionization process becomes short-distance dominated, and a relativistic atomic treatment becomes necessary [60].
The minimum incoming momentum to produce an electron recoil E R,e is obtained from the δ-function in (9) when q and p χ are parallel, i.e. cos θ qpχ = 1, This expression is exact and is used in the integral that computes the average over the incoming energy spectrum (see main text).In the non-relativistic limit, expanding in ∆E e /q 1 one recovers the expression for v min = p min χ /m χ given in previous works [19].

MODELING OF LXE DETECTOR RESPONSE
Here we describe a simple procedure that is aimed at capturing the dominating factors in the detection of scintillation (S1) and ionization (S2) signals following an electron recoil in a LXE detector.Our treatment largely follows previous experimental, theoretical, and joint theory-experiment studies [6,61,62].With the results of the paper together with the details of the MC simulation provided in the previous section of the supplement, we invite the experimental collaborations to perform their own dedicated analysis of the reflected dark matter signal.
Given a total energy deposition E dep.= E χ − E χ from DM scattering on an atomic electron, the average number of produced quanta at the interaction point is, with W = 13.7 eV [28,29].The quanta are partitioned into n e ionized electrons escaping the interaction point and n γ scintillation photons; Q y and L y denote the energy-dependent charge and light yields, respectively.For the purpose of setting limits we only use data above E dep.= 0.19 keV in the computation of n e , corresponding to the lowest energy at which Q y was measured [30].[63] In addition, data on Q y from [31] is used; see also [32,33].The light output is then obtained self-consistently by demanding conservation of N Q , For E dep. 10 keV the charge and scintillation yields depend only mildly on the applied drift voltage [32,64,65]; for given E dep., Q y (L y ) varies by about 10% in the range 120 − 730 V/cm, i.e. in the range of applied drift-fields accross the various experiments.We hence neglect this experiment-specific detail as it is expected to have only relatively minor impact on the resulting bounds.Finally, heat losses lead to a quenching in the nuclear recoil signal, and subsequently to a fluctuation in N Q .For electron recoils, heat losses are negligible [32,64] and we take the number of produced quanta for a given deposited energy as a constant, N Q = N Q .

Fluctuations at the interaction point
Recombination at the interaction point shifts the partition of n γ and n e while holding their sum N Q = n γ +n e fixed.The quantities n γ and n e are the end products after an initial number of ions n i (yielding electrons) and excitons n ex (yielding scintillation) had been created but were redistributed because of recombination described by parameter r, n e = n i (1 − r), n γ = n i (r + α), with α ≡ n ex /n i .Note that N Q = n γ + n e = n i + n ex .Fluctuations in r itself, leading to an observed variance that is in excess from one that is expected from a binomial process, are most important for larger energy depositions E R 2 keV [29,31]; see also [65].Since the bulk of the solar flux is energetically lower, we neglect this complication.Thus, we follow [62] and approximate the primary signal formation by directly computing the probabilities of producing either n γ or n e as, where,

Detector-specific fluctuations
The S1 signal is obtained from the probability of detecting n PE photons from n γ produced, and is modelled by a binomial distribution with overall light collection efficiency g 1 (the experiment-specific values of g 1 are found in the main text), For computing the S2 signal, one needs to account for the survival probability p surv of electrons when they are drifted by a distance ∆z until the liquid-gas interface with a ballpark velocity v d ∼ 1.7mm/µ sec [61,66], The electron lifetime varies across experiments.We follow the reasoning presented in [62] and assume that the electron lifetime is distributed uniformly over [0, 2/3] mm/µs.The probability of n surv e electrons reaching the gas-phase of the detector is then found from the compound distribution function, where one marginalizes over the production location, To relax the numerical demand, in the actual analysis we assign an average electron survival probability P (n surv e |n e ) p surv = 0.8 obtained from the expectation value of (21) which then allows to perform the sum in (27) below.
In a final step we account for the PMT resolution, using as a representative value σ PMT 0.4 √ n PE for the lowphoton count in S1.For S2, once the electron reaches the gas phase it receives a gain factor g 2 in the number of produced photo-electrons that are detected (for the experiment-specific values of g 2 see main text).( Such expression for the joint PDF in S1 and S2 makes it amenable to a likelihood analysis of the experimental data.
For the purpose of this paper, where we primarily explore the principal experimental sensitivity and do not intend to forestall a dedicated experimental analysis, we consider the signals S1 and S2 separately, leaving an analysis in the correlated signal for more in-depth work.The pdf for S1 is then obtained from, where in the last equality we have performed the sum over n γ .For obtaining the PDF in S2 only, one evaluates, Of course, the PDFs ( 26) and ( 27) also follow directly from (25) by marginalizing over S2 and S1, respectively.

DETAILS ON THE EXEMPLARY DM MODEL
The nonrelativistic scattering cross section on free electrons obtained from the Lagrangian given in the main text is found to be where µ is the reduced mass.In turn, the total annihliation cross section to electrons in terms of the squared center-of-mass energy s reads, In the non-relativistic expansion, s = 4m 2 χ + m 2 χ v, and one finds the velocity scaling corresponding to p-wave annihilation, We have calculated the relic abundance with three different methods: semi-analytically following [67], numerically following [48], and through an implementation of the publicly available computer code Micromegas [49] with mutually agreeing results.Imposing the relic density condition, allows one to obtain the expectation for σ e as a function of m χ as provided in the main text.Within this model, the area below the line is disfavored by DM overproduction, while the area above the line makes χ to form only a fraction of the total DM energy density.

FIG. 1 .
FIG. 1.A schematic illustration of the reflected dark matter flux generated through solar scattering.For bound solar electrons with energy Ee ∼ kT , the DM recoil energy is bounded by the expression in Eq. (1) and can be ∼ keV.

38 FIG. 2 .
FIG. 2. Exclusion contours for 'reflected DM' from a rangeof experiments are shown in comparison to previous limits from XENON10 and XENON100 on scattering from the 'galactic DM' halo population[20,23].Filled contours reflect current limits, while dashed contours denote future projections.The thick gray relic density contour is for the DM model in Eq. (5).A vertical line at 100 keV indicates a schematic lower limit on m DM from stellar energy loss while the more model-dependent cosmological N eff constraint is not shown (see text).

FIG. 3 .
FIG. 3. Normalized energy distributionsF Aρ=16πR 2 (E) (in eV ),are shown for reflected DM with a mass of 3 MeV and the range of scattering cross sections indicated.The initial velocity is assumed to follow a Maxwell-Boltzmann distribution with an expectation value of 10 −3 , and an escape velocity cut-off at 2 × 10 −3 .It is apparent that the distributions below 5-7 eV tend to that of the background halo.

m χ = 1 FIG. 4 .
FIG. 4.Exemplary electron scattering event rates as a function of S2 in XENON100 (upper panel) and as a function of S1 for XENON1T (lower panel).When setting limits we require a minimum deposited energy of E dep.> 0.19 keV (dashed curve).

37 FIG. 5 .
FIG.5.Impact of various systematic theoretical uncertainties on the direct detection limits originating from the solar reflection process as well as from the evaluation of the atomic scattering section.
photon count in S1.For S2, once the electron reaches the gas phase it receives a gain factor g 2 in the number of produced photo-electrons that are detected (for the experiment-specific values of g 2 see main text).The process is modeled by a Gaussian with representative standard deviation σ S2 = 7√ n surv e[40].The respective probabilities of detection read, P (S1|n PE ) = gauss(S1|n PE , σ PMT ), (23)P (S2|n surv e ) = gauss(S2|g 2 n e surv , σ S2 ).(24)Collecting all factors allows one to estimate the correlated PDF for observing S1 and S2, given an energy deposition of E dep., P (S1, S2|E R ) = e |n e )P (n e | n e ).
1/2 (dΦ refl /dE)Φ −1 halo .Multiplying it by the flux and target density N T , we arrive at the total rate from the (n, l) state, dR nl /d ln E R,e = N T Φ halo d σ nl v /d ln E R,e , where E min is the minimum DM energy required to produce an electron with E R,e recoil energy.The resulting electron recoil energy spectrum is converted into scintillation (S1) and ionization (S2) responses in liquid xenon experiments, dR nl /dSi = ε(Si) dE R,e pdf(Si|E nl dep.)dR/d ln E R,e .Here, ε(Si) is the Si detection efficiency and pdf(Si|E dep. ) is the probability to produce Si given a deposited energy E nl dep.= E R,e + |E nl B |.For the purpose of this work, we consider the signals in S1 and S2 separately, and model pdf(Si|E nl dep.