Electron-interacting dark matter and implications from DAMA/LIBRA-phase2

We investigate the possibility for the direct detection of low mass (GeV scale) WIMP dark matter in scintillation experiments. Such WIMPs are typically too light to leave appreciable nuclear recoils, but may be detected via their scattering off atomic electrons. In particular, the DAMA Collaboration [R. Bernabei et al., Nucl. Phys. At. Energy 19, 307 (2018)] has recently presented strong evidence of an annual modulation in the scintillation rate observed at energies as low as 1 keV. Despite a strong enhancement in the calculated event rate at low energies, we find that an interpretation in terms of electron-interacting WIMPs cannot be consistent with existing constraints. We also demonstrate the importance of correct treatment of the atomic wavefunctions, and show the resulting event rate is very sensitive to the low-energy performance of the detectors, meaning it is crucial that the detector uncertainties be taken into account. Finally, we demonstrate that the potential scintillation event rate can be much larger than may otherwise be expected, meaning that competitive searches can be performed for m_chi ~ 1 GeV scale WIMPs using the conventional prompt (S1) scintillation signals. This is important given the recent and upcoming very large liquid xenon detectors.


I. INTRODUCTION
The identity and nature of dark matter (DM) remains one of the most important outstanding problems in modern physics. Despite the overwhelming astrophysical evidence for its existence, no conclusive terrestrial observation of DM has yet been reported [1,2]. Currently, most of the effort in the search for DM has focussed on weakly interacting massive particles (WIMPs) with masses 10 − 100 GeV through their hypothesised nongravitational interactions with Standard Model particles. In this work, we consider low-mass WIMP DM with masses on the order of 1 GeV.
One long standing claim of a potential DM detection was made by the DAMA Collaboration, which uses a NaI-based scintillation detector to search for possible DM interactions within the crystal in the underground laboratory at the Gran Sasso National Laboratory, INFN, Italy [3] (see also Refs. [4,5]). The results from the combination of the DAMA/LIBRA and DAMA/NaI experiments indicated an annual modulation in the event rate at around 3 keV electron-equivalent energy deposition (with a low-energy threshold of ∼ 2 keV) with a 9.3σ significance [3]. The phase of this modulation agrees very well with the assumption that the signal is due to the scattering of WIMP DM present in the galactic halo. An annual modulation in the observed WIMP scattering event rate is expected due to the motion of the earth around the sun, which results in an annual variation of the DM flux through a detector (and the mean DM kinetic energy); see, e.g., Refs. [6,7].
Despite the significant signal, there is strong doubt that the DAMA Collaboration result can be due to WIMPs, since it is seemingly in conflict with the null results of many other direct detection experiments, e.g., Refs. [8][9][10][11][12]. There are also several works which offer explanations the DAMA result in terms of non-DM origins, e.g., in Refs. [13,14]. However, it is not always possible to compare different experiments in a model-independent way, meaning it is difficult to make general statements to this effect.
For example, one possibility that has been considered in the literature is that the DAMA modulation signal may be caused by WIMPs that scatter off the atomic electrons [15][16][17], as opposed to nuclear scattering as is assumed in typical experiments. This is particularly applicable for lighter WIMPs ( 10 GeV), which will not leave appreciable nuclear recoils. Most direct detection experiments try to reject pure electron scattering events, in order to perform nuclear recoil searches with as low as possible background. Conversely, the DAMA experiment is sensitive to WIMPs which scatter off either electrons or nuclei, potentially allowing electron-interacting DM to explain the DAMA modulation while avoiding the tight constraints from other experiments. In a recent work [18], however, we used scintillation and ionisation signals from the XENON100 [19] and XENON10 [20] experiments to rule out this possibility for the observed signal above 2 keV; see also Refs. [21][22][23].
Recently, newer results from the DAMA/LIBRA-phase2 experiment have become available [24] (see also [25,26]). These results strengthen the claim for a detected signal, with the significance of the annual modulation in the 2-6 keV energy window rising to 12.9σ. Importantly, the low-energy threshold has been lowered in the new experiment to 1 keV, and the annual modulation is also clearly present in this region (9.5σ significance). This may be of particular significance for the interpretation in terms of electron interacting DM. In our previous work [18] we showed that there would be an almost exponential increase in the potential event rate at lower energies for such models of light (∼ 1 GeV) WIMPs.
For the ∼ keV energy depositions of interest to this work, the relevant process for electron-scattering DM is atomic ionisation. Such processes are kinematically disfavoured at these energy scales, and therefore the scattering probes deep inside the bound-state wave function, with the main contribution coming from wavefunction at distances much smaller than the characteristic Bohr radius of an atom. In such a situation, incorrect smalldistance scaling of the wavefunctions (for example, by using an "effective Z" model, or assuming plane waves for the outgoing ionisation electron) can lead to large errors in the predicted ionisation rates [18]. Further, the relativistic effects for the electron wavefunction are crucial and must be taken into account [27]. As such, interpretation in terms of light WIMPs requires non-trivial calculations of the atomic structure and ionisation rates. Finally, we note that there are several ongoing experiments [28][29][30][31] and proposals [32][33][34][35] to search for light WIMPs in direct detection experiments. We also note that weak evidence for annual modulation at 2 keV from the COSINE collaboration has been recently made public [36] (see also Ref. [37]).

