The Migdal effect in semiconductors

When a nucleus in an atom undergoes a collision, there is a small probability to inelastically excite an electron as a result of the Migdal effect. In this Letter, we present a first complete derivation of the Migdal effect from dark matter-nucleus scattering in semiconductors, which also accounts for multiphonon production. The rate can be expressed in terms of the energy loss function of the material, which we calculate with density functional theory (DFT) methods. Because of the smaller gap for electron excitations, we find that the rate for the Migdal effect is much higher in semiconductors than in atomic targets. Accounting for the Migdal effect in semiconductors can therefore significantly improve the sensitivity of experiments such as DAMIC, SENSEI and SuperCDMS to sub-GeV dark matter.

Introduction -Direct detection of nuclear recoils from sub-GeV dark matter (DM) is challenging because the typical energy deposited in an elastic nuclear recoil scales as E N ∼ m 2 χ v 2 χ /m N , which is exceedingly difficult to detect for lighter dark matter candidates. Furthermore, low-energy nuclear recoils do not deposit much energy in the form of charge (electrons) or scintillation light, which are primary detection channels in many experiments. These considerations are drivers for new technologies and experiments capable of lower thresholds and phonon-based detection [1,2].
In the meantime, existing experiments can significantly extend their reach by searching for inelastic scattering processes, during which additional excitations are created, such as ionization electrons [3], photons [4], or plasmons [5,6]. While there is generally a rate penalty for such processes, there are two key advantages: the unfavorable energy relation in the previous paragraph can be broken and the additional excitation gives rise to charge signals that are more readily detected.
The Migdal effect [7,8] describes inelastic collisions where atoms are excited or ionized during the initial, hard nuclear recoil. For m χ 70 MeV, this process takes place on time scales much shorter than the time for the recoiling nucleus to travel a distance comparable to the interatomic spacing. It can therefore be factorized from secondary ionizations that can occur as the recoiling nucleus interacts with surrounding atoms, as described by Lindhard's theory [9]. As we will show, in crystals the underlying physics of the primary (Migdal) and secondary (Lindhard) ionizations is closely related, but they take place on different time scales.
So far, the Migdal effect and other inelastic processes have primarily been studied in atomic targets. As applied to dark matter direct detection, the most complete derivation of the Migdal effect can be found in Ref. [3], while additional discussion can be found in Refs. [10][11][12][13][14][15][16][17][18]. * simon.knapen@cern.ch † jkozaczuk@physics.ucsd.edu ‡ tongyan@physics.ucsd.edu The effect has been applied to set strong experimental limits in noble liquid targets [17][18][19][20]. Charge thresholds in low-threshold semiconductor experiments are even lower [21][22][23], which makes them even better suited to exploit the Migdal effect. However, inelastic processes are less well-studied in crystal targets, in part due to the more complicated spectrum of excitations. Moreover, the existing calculation for atomic targets [3] relies on boosting the system to the rest frame of the recoiling nucleus. While convenient for atomic calculations, this method breaks down for semiconductors, as the rest frame of the crystal is a preferred frame. So far, Refs. [15][16][17] presented estimates of the Migdal effect in solid state targets, but as of now a fully robust, first principles derivation is still lacking.
In this Letter, we show that the Migdal effect in semiconductors can be described by the 2 → 3 process of DM-nucleus scattering in association with a nucleuselectron Coulomb interaction. Furthermore, the details of the electronic band structure can be packaged into the well-studied energy loss function, Im(−1/ ), where is the dielectric function of the material. When the energy deposited into electronic excitations is close to the plasma frequency, plasmons are resonantly excited. (A plasmon resonance can be thought of as a collective electron excitation or as a longitudinal in-medium photon mode.) A first calculation of the rate for inelastic DM-nucleus scattering with associated plasmon production was presented in Ref. [6]. By directly calculating the non-resonant contributions, we show that plasmon production is simply the resonant component of the Migdal effect.
We will begin with a general description of the process in semiconductors and lay out the assumptions in our calculation. We present our main results of the Migdal rate in semiconductors, and then clarify its relation to the atomic Migdal effect. Finally, we present sensitivity estimates for experiments using Si and Ge targets. We leave the bulk of our calculations for the Appendices.
Description of process -For an elastic recoil off a free nucleus of mass m N , the typical momentum deposited by sub-GeV DM is q N ≈ m χ v χ ≈ MeV × (m χ /GeV) and the typical recoil energy is E N ≈ arXiv:2011.09496v2 [hep-ph] 3 Sep 2021 FIG. 1. We compute DM scattering off a nucleus in a harmonic crystal using the impulse approximation, which is valid when the time scale of the initial collision (t0) is short compared to the time scale to traverse its potential well (t1), set by the phonon frequency 1/ω. This holds for mχ 70 MeV. m 2 χ v 2 χ /m N ≈ 35 eV × (m χ /GeV) 2 , taking the example of a Si target. For sub-GeV DM, the energy and momentum scales are then comparable to various scales inherent to the crystal. Some care is therefore needed with respect to the regime of validity of our approximations.
Concretely, in a typical crystal, each nucleus sits in an approximately harmonic potential with size of ∼Å and frequencyω ≈ 30 − 50 meV (see Fig. 1). We will deal with m χ 10 MeV such that the inverse momentum transfer 1/q N Å . Then we can consider the interaction of the DM with a single nucleus, the so-called incoherent approximation 1 . We will thus compute DM scattering off a nucleus in the ground state of the potential, with associated nucleus-electron interaction.
To treat the excited states of the nucleus, we will rely on the impulse approximation, which is valid if the collision happens quickly relative to the time scale set by the potential well, 1/ω. (See e.g. [31]) The initial DMnucleus collision and the emission of the Migdal electrons take place on a time scale ∼ 1/E N , during which the nucleus remains near the minimum of the potential well (see Fig. 1). Only at a much later time 1/ω, the nucleus reaches the edge of the unit cell and loses its residual kinetic energy to phonons, or becomes unbound, depending on the initial collision energy. The impulse approximation thus allows us to model the excited states as plane waves for the duration of the hard collision 2 . We show this in Appendix B by calculating the all orders, multiphonon response of a harmonic crystal and find that the approximation holds if the momentum transfer satisfies q N √ 2m Nω . For DM with the standard velocity profile, this implies m χ 70 MeV. The incoherent approximation is then satisfied automatically as well. In practice, we will apply a set of cuts to exclude regions of the phase space where the approximations start to break down.
In this discussion we have focused on the DM momentum transfer to the nucleus, q N , because the momentum deposited in electrons will be much smaller. The Migdal rate then factorizes as: where ω is the energy deposited into electronic excitations, dσ el /dE N the elastic DM-nucleus cross section and dP (E N )/dω is the differential ionization probability. This is identical to the expansion made in the bremsstrahlung of a soft photon from a heavy charged particle, and we refer to it as the soft limit. The soft limit holds as long as |q N · k| m N ω and k q N . Estimating q N ∼ vm χ and k ∼ 1 − 10 keV, this translates to 50 MeV m χ 1 GeV and ω eV, which covers the most relevant parameter space. While our formal result does not rely on the soft limit, it is a useful technical and conceptual simplification when performing the phase space integrals, and is valid whenever the impulse approximation holds.
Finally, we will treat the nuclei and tightly-bound core electrons together as a particle with charge Z ion , and only consider excitations of the valence electrons. In other words, we assume an ion potential which behaves as Z ion e/r for large r compared to the wavefunctions of the inner shell electrons. To account for the effective ion charge at shorter distances, we include a momentumdependent Z ion (k) in the Fourier-transformed potential, which we obtain using tabulated ionic form factors [32]. While the electron-ion momentum exchange will be 10 keV such that the long-range behavior of the potential is most important, including Z ion (k) leads to O(1) rate increases.
Calculation -The Migdal rate is given by the leading order expansion in both the DM-nucleus and the electromagnetic interactions, analogous to bremsstrahlung. We assume a contact interaction between the DM and the nuclei, given by the Hamiltonian H χ = (2πb χ /m χ )δ(r χ − r N ) for m χ m N , with b χ the DM-nucleus scattering length and r χ , r N the position operators for the DM and nucleus. (The DM-nucleus elastic cross section is therefore given by σ N = 4πb 2 χ = A 2 σ n , with A and σ n respectively the atomic mass number and DM-nucleon elastic cross section.) The electron-nucleus interaction is H e = d 3 r −1 (r , r, ω) Z ion α/|r − r N | with r the position operator for the electron. α is the electromagnetic fine structure constant, and is the frequency-dependent, microscopic dielectric function, which encodes the screening by the spectator valence electrons.
In a crystal, the linear response depends on both r and r, up to the lattice periodicity. The Fourier transform of the response −1 (r , r, ω) can then be written as −1 KK (k, ω) ≡ −1 (k + K, k + K , ω) where k is in the first Brillouin Zone (BZ) and K, K are reciprocal lattice vec-tors. −1 KK can be regarded as a matrix in the reciprocal lattice vectors, but for Si and Ge we find the contribution of the off-diagonal pieces to be subleading. Here we just present results assuming a diagonal response matrix −1 KK and reserve the full result for Appendix A. Including the momentum-dependent ion charge, H e can then be written as We can apply Fermi's golden rule with second-order perturbation theory to compute the cross section for DM-nucleus inelastic scattering. We take the initial ions to be in a ground state of a harmonic crystal potential. Following the impulse approximation, we use plane waves for intermediate and final states. Meanwhile, the electron states are treated as Bloch states. Though the computation itself is a straightforward application of second order perturbation theory, the formulas and derivation are somewhat lengthy due to the appearance of the reciprocal lattice vectors and a form factor for the recoiling ion. We refer the reader to Appendix A for further details, and only present the final result here: where q N and p f are the final ion and DM momentum, respectively, and k + K is the momentum deposited to the electrons. In the first part of (3), we see the same factors and phase space integral that appear for elastic DM-nucleus scattering, except with the additional energy ω being deposited in electronic excitation. While for free nucleus scattering there would be a momentum conservation delta function, here we have a form factor F which encodes the details of the ion ground state, and for a harmonic crystal it is given by whereω is an oscillator frequency, averaged with respect to the density of states D(ω) and the thermal Bose factor. The remaining pieces of (3) contain the probability for exciting the electron. We sum over all initial and final electron states p e and p e + k, weighted by the occupation numbers f , and where band indices have been suppressed. The electronic wavefunction overlaps [p e + k|e ir·K |p e ] Ω are performed over the unit cell, and V is the volume of the crystal. In (3), the bracketed quantity can be rewritten in terms of the imaginary part of the dielectric function in the random phase approximation, Im [ KK (k, ω)]. Then we can write Im [ KK (k, ω)]/| KK (k, ω)| 2 = Im [−1/ KK (k, ω)], which is the energy loss function (ELF) governing energy loss of charged particles in a material. Physically, the ion-electron interaction in the inelastic process can be encapsulated in the same ELF as ions passing through a material. Since the ELF is a well-measured and calcu-lated quantity in many materials, this provides a useful starting point for numerical evaluations of (3).
In the soft limit |k + K| |q N |, the cross section factorizes as in (1), and the form factor F only modifies the elastic recoil cross section. Then the differential ionization probability is with v N ≡ q N /m N . This simplified formula is only valid for k in the first Brillouin zone, see Appendix A for the full expressions used in our numerical results. Eq. (6) was also derived in [6], but that work did not account for the ion ground state or electron momentum transfers outside of the first BZ, since it was focused on long-wavelength plasmons. Furthermore, [6] used an analytic approximation for (k, ω) near the plasmon pole. In the results below, we will study the impact of accounting for the ion ground state and use numerical calculations of (k, ω) valid away from the plasmon resonance. Before doing so, we clarify the relation of this process with the atomic Migdal effect.
Comparison with atomic Migdal effect -In Migdal's original calculation [7,8] for an atomic target, the ground state of the electron cloud (|i ) is first boosted to the rest frame of the moving nucleus |i → e imev N · β r β |i . He then computes the overlap with the where β runs over all the electrons in the atom. The transition probabilities |M if | 2 can then be evaluated with known atomic wave functions, and it was found that single ionizations dominate for sub-GeV dark matter [3].
To demonstrate the connection with the semiconductor Migdal effect derived above, we instead rewrite (7) using the following operator identity: is the total energy deposited and H 0 the electron Hamiltonian. We assume a non-relativistic 3 Hamiltonian such that the H 0 is a sum of kinetic terms, Coulomb interaction terms between electrons, and the Coulomb interaction of the electrons with the nucleus. Then the commutator β [p β , H 0 ] will be proportional to the total force from the nucleus, since the electronelectron forces cancel out. Contracting the matrix element with v N , we find (see also [33]) with r N the position operator of the nucleus. In the right-hand-side matrix element above, we see the time derivative of the dipole potential from a nucleus which has been displaced by |r N | |r β |. This is already very similar to the Coulomb potential in (2) and suggestive of the same physical interpretation as in the semiconductor case. One can now take the Fourier transform and evaluate the transition probability for single ionizations: where we have dropped the e −ik·r N phase factor in the soft limit. We have thus shown that the atomic Migdal effect has a form nearly identical to (5) for semiconductors, up to a few differences reflecting the different physical systems. In (5), the integral over k appears outside the amplitude squared; this reflects conservation of crystal momentum in the semiconductor, which requires that the final state have momentum p e + k, whereas the atomic states are not momentum eigenstates. The nucleus charge Z N appears here, since we are considering the all-electron wavefunctions, whereas in the semiconductor case we effectively integrated out the core electrons and treated the ion with effective charge Z ion (k). Finally, in (5), we accounted for the in-medium dielectric screening 1/| (k, ω)| 2 due to all the other electrons, which is neglected in the atomic case.
With this result, we find an equivalent formulation of the atomic Migdal effect, (9), but which has a physical interpretation that applies also in semiconductors. Previously, Ref. [34] assumed the atomic formulation in (7) could be generalized directly to semiconductors. However, the boosting argument used to obtain (7) does not apply in this case, since in a crystal there is a preferred coordinate frame. In other words, applying a boost operator would boost all nuclei in the lattice. Specifically, the two approaches are not equivalent because the operator relation (8) only applies for an individual atom: in the presence of a lattice of nuclei, we would have contributions from all nuclei on the right-hand-side of (8). Ref. [15] attempted to address this subtlety by using atom-centered localized Wannier functions in (7), to mitigate the contribution from the remaining nuclei in the crystal. For Si and Ge in particular, the Wannier approach is however found to be computationally challenging due to slow numerical convergence [15].
Between these different starting points, (9) has a clear physical interpretation as the in-medium analogue of bremsstrahlung, which nicely generalizes to the semiconductor case. This interpretation is discussed further in Appendix D. The final result in semiconductors can moreover be expressed in terms of the ELF, which can be calculated with a number of public codes. We therefore argue for this approach in generalizing the atomic Migdal effect.
Results -The Migdal rate is given by where N T is number of target nuclei per kilogram and n χ = ρ χ /m χ is the DM number density. We take ρ χ = 0.4 GeV/cm 3 and assume the standard halo model for the DM velocity distribution f (v) with escape velocity v esc = 500 km/s [41,42] and Earth velocity v e = 240 km/s. We calculate the dielectric function with the public code GPAW [43][44][45][46]. The wavefunctions are computed on an 8 × 8 × 8 k-space grid for Si, and a 12 × 12 × 12 grid for Ge. The TB09 exchangecorrelation functional [47] is used and local field effects are incorporated [48]. A scissor correction is applied to match on to the experimentally-measured band gap in both materials. We have averaged the ELF over crystal directions, and for computational reasons, we currently only include the diagonal components of the loss function, although we verified that the off-diagonal terms do not contribute more than an O(1) amount to the total rate. Further details, more refined numerical studies and other applications will be presented in a forthcoming publication [49].
The resulting rates and sensitivity estimates are shown in Fig. 2 for Si and Ge. The bands indicate an estimate of the theory uncertainty due the breakdown of the impulse approximation for low energy nuclear recoils. We find that our calculation starts to break down for m χ 50 MeV, at which point one has to go beyond the impulse approximation by matching onto the phonon regime. The dashed lines in the right-hand panel indicate the free ion approximation, where the F form factor is replaced with The differential rate is translated into total number of electron-hole pairs created (Q) using Eq. 5.1 of Ref. [35]. (Right) Expected 90% CL sensitivity to DM-nucleon cross section σn assuming a heavy mediator and 1 kg-year of exposure. We take Q ≥ 2 and zero background, corresponding to an upper limit of 2.4 events. The red line is a 90% CL limit obtained using the recent upper limit on the 2-electron rate from SENSEI [22], while the shaded region includes bounds from XENON1T [20], LUX [19], CRESST III [36] and CDEX [37], as well as a recast of XENON10 [38], XENON100 [39], and XENON1T [40] data in terms of the Migdal effect [17]. For comparison, we show a projection for the Migdal effect in Xenon (dotted line) based on the atomic ionization signal [17]. In both panels, the shaded bands are an estimate of the theoretical uncertainty due to the impulse approximation, obtained by varying the threshold on EN from 4ω to 9ω. See Appendix C for more details.
a δ-function. We find that the free ion approximation is excellent in the regime where the impulse approximation applies, further validating the approach in [6]. Compared to [6], we however find significantly stronger sensitivity by including the contributions away from the plasmon resonance. Most importantly, we confirm that the Migdal rate for a low threshold detector such as SENSEI, DAMIC or SuperCDMS is much higher than in noble liquid detectors. This is due to the lower ionization gap, the ω −4 scaling in (6), as well as the possibility of detecting all electronic excitations rather than only atomic ionizations.
Note added -In the final stages of preparing this work, [50] appeared, which studies the multiphonon response and has some overlap with the calculation in Appendix B.
Acknowledgements -We are grateful to Noah Kurinsky, Daniel Baxter, Gordan Krnjaic, Toby Opferkuch, and Chih-Pan Wu for useful discussions and to Yonatan Kahn for useful discussions and valuable comments on the draft. We further thank Diego Redigolo for useful discussions, proofreading the manuscript and collaboration on related work. TL and JK are supported by the Department of Energy under grant DE-SC0019195 and a UC Hellman fellowship. TL is also supported by an Alfred P. Sloan foundation fellowship.

