Dark Matter Detection With Bound Nuclear Targets: The Poisson Phonon Tail

Dark matter (DM) scattering with nuclei in solid-state systems may produce elastic nuclear recoil at high energies and single-phonon excitation at low energies. When the dark matter momentum is comparable to the momentum spread of nuclei bound in a lattice, $q_0 = \sqrt{2 m_N \omega_0}$ where $m_N$ is the mass of the nucleus and $\omega_0$ is the optical phonon energy, an intermediate scattering regime characterized by multi-phonon excitations emerges. We study a greatly simplified model of a single nucleus in a harmonic potential and show that, while the mean energy deposited for a given momentum transfer $q$ is equal to the elastic value $q^2/(2m_N)$, the phonon occupation number follows a Poisson distribution and thus the energy spread is $\Delta E = q\sqrt{\omega_0/(2m_N)}$. This observation suggests that low-threshold calorimetric detectors may have significantly increased sensitivity to sub-GeV DM compared to the expectation from elastic scattering, even when the energy threshold is above the single-phonon energy, by exploiting the tail of the Poisson distribution for phonons above the elastic energy. We use a simple model of electronic excitations to argue that this multi-phonon signal will also accompany ionization signals induced from DM-electron scattering or the Migdal effect. In well-motivated models where DM couples to a heavy, kinetically-mixed dark photon, we show that these signals can probe experimental milestones for cosmological DM production via thermal freeze-out, including the thermal target for Majorana fermion DM.

The strategy of searching for dark matter (DM) by detecting the energy imparted to a nucleus during a scattering process dates back to the original proposals of Goodman and Witten [1] and Drukier, Freese, and Spergel [2] and forms the basis for the multi-ton liquid noble experiments XENONnT [3] and LZ [4]. In recent years, a vast landscape of plausible DM models has opened up at much lighter DM masses compared to the GeV-scale WIMPs these experiments were originally designed to search for (see [5] and references therein), necessitating a rethinking of the basic assumptions that went into the calculations of the nuclear scattering rates. In particular, for elastic scattering of free nuclei, the energy transfer between the DM and nucleus is inefficient if the DM mass is much below the nucleus mass: for a xenon nucleus of mass m N = 122 GeV and DM of mass m χ = 100 MeV, the typical momentum transfer is q ∼ m χ v ∼ 150 keV, and the recoiling nucleus has energy q 2 /(2m N ) ∼ 0.1 eV, well below the thresholds of typical noble liquid detectors.
However, as emphasized in recent work [6,7], a nucleus in a lattice is not a free particle, and therefore the elastic scattering argument does not necessarily apply for solid-state detectors. Indeed, an individual nucleus in a solid is subject to a harmonic oscillator potential from its neighboring atoms where ω 0 ∼ 60 meV is a typical optical phonon energy for silicon (for other materials, ω 0 ∈ [20 meV, 140 meV]).
This potential is short-range, acting only over interatomic distances a ∼ (keV) −1 , and thus the displacement energy required to remove a nucleus from its lattice site is E d ∼ 1 2 m N ω 2 0 a 2 ∼ O(10 eV). Only for recoil energies E R well above E d may the final-state nucleus be wellapproximated as a free particle; below E d , the final-state energy spectrum is that of a harmonic oscillator, and of course the initial state is the ground state of a harmonic oscillator rather than a zero-momentum plane wave. This observation has previously been made in the context of single-phonon excitations in solids [8,9], where DM with de Broglie wavelength exceeding the interatomic spacing interacts primarily with collective modes of many oscillating nuclei, which are quantized into a phonon spectrum.
In this Letter, we study the intermediate regime where the momentum transfer satisfies q 1/a, such that the interaction is localized to a single lattice site, but E R E d so that both the initial and final nuclear states belong to the harmonic oscillator spectrum. We construct a toy model of a solid-state detector by considering the non-relativistic single-particle quantum mechanics of a nucleus subject to the potential (1). The characteristic momentum spread of the ground state is where q 0 ≈ 56 keV for silicon. We will show that when a momentum transfer q > q 0 is kinematically allowed, the typical energy deposited is This is the energy expected from elastic scattering, but here it may be interpreted as the production of a multiphonon state with occupation number n = (q/q 0 ) 2 . We will show that in this greatly simplified model of a single oscillator, the distribution of phonon number is exactly Poissonian, so that unlike elastic scattering where only a single energy E is allowed for a given q, there is a spread ∆n = q/q 0 corresponding to an energy spread This observation extends the reach of solid-state nuclear recoil detectors to sub-GeV DM compared to previous analyses which assumed the kinematics of elastic scattering. As a cross-check, we show that our results parametrically reduce to the previously-calculated single-phonon rate when q q 0 , offering a pleasing interpretation of the single-phonon rate as an upward Poisson fluctuation when the elastic energy is well below the phonon energy; our results converge exactly on the expected elastic recoil spectrum when q q 0 . Our work complements that of [6,7] by explicitly showing the transition between the single-phonon and elastic regimes in a simple model (see also [10] which considers production of multiple acoustic phonons for q q 0 ). We will show that when the DMnucleus interaction is mediated by a heavy dark photon, the full single-and multi-phonon spectrum can probe a wide range of thermal relic targets, including Majorana DM which suffers a velocity-suppressed cross section.
Furthermore, we argue that when a nuclear recoil induces secondary ionization as in the Migdal effect [11][12][13][14][15][16][17][18][19] (see Fig. 1 for an illustration), the electron spectrum is independent of the nuclear recoil spectrum except at the very largest kinematically-allowed ionization energies. This result was previously known for isolated atoms [15], but we show that it persists in a simple model of a bound nucleus. For similar reasons, we also find that DM-electron scattering also yields secondary phonon excitations. Therefore, the multi-phonon spectrum we compute should be an irreducible component of an ionization signal in solid state detectors, and due to the Poisson tail it may lie above the threshold of next-generation calorimetric detectors and would help to distinguish it from other non-particle backgrounds that create charge pairs without accompanying recoil. Finally, our results suggest a qualitatively different figure of merit for detector materials to detect sub-GeV DM-nuclear scattering: to maximize the sensitivity at low DM masses, large phonon energies are preferred to maximize the Poisson fluctuations when the elastic energy deposit is close to threshold.

