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

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. In this work 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 $Z^2 (m_e/m_N)^2 \mathbf{q}^2$ where $m_N$ is the mass of the target nucleus and $\mathbf{q}$ is the 3-momentum transferred to the atom. 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.


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 [1]. In recent years, there have been several novel experimental techniques introduced to directly detect hitherto inaccessible DM candidates below the GeV scale [2]. One particularly promising strategy involves detecting single electron ionization from DM-electron scattering .
Since momentum transfer from the DM to a target particle T is most efficient when m DM > m T , a bound atomic electron can capture an order-one fraction of the DM kinetic energy for m 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 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 [28][29][30][31][32][33][34] 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 WIMP models [35][36][37][38].
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 /dq due to the Migdal effect (which we will refer to as "Migdal scattering" for brevity) to the corresponding direct electron scattering rate dR e /dq scales as where m N is the mass of the target nucleus, Z its atomic number, and q the 3-momentum transferred from the DM to the atom. 1 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 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 [39,40], as well as the recent XENON1T results on light DM scattering [41]. 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 χ, which scatters off both electrons and protons through the exchange of a massive dark photon A [42,43]. The A kinetically mixes with the visible photon, and after rotating away the kinetic mixing term 2 F µν F µν , the Lagrangian for this scenario contains where is the kinetic mixing parameter, m A is the A mass, J µ EM = f Q ff γ µ f is the electromagnetic current of all Standard Model fermions f with charges Q f , g D ≡ √ 4πα D is the dark photon gauge coupling, and J µ D = χγ µ χ or iχ * ∂ µ χ is the DM current for a Dirac fermion or complex scalar DM candidate, respectively. The fiducial non-relativistic cross section for χ scattering off a freeparticle target T is defined at a reference 3-momentum transfer q 0 as where µ χT is the χ-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 [44]; it is also generic for (pseudo)scalar-mediated DM scattering to feature comparable electron and proton couplings [45,46].

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 |e|, but our results are applicable to any model where the momentum dependence of DMelectron and DM-proton scattering (that is, the form factor F DM (q) defined below Eq. (15)) is identical.

A. 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 p A = p N + Z i=1 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. [32], we will take the second perspective where energy-momentum conservation is more transparent, in which case the entire atom recoils with velocity v A and has momentum p A = m A v A , where 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 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 χ , incoming velocity v, and outgoing momentum p χ , momentum conservation requires and energy conservation requires where µ χN = m χ m N /(m χ + m N ) is the DM-nucleus reduced mass and ∆E e ≡ 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 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 q as the nuclear recoil momentum, as shown in Eq. (4).

B. 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 where M eχ (q) is the Lorentz-invariant matrix element for DM scattering off a free electron through 4-momentum transfer q ≈ (0, q), and x χ and 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 [3] where we have made the approximation that the initialand final-state electron cloud wavefunctions (Ψ i and Ψ f , respectively) factorize such that only a single electron (with position operator x) participates in a transition between the single-electron states ψ i and ψ f . For Migdal scattering, where the fundamental DMatom interaction is with the nucleus, the interaction Hamiltonian is where x N is the position operator for the nucleus. H int,N does not involve the electron position operators x i , so by itself, it cannot induce electronic transitions. However, the light-crossing time of the electron cloud is ∼ nm/c ∼ 5 keV −1 , so the timescale of momentum transfer to the nucleus is "fast" as long as the mediator mass m A satisfies m A keV. 2 In this regime, the entire atom suddenly acquires velocity 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 int,N .
Rather, following Ref. [32], we construct approximate energy eigenstates of the moving atom by applying a Galilean transformation with velocity parameter v A , which results in a final-state atomic wavefunction containing a phase exp(i Z i=1 q e · x i ) multiplying the full wavefunction of the atom at rest, where q e ≡ m e v A . Consequently, as shown in Ref. [32], the matrix element of the nuclear wavefunction with H int,N results in a factor of M eN (q) times the overlap of the electronic wavefunctions, and thus where Ψ v A is the Galilean transformation of the initial state Ψ i of the electron cloud, with velocity parameter v A . We have made the same approximation as in Eq. (7) that only a single electronic transition contributes; note that q e instead of q appears in the exponent.
The key observation of this paper is that by momentum conservation from Eq. (4), For sub-GeV DM, |q| is typically of order m χ v, which is at least ∼ keV for m χ > MeV and hence q · x O(1) where the atomic wavefunctions have support (within a few Bohr radii a 0 of the nucleus). Thus we may expect the electronic matrix element in Eq. (7) to be as large as O(1) when this momentum transfer is kinematically allowed. On the other hand, q e · x 1 because m e /m N 1. As q e → 0, the matrix element in Eq. (9) must vanish because ψ f and ψ 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 q e , and R M scales as q 2 e = q 2 (m e /m N ) 2 . There are also additional selection rules now that the matrix element has a dipole form, ψ f |x|ψ i , but in general, for a given q and a choice of initial-and final-state wavefunctions, the ratio of Migdal and electron scattering spectra for each i → f transition scales as where the Bohr radius, a 0 , arises from the expectation value of position when expanding Eq. (9). We have written the above relation as an inequality because, for sufficiently large momentum transfers (|q| keV), the exponential in Eq. (7) will oscillate rapidly and R e will become suppressed, thereby enhancing the Migdal rate relative to the electron scattering rate.

