Dark matter direct detection from the single phonon to the nuclear recoil regime

In most direct detection experiments, the free nuclear recoil description of dark matter scattering breaks down for masses $\lesssim$ 100 MeV, or when the recoil energy is comparable to a few times the typical phonon energy. For dark matter lighter than 1 MeV, scattering via excitation of a single phonon dominates and has been computed previously, but for the intermediate mass range or higher detector thresholds, multiphonon processes dominate. We perform the first calculation of the scattering rate via multiphonon production for the entire keV-GeV dark matter mass range, assuming a harmonic crystal target. We provide an analytic description that connects the single phonon, multiphonon, and the nuclear recoil regimes. Our results are implemented in the public package $\texttt{DarkELF}$.


I. INTRODUCTION
The effort to directly detect dark matter (DM) is entering the sub-GeV mass regime, thanks to experimental innovations which allow for ever lower energy thresholds. For kinematic reasons, this regime is especially challenging for DM which primarily couples to hadronic matter. For a DM mass (m χ ) below 1 GeV, the energy that the DM can deposit in an elastic collision with a nucleus of mass m N is bounded by For m χ ≪ m N this is only a small fraction of the total available DM kinetic energy, which can make it very difficult to detect. This problem can be mitigated to some extent by choosing light element targets such as H [1], He [2][3][4], or diamond [5] and by pushing for lower thresholds. Alternatively, one may leverage inelastic processes such as the Migdal effect [6][7][8] or bremsstrahlung [9]. Inelastic processes occur at substantially lower rate, but are not subject to the constraint in (1) and can also yield signals that are more easily detected than a nuclear recoil, such as electronic excitations, ionizations or X-rays. Which approach is preferable depends on the characteristics of the detector. At sufficiently low energy and momentum scales, DMnucleus scattering is also not subject to (1) because atomatom interactions become important. In particular, the relevant excitations in a crystal target are phonons instead of elastic nuclear recoils. For m χ ≲ MeV, the momentum transfer from DM scattering corresponds to wavelengths comparable or larger than the interatomic spacing of a typical target. In this regime, the dominant process will be coherent scattering off multiple atoms, with creation of a single phonon. For crystalline targets with phonon energies as high as ∼100 meV, the energy deposited from DM can be well above the naive estimate in (1). Single phonon excitation has been studied extensively for sub-MeV dark matter, where numerical and analytic calculations by different groups are in good agreement [10][11][12][13][14][15][16][17]. These calculations have also been extended to diphonon production from sub-MeV dark matter 1 [25] as well as to single phonon production from MeV-GeV dark matter by including Umklapp processes [12,14]. However, so far there has not been a complete description of DM scattering for intermediate energy and momentum transfers, where multiphonon processes are expected to dominate.
In this work, we develop an analytic treatment of DM scattering that interpolates between the single phonon and nuclear recoil regimes. The relevant approximations are set primarily by the momentum transfer q. For single phonon excitations and q < 2π/a, where a is typical atomic lattice spacing, we use a long-wavelength approximation used earlier in the literature [10,11,13,25]. For q > 2π/a, we employ the incoherent approximation, which neglects interference effects between the response of neighboring atoms. This allows us to organize the calculation as a systematic expansion in the number of final state phonons, where each additional phonon comes with a factor of q/ √ 2m dωd . Here, m d andω d are the mass and average oscillation frequency of the atom in the position indexed by d. For q < √ 2m dωd it is numerically practical to compute the rate order-by-order in terms of the phonon density of states of the material. For q ≫ √ 2m dωd , scattering into many phonons dominates and the perturbation series requires increasingly large orders in q/ √ 2m dωd to converge. It can however be resummed by making use of the impulse approximation, which in turn smoothly matches onto the free nuclear recoil regime. A similar expansion in number of modes has been performed previously for the integrable toy model that is the harmonic oscillator [26]. Here we have generalized the approach to a harmonic crystal, analogous to the procedure followed in [27] and [28], in calculations of the Migdal effect and X-ray backgrounds, respectively. Fig. 1 illustrates our results from applying these approximations. All of our calculations are implemented as part 1 Analogous calculations were performed for superfluid He [18][19][20][21][22][23][24], for which diphonon production is the leading observable process for mχ ≲ 1 MeV. of the DarkELF public code [29]. 2 The remainder of this paper is organized as follows: In Sec. II we introduce the dynamic structure factor, which captures the material-dependence of the DM scattering cross section, and motivate the incoherent approximation for the structure factor. In Sec. III, we describe our analytic approximations in detail across the different regimes in energy and momentum transfer. We perform checks on our use of the incoherent approximation by comparing with previous calculations for single-phonon production and analytic calculations for diphonon production. Our results for GaAs are discussed in detail in Sec. IV and we conclude in Sec. V. Appendix A contains the formulas for diphonon production and Appendix B provides details on the impulse approximation. The implementation in DarkELF is documented in Appendix C. We further provide numerical results for Ge, Si and diamond in Appendix D.

II. DYNAMIC STRUCTURE FACTOR
Our starting point will be a general potential for spinindependent DM-nucleus interactions, although the formalism below could also be applied to spin-dependent interactions. For a DM particle of mass m χ incident on a crystal with N unit cells and n ions per unit cell, the potential in Fourier space is given bỹ Here, we sum over the N unit cells, labeled by lattice vectors ℓ, and atoms within the unit cell, labeled with the index d, such that all atoms in the crystal with positions r ℓd are summed over. The DM-proton scattering length b p is defined by the DM-proton scattering cross section σ p ≡ 4πb 2 p at some reference momentum, and µ χ is the DM-proton reduced mass. We first consider a general coupling strength f ℓd of the nucleus labeled by ℓ, d relative to that of a single proton. f ℓd is specified for various interactions in Section IV, such as nucleon number for scalar mediators and the effective electric charge for scattering via a dark photon mediator. In the latter case f ℓd is q dependent when accounting for screening effects.
We consider two form factors in (2) representing limiting cases of interactions: scattering via a heavy mediator, whereF (q) = 1; and scattering via a massless mediator, whereF (q) = q 2 0 /q 2 with a model-dependent reference momentum q 0 .
Collecting the overall factor 2πb pF (q)/µ χ , we define the differential cross section as where v is the initial velocity of the dark matter (incident on a target at rest), Ω c = V /N is the volume of the unit cell in the crystal, and ω q = q · v − q 2 /2m χ is the kinematic constraint on the momentum and energy transfers to the crystal q and ω. We have in turn also defined the dynamic structure factor (4) Note that the convention for S(q, ω) varies across the literature; here we use the convention that gives a similar S(q, ω) definition for both phonon interactions and DMelectron interactions [12,30]. We also assume the system is initially in its ground state |0⟩ prior to the collision, corresponding to a zero temperature system. We sum over final states with energies E f , such that each term represents the probability to excite the final state |Φ f ⟩.