A. Atomic ionisation
Throughout the text and in the figures, we use relativistic units ( = c = 1), with masses, energies and momenta presented in eV, as is standard in the field. However, it is also customary, e.g., to present cross sections in cm 2 , event rates in counts/kg/keV/day. Further, for the calculations of atomic ionisation, it is convenient and common to use atomic units ( = m e = 1, c = 1/α). Therefore, to avoid any possible confusion, we leave all factors and c in the equations.
We consider DM particles that have electron interactions of the form where µ is the inverse of the length-scale for the interaction, set by the mediator mass (e.g., µ = m v c/ ), and α χ is the effective DM-electron coupling strength. Such effective interaction Hamiltonians arise generally in the case of either scalar or vector interactions (e.g., via the exchange of a dark photon). The coefficient in (1) is chosen so that in the case of a massless mediator (long-range interaction, µ = 0), this reduces to a Coulomb(-like) potential (with α → α χ ). In the limit of a very heavy mediator, the above reduces to the contact interaction: V (r) = 4π c(α χ /µ 2 )δ(r). The differential cross section (for fixed velocity v) for the excitation of an electron in the initial state njl is where q is the magnitude of the momentum transfer, E is the energy deposition, and K is the atomic excitation factor, defined [27] Here, f is the density of final states, m = j z , and the total cross section is to be summed over all electrons dσ = dσ njl . The factor of the Hartree energy unit (E H ≡ m e c 2 α 2 27.2 eV) is included in Eqs.
(2) and (3) in order to make the K factor dimensionless (q and µ have dimensions of inverse length). Since we are considering ionisation processes, the final state is a continuum electron with energy ε = E − I njl (I njl is the ionisation energy). Formulas for calculating the atomic excitation factor (3) are given in Appendix A.
Equation (2) is to be integrated over all possible values for the momentum transfer. From conservation of momentum, the allowed values fall in the range between where both the DM particle and ejected electron are assumed to be non-relativistic. For ∼ GeV WIMPs leaving ∼ keV energy depositions, the typical momentum transfer is q ∼ 2m χ E ∼ MeV, which is very large on atomic scales [18]. The resulting differential event rate (per unit mass of target material) is proportional to the cross section (2) averaged over incident DM velocities: where n T is the number density of target atoms (per unit mass), and ρ DM ∼ 0.3 − 0.4 GeV cm −3 is the local DM energy density [38]. We follow Ref. [39] and parameterise the velocity-averaged cross section as whereσ e is the free electron cross section at fixed momentum transfer of q = a −1 0 , a 0 = /(m e cα) is the Bohr radius, α ≈ 1/137 is the fine-structure constant, and f is the DM speed distribution (in the lab frame, normalised to f (v) dv = 1). In the case of a vector or scalar mediated interaction such as (1), these are expressed We have assumed here that m χ m e , which is valid for the considered mass ranges. In the limit of a heavy mediator (contact-like interaction), the DM form factor reduces simply to F χ = 1, while for an ultra-light mediator (Coulomb-like interaction), it reduces to F χ = (a 0 q) −2 . This is a convenient way to parameterise the calculations, and allows for easy model independent comparison between different results. There is no contribution to the event rate stemming from DM velocities below v min = 2E/m χ , the minimum required to deposit energy E. If the majority of the target momentum-space wavefunction density lies inside the (q − , q + ) region, then the integration over q in Eq. (2) is essentially independent of the integration limits, so one may write where η(v min ) is the mean inverse speed of DM particles fast enough to cause the ionisation (v > v min ) for the given velocity distribution. This is a common way to calculate DM direct detection event rates, particularly for nuclear recoils, however, we note that in the case of electron scattering it is not valid. The cross section depends strongly on q − and hence the DM velocity, since, in many cases, the bulk of the electron momentum-space wavefunction lies below the allowed region for momentum transfer. This means that a careful treatment of the DM speed distribution, including uncertainties, is required for the analysis (see also Ref. [44]).