Harmonic Oscillator Model
We consider a single nucleus of mass m N subject to a 3-dimensional isotropic harmonic oscillator potential: The energy eigenstates are | n , where n = {n x , n y , n z }, with energies E n = n + 1 2 ω 0 , where n = n x + n y + n z . To model galactic DM particles χ scattering from this nucleus, we follow [20] and write the total scattering rate for a detector with N T nuclear targets as where ρ χ /m χ = 0.3 cm −3 (GeV/m χ ) is the local DM number density [21], Z is the atomic number of the target and σ p is a fiducial DM-proton cross section; for concreteness we restrict our attention to a contact interaction between DM and protons, but our results are easily generalizable to other interactions. We do not include a nuclear form factor because the momentum transfers we consider for sub-GeV DM will always be smaller than the inverse nuclear size, q 1/R 0 ∼ O(10 MeV). The dimensionless functions appearing in the rate integral are the inelastic form factor and the DM inverse mean speed where f χ (v) is the DM velocity distribution, which we assume to be the isotropic Standard Halo Model. Here is the minimum DM speed needed to excite the nucleus to harmonic oscillator level n. We can calculate the form factor analytically using either momentum-space wavefunctions or the harmonic oscillator algebra. Due to the isotropy of the harmonic oscillator, the angular integrals can also be performed analytically. Defining an angular average |f (n, q)| 2 ≡ 1 4π dΩ q |f (n, q)| 2 , we find While to our knowledge this result has not appeared before in the DM literature, it is known in the case of neutron scattering, see for example [22]. Performing the angular average in Eq. (6) and using in Eq. (10) yields a differential spectrum (11) and in Fig. 2 we show some representative spectra for different DM masses. See the Supplementary Material for derivations of the main results in this section. The form factor in (10) is nothing but a Poisson distribution in n with meann = q 2 /q 2 0 = q 2 /(2m N ω 0 ). Therefore, the mean energy deposited in a DM-nuclear scattering isnω 0 = q 2 /(2m N ), the elastic value. In the limit q 2 q 2 0 , the Poisson distribution approaches a delta function δ(n −n); taking the continuum limit n → E R /ω 0 , n → dn in Eq. (6), and replacing qdq = m N dE R , we recover the usual elastic recoil spectrum with v . This matching justifies the standard approximation of treating the nucleus as a free particle with the elastic dispersion relation, at least for ω 0 E R < E d such that the nucleus remains bound in the harmonic oscillator potential.
Consider now the opposite limit, q 2 q 2 0 . From [6][7][8][9], the form factor for production of a single optical phonon is parametrically (q/q 0 ) 2 e −q 2 /q 2 0 , where in the context of condensed matter physics the exponential is known as the (zero-temperature) Debye-Waller factor. We can now interpret this as the Poisson probability for n = 1 phonons whenn 1. While the highest-probability outcome is producing no phonons, n = 0, the most likely excitation above the ground state is n = 1, with larger phonon numbers strongly suppressed by powers of q 2 /q 2 0 1. The advantage of our model is that we can now seamlessly interpolate between the single-phonon regime q 2 q 2 0 and the elastic regime q 2 q 2 0 . When q 2 ∼ q 2 0 , as is the case for kinematics of 50 MeV DM in the Standard Halo Model, the Poissonian fluctuations in phonon number become important. In the right panel of Fig. 2 we show representative quantized spectra that illustrate this behavior. Indeed, consider a detector with threshold close to ω 0 . If the scattering were purely elastic, DM with maximum momentum q 0 could not produce a detectable nuclear recoil. However, for q = q 0 , the Poisson probability for n = 2 phonons is 0.18, compared to 0.37 for n = 1. Thus the true reach of the detector extends to lower DM masses, because the probability to deposit energy well above threshold is comparable to the probability to deposit energy at threshold. A similar version of this argument applies to detectors with thresholds some-  6) with n ≤ E thresh /ω0 (solid), along with the equivalent curves for elastic scattering (long dashed) and single-phonon production (short dashed magenta). The gray shaded regions represent the accelerator bounds on invisibly decaying A from the LSND [23], NA64 [24], E137 [25,26], BABAR [27][28][29], and MiniBooNE [30] experiments. Also shown are limits from XENON1T [31] (limits on electron scattering interpreted as σp in the dark photon model), CRESST II [32], and CRESST III [33]. We also show the neutrino floor background from [5].
what above the optical phonon energy: the Poisson tail of events with large n permits sensitivity to smaller DM masses than would be expected based on elastic kinematics alone.

