Constraints from a many-body method on spin-independent dark matter scattering off electrons using data from germanium and xenon detectors

Scattering of light dark matter (LDM) particles with atomic electrons is studied in the context of effective field theory. Contact and long-range interactions between dark matter and an electron are both considered. A state-of-the-art many-body method is used to evaluate the spin-independent atomic ionization cross sections of LDM-electron scattering, with an estimated error about 20%. New upper limits are derived on parameter space spanned by LDM mass and effective coupling strengths using data from the CDMSlite, XENON10, XENON100, and XENON1T experiments. Comparison with existing calculations shows the importance of atomic structure. Two aspects particularly important are relativistic effect for inner-shell ionization and final-state free electron wave function which sensitively depends on the underlying atomic approaches.


I. INTRODUCTION.
Astronomical and cosmological observations not only provide evidences of dark matter (DM) but also point out its properties such as nonrelativistic, non-baryonic, stable with respect to cosmological time scale, and interacting weakly, if any, with the standard model (SM) particles. Its non-gravitational interactions with normal matter are still unknown. A generic class of cold dark matter candidates, the so-called Weakly-Interacting Massive Particles (WIMPs), receive most attention, as they lead to predictions of DM's relic abundance comparable to the measured value and have coupling strengths of weak interaction scales to SM particles, which can be experimentally tested. Also, the existence of such particles are predicted in many extensions of the SM (see, e.g., Refs. [1] for review). Recently there has been remarkable progress made in direct WIMP searches, thanks to novel innovations in detector technologies and increment of detector size and exposure time. As a result, a substantial portion of the favored WIMP parameter space has now been ruled out. For example, the stringent bounds on the spinindependent WIMP-nucleon cross section are currently set by the XENON1T [2]: 4.1 × 10 −47 cm 2 at 30 GeV dark matter mass, 1 , PandaX-II [3]: 8.6 × 10 −47 cm 2 at 40 GeV, and DarkSide-50 [4]: 1 × 10 −41 cm 2 at 1.8 GeV, respectively.
In spite of tremendous efforts in experiment, no concrete evidence of WIMPs has been found to date, di- * jwc@phys.ntu.edu.tw; corresponding author † cpliu@mail.ndhu.edu.tw; corresponding author 1 We use the natural units = c = 1.
rectly or indirectly. This motivates searches of DM particles with masses lighter than generic WIMPs, i.e., 10 GeV/c 2 . Theoretically, such light dark matter (LDM) candidates arise in many well-motivated models, and to account for the relic DM abundance, there are mechanisms suggesting LDM interacts with SM particles through light or heavy mediators with coupling strengths smaller than the weak scale (see Ref. [5] for review). Moreover, annihilations or decays of LDM candidates are possible sources of the anomalous 511 keV [6,7] and 3.5 keV [8,9] emission lines recently found in the sky. Consequently, new ideas to search for LDM flourish and have good discovery potential (see Ref. [10] for general survey).
For energy deposition in the sub-keV region, electron recoil becomes an important subject, no matter being taken as a signal or background, because LDM particles transfer their kinetic energy more efficiently to target electrons than nuclei. Furthermore, electron recoil signals can be used to directly constrain LDMelectron interactions; this complements the study of LDM-nucleon interactions through nuclear recoil and extends a direct detector's scientific reach. Constraints of LDM-electron scattering by direct detection experiments emerged recently, e.g., DAMA/LIBRA [27,28], DarkSide-50 [29], SuperCDMS [25], XENON10 [30,31], XENON100 [31,32], and XENON1T [33]; and much improvement will certainly be expected in next-generation sub-keV detectors.
While electron recoil at sub-keV energies opens a new, exciting window for LDM searches, the scattering processes of LDM particles in detectors pose a fundamental theoretical challenge: The typical energy and momentum of a bound electron is on the order of Z eff m e α and Z 2 eff m e α 2 /2, respectively, where Z eff is the effective nuclear charge felt by an electron of mass m e in a certain shell and α is the fine structure constant with m e α ≈ 3.7 keV. Consequently a sub-keV scattering event strongly overlaps with the atomic scales. This implies a reliable calculation of LDM-electron scattering cross section, which is needed for data analysis, should properly take into account not only the bound nature of atomic electrons but also the electron-electron correlation.
In this work, we applied a state-of-the-art many-body method to evaluate the atomic ionization cross sections of germanium (Ge) and xenon (Xe) by spin-independent LDM-electron scattering. New upper limits on parameter space spanned by the effective coupling strengths and mass of LDM are derived with data from CDMSlite [25], XENON10 [34], XENON100 [35], and XENON1T [33]. The results are also compared with existing calculations.