C. Spectra and rates
To compute the ionization rate for both processes, we must integrate over the momentum transfer q and the DM velocity v, and sum over the final electronic states ψ f , weighted by a delta function enforcing energy conservation (3-momentum conservation has already been enforced in the definitions of q and q e above). 3 We perform the integral over v by approximating the DM velocity distribution as spherically symmetric, f (v) = f (v), such that the total rate between initial state i and final state f is where ρ χ is the local DM density. For the sum over final states, we choose the normalization [3,16] which is appropriate for scattering states in a sphericallysymmetric potential which have asymptotic momentum k = √ 2m e E e and angular momentum quantum numbers l and m . 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 σv i→f .
To perform the integral over q and compute σv i→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 χ scattering off an isolated target T with charge |e| as in Eq. (3). For T = p, e, these fiducial cross sections satisfy 3 Note that integrating over q is equivalent to integrating over the nuclear recoil energy E R ≈ q 2 /(2m N ), since in this paper we are concerned only with the electronic energy spectrum.
so σ e and σ p are proportional by a factor which only depends on the DM mass χ. To emphasize the point that σ e and σ p are related in this model, we shall refer to σ e as simply σ.
The fiducial cross section defined in Eq. (3) is related to the free-particle scattering matrix element as where we have assumed that the appropriate electron and DM spins have been summed and/or averaged. Here, F DM (q) is the DM form factor which parametrizes all momentum dependence in the free-particle matrix element: Note that in the dark photon model with equal proton and electron couplings, in Eq. (11), 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 where N T is the number of atomic targets and Here, we have solved the delta function for energy conservation, δ(E e −E nl + q 2 2µ χN −q · v), to perform the integral over the DM velocity distribution, resulting in a factor of η(v min ) ≡ v −1 θ(v − v min ) , the mean inverse DM speed in the lab frame, as a function of which is the minimum DM velocity required to ionize the target electron through a momentum transfer |q|. The lab frame velocity distribution is cut off at where v E ∼ 240 km/s is the average speed of the Earth relative to the DM halo, and v esc. = 544 km/s is the galactic escape velocity (these parameters are chosen to facilitate comparisons with Ref. [16]). The differences between the Migdal and electron scattering processes are entirely contained in the ionization form factors |f e,M (E e , 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 B the issues with accurately computing the atomic wavefunctions required for these ionization form factors. For electron scattering, where ψ i E nl is a bound orbital of energy E nl with unit norm and ψ f Ee is an unbound electronic state of energy E e (indexed by the continuously-valued energy E e and the angular momentum quantum numbers l and m ), normalized to ψ E1,l,m |ψ E2,l ,m = 1 where k 1,2 = 2m e E 1,2 . For Migdal scattering, the analogous ionization form factor is 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 q e instead of q in the matrix element between initial and final states.
The key quantity controlling the relative size of Migdal and electron scattering rates is q. From Eq. (19), the smallest allowed |q| is where E b is the first ionization energy (positive by convention) of the atom in question, and v max is the largest possible DM speed, which is the Galactic escape velocity in the lab frame. Note that |q| min is independent of the DM mass: for xenon with E b ∼ 12 eV, and v max ∼ 770 km/s, |q| min ∼ 5 keV. Thus for all kinematically-allowed momentum transfers, |q|a 0 > 1, and electron scattering is dominated by the smallest possible q before f e (E e , q) is suppressed by the quicklyoscillating exponential in the matrix element. On the other hand, the largest allowed |q| is which grows with DM mass and can be as large as hundreds of keV for m χ = O(100 MeV). For these momentum transfers, |f M (E e , q)| 2 still does not feel any suppression from the nuclear form factor F N , which is still ∼ 1 for |q| MeV, and likewise is still in the regime of small q e and so grows with q 2 . Thus the Migdal ionization form factor is largest when q is the largest, and Migdal scattering is dominated by the largest kinematically-allowed momentum transfers. 4 We illustrate this behavior in Fig. 1 for m χ = 300 MeV, E e = 10 eV, and F DM = 1. The left plot shows |f e,M (E e , q)| 2 evaluated at E e = 10 eV, and the right plot shows the integrand of Eq. (18) which is weighted by η(v min ). The q values plotted span the kinematicallyallowed range between |q| min and |q| max , as can be seen MeV (blue) and 300 MeV (red). For both spectra, we show the inclusive rates summed over all E nl → Ee transitions in xenon, where contributions from nl = 5p, 5s, and 4d dominate. Migdal spectra are computed using the wavefunctions and binding energies from Ref. [32], while electron spectra are computed using wavefunctions and binding energies from Ref. [16]. The differences between these choices are irrelevant for our qualitative argument here and are discussed further in Sec. IV B. from the right plot where the velocity distribution cuts off the integrand at small and large q. To compute f e,M for both electron and Migdal scattering from the same set of wavefunctions, the orthogonality of ψ f and ψ 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 B. 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 ∼ 1/q 4 which heavily suppresses the Migdal spectrum compared to the electron spectrum for all DM masses. 5 For these form factors, electron scattering always dominates over Migdal 5 As 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.
scattering by several orders of magnitude, and as such, for the remainder of this paper we will focus on the case F 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. [16], and the Migdal spectra are calculated using f M as tabulated in Ref. [32]. Despite some differences in these wavefunctions and binding energies (which we discuss further in Sec. IV B), 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 [16] of the XENON10 [47,48] and XENON100 [49] detectors, as well as the newlyreleased S2-only analysis [41] from the XENON1T detector [50]. 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 (TPCs) 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 . 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 ∼300 V/cm electric field to a liquidgas interface where a stronger (∼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 [51][52][53][54]. 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) [16] and a 5 electron threshold for XENON1T [41].

A. 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. [3,16] to determine the number of electrons produced (n e ) from an initial energy transfer ∆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 [51,55] in good agreement with measurement [56] and is assumed to be zero for this analysis. We can thus write the probability of observing an initially produced electron as where the ratio of initially excited atoms to initially ionized atoms satisfies N /N i ≈ 0.2 at high energies [57,58]. At high energies, the average energy W required to produce one charge quantum in xenon is measured to be W = 13.7 eV [51]. 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 where |E nl | − E b is the available de-excitation energy. Thus, we can write 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 ∼100% [59,60]. Example quantized spectra for the F 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. [16] 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 [41] are not given directly. Instead, an upper bound of 22.5 events is reported in the range 165-275 photoelectrons in 15 tonne-days of exposure 6 . We take the measured ratio of photoelectrons detected per n e (single electron gain) to be ∼33 [41] 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.

B. Electron Binding Energies and Wavefunctions
The dominant quantity controlling the sensitivity for small DM masses ( 100 MeV) is the outer-shell binding energy of xenon, or equivalently its first ionization  I. Comparison of the ionization energy E b and electron binding energies |E nl | (eV) of the 5p, 5s, and 4d electron shells of xenon from calculations using the formalisms of Ref. [32] and Ref. [16]. The measured values [62,63] tend to fall somewhere in the middle.
FIG. 4. Effect of wavefunction and binding energy choices on 90% CL limits on σ for direct DM-electron scattering in XENON10 with ne = 1 (blue), 2 (green), and 3 (purple) and in XENON100 with ne = 4 (red). Solid curves are published limits from Ref. [16] and dashed curves use binding energies from Ref. [32] to construct hydrogenic final-state wavefunctions.
energy. This quantity is well-measured experimentally with a value of 12.1 eV [62]. Following Ref. [32], we define the electron shell binding energies E nl as an average over angular momentum states κ, 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 5p |, but the measured values show a ∼ 5% difference between the two, as can be seen in Table I. 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 Refs. [16] and [32] are compared in Table I. As can be seen from this table, obtaining accurate binding energies is not entirely trivial, as the FAC code [64] used in Ref. [32] 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. [16] takes the calculated binding energy from the Roothaan-Hartree-Fock atomic wavefunctions tabulated in Ref. [61] 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. [61] 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 q, as discussed above, so it is not possible to compute Migdal scattering rates using these wavefunctions. However, as noted in Sec. III C, electron scattering is dominated by the region where q · x 1, so the form factor never probes the region where the matrix element must vanish as q → 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. [16]), but constructed from the systematically-lower binding energies used in Ref. [32]. The main difference from the illustrative example computed in Sec. III C is that the initial-state wavefunctions are now taken from Ref. [61] 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. [16]. 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. [16] 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 phtoabsorption data as a figure of merit for the quality of the wavefunctions [65]; 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  14) and equivalent to σe for electron scattering) for heavy A (FDM = 1) mediated DM-electron (dashed purple) and DM-Migdal (dashed red) scattering are shown for ne = 4 for XENON100 (left) and for ne = 5 for XENON1T (right). In the mass and coupling range shown on the plot, the XENON10 limits are sub-dominant. For comparison, we show the published electron scattering limits [16,41] computed with hydrogenic final-state wavefunctions and binding energies from Ref. [16] (solid purple); our electron scattering results use the smaller (unphysical) binding energies from [32] to facilitate a comparison with Migdal scattering using the same binding energies (see Sec. IV B). The thick blue curve is the complex scalar DM freeze-out target (particle-antiparticle symmetric DM population). Points along this curve account for the full DM abundance as long as m A mχ; near resonance at m A ≈ 2mχ this target moves down in the parameter space, but is otherwise robust [66,67]. The thin blue curve is the looser asymmetric Dirac fermion DM target. Any points above this line can account for the full DM abundance, but with different particle-antiparticle asymmetries [2,66]; points below this curve are excluded by Planck limits on CMB energy injection from the annihilation of the symmetric component [68]. The dotted blue curve taken from Ref. [2] represents sensitivity targets for ELDER DM [69]; points above this curve correspond to SIMP DM models with the same A mediator considered here [70]. Shaded regions represent an envelope of exclusions from beam dump searches (LSND [71], E137 [72,73], and MiniBooNE [74,75]), nuclear recoil direct detection limits from CRESST II [76], and the BaBar monophoton search for invisibly decaying dark photons [77][78][79]. 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 DM = 1), the Migdal rate can dominate for masses above ∼100 MeV. Thus, all existing limits for electron scattering in such models (such as dark photonmediated scenarios), including those from XENON10, XENON100, and XENON1T [16,41], 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 σ 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 com-plex scalar DM freezing out through a heavy dark photon, m A m χ , 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 electronionization direct detection. This additional signal channel can be probed by numerous future and ongoing experiments, including DAMIC at SNOLAB [80], SEN-SEI [26], SuperCDMS [24], and DAMIC-M [81]. 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 DMelectron 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. [38] 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 [65], 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.