Nuclear Scattering Reach
Now we apply our observations to estimate the reach of a next-generation calorimetric detector whose sensitivity to light DM can be greatly enhanced by including irreducible multi-phonon contributions from the Poisson tail. For concreteness, we consider a DM-nucleus interaction mediated by a massive, kinetically-mixed dark photon A with Lagrangian where is a kinetic mixing parameter, g D is the A -χ coupling constant, J µ EM is the electromagnetic current, and J µ D is the dark matter current. In the contact limit m A q, for complex scalar χ, the single proton cross section is which enters into Eq. (6). For m A > m χ , this model can realize thermal freeze-out via χχ * → f + f − annihilation where f is a Standard Model fermion, and there is a oneto-one correspondence between σ p and the early universe annihilation rate [34,35]. In Fig. 3 we show how including the additional phonons from the Poisson tail can enhance the DM signal yield for detectors with various energy thresholds, from 1 eV down to the single-phonon energy ω 0 . Each solid curve represents a different threshold for a silicon crystal target. We only include transitions into bound oscillator final states 0 → n with nω 0 < E d , so this simple model represents a lower bound on the total signal rate. The dashed curves show the would-be sensitivity of a low-threshold detector in the elastic regime, demonstrating that our result continuously interpolates between the discrete phonon regime and elastic DM-nuclear scattering which takes over forn 1, and that multi-phonon production provides increased sensitivity at lower DM masses. Note that, as in [7], the dashed elastic curves in this figure flatten towards higher masses reflecting the µ χp dependence in Eq. (12); the inelastic phonon curves shift upwards as a greater fraction of halo particles can deposit energy above E d to displace the nucleus, which we do not consider in our model.
Here, the A couples universally to charge, so we also include an atomic form factor in Eq. (6) to account for screening [36][37][38]: where λ TF ≈ 0.89/Z 1/3 a 0 ≈ 0.37a 0 is the Thomas-Fermi screening length for silicon and a 0 is the Bohr radius. At low momentum q 1/λ TF , |F A (q)| 2 → 0, reflecting complete charge screening by the neutral atom. This is also a manifestation of the well-known fact that dark photons do not couple efficiently to optical phonons at low momentum in non-polar crystals, since out-of-phase oscillations are suppressed compared to in-phase oscillations by powers of q. However, because λ TF q 0 ≈ 5.2 in silicon, |F A (q 0 )| 2 = 0.93, so this screening only affects the kinematic regime with q q 0 . In practice, screening slightly suppresses the single-phonon rate and the total rate below m χ 10 MeV compared to a generic heavy mediator coupling only to protons. For larger DM masses, our single-phonon curve roughly reproduces the analysis of [7] for a heavy hadrophilic mediator, which is an important check on the validity of this simple model.
The blue curves in Fig. 3 correspond to thermal freezeout targets for complex scalar and Majorana fermion DM candidates coupled to A . The Majorana cross section is proportional to σ p in Eq. (14) (which is defined for scalar χ), and further suppressed by v 2 ∼ 10 −6 , where v is the DM velocity in the earth frame (see [35] for details). Both models here feature p-wave DM annihilation and are, therefore, safe from CMB bounds, which exclude freeze-out for s-wave candidates with m χ 10 GeV [39]. It is notable that calorimeters with sufficiently low thresholds in the few-phonon range may soon begin to explore the Majorana thermal target.