II. FORMALISM.
A general framework for dark matter interaction with normal matter has recently been developed using effective field theory (EFT). This framework accommodates scalar, fermionic, and vector nonrelativistic (NR) DM particles interacting with NR nucleons via scalar and vector mediators [36]. All leading-order and next-to-leadingorder operators in the effective DM-nucleon interaction are identified [37]. The DM-electron interaction can be formulated similarly with the electron being treated relativistically, as it is essential for atomic structure of Ge and Xe [38]. At leading order (LO), the spin-independent (SI) part is parametrized by two terms: where χ and e denote DM and electron fields, respectively, and q = | q| is the magnitude of 3-momentum transfer, which can be determined by the NR DM particle's energy transfer T and scattering angle θ. The low-energy constants c 1 and d 1 characterize the strengths of the short-range (for heavy mediators) and long-range (for light mediators) interactions, respectively. While the masses of the mediators can vary in broad ranges, it is customary to consider the two extremes: the extremely massive and the massless, which give rise to the contact (or zero-range) and the (infinitely) long-range interactions, respectively. The main scattering process that yields electron recoil is atomic ionization: and the energy deposition by DM is reconstructed by subsequent secondary particles, such as photons and more ionized electrons, recorded in a detector. The differential DM-atom ionization cross section in the laboratory frame through the LO, SI DM-electron interaction is derived in Ref. [38] dσ dT where m χ and v χ are the mass and velocity of the DM particle, respectively. The full information of how the detector atom responds to the incident DM particle is encoded in the response function where |A and |A + , e − denote the many-body initial (bound) and final (ionized) state; M and µ the total and reduced mass of the ion plus free electron system, respectively, with µ ≈ m e . The summation is over all electrons, and the ith electron has its binding energy E Bi , relative coordinate r i , and relative momentum p i . The Dirac delta function imposes energy conservation and constrains the kinematics of the ejected electron, whose energy is NR in the kinematic range of our study but wave function still in the fully relativistic form. Evaluation of R(T, q) is non-trivial. In this work, we use a procedure benchmarked by photoabsorption to a few percent accuracy based on an ab inito method, the (multi-configuration) relativistic random phase approximation, (MC)RRPA [39][40][41][42][43]. Details of how the theory was applied to the responses of Ge and Xe detectors in cases of neutrino scattering are documented in Refs. [44][45][46][47]; here we only give a brief outline and focus on the points that are new in the case of LDM scattering.
First, the ground-state atomic wave functions are calculated by (MC) Dirac-Fock (DF) theory. The multiconfiguration feature is needed for open-shell atoms like Ge, but not for noble gas atoms like Xe. Quality of the initial wave function is benchmarked by the ionization energies of all atomic shells, which can be determined by edges in photoabsorption data.
Second, (MC)RRPA is applied to calculate the transition matrix elements of photoionization. Quality of the final state wave function is benchmarked by how good the calculated photoionization cross section is compared with experiment data. For Ge and Xe, the atomic numbers Z = 32 and 54 are not small, so relativistic corrections to inner-shell electrons or at large 3-momentum transfer are not negligible (a detailed study can be found in Refs. [27,28]). Furthermore, the residual electron correlation is important for excited states, as a result, its (partial) inclusion by RPA makes our calculated photoionization cross sections of Ge and Xe in excellent agreement with experiments. The only exception is T < 80 eV for Ge, where the crystal structure of outer-shell electrons in Ge semiconductor can not be described by our pure atomic calculations [44,46]. Taking the well-benchmarked initial and final state wave functions, response functions for DM-atom scattering, Eq. (4), can in principle be computed in a similar procedure as we previously did for neutrino-atom scattering. However, the high 3-momentum transfer associated with the DM-atom scattering, q m e α, dramatically slows down the (MC)RRPA computation. Therefore, we resorted to a conventional frozen-core approximation (FCA, see, e.g., Refs. [48,49]), so that the finalstate continuum wave function of the ionized electron can be efficiently computed with an electrostatic mean field determined from the ionic state prescribed by (MC)DF. The details of our FCA scheme and its benchmark against (MC)RRPA are given in Appendix A. Except when energy transfer is close to ionization edges, the FCA results generally agree with (MC)RRPA within 20%.
At a direct detector, the measured event rate is where ρ χ = 0.4 GeV/cm 3 is the local DM density [50], and N T is the number of target atoms. The averaged velocity-weighted differential cross section is folded to the conventional Maxwell-Boltzmann velocity distribution f ( v χ ) [51], with escape velocity v esc = 544 km/s [52], circular velocity v 0 = 220 km/s, and averaged Earth relative velocity v E = 232 km/s. The maximum DM velocity seen from the Earth is v max = v esc +v E , and the minimum v min = 2T /m χ is to guarantee enough kinetic energy. This velocity average is timeconsuming because dσ/dT needs to be computed on every grid point of v χ . As commonly seen in literature, e.g., Refs. [53,54], the procedure is simplified by an interchange of integration order: first on v χ of Eq. (6) with [51,55], then the remaining q integration. The details is given in Appendix BB, and we numerically checked that the two average schemes agree very well for all kinematics considered in this work.
In Fig. 1, we show some results of d σv χ /dT for Ge and Xe targets with selected LDM masses. There are several noticeable features. First, the sharp edges correspond to ionization thresholds of specific atomic shells. They clearly indicate the effect of atomic structure, and the peak values sensitively depend on atomic calculations. If direct DM detectors have good enough energy resolution, these peaks can serve as powerful statistical hot spots. Second, away from these edges, the comparison between Ge and Xe cases do point out that the latter has a larger cross section, but the enhancement is not as strong as Z 2 for coherent scattering nor Z for incoherent sum of free electrons. In other words, a heavier target atom does not enjoy much advantage in constraining the SI DM-electron interaction, in opposition to the SI DM-nucleon case. Third, the long-range interaction has a larger inverse energy dependence than the short-range one. As a result, lowering threshold can effectively boost a detector's sensitivity to the long-range DM-electron interaction.  [25], XENON10 (black) [34], XENON100 (red) [35], and XENON1T (magenta and green) [33] data. Superimposed are constraints from Ref. [31] with XENON10 (black-dotted) and XENON100 (red-dotted), and from Ref. [33] with XENON1T (magenta-dotted and green-dotted, short-range only), in which the choices of FDM = 1 (left) and FDM = 1/q 2 (right) correspond to c1-and d1-type interactions, respectively, in this work.