B. Calculation of the atomic ionisation factor
In Fig. 1 we show a comparison of the atomic ionisation factor as calculated using a number of different approximations, as a function of the momentum transfer q (for fixed E). For the values relevant to this work, around q ∼ 1 MeV, there is almost four orders of magnitude difference between the various approximations. Also note that the relativistic effects are very important for large q, and the corrections continue grow with increasing q [27].
Since the typical kinetic energy of a ∼ GeV mass WIMP is large compared to typical atomic transition energies, the minimum momentum transfer is given q min ∼ m χ v χ ∼ E/v χ (4). Therefore, we see that q min with v e ∼ αc − Zαc the typical velocity of an atomic electron. The consequence is that only the very high-momentum tail of the wavefunctions (in momentum space) can contribute to such processes. In position space, this part of the wavefunction comes from distances very close to the nucleus. For a detailed discussion, see Ref. [27].
Therefore, care must be taken to perform the calculations of such processes correctly. For example, it is common to calculate such processes using analytic hydrogen-like wavefunctions, with an effective nuclear charge, which is chosen to reproduce experimental binding energies: Z eff = n 2I njl /E H . While such functions give a reasonable approximation for low q, for the large q values important for this work, they drastically underestimate the cross-section. This is because such functions have incorrect scaling at distances close to the nucleus, which is the only part of the electron wavefunction that can contribute enough momentum transfer.
Another common approach is to approximate the outgoing ionisation electron wavefunction as a plane wave state. Such functions also have the incorrect scaling at small distances, and underestimate the cross section by orders of magnitude for large q. (This is mostly due to the missing Sommerfield enhancement as discussed in Ref. [40].) More details regarding this point are given in Appendix A.
Therefore, to perform accurate calculations one must employ a technique known to accurately reproduce the electron orbitals. Namely, the relativistic Hartree-Fock method, including finite-nuclear size, and using continuum energy eigenstates as the outgoing electron orbitals. Detailed calculations and discussion was presented in Ref. [18]. Formulas are given in the appendix A.
Given the extreme dependence on the atomic physics seen in Fig. 1, it is important to estimate the uncertainty in the calculations. To gauge this, we also calculation the cross-section using other (simpler) methods. Namely, we exclude the effect of the exchange potential from the Hartree-Fock method, and also solve the Dirac equations using only a local parametric potential (chosen to reproduce the ionisation energies) instead of the Hartree-Fock equations. The effect this has on the calculations is very small, with the main difference coming from small changes in the calculated values for the ionisation energies. This is as expected, since the cross-section is due mainly to the value of the wavefunctions on small distances, close to the nucleus, where many-body electron effects are less important (but the correct scaling is crucial). All of these methods (unlike the effective Z method, or plane-wave assumption) give the correct small-r scaling of the bound and continuum electron orbitals.
The finite-nuclear size correction is important for large values of q, but is small compared to the relativistic cor- Velocity-averaged differential cross section for a single Xe atom, withσe = 10 −37 cm 2 . The left panel is for a contact (heavy mediator) interaction, and the right is for a long-range (Coulomb-like) interaction. The kinks in the plots are due to the opening up of deeper electron shells. There is no signal above Emax = mχv 2 max /2, where vmax is the maximum DM velocity.
rections, and ultimately is not a leading source of error. In any case, we include this in an ab initio manner, by directly solving the electron Dirac equation in the field created by the nuclear charge density, which we take to be given by a Fermi distribution, ρ(r) Here t ≡ 4a ln 3 2.3 fm and c 1.1A 1/3 fm are the nuclear skin-thickness and halfdensity radius, respectively, e.g., [41], and ρ 0 is the normalisation factor. We note that the uncertainties stemming from the atomic physics errors are small compared to those coming from the assumed dark matter velocity distribution and detector performance, as discussed in the following sections.
Plots of the velocity averaged differential cross-sections for several WIMP masses and mediator types are presented in Fig. 2. We find very good agreement with similar recent calculations for Xe atoms in Ref. [42]. We present these plots for the xenon atom, since it is the most common target material. For DAMA/LIBRA experiment, the cross section is dominated by scattering off iodine (Z = 53), which has very similar electron structure to xenon (Z = 54).