A. Coherent and incoherent structure factors
For a given crystal there are many possible configurations of interaction strengths f ℓd which may vary even for different samples of the same material, e.g. the exact distribution of spins or isotopes in the material for spin-dependent 3 or mass-dependent interactions, respectively. This can be accounted for by averaging over a large collection of target samples. With a large number of nuclei in the crystal, we expect the exact distribution of interaction strengths in a given sample to be inconsequential relative to the result averaged over many samples. We can keep track of fluctuations away from the average configuration by splitting the scattering rate into a coherent and incoherent contribution, as explained below.
We follow the procedure of Refs. [31,32] and first reexpress (4) by expanding the square and Fourier transforming the δ-function, giving where C ℓ ′ d ′ ℓd is the time-dependent two-point function: In the second line we used the completeness of the basis of states. It will also be advantageous to define a shorthand notation for the auto-correlation function for an atom with itself as We assume that the f ℓd are random throughout the crystal. Under this assumption, the average of f ℓd f * ℓ ′ d ′ over target configurations, f d f * d ′ , must be independent of the lattice sites ℓ, ℓ ′ . Making this replacement in (5) gives where the averages may be written as For the d ̸ = d ′ case we assumed that the expectation values of the f d for different atoms in the unit cell are uncorrelated. This allows one to split the structure factor into two contributions: where the second line is obtained by adding and subtracting the term proportional to (f d ) 2 and regrouping. The first and second term in (11) are usually referred to as the coherent and incoherent structure factors in the neutron scattering literature. The coherent structure factor relays the scattering rate if the interaction strengths of all atoms in equivalent lattice sites are equal to a common value f d . For example, one can consider low energy, spin-independent neutron scattering in a very pure crystal with only a single isotope per atom type. This implies f d = f ℓd = A d , with A d the atomic mass number, such that the incoherent contribution in (11) vanishes exactly. The sum in (10) then crucially includes position correlators between differing nuclei, which capture the interference between different lattice sites. In practice, this interference leads to a coherence condition, which demands that the momentum in the scattering process must be conserved up to a reciprocal lattice vector. In particular, the 0th order term in a low q expansion of (6) corresponds to Bragg diffraction.
The incoherent structure factor on the other hand accounts for the statistical variations in interaction strengths between different scattering centers in the lattice. The second sum in (10) contains no cross terms and thus does not include interference between different lattice sites. There is therefore no corresponding coherence condition and the incoherent structure factor does not enforce momentum conservation. 4 For most earlier DM direct detection calculations the focus has been on spin-independent scattering in high purity crystals with little isotopic variation. In this scenario, we take the single isotope approximation f 2 d − (f d ) 2 = 0, implying that only the coherent scattering con-tributes. For spin-dependent dark matter scattering, the average will be the quantum expectation value of the spin operator, resulting in f 2 d ̸ = (f d ) 2 . We therefore expect the incoherent piece in (11) to be important in this case. In this paper we focus exclusively on spin-independent scattering in the single isotope limit and the corresponding coherent structure factors. The coherent structure factors are however more difficult to evaluate, due to the conservation of crystal momentum that is built into (6). This results in increasingly complicated phase space integrals for multiphonon processes [25]. For our purposes, the utility of studying the incoherent structure factor will be that the auto-correlation function can be used to obtain a reasonable and more manageable approximation of the coherent structure factor at sufficiently high momenta. Our results can also be extended to the case of spin-dependent scattering, but we leave this for future work.
Before venturing further into this approximation and its validity, we must first develop the structure factors into a form which lends itself to a direct calculation. In order to evaluate the structure factors in (4)-(8), the position vector of each atom may be decomposed in terms of the equilibrium lattice positions and displacement vectors, r ℓd = ℓ + r 0 d + u ℓd . Here r 0 d is the equilibrium location of atom d relative to the origin of the primitive cell and u ℓd is the displacement relative to that equilibrium. Following this decomposition, we quantize the displacement vector in the harmonic approximation with a phonon mode expansion The index ν denotes the phonon branches, of which there are 3n, and k labels the phonon momentum in the first Brillouin Zone (BZ). Theâ † ν,k andâ ν,k are the creation and annihilation operators for the phonons, ω ν,k is the energy of the phonon, e ν,d,k is the phonon eigenvector for atom d normalized within a unit cell, d e * ν,d,k · e µ,d,k ′ = δ µν δ k,k ′ , and m d is the mass of atom d.
The structure factor in (8) can then be explicitly evaluated by applying (12) to (6). For a pure single isotopic crystal with f 2 d = (f d ) 2 , this is given by [25] is the matrix element for scattering into the final state of the crystal denoted by f . The Debye-Waller factor e −W d (q) is given in terms of the function W d (q) ≡ 1 2 ⟨(q · u ℓd (0)) 2 ⟩. We may Taylor expand the inner exponential in powers of q where the nth term can excite a final state consisting of n phonons. The phonon eigenvectors and energies may be obtained numerically using Density Functional Theory (DFT) (see e.g. [33]); using these, the leading single phonon structure factor has been calculated [11,14,17]. These DFT-based calculations quickly become cumbersome, however, and have not yet been performed for generic n-phonon terms. Analytic calculations may be performed more easily, and have been carried out for the single-and two-phonon terms [25], but are only tractable when assuming an isotropic crystal and that |q| is small relative to the size of the first Brillouin zone. Such analytic calculations likewise lack scalability for higher order phonon terms. In summary, since the direct evaluation of (13) is very tedious and not always possible, we will rely instead on an approximate form of S (coh) (q, ω), bypassing the need to deal with (13). This is described in the next section.

B. Incoherent approximation
The incoherent approximation amounts to dropping the cross terms in (ℓ ̸ = ℓ ′ or d ̸ = d ′ ) from the sum in (10), thus neglecting the interference between nonidentical atoms. In other words, one approximates the coherent structure factor by The incoherent structure factor remains unchanged, and the total structure factor is then given by S (tot) (q, ω) ≈ N ℓ n d f 2 d C ℓd . In this work we will focus only on pure crystals with a single isotope for each type of atom, so that the total structure factor can be computed with (15). The incoherent approximation is expected to be a good approximation when the momentum transfer is larger than 2π/a with a the inter-particle spacing. Then the phase factors associated with the interference terms are expected to add up to a small correction compared to the ℓ = ℓ ′ , d = d ′ terms in the sum. For an argument justifying (15) we refer to [31,34].
For momentum transfers within the first Brillouin zone, single phonon scattering always dominates the inclusive scattering rate. It is however possible that the detector threshold is such that single phonon processes cannot be accessed but the double or multiphonon processes can. In this case the incoherent approximation cannot a priori be taken for granted. We nevertheless use it, but verify the results against our earlier two-phonon calculations [25] whenever possible (Sec. III B), finding satisfactory agreement. The accuracy of the calculations in this part of phase space is however less well understood and further work is needed.
To evaluate the auto-correlation function, we first replace the atomic positions r ℓd in (11) with their displacement operator decomposition, noting that the ℓ+r 0 d constant cancels, amounting to a simple substitution of r ℓd → u ℓd : The expectation value may be rewritten with an application of the Baker-Campbell-Hausdorff formula, Bloch's identity ⟨eÂ⟩ = e 1 2 ⟨Â 2 ⟩ , and some matrix algebra [11] giving: When we deployed Bloch's identity, we implicitly used the harmonic approximation, by only considering displacement operators of the form in (12).
The correlator ⟨q · u ℓd (0) q · u ℓd (t)⟩ may be evaluated with the form of the displacement operator in (12), wherein the ℓ dependence cancels. This gives which can be simplified further by averaging over the direction of momentum vector q where we defined the partial density of states for each atom in the primitive cell as The partial density of states was normalized to satisfy +∞ −∞ dωD d (ω) = 1. This can be shown by using the eigenvector completeness condition, which imposes ν e * ν,k,d,i e ν,k,d,j = δ ij for fixed k, d, where i, j are spatial indices. In addition, the total density of states of the material is defined by which satisfies +∞ −∞ dωD(ω) = n with n the number of atoms in the unit cell. 5 In materials such as Ge, Si, or GaAs all atoms in the primitive cell have the same or similar mass and as such contribute roughly equally to the density of states, see Fig. 2. One could therefore approximate D d (ω) ≈ D(ω)/n in (20) for these materials. We however choose to keep track of the partial density of states, to keep the calculations as general as possible.
For mono-atomic lattices, the density of states can be extracted directly from neutron scattering data through the incoherent structure factor. This is not always possible for multi-atomic lattices, since the scattering is only sensitive to the combination d |f d | 2 D d (ω)/m d . To infer the individual D d (ω) as well as D(ω), one therefore needs a set of scattering techniques which allows one to effectively vary the f d . This is not available for all materials, and it is therefore often most convenient to extract the D d (ω) from DFT calculations. A comprehensive library of results has been made available by the materials project [35]. 5 In the literature, the density of states is also sometimes normalized to 3na, where na is the atomic density.  [35]. Labels indicate the regions in which a particular phonon branch dominates.
Returning now to the calculation of the autocorrelation function, we can expand the exponential term in (17) using the form of the correlator in (20). This yields an explicit representation of C ℓd as an expansion in number of phonons n being excited: where the delta function arises from the time integral 1 2π dt e i( ωi)t e −iωt and ensures energy conservation. Here, by using (20), the Debye-Waller function takes the form of Thus, in comparison to the difficulties discussed surrounding (13), inputting this form of the correlator into (15) gives an analytic approximation for all phonon terms in the appropriate regime of validity.
In this paper, we utilize the incoherent approximation to calculate the contributions from higher-order phonon terms to an arbitrary degree in a simple and fast manner. This allow us to make rate predictions for the entire relevant mass range, going from the low-mass (m χ ≳ keV) single phonon regime to the high-mass (m χ ≳ 50 MeV) nuclear recoil regime.