III. RESULTS AND DISCUSSIONS
The CDMSlite experiment with Ge-crystals as target has recently demonstrated the novel mechanism of bolometric amplification [56] and achieved low ionization threshold making it sensitive to LDM searches. A data set of 70.1 kg-day exposure [25] and threshold of 80 eV is adopted for this analysis. The combined trigger and pulse-shape analysis efficiency is more than 80%. The fiducial-volume cut significantly reduced the background and the total combined efficiency including fiducial-volume is adopted from Fig. 4. of Ref. [25]. Limits on LDM-electron scattering are derived without background subtraction with optimum interval method [57]. The derived 90% C.L. limits for both short-and longrange coupling are depicted in Fig. 2.
Dual-phase liquid Xe detectors have demonstrated the sensitivity to ionization of a single electron with their "S2-only" signals [58]. Constraints have been placed in Refs. [30,31] with XENON10 [34] , XENON100 [35], and XENON1T [33] data on LDM-electron scattering using an alternative theoretical framework with different treatment to the atomic physics from this work.
Efficiency-corrected data of XENON10 [30,34] with 15 kg-day exposure, XENON100 [35] with 30 kg-year exposure, and XENON1T [33] with 1 ton-year exposure are extracted from the literature. We follow the same procedure of Refs. [30,31] to convert energy transfer T first to the number of secondary electrons, n e , and then to the photoelectron (PE) yield. Under a conservative assumption that all observed events are from potential LDMelectron scattering, upper limits at 90% C.L. on both short-and long-range interactions are derived and displayed in Fig. 2. The electron recoil charge yield (Q y ) cutoff can change the exclusion region. For the analy-sis of XENON1T data [33], we present both results with and without a cutoff at 12 produced electrons. For the latter, events of smaller ionized electrons enter through Gaussian smearing of PEs.
Comparing the various exclusion curves in Fig. 2 there are several important observations to note. First, the lowest reach of a direct search experiment in LDM mass is determined by its energy threshold. According to what we set for CDMSlite, XENON10, XENON100, and XENON1T: 80, 13.8, 56, and 186 eV, the lightest DM masses can be probed are ∼ 30, 10, 20, and 50 (25 without Q y cutoff) MeV, respectively.
The exclusion limits on DM-electron interaction strengths depend on several factors: experimentally, detector species, energy resolution, background, and exposure mass-time [59]; and theoretically, the DMelectron interaction type and the atomic structure. For the contact interaction, XENON1T yields the best limit when m χ 50 (25 without Q y cutoff) MeV for its overwhelmingly-large exposure mass-time. However, in the lower-mass region, it is the extremely-low threshold of XENON10 that makes it more competitive, despite a much smaller exposure mass-time. The exclusion limit on the long-range interaction is more subtle, because the differential cross section has a sharper energy dependence and tends to weight more at low T . This explains why XENON10's constraint is better than others (less so when XENON1T has no Q y cutoff) in the plot. Among other three experiments, the dominance of XENON1T is shrinking for the same reason. More surprisingly, the finer energy resolution and lower background of CDM-Slite achieves a better constraint than XENON100 when m χ 80 MeV, despite the latter has a big exposure advantage.
In Fig. 2, the exclusion limits derived in Refs. [31,33], using the same xenon data sets, are compared. 2 The differences in the overall exclusion curves are obvious and most likely of theoretical origins. In Fig. 3, we use a sample case to illustrate the difference in predicted event numbers as a function of n e . For both types of interactions, our results are comparatively smaller at small n e but bigger at large n e . This provides a qualitative explanation for the overall differences observed in the exclusion curves: The larger the DM mass m χ , the larger its kinetic energy and hence the increasing chance of higher energy scattering that produces more n e . Therefore, our calculations yield tighter constraints on c 1 for heavier DM particles, but looser for lighter DM particles. As for the long-range interaction, the low-energy cross section is so dominant that the derivation of exclusion limit with the XENON10 data is dictated by the one-electron event, i.e., the first bin. Consequently, the larger event number predicted in Ref. [31] leads to a better constraint on d 1 by a similar size. As for XENON100, its higher energy threshold cut out most sensitivity to low n e events, so both exclusion curves become similar. While there is no exclusion from Ref. [33] to compare, we note that the analysis without Q y cutoff does improve the constraint substantially.
proach, a continuum final state is simply the Coulomb wave function with an effective nuclear charge determined from the binding energy of the atomic shell from which the electron is ionized. In Appendix C, we discuss and compare these atomic approaches in more detail.
As the comparison of FCA vs. NRFCA shows in Fig. 3, for n e ≥ 6, the relativistic effect becomes increasingly important, because the inner shell (in this case, 4d) ionization dominates the cross section. This is consistent with what has been reported in Refs. [27,28]. However, this effect is not enough to explain the discrepancy with Ref. [31]. The difference between the NRFCA and Ref. [31] is most likely due to different formulations of the effective Coulomb potential felt by an ionized electron. However, no further comment can be made as the detail is not explicitly given in Refs. [30,31]. On the other hand, we did find the results of Ref. [31] falls in between NRFCA and HLA, so perhaps is the reconstructed Coulomb potentials.
We note that a later relativistic calculation by Roberts et al. [63] is similar to our framework, but differs in the formulation of the frozen core potential. As a result, even though the agreement of both calculations are generally good, there is still difference due to atomic treatments.