C. Annual modulation
We assume the DM velocity distribution is described by the standard halo model, with a cut-off (in the galactic rest frame) of v esc = 550(55) km/s, and a circular velocity of v 0 = 220(20) km/s, see, e.g., Ref. [6,43]. The numbers in the parenthesis above represent estimates for the uncertainties in the values. This is important, due to the strong velocity dependence of the cross section (see also, e.g., Ref. [44]). We use these uncertainties to estimate the resulting uncertainty in the calculated event rates.
For the calculations, the velocity distribution is boosted into the Earth frame, which has a speed of v E (t) ≈ v L + v orb cos β cos(2π · yr −1 t + φ).
Here, v L = v 0 + 13 km/s is the average local rest frame velocity, accounting for the peculiar motion of the Sun, v orb is the Earth's orbital velocity, and cos β ≈ 0.49 accounts for the inclination of Earth's orbit to the galactic plane. The cos(ωt) term accounts of the annual change in the local frame velocity due to the orbital motion around the sun, with phase φ chosen such that v E is maximum at June 2, when the Earth and sun velocities add maximally in the galactic halo frame. Due to the strong velocity dependence of the crosssection, the resulting event rates are not perfectly sinusoidal, particularly at higher energies and lower WIMP masses [18]. However, the general sinusoidal feature remains a reasonable approximation. We define the modulation amplitude as (R max − R min )/2.