Appendix A: Derivation of Migdal rate in semiconductors
The derivation of the Migdal effect for semiconductors is complicated by the spatial delocalization of the valence electrons. As a consequence, each valence electron feels the presence of a large number of ions and their fellow electrons in the crystal. The system is often described in the single electron approximation, given by the Hamiltonian where U (r) is an effective, periodic potential felt by the electron, due to the presence of the ions and the remaining electrons. In general U (r) is very complicated, and its eigenstates are typically obtained with specialized numerical methods in the realm of Density Functional Theory (DFT). For now, we can however keep U (r) as an abstract operator, and just work in terms of its eigenstates as long as possible. Concretely, the eigenstates of (A1) are Bloch wave functions where j indexes the electronic bands. Going forward, we use |k and |k] as shorthand for the full (ψ k ) and Bloch functions (u k ) respectively. To keep the notation manageable, the band indices j will be often be suppressed in what follows. We quantize the system over a finite volume V with periodic boundary conditions. The unit cell's volume is Ω and the number of cells in the crystal is therefore N = V /Ω. We will take the infinite volume limit only at the end of the calculation. This means that the sampling over the first Brillouin zone (BZ) will be a finite, discrete sum with N terms until we take the continuum limit by sending V, N → ∞.
To treat the ion-electron interaction, we must account for the screening from the spectator valence electrons, which is parameterized by , the frequency-dependent, microscopic dielectric function. The dielectric function of a material is defined by the relation with E and E ext the total and external electric fields, respectively. For a homogeneous material at distances much larger than the lattice spacing, the dielectric function is invariant under translations, and its Fourier transform is a function of a single momentum vector, such that we can write it as (q, ω). At short distances, however, there can be large oscillations of the electric field inside the primitive cell, and translation invariance is broken down to the discrete translation vectors of the crystal. Using this residual translation invariance, one can show that the Fourier transform of (A3) can be written as (A4) where k is in the first BZ (1BZ), and K and K are reciprocal lattice vectors. We will be using the common shorthand notation −1 KK (k, ω) ≡ −1 (k + K, k + K , ω). Note that −1 KK (k, ω) can be thought of as a matrix in the space of reciprocal lattice vectors, such that KK (k, ω) is defined as the matrix inverse of −1 KK (k, ω). Eventually we will drop terms for which K = K , as they tend to be numerically subleading, and this distinction will not be important.
The Migdal rate is given by the leading order expansion in both the dark matter nucleus interaction (H χ ) and the electromagnetic interaction (H e ), very analogously to a bremsstrahlung calculation. This is illustrated by the diagram in Fig. 3, which also serves to define the with b χ the DM-nucleus scattering length and r N , r χ and r the position operators for the dark matter, nucleus and electronic wave functions respectively. Here we have assumed sub-GeV DM, so that the scattering length is related to the DM-nucleus elastic cross section by σ N = 4πb 2 χ = A 2 σ n , where σ n is the DM-nucleon cross section and A is mass number. α is the electromagnetic fine structure constant, and we have included the charges of the tightly bound core electrons in the charge of the ion, denoted by Z ion . In order to simplify the derivation below, we will simply keep a constant Z ion , and only restore the momentum-dependent Z ion (k) in the final steps.
The Fourier transform of H e is where the sum over k samples the first Brillouin zone and K, K run over the reciprocal lattice vectors. We used the symmetrized dielectric matrix˜ −1 KK , which has the property that˜ −1 KK =˜ −1 K K . In what follows we will drop the˜, and the −1 KK will always refer to the symmetrized dielectric matrix.
There are two terms in the matrix element, corresponding to the advanced and retarded contributions, as in a standard bremsstrahlung calculation where the sum runs over the full momentum space for the ions. As laid out in the main text and Appendix B, our use of the impulse approximation means that we approximate both the final and the intermediate state nucleus wave function as plane waves, represented by |q N and |q respectively. The plane waves are normalized as r| q = e iq·r / √ V , etc. |ψ 0 indicates the ground state of the ion, which is in general not described as a plane wave, as it is bound in the crystal. We have neglected the small ground state energy in the propagator. We now calculate both terms separately.
First we calculate the DM-nucleus matrix element in (A8), accounting for the fact that the nucleus is bound in the crystal. Modeling both the incoming and outgoing the dark matter particle as a plane wave as well, the matrix element is with |0 the ground state of the free theory. The position space representation of the ground state wave function of a nucleus in a harmonic crystal is withω the oscillator frequency, averaged with respect to the density of states D(ω) and the thermal Bose factor, or ω ≡ ∞ 0 dω ω D(ω ) coth( ω 2T ). We prove equation (A11) in Appendix B. Evaluating the matrix element, one finds In the limit of a free nucleus, we recover as expected.
Next we calculate the electronic matrix element in (A8) by decomposing the electronic wave functions into the Bloch functions |p e ]: where in the second line we used the plane wave forms to evaluate the matrix element corresponding to the nucleus. Note that we use the Kronecker delta function for the ion wavefunction overlap here, and we will take δ q,q N +k +K → ((2π) 3 /V )δ(q N + k + K − q) in the continuum limit at the end. Next we can compute the Bloch amplitude, where we rewrite the amplitude as an integral just over the primitive cell, with volume Ω: where the R are the lattice vectors of the Bravais lattice. The Ω-subscript in (A20) serves as a reminder that this matrix element is only defined over the volume of the primitive cell. In the second step we used the periodicity of the Bloch functions under translations along the Bravais lattice (u j,p (r + R ) = u j,p (r)). In (A19) we furthermore relied on the identities e iR ·(k −k) = V Ω L δ k,k +L and e iR ·K = 1. Since both k and k are restricted to the first Brillouin zone, the only possible contributions to the sum come from L = 0 and from the L in the second Brillouin zone. However, the latter configurations in k and k are a set of measure zero in the N → ∞ limit, so we drop them from the calculation. Inserting (A20) back into (A16) results in and thus The calculation of M b is analogous: The H χ element in (A9) is trivial, as all incoming and outgoing states are plane waves Using the same methods as above, the H e matrix element is Adding both contributions in (A22) and (A25), the total matrix element is where we have dropped the negligible |k + K | 2 /2m N terms from the propagators. This is justified, as the form factor F exponentially suppresses contributions from terms with K p i − p f − q N . The cross section is obtained as usual by invoking Fermi's golden rule with v χ is the incoming DM velocity and f j (p e ) the temperature-dependent occupation numbers for the electronic bands, which are indexed by j, j . We have also introduced an integral over ω with a corresponding δ-function, such that we can continue to use the shorthand notation ω = ω j ,pe+k − ω j,pe . In the second line we rewrote the Pauli blocking factor by identifying the f j (p e ) as the Fermi-Dirac distribution with temperature T . CCD detectors such as DAMIC and SENSEI are run at T ≈ 140 K, while SuperCDMS is operated at even lower temperatures. This means that for ω ∼ eV we can safely neglect the e −ω/T factor in (A28). Next we take the infinite volume limit, which in momentum space corresponds to the continuum limit, by replacing 3 , etc. Plugging in (A26), the differential cross section is thus where the last line in (A29) equals Im[ KL ] in the random phase approximation (RPA). This follows from the generalization of Lindhard's formula [51] to the reciprocal lattice [52] 4 Note that the transition matrix elements are real for systems respecting time reversal and parity symmetry [53], such that the only contribution to Im[ KL ] is from the pole in the propagator. One can moreover show that The Im[− −1 KL ] object is a positive function and corresponds to the well known energy loss function (ELF), which governs the energy loss of charged particles passing through a material. To reduce the computational burden, we currently set Im[− −1 KL ] = 0 for K = L, although we verified that the combined, neglected contributions are at most comparable to the diagonal terms. With this additional assumption, we arrive at Note that in this expression, we have finally have restored the effective ion charge Z ion (|k + K|) associated with the momentum transfer from the ion to the electrons, k + K. In (A32), all information regarding the material and its properties is encoded in the reciprocal lattice vectors, the ELF and the form factor F (q), which is sensitive to the phonon density of states (see (A13)). Formulating the rate for the Migdal process in terms of the ELF and the phonon density of states has a number of advantages. Both can be measured directly for individual materials under various experimental conditions. A large body of experimental data is already available, though primarily for k ≈ 0 in the case of the ELF. It is therefore often necessary to extrapolate the experimental data on the ELF to finite k, which is can be done with Lindhard or Mermin oscillator models [54,55]. In fact, Ref. [16] already proposed to use optical data to estimate the rate for the Migdal effect in various materials. Our calculations put this approach on a more firm theoretical footing. In particular, [16] is relying on (7), which we have shown to be invalid for semiconductors. In addition, both the ELF and the density of states can be directly computed from first principles with well established, public codes such as GPAW and Quantum ESPRESSO. This has the advantage that no extrapolation to finite k is needed, though this method is substantially more computationally intensive. We explore both techniques further in a companion paper [49], but illustrate the advantage of this formalism in Fig. 4 which compares the ELF in Si computed in the DFT code GPAW to that measured from data. The ELF for k outside the 1BZ is defined as Im[− −1 KK (k, ω)] where k = k + K with k ∈ 1BZ and K a reciprocal lattice vector. The results shown are for momentum transfer along the [111] direction. The DFT calculation reproduces the data well, and comparison with data allows for more control over theoretical uncertainties in the Migdal rate calculation.
It is possible to further simplify the expression above by making a number of additional approximations, which provide some conceptual insight and simplify the numerical evaluation of the integrals in (A32). First, we observe that in most of the phase space q N · (k + K)/m N ω and |k + K| |q N |, which corresponds to the so-called soft limit. This allows us to factorize the various integrals as follows This is the formula we use to obtain the numerical results shown in Fig. 2. The free ion approximation is obtained by replacing the |F | 2 form factor with a δ-function and subsequently eliminating the p f integral  [56] and inelastic x-ray scattering measurements (k > 0) [57]. The DFT calculation reproduces the experimental results well. For low k, the plasmon peak is clearly visible, while it disappears for high k.
To make contact with earlier analytic work, we can further restrict the phase space to K = 0, such that this expression reduces to which matches the semi-classical result derived in [6] with fixed Z ion . Now recall that dσ dωE N ≈ dσ el dE N × dP dω in the soft limit, with and for sub-GeV dark matter we can take µ χn , µ χN → m χ . After carrying out the q N integral in (A35) we can read off the electronic excitation probability To make the connection with the Migdal effect in atoms in (9) as manifest as possible, we can substitute the matrix element back for Im[ (k, ω)] where (k, ω) is the macroscopic dielectric function.

