Electron Ionization via Dark Matter-Electron Scattering and the Migdal Effect

TL;DR

This paper compares electron ionization signals from direct DM-electron scattering and the Migdal effect, showing the latter can dominate for heavy mediators and certain DM masses, impacting experimental interpretations.

Contribution

It provides a unified theoretical framework for both ionization processes and clarifies when the Migdal effect becomes significant in dark matter detection.

Findings

Migdal effect scales with Z^2 and q^2, making it more relevant for heavy mediators.

For light mediators, electron scattering dominates ionization signals.

Heavy mediators and DM masses in hundreds of MeV can lead to Migdal-induced ionization surpassing direct electron scattering.

Abstract

There are currently several existing and proposed experiments designed to probe sub-GeV dark matter (DM) using electron ionization in various materials. The projected signal rates for these experiments assume that this ionization yield arises only from DM scattering directly off electron targets, ignoring secondary ionization contributions from DM scattering off nuclear targets. We investigate the validity of this assumption and show that if sub-GeV DM couples with comparable strength to both protons and electrons, as would be the case for a dark photon mediator, the ionization signal from atomic scattering via the Migdal effect scales with the atomic number $Z$ and 3-momentum transfer $\mathbf{q}$ as $Z^2 \mathbf{q}^2$. The result is that the Migdal effect is always subdominant to electron scattering when the mediator is light, but that Migdal-induced ionization can dominate over electron scattering for heavy mediators and DM masses in the hundreds of MeV range. We put these two ionization processes on identical theoretical footing, address some theoretical uncertainties in the choice of atomic wavefunctions used to compute rates, and discuss the implications for DM scenarios where the Migdal process dominates, including for XENON10, XENON100, and the recent XENON1T results on light DM scattering.

Tables

Table 1. Table 1: Comparison of the ionization energy E b subscript 𝐸 𝑏 E_{b} and electron binding energies | E n ​ l | subscript 𝐸 𝑛 𝑙 |E_{nl}| (eV) of the 5p, 5s, and 4d electron shells of xenon from calculations using the formalisms of Refs.  Ibe  et al. ( 2018 ); Essig  et al. ( 2017b ) . The measured values  Brandi  et al. ( 2001 ); Cardona and Ley ( 1978 ) tend to fall somewhere in the middle.

Method		$E_b$		$\left|E_{5p}\right|$	$\left|E_{5s}\right|$	$\left|E_{4d}\right|$
Ibe et al.  Ibe et al. (2018)		9.8		9.8	21	61
Essig et al. Essig et al. (2017b); Bunge et al. (1993)		12.4		12.4	25.7	75.6
Measured Brandi et al. (2001); Cardona and Ley (1978)		12.1		12.8	23.3	68.5

Equations

$$\frac{d{R}_M/d\mathbf{q}}{d{R}_e/d\mathbf{q}}>Z^2{\left(\frac{m_e}{m_N}\right)}^2(\mathbf{q}{r}_a{)}^2,$$

$$\mathcal{L}=-\frac{1}{4}{F}_{\mu \nu }^{\prime }{F}^{\prime \mu \nu }+\frac{m_{A^{\prime }}^2}{2}{A}_{\mu }^{\prime }{A}^{\prime \mu }+A_{\mu }^{\prime }\left(g_D{J}_D^{\mu }+ϵe{J}_{\mathrm{EM}}^{\mu }\right),$$

$${\overline{\sigma }}_T=\frac{16\pi {ϵ}^2\alpha {\alpha }_D{\mu }_{\chi T}^2}{(m_{A^{\prime }}^2+|{\mathbf{q}}_0{|}^2{)}^2},$$

$$\frac{d{R}_M/d\mathbf{q}}{d{R}_e/d\mathbf{q}}\gtrsim {\left(\frac{Z{g}_p+(A-Z)g_n}{g_e}\right)}^2{\left(\frac{m_e}{m_N}\right)}^2(\mathbf{q}{r}_a{)}^2,$$

$$\mathbf{q}$$

$$\Delta E_e=\frac{1}{2}{m}_{\chi }{v}^2-\frac{|m_{\chi }\mathbf{v}-\mathbf{q}{|}^2}{2m_{\chi }}-\frac{{\mathbf{q}}^2}{2m_N}=\mathbf{q}\cdot \mathbf{v}-\frac{{\mathbf{q}}^2}{2{\mu }_{\chi N}}$$

$$H_{\mathrm{int},e}=-\int \frac{d^3\mathbf{q}}{(2\pi {)}^3}{e}^{i\mathbf{q}\cdot (-{\mathbf{x}}_{\chi }+{\sum }_{i=1}^Z{\mathbf{x}}_i)}\frac{{\mathcal{M}}_{e\chi }(q)}{4m_{\chi }{m}_e},$$

$$R_e\propto |⟨{\Psi }_f|H_{\mathrm{int}}|{\Psi }_i⟩{|}^2\sim |⟨{\psi }_f|e^{i\mathbf{q}\cdot \mathbf{x}}|{\psi }_i⟩{|}^2,$$

$$H_{\mathrm{int},N}=-\int \frac{d^3\mathbf{q}}{(2\pi {)}^3}{e}^{i\mathbf{q}\cdot ({\mathbf{x}}_N-{\mathbf{x}}_{\chi })}\frac{{\mathcal{M}}_{N\chi }(q)}{4m_{\chi }{m}_N},$$

$$R_M\propto |⟨{\Psi }_{{\mathbf{v}}_A}|{\Psi }_i⟩{|}^2\sim |⟨{\psi }_f|e^{i{\mathbf{q}}_e\cdot \mathbf{x}}|{\psi }_i⟩{|}^2,$$

$${\mathbf{q}}_e=\frac{m_e}{m_N}\mathbf{q}.$$

$$\frac{d{R}_M/d\mathbf{q}}{d{R}_e/d\mathbf{q}}>{\mathbf{q}}^2{r}_a^2{\left(\frac{m_e}{m_N}\right)}^2\frac{|{\mathcal{M}}_{\chi N}(q)/m_N{|}^2}{|{\mathcal{M}}_{\chi e}(q)/m_e{|}^2},$$

$$R_{i\to f}=\frac{{\rho }_{\chi }}{m_{\chi }}\int d^{3}vf(v)\sigma v_{i\to f},$$