III. RESULTS
In order to calculate the number of events detected within a particular energy range, the energy resolution of the detectors must be taken into account. To do this, we follow the procedure from Ref. [15], and take the detector resolution to be described by a Gaussian with standard deviation which is measured at low energy to be given by α LE = 0.45 (4), and β LE = 9(5) × 10 −3 [45]. The calculated rate, R, is integrated with the Gaussian profile to determine the observable event rate, S: where E HW is the hardware threshold, which for DAMA is 1 photoelectron [45]. The extracted number of photoelectrons is measured by the DAMA collaboration to be 5.5-7.5 photoelectrons/keV, depending on the detector [45]. We take an average value of 6.5, with ±1 as an FIG. 3. Calculated modulation amplitude for NaI, accounting for the DAMA detector resolution, withσe chosen to reproduce the observed DAMA/LIBRA-phase2 modulation amplitude averaged across the 1-2 keV energy bin. The black points are the combined DAMA/LIBRA data (the points below 2 keV from phase2 alone) [24]. The shaded blue region shows the uncertainties from the calculations (this work), which are mostly due to uncertainties in the DAMA energy resolution. The plot is drawn assuming mχ = 1 GeV for a contact-like interaction (left), and a long-range interaction (right). Clearly, the fit is poor.
error term, so that E HW = 0.15(3) keV. We don't take the detector efficiency into account, because the DAMA collaboration present their results corrected for this [24]. The effect of the finite detector resolution is that it allows events that originally occur at lower energies to be visible in the observed data above the threshold. This is particularly important due to the strong enhancement in the cross section at low energies (see Fig. 2).
Due to the strong atomic number Z dependence, the cross section for scattering off sodium electrons is negligible [18]. So, for the DAMA NaI crystals, it is sufficient to calculate the rate due to scattering of the iodine electrons. We have treated iodine as though it were a free atom, whereas in fact, it is bound in the NaI solid. Only the outermost 5p orbitals are involved in binding. However, even after accounting for the detector resolution, the 5p (and 5s) orbitals contribute negligibly, with the dominant contribution at ∼ 1-2 keV coming from the inner 3s shell, which is very well described by atomic wavefunctions.
Using this approach, we calculate the expected event rate and annual modulation amplitude for DAMA, as a function of the incident WIMP mass, assuming both an ultra-light and super-heavy mediator. Due to the very large enhancement in the expected event rate at smaller energies, the calculated modulation amplitude is a poor fit to the observed DAMA spectrum. In Fig. 3, we present the calculated spectrum along-side the DAMA/LIBRA-phase2 data [24]. For the coupling strength (parameterised in terms ofσ e ), we have fitted the expected event rate to the observed DAMA modulation signal only for the lowest 1 − 2 keV bins. Taking the higher energy bins into account can only increase the best-fit value forσ e , so (as discussed below) this is the most conservative choice.
In Fig. 4 we plot the best-fit regions for the lowest energy DAMA/LIBRA-phase2 modulation signal, as a function of possible DM masses m χ and coupling strengthsσ e . Despite the large enhancement in the ex-pected event rate at the lower energies, and the conservative assumptions made for extracting the best-fit, the interpretation of the observed modulation amplitude in terms of electron-interacting dark matter is inconsistent with existing bounds. All regions of parameter space that could possibly explain the observed DAMA signal are excluded by constraints derived in Ref. [32], using S2 "ionisation-only" results from the XENON10 [20] and XENON100 [46] experiments.
Note the large uncertainties visible in the plots in Figs. 3 and 4. The dominating source of error comes from the uncertainties in the detector response and energy resolution. Sizable errors also arise due to uncertainties in the standard halo model DM velocity distribution. Uncertainties coming from the atomic physics calculations are also included, but are negligible. The uncertainties in the detector resolution and DM velocities themselves are not so large (∼ 10%) -but they lead to very large uncertainties (up to an order of magnitude) in the observable event rate. This is due to the very strong enhancement in the event rate at low energies, which makes the observed rate very sensitive to the detector cut-offs and energy resolution. Clearly, taking these uncertainties into account is crucial.