IV. SUMMARY
In summary, we conclude the scattering cross section of sub-GeV dark matter off atoms depends sensitively on atomic structure. Two aspects particularly important are: the relativistic effect, which is sizable when innershell electrons are ionized, and the final-state wave functions on which electron correlation plays an important role. The atomic approached used in our work is fully relativistic, and as the benchmark with (MC)RRPA (a truly many-body approach) validates that our adopted frozen-core mean field is quite robust, and the theoretical uncertainty of our results is estimated to be about 20%. In this appendix, we give an outline of our formulation of the frozen core potential, by which the final ionized electron wave function is computed.
Starting from the self-consistent, relativistic mean field theory, the (MC)DF routine yields a set of single particle orbitals which take a 2-spinor form where a = (n a , κ a ) is a collective label for the principle (n a ) and relativistic orbital (κ a ) quantum numbers of the ath atomic shell; g a and f a the reduced radial wave functions of the large and small components, respectively; and Ω κa,ma (r) the spin-angular function that depends on κ a , the solid angler, and the total magnetic quantum number m a .
In an ionization process where the initial state is unpolarized and the final angular distribution of the free electron is summed, it is more convenient to consider the process isotropic so the frozen core potential has only the radial dependence, and the quantum label a denotes the ionized shell, and the summation over b is performed to all occupied shells. The single-electron potential where r > is the larger one of r and r . The physical meaning of V (a) FCA (r) is clear: In a mean field approach where all electrons are treated as independent particles, the Coulomb potential felt by an ionized a-shell electron is the sum of contributions from the nucleus (of charge Z) and all electrons (the summation b with a degeneracy factor 2j b + 1) except for itself (v a ). As Eq. (A2) looks very similar to the direct term in a typical Dirac-Hatree-Fock approximation, the potential is also called the relativistic Hartree potential (see, e.g., Refs. [48,49]).
We remark that there are other ways of implementing the frozen core approximation. In some approaches, it is even used to replace the two-body Coulomb potential, so the self-consistent mean field equation can be solved more efficiently. When it comes to comparisons of difference atomic approaches, these details can be important.
The procedure we just introduce can be easily implemented in typical non-relativistic Hartree-Fock schemes, which solve the Schrödinger equation instead. The only change in formulae above is the wave function now takes a 1-spinor form, so there is no small component f a , and the quantum label a = (n a , l a ), where l a is the familiar orbital quantum number.
In Fig. 4, the FCA results (in lines) are compared with the (MC)RRPA (in dots). The agreement is generally good at the level of 20%, except for regions close to the ionization edges. Combining Eqs. (3 and 6), the computation of averaged velocity-weighted differential cross section involves a double integration R(T, q)) .