$$\sum _f=\frac{1}{2}\sum _{l^{\prime }{m}^{\prime }}\int \frac{k^{\prime 3}d\ln E_e}{(2\pi {)}^3},$$

$$\frac{{\overline{\sigma }}_e}{{\mu }_{\chi e}^2}=\frac{{\overline{\sigma }}_p}{{\mu }_{\chi p}^2},$$

$$|\mathcal{M}(q){|}^2=\frac{16\pi m_{\chi }^2{m}_T^2{\overline{\sigma }}_T}{{\mu }_{\chi T}^2}|F_{\mathrm{DM}}(q){|}^2,$$

$$\frac{|{\mathcal{M}}_{\chi N}(q){|}^2}{m_N^2}=Z^2|F_N(q){|}^2\frac{|{\mathcal{M}}_{\chi e}(q){|}^2}{m_e^2}$$

$$\frac{d{R}_{e,M}}{d\ln E_e}=N_T\frac{{\rho }_{\chi }}{m_{\chi }}\frac{\overline{\sigma }}{8{\mu }_{\chi e}^2}{I}_{e,M}(E_e),$$

$$I_{e,M}(E_e)=\int d|\mathbf{q}||\mathbf{q}||F_{\mathrm{DM}}(q){|}^2\eta (v_{\min })|f_{e,M}(E_e,\mathbf{q}){|}^2.$$

$$v_{\min }=\frac{\Delta E_e}{|\mathbf{q}|}+\frac{|\mathbf{q}|}{2{\mu }_{\chi N}}=\frac{|E_{nl}|+E_e}{|\mathbf{q}|}+\frac{|\mathbf{q}|}{2{\mu }_{\chi N}},$$

$$|f_e(E_e,\mathbf{q}){|}^2=\frac{k^{\prime 3}}{4{\pi }^3}\times 2\sum _{n,l,l^{\prime },m^{\prime }}|⟨{\psi }_{E_e}^f|e^{i\mathbf{q}\cdot \mathbf{x}}|{\psi }_{E_{nl}}^i⟩{|}^2,$$

$$⟨{\psi }_{E_1,l,m}|{\psi }_{E_2,l^{\prime },m^{\prime }}⟩=\frac{(2\pi {)}^3}{k_1^2}\delta (k_1-k_2){\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }},$$

$$|f_M(E_e,\mathbf{q}){|}^2$$

$$\times 2\sum _{n,l,l^{\prime },m^{\prime }}|⟨{\psi }_{E_e}^f|e^{i{\mathbf{q}}_e\cdot \mathbf{x}}|{\psi }_{E_{nl}}^i⟩{|}^2.$$

$$|\mathbf{q}{|}_{\min }=\frac{E_b}{v_{\max }},$$

$$|\mathbf{q}{|}_{\max }=2{\mu }_{\chi N}{v}_{\max }\sim 5\mathrm{keV}\left(\frac{m_{\chi }}{\mathrm{MeV}}\right),$$

$$f_0=\frac{1-f_R}{1+(N_{\star }/N_i)}\approx 0.83,$$

$$n_t=\mathrm{floor}\left(\frac{E_e}{W}\right)+\mathrm{floor}\left(\frac{|E_{nl}|-E_b}{W}\right),$$

$$n_e=(1-f_R)+\left(\genfrac{}{}{0ex}{}{n_t}{f_0}\right),$$

$$E_{nl}=\frac{1}{2}\sum _{\kappa }{\delta }_{l,|\kappa +1/2|-1/2}{E}_{n\kappa },$$

Full text

††preprint: FERMILAB-PUB-19-257-A

I Introduction

Although the evidence for dark matter (DM) is overwhelming, its microscopic properties remain unknown and motivate various experimental techniques to uncover its possible non-gravitational interactions Bertone and Hooper (2018). In recent years, there have been several novel experimental techniques introduced to directly detect hitherto inaccessible DM candidates below the GeV scale Battaglieri et al. (2017). One particularly promising strategy involves detecting single electron ionization from DM-electron scattering Essig et al. (2012a, b); Graham et al. (2012); An et al. (2015); Lee et al. (2015); Essig et al. (2016); Hochberg et al. (2016a, b); Cavoto et al. (2016); Derenzo et al. (2017); Essig et al. (2017a); Hochberg et al. (2017); Kouvaris and Pradler (2017); Essig et al. (2017b); Cavoto et al. (2018); Fichet (2018); Tiffenberg et al. (2017); Hochberg et al. (2018); Romani et al. (2018); Agnes et al. (2018); Crisler et al. (2018); Agnese et al. (2018); Essig et al. (2018); Abramoff et al. (2019); Emken et al. (2019).

Since momentum transfer from the DM to a target particle $T$ is most efficient when $m_{\rm DM}>m_{T}$, a bound atomic electron can capture an order-1 fraction of the DM kinetic energy for $m_{\rm DM}>m_{e}$ and be ionized. Similar reasoning would suggest that DM lighter than an atomic nucleus cannot efficiently transfer momentum to the nucleus, which is why experiments searching for nuclear recoils are typically insensitive to $m_{\rm DM}<m_{p}$. However, in a bound atomic system, momentum transfer to the entire atom will be redistributed among electrons and the nucleus through the electronic binding energy. This is known as the Migdal effect Migdal (1941); Vergados and Ejiri (2005); Moustakidis et al. (2005); Bernabei et al. (2007); Ibe et al. (2018); Dolan et al. (2018); Bell et al. (2019) and can also result in a final state with an ionized electron and a recoiling atom. Until now, the Migdal effect has solely been used to set limits on nucleon coupling for low mass weakly interacting massive particle models Akerib et al. (2019); Armengaud et al. (2019); Liu et al. (2019); Aprile et al. (2019a).

The main result of this paper is the following: in models where sub-GeV DM couples comparably to electrons and protons, the ratio of the differential ionization rate $dR_{M}/d\mathbf{q}$ due to the Migdal effect (which we will refer to as “Migdal scattering” for brevity) to the corresponding direct electron scattering rate $dR_{e}/d\mathbf{q}$ satisfies

[TABLE]