III. PROCESSES
Using the autocorrelation function, (23), we can estimate the scale at which a generic n-phonon term starts becoming a relevant contribution to scattering. To organize the multiphonon expansion, it is useful to define an average phonon energȳ Whileω d technically depends on the atom d, this just gives an O(1) dependence in the phonon scale. Since n! ∝ n n at large n, we see that the nth term of the series (23) will roughly begin giving an O(1) contribution when This means that for a given q (or consequently, m χ ) one can determine the dominant scattering processes. When q 2 2m dωd ≲ 1, single phonon excitations will be the primary channel; for m d ∼ 30 GeV andω d ∼ 30 meV, this corresponds to q ≲ 30 keV. Conversely, when q 2 2m dωd ≫ 1, phonons are no longer a suitable description and the scattering is instead well modeled by the recoil of a single nucleus. This transition occurs roughly at q ≳ 2 √ 2m dωd . In between these two extremes, we have n ∼ few, indicating multiphonon excitations as the primary process. The precise nature of the dominant process for a given m χ will vary based on the mediator mass and experimental threshold.
In this section, we describe analytic approaches for characterizing the structure factor in crystal targets, broken into subsections corresponding to the previously mentioned processes. Secs. III A and III B deal with single phonon and two phonon excitations. Here we can also compare calculations of the full structure factor with the incoherent approximation. Sec. III C deals with many phonon excitations, and Sec. III D describes the impulse approximation, which gives a good approximation to the structure factor for momenta approaching the nuclear recoil limit. For all numerical results in this section, we will assume a coupling to nucleons (replacing the generic average interaction strength f d with the nucleon number A d ) for both massive and massless mediators, and take a GaAs target as a typical example of a simple cubic crystal of interest.