Migdal Effect and Electron Scattering
Nearly a century ago it was pointed out by Migdal [11] that a sudden impulse delivered to the nucleus could result in electronic transitions in atoms. In recent years this observation has become increasingly relevant for light DM searches because the ionization energy may greatly exceed the elastic energy (for calorimetric detectors) or quenched elastic energy (for ionization detectors) from nuclear recoil, improving the possibility of detecting nuclear recoil events through secondary ionization even if the nuclear recoil energy is below the detector threshold.
However, as noted in [18], the theoretical formalism of the Migdal effect only applies to isolated atoms, where the energy eigenstates are free-particle plane waves and thus a unitary transformation may be applied to transform to the frame of the recoiling nucleus. The overlap between this frame and the stationary frame leads to a factor of exp(iq e · r e ) in matrix elements between electronic states, where q e ≡ (m e /m N )q is suppressed from the physical momentum transfer q by the ratio of electron to nuclear masses. Despite the rather limited validity of the formalism, theoretical estimates of Migdal rates have been performed for semiconductors [15], including for the valence bands [19], where the effective mass of the electron may differ considerably from the vacuum mass m e used in the standard Migdal calculation.
While we do not attempt to provide a rigorous derivation of the Migdal effect in solid-state systems in this paper (especially because many-body electron states such as the plasmon may contribute significantly [40,41]), we consider a simplified model in the spirit of our nuclear scattering analysis, a nucleus in a harmonic potential which exerts a potential V e on a single electron: Transforming to relative and center-of-mass coordinates, r ≡r N −r e ,R = (m NrN +m ere )/(m N +m e ), the Hamiltonian becomes separable up to small perturbations (see the Supplementary Material), and the eigenstates can be written |Ψ = | n; ψ e , where the first label is the harmonic oscillator state and ψ e (r) is the electronic wavefunction. The generalized transition amplitudes from Eq. (7) now factorize into oscillator and ionization terms where q e ≡ (m e /m N )q is the familiar Migdal ionization factor [15]. If the ionization matrix element depends only on the magnitude |q e |, the differential scattering rate factorizes where f (n, q) is the Poisson distribution in Eq. (10) and is the ionization probability as in [15]; if Z f does not factorize, the integration is a convolution of the ionization and oscillator form factors and may deviate from the Poisson form (see Supplementary Material). We see that in this model, the electronic spectrum is unmodified to leading order in m e /m N compared to the free-nucleus picture, and the only changes come in replacing the initial and final nuclear states with harmonic oscillator states, the matrix elements of which we have already calculated. Since the matrix element factorizes, the only coupling between the electronic and nuclear excitation energies in Eq. (18) comes from energy conservation, which after integrating over the DM velocity distribution modifies v As our earlier analysis has shown, the mean nuclear energy is the elastic energy, which is typically much less than the electronic excitation energy, and thus the only effect of the harmonic oscillator spectrum in this model is to truncate the electronic spectrum at slightly smaller energies than would be expected from elastic nuclear scattering due to the Poisson tail. Compared to previous results on isolated atoms [15], our new result is the Poisson spectrum of phonons which replaces the continuum of elastic recoil energies for a free nucleus. On the other hand, for DM scattering directly off an electron [20,42], the matrix element is proportional to Ψ |e iq·re |Ψ , which now involes r e = R − µ me r, where µ ≡ m e m N /(m e + m N ) ≈ m e . In the limit m N m e , the above analysis yields where the electron matrix element, ψ e |e −iq·r |ψ e , is familiar from previous analyses of DM-electron scattering, but the nuclear matrix element is identical to that of the the Migdal effect. Note that unlike in the case of the Migdal effect, since the coefficient of r in r e is unity up to O(m e /m N ), the electron scattering matrix element is independent of the vacuum electron mass to leading order and can be computed entirely with band-structure wavefunctions which fully take into account the effective mass.
Our arguments about electronic excitations are fairly general but rely on an electronic potential which is a function only of r N − r e and a system which is separable enough that it can be treated as a two-body problem and decomposed into relative and center-of-mass coordinates. At this point it is unclear to us how well these assumptions hold in realistic solid-state systems (at the very least, this analysis completely ignores electron-phonon interactions and anisotropies due to the lattice structure which are likely to couple the two spectra at some level beyond simple kinematics), but regardless, our analysis of the nuclear scattering matrix element suggests that the Poisson tail may help push the nuclear recoil energies above threshold. A detector capable of converting ionization energy to phonons, such as the low-threshold calorimetric detectors used by SuperCDMS [43][44][45] and EDEL-WEISS [46,47], could observe the Poisson phonon spectrum of nuclear recoil simultaneously with the electronhole pairs created by ionization, which could help distinguish between true scattering events and low-momentumtransfer background processes such as charge leakage.