where $m_{N}$ is the mass of the target nucleus, $Z$ its atomic number, $\mathbf{q}$ the 3-momentum transferred from the DM to the atom, and $r_{a}$ is an effective atomic radius which we will define more precisely in Sec. III.2.111If DM couples equally to protons and neutrons, $Z$ should be replaced by the mass number $A$. As we will show, the rate computation for both electron scattering and Migdal scattering involves identical atomic ingredients because the initial hard scatter factorizes from the dynamics of ionization in bound atoms. However, the scattering probability in the latter case is enhanced by $Z^{2}$ due to coherent scattering off the nucleus, and simultaneously suppressed by the small electron mass compared to the heavy nucleus, though this suppression is mitigated somewhat when the momentum transfer to the atom is large.

Due to these competing effects, the Migdal scattering rate in heavy atoms such as Xe is generically dominated by the largest kinematically-permitted momentum transfers, typically hundreds of keV, which are small on nuclear scales but large on electron scales; by contrast, direct DM-electron scattering is dominated by the smallest momentum transfers. Thus, the direct DM-electron scattering rate generically dominates over Migdal scattering for light-mediator exchange, which favors small momentum transfers, whereas Migdal scattering can dominate for heavier mediators and heavier DM which imparts larger momentum transfers to the target system. When Eq. (1) is integrated over the momentum transfer $\mathbf{q}$, the total rate $R_{M}$ will then dominate over $R_{e}$ for sufficiently heavy DM.

This paper is organized as follows. In Sec. II we define our reference model of DM coupling to both electrons and protons. In Sec. III we develop in parallel the formalisms for electron and Migdal scattering, illustrating their similarities and differences. In Sec. IV we discuss the conversion from electron recoil spectra to observed ionization spectra in Xe and highlight the importance of the electron binding energies and wavefunctions in obtaining accurate exclusion curves. We conclude in Sec. V with a comparison of Migdal and electron exclusion curves for XENON10 and XENON100 data Angle et al. (2011a); Aprile et al. (2014), as well as the recent XENON1T results on light DM scattering Aprile et al. (2019b). We emphasize throughout that considerable theoretical uncertainty exists as to the correct choice of wavefunctions to use in computing limits on both Migdal and electron scattering. Consequently, our limits presented here should be considered provisional pending a dedicated analysis of relativistic and electron correlation effects in heavy atomic systems.

II Reference Model

Our benchmark model consists of a DM candidate $\chi$, which scatters off both electrons and protons through the exchange of a massive dark photon ${A^{\prime}}$ Holdom (1986); Okun (1982). The ${A^{\prime}}$ kinetically mixes with the visible photon, and after rotating away the kinetic mixing term $\frac{\epsilon}{2}F^{\mu\nu}F^{\prime}_{\mu\nu}$, the Lagrangian for this scenario contains

[TABLE]

where $\epsilon$ is the kinetic mixing parameter, $m_{{A^{\prime}}}$ is the ${A^{\prime}}$ mass, $J_{\rm EM}^{\mu}=\sum_{f}Q_{f}\bar{f}\gamma^{\mu}f$ is the electromagnetic current of all Standard Model fermions $f$ with charges $Q_{f}$, $g_{D}~{}\equiv~{}\sqrt{4\pi\alpha_{D}}$ is the dark photon gauge coupling, and $J_{D}^{\mu}~{}=\bar{\chi}\gamma^{\mu}\chi$ or $i\chi^{*}\partial_{\mu}\chi$ is the DM current for a Dirac fermion or complex scalar DM candidate, respectively. The fiducial non-relativistic cross section for $\chi$ scattering off a free-particle target $T$ is defined at a reference 3-momentum transfer $\mathbf{q}_{0}$ as

[TABLE]

where $\mu_{\chi T}$ is the $\chi$-$T$ reduced mass; by coincidence this same parametric expression hold for both complex scalar and Dirac fermion DM candidates.

This popular scenario features comparable mediator couplings to electrons and protons, so it serves as a good benchmark for comparing DM-induced ionization from direct electron scattering and Migdal scattering. Indeed, both processes will always be present, so for the remainder of this paper, we will only consider this model. However, comparable quark and lepton couplings are by no means unique to dark photons. This feature applies to most anomaly-free $U(1)$ extensions to the SM (e.g. gauged $B-L$) whose gauge bosons couple to DM Bauer et al. (2018); it is also generic for (pseudo)scalar-mediated DM scattering to feature comparable electron and proton couplings Krnjaic (2016); Knapen et al. (2017).

We note for completeness that, in the general case with arbitrary mediator couplings to different particle species, the relationship between Migdal and electron scattering in Eq. (1) becomes

[TABLE]

where $g_{e}$ and $g_{p,n}$ are the mediator’s couplings to electrons, protons, and neutrons respectively. Thus, unless $g_{p,n}\ll g_{e}$, there is a kinematic regime for which Migdal scattering can predominate over electron scattering.

III Comparison of Electron Scattering and Migdal Scattering

In this section, we carefully define the kinematics and dynamics relevant for sub-GeV DM interacting with atoms through electron and Migdal scattering. We will work in the framework of the dark photon model described above in Sec. II, where the dark photon mediates DM-SM interactions and couples equally to electrons and protons with strength $\epsilon|e|$, but our results are applicable to any model where the momentum dependence of DM-electron and DM-proton scattering (that is, the form factor $F_{\rm DM}(q)$ defined below Eq. (16)) is identical.

III.1 Kinematics

In both Migdal and electron scattering, the incoming and outgoing states are the same: a DM particle plus a bound atom and a DM particle plus an ionized atom plus an unbound electron, respectively. The incoming DM is assumed to be a plane wave, which is both an energy eigenstate and a momentum eigenstate. The incoming atom (at rest in the lab frame) is an energy eigenstate, and is also a momentum eigenstate for the total momentum of the atom $\mathbf{p}_{A}=\mathbf{p}_{N}+\sum_{i=1}^{Z}\mathbf{p}_{i}$, where the sum runs over the $Z$ electrons in the electron cloud of the (neutral) atom. The outgoing DM is also a plane wave, but the outgoing atom can either be treated as an ionized atom with a separate ionized electron, or an atom in an excited state where the ionized electron belongs to the continuum spectrum of the atomic Hamiltonian. Following Ref. Ibe et al. (2018), we will take the second perspective where energy-momentum conservation is more transparent, in which case the entire atom recoils with velocity $\mathbf{v}_{A}$ and has momentum $\mathbf{p}_{A}=\overline{m}_{A}\mathbf{v}_{A}$, where $\overline{m}_{A}=m_{N}+Zm_{e}$ is the nominal mass of the atom neglecting binding energy. In all cases we will consider, it is appropriate to approximate $\overline{m}_{A}$ by $m_{N}$ since the nucleus is so much heavier than the electron cloud. The energetics of the ionized electron are accounted for by treating it as an excited state of the electron cloud.