A. Single phonon production
If the unit cell contains at least two atoms, there are two types of phonons that can be produced: acoustic and optical phonons. As discussed in Sec. II, DFT-based calculations for both single acoustic and single optical phonon excitations have been performed across a large dark matter mass range (∼keV to GeV) [11,14,17]. Meanwhile analytic calculations so far have been limited q ≲ 1 keV, which corresponds to m χ ≲ MeV [10,25]. Although the DFT-based calculations span the entire mass range of interest and can provide information such as directional dependence, the numerics are more intensive; the phonon band structure, eigenvectors and structure factors must be calculated from first principles for each material. For high q, the sum over the reciprocal lattice must also be accounted for [12,16]. Here we extend the analytic calculations to the high q regime by using the incoherent approximation. The comparison with the DFT results of [11] will also serve as a validation of the incoherent approximation.
To organize the calculations, it is useful to define a momentum scale (q BZ ) which approximately reflects the size of the first Brillouin zone. We take q BZ = 2π a ≈ 2 keV, where a is the lattice constant. We first review the single phonon response for q < q BZ . In this regime, we compute the structure factors in the isotropic approximation and in the limit q ≪ q BZ . For this purpose we assume linear dispersions ω = c s q for the longitudinal acoustic (LA) and transverse accoustic (TA) modes, with c s replaced by c LA and c TA for the longitudinal and transverse sound speeds, respectively. The optical modes are assumed to have flat (constant) dispersions for the longitudinal optical (LO) and transverse optical (TO) phonon energies ω LO and ω TO . The sound speeds and optical phonon energies are taken to be their long-wavelength values (q = 0). We will refer to this set of assumptions as the long-wavelength approximation.
The matrix element is given by the leading non-trivial term in the small q expansion of (14). The only relevant contributions for q ≪ q BZ are those of the single LA and LO phonons. We approximate the long-wavelength acoustic eigenvectors as note that this form is valid for generic crystal targets and not limited to GaAs. For the LO phonon, we use the following eigenvectors, which are only valid for diatomic lattices [25] e LO,k,1 ≈ √ where the first atom is taken to be at the origin of the primitive cell, and the second atom is taken to be at the coordinate r 0 2 = (a/4, a/4, a/4) for GaAs. The acoustic and optical transverse eigenvectors are orthogonal to these, but do not contribute to the scattering into a single phonon. With these approximations and taking f d = A d , the analytic expressions for the single phonon contributions to the structure factor are [25] with Ω c the volume of the primitive cell. Here we have introduced a cut-off of ω = ω LO to the longitudinal acoustic branch to avoid overestimating the scattering rate with the LA mode near the edge of the Brillouin zone. The q 4 scaling and appearance of the lattice constant a in the optical structure factor comes from averaging over angles with the eigenvectors, giving (q · r 0 2 ) 2 ≈ q 2 a 2 /16 [13]. For dark matter with a standard velocity dispersion v ∼ 10 −3 , the typical momentum transfer begins to fall outside the first Brillouin zone for m χ ≳ 1 MeV. Physically, this corresponds to the wavelength becoming smaller than the interatomic spacing, and the long-wave length formulas from (27) to (31) are no longer valid. We can however utilize the incoherent approximation in (15) and (23), which yields The forms of the structure factor are qualitatively quite different in the two q regimes. In the coherent regime q < q BZ , summing over the response of multiple atoms with constructive interference leads to a resonant response in (32). The impact of the interference is greatly reduced for q > q BZ , such that the incoherent approximation becomes a viable description. While the sharp transition in the structure factor is an artifact of our approximations, (32)-(33) can accurately describe the integrated structure factor above or below q BZ . Fig. 3 compares our combined analytic single phonon description with numerical DFT calculations. For the DFT result we follow [11], computing the dynamical matrix and phonon dispersions with respectively VASP [36] and phonopy [33] (see also [14]), and take the angular average of S(q, ω) over all q directions for comparison with the isotropic approximation. The top panels show the structure factors in (32) as a function of q, integrated over ω. The top left panel shows S(q, ω) integrated over ω ∈ [1 meV, 27 meV] to select the acoustic phonon branches only and the top right panel shows the integral over ω ∈ [27 meV, 40 meV] for optical phonon branches. The analytic approximations are in good agreement with the DFT result in their respective regimes of validity. For q < q BZ , integrating (32) leads to respectively ∼ q and ∼ q 4 scaling, while the incoherent approximation in (33) always scales as ∼ q 2 . As discussed above, the ω-dependence of the analytic structure factors is quite different in the two regimes, with the coherent structure factor giving a resonant response around the single-phonon dispersion while the incoherent approximation is continuous in ω. However, the integrated result matches the full DFT calculation of the coherent structure factor well, indicating that the analytic approach will be useful in calculating integrated quantities such as rates. Furthermore, the analytic approach provides physical insight into the change in the q-scaling of the structure factor in Fig. 3a.
The plots in Fig. 3b show single phonon integrated rates for both massive and massless scalar mediators. For the massless mediator, scattering into the acoustic phonon specifically favors small q due to the ∝ q −4 contribution of the mediator form factor. The analytic result of (30) therefore applies across the entire DM mass range, as the large q contributions are negligible. For all other cases the structure factor scales with a positive power of q so that large q contributions are the most important. We therefore see a change in slope of the σ p reach around m χ ∼ MeV, when q ≳ q BZ becomes kinematically accessible. These features are captured by the q > q BZ analytic description from the incoherent approximation, and again agree with the DFT results.
B. Two-phonon production (q < qBZ) We next turn to the use and accuracy of the incoherent approximation for two-phonon production, in particular for q < q BZ . Single phonon production always dominates in this regime if above threshold [25]. It is however expected that there will be a phase in the experimental program for which the energy threshold will still be too high to access single optical and accoustic phonons, such that the formally subleading double phonon production can be relevant.
While the incoherent approximation is expected to be the least accurate for q < q BZ , it is still useful to compare  it with existing analytical results for the structure factor. The analytic results are obtained in the long-wavelength approximation, as defined in Sec. III A. In this limit, the Wilson coefficients of the self-interaction operators for the acoustic modes can be extracted from the measured or calculated elasticity parameters. With these assumptions, one can explicitly evaluate (13) to second order in q/ √ m d ω [25].
In this work, we will extend the long-wavelength calculations to all final states (see Appendix A) and compare them with the incoherent approximation. For this purpose we extrapolate the results of Ref. [25] to higher q values and make a number of additional assumptions to model the self-interactions of the optical modes, thus giving the complete structure factor. For these reasons the calculations in this section should however be considered only a toy model of a GaAs-like crystal. We will show below that for this toy model and in the limit of small momentum transfer, the incoherent and long-wavelength approximations give qualitatively similar DM scattering rates.
From Ref. [25], the two-phonon structure factor can be written as in the long-wavelength limit. The first term is the structure factor in the harmonic limit (also referred to as the contact piece in [25]), where anharmonic corrections to the atomic potentials are neglected. It can be obtained by expanding (14) to second order, and evaluated analytically in the long-wavelength limit. The second term contains contributions to the structure factor from anharmonic interactions. In order to evaluate this, one needs to include a phonon self-interaction Hamiltonian in computing (14), as described in detail in [25]. The interactions of acoustic phonons are based on an effective three-phonon Hamiltonian valid in the long-wavelength limit, but to obtain a more complete picture we include a highly approximate three-phonon Hamiltonian for interactions involving optical phonons. These calculations are summarized in Appendix A.
To perform the most meaningful comparison between the incoherent and long-wavelength approximations, we assume the following Debye model for the partial density of states for a diatomic crystal which is derived from the long-wavelength approximation as described in Sec. III A. 6 The explicit structure factor from using this toy density of states in (23) is given in Appendix A, which for simplicity we evaluate with A 1 = A 2 for GaAs. The top panel of Fig. 4 compares the calculations of the two-phonon structure factor in the incoherent and long-wavelength approximations. For the incoherent approximation, we show the result with the toy density of states in (35) as well as with the true density of states 6 Here the maximum momentum of the modes is determined by requiring that the sum over all modes is equal to the total number of degrees of freedom. For GaAs and in the isotropic approximation, the exact momentum cutoff is about 2% different from q BZ = 2π/a. This error is negligible compared to the uncertainties on the other assumptions made in this section.
from Fig. 2. The dashed line shows the harmonic limit, meaning that S (anh) is neglected. This is the case that is most directly comparable to the incoherent approximation, which assumes the harmonic mode expansion in (12). For the dotted line, the leading phonon selfinteractions were included.
In the harmonic limit, all modes scale as ∼ q 4 except for optical-acoustic final state, which scales as ∼ q 6 . The incoherent approximation naturally misses these more subtle destructive interference effects, but still captures the correct q 4 scaling for most of the modes. We see in Fig. 4 that the incoherent approximation is within a factor of ∼ 5 of the long-wavelength approximation for all ω > ω LO , for both the toy model and true density of states. The difference at smaller ω is not experimentally relevant, as the single phonon rate will completely dominate in this region. There are also delta-function terms from the optical-optical branches which do not appear in the plot; their contributions to the overall scattering rate are comparable for the incoherent and long-wavelength approximations as well. See Appendix A for details. These terms dominate the scattering rate at higher energies, and overall we see in Fig. 4 that the incoherent approximation reproduces the structure factor in the harmonic limit to within a factor of few.
When anharmonic interactions are included, the difference becomes larger and the incoherent approximation may under-predict the rate by up to an order of magnitude in our estimate. However, as discussed above, the anharmonic Hamiltonian used is itself also only valid at the order of magnitude level, particularly for optical modes. We expect that our approach can model the rate in this regime at the order-of-magnitude level, but a proper DFT calculation is needed for it to be rigorously validated.
Finally, we show in the bottom panel of Fig. 4 a comparison of the cross sections corresponding to a rate of 3 events/kg year, with the different approximations for the two-phonon structure factor. We assume ω > 40 meV, since for lower thresholds the rate is dominated by single-phonon production [25]. We emphasize that here we are only illustrating that the incoherent approximation is within a factor of few of the full structure factor, as long as the same assumptions are made for the phonon dispersion relations. Therefore, we restrict our comparison to m χ < MeV such that we can restrict to q < q BZ . The incoherent approximation underestimates the rate by a factor of few in the harmonic limit, and up to an order of magnitude when anharmonic interactions are included. Using the true density of states slightly improves the agreement. Though this compari- Optical-optical channels give a δ-function and are not plotted. Bottom: Cross sections for producing two phonons at a rate of 3 events/kg-year using the same approximations as above. We restrict the mass range to mχ ≲ 1 MeV so that typical q values are below qBZ, where our long-wavelength approximations are valid. The energy threshold is taken to be 40 meV, above the single phonon energies.
son only applies to a limited q range, our result suggests that the incoherent approximation should give a reasonable, order-of-magnitude estimate for multiphonon production even at low q. We expect this uncertainty to decrease for larger q where the incoherent approximation is most justified, and in particular we will see that the incoherent approximation reproduces the expected rate in the free nuclear recoil limit, as discussed in the next sections.

C. Multiphonon production
In the previous section, where we dealt with q < q BZ , the incoherent approximation should be viewed as an order-of-magnitude estimate only. For q > q BZ , it is however on firm ground [31,34] and is used routinely to measure the density of states from neutron scattering data [31]. Moreover, in the q ≫ q BZ regime multiphonon processes become important. This follows from the form of the structure factor, obtained by inserting (23) into the incoherent approximation (15): From the discussion around (26), the typical number of phonons is n ∼ q 2 2m dωd . Withω d ≳ 30 meV and m d ≳ 30 GeV for most crystals, the self-consistency condition for the incoherent approximation (q ≳ q BZ ) is therefore always satisfied for n > 2 processes. The evolution of (36) for increasingly large q is shown in Fig. 5a.
We can obtain an approximate scaling for (36) by separating each term in the sum over n into q-dependent and ω-dependent parts. The ω-dependent part is given by the second line of the equation, which is only non-zero at ω ≲ n ω LO in order to satisfy the delta function. This part of the structure factor can be estimated to have at most the value of 1/(n!ω n+1 d ); this is illustrated in Fig. 11 of Appendix C, where we plot the numerical result. For q ≲ √ 2m dωd (left and center panels of Fig. 5a), the Debye-Waller factor can be neglected and the structure factor then scales as S(q, ω) ∝ n 1 n! q 2 2m dωd n . For q 2 /(2m dωd ) ≲ 1, the structure factor therefore scales as S(q, ω) ∼ q 2m , with m the lowest number of phonons that is kinematically allowed. This scaling will be useful in Sec. IV, where we use it to extract the approximate scaling behavior of the DM cross section curves. It no longer holds for q ≳ √ 2m dωd (right-hand panel of Fig. 5a), where many modes contribute equally. This regime however can be understood in the impulse approximation, which is the subject of the next section.