Appendix B: Multiphonon response and impulse approximation
In this Appendix we compute the response associated with dark matter scattering off a nucleus which is embedded in an isotropic, harmonic crystal, to all orders in the phonon expansion. This is the generalization of the elastic nuclear recoil calculation, without the Migdal effect. Analogous to the analysis in [58], our calculation allows us to interpolate between the single phonon regime at low momentum transfers, and the free ion regime at high momentum transfers. Our goals hereby are two-fold: (i) We explicitly demonstrate that the impulse approximation is an excellent approximation to the full multiphonon response when the dimensionless parameter q/ √ 2m Nω 1, where q is the momentum imparted on the nucleus andω the typical frequency of the phonons.
(ii) We show that the response in this regime is well modeled by using the ground state of the harmonic system for the initial state of the nucleus, along with the plane wave approximation for its final state, justifying the approach in Appendix A.
For a short range interaction, the interaction between the nucleus and the DM is simply withṼ the Fourier transform of the potential. Here we are only concerned with establishing the momentum transfer for which the free ion approximation breaks down, rather than the overall normalization. We will therefore set the constant prefactorṼ 0 = 1 going forward. The dynamical structure factor describes the response of the lattice, and can be defined by Fermi's Golden rule: with λ i,f | the initial and final states respectively. The expectation value includes a thermal average over initial states. Upon Fourier transforming the energy δ-function, this can be written as Given that the λ i,f form a complete basis of eigenstates of the Hamiltonian, this can be expressed in terms of a time-dependent two-point correlation function where in the last line we dropped the explicit reference to the initial states λ i |. With Bloch's identity and a moderate amount of matrix algebra, the expectation value can be taken into the exponent (see e.g. Appendix B of [25] or chapter 9 of [59]) with W (q) ≡ 1 2 (q · r N (0)) 2 the Debye-Waller factor. In the small q limit, one can expand the timedependent exponential to extract the DM scattering rate to a single phonon [25,27] and two phonons [29]. Here we are however interested in the full momentum range, and cannot expand the exponential. Instead, we first compute the correlation function exactly for the special case of an isotropic harmonic oscillator with frequency ω 0 . In particular, the second quantization of the displacement operator is given by where the e j are a set of three unit vectors in position space and a (a † ) are creation (annihilation) operators. All directions are weighted equally since we assumed the isotropic limit. Here we also implicitly assumed that the equilibrium location of the ion is at the origin of the coordinate frame; more general expressions can be found in Appendix B of [25] or in [59]. With (B6) in hand, we can calculate with n(ω 0 ) ≡ 1/(e ω0/T −1) the Bose-Einstein distribution for a temperature T . The first and second term respectively correspond to the absorption and production of an oscillatory mode. This can be further rewritten as The Debye-Waller factor is obtained by setting t = 0 In a realistic crystal, there is of course a continuum of frequencies, rather than a single frequency ω 0 . 5 This can be accounted for by including a density of states D(ω), such that (B8) generalizes to and similarly for the Debye-Waller function. While this system is in general no longer integrable, the frequency and time integrals in (B5) and (B10) can be carried out numerically for a given D(ω). To illustrate the qualitative features of the transition to the free ion regime, we work with a simple, isotropic Debye density of states defined by where ω ph is the frequency of the acoustic phonons near the edge of the first Brillouin zone. 6 For Si and Ge, ω ph can taken to be roughly 40 meV and 25 meV respectively. In the high q regime, the response will be dominated by short time scales. An asymptotic expression can be derived through a steepest-descent expansion of (B5) around small t [58]. This corresponds to the impulse approximation, as the superscripts indicate whereω is the phonon frequency, averaged over the density of states and the thermal Bose-factor. (Negative frequencies correspond to the dark matter absorbing energy from the crystal.) When T = 0,ω = 3ω ph /4 for the density of states in (B11). In the free limit ofω → 0, or equivalently q → ∞, (B12) reduces to the familiar δfunction, enforcing the on-shell condition of an unbound nucleus. In Fig. 5, our numerical evaluation of (B5) is compared with the impulse approximation in (B12), for various values of the momentum transfer. For q/ √ 2m Nω = 0.5, the rate is dominated by the production of a single phonon, and the structure factor simply reflects the phonon density of states D(ω). Though the impulse approximation is clearly not valid, expanding the exponential factor in (B5) in low q is however an excellent approximation. This regime was studied in [24][25][26][27][28][29][30]. For q/ √ 2m Nω ≈ 1, a bump at higher frequencies starts to emerge, which corresponds to two phonon production. For even higher q, the peaks associated with single, double and triple phonon production are clearly visible, as well as a hint of the four phonon bump. Finally at q/ √ 2m Nω 3, the individual phonon peaks are no longer discernible, and the structure function has converged to the Gaussian form in (B12), associated with the impulse approximation. The slight offset of the maximum between both results is due to the fact that at high q, the true saddle point in (B5) is slightly offset from t = 0, an effect which can be corrected for analytically [58], though here we choose not to do so for simplicity. For even higher momentum, the central value of the Gaussian continues to shift to higher energies, as is apparent from the bottom right panel in Fig. 5. From (B12) one can furthermore see that the structure factor asymptotes to δ-function in the q → ∞ limit, correctly reproducing the free particle limit.
It then remains to be shown that the asymptotic form in (B12) is equivalent to using the plane wave approximation for the recoiling nucleus. We start from the position space wave function of a nucleus in an isotropic harmonic potential with characteristic frequency ω The ground state wave function of a harmonic crystal is the direct product of the ground states of a set of harmonic oscillator, weighted by the density of states of the crystal where we are now working in the T ≈ 0 limit. In the second line we used the definition (B14) and normalized the wave function such that ψ 0 |ψ 0 = 1. Now returning to the definition in (B2), we insert a set of plane waves |q N for the final states with E 0 the ground state energy. Plugging in the plane wave form for |q N and shifting the integration variable q N → q N + q, this yields √ 2ωmN , the structure factor resembles a Gaussian which is centered at the free-particle recoil energy q 2 /(2mN ) and with width ωq 2 /(2mN ).