To summarize, when treating the atom as a composite system of electrons and nucleus with a spectrum of internal energy levels, both energy and momentum are conserved in DM-atom interactions. For DM with mass $m_{\chi}$, incoming velocity $\mathbf{v}$, and outgoing momentum $\mathbf{p}^{\prime}_{\chi}$, momentum conservation requires

[TABLE]

and energy conservation requires

[TABLE]

where $\mu_{\chi N}=m_{\chi}m_{N}/(m_{\chi}+m_{N})$ is the DM-nucleus reduced mass and $\Delta E_{e}\equiv E_{e,f}-E_{e,i}$ is the energy transferred to the scattered electron.

We emphasize that these kinematics are identical for electron scattering and Migdal scattering, provided the ionized electron is treated as a scattering state of the electron cloud Hamiltonian. In thinking of the nucleus and electrons as a single many-particle system in this formalism, it helps to regard $\mathbf{q}$ as simply the momentum transferred from the DM, rather than as a momentum transferred to any particular constituent of the target system. However, since the nucleus makes up the vast majority of the mass of the atom, one may think of $\mathbf{q}$ as the nuclear recoil momentum, as shown in Eq. (5).

III.2 Dynamics

While the kinematics of Migdal and electron scattering are identical, their dynamics differ in a crucial way depending on whether DM interacts directly with electrons or nuclei. In the language of nonrelativistic quantum mechanics, the perturbing Hamiltonian for the DM-atom interaction in the case of electron scattering is

[TABLE]

where $\mathcal{M}_{e\chi}(q)$ is the Lorentz-invariant matrix element for DM scattering off a free electron through 4-momentum transfer $q\approx(0,\mathbf{q})$, and $\mathbf{x}_{\chi}$ and $\mathbf{x}_{i}$ are the position operators for the DM and electrons, respectively. Because the DM interacts directly with electrons, we can ignore the nuclear part of the atomic Hamiltonian, and the rate will be proportional to Essig et al. (2012a)

[TABLE]

where we have made the approximation that the initial- and final-state electron cloud wavefunctions ($\Psi_{i}$ and $\Psi_{f}$, respectively) factorize such that only a single electron (with position operator $\mathbf{x}$) participates in a transition between the single-electron states $\psi_{i}$ and $\psi_{f}$. The scale of atomic wavefunctions is parametrically the size of the atom: to make this precise, we define an effective atomic radius $r_{a}$ as $r_{a}\equiv 1/|\mathbf{q}|_{\rm max}$, where $|\mathbf{q}|_{\rm max}$ is the momentum transfer at which the matrix element in Eq. (8) is maximized, for a given choice of initial and final states. Thus we may expect the electronic matrix element in Eq. (8) to be as large as $\mathcal{O}(1)$ when $|\mathbf{q}|r_{a}=\mathcal{O}(1)$. Typical $r_{a}$ are on the scale of the Bohr radius $a_{0}$, which is $\sim 1/{\rm keV}$ in natural units, and hence $\mathbf{q}\cdot\mathbf{x}\gtrsim\mathcal{O}(1)$ where the atomic wavefunctions have support. Note that for DM heavier than 1 MeV, a momentum transfer of $1/r_{a}$ is always kinematically allowed, since the DM carries momentum of at least 1 keV.

For Migdal scattering, where the fundamental DM-atom interaction is with the nucleus, the interaction Hamiltonian is

[TABLE]

where $\mathbf{x}_{N}$ is the position operator for the nucleus. $H_{{\rm int},N}$ does not involve the electron position operators $\mathbf{x}_{i}$, so by itself, it cannot induce electronic transitions. However, the light-crossing time of the electron cloud is $\sim{\rm nm}/c\sim 5\ {\rm keV}^{-1}$, so the timescale of momentum transfer to the nucleus is “fast” as long as the mediator mass $m_{A^{\prime}}$ satisfies $m_{A^{\prime}}\gg{\rm keV}$.222By construction, this always holds for contact interactions, $F_{\rm DM}=1$, but it is not clear to us if the formalism developed in Ref. Ibe et al. (2018) remains valid for ultralight mediators with $F_{\rm DM}\propto q^{-2}$, where for sufficiently small momentum transfers, the timescale for momentum transfer can be slow enough that the atomic state changes adiabatically. In this regime, the entire atom suddenly acquires velocity $\mathbf{v}_{A}$ but leaves behind its stationary electrostatic potential; the electrons of the moving atom are no longer in energy eigenstates of the old electron cloud Hamiltonian. As a result, electronic transitions can arise, but not through a matrix element with the perturbing Hamiltonian $H_{{\rm int},N}$.

Rather, following Ref. Ibe et al. (2018), we construct approximate energy eigenstates of the moving atom by applying a Galilean transformation with velocity parameter $\mathbf{v}_{A}$, which results in a final-state atomic wavefunction containing a phase $\exp(i\sum_{i=1}^{Z}\mathbf{q}_{e}\cdot\mathbf{x}_{i})$ multiplying the full wavefunction of the atom at rest, where $\mathbf{q}_{e}\equiv m_{e}\mathbf{v}_{A}$.333It is important to note that $\mathbf{q}_{e}$ is not to be interpreted as the momentum of the outgoing electron, but rather the effective momentum which appears in the matrix element due to the Galilean transformation to account for the atomic recoil velocity $\mathbf{v}_{A}$. The advantage of defining this quantity is to explicitly enable the dipole approximation of the matrix element for the masses considered here. Care should be taken not to place too much physical interpretation in this quantity, and so we translate our final result in Eq. (12) below back in terms of the physical momentum of the atomic system q. Consequently, as shown in Ref. Ibe et al. (2018), the matrix element of the nuclear wavefunction with $H_{{\rm int},N}$ results in a factor of $\mathcal{M}_{eN}(q)$ times the overlap of the electronic wavefunctions, and thus

[TABLE]

where $\Psi_{{\mathbf{v}_{A}}}$ is the Galilean transformation of the initial state $\Psi_{i}$ of the electron cloud, with velocity parameter $\mathbf{v}_{A}$. We have made the same approximation as in Eq. (8) that only a single electronic transition contributes; note that $\mathbf{q}_{e}$ instead of $\mathbf{q}$ appears in the exponent.