D. The impulse approximation (q ≫ qBZ)
For q ≫ q BZ the sum of the multiphonon terms asymptotes to an approximately Gaussian envelope, as can be seen most clearly from the rightmost panel in Fig. 5a. This asymptotic form can be derived directly with a steepest descent approximation, also known as the impulse approximation. It is valid whenever the interaction with the probe particle happens on a time scale short compared to that of the phonon modes.
To derive this, it is most insightful to take a step back from (36) and return to using (20) in (17). The autocorrelation function is then When q ≫ √ 2m dωd , the exponent involving the density of states integral will be highly oscillatory in t, and the integral may be approximated by expanding about t = 0 through a steepest descent method. (See Appendix B). Doing so gives where ∆ 2 d ≡ q 2ω d 2m d . This approximation is referred to as the impulse approximation since the saddle-point around t = 0 dominates the rate. The true peak is shifted slightly from the result (38), which can be corrected by including higher orders in the expansion [37]. Including these additional terms has negligible impact on later results.
From (38), we see that the structure factor in the impulse approximation is which is a sum of Gaussians peaked around q = √ 2m d ω, one for each atom in the unit cell. In Fig. 5a we see that (39) is a reasonable approximation for q ≈ √ 2m dωd and converges rapidly to the full result in (36) for q ≳ 2 √ 2m dωd . As expected, it does not capture the features in the structure factor for q ≲ √ 2m dωd . In our final results, we use (39) for q > 2 √ 2m dωd , as it is numerically much faster than (36). For crystals composed of multiple atoms, we define the boundary as max d 2 √ 2m dωd . At this scale, the average number of phonons is about four, and it is sufficient to truncate the sum at n = 10 for all smaller q.
As we consider larger DM masses which access larger q and ω, the Gaussian becomes more sharply peaked. This can be seen by comparing the width ∆ d to the peak value ω = q 2 /2m d . In the large-q limit, we have so the Gaussian becomes narrow for ω well above the typical phonon energy. Then the narrow width limit exactly reproduces the expected free nuclear recoil delta function response: We therefore recover the familiar free nuclear recoil response for each individual atom in the unit cell.
In Fig. 5b we show cross section curves with a GaAs target, for both massive and massless scalar mediators. We compare the reach obtained with the full structure factor (in the incoherent approximation), the impulse approximation, and the free nuclear recoil limit. For m χ ≲ 20−40 MeV, the full structure factor must be used to capture the rate, depending on the mediator mass and threshold. For m χ ≳ 20 − 40 MeV, the q values compatible with the impulse approximation start to dominate,

GaAs, Multiphonon Response
(a) The first ten phonon structure factors in the incoherent approximation for GaAs, plotted for various fixed q. At sufficiently large q > √ 2m dωd , the total structure factor converges to the impulse approximation (IA, dashed line). In the right panel, there is a slight difference between the peak of the true structure factor and the impulse approximation. This can be accounted for in the impulse approximation by including higher orders in the steepest descent expansion [37]. (b) Cross sections for 3 events/kg-yr in GaAs for a hadrophilic mediator. Rates are computed with the n ≤ 10 phonon terms in the incoherent approximation (solid lines), the impulse approximation (IA; dashed), and the analytic free nuclear recoil result (NR; dotted). We see that at sufficiently high masses-and hence momentum transfers-the impulse approximation sufficiently recovers the result of summing the phonon terms. Likewise, for yet larger momenta the impulse approximation merges onto the free nuclear recoil result, as discussed in Sec. III D.

FIG. 5. Multiphonon transition into the nuclear recoil regime.
and we see that it reproduces the full result very closely. At even higher masses, the free nuclear recoil response becomes an excellent approximation, as expected.
A particular feature to notice from Fig. 5b is that the free nuclear recoil rate agrees with the impulse approximation result even in regions of the q, ω phase space where the Gaussian is not narrow. For example, for the massive mediator and m χ = 50 MeV, the rate will be dominated by momentum transfers q ∼ 2m χ v ∼ 100 keV, corresponding most closely to the rightmost panel of Fig. 5a. From (40) this gives ∆ d /ω ≈ 0.5 which is not particularly small. The nuclear recoil approximation nevertheless works remarkably well. The reason is that phase space integral in (3) has a trivial ω dependence aside from the S(q, ω) factor, since the delta function in ω just determines the region of phase space that is integrated over. Therefore, as long as the energy threshold is small compared to the peak in ω, the phase space integral over (39) and (42)  The "1-ph long wavelength" regime is discussed in Sec. III A, the "n-ph incoherent approximation" regime in Sec. III B and III C and the "Impulse approximation" region in Sec. III D.
E. Summary   Fig. 6 schematically illustrates the various approximations for the structure factor discussed in this section. The boundaries reflect only our choice of approximation and not a sharp transition in the behavior of the structure factor. The dotted gray parabola represents the phase space boundary for a given m χ and v (see Sec. IV). This parabola extends upwards and rightwards as m χ is increased, such that multiple different regimes are sampled for high enough m χ .
For the single phonon excitations (n = 1) described in Sec. III A, we use the long-wavelength and incoherent approximations for q < q BZ and q > q BZ , respectively. This combination gives good agreement with a full DFT calculation of the scattering rate, at least for a cubic crystal such as GaAs.
For multiphonon excitations (n ≥ 2), we use the incoherent approximation for the structure factor for all q below max d [2 √ 2m dωd ]. This is motivated by Sec. III B, where we argued that the incoherent approximation can serve as an order-of-magnitude estimate even for q ≪ q BZ . Given the limitations of the long-wavelength approximation, a dedicated DFT calculation is needed in this regime. For multiphonon excitations, we sum terms in (36) until we achieve convergence, as explained in Sec. III C. Finally, for q ≥ max d [2 √ 2m dωd ] we make use of the impulse approximation, which ultimately transitions into the well-known free nuclear recoil regime. This was explained in Sec. III D. Fig. 7 shows our full calculation of the structure factor for GaAs, overlaid with the phase space boundaries for a few representative DM masses. In the low q, single phonon regime, the response is given by a set of δfunctions on the LO and LA phonon dispersions, represented by the orange curves. At intermediate and high q, the structure function is modeled by a continuous function, where the layered structure for ω ≲ 50 meV reflects the various single and multiphonon contributions. At higher q and ω the individual resonances cease to be visible and one transitions into the smooth S(q, ω) predicted by the impulse approximation. At very high ω the structure function converges towards its free nuclear recoil form, which is represented by the black dashed line.

IV. RESULTS
In this section we convert our newly-gained understanding of the structure factor into concrete predictions for the DM scattering rate in a crystal target. The event  Fig. 6 for each corresponding regime in the (q, ω) phase space. For the massive mediator, we see the dominance of the single acoustic phonon at low masses and low thresholds, and of the optical phonon for intermediate thresholds. Eventually, for sufficiently high masses the process becomes dominated by the free nuclear recoil response. For the massless mediator, the q −4 form factor favors small momenta, and the rate is dominated by the lowest accessible mode for a given threshold.
rate per unit of target mass is where the experimental energy threshold is implicit in the boundary of the ω integral. f (v) is the DM velocity distribution, which we take to be with v 0 = 220 km/s, the Earth's average velocity v e = 240 km/s, and v esc = 500 km/s the approximate local escape velocity of the Milky Way. The scattering rate can be further simplified in the isotropic limit; using (3), where ω th is the energy threshold of the experiment, and the other integration limits 7 are Note (47) defines the phase space boundary shown in Fig. 6 for a given m χ and v. Finally, ρ T is the mass density of the target material and we have recast the rate in terms of the DM-proton scattering cross section σ p ≡ 4πb 2 p .