Conclusions
In this Letter we have constructed a simple quantummechanical model which describes the multi-phonon regime of DM-nuclear scattering in solid-state systems, interpolating between the elastic regime and singlephonon production when the energy deposit is less than the displacement energy E d . Our key finding is that when the elastic energy q 2 /2m N is close to the optical phonon energy ω 0 , there is an order-1 variance in the number of phonons produced, such that the sensitivity of low-threshold detectors to low-mass DM is stronger than previously expected from elastic kinematics alone. This observation suggests that detectors with an especially large optical phonon energy, such as diamond [48], may be able to take advantage of these Poisson fluctuations even when the detector threshold is somewhat above the single-phonon energy.
As an illustrative example, we have shown that plausible next-generation calorimeter detectors with ∼ eV scale thresholds can exploit this irreducible effect to greatly enhance their sensitivity to sub-GeV DM. Remarkably, the potential gains identified here could enable such detectors to probe the full thermal relic parameter space for Majorana DM candidates with mass between 10 MeV and 1 GeV freezing out via kinetically-mixed dark photons, corresponding to parameter space which is far below current direct detection limits.
Finally, according to [22], a more realistic model would incorporate a nontrivial phonon density of states and anharmonicities which broaden the phonon spectrum, but we expect that the parametric scaling of our result would persist, at least for monoatomic crystals. Crystals with more than one atom per unit cell, in particular polar materials where dark photons couple to optical phonons even at low q [6-9], may exhibit similar behavior but deserve a dedicated analysis, especially because of their directional detection capabilities. We emphasize here that the rates we have computed are only lower bounds on the total nuclear scattering rates for sub-GeV DM. Indeed, DM heavier than about 10 MeV has sufficient kinetic energy to displace a nucleus from its lattice site, and the fact that neither the initial nor the final states are free plane waves may allow for the possibility of inelastic scattering when the elastic rate below threshold is zero. This would give a spectrum of recoil events with E N > E d additive to the one we consider here. We plan to investigate this possibility in future work.