By momentum conservation from Eq. (5),

[TABLE]

Unlike the case for electron scattering, where $\mathbf{q}\cdot\mathbf{x}\gtrsim 1$ for all DM masses greater than an MeV, $\mathbf{q}_{e}\cdot\mathbf{x}\ll 1$ for all sub-GeV DM because $m_{e}/m_{N}\ll 1$. Indeed, $r_{a}\sim 4a_{0}$ for xenon, and $\mathbf{q}_{e}<1/r_{a}$ as long as $m_{\chi}<100\ {\rm GeV}$. As $\mathbf{q}_{e}\to 0$, the matrix element in Eq. (10) must vanish because $\psi_{f}$ and $\psi_{i}$ are energy eigenstates of the same Hamiltonian with different energy eigenvalues, by assumption. Hence the leading order term in the Taylor expansion of the exponential is linear in $\mathbf{q}_{e}$, and $R_{M}$ scales as $\mathbf{q}_{e}^{2}=\mathbf{q}^{2}(m_{e}/m_{N})^{2}$. There are also additional selection rules now that the matrix element has a dipole form, $\langle\psi_{f}|\mathbf{x}|\psi_{i}\rangle$, but in general, for a given $\mathbf{q}$ and a choice of initial- and final-state wavefunctions, the ratio of Migdal and electron scattering spectra for each $i\to f$ transition scales as

[TABLE]

where the appearance of $r_{a}$ arises from the expectation value of $\mathbf{x}$, which is parametrically of order $r_{a}$. We have written the above relation as an inequality because, for sufficiently large momentum transfers ($|\mathbf{q}|\gg{\rm keV}$), the exponential in Eq. (8) will oscillate rapidly and $R_{e}$ will become suppressed, thereby enhancing the Migdal rate relative to the electron scattering rate. On the other hand, the Migdal spectrum may always be approximated by a dipole matrix element which scales with $\mathbf{q}^{2}r_{a}^{2}$, since $\mathbf{q}_{e}<1/r_{a}$ for all targets relevant for sub-GeV DM.

III.3 Spectra and rates

To compute the ionization rate for both processes, we must integrate over the momentum transfer $\mathbf{q}$ and the DM velocity $\mathbf{v}$, and sum over the final electronic states $\psi_{f}$, weighted by a delta function enforcing energy conservation (3-momentum conservation has already been enforced in the definitions of $\mathbf{q}$ and $\mathbf{q}_{e}$ above).444Note that integrating over $\mathbf{q}$ is equivalent to integrating over the nuclear recoil energy $E_{R}\approx\mathbf{q}^{2}/(2m_{N})$, since in this paper we are concerned only with the electronic energy spectrum. We perform the integral over $\mathbf{v}$ by approximating the DM velocity distribution as spherically symmetric, $f(\mathbf{v})=f(v)$, such that the total rate between initial state $i$ and final state $f$ is

[TABLE]

where $\rho_{\chi}$ is the local DM density. For the sum over final states, we choose the normalization Essig et al. (2012a, 2017b)

[TABLE]

which is appropriate for scattering states in a spherically-symmetric potential which have asymptotic momentum $k^{\prime}=\sqrt{2m_{e}E_{e}}$ and angular momentum quantum numbers $l^{\prime}$ and $m^{\prime}$. Here, $E_{e}$ is the recoil energy of the ionized electron asymptotically far away from the ionized atom; from now on our final state $f$ will always be a scattered electron of energy $E_{e}$, and the initial state $i$ will be a bound state of (negative) energy $E_{nl}$ indexed by principal quantum number $n$ and angular momentum quantum number $l$, appropriate for a spherically-symmetric atom ignoring spin-orbit coupling and relativistic effects. The only difference between Migdal and electron scattering in the above procedure is the expression for $\sigma v_{i\to f}$.

To perform the integral over $\mathbf{q}$ and compute $\sigma v_{i\to f}$ we must specify the free-particle matrix elements. In the dark photon model, we can define a spin-averaged fiducial cross section for DM $\chi$ scattering off an isolated target $T$ with charge $|e|$ as in Eq. (3). For $T=p,e$, these fiducial cross sections satisfy

[TABLE]

so $\overline{\sigma}_{e}$ and $\overline{\sigma}_{p}$ are proportional by a factor which only depends on the DM mass $\chi$. To emphasize the point that $\overline{\sigma}_{e}$ and $\overline{\sigma}_{p}$ are related in this model, we shall refer to $\overline{\sigma}_{e}$ as simply $\overline{\sigma}$.

The fiducial cross section defined in Eq. (3) is related to the free-particle scattering matrix element as

[TABLE]

where we have assumed that the appropriate electron and DM spins have been summed and/or averaged. Here, $F_{\rm DM}(q)$ is the DM form factor which parametrizes all momentum dependence in the free-particle matrix element: if $m_{A^{\prime}}\ll m_{\chi}v$, $F_{\rm DM}(q)\propto 1/q^{2}$, while if $m_{A^{\prime}}\gg m_{\chi}v$, $F_{\rm DM}(q)=1$. Note that in the dark photon model with equal proton and electron couplings,

[TABLE]

in Eq. (12), where $F_{N}$ is the form factor of the nucleus; this relation between the matrix elements gives Eq. (1).

Putting all the pieces together, the electron recoil spectrum per unit detector mass for both electron and Migdal scattering is

[TABLE]

where $N_{T}$ is the number of atomic targets and

[TABLE]

Here, we have solved the delta function for energy conservation, $\delta(E_{e}-E_{nl}+\frac{\mathbf{q}^{2}}{2\mu_{\chi N}}-\mathbf{q}\cdot\mathbf{v})$, to perform the integral over the DM velocity distribution, resulting in a factor of $\eta(v_{\rm min})\equiv\langle v^{-1}\theta(v-v_{\rm min})\rangle$, the mean inverse DM speed in the lab frame, as a function of

[TABLE]

which is the minimum DM velocity required to ionize the target electron through a momentum transfer $|\mathbf{q}|$. The lab frame velocity distribution is cut off at $v_{E}+v_{\rm esc.}$, where $v_{E}\sim 240\ {\rm km/s}$ is the average speed of the Earth relative to the DM halo, and $v_{\rm esc.}=544\ {\rm km/s}$ is the galactic escape velocity (these parameters are chosen to facilitate comparisons with Ref. Essig et al. (2017b)).