A. Massive hadrophilic mediator
In the case of a massive mediator coupling to baryon number, we calculate the scattering rate by taking f d = A d andF (q) = 1. The cross sections corresponding to a rate of 3 events/kg-year exposure are shown in the left panel of Fig. 8, assuming a GaAs target and for different energy thresholds. The same figures for Si, Ge and diamond can be found in Appendix D.
We can understand the numerical results in Fig. 8 analytically using the scaling of the structure factor discussed 7 In numerical implementations of (45), as done in DarkELF, it is beneficial to change the order of integration by first integrating over v, then q, and finally over ω.
in Secs. III A-III D. First, from (45), the m χ dependence of the rate is contained in The structure factor only contains positive powers of q across the entire phase space, so for a massive mediator, the integral (48) will be dominated by the largest kinematically accessible momentum transfers. For m χ ≫ 30 MeV, the kinematically allowed phase space is extended to q and ω where the free nuclear recoil approximation can be used. The rate therefore approximately scales as R ∼ 1/m χ for m p ≳ m χ ≫ 30 MeV. For low enough thresholds, this scaling holds even as the dark matter mass comes within O(few) of 30 MeV, where the structure factor is relatively broad in ω. The reason is that the kinematically allowed phase space is wide enough in ω that the integral over the Gaussian in the impulse approximation gives within a factor of few of the integral over the delta function in (42), as discussed earlier in Sec. III D.
For dark matter masses of 1 to 30 MeV, the allowed phase space is restricted to values of q < √ 2m dω . Here the structure factor can be expanded in powers of q/ √ 2m dω and favors small ω. As noted in Sec. III C the structure factor scales as ∼ q 2m , with m the smallest number of phonons whose total energy is above the energy threshold. We see there is significant threshold dependence: the single phonon final state strongly dominates the rate if it is above the energy threshold, while for higher thresholds only multiphonons contribute. The rate integral now scales as where q was evaluated at its maximum q ∼ 2m χ v. The ω integral does not contribute to the m χ scaling of the rate, since the integrand is peaked in ω somewhere near the energy threshold ω th . This expression then gives the approximate scaling R ∝ m 2m−1 χ . Since m is dependent on the energy threshold, this explains why different thresholds in Fig. 8 result in a different scaling as a function of m χ .
At even lower dark matter masses (m χ < 1 MeV), the phase space is restricted to q values within the first Brillouin zone, which is dominated by single phonon production in the long wavelength regime. If the threshold is low enough to access a single phonon, the scaling further depends on whether the threshold captures an appreciable part of the LA branch. If so, the leading contribution comes from the acoustic mode (30), which gives approximately independent of m χ . This behavior is clearly reproduced in Fig. 8 for the 1 meV threshold, for which the acoustic branch is always accessible. If the threshold is too high to access the acoustic branch, but can detect the optical branch, the structure factor has an extra q 3 scaling and we find R ∝ m 3 χ . This case occurs for m χ ≲ 0.3 MeV on the 20 meV curve in Fig. 8. For m χ ≳ 0.3 MeV the DM can excite the acoustic branch, resulting in a sharp enhancement of the rate.

B. Massless hadrophilic mediator
If we instead have a massless mediator that couples to baryon number, then by convention, the mediator form factor is taken to be |F (q)| 2 = mχv0 q 4 with v 0 = 220 km/s. The cross section curves for this scenario are given in the right panel of Fig. 8 again for different thresholds. As in Sec. IV A, we can analytically explain the scaling of the different curves across the DM mass range. The main difference with the massive mediator case is that for a massless mediator, there is a 1/q 4 scaling in the form factor, which leads to a scattering rate that generally favors low q and ω. The main contribution to the rate will therefore be much more threshold dependent across all DM masses.
If the threshold is small enough to access single acoustic phonon excitations, then this will be the dominant contribution to the rate at all masses. Again from (45) and using the analytic acoustic structure factor, the rate for thresholds that are sensitive to a single acoustic phonon scales as The integrand is largest at the smallest q, so we estimate the q integral by evaluating the integrand at q ≈ ω th /c LA in (46). The integrand therefore has no m χ dependence and gives the scaling R ∝ m χ for the ω > 1 meV curve in Fig. 8. Note however that this scaling behavior is sensitive to our convention for the reference momentum iñ F (q). For example, in models with both electron and nucleon couplings one often chooses to normalize the form factor with the reference momentum q 0 = αm e , which would yield R ∝ m −3 χ . If the LA branch is not accessible but the LO branch is, the production of a single LO mode will generally dominate. This introduces a different m χ dependence, which can be seen in Fig. 8 by comparing the 1 meV and 20 meV curves in the region with m χ ≲ 30 MeV. If m χ < 1 MeV, using the expression in (31) gives Unlike for the acoustic phonon, the structure factor favors high q so that the largest contribution is near q ∼ 2m χ v, giving R ∝ m 3 χ . If m χ > 1 MeV, the rate integrand is dominated by momentum transfers q ∼ q BZ . This is because when q > q BZ and ω ≤ ω LO we are using the incoherent approximation for single phonon production, where the q integrand drops as q −1 . Thus, we estimate the rate by integrating up to q BZ only: and find that R ∝ m χ . This is the reason why the 20 meV curve in Fig. 8 changes slope around m χ ∼ 1 MeV. We next turn to the intermediate mass range (1 − 30 MeV) with ω th > ω LO , such that n ≥ 2 phonons. In Fig. 8 this corresponds to the curves with thresholds of 40 meV and above. As in Sec. IV A, we again notice that the leading contribution to the structure factor will be given by the smallest number of phonons, m, that can exceed the threshold energy. In this regime, the integrand ∝ S(q, ω)/q 3 scales with positive powers of q for m ≥ 2 phonons, since (23) grows faster than q 3 . The analysis for multiphonons then follows exactly the same logic as the discussion in the previous section and we find that For large dark matter masses (≫ 30 MeV), again if the threshold is well above the single phonon energy, we can apply the free nuclear recoil approximation to obtain the scaling. Using the free nuclear structure factor gives The q-integral is dominated by low-momentum transfers along the free nuclear recoil dispersion, so we evaluate the integral at the intersection of ω = ω th and ω = q 2 2m d , or q = √ 2m d ω th . Then, the approximate scaling in this regime is R ∝ m χ /ω th , which we verify numerically in Fig. 8.