IV. FUTURE PROSPECTS
We also calculate the potentially observable S1 (prompt scintillation signal) event rate and modulation amplitude for a hypothetical future liquid xenon detector. We model this detector after that of XENON100, and follow Ref. [47] for the conversion from the energy deposition to the observable photoelectron (PE) count (see also Ref. [48]). In this case, the relevant quantity is a counted rate as a function of observable photoelectrons, denoted s1.
The calculated event rate for the production of n pho-  , derived from the S2 ionization signals from the XENON10 and XENON100 experiments, respectively. The "DAMAallowed" regions are excluded for all relevant WIMP masses by these bounds. The fit for DAMA was performed by averaging over just the lowest energy bins without regard to the shape of the spectrum. Taking the higher energy bins into account pushes the DAMA region higher, strengthening this conclusion (see Fig. 3).
FIG. 5. Hypothetically observable WIMP electron recoil event count expected for a 1 tonne·year exposure of a liquid xenon detector (based on XENON100) using the prompt scintillation (S1) signal; (left) for a contact interaction (withσe = 10 −38 cm −1 ), (right) for a contact interaction (withσe = 4 × 10 −35 cm −1 ). Theσe values are chosen to be below the present constraints, which are most stringent around ∼ 0.1 GeV (see Fig. 4). Note that for larger masses, the constraints are less stringent, and larger events rates are not ruled out.
toelectrons is obtained by applying Poisson smearing to the calculated differential rate [47]. We do this according to a Poisson distribution, P n (N ) = e −N (N n /n!), where N = N (E) is the expected/average number of photoelectrons produced for a given energy deposition E, and n is the actual number of photoelectrons produced. The relation between the deposited energy (electron recoil energy) and the produced number of photoelectrons is given in Fig. 2 of Ref. [47]. We model this as a power law: N (E) = aE b , with a = 1.00(25) and a = 1.53 (10), which give the best fit for the lower energies applicable for this work, accounting for the uncertainties from Ref. [47]. Further, to account for the photomultiplier tube (PMT) detector resolution, we convolve the calculated rate with a Gaussian of standard deviation σ = σ PMT √ n, with σ PMT = 0.5 PE [48]. We do not include uncertainty contributions from the PMT resolution, though note that we have checked, and error in the N (E) conversion is by far the dominant source of uncertainty in this step. Finally, the detection acceptance is taken into account as (s1) = 0.88(1 − e s1/3 ) [47], though we note that this has an insignificant impact on the results. The final expression for the observable event rate, S, as a function of counted PEs s1 is dS ds1 We calculate potential event rates, assuming a value forσ e that is not excluded by current experiment, for a one tonne-year exposure in Fig. 5  PE, as in Ref. [22] (see also Refs. [19,21]). This roughly corresponds to the 2 − 6 keV energy window. The rate is strongly dominated by the lower PE contribution, so it doesn't matter where the higher PE cut is taken. We also present the expected rates for the ranges including 1 and 2 PE. The larger and less-well understood background at these lower energies make experiments more difficult to interpret. However, the much enhanced event rate, and large annual modulation amplitudes, may make these regions interesting for future experiments.
In Fig. 6, we show the expected annual modulation fraction for the same type of experiment. Due to the strong velocity dependence of the cross-section, the fractional modulation amplitude is large. For example, for a ∼ 0.1 GeV WIMP, where the event rate may be expected to be high, it is ∼ 15-20%. The peaks in the annual modulation curves (Fig. 6) at around 0.04 and 1 GeV are due to the opening of the n = 3 and n = 4 shells in Xe. Electrons may only become ionised if their binding energies are lower than the maximum kinetic energy of the incident WIMPs: The ionisation rate for shells with energies close to this number (i.e., that are "only just" accessible) will be sensitive to small changes in the velocity distribution. For Xe, these occur for the n = 3 shell just above ∼ 1 keV, and the n = 4 shell just below ∼ 0.1 keV (see Fig. 2).

V. CONCLUSION
We have calculated the expected event rate for atomic ionization by ∼ GeV scale WIMPs that scatter off atomic electrons, relevant to the DAMA/LIBRA direct detection experiment. Though the calculated event rate and annual modulation amplitude is much larger than may be expected, we show that such WIMP models cannot explain the observed DAMA modulation signal without conflicting with existing bounds, even when just the lowest energy 1 − 2 keV bin is fitted. Taking higher bins into account strengthens this conclusion. Further, we demonstrate explicitly the importance of treating the electron wavefunctions correctly, and note that the expected event rates are extremely sensitive to the detector resolution and low-energy performance, and the assumed dark matter velocity distribution. Uncertainties in these quantities lead to large uncertainties in the calculated rates, and therefore must be taken into account. Finally, we calculate the potentially observable event rate for the prompt scintillation signal of future liquid xenon detectors. Large event rates would be expected for dark matter parameters which are not excluded by current experimental bounds, making this an important avenue for potential future discovery.
For the electron wavefunctions, we employ the Dirac basis, in which single-particle orbitals can be expressed where Ω is a spherical spinor, with j 1 m 1 j 2 m 2 |JM a Clebsch-Gordon coefficient, Y lm a spherical harmonic,l = l ± 1 for j = l ± 1/2, and χ σ a spin eigenstate with σ = ±1/2 being the spin orientation. Note, that in the non-relativistic limit, the small component g → 0, and f → P , where P/r is the radial solution to the non-relativistic Schrödinger equation. To reach non-relativistic limit in the calculations, we allow the speed of light c → ∞ inside the code before the Dirac equation is solved. We also note that the relativistic enhancement discussed in the main text does not stem from the lower g component, whose contributions to R scale as (Zα) 2 (except in the case of pseudo-scalar/pseudo-vector interactions, where the g functions contribute at leading order [18]). Instead, they come from differences in the radial dependence of the upper f component [27].  7. The left panel shows the dominating contributions to the atomic ionisation factor K for E = 0.5 keV. For a given energy, K is dominated by the deepest accessible shell (lowest principal quantum number n). The solid black line is the total sum (including states with higher n that are not shown explicitly). The main contribution at low q comes from the states with highest total angular momentum j, while the main large q contribution comes from states with the lowest orbital momentum l. The right panel shows the total K (summed over all accessible atomic electrons) for E = 0.03, 0.5, and 2 keV, which are dominated by the n = 5, 4, and 3 shells, respectively.
with energy normalisation [49], so that:  [50]). Then, from the standard rules for angular momentum, the atomic factor can be expressed as Here, R is the radial integral where ε = I njl − E, j l is a spherical Bessel function, and C is an angular coefficient given by (for closed shells) with [J] ≡ 2J + 1, (. . .) being a Wigner 3j symbol, and Π L ll = 1 if l+l +L is even and is 0 otherwise. The primed quantities refer to the angular momentum state of the ejected ionisation electron (final state). For q 1 MeV, only the L = 0 term contributes significantly, while for q 0.01 MeV, only the L = 1 term is important. For the intermediate region ∼ 0.1 MeV, the sum saturates reasonably rapidly, and convergence is reached by L = 4. These equations are valid for the case of DM that interacts via vector and scalar mediators; similar expressions for the case of pseudo-vector and pseudo-scalar mediators are given in Ref. [18].
Plots of the atomic ionisation factors showing the energy and momentum-transfer dependence, as well as the contributions from different atomic orbitals, are shown in Fig. 7.

Plane wave final states
Here we present the formulas for calculating the ionisation assuming a plane-wave final state. This is done only as a demonstration; we stress, as discussed above, that this is not a reasonable approximation for the processes considered in this work. Take the final ionisation electron state as a plane wave with |p| = √ 2m e ε, subject to the normalisation d 3 p (2π ) 3 p|p = 1.
The relativistic corrections to (A8) are suppressed as ε/(2m e c 2 ), and can be safely excluded. If k ( k ) is the initial (final) WIMP momentum, the minimum allowable momentum transfer can be expressed q − = |k − k| ≈ E/v m e v 2 e /v = p e (v e /v). Since for innershell electrons v e ∼ Zαc, while the DM speed v ∼ 10 −3 c, this implies q p e . Then, the form factor can be expressed as This method of calculating the event rate is used widely in the literature, however, for large values of q it can drastically underestimate the cross section by orders of magnitude (see main text). This is due partly to the missing Sommerfield enhancement, which was also discussed in the context of DM-induced ionisations in Ref. [40]. The effect arises due to the attractive potential of the nucleus, which enhances the value of the unbound electron wavefunction near the nucleus. As discussed above in the main text, it is only the portion of the wavefunction close to the nucleus that contributes to the cross section.
The size of the Sommerfield enhancement can be esti-mated for hydrogen-like s-states as (in atomic units) where each of the K terms is calculated using only the leading small-r terms in the expansion of the wavefunction, and p = √ 2m e ε is the momentum of the outgoing electron. In the non-relativistic limit, the contribution to K coming from the first-order term in the small-r expansion of the electron wavefunctions is identically zero, and the leading non-zero contribution only arises at second order. In contrast, using relativistic functions, one finds that the lowest-order term survives, leading to significant relative enhancement due to relativistic electron effects [27]. Therefore, the non-relativistic equation (A11) also underestimates the relative enhancement. Scaling of inner-shell electron wavefunctions near the nucleus goes as |ψ(0)| 2 ∼ Z 3 , as for hydrogen-like functions. However, the scaling for outer-shell electrons is not as simple [51], so it's important to use electron wavefunctions that correctly reproduce the low-r behaviour, including the correct screening and electron relativistic effects (e.g., the relativistic Hartree-Fock method, as employed here).