The differences between the Migdal and electron scattering processes are entirely contained in the ionization form factors $|f_{e,M}(E_{e},\mathbf{q})|^{2}$, which are independent of all DM properties and depend only on the electronic and nuclear structure of the target. We will discuss in some detail in Sec. IV.2 the issues with accurately computing the atomic wavefunctions required for these ionization form factors. For electron scattering,

[TABLE]

where $\psi^{i}_{E_{nl}}$ is a bound orbital of energy $E_{nl}$ with unit norm (the factor of 2 accounts for the approximate spin degeneracy of occupied states) and $\psi^{f}_{E_{e}}$ is an unbound electronic state of energy $E_{e}$ (indexed by the continuously-valued energy $E_{e}$ and the angular momentum quantum numbers $l^{\prime}$ and $m^{\prime}$), normalized to

[TABLE]

where $k_{1,2}=\sqrt{2m_{e}E_{1,2}}$. For Migdal scattering, the analogous ionization form factor is

[TABLE]

The differences with respect to electron scattering are the appearance of $Z^{2}$ from coherent scattering off all the protons in the nucleus, a nuclear form factor $F_{N}$ parametrizing loss of coherence at large momentum transfers (which is largely irrelevant for the sub-MeV momentum transfers typical of sub-GeV DM), and the appearance of $\mathbf{q}_{e}$ instead of $\mathbf{q}$ in the matrix element between initial and final states.

The key quantity controlling the relative size of Migdal and electron scattering rates is $\mathbf{q}$. From Eq. (20), the smallest allowed $|\mathbf{q}|$ is

[TABLE]

where $E_{b}$ is the first ionization energy (positive by convention) of the atom in question, and $v_{\rm max}$ is the largest possible DM speed, which is the Galactic escape velocity in the lab frame. Note that $|\mathbf{q}|_{\rm min}$ is independent of the DM mass: for xenon with $E_{b}\sim 12\ {\rm eV}$, and $v_{\rm max}\sim 770\ {\rm km}/{\rm s}$, $|\mathbf{q}|_{\rm min}\sim 5\ {\rm keV}$. Thus for all kinematically-allowed momentum transfers, $|\mathbf{q}|r_{a}>1$, and electron scattering is dominated by the smallest possible $\mathbf{q}$ before $f_{e}(E_{e},\mathbf{q})$ is suppressed by the quickly-oscillating exponential in the matrix element. On the other hand, the largest allowed $|\mathbf{q}|$ is

[TABLE]

which grows with DM mass and can be as large as hundreds of keV for $m_{\chi}=\mathcal{O}(100\ {\rm MeV})$. For these momentum transfers, $|f_{M}(E_{e},\mathbf{q})|^{2}$ still does not feel any suppression from the nuclear form factor $F_{N}$, which is still $\sim 1$ for $|\mathbf{q}|\lesssim{\rm MeV}$, and likewise is still in the regime of small $\mathbf{q}_{e}$ and so grows with $\mathbf{q}^{2}$. Thus the Migdal ionization form factor is largest when $\mathbf{q}$ is the largest, and Migdal scattering is dominated by the largest kinematically-allowed momentum transfers.555In principle, there should be interference between electron and Migdal scattering, which have identical final states, but the distinct kinematics of these two processes should minimize these effects.

We illustrate this behavior in Fig. 1 for $m_{\chi}=300\ {\rm MeV}$, $E_{e}=10\ {\rm eV}$, and $F_{\rm DM}=1$. The left plot shows $|f_{e,M}(E_{e},\mathbf{q})|^{2}$ evaluated at $E_{e}=10\ {\rm eV}$, and the right plot shows the integrand of Eq. (19) which is weighted by $\eta(v_{\rm min})$. The $\mathbf{q}$ values plotted span the kinematically-allowed range between $|\mathbf{q}|_{\rm min}$ and $|\mathbf{q}|_{\rm max}$, as can be seen from the right plot where the velocity distribution cuts off the integrand at small and large $\mathbf{q}$. To compute $f_{e,M}$ for both electron and Migdal scattering from the same set of wavefunctions, the orthogonality of $\psi^{f}$ and $\psi^{i}$ is crucial, and to ensure this, both bound and free wavefunctions must be constructed from the same atomic Hamiltonian. A complete treatment would require a full numerical solution to the many-body Schrödinger equation for the atom in question, but here we capture the essential features by using hydrogenic wavefunctions for the 5p shell of Xe and matching the effective nuclear charge to the binding energy of the 5p state, with scattering states constructed from the same hydrogenic potential. This unphysical choice of wavefunctions is for illustrative purposes only; the wavefunctions used in the remainder of this paper are discussed in detail in Sec. IV.2. We note that for a DM form factor proportional to $1/q^{2}$, as would be the case for an ultralight dark photon mediator, the spectrum integrand is weighted by $1/q^{4}\sim 1/\mathbf{q}^{4}$ which heavily suppresses the Migdal spectrum compared to the electron spectrum for all DM masses.666As noted in Footnote 2, long-range interactions may result in adiabatic rather than sudden changes in atomic states during Migdal scattering, which would further suppress electronic transitions. For these form factors, electron scattering always dominates over Migdal scattering by several orders of magnitude, and as such, for the remainder of this paper we will focus on the case $F_{\rm DM}=1$.

We can confirm the relative strength of Migdal and electron scattering by considering the full electron recoil spectrum, as shown in Fig. 2. Here, to facilitate comparison with the literature, the electron spectra are calculated using the wavefunctions and binding energies of Ref. Essig et al. (2017b), and the Migdal spectra are calculated using $f_{M}$ as tabulated in Ref. Ibe et al. (2018). Despite some differences in these wavefunctions and binding energies (which we discuss further in Sec. IV.2), the intuition developed above holds well: for sufficiently large DM masses and equal couplings to protons and electrons, the Migdal spectrum dominates over the electron spectrum for all electron recoil energies.

IV Numerical modeling and systematic uncertainties