Dark Matter Direct Detection With Bound Nuclear Targets: The Poisson Phonon Tail Supplementary Material
Yonatan Kahn, Gordan Krnjaic, Bashi Mandava

DERIVING THE SCATTERING RATE
The transition rate between nuclear states induced by DM-nuclear scattering can be computed using Fermi's Golden Rule. Following [20], the scattering cross section times velocity for the inelastic χ(p)N (0) → χ(p )N ( n) transition between individual oscillator levels 0 → ( n) can be written in relativistic normalization as where primes denote final-state quantities. To calculate the momentum-space matrix element M(q), we begin by postulating a contact interaction between dark matter and nuclei of the form where r χ , r N are the DM and nuclear coordinates, respectively, and M N is the dimensionless relativistic scattering matrix element between DM and a free nucleus, which is simply a constant for a contact interaction. For simplicity we consider the case of scattering mediated through a heavy dark photon, such that the fundamental interaction is between DM and protons; the matrix elements are related by From this we can define the fiducial single proton cross section Using relativistic state normalization, the initial/final state wave functions in position space are where we treat the dark matter as a free-particle plane wave and the ψ i (r N ) are harmonic oscillator wavefunctions. The matrix element can then be written where q ≡ p − p is the momentum transferred from the dark matter to the harmonic oscillator system. The harmonic oscillator matrix element can be written in momentum space as follows: whereφ n are harmonic oscillator wave functions in momentum space. Note that, for comparison with the derivation in [20], here we adopt a normalization convention for which Putting it all together, we insert M(q) back into Eq. (S1) and replace E χ , E χ → m χ , E N , E N → m N in the nonrelativistic limit. Squaring and summing over all oscillator levels such that n x + n y + n z = n yields the form factor |f (q, n)| 2 defined in Eq. (7). Summing over all allowed energies (labeled by n) and scattering targets and integrating over the DM velocity distribution, we obtain the usual rate formula [20] stated in Eq. (6): is the inverse mean speed and v min is the minimum DM velocity required to upscatter the nucleus into a state with energy E n = (n + 1/2)ω 0 . Note that our derivation here is essentially identical to that of [20], replacing DM scattering from bound electrons with DM scattering from bound nuclei. The only essential difference is the factor of Z 2 arising from coherent scattering over all protons in the nucleus.

A. Momentum Space Wavefunctions
We begin by constructing the harmonic oscillator wavefunctions in momentum space, which is convenient for evaluating the matrix element where the argument of one wavefunction is translated by the momentum transfer q. In momentum space, the position and momentum operators are represented asp → p,x → i∇ p , so the time independent Schrödinger equation with = 1 becomes which we have written in terms of separable solutions where denotes differentiation with respect top. For each i = x, y, z, the normalized solutions satisfy where H n is an Hermite polynomial. In particular, the ground state has a Gaussian profile, with momentum spread of order q 0 = √ 2m N ω 0 as claimed. Note that the normalized wavefunctions (S12) satisfy the usual non-relativistic normalization convention which differs from the convention in Eq. (S9) byφ n = (2π) 3/2 φ n which is more common in relativistic treatments. Throughout this paper (including the remaining Supplementary Material), we use the convention in Eq. (S14).