(B1)
The integrand has an apparent v χ dependence which is v χ × v −2 χ = 1 vχ , and an implicit one in q ± , which determine the allowed range of 3-momentum transfer under the given kinematics. Therefore, an interchange of integration order needs to preserve the phase space for it to be exact; or at least the part where most strength lies so it is a good approximation.
A standard procedure of interchanging integration order (see, e.g., Refs. [53,54]) is the following: From energy and momentum conservation, the threshold velocity that a LDM can instigate an energy transfer T and 3momentum transfer q is Fixing the q integration range by equating v min = T /q− + q− /2mχ and q + = 2m χ v max so that q ± no longer have v χ dependence, then the integration of v χ can simply be reduced to a function where the step function Θ(v χ −ṽ min ) imposes the velocity requirement. The analytic functional form of η(ṽ min ) can be found, e.g., in Ref. [51,55]. We have verified within the FCA that the η function approach yields the same result as the more cumbersome double integration of Eq. (B1) for kinematics relevant in this work. This is contrary to the assertion of Ref. [63]. The existing works, including our FCA approach, on DM scattering off atomic electrons are all based on mean field approaches, either non-relativistic [4,30,31,53,54,62] or relativistic [27,28,63]. The first common feature of all is solving the single-electron orbital functions self-consistently from the Hartree-Fock or Dirac-Fock equation. With these orbitals, a good description of the atomic ground state is constructed from one or a linear combination of Slater determinants by filling Z electrons. While relativity does introduce slight changes of the eigenenergies and eigenfunctions of occupied orbitals, it is fair to say that all approaches have very similar atomic initial states to start with, and it is the atomic final states that different approaches diverge.
So far, except for the limited data we carried out with (MC)RRPA, the second common feature of all is the ionized electron state is obtained by solving some mean field equation. Except Refs. [28,63], it simply takes a Schrödinger or Dirac form with an effective, isotropic, electrostatic potential V (a) (r) = Z (a) (r) r , where a is the quantum label of the ionized atomic shell. A few prescriptions of the effective charge Z (a) (r) are listed in the following.
• Plane wave approximation (PWA) [53,54]: This is done by simply taking Z (a) = 0, i.e., the ionized electron is completely free. Calculation-wise, the transition matrix element is simplified to the Fourier transform of the bound state wave function. However, unless the scattering energy is much higher than the atomic scale, this assumption can hardly be justified.
• Hydrogen-like approximation (HLA) [4,54,62]: This approximation assumes the ionized electron behaves like a hydrogen-like electron. The effective charge is simply a unscreened point source at the atomic center with its magnitude fixed by the orbital binding energy E B : where n a is the principle quantum number (in the NR case, n a is an integer; in the relativistic case, there is some correction); and Ry = 13.6 eV. While this approximation is tuned to reproduce the correct binding energy, this point charge picture is not realistic, either.
• Frozen core approximation (FCA): This approximation has been explained in Appendix A. We additionally comment here that the resulting effective charge has the correct asymptotic behaviors: i.e. Z  FCA (r → ∞) → 1. Therefore, it is more realistic than PWA and HLA.
In Fig. 5, we use the xenon 5p shell(s) as an example to illustrate the differences in Z is keV, i.e., the inverse of atomic size, the scattering amplitude depends sensitively on the effective charge distribution. On the other hand, for much smaller or larger q, atomic physics still plays an important role in getting reliable results. For the former, the correct asymptotic requirement Z (a) FCA (r → ∞) → 1 is still crucial (otherwise phase shifts would not be correctly inferred). For the latter, the relativistic and nuclear finite-size effect are both significant, as already pointed out in Refs. [27,28].