To apply the formalism of the previous section to experimental data, we must choose a model for generating ionization spectra from recoil spectra, as well as a set of atomic wavefunctions. We choose the previously published low-threshold analysis Essig et al. (2017b) of the XENON10 Angle et al. (2011b); Essig et al. (2012c) and XENON100 Aprile et al. (2016) detectors, as well as the newly-released S2-only analysis Aprile et al. (2019b) from the XENON1T detector Aprile et al. (2019c). These data are chosen for containing a relatively large exposure (for this mass range) at extremely low thresholds. All of these experiments use xenon time projection chambers to measure charge and light produced from energy deposited in liquid xenon.

An interaction in the xenon will create some number of xenon ions $N_{i}$ and initially excited xenon atoms $N_{\star}$. These atoms will form dimer states in the xenon which will release energy in the form of charge and UV scintillation photons. Scintillation photons produced at this step are immediately detected in what is referred to as the S1 signal. Before they can recombine, emitted electrons are drifted in a $\sim$300 V/cm electric field to a liquid-gas interface where a stronger ($\sim$10 kV/cm) extraction field is used to accelerate the electrons into the gas, producing a second (amplified) burst of light, referred to as the S2 signal. It has been well measured that (relatively speaking) interactions with xenon nuclei preferentially deposit energy via S1, whereas interactions with electrons in the detector preferentially deposit energy via S2 Dahl (2009); Akerib et al. (2016a, b); Aprile et al. (2018). In the case of sub-GeV DM interacting with a xenon atom through either electron or Migdal scattering, the momentum transfer to the recoiling atom is sufficiently small that the S1 signal is expected to be effectively zero. Thus, we consider an S2-only analysis using a 1(4) electron threshold for XENON10(100) Essig et al. (2017b) and a 5 electron threshold for XENON1T Aprile et al. (2019b).

IV.1 Ionization model and quantization

To compare with data, we must quantize the calculated recoil spectra in terms of the number of electrons extracted. For this, we adopt the ionization model from Refs. Essig et al. (2012a, 2017b) to determine the number of electrons produced ($n_{e}$) from an initial energy transfer $\Delta E_{e}$ which ionizes an electron from a specified electron shell with binding energy $E_{nl}$ to the continuum with energy $E_{e}$. We begin by considering the ejected electron, which has a probability $f_{R}$ to recombine (avoiding detection). According to the Thomas-Imel recombination model, $f_{R}$ is determined to be very small at low energies Dahl (2009); Sorensen and Dahl (2011) in good agreement with measurement Akerib et al. (2017) and is assumed to be zero for this analysis. We can thus write the probability of observing an initially produced electron as

[TABLE]

where the ratio of initially excited atoms to initially ionized atoms satisfies $N_{\star}/N_{i}\approx 0.2$ at high energies Doke et al. (2002); Aprile et al. (2007).

At high energies, the average energy $W$ required to produce one charge quantum in xenon is measured to be $W=13.7$ eV Dahl (2009). To convert $E_{e}$ and $E_{nl}$ into an expected quantized signal, we consider $n_{t}$ trials of a binomial process with probability of success, $f_{0}$, which satisfy

[TABLE]

where $|E_{nl}|-E_{b}$ is the available de-excitation energy. Thus, we can write

[TABLE]

where we assume that the number $n_{e}$ of quanta produced is equivalent to the number of electrons extracted (observed), as the extraction efficiency should be $\sim$100$\%$ Gushchin and Obodovskii (1982); Xu et al. (2019). Examples of quantized spectra for the $F_{\rm DM}=1$ case and different DM masses are shown in Fig. 3. This simple ionization model is sufficient for comparing electron and Migdal scattering here, but a more robust model would be needed to correctly interpret a signal through either channel.

Upper limits on the number of events for each value of $n_{e}$ in XENON10(100) have been determined in Ref. Essig et al. (2017b) to be $r_{1}<15.18$, $r_{2}<3.37$, $r_{3}<0.95$, and $r_{4}<0.35(0.17)$ counts kg*-1* day*-1*. Upper limits from XENON1T Aprile et al. (2019b) are not given directly. Instead, an upper bound of 22.5 events is reported in the range 165-275 photoelectrons in 15 ton-days of exposure777The average exposure for the range 165-198 photoelectrons has been determined from Figs. 1 and 4 of Ref. Aprile et al. (2019b).. We take the measured ratio of photoelectrons detected per $n_{e}$ (single electron gain) to be $\sim$33 Aprile et al. (2019b) and conservatively assume that all events in the reported range are for the lowest encompassed bin $n_{e}=5$, to obtain $r_{5}<0.0015$ counts kg*-1* day*-1*. These rates already account for analysis and detector efficiencies and thus can be directly compared against our quantized spectra.

IV.2 Electron Binding Energies and Wavefunctions

The dominant quantity controlling the sensitivity for small DM masses ($\lesssim 100\ {\rm MeV}$) is the outer-shell binding energy of xenon, or equivalently its first ionization energy. This quantity is well-measured experimentally with a value of 12.1 eV Brandi et al. (2001). Following Ref. Ibe et al. (2018), we define the electron shell binding energies $E_{nl}$ as an average over angular momentum states $\kappa$,

[TABLE]

but define the ionization energy $E_{b}$ as the minimum binding energy for all spin states in the atom. As a result, if spin-orbit coupling is ignored, $E_{b}=|E_{\rm 5p}|$, but the measured values show a $\sim 5\%$ difference between the two, as can be seen in Table 1.

Ideally, one wishes to compute the ionization spectrum using atomic wavefunctions with energy eigenvalues matching the observed binding energies. The binding energies used in the two analyses in Refs. Essig et al. (2017b); Ibe et al. (2018) are compared in Table 1. As can be seen from this table, obtaining accurate binding energies is not entirely trivial, as the FAC code Gu (2008) used in Ref. Ibe et al. (2018) gives an outer-shell binding energy of 9.8 eV, a significant difference of 20% from the observed value. On the other hand, the procedure used in Ref. Essig et al. (2017b) takes the calculated binding energy from the Roothaan-Hartree-Fock atomic wavefunctions tabulated in Ref. Bunge et al. (1993) and constructs outgoing wavefunctions from a hydrogenic potential with a shell-dependent effective nuclear charge which reproduces the appropriate binding energy for each electron shell. As the binding energies from Ref. Bunge et al. (1993) are closer to the observed values, this procedure retains the physical binding energies at the cost of losing orthogonality between initial and final electronic states, as well as any electron correlation effects. This orthogonality is crucial in order to obtain the behavior of $f_{M}$ as a function of $\mathbf{q}$, as discussed above, so it is not possible to compute Migdal scattering rates using these wavefunctions. However, as noted in Sec. III.3, electron scattering is dominated by the region where $\mathbf{q}\cdot\mathbf{x}\gtrsim 1$, so the form factor never probes the region where the matrix element must vanish as $\mathbf{q}\to 0$; thus, the overall kinematic features of electron scattering are sufficiently captured by this formalism.