B. Poisson Distribution
According to Eq. (S8), we may write the matrix element between harmonic oscillator states in momentum space as n|e iq·r N |0 = dp x dp y dp z φ * n (p)φ 0 (p − q) . (S15) Since the wave functions are separable in Cartesian coordinates, we need only evaluate the integral for a single component. For ease of notation, we will do this for the x-coordinate and write n ≡ n x : Defining the dimensionless variables a ≡ p x / √ m N ω 0 and b ≡ q x / √ m N ω 0 , the right-hand side of Eq. (S16) is so our task reduces to evaluating the expression (S18) Note that the integrand is related to the generating functions for Hermite polynomials: By performing a summation over all I n (b) as follows: and using the Taylor expansion of exp(bt) on the right hand side, we can read off the expression for individual I n (b) by matching terms of equal n in the summation: Using this result to evaluate our original Cartesian integral, we obtain (now restoring the index n x ) Since the integrals in Eq. (S15) are identical in p x , p y , p z , we can use Eq. (S12) to write their product as where n = n x + n y + n z and we have written the momentum components q i in spherical coordinates for future convenience. Squaring and summing over degenerate states, the form factor for exciting the n th energy level of the harmonic oscillator is given by which resembles a Poisson distribution up to the angular factors. We now turn to performing an angular average of this expression to justify the full Poissonian form shown in Eq. (10).

C. Angular Integrals
Since phonon directions are not observable, we are interested in computing the angular average 1 4π dΩ q |f (n, q)| 2 , summing over all occupation numbers that satisfy the constraint n = n x + n y + n x . To average the angular part of Eq. (S24), we need to compute Let (sin θ cos φ) 2 = A, (sin θ sin φ) 2 = B and (cos θ) 2 = C. The integrand with the summation is Multiplying and dividing by n!, we have where the summation is simply the multinomial expansion, With this reorganization of terms, our integrand of interest simplifies considerably: Using this form, Eq. (S25) becomes as claimed in Eq. (10).

D. Operator Algebra
We can obtain the same result by using the algebra of creation and annihilation operators. The three-dimensional harmonic oscillator can be separated into three mutually-commuting sets of creation and annihilation operators, in terms of which the position and momentum operatorsx i andp i can be written Since q ·x N = q xxN + q yŷN + q zẑN and the spatial operators commute with each other, without loss of generality we can simply compute and copy the result for n y and n z ; multiplying these together gives the desired result for arbitrary n. First let's write the exponential operator in terms of creation and annihilation operators: Identifying A →â † and B →â, and noting that [â,â † ] = 1 is a c-number, the series in Eq. (S34) truncates after the third term. Indeed, the third term is just a number, so we have the exact result and since the last factor is just a c-number, this can be written as Note that this is the same argument presented in [6] to compute the amplitude for single-phonon production, only restricted here to the greatly simplified context of a single 3-dimensional oscillator.
Since we are interested in simplifying e iqxx N = e iκ(â+â † ) , where κ = iq x / √ 2m N ω 0 , the relevant commutator is [κa, κa † ] = κ 2 . Following the above argument yields Consider taking the matrix element of this operator between the states 0| and |n x . Acting on |n x on the left, we have to get to the state |0 to have a nonzero matrix element with 0|. The only way to get there is to act n x times withâ and zero times withâ † . This means we take the n th x term from the right-most exponential series, and the 0 th term from the middle exponential. Since these operators satisfŷ a n x |n x = n x !|0 , we can act on our initial and final states to obtain Taking the modulus squared of this expression gives and multiplying identical expression for the n y and n z contributions recovers (S24). In fact, from Eq. (S41) we can see that the Cartesian occupation numbers n x , n y , and n z are also Poisson-distributed, so their sum n = n x + n y + n z will also be Poisson, confirming our more detailed calculation.