C. Dark photon mediators
The defining feature of a dark photon mediator is that it couples to the electric charge of the SM particles. In the regime where phonons are the relevant degrees of freedom, the charge of the nucleus is (partially) screened by the electrons. This means that we need a notion of an effective charge, as seen by the DM, which is momentum dependent. For individual atoms, this effective charge interpolates between zero in the low momentum, fully screened regime and the nuclear charge in the high momentum regime. We use the calculations from Brown et. al. [38] of the effective charge for individual atoms, as shown in Fig. 9. We expect this approximation to hold only for q ≳ q BZ , since additional many-body effects should be relevant for q < q BZ . This is particularly true for a polar material such as GaAs, where the Born effective charge of the Ga and As atoms is non-zero in the q → 0 limit. In this regime a full DFT calculation of the momentum dependence of the effective charge is needed, which we do not attempt here.
In this work, we will therefore focus on the momentum regime q ≳ q BZ , which corresponds to m χ ≳ MeV. In this case we can use the incoherent approximation and take f d = Z d (q), with Z d (q) the atomic effective charges in Fig. 9. This allows us to compute scattering rates with dark photon mediators for the production of two or more phonons, which is dominated by the highest kinematically accessible momentum transfers.
The regime q < q BZ is relevant primarily for massless dark photon mediators. (For massive dark photon mediators, there are strong BBN constraints that severely limit the scattering rate for sub-MeV dark matter, see e.g. [40].) In this regime, there are substantial devia- tions from the atomic effective charges due to the delocalized nature of the valence electrons. For instance, a polar material such as GaAs, SiC and sapphire can have a residual dipole moment associated with atomic displacements even for q → 0. The effective couplingsf d in this limit are given by Z * d /ϵ ∞ , where Z * d is the Born effective charge and ϵ ∞ is a screening due to valence electrons; the Born effective charges can be calculated with DFT methods [11,14,16]. This was treated in previous studies of single-phonon production through a massless dark photon mediator [10][11][12][13][14][15][16][17]. For non-polar materials such as Si, Ge and diamond, the Born effective charges vanish and instead multiphonon production is expected to dominate. This can be estimated with the energy loss function [29], at least for sub-MeV dark matter. Since this q < q BZ regime is already included in DarkELF [29], we restrict our results here to multiphonon processes with q > q BZ and ω > ω LO .
Our results are shown in Fig. 10 for GaAs; the results for Ge, Si and diamond are deferred to Appendix D. As is conventional for dark photon mediators, we choose the reference momentum for the massless mediator to be q 0 = αm e and present the results in terms of the effective DMelectron cross sectionσ e [41], with and µ χe the DM-electron reduced mass. In our calculations using the atomic effective charges, we impose q > q BZ to ensure we are not sampling the area of phase space for which these charges are clearly invalid. This means that our rate calculations for m χ ≲ 10 MeV are a slight underestimate of the true result.

V. CONCLUSIONS AND OUTLOOK
It is well-known that DM scattering in crystals can lead to one or more phonons being produced if DM has MeV-scale mass, as well as a recoiling nucleus if DM has GeV or higher mass. These processes are two sides of the same coin, depending on whether the momentum transfer is comparable to the inverse of the interparticle spacing and whether the energy deposition is comparable to the typical phonon energy ∼ω. When both momentum and energy scales are small, single phonon production dominates, and when both are large, nuclear recoils dominate. Here we studied the intermediate regime which is dominated by many phonons, which allows us to smoothly interpolate between single phonon production and nuclear recoils (see Fig. 8).
To make the multiphonon calculation tractable, we relied on the isotropic, incoherent, and harmonic crystal approximations. This allowed us to obtain analytic results for the scattering rate in terms of the phonon density of states in the crystal. These approximations are expected to be very good for q ≫ q BZ (m χ ≫ 1 MeV), as they explicitly reproduce the nuclear recoil limit when q ≫ √ 2m Nω . For q ≲ q BZ (m χ ≲ 1 MeV) the experimental threshold determines which theoretical treatment is most appropriate: for single phonon production, one can obtain analytic formulas by instead using a long wavelength, isotropic approximation. These results are currently only valid for cubic crystals such GaAs, Si, Ge and diamond. For strongly anisotropic materials such as sapphire, one must find a way to generalize them further or rely on DFT calculations. For multiphonon production and q ≲ q BZ , the situation is more complicated: in this case it cannot be taken for granted that anharmonic corrections to the various multiphonon channels can be neglected. The anharmonic multiphonon contributions involving optical modes are particularly difficult to model analytically, and at the moment we perform a simple estimate in a toy model to justify extrapolating the incoherent and harmonic approximations to q ≲ q BZ . A dedicated DFT calculation is needed to improve their accuracy.
Our approach provides a smooth description of sub-GeV dark matter scattering down to keV masses for hadrophilic mediators. For dark photon mediators, a DFT calculation of the momentum-dependent couplings in the q ∼ q BZ regime is needed to complete the interpolation. For both mediators, we have provided results for multiple direct detection materials of interest, and also included our calculation as part of the DarkELF public code package. These will be essential to interpret direct detection results as experimental thresholds for calorimetric detectors reach the eV scale and lower. In Sec. III B we compared the long-wavelength and incoherent approximations for the two-phonon final states, for q within the first BZ. In this appendix we provide the analytic expressions for both approximations.

Long-wavelength approximation
Here we discuss how we extend the analytic calculations from [25] for the coherent two-phonon structure factor to additional combinations of final state phonon pairs. As in Sec. III B, we assume a hadrophilic mediator with f d = A d throughout this appendix. It was shown in [25] that the structure factor separates into harmonic and anharmonic contributions S(q, ω) = S (harm) (q, ω) + S (anh) (q, ω) (A1) which do not interfere at leading order in the long wavelength limit. The first term involves expanding (13) to second order; note that it was referred to as the contact term in [25]. The anharmonic term is computed using an anharmonic phonon interaction Hamiltonian to first order. The specific matrix elements to be used are given in equations (12) and (13) of [25]. We take the longwavelength approximation for the phonon modes, as described in Sec. III A. For a crystal with two atoms in the unit cell, the longitudinal eigenvectors can be approximated by withk the unit vector along the phonon propagation direction. Note that the r 0 2 dependence was neglected in the LA eigenvector in (27) and in [25]; here we have kept this additional phase so that the acoustic and optical eigenvectors are explicitly orthogonal across a unit cell. This additional phase factor will only be relevant in cases where there is a destructive interference in the leading coupling to acoustic phonons, which occurs for some final states [13]. The transverse eigenvectors lay in the plane perpendicular tok and have analogous normalizations.
Analytic expressions for the harmonic structure factor were provided in Ref. [25] for acoustic-acoustic final states only. We require expressions for the opticaloptical and optical-acoustic final states as well to perform the comparison with the incoherent approximation. A straightforward application of (16) in [25] to the lowest order in q gives for the optical-optical modes.
For the optical-acoustic modes, the harmonic structure factors are of the form where x ≡ cLAq ω−ωLO . The other structure factors for optical-acoustic final states are given by relabelings LO → TO, LA → TA, where the expressions g expanded at small q are g (harm) TOLA (x ≪ 1) ≈ We see that at leading order in small q, the opticalacoustic structure factors are all suppressed by an additional factor of q 2 relative to the optical-optical modes, which is due to destructive interference. Since we will be comparing with the incoherent approximation at small q, we can effectively neglect these final states. We would also like to compute the anharmonic contributions to the 2-phonon structure factor, which we do with the inclusion of an anharmonic interaction Hamiltonian. For acoustic phonons in the long-wavelength limit, we have an effective Hamiltonian for acoustic phonons where the interactions are given in terms of macroscopic properties of the crystal through the Lamé parameters, as described in [25]. For the interactions of optical phonons, however, it is more difficult to write down a reliable analytic Hamiltonian. In this case we use (45) of Ref. [25], which comes from [42]. This Hamiltonian should be taken only at the order-of-magnitude level. We restrict the use of both effective Hamiltonians to the first BZ. The analytic expressions for the acoustic-acoustic and acousticoptical final states are given already, so we complete this by calculating the optical-optical terms. At leading order in q, this gives where c ≡ (c LA + c TA )/2. We have also assumed that the Grüneisen constant γ G ≈ 1.