While our interest in this paper is primarily a qualitative comparison of electron and Migdal scattering, we estimate one source of systematic error which can affect both processes by computing electron scattering rates using hydrogenic final-state wavefunctions (as in Ref. Essig et al. (2017b)), but constructed from the systematically-lower binding energies used in Ref.  Ibe et al. (2018). The main difference from the illustrative example computed in Sec. III.3 is that the initial-state wavefunctions are now taken from Ref. Bunge et al. (1993) instead of using a crude hydrogenic potential, in order to isolate the effects of binding energies and final-state wavefunctions. Fig. 4 compares the cross section limits obtained from these binding energies to the published electron scattering limits from Ref. Essig et al. (2017b). The systematic error on the electron scattering case associated with this procedure is less than an order of magnitude over the full DM mass range, but as expected, smaller binding energies lead to stronger cross section limits. The same procedure cannot be directly applied to Migdal scattering as the wavefunctions from Ref. Essig et al. (2017b) are not orthogonal, but the magnitude of the difference should be comparable (about a factor of 2 for masses above 100 MeV where Migdal scattering dominates).

We leave to future work a precise determination of experimental limits on Migdal scattering using more sophisticated quantum chemistry codes which correctly reproduce the observed binding energies with orthogonal wavefunctions. Indeed, recent progress for electron scattering has already been made by incorporating relativistic effects and electron-electron interactions into a many-body calculation, using the observed ionization energies from photoabsorption data as a figure of merit for the quality of the wavefunctions Pandey et al. (2018); the result is that at large $n_{e}$, the spectrum differs significantly from that obtained with hydrogenic final-state wavefunctions, potentially affecting the limits at large DM masses by an order of magnitude. In particular, the size of relativistic effects will grow as the atomic number of the atom increases, so this may be a significant source of systematic uncertainty for Migdal scattering in xenon.

V Results and conclusion

In this paper we have placed sub-GeV DM detection via electron and Migdal scattering on equivalent theoretical footing. Intriguingly, we have found that if DM couples comparably to electrons and protons through a contact interaction ($F_{\rm DM}=1)$, the Migdal rate can dominate for masses above $\sim$100 MeV. Thus, all existing limits for electron scattering in such models (such as dark photon-mediated scenarios), including those from XENON10, XENON100, and XENON1T Essig et al. (2017b); Aprile et al. (2019b), have omitted the dominant signal component at higher DM masses. In Fig. 5, we recalculate the full signal for DM-xenon scattering in XENON100 and XENON1T and extract upper bounds on $\overline{\sigma}$ which include both electron and Migdal scattering. It is clear that by exploiting the combination of electron and Migdal scattering, experiments with exposures and background rates comparable to XENON1T can start to probe the target parameter space for complex scalar DM freezing out through a heavy dark photon, $m_{A^{\prime}}\gg m_{\chi}$, but as we have emphasized, a definitive conclusion requires a more careful treatment of the atomic wavefunctions than has been used previously in the literature. In the event of a signal in the DM mass range where electron and Migdal rates are within a few orders of magnitude of each other, the unique spectral shapes can be used as a discriminant, though interference effects should be carefully considered.

Although our treatment here has focused on scattering from isolated atoms, we note that additional ionization from Migdal scattering should also contribute in semiconductor targets (mainly Si and Ge), for which low electronic band gaps represent the next frontier in electron-ionization direct detection. This additional signal channel can be probed by numerous future and ongoing experiments, including DAMIC at SNOLAB Aguilar-Arevalo et al. (2019), SENSEI Abramoff et al. (2019), SuperCDMS Agnese et al. (2018), and DAMIC-M Settimo (2018). However, a proper comparison of electron scattering and Migdal scattering in such materials is beyond the scope of the present work and deserves a dedicated study. At a minimum, the formalism for Migdal scattering must incorporate the nontrivial harmonic potential between neighboring ions, which may result in some portion of the DM energy loss appearing as phonons rather than electronic excitations.

Despite its robust theoretical underpinnings, Migdal scattering has not yet been experimentally observed. We emphasize that the same theoretical uncertainties which apply to Migdal scattering are present for the case of DM-electron scattering. We take advantage of this fact to show that the systematic uncertainty in the rate calculations due to different computations of the xenon binding energies can be estimated to be a factor of a few at high masses, where the Migdal scattering rate is dominant, but other sources of theoretical uncertainty due to relativistic effects and electron correlations may be equally important. The uncertainty due to the ionization model does not significantly affect the relative comparison of our calculated electron and Migdal scattering limits, but does explain the bulk of the difference between our limits and those in Ref. Aprile et al. (2019a) when interpreted in the dark photon model.

This analysis further highlights the importance of developing low-energy calibration techniques. We have shown that Migdal and electron scattering processes probe the same atomic wavefunctions, but in different kinematic regimes. As noted in Pandey et al. (2018), atomic many-body effects are crucial for understanding DM-atom interactions. Calibrations of both ionization processes, Migdal scattering and electron scattering, would help to resolve the theoretical uncertainty in the wavefunctions, which we believe has been substantially underappreciated by the sub-GeV DM community.

Acknowledgments. We thank Rouven Essig, Noah Kurinsky, Masahiro Ibe, Josef Pradler, Lucas Wagner, and Tien-Tien Yu for helpful conversations. We additionally thank Evan Shockley for helpful clarification regarding the recent XENON1T results. Fermilab is operated by Fermi Research Alliance, LLC, under Contract No. DE-AC02-07CH11359 with the US Department of Energy. This work was supported in part by the Kavli Institute for Cosmological Physics at the University of Chicago through an endowment from the Kavli Foundation and its founder Fred Kavli. We thank the KICP pheno journal club for providing a stimulating environment for discussion where this work originated. This research was supported by the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence “Origin and Structure of the Universe.” GK thanks MIAPP for hospitality while part of this work was completed.