GENERALIZING THE MIGDAL EFFECT
Dark matter scattering off bound nuclear targets can also yield electronic energy in the form of Migdal ionization. However, unlike previous studies of this effect [15], here the nucleus is not a free particle, so here we revisit and generalize this result with a harmonic oscillator dispersion relation for the nuclear target. Rather than considering the final state in the boosted frame of the recoiling nucleus, we will decompose the problem into relative coordinates in the lab frame, which more easily generalizes for a bound nucleus. 1

A. Choosing Coordinates
We begin by first considering a simple atomic system in which a single electron and nucleus are bound by a potential V e , and the nucleus is held in place by a harmonic oscillator potential. The Hamiltonian for this system can be written (S42) 1 We thank Gordon Baym for suggesting this perspective on the problem.
Transforming to relative and center-of-mass coordinates, the Hamiltonian becomesĤ =Ĥ 0 + ∆Ĥ, wherê where µ = m e m N (m e + m N ) ≈ m e is the electron-nucleus reduced mass, and expanding the original harmonic oscillator term gives ∆Ĥ = −µω 2 0R ·r + Written in this way,Ĥ 0 is separable and can be solved by Ψ(R, r) = ψ N (R)ψ e (r), where ψ N is a simple harmonic oscillator wavefunction for the nucleus N and ψ e an electronic wavefunction. The terms in ∆Ĥ are suppressed by powers of m e /m N 1 and can be treated as small perturbations. In particular, the first-order energy shift is parametrically where we have assumed that the typical spread in position space of the electronic wavefunction is of order the lattice spacing a; note that ∆E independent of the harmonic oscillator level N because R = 0 in any stationary state, so only the second term in ∆Ĥ contributes. For ω 0 ∼ 50 meV and m N = 28 GeV as for silicon, we have ∆E ∼ 25 neV ω 0 and thus we are justified in ignoring the perturbation and treating the nuclear spectrum as purely a harmonic oscillator spectrum.

B. Including Dark Matter
Equipped with this formalism, we can now include the contact interaction from Eq. (S2) that couples the DM to the nucleus and repeat the argument that culminates in Eq. (S10) in Sec. 1 of this supplement with the initial/final states |Ψ = 2m χ √ 2m N ψ 0 (r N )ψ e (r e ) e ip·rχ (S47) |Ψ = 2m χ √ 2m N ψ n (r N )ψ e (r e ) e ip ·rχ , where we have merely extended Eq. (S5) to include electron wave functions which have non-relativistic normalization to match the conventions of [15]. Since the DM potential is only a function of r N − r χ , all the steps leading up to Eq. (S7) are identical and the ψ e states are spectators up until the last step where, instead of n|e iq·r N |0 , we get Ψ |e iq·r N |Ψ = n|e iq·R |0 ψ e |e iqe·r |ψ e , q e ≡ m e m N q .
Upon squaring this result, the first factor recovers Eq. (S24) and the second factor is the same ionization probability found in Ref. [15]: Although here we have only considered single-electron atoms, the calculation straightforwardly generalizes to atomic systems with multiple electrons at relative coordinates r j : Z f (q e ) → j d 3 r j ψ * e (r j ) exp i j q e · r j ψ e (r j ).
Thus, including the electronic matrix elements, the total scattering rate for the Migdal effect now becomes where the sum on f includes all allowed electron final states and the additional electronic energy E e released in this inelastic process shifts the minimum velocity required to undergo a given transition v (n,e) min = E e + nω 0 q + q 2m χ . (S53) For systems with spherically symmetric ionization probabilities, Z f depends only on E e (as in [15] with free atoms), so the electron term factorizes from the 1 4π dΩ q angular average that yields the Poisson distribution in Eq. (10). In this case, the differential scattering rate becomes which recovers the expression in Eq. (18). However, in general, it need not be the case that these terms factorize (in particular, the valence electrons in semiconductors are not necessarily in spherically-symmetric states) and the final distribution will involve an integration over the combined oscillator and ionization probabilities. For a detailed discussion of the Migdal effect in semiconductors see [49].