S(q, ω) = (πm 3
Nω The impulse approximation is a good approximation when q √ m Nω , while the integrand is exponentially suppressed unless q N √ m Nω . We can therefore drop the q 2 N /2m N term from the δ-function, as well as the ground state energy E 0 . Carrying out the residual integrals then delivers the same expression as in (B12). We therefore conclude that at low T , the usage of plane wave wavefunctions for the excited states for the recoiling nucleus is justified in the limit q √ m Nω , as long as the bound state nature of the initial state is accounted for.
Appendix C: Details of rate calculation In this appendix, we provide expressions for the Migdal rate, dR/dω, in the free-particle and impulse approximations in the soft limit. We also illustrate the theoretical uncertainties arising from our various approximations. The differential rate is given by dR/dω = N T n χ d 3 v χ f (v χ ) v χ dσ/dω, where N T is number of target nuclei per kilogram and n χ = ρ χ /m χ is the DM number density. Consider first the free particle approximation in the soft limit. The differential cross-section entering the rate is then given by (A34). We approximate the ELF by defining Im[− (k , ω)] ≡ Im[− KK (k, ω)] where k ≡ |K + k| can lie outside the 1BZ and the overline corresponds to the angular average for fixed k and interpolation between the discrete k mesh points sampled by the DFT code GPAW. (We neglect the offdiagonal parts of the dielectric matrix.) The d 3 k integral can then be extended outside the 1BZ. In the soft limit, the k-dependence can be factorized and isolated into the quantity Due to this factorization, we can use the delta function to perform the angular q N integral. We then convert the q N integral to an integral over E N , the nuclear recoil energy, which can also be performed straightforwardly.
The result for the differential rate becomes where f (v) is the standard halo model for the DM velocity distribution in the Earth frame, and v min = 2ω µ χN , v max = v esc + v e are limits on the |v| integration. We have also defined the minimum and maximum nuclear recoil energies N cuts off the E N integration to account for a possible energy threshold. µ χN is the DM-nucleus reduced mass. The DM velocity integral is straightforward to perform numerically.
The corresponding derivation in the impulse approximation proceeds similarly, except we use (A33) for dσ/dω. In the soft limit, the k integral again factorizes into I(ω) given above. There is an additional integral over p f in this case due to the appearance of the form factor F rather than the usual momentum delta function. We shift the p f integration to be over q = p i − p f . We then first perform the angular q N integral, and use the energy delta function to perform the angular q integral (the azimuthal integrals are trivial). The |q| integration can then be performed, followed by the integral over q N , which we trade for an integral over E N . In doing so, we again cut off the E N integration at a threshold value E th N , which we will vary to illustrate the theoretical uncertainties associated with the impulse approximation. We could have instead swapped the order of q and q N integration and imposed a cutoff on q, but we find that the two procedures yield similar results for m χ 50 MeV. The DM velocity integral proceeds as before, only with v min = 2(ω + E th N )/m χ . The result for the rate is then In Fig. 5 we showed that the impulse approximation starts to break down for q/ √ 2ωm N 3. To estimate the size of the uncertainty associated with this approximation, we calculate the rate for both E th N = 4ω and E th N = 9ω, which is equivalent to restricting the phase space to respectively q/ √ 2ωm N > 2 and q/ √ 2ωm N > 3 in the free particle limit. The difference is shown by the shaded bands in Fig. 2. For m χ 50 MeV, the rate is dominated by the phase space corresponding to E th N < 9ω, as is evident by the diverging uncertainty bands. Here the impulse approximation ceases to be reliable and we chose not to extrapolate our results to this regime.
In Fig. 6 we show the differential rate for m χ = 100 MeV for Si and Ge under different cuts and assumptions. First, one observes that the free ion and impulse approximations are essentially identical for this mass point. To isolate the contribution from the plasmon pole in particular, we also plot the rate with a k < 2.5 keV cut, which reveals the plasmon enhancement around ω ≈ 20 eV. The separation between the blue and the green curves however shows that this contribution is highly subleading as compared to the high k, off-resonance part of the ELF. Finally, in orange we show the rate for the Migdal effect in atomic Si or Ge, using the results of [3]; this case corresponds to the collision of the DM with a single, isolated Si or Ge atom, and includes only atomic ionizations as a possible charge signal. Clearly, this gives a poor approximation for the rate in semiconductors, which is much larger due to the reduced energies for valence states, as well as the possibility of detecting all electronic excitations between bands. In contrast, for the atomic Migdal effect, transitions between bound electronic states do not give an ionization signal.
In Fig. 7, we show the impact of using momentumdependent Z ion (k) (the default for results in this paper) as opposed to constant Z ion . The momentum-dependent ion charges were obtained using the results of Ref. [32], which provides tabulated form factors for Si 4+ and Ge 4+ as a function of k/(4π) inÅ −1 , where Z ion (k) starts deviating from Z ion = 4 at k ∼ 5 − 10 keV. The rate with Germanium, m χ = 100 MeV, σ n = 10 −38 cm 2 Impulse Impulse, k < 2.5 keV Free Free, k < 2.5 keV Atomic target FIG. 6. Differential Migdal ionization rate in Si (left) and Ge (right) for mχ = 100 MeV and σn = 10 −38 cm 2 . The total rate in the impulse (free particle) approximation with EN > 4ω is given by the solid (dashed) curves. The corresponding result cutting off the k integral at 2.5 keV is shown in light blue, where the plasmon resonance (which is only present for small k) is clearly visible. For comparison, the dotted curves show the atomic ionization rate due to the Migdal effect, computed in Ref. [3] for Si and Ge atoms. momentum-dependent charges is larger by a factor of 3 for Ge and 1.6 for Si, where the effect is larger for Ge because the semi-core electrons are not as tightly bound. The difference with constant Z ion also becomes larger at higher ω, where the differential rate is also weighted at higher momentum transfers.
The first condition v N t r Bohr justifies the expansion in small R N compared to the typical extent of electronic wavefunctions, while the second condition r Bohr t justifies the use of the electrostatic potential.
For small t, we can expand the potential about small R N (t): which is just the dipole potential for a recoiling nucleus with dipole moment Z N R N (t). In time-dependent perturbation theory, we can now compute the transition probability as where ω = E f − E i and η is a small positive number that accounts for dissipation and turns off the potential at large times. In the last line, we took η → 0. This is precisely the matrix element we obtained in the main text by rewriting (7) in terms of (8) with operator identities. This gives a physical interpretation for the Migdal effect that is more difficult to see from the boosting argument. To further make contact with our semiconductor calculation, note that we can rewrite Then the transition probability P i→f for single ionizations is given by and we see the same form of the matrix element as in the semiconductor calculation, (A26), up to overall factors that account for the differences in the physical systems and for the DM-nucleus scattering rate. The situation is analogous to bremsstrahlung, which can be computed either in a semi-classical approximation (as in this section) or with quantum mechanics (as we did in the semiconductor case, in Appendix A). Summing over initial and final states and taking the soft limit v N · k ω of this transition probability, we can then obtain (9).