Incoherent approximation
The second result needed for the comparison in Sec. III B is the two-phonon structure factor for GaAs in the incoherent approximation. To calculate this, we use the simplified density of states in (35) corresponding to the long-wavelength limit. Performing the n = 2 integral in (23) gives where each S is a contribution to the n = 2 structure factor from the part of the density of states associated with the subscripted modes, and the ellipsis indicates we sum over all combinations of modes. The first term of the sum in (A10) is and S TATA is given by S LALA with the replacement LA → TA and an additional overal factor of 4. The same procedure gives the LATA term as as well as the LOLA term, Again we may find S LOTA , S TOLA , and S TOTA by relabelings and inserting relevant factors of two for polarizations. Note that, since the incoherent approximation does not recover the q 6 scaling resulting from interference, we have written the structure factor here using q BZ = 2π/a to make the comparison more explicit. At lowest order in x and for A 1 ≈ A 2 , such a comparison of (A7) and (A13) shows a relative factor of 40/π 3 ≈ 1 for the LOLA channel. Lastly, for the remaining opticaloptical channels we find A comparison now of (A6) and (A14) shows the incoherent approximation gives a smaller structure factor by factors of 2π/5 -6π/5 ≈ 2 -4.

Appendix B: Impulse approximation
In this section we discuss how to obtain the impulse approximation form of the structure factor, (39) in Sec. III D. To achieve this we must approximate the t integral in (37) for large q. The expression in (37) can be written as with From this, we see there is a global maximum in the real part and a global minimum in the modulus of the imaginary part at t = 0. This allows us to perform a steepestdescent expansion about t = 0, giving where againω d = dω ′ ω ′ D d (ω ′ ). Note that the leading term in the expansion about t = 0 cancelled the Debye Waller factor, assuming the form given in (24). Evaluating the above gives which is the impulse approximation result.
In obtaining this form, we have assumed that any other local maxima in t gives a subdominant contribution to the t = 0 maximum. In particular, aside from the t = 0 point, which is a global maximum in Re[f (t)], there are local maxima in the real part which will generally be near integer multiples of 2π/ω d . The leading order contribution from each additional maxima t max is given by evaluating the real part in the exponential at the location of the maxima. This must necessarily be smaller than the t = 0 contribution since the following inequality is always satisfied Since t max ∼ 2π/ω d , the left hand side will be suppressed by an O(1) amount due to presence of the cos(ω ′ t max ). Then, the contribution from the local maxima will be exponentially suppressed: as long as the following condition is satisfied Here we have taken dω ′ D d (ω ′ ) ω ′ ∼ 1/ω d as a typical scale for this integral, although it will differ by an O(1) factor. Therefore, as long as the free nuclear recoil energy ω = q 2 /(2m d ) is well above the typical phonon energyω d for a scattering off of atom d, the t = 0 maximum is dominant and the impulse approximation should be accurate.
In the regime where q 2 /2m d is comparable toω d , the contributions from the additional maxima in t can become important. Nevertheless, the impulse approximation is still accurate at large ω even in this case because of cancellations from the rapidly changing phase in Im[f (t)]. When ω ≫ω d , then Im[f (t)] ≈ −ωt for t around t max ∼ 2π/ω d . This implies large oscillations of f (t) around t max , which suppresses the contribution from these local maxima. On the other hand, if ω ≲ω d , there may be large corrections to the impulse approximation due to these additional maxima.
These effects were shown in Fig. 5a when comparing the multiphonon expansion result to the impulse approximation. The middle panel showed the result if q = √ 2m dωd , in the m Ga ≈ m As approximation. For ω ≳ω d the structure factor falls smoothly and can be reasonably captured by the impulse approximation, while for ω ≲ω d ≈ 22 meV or at the optical phonon energies 31 and 33 meV there are sharp peaks in the multiphonon response that are not captured by the impulse approximation. For q = 2 √ 2m dωd the many multiphonon peaks merge and add up to a shape similar to the impulse approximation over the whole ω range. Practically, for our calculations, we use the impulse approximation for the structure factor at q > 2 √ 2m dωd . Though the approximation has small differences with the exact result when q ∼ 2 √ 2m dωd , integrating over the allowed phase space for the rate largely washes out these differences.

Appendix C: Implementation in DarkELF
In the main text, we presented the formulas in the manner which is most clear from the point of view of the various approximations and their regimes of validity. These formulas were not always suitable however for an efficient numerical implementation, which we address in this section. We also provide details on their implementation in the DarkELF package [29].
In the main text we gave the rate in the isotropic limit, (45). In order to calculate the rate for any mediator and to obtain the differential rate dR/dω, it is convenient to perform the v-integral first and rewrite the rate as: where now the integration limits are given by with v max = v esc + v e the maximum DM speed in the lab frame. The η function is given by with v min (q, ω) = q 2mχ + ω q . To evaluate the rate using incoherent approximation, we provide look-up tables for the structure factor. At each n for the sum in (36), the q and ω parts of the integral are separable, so we can capture the ω-dependent part with the family of functions and calculate the rate in terms of functions F n,d . These functions are simple to calculate numerically up to n ≤ 10, which we have tabulated and provided in DarkELF as look-up tables to speed up the calculation. The combinationω n F n (ω) is shown in Fig. 11 for GaAs in the m Ga ≈ m As approximation. For increasingly high n, the F n,d become increasingly smooth.
We have added several additional functions to DarkELF for the differential and integrated rate calculations from the single phonon to the nuclear recoil regime. Tab. I describes some of the new relevant functions. These functions currently work for materials with up to two atoms per unit cell. We have included the necessary data files for the multiphonon calculation for GaN, Al, ZnS, GaAs, Si, and Ge from a combination of DFT and experimental sources. We also allow the user to input their own calculations or extractions of the (partial) density of states, as well as momentum-dependent dark matternucleon couplings. Before calculating multiphonon scattering rates in DarkELF, it is necessary to tabulate the auxiliary function (C5) for each atom. This is done using the DarkELF function create Fn omega. This step is the most time consuming part of the calculation, so we provide these pre-tabulated for the aforementioned materials. For calculations with a user-supplied (partial) density of states, these tables must first be updated by running create Fn omega. DarkELF will save these new look-up tables for future computations, such that this step only need to be performed once. Next we describe DM-multiphonon scattering function description available for dRdomega_multiphonons_no_single(omega) Differential rate dR/dω in 1/kg/yr/eV all except SiO2, Al2O3 excluding long-wavelength single phonons

sigma_multiphonons(omega)
Nucleon cross section to produce 3 events/kg-yr all except SiO2, Al2O3 the functions that return important results. All of the following straightforwardly apply equations (C1-C4). R single phonon: This function takes the energy threshold and DM-nucleon cross sections and outputs the rate in the long-wavelength single phonon regime using the analytic functions (30)(31).
R multiphonons no single: This function takes the energy threshold and DM-nucleon cross section as inputs and calculates the total integrated rate, excluding the single phonon processes at long wavelengths q < q BZ . In other words, this calculation includes only the purple (multiphonon expansion) and red (impulse approximation) phase space regions in Fig. 6.
sigma multiphonons: This takes the energy threshold as input and returns the necessary DM-nucleon cross section to produce three events per kg-year for any number of phonons.
In order to return this cross section, this function first calculates the total rate by summing the outputs of R single phonon and R multiphonons no single, so it includes the entire calculation scheme depicted in Fig. 6.
dR domega multiphonons no single: This function takes the energy transfer ω and DM-nucleon cross section and returns the differential rate dR dω at that energy excluding single phonons in the long wavelength regime. This comes from equation (C1) without evaluating the ω integral. We exclude the single coherent phonon here since the long-wavelength approximation has delta functions in energy in the differential rate.

Appendix D: Additional results
Here, we provide additional results for Ge, Si, and diamond. Concretely, Fig. 12 shows the density of states for these three materials, as extracted from [35]. Fig. 13 shows the differential scattering rate via a massive scalar mediator for two example DM masses in GaAs, Ge and Si targets. Finally, Figs. 14, 15, and 16 are the cross section Densities of states for germanium, silicon, and diamond [35].
plots corresponding to an integrated rate of 3 events/kgyear for Ge, Si, and diamond, respectively. The electron recoil cross sections shown (dashed black lines) are based on calculations in [39] for Ge, Si and in [43] for diamond.