Electron-interacting dark matter: Implications from DAMA/LIBRA-phase2 and prospects for liquid xenon detectors and NaI detectors

TL;DR

This paper explores the detection prospects of low-mass WIMP dark matter via electron scattering in scintillation detectors, emphasizing the importance of atomic physics and detector performance, and finds current constraints challenge such interpretations.

Contribution

It highlights the significance of atomic wavefunction treatment and detector uncertainties in interpreting low-energy scintillation signals for GeV-scale WIMPs.

Findings

Electron-interacting WIMPs cannot explain DAMA/LIBRA results within current constraints.

Correct atomic wavefunction modeling significantly affects event rate predictions.

Large liquid xenon detectors can effectively search for GeV-scale WIMPs using scintillation signals.

Abstract

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 GeV scale WIMPs using the conventional prompt S1 scintillation signals. This is important given the recent and upcoming very large liquid xenon detectors.

Figures

Equations

$$V(\mathbf{r})=\hslash c{\alpha }_{\chi }\frac{e^{-\mu r}}{r},$$

$$\frac{\mathrm{d}{\sigma }_{njl}}{\mathrm{d}E}=8\pi {\alpha }_{\chi }^2{\left(\frac{c}{v}\right)}^2{\int }_{q_{-}}^{q_{+}}\frac{q\mathrm{d}q}{(q^2+{\mu }^2{)}^2}\frac{K_{njl}(E,q)}{E_H},$$

$$K_{njl}\equiv E_H\sum _m\sum _f{\left|⟨f|e^{i\mathbf{q}\cdot \mathbf{r}}|njlm⟩\right|}^2{\varrho }_f(E).$$

$$\hslash q_{\pm }=m_{\chi }v\pm \sqrt{m_{\chi }^2{v}^2-2m_{\chi }E},$$

$$dR=\frac{n_T{\rho }_{\mathrm{DM}}}{m_{\chi }{c}^2}\frac{\mathrm{d}⟨{\sigma }_{njl}{v}_{\chi }⟩}{\mathrm{d}E}\mathrm{d}E,$$

$$\frac{⟨\mathrm{d}\sigma v⟩}{\mathrm{d}E}=\frac{\stackrel{ˉ}{}{\sigma }_e}{2m_{e}c}\int \mathrm{d}v\frac{f(v)}{v/c}\underset{q_{-}}{\overset{q_{+}}{\int }}{a}_0^{2}q\mathrm{d}q|F_{\chi }^{\mu }(q){|}^{2}K(E,q),$$

$${\bar{\sigma }}_e$$

$$F_{\chi }(q)$$

$$⟨\mathrm{d}\sigma v⟩\propto \eta (v_{\min })\int q\mathrm{d}q|F_{\mathrm{DM}}(q){|}^{2}K(E,q)\mathrm{d}E,$$

$$v_E(t)\approx v_L+v_{\mathrm{orb}}\cos \beta \cos (2\pi \cdot {\mathrm{yr}}^{-1}t+\varphi ).$$

$${\sigma }_{\mathrm{LE}}/E={\alpha }_{\mathrm{LE}}/\sqrt{E/\mathrm{keV}}+{\beta }_{\mathrm{LE}},$$

$$\frac{\mathrm{d}S}{\mathrm{d}E}={\int }_{E_{\mathrm{HW}}}^{\infty }{g}_{{\sigma }_{\mathrm{LE}}}(E^{\prime }-E)\frac{\mathrm{d}R(E^{\prime })}{\mathrm{d}{E}^{\prime }}\mathrm{d}{E}^{\prime },$$

$$\frac{\mathrm{d}S}{\mathrm{d}s1}=ϵ(s1)\sum _n{g}_{\sigma }(n-s1){\int }_0^{\infty }\frac{\mathrm{d}R(E)}{\mathrm{d}E}{P}_n(N)\mathrm{d}E.$$

$$I_{njl}<K_{\max }\simeq \frac{1}{2}{m}_{\chi }{v}_{\max }^2\sim 4\left(\frac{m_{\chi }{c}^2}{\mathrm{GeV}}\right)\mathrm{keV}.$$

$$⟨\mathbf{r}|njlm⟩=\frac{1}{r}\left(\begin{array}{c}{f}_{njl}(r){\Omega }_{jlm}(\hat{\mathbf{n}})\\ i{g}_{njl}(r){\Omega }_{j\tilde{l}m}(\hat{\mathbf{n}})\end{array}\right),$$

$${\Omega }_{jlm}\equiv \sum _{\sigma }⟨l,m-\sigma ,1/2,\sigma |j,m⟩Y_{l,m-\sigma }(\hat{\mathbf{n}}){\chi }_{\sigma },$$

$$⟨\mathbf{r}|\epsilon jlm⟩=\frac{1}{r}\left(\begin{array}{c}{f}_{\epsilon jl}(r){\Omega }_{jlm}(\hat{\mathbf{n}})\\ i{g}_{\epsilon jl}(r){\Omega }_{j\tilde{l}m}(\hat{\mathbf{n}})\end{array}\right)$$

$${\int }_{\epsilon -\delta \epsilon }^{\epsilon +\delta \epsilon }⟨{\epsilon }^{\prime }jlm|\epsilon jlm⟩\mathrm{d}{\epsilon }^{\prime }=1.$$

$$K_{njl}(E,q)=E_H\sum _L\sum _{j^{\prime }{l}^{\prime }}|R_{njl}^{j^{\prime }{l}^{\prime }L}(\epsilon ,q){|}^2{C}_{jl}^{j^{\prime }{l}^{\prime }L}.$$

$$R_{njl}^{j^{\prime}l^{\prime}L}(\varepsilon,q)=\int_{0}^{\infty}\big{[}f_{njl}(r)f_{\varepsilon j^{\prime}l^{\prime}}(r)\\ +g_{njl}(r)g_{\varepsilon j^{\prime}l^{\prime}}(r)\big{]}j_{L}(qr)\,{\rm d}r,$$

$$C_{jl}^{j^{\prime }{l}^{\prime }L}=[j][j^{\prime }][L]{\left(\begin{array}{ccc}j& j^{\prime }& L\\ -1/2& 1/2& 0\end{array}\right)}^2{\Pi }_{l{l}^{\prime }}^L,$$

$$⟨\mathbf{r}|\mathbf{p}⟩=e^{i\mathbf{p}\cdot \mathbf{r}/\hslash },$$

$$\int \frac{{\mathrm{d}}^3\mathbf{p}}{(2\pi \hslash {)}^3}⟨\mathbf{p}|\mathbf{p}⟩=1.$$

$$K_{njl}^{\mathrm{pw}}=\frac{2E_H}{\pi }|{\Phi }_{njl}(q){|}^2\sqrt{\frac{2m_e^3(E-I_{njl})}{{\hslash }^3}}(2j+1),$$

$${\Phi }_{njl}(q)={\int }_0^{\infty }{f}_{njl}(r)j_l(qr)r\mathrm{d}r.$$

$$\frac{K_{ns_{1/2}}}{K^{\rm pw}_{ns_{1/2}}}\Bigg{|}_{r\to 0}\approx\frac{8\pi Z}{\left[1-\exp(-\frac{2\pi Z}{|p^{\prime}|})\right]n^{3}|p^{\prime}|},$$

Full text

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 Liu et al. (2017); Bertone and Tait (2018). Currently, most of the effort in the search for DM has focussed on weakly interacting massive particles (WIMPs) with masses $\gtrsim 10-100\,$GeV through their hypothesised non-gravitational 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 Bernabei et al. (2013) (see also Refs. Bernabei et al. (2008a, 2012)). 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 $\sim 2$ keV) with a 9.3$\sigma$ significance Bernabei et al. (2013). 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. Freese et al. (2013); Lee et al. (2015).

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. The XENON Collaboration (2019); Aprile and The XENON Collaboration (2018); The LUX Collaboration (2017a); Ren and The PandaX-II Collaboration (2018); Agnes and The DarkSide Collaboration (2018). There are also several works which offer explanations the DAMA result in terms of non-DM origins, e.g., in Ref. Pradler et al. (2013). 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 Bernabei et al. (2008b); Foot (2018); Savage et al. (2009), as opposed to nuclear scattering as is assumed in typical experiments. This is particularly applicable for lighter WIMPs ($\lesssim$ $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 Roberts et al. (2016a), however, we used scintillation and ionisation signals from the XENON100 The XENON Collaboration (2015a) and XENON10 The XENON10 Collaboration (2011) experiments to rule out this possibility for the observed signal above 2 keV; see also Refs. The XENON Collaboration (2017); The XENON Collaboration (2015b); Kopp et al. (2009).

Recently, newer results from the DAMA/LIBRA-phase2 experiment have become available Bernabei et al. (2018) (see also Bernabei and Incicchitti (2017); Bernabei et al. (2015)). 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$\sigma$. 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$\sigma$ significance). This may be of particular significance for the interpretation in terms of electron interacting DM. In our previous work Roberts et al. (2016a) we showed that there would be an almost exponential increase in the potential event rate at lower energies for such models of light ($\sim$ 1 GeV) WIMPs.

For the $\sim$ 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 small-distance 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 Roberts et al. (2016a). Further, the relativistic effects for the electron wavefunction are crucial and must be taken into account Roberts et al. (2016b). 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 Akerib and The LUX Collaboration (2018a, b); The LUX Collaboration (2017b); Kobayashi et al. (2018); Ema et al. (2019); Cappiello and Beacom (2019); Antonello et al. (2019); The COSINE-100 Collaboration (2019a) and proposals Essig et al. (2017); Kouvaris and Pradler (2017); McCabe (2017); Wolf et al. (2019); Emken et al. (2019); Bringmann and Pospelov (2018); Kurinsky et al. (2019) 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 The COSINE-100 Collaboration (2019b) (see also Ref. The Cosine-100 Collaboration (2018)).

II Theory

II.1 Atomic ionisation

Throughout the text and in the figures, we use relativistic units ($\hbar=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 $\,{\rm{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 ($\hbar=m_{e}=1$, $c=1/\alpha$). Therefore, to avoid any possible confusion, we leave all factors $\hbar$ and $c$ in the equations.

We consider DM particles that have electron interactions of the form

[TABLE]

where $\mu$ is the inverse of the length-scale for the interaction, set by the mediator mass (e.g., $\mu=m_{v}c/\hbar$), and $\alpha_{\chi}$ 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, $\mu=0$), this reduces to a Coulomb(-like) potential (with $\alpha\to\alpha_{\chi}$). In the limit of a very heavy mediator, the above reduces to the contact interaction: $V(\boldsymbol{r})=4\pi\hbar c({\alpha_{\chi}}/{\mu^{2}})\delta(\boldsymbol{r})$.

The differential cross section (for fixed velocity $v$) for the excitation of an electron in the initial state $njl$ is

[TABLE]

where $\hbar q$ is the magnitude of the momentum transfer, $E$ is the energy deposition, and $K$ is the atomic excitation factor, defined Roberts et al. (2016b)

[TABLE]

Here, $\varrho_{f}$ is the density of final states, $m=j_{z}$, and the total cross section is to be summed over all electrons ${\rm d}\sigma=\sum{\rm d}\sigma_{njl}$. The factor of the Hartree energy unit ($E_{H}\equiv m_{e}c^{2}\alpha^{2}\simeq 27.2\,{\rm{eV}}$) is included in Eqs. (2) and (3) in order to make the $K$ factor dimensionless ($q$ and $\mu$ have dimensions of inverse length). Since we are considering ionisation processes, the final state is a continuum electron with energy $\varepsilon=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

[TABLE]

where both the DM particle and ejected electron are assumed to be non-relativistic. For $\sim$ GeV WIMPs leaving $\sim$ keV energy depositions, the typical momentum transfer is $q\sim\sqrt{2m_{\chi}E}\sim{\rm MeV}$, which is very large on atomic scales Roberts et al. (2016a).

The resulting differential event rate (per unit mass of target material) is proportional to the cross section (2) averaged over incident DM velocities:

[TABLE]

where $n_{T}$ is the number density of target atoms (per unit mass), and $\rho_{\rm DM}\sim 0.3-0.4\,{\rm{GeV}}\,{\rm{cm}}^{-3}$ is the local DM energy density Bovy and Tremaine (2012). We follow Ref. Essig et al. (2012a) and parameterise the velocity-averaged cross section as

[TABLE]

where $\bar{\sigma}_{e}$ is the free electron cross section at fixed momentum transfer of $q=a_{0}^{-1}$, $a_{0}=\hbar/(m_{e}c\alpha)$ is the Bohr radius, $\alpha\approx 1/137$ is the fine-structure constant, and $f$ is the DM speed distribution (in the lab frame, normalised to $\int f(v)\,{\rm d}v=1$). In the case of a vector or scalar mediated interaction such as (1), these are expressed

[TABLE]

We have assumed here that $m_{\chi}\gg 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_{\chi}=1$, while for an ultra-light mediator (Coulomb-like interaction), it reduces to $F_{\chi}=(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_{\text{min}}=\sqrt{2E/m_{\chi}}$, 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

[TABLE]

where $\eta(v_{\rm min})$ is the mean inverse speed of DM particles fast enough to cause the ionisation ($v>v_{\rm 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. Wu et al. (2019)).

II.2 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\sim 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$ Roberts et al. (2016b).

Since the typical kinetic energy of a $\sim\,$GeV mass WIMP is large compared to typical atomic transition energies, the minimum momentum transfer is given $\hbar q_{\rm min}\sim m_{\chi}v_{\chi}\sim E/v_{\chi}$ (4). Therefore, we see that $\hbar q_{\rm min}\gtrsim m_{e}v_{e}^{2}/v_{\chi}=p_{e}(v_{e}/v_{\chi})\gg p_{e}$, with $v_{e}\sim\alpha c-Z\alpha 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. Roberts et al. (2016b).

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_{\rm eff}=n\sqrt{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 Sommerfeld enhancement as discussed in Ref. Essig et al. (2012b).) 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. Roberts et al. (2016a). 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 corrections, 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, $\rho(r)=\rho_{0}[1+\exp(r-c/a)]^{-1}$. Here $t\equiv 4a\ln 3\simeq 2.3\,{\rm fm}$ and $c\simeq 1.1A^{1/3}\,{\rm fm}$ are the nuclear skin-thickness and half-density radius, respectively, e.g., Fricke et al. (1995), and $\rho_{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. Pandey et al. (2018). 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$).

II.3 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_{\rm esc}=550(55)\,{\rm{km/s}}$, and a circular velocity of $v_{0}=220(20)\,{\rm{km/s}}$, see, e.g., Ref. Freese et al. (2013); Baum et al. (2019). The numbers in the parentheses 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. Wu et al. (2019)). 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

[TABLE]

Here, $v_{L}=v_{0}+13\,{\rm{km/s}}$ is the average local rest frame velocity, accounting for the peculiar motion of the Sun, $v_{\rm orb}$ is the Earth’s orbital velocity, and $\cos\beta\approx 0.49$ accounts for the inclination of Earth’s orbit to the galactic plane. The $\cos(\omega t)$ term accounts of the annual change in the local frame velocity due to the orbital motion around the sun, with phase $\phi$ 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 cross-section, the resulting event rates are not perfectly sinusoidal, particularly at higher energies and lower WIMP masses Roberts et al. (2016a). However, the general sinusoidal feature remains a reasonable approximation. We define the modulation amplitude as $(R_{\rm max}-R_{\rm min})/2$.

III Implications from DAMA/LIBRA-phase2

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. Bernabei et al. (2008b), and take the detector resolution to be described by a Gaussian with standard deviation

[TABLE]

which is measured at low energy to be given by $\alpha_{\rm LE}=0.45(4)$, and $\beta_{\rm LE}=9(5)\times 10^{-3}$ R. Bernabei et al. (2008). The calculated rate, $R$, is integrated with the Gaussian profile to determine the observable event rate, $S$:

[TABLE]

where $E_{\rm HW}$ is the hardware threshold, which for DAMA is 1 photoelectron R. Bernabei et al. (2008). The extracted number of photoelectrons is measured by the DAMA collaboration to be 5.5–7.5 photoelectrons/keV, depending on the detector R. Bernabei et al. (2008). We take an average value of 6.5, with $\pm 1$ as an error term, so that $E_{\rm HW}=0.15(3)\,{\rm{keV}}$. We don’t take the detector efficiency into account, because the DAMA collaboration present their results corrected for this Bernabei et al. (2018).

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 Roberts et al. (2016a). 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 $\sim$ 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 Bernabei et al. (2018). For the coupling strength (parameterised in terms of $\bar{\sigma}_{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 $\bar{\sigma}_{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_{\chi}$ and coupling strengths $\bar{\sigma}_{e}$. Despite the large enhancement in the expected 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. Essig et al. (2017), using $S2$ “ionisation-only” results from the XENON10 The XENON10 Collaboration (2011) and XENON100 The XENON100 Collaboration (2016) 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 ($\sim$ 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 Propsects for liquid xenon detectors

In this section, we discuss the prospects for the detection of light ($\sim$ GeV) WIMPs using xenon duel-phase time projection chambers. (We base our discussion here on XENON Collaboration detectors, see e.g., Ref. Aprile and The XENON100 Collaboration (2012); similar principals apply for other experiments.) When a scattering event occurs in the liquid xenon bulk of such a detector, a prompt $S1$ scintillation signal is induced, which is proportional to the total energy deposited in the detector. Then, any ionized electrons are drifted upwards through the liquid/gas boundary (via an applied electric field), where a secondary scintillation signal ($S2$) that is proportional to the number of ionized electrons may be observed. Combining the $x,y$ spatial resolution of the top and bottom photodetectors with the $z$-resolution from the time between the $S1$ and $S2$ signals allows three-dimensional reconstruction of the event geometry. This allows for the “fiducialization” of the target material, where only scattering events occurring within the inner volume of the detector are included in the analysis. This is an important stage of background rejection, since charged particles are much more likely to scatter quickly, i.e., at the outer regions of the xenon chamber, whereas feebly interacting particles such as WIMPs are equally likely to scatter anywhere within the detector volume. Further, the ratio between the relative strengths of the $S1$ and $S2$ signals may be used to distinguish between nuclear and electronic scattering events. The combination of both the $S1$ and $S2$ signals is thus key to understanding the source of any scattering events.

Proposals to use the ionization-only ($S2$) signals to search for sub-GeV WIMPs have been made previously Essig et al. (2012b, a), and limits from $S2$ observations using XENON10 and XENON100 experiments have been setEssig et al. (2017) (as discussed in the previous section). It is worth noting that the best constraints actually come from the older XENON10 experiment (finished in year 2011), despite its much smaller detector mass, older generation of detectors, and much smaller total exposure. This is due to the detection strategy of the modern experiments, which rely on the combination of $S1$ and $S2$ signals.

The reason $S2$ ionisation-only signals were considered, is because for low-mass WIMPs, the typical energy deposited in the detectors is much smaller than the $\sim$ few keV effective low-energy threshold for $S1$ signals. It was thus believed that the $S1$ scintillation signal produced from such events would be negligible. In this work, we demonstrate that, due to the large enhancement from lower energies and the finite detector resolution, the prompt $S1$ scintillation signal can be many times larger than otherwise expected, and that it therefore can be a promising WIMP direct detection observable. Thus, it would be possible to perform a low-mass WIMP search with modern liquid xenon detectors using the combined $S1$ and $S2$ signals. Detailed calculations of the observable $S2$ spectrum from low mass WIMPs was presented recently in Ref. Essig et al. (2017); here we present calculations for the corresponding $S1$ signal.

We 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. The XENON100 Collaboration (2014a) for the conversion from the energy deposition to the observable photoelectron (PE) count (see also Ref. The XENON100 Collaboration (2014b)). 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$ photoelectrons is obtained by applying Poisson smearing to the calculated differential rate The XENON100 Collaboration (2014a). 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 photo-electrons 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. The XENON100 Collaboration (2014a). 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. The XENON100 Collaboration (2014a).

Further, to account for the photomultiplier tube (PMT) detector resolution, we convolve the calculated rate with a Gaussian of standard deviation $\sigma=\sigma_{\rm PMT}\sqrt{n}$, with $\sigma_{\rm PMT}=0.5\,$PE The XENON100 Collaboration (2014b). 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 $\epsilon(s1)=0.88(1-e^{s1/3})$ The XENON100 Collaboration (2014a), 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

[TABLE]

We calculate potential event rates, assuming a value for $\bar{\sigma}_{e}$ that is not excluded by current experiment, for a one tonne-year exposure in Fig. 5 as a function of the WIMP mass, for both a contact and long range interaction. We show the rate integrated between $3$ and 14 PE, as in Ref. The XENON Collaboration (2015b) (see also Refs. The XENON Collaboration (2017); The XENON Collaboration (2015a)). 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 $\sim$ $0.1\,$GeV WIMP, where the event rate may be expected to be high, it is $\sim$ 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:

[TABLE]

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 $\sim$ $1\,{\rm{keV}}$, and the $n=4$ shell just below $\sim$ $0.1\,{\rm{keV}}$ (see Fig. 2).

V Conclusion

We have calculated the expected event rate for atomic ionization by $\sim$ 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.

Acknowledgements.

BMR gratefully acknowledges financial support from Labex FIRST-TF. This work was also supported by the Australian Research Council and the Gutenberg Research College fellowship.

Appendix A Atomic ionisation factor

A.1 Continuum final states (energy eigenstates)

For the electron wavefunctions, we employ the Dirac basis, in which single-particle orbitals are expressed

[TABLE]

where $\Omega$ is a spherical spinor,

[TABLE]

with $\langle j_{1}m_{1}j_{2}m_{2}|JM\rangle$ a Clebsch-Gordon coefficient, $Y_{lm}$ a spherical harmonic, $\tilde{l}=l\pm 1$ for $j=l\pm{1}/{2}$, and $\chi_{\sigma}$ a spin eigenstate with $\sigma=\pm{1}/{2}$ being the spin orientation. Note, that in the non-relativistic limit, the small component $g\to 0$, and $f\to 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\to\infty$ 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\alpha)^{2}$ (except in the case of pseudo-scalar/pseudo-vector interactions, where the $g$ functions contribute at leading order Roberts et al. (2016a)). Instead, they come from differences in the radial dependence of the upper $f$ component Roberts et al. (2016b).

The continuum state orbitals are defined similarly,

[TABLE]

with energy normalisation Bethe and Salpeter (1977), so that:

[TABLE]

In practice, the normalisation is achieved by a comparison with analytic Coulomb functions at large $r$ Bethe and Salpeter (1977). Note, this formalism means that $\varrho_{f}$ (3) is included already in the definition of the orbitals $\langle\boldsymbol{r}|\varepsilon jlm\rangle$, which have dimension $[{\rm Length}]^{-3/2}[{\rm Energy}]^{-1/2}$.

To calculate the atomic ionisation factor, we first expand the exponential in Eq. (3) as a sum over irreducible spherical tensors: $e^{i\boldsymbol{q}\cdot\boldsymbol{r}}=\sum_{L=0}^{\infty}\sum_{M=-L}^{L}T_{LM}$ (see, e.g., Varshalovich et al. (1988)). Then, from the standard rules for angular momentum, the atomic factor can be expressed as

[TABLE]

Here, $R$ is the radial integral

[TABLE]

where $\varepsilon=I_{njl}-E$, $j_{l}$ is a spherical Bessel function, and $C$ is an angular coefficient given by (for closed shells)

[TABLE]

with $[J]\equiv 2J+1$, $\left(\ldots\right)$ being a Wigner $3j$ symbol, and $\Pi_{ll^{\prime}}^{L}=1$ if $l+l^{\prime}+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\gtrsim 1\,{\rm{MeV}}$, only the $L=0$ term contributes significantly, while for $q\lesssim 0.01\,{\rm{MeV}}$, only the $L=1$ term is important. For the intermediate region $\sim$ $0.1\,{\rm{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. Roberts et al. (2016a).

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.

A.2 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

[TABLE]

with $|p|=\sqrt{2m_{e}\varepsilon}$, subject to the normalisation

[TABLE]

The relativistic corrections to (21) are suppressed as $\sqrt{\varepsilon/(2m_{e}c^{2})}$, and can be safely excluded. If $\hbar\boldsymbol{k}$ ($\hbar\boldsymbol{k^{\prime}}$) is the initial (final) WIMP momentum, the minimum allowable momentum transfer can be expressed $\hbar q_{-}=\hbar|\boldsymbol{k}^{\prime}-\boldsymbol{k}|\approx E/v\gtrsim m_{e}v_{e}^{2}/v=p_{e}(v_{e}/v)$. Since for inner-shell electrons $v_{e}\sim Z\alpha c$, while the DM speed $v\sim 10^{-3}c$, this implies $\hbar q\gg p_{e}$. Then, the form factor can be expressed as

[TABLE]

with the radial integral defined

[TABLE]

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 Sommerfeld enhancement, which was also discussed in the context of DM-induced ionisations in Ref. Essig et al. (2012b). 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 Sommerfeld enhancement can be estimated for hydrogen-like $s$-states as (in atomic units)

[TABLE]

where each of the $K$ terms is calculated using only the leading small-$r$ terms in the expansion of the wavefunction, and $p^{\prime}=\sqrt{2m_{e}\varepsilon}$ 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 Roberts et al. (2016b). Therefore, the non-relativistic equation (24) also underestimates the relative enhancement. Scaling of inner-shell electron wavefunctions near the nucleus goes as $|\psi(0)|^{2}\sim Z^{3}$, as for hydrogen-like functions. However, the scaling for outer-shell electrons is not as simple Landau and Lifshitz (1977), 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).