Anharmonic effects in nuclear recoils from sub-GeV dark matter

Direct detection experiments are looking for nuclear recoils from scattering of sub-GeV dark matter (DM) in crystals, and have thresholds as low as ~ 10 eV or DM masses of ~ 100 MeV. Future experiments are aiming for even lower thresholds. At such low energies, the free nuclear recoil prescription breaks down, and the relevant final states are phonons in the crystal. Scattering rates into single as well as multiple phonons have already been computed for a harmonic crystal. However, crystals typically exhibit some anharmonicity, which can significantly impact scattering rates in certain kinematic regimes. In this work, we estimate the impact of anharmonic effects on scattering rates for DM in the mass range ~ 1-10 MeV, where the details of multiphonon production are most important. Using a simple model of a nucleus in a bound potential, we find that anharmonicity can modify the scattering rates by up to two orders of magnitude for DM masses of O(MeV). However, such effects are primarily present at high energies where the rates are suppressed, and thus only relevant for very large DM cross sections. We show that anharmonic effects are negligible for masses larger than ~ 10 MeV.

Over the past few decades, a significant theoretical and experimental effort has been dedicated to detect dark matter (DM), but the particle nature of DM still remains a mystery.Direct detection experiments look for the direct signatures left by halo DM depositing energy inside the detectors.Traditionally, such experiments have looked for elastic nuclear recoils induced by DM particles in detectors [1].This strategy has had tremendous sensitivity for DM particles with masses higher than the GeVscale that interact with nuclei [2][3][4].However, in recent years it has also been recognized that sub-GeV dark matter models are also compelling and motivated dark matter candidates [5][6][7][8][9][10][11].These DM particles would leave much lower energy nuclear recoils, motivating experimental efforts to lower the detector thresholds for nuclear recoils.Inelastic processes like the Migdal effect [12][13][14][15][16] or bremsstrahlung [17] provide alternative channels to detect nuclear scattering in the sub-GeV DM regime.
The majority of experiments achieving lower thresholds in nuclear recoils (down to ∼ 10 eV) are doing so with crystal targets [18][19][20][21], although there is also progress in using liquid helium [22].Future experiments like SPICE [23] will reach even lower thresholds by measuring athermal phonons produced in crystals like GaAs and Sapphire (i.e.Al 2 O 3 ).In crystal targets, DMnucleus scattering can deviate substantially from the picture of a free nucleus undergoing elastic recoils.Nuclei (or atoms) are subject to forces from the rest of the lat-arXiv:2309.10839v1[hep-ph] 19 Sep 2023 tice, which play a role at the lower energies relevant for sub-GeV DM.For recoil energies below the typical binding energy of the atom to the lattice (O(10 eV)), the atoms are instead treated as being bound in a potential well.At even lower energies, the relevant degrees of freedom are the collective excitations of the lattice, known as phonons.In this regime, single phonon excitations with typical energies ≲ 0.1 eV are possible.
In the DM scattering rate, crystal scattering effects are all encoded within a quantity known as the dynamic structure factor, S(q, ω).The differential cross section for a DM particle of velocity v and mass m χ to scatter with energy deposition ω and momentum transfer q can be written in terms of S(q, ω) as: Here b p is the scattering length of the DM with a proton, µ χ is the reduced DM-proton mass, Ω c ≡ V /N is the volume of the unit cell in the crystal with total volume V and N unit cells, and ω q = q • v − q 2 /2m χ is equal to the energy ω lost by the DM particle when it transfers momentum q to the lattice.The q-dependence of the DM-nucleus interaction is encapsulated in the DM form factor F (q). S(q, ω) can thus be viewed as a form factor for the crystal response.For a recent review, see Ref. [24].Understanding S(q, ω) in crystals is critical to direct detection of sub-GeV dark matter.Thus far, the limiting behavior of S(q, ω) is well understood [25].In the limit of large ω and q (ω ≳ eV and q ∼ √ 2m N ω for nucleus of mass m N ), the structure factor behaves as S(q, ω) ∝ δ q 2 /(2m N ) − ω , reproducing the cross section for free elastic recoils.At low ω comparable to the typical phonon energy ω 0 and q comparable to the inverse lattice spacing, S(q, ω) instead is dominated by single phonon production.The intermediate regime, particularly q ∼ √ 2m N ω 0 , is dominated by multiphonon production.For a large number of phonons being produced, this should merge into the free nuclear recoil limit.
For DM masses below ∼MeV, the momentum-transfers are smaller than the typical inverse lattice spacing of crystals, q < 2π/a ∼ O(keV), where a is the lattice spacing.The dominant process is the production of a single phonon.In recent years, the single phonon contribution to S(q, ω) has been computed extensively in a variety of materials, often using first-principles approaches for the phonons [26][27][28][29][30][31][32][33][34][35][36].In most of the crystals, single phonons have a maximum energy of O(100 meV), however, requiring extremely low experimental thresholds to detect them.
Production of multiphonons is an enticing channel to look for sub-GeV DM with detectors having thresholds higher than O(100 meV).They are also important to understand in the near term as experiments lower their thresholds.However, multiphonon production has been more challenging to compute.The numerical firstprinciples approach taken for single phonon production does not scale well with number of phonons being produced, where even the two-phonon rate becomes very complicated.Alternate analytic methods are thus valuable.In Fig. 1, we show a classification of the different regimes in which a multiphonon calculation has been performed, including this work.We discuss the details of these regimes and calculations below.
One analytic approach was taken in Ref. [37], which calculated the two-phonon rate in the long-wavelength limit, but this study was limited to the regime q < 2π/a and focused on acoustic phonons only.For q > 2π/a, a different approximation is possible, the incoherent approximation, which drops interference terms between different atoms of the crystal in calculating S(q, ω).Then scattering is dominated by recoiling off of individual atoms.This approach was taken in [25], which found a general n-phonon production rate scaling as (q 2 /(2m N ω 0 )) n .This result also showed how the free nuclear recoil cross section was reproduced in the multiphonon structure factor as q ≫ √ 2m N ω 0 .However, one limitation of the multiphonon production rate in Ref. [25] was that it worked in the harmonic approximation, where higher order phonon interactions like the three-phonon interaction are neglected.Typical crystals have some anharmonicity which introduces phonon self-interactions, leading to various observable effects like phonon decays, thermal expansion, and thermal conductivity of crystals [38][39][40].Using a simplified model of anharmonic phonon interactions, Ref. [25] estimated that anharmonic three-phonon interactions may give the dominant contribution to the two-phonon rate q < 2π/a, and are larger than the harmonic piece by almost an order of magnitude in the regime.On the other hand, we do not expect anharmonic effects to be important in the opposite limit of large q (q ≫ √ 2m N ω 0 ), where the nucleus can be treated as a free particle.It is thus necessary to bridge these two extremes and estimate the anharmonic effects in the intermediate regime where multi-phonons dominate the scattering.
In this work, we estimate the anharmonic effects on the rate of multiphonon production by working in the incoherent approximation and q > 2π/a.In this limit, the multiphonon scattering rate looks similar to that of an atom in a potential [41], although the spectrum of states is smeared out due to interactions between neigh- Here we show a classification of regimes in which a multiphonon calculation has been performed, as well as approximations made in each case.In this work, we show that anharmonic corrections can be significant for q ≲ √ 2mN ω0 (Sec.III B) but are negligible when q ≫ √ 2mN ω0 (Sec.III C).We obtain results for all q using numerical calculations (Sec.IV A). (Right) To estimate anharmonic effects, we take a toy model of dark matter scattering off an atom in a 1D anharmonic potential.We obtain the anharmonicity by fitting to empirical models of interatomic potentials.
boring atoms.Given this similarity, we will take a toy model of an atom in a 1D potential.This gives a simple approach to including anharmonic effects, which is also illustrated in the right panel of Fig. 1.The anharmonic corrections to the atomic potential only capture a part of the contributions to anharmonic phonon interactions, but they have a similar size (in the appropriate dimensionless units) and should give a reasonable estimate of the size of the effect.We can therefore use this approach to estimate theoretical uncertainties and gain analytic understanding for the multiphonon production rate.However, the result should not be taken as a definitive calculation of the anharmonic corrections.Fortunately, we will find that anharmonic corrections are large only in certain parts of the phase space which are more challenging to observe, and that the multiphonon rate quickly converges to the harmonic result for DM masses above a few MeV.
The outline of this paper is as follows: In Sec.II, we discuss the formalism of DM scattering in a crystal and the dynamic structure factor, which encodes the information about the crystal response.We consider the calculation of the structure factor under the incoherent approximation, and motivate the anharmonic 1D toy potentials we use in this paper.In Sec.III, we study the behavior of the dynamic structure factor analytically for the anharmonic 1D potentials.Using perturbation theory, we show that anharmonic corrections can dominate for q ≪ √ 2m N ω 0 and become more important for higher phonon number.In the opposite limit q ≫ √ 2m N ω 0 , we use the impulse approximation to show that anharmonic corrections are negligible and that the structure factor indeed approaches that of an elastic recoil.In Sec.IV, we present numerical results for the structure factor in anharmonic 1D potentials obtained from realistic atomic potentials in various crystals.In Sec.IV A, we calculate the impacts of including anharmonic effects on DM scattering rates.We conclude in Sec.V.
Appendix A gives the details of the modeling of the interatomic forces on the lattice, used to extract 1D single atom potentials.Appendix B gives additional details of the analytic perturbation theory estimates of the anharmonic structure factor.Appendix C includes additional details relevant to the impulse approximation calculation.Appendix D summarizes the exactly solveable Morse potential model, which further validates the results in the main text.

II. DARK MATTER SCATTERING IN A CRYSTAL
Consider DM that interacts with nuclei in the crystal.We will parameterize the interaction with the lattice by a coupling strength f ℓd relative to that of a single proton, where ℓ denotes the lattice vector of a unit cell and d denotes the atoms in the unit cell.In the DM scattering cross section, (1), the material properties of the crystal are encoded in the structure factor S(q, ω) which is defined as, where |Φ f ⟩ is the final excited state of the crystal with energy E f and r ℓd denotes the position of the scattered nucleus.The crystal is considered to be in the ground state |0⟩ initially.Note for simplicity we assume a pure crystal where each atom has a unique coupling strength; the scattering is modified if there is a statistical distribution for the interaction strengths at each lattice site, for instance if different isotopes are present [25].
The states |Φ f ⟩ are the phonon eigenstates of the lattice Hamiltonian, where the first term is the kinetic energy of the atoms in the lattice and the lattice potential V lattice in general is given by, where the u α (ℓd) is the displacement from the equilibrium position in the Cartesian direction α for the atom at the position d in the unit cell located at ℓ, and k αβ , k αβγ are the second-, and third-order force constants respectively.Note that as the displacements are considered around equilibrium, we do not have a term in the potential which is linear in the displacements.
A number of approximations are useful in evaluating S(q, ω).The first is the harmonic approximation, which amounts to keeping the terms up to second-order force constants and neglecting the higher order terms (k (3) αβγ = 0).This vastly simplifies the Hamiltonian into a harmonic oscillator system, and has been used in most previous calculations of DM scattering in crystals.While this is generally an excellent approximation in crystals, including higher order terms in the Hamiltonian (anharmonicity) is necessary to explain a number of observable effects, as we will discuss further below.
The second approximation is the incoherent approximation, used for scattering with momentum transfers much bigger than the inverse lattice spacing of the crystal, q ≫ 2π/a.In this limit, we drop the interference terms between different atoms in the crystal in (2).This amounts to summing over the squared matrix elements of individual atoms in the structure factor in (2), The calculation of the structure factor then simplifies to computing matrix elements ⟨Φ f |e iq•r ℓd |0⟩ 2 which are identical for the atoms in all unit cells ℓ.
Below, we will first discuss this calculation under the approximation of a harmonic crystal, before going on to setting up a model that accounts for anharmonicity in crystals.

A. Harmonic approximation
In the harmonic approximation, the lattice Hamiltonian can be written as a sum of harmonic oscillators in Fourier space [42], where the phonon eigenmodes of the lattice are labelled by the momentum q and the 3n branches ν with n being the number of atoms in the unit cell.The â † q,ν (â q,ν ) are the creation (annihilation) operators, and ω q,ν are the energies of the phonons.The lattice eigenstates that appear in (2) can then be written as, where |Φ n ⟩ is an n-phonon state.The displacement operators in this harmonic approximation are given by, where the e q,ν (d) indicates the eigenvector of the displacement of atom d for that phonon.The equilibrium position of the atom is denoted by r 0 ℓd .Using r ℓd = r 0 ℓd +u(ℓd) inside (2), the dynamic structure factor can be calculated in the harmonic approximation.This approach has been applied to calculate single-phonon excitations using numerical results for phonon energies and eigenvectors [27][28][29][30][32][33][34], but becomes computationally much more burdensome for multi-phonons in the final state.
Under both the incoherent and harmonic approximations, it is possible to compute the multiphonon structure factor in (5).This was given in Ref. [25] as an expansion in the number of phonons produced n, S(q, ω) ≈ 2π where D d (ω) is the partial density of states in the crystal, normalized to dωD d (ω) = 1.W d (q) is the Debye-Waller factor defined as, (9) shows that with higher momentum q, there is an increased rate of multiphonons; the typical phonon number is n ∼ q 2 2mω with ω a typical phonon energy.In the limit of n ≫ 1, this reproduces the nuclear recoil limit.
In the incoherent approximation above, we still assumed the final states |Φ f ⟩ are the phonon eigenstates of the harmonic lattice Hamiltonian in (6).Let us now make a further approximation that the final states are isolated atomic states, where each atom is bound in a potential.Assuming an isotropic potential, and a single frequency ω 0 for the oscillators, a toy atomic Hamiltonian for atom d in the lattice can be written as, where r d is the displacement of the atom d from its equilibrium position.Following (5), the dynamic structure factor can be written as, where |⃗ n⟩ are the energy eigenstates of the toy harmonic Hamiltonian considered for atom d, with ⃗ n = {n x , n y , n z }.The energies with respect to the ground state equilibrium are given by E n − E 0 = nω 0 with n = n x + n y + n z .We have also absorbed the sum over the lattice vector ℓ and the volume V into the density n d of atom d in the lattice.As shown in [41], this structure factor is given by, where the Debye-Waller factor in the toy model is given by, W toy d (q) = q 2 /4m d ω 0 .This picture can be simplified even further by considering a toy one-dimensional harmonic potential for the atom d given by Note that in general ω 0 will depend on the atom d within the unit cell, but we suppress this dependence for simplicity.The structure factor in this 1D case is exactly the same expression as the toy three-dimensional case in (13), as expected given the isotropic 3D potential assumed.A derivation of the 1D result is given in Sec.III A.
The toy model of DM scattering off a 1D harmonic potential gives a simple intuitive picture for the result in (9).We see a very similar form of the structure factor in (13), but with a discrete spectrum of states for the isolated oscillator of the toy model.By assuming that the final states are isolated atomic states, we have effectively neglected the interactions between atoms, and the excited states of all the atoms are discrete and degenerate.In a real material, the interaction with neighboring atoms will lead to a splitting of the degenerate levels, and give a broad spectrum of allowed energy levels (the phonon spectrum).The interpretation for the structure factor is therefore also somewhat different in the two cases, as it gives a probability for exciting the nth excited state in an isolated oscillator.But we will still continue to refer the nth excited state as the n-phonon state to make the connection with the full incoherent structure factor in (9).
The similarity in the structure factor gives a route forward to including anharmonic effects, which is much easier to understand in the toy model.We can proceed by including anharmonic corrections to the 1D potential in (14), and in some cases obtain analytic results that illustrate their importance.In order to quantitatively estimate the impact on dark matter scattering rates, a few remaining ingredients are needed.In practice, the toy Scattering in Harmonic Crystal and 1D Oscillator FIG. 2. Comparison of scattering in a harmonic crystal to 1D harmonic oscillator.The dotted lines show the DM cross section reach computed using the multiphonon structure factor in a harmonic crystal, (9), and assuming the incoherent approximation [25].Using the structure factor of the toy 1D harmonic oscillator in (13) combined with the energy smearing prescription in (16) gives a very similar result (solid lines).There are some small deviations at low momentum since we place a hard cut on the allowed momentum transfer q > 2π/a ≈ 2 keV for the 1D oscillator.
model can give very different results in certain parts of parameter space due to the discrete spectrum assumed and depending on the choice of ω 0 .We therefore need a prescription to identify the appropriate ω 0 for the isolated oscillator, and to smear it out appropriately to mimic a real material.Comparing Eqs. 9 and 13, we see that the complete structure factor can be attained by making a replacement In this expression, we can identify D(ω)/(ωω −1 ) as a normalized probability distribution for ω, where ω −1 = dω ′ D(ω ′ )/ω ′ .This distribution yields a mean value for ω of (ω −1 ) −1 .The right hand side of ( 15) is proportional to the joint probability distribution for total energy ω, and we can simplify it when n ≫ 1 by applying the Central Limit Theorem.This allows us to replace the right hand side with a Gaussian, which simplifies computations: Θ(ω max − ω).
(16) Note we have included a cutoff at multiples of the max-imum allowed energy in the density of states, ω max = n × (min(ω)|D(ω) = 0) so that we do not include the region where D(ω i ) = 0 on the right hand side of (15).The width of the Gaussian for n = 1 is given by and ω = dω ′ D(ω ′ )ω ′ .This discussion therefore makes it clear that we should identify the frequency of the 1D toy model as ω 0 = 1/ω −1 , which can be calculated numerically given the phonon density of states.This approach is validated in Fig. 2, where we compare our previous result using the full density of states [25] to the prescription described above.Note that small deviations at low mass arise from the lack of a cutoff at the Brillouin zone momentum in the previous density of states result.We reiterate that in this work, we shall include this Brillouin zone cutoff across all rate calculations since the incoherent approximation and subsequent approximations are only valid in this regime.We will utilize this prescription to extend the multiphonon calculations for an anharmonic potential.To set up toy 1D anharmonic potentials, we first need to understand the anharmonic properties of typical crystals to extract the behavior of the potentials.We do this in the following subsection.

B. Anharmonic crystal properties
In general, a crystal lattice will exhibit some anharmonicity.Anharmonicity technically refers to the presence of non-zero force constants which are higher than second-order in the lattice potential in (4).For example, cubic anharmonicity in the crystal is parameterized by the third-order force constants k (4).Such force constants can be computed with DFT methods, similar to the harmonic case [43].In the presence of such terms, the phonon eigenstates are no longer the harmonic phonon eigenstates of the crystal, and higher order phonon interactions, such as a three-phonon interaction, will be present.Calculating the full dynamic structure factor in (5) for a crystal with such anharmonicity would require accounting for these higher order force tensors in both the matrix elements and in the final states, which quickly becomes a very challenging numerical problem.The rough size of the anharmonic force constants can be inferred from measurable crystal properties, however.We will briefly discuss some of the anharmonic effects below, and use them to justify our estimate of anharmonic effects.
An important effect of keeping cubic or higher order terms in ( 4) is to introduce interactions between the phonon modes which are the eigenstates of the harmonic Hamiltonian.For example, from (8), we can see that a cubic term in the displacements u(ℓd) will introduce three-phonon interactions like â † q,ν âq ′ ,ν ′ âq ′′ ,ν ′′ (i.e.annihilation of two phonons to create a single phonon) or â † q,ν â † q ′ ,ν ′ âq ′′ ,ν ′′ (i.e.decay of a single phonon into two phonons) in the Hamiltonian at the first order in the anharmonic force constant k (3) .Phonon lifetimes in crystals are thus directly related to the anharmonic force constants, and can be measured to estimate the size of the anharmonicity [40,44,45].
Anharmonicity is also necessary to explain thermal expansion and conductivity in crystals.In particular, the linear volume expansion coefficient of crystals can be directly written in terms of the mode Gruneisen constants γ qν which is defined for phonon modes labelled by the momentum q and branch index ν as [46], Note that the change in volume in the equation above is at a fixed temperature.In a purely harmonic crystal, the phonon frequencies are determined by the secondorder force constants which do not get modified with changes in volume, thus leading to zero Gruneisen constant.However, in the presence of cubic anharmonic-ity, the phonon frequencies are determined by the effective second-order force constants, which receive corrections depending on both the third-order force constants k (3) and the changes in volume, thus giving a non-zero Gruneisen constant [47].An increase in volume leads to larger displacements of atoms, which typically makes the effective second order constants and the phonon frequencies smaller, providing a positive Gruneisen constant.In the case of a non-zero Gruneisen constant, the free energy of the crystal, which has a harmonic contribution ∝ ∆V 2 , receives a volume-dependent correction ∝ −∆V γ qν Ēqν , where Ēqν is the mean energy in the phonon mode qν at a particular temperature [38].As the temperature increases, the mean energy Ēqν goes up, and thus this leads to a new equilibrium volume which minimizes the free energy.For a positive Gruneisen constant, this leads to thermal volume expansion.The Gruneisen constants are thus directly related to the cubic force constants of the material, and have also been used to extract them [47].Concretely, the relationship between the mode Gruneisen constants and the anharmonic force constants for weak anharmonicity can be shown to be [48], where the e β q,ν (d) indicates the displacement of atom d in the Cartesian direction β for the phonon qν, and r 0,α 0d is the equilibrium position of atom d in the Cartesian direction α for the unit cell at the origin.To get a rough estimate of the maximum anharmonicity strength in the crystal, the relation in (19) can be inverted and written in terms of the maximal mode Gruneisen constant γ max found in a crystal, where ω 0 is the typical phonon energy of the lattice and l is the nearest neighbor distance.Now consider a typical displacements ∼ ( √ 2m d ω 0 ) −1 of an atom in the crystal; the change in the potential energy δV anh due to anharmonic force constant estimated above is given by, where in the second line we use parameters for Si.We use an estimate for the maximal value of the mode Gruneisen constant in Si from [38] at 0K.In Ge, the maximal Gruneisen constant is similar to that in Si, while in GaAs, it could be as high as 3.5 for certain phonon modes [38].
The Gruneisen constant thus provides a rough estimate of the overall anharmonicity in the crystal, including the cubic terms which depend on displacements of multiple atoms.
In this paper, we will work with a toy model of anharmonic interactions similar to the 1D oscillator model in Sec.II A. In particular, we consider excitations for an isolated atom in a 1D anharmonic potential.The anharmonicity is controlled by force constant terms like k terize the modification to the potential of a single atom in a lattice.Since the Gruneisen constants involve a sum over many cubic force terms, we instead directly obtain the single-atom anharmonic force constants with an empirical model of the lattice.
We model the lattice assuming empirical interatomic potentials, which have been shown to accurately reproduce phonon dispersions and transport properties [49].Concretely, we assume the Tersoff-Buckingham-Coulomb interatomic potential with the parameter set given in Ref. [49] (see Appendix A for details).We then fix all atoms at their equilibrium positions except for one atom denoted by ℓd, which is displaced by a small distance in different directions.The single atom potential calculated from this procedure is shown in Fig. 3 for Si, with deviations from the harmonic potential that depend on the direction of displacement.The maximum anharmonicity is along the direction of the nearest neighbor atom.Along this direction, we find that the typical change in the potential energy for an atom displaced by r ∼ ( Comparing this estimate with ( 21), we see that the anharmonicity strength inferred from the potential of a single atom is roughly of the same size as the overall anharmonicity strength of the lattice inferred from the Gruneisen constant.Thus, even though we do not perform a full calculation of the structure factor for an anharmonic crystal including the modification of the phonon spectrum and the lattice states, the comparison above suggests that the effects in a full calculation are expected to be similar in magnitude to the effects we estimate in this work using single atom potentials.
Single atomic potential: Potential of a single atom displaced along various directions with all other atoms at their equilibrium positions.In zincblende Si, the largest anharmonicity is in the direction of the nearest-neighbor atom, while the smallest anharmonicity is in the direction of the next-nearest-neighbor.We have also included a third direction orthogonal to the other two, with intermediate anharmonicity strength.

C. Toy anharmonic potential
As shown in Sec.II A for the harmonic crystal, the features of the dynamic structure factor under the incoherent approximation can be well-approximated with just a 1D toy potential for an individual atom.This gives a much simpler path to calculating DM scattering in anharmonic crystals for q ≫ 2π/a, where many phonons may be produced.In contrast, prior work including anharmonicity focused on the limit q ≪ 2π/a, restricted to two phonons [37], and does not scale well to large number of phonons.We can then stitch together the two approaches to gain a more complete understanding of anharmonic effects.
In this work, we take a 1D anharmonic potential and calculate the 1D structure factor, in order to simplify the problem as much as possible.Taking the 1D approximation is more subtle in the presence of anharmonicity since a generic potential in 3D is not separable, unlike the harmonic case.Denoting the small displacement around equilibrium by r, and the polar and azimuthal directions by θ and ϕ respectively, the potential energy for atom d in the lattice can be expanded in powers of r as, where λ k are dimensionless constants parameterizing the degree of anharmonicity at k th order, and f k (θ, ϕ) are functions which specify the angular dependence and whose range is [−1, 1].Solving the full 3D problem would require numerically finding the eigenstates of this general potential, while in the 1D case we can make much more progress analytically.We will therefore select directions of maximum anharmonicity and use this for our simplified 1D problem.Our expectation is that this gives a conservative estimate of the importance of anharmonic couplings, in that the full 3D calculation would give somewhat reduced effects.
As discussed in Sec.II B, we can extract realistic single atom potentials by modeling the interatomic potentials on the lattice and displacing a single atom (see Appendix A for details).We typically find that, for small displacements around equilibrium, the anharmonicity is dominated by the cubic and quartic terms parametrized by λ 3 and λ 4 , respectively.Motivated by these observations, we consider the following forms of toy potentials in our study: • Single cubic or quartic perturbations: We first consider a harmonic potential with a single perturbation, where k = 3 or 4.This case is amenable to perturbation theory, and in Sec.III B, we apply it to discuss the power counting of anharmonic corrections.
• Morse potential: It is possible to obtain exact (non-perturbative) analytic results for the Morse potential defined by, where a is a parameter controlling the width of the potential and B is the normalization.We fit these two parameters to the cubic anharmonicity estimated from the single atom potentials discussed earlier, and calculate the dynamic structure factor for this potential in App.D.
• Fit to realistic atomic potentials: We numerically calculate the structure factor in a potential with both cubic and quartic terms, where the dimensionless anharmonic couplings are obtained by fitting to the actual single atom potential.The potential in this case is given by We find that typically, λ 3 ∼ 0.01, and λ 4 ∼ 10 −4 .
For the 1D toy potentials discussed above, we compute the 1D dynamic structure factor in the incoherent approximation (q ≫ 2π/a): Again, we have summed over all atoms of type d in the lattice and defined the number density of atom d by n d .The wavefunctions |Φ⟩ are the eigenfunctions of the Hamiltonian, The computation of the dynamic structure factor then boils down to computing the ground state |0⟩ and the excited eigenstates |Φ f ⟩ for this Hamiltonian, and calculating the structure factor under the incoherent approximation as in Eq. ( 27).
As discussed in Sec.II A, for a 1D toy model the phonon levels are discrete and in a real crystal there is a broad spectrum of energy levels.Similar to the harmonic case, we need a prescription to account for this smearing of energies.In the case with anharmonicity, the spectrum is shifted.The 1D toy model will instead give a modified energy-conserving delta function: where f (n)ω 0 is the energy difference between the nth excited state and the ground state.f (n) will depend on the exact form of the potential.Guided by the harmonic result, we again shall fix ω 0 = 1/ω −1 and introduce a width to the delta function in a similar fashion: This is in the 1D approximation, and that including the full 3D anharmonic potential would be expected to have an additional effect on the spectrum of states.However, in practice the anharmonicity is sufficiently small that the shift of the spectrum is subdominant to the other anharmonic effects in the structure factor.This forms the basis of the toy model we consider in this paper.Focusing on the high q regime where the incoherent approximation applies, we consider independent lattice sites and calculate scattering in them with 1D toy anharmonic potentials.We now describe different approaches to understand the dynamic structure factor in this setting.

III. ANALYTIC RESULTS FOR STRUCTURE FACTOR
In this section, we study the features of the structure factor for a 1D anharmonic potential with analytic methods.This will allow us illustrate the general behavior for the limits q ≪ √ 2m d ω 0 and q ≫ √ 2m d ω 0 .First, we review the derivation of the structure factor for a 1D harmonic potential.For n-phonon production in the harmonic limit, the structure factor in the regime Treating the anharmonic 1D potential as a perturbation, we then show that the q dependence of the n-phonon term can be substantially modified in the regime q ≪ √ 2m d ω 0 , leading to large anharmonic corrections.In particular, we obtain the power counting of the structure factor in powers of q and the anharmonicity parameter λ k , which allows us to roughly identify the regime of q where we expect the anharmonic effects to be dominant.As we will see later, this proves useful to explain the numerical results for realistic potentials.
Finally, we will also use the impulse approximation to perform an analytic estimate of the structure factor in the regime q > √ 2m d ω 0 .We show that the nuclear recoil limit is reproduced, with the structure factor approximated by a Gaussian envelope similar to the harmonic case.Anharmonic terms give rise to slightly modified shape of the Gaussian, which have negligible impact on scattering rates.

A. Harmonic oscillator
First, we briefly review the calculation of the dynamic structure factor in the harmonic approximation.In this case the potential V d (x) is given by The energy E n of the n-th excited state |n⟩ of this simple harmonic oscillator is given by, The structure factor in Eq. ( 27) thus becomes, The matrix element can be evaluated in the following way, where we use e −iqx ae iqx = a + iq √ 2m d ω0 in the second equality.Plugging the above matrix element to the structure factor in (33) becomes, where W toy d (q) = q 2 /(4m d ω 0 ) is the Debye-Waller factor in the toy model.The structure factor follows a Poisson distribution with mean number of phonons µ = q 2 /(2m d ω 0 ), as also shown in the case of the 3-dimensional harmonic oscillator in [41].

B. Perturbation theory for anharmonic oscillator:
q ≪ √ 2mω0 We now turn to more general case where small anharmonic terms are included in the 1D toy potential.An exact solution is no longer possible.But as we will see, in the kinematic regime q ≪ √ 2m d ω 0 , we can use perturbation theory to obtain the behavior of the structure factor and illustrate the importance of the anharmonic corrections as a function of momentum and energy deposition.Our goal in this section then is to obtain the power counting of the anharmonic contributions to the structure factor.
The toy Hamiltonian we consider is given by, We will concretely consider k equal to 3 and 4, corresponding to a leading cubic and quartic anharmonicity, respectively.Treating the dimensionless anharmonicity parameter λ k as a perturbation, the eigenstates |Φ n ⟩ are given by and E ′ n are the perturbed energies, With time-independent perturbation theory, the dynamic structure factor can be explicitly computed at different orders in λ k using (27).We defer the details of the explicit calculation to Appendix B. Instead, from the structure of the expansion we can already learn about the relevant corrections.In general, we can express the dynamic structure factor as an expansion in both λ k and q 2 /(2m d ω 0 ).At zeroth order in λ k , we see from (35) that the n-phonon term appears with a q-scaling of q 2n /(2m d ω 0 ) n .As we will show below, anharmonicity introduces departures from this q-scaling at higher orders of λ k .In the kinematic regime under consideration (q ≪ √ 2m d ω 0 ), powers of q 2 /(2m d ω 0 ) smaller than n can lead to large anharmonic corrections to the n-phonon term in the structure factor. 1 The aim of this section is thus to illustrate the behavior of the q-scaling at different orders of λ k .
The general expression for the dynamic structure factor in the toy model can be written as, For each n, the harmonic contribution appears at O((q 2 /(2m d ω 0 )) n ) as seen in (13); note that we do not 1 Perturbation theory in λ k is still valid.For instance, the expansion in (37) still holds.But the harmonic contribution in the structure function could be suppressed by small q for multiphonon states.
include the Debye-Waller factor in this power counting discussion since it always appears as an overall factor.The anharmonic corrections are included here as an expansion in powers of q 2 /(2m d ω 0 ) which are denoted by i.From the orthogonality of the states |Φ n ⟩ with the ground states, we see that the dynamic structure factor should vanish for q → 0, which in turn implies that i ≥ 1.Each power i of q 2 /(2m d ω 0 ) appears with non-zero powers of λ k , denoted by ν(n, i).Here the power ν(n, i) is the smallest allowed power of λ k for a given phonon number n and the power i of q 2 /(2m d ω 0 ).However, numerical cancellations can sometimes force this leading behavior to vanish.Typically, the bigger the difference in i and n, the larger the power of λ k that is required.We will explicitly see the behavior of the powers ν(n, i) for k equal to 3 and 4 below, but we first discuss the implications of this form.
For the single phonon structure factor (i.e. for n = 1), the anharmonic terms are always suppressed compared to the harmonic term because of the additional powers of λ k and q 2 /(2m d ω 0 ).But for phonon numbers n > 1, it is possible for anharmonic contributions to dominate for q ≪ √ 2m d ω 0 .As a simple example, in the 3-phonon state, the harmonic contribution to the structure factor is proportional to q 6 /(2m d ω 0 ) 3 , while the aharmonic result contains λ 2 3 q 4 /(2m d ω 0 ) 2 .So when q ≪ √ 2m d ω 0 , the anharmonic effect can lead to a large correction to the dynamic structure factor.
In a generic n-phonon state, the harmonic piece scales as (q 2 /(2m d ω 0 )) n .Comparing this with the anharmonic term ∝ λ ν(n,i) k q 2i /(2m d ω 0 ) i , we note that the anharmonic term dominates the harmonic term for q ≪ √ 2m d ω 0 λ ν(n,i)/(2(n−i)) k . For small enough q, the behavior is governed by the anharmonic effects.Of course, at even smaller q ∼ q BZ one would expect the incoherent approximation to break down.For the values of λ in realistic materials, we find that the dominance of the anharmonic terms can happen for q above q BZ , particularly for larger n.These corrections become larger with n since the harmonic piece is progressively more suppressed in q 2 /(2m d ω 0 ).
We now illustrate the origin of the λ k powers ν(n, i) with an example in the case of k = 3.In this case, the perturbation x 3 ∼ (a+a † ) 3 implies the leading correction to the state can change the oscillator number by ±1 or ±3.Then the perturbed eigenstates have the schematic form: We neglect the numerical prefactor in front of each state.
Note that the terms are only present if the integer labelling the state is non-negative, for example for the ground state |Φ 0 ⟩ ∼ |0⟩+ λ 3 (|1⟩+|3⟩)+O(λ 2 3 ).The matrix element appearing in the n-phonon structure factor can be expressed as, where the coefficients are schematically given by, In order for given term in the coefficient to be nonzero, a minimum number of powers of iqx are required in the series expansion for e iqx .This therefore links the powers of q with powers of λ 3 .
Taking n = 3 as an example, then b 0 ∝ (iq) 3 at leading order in the q expansion.Meanwhile, b 1 ∝ (iq) 2 +(iq) 4 +....Note that the matrix elements ⟨0|e iqx |0⟩ and ⟨3|e iqx |3⟩ in b 1 contain terms proportional to (iq) 0 , but they cancel each other, consistent with a matrix element that always vanishes as q → 0. Also note that the coefficients b 0 , b 1 always alternate in even or odd powers of (iqx) and therefore alternate in being purely real or imaginary.The resulting matrix element squared thus goes as For the cubic interaction, only even powers of λ 3 appear in the matrix element squared due to the alternating even and odd powers of (iqx) in the b coefficients.In this example, in order to achieve the minimum q scaling of q 2 , higher powers of λ 3 are required, which will introduce more terms in the expansion.Here we see a correction to the matrix element squared at O(q 2 λ 4 3 ).The explicit derivation of ν(n, i) is given in Appendix B. The minimum power of λ 3 required to get the leading behavior ∝ q 2 /(2m d ω 0 ) in the anharmonic terms is given by, The minimum power of λ 3 as a function of the phonon number n and the power i of q 2 /(2m d ω 0 ) for i > 1 is given by, We show the expansion of the structure factor in the powers of λ 3 and q 2 /(2m d ω 0 ) schematically in Fig. 4, where we drop the numerical coefficients for all the terms and only illustrate the behavior of the powers of λ 3 and q 2 /(2m d ω 0 ).In the right part of the schematic, we show the behavior of the n-phonon term for n > 3, and in the left part of the schematic, we show the expansion for n = 1, 2, and 3.
The relationship between the powers in λ 3 and the powers of q 2 /(2m d ω 0 ) in (46) can also be understood in the following way.The powers of q 2 /(2m d ω 0 ) that appear at O(λ ν 3 ) can range from n − 3ν/2 to n + 3ν/2, with the minimum power allowed being 1, and ν being an even positive integer.Contributions from powers larger than n are suppressed in the kinematic regime q ≪ √ 2m d ω 0 .But powers smaller than n can lead to significant corrections in the same regime.
For example, the anharmonic contribution to the 2-phonon structure factor has a leading behavior ∝ λ 2 3 q 2 /(2m d ω 0 ), which is expected to dominate the harmonic behavior ∝ q 4 /(2m d ω 0 ) 2 for small enough q (explicitly for q ≲ √ 2m d ω 0 λ 3 ).Assuming m d ∼ 28 GeV, ω 0 ∼ 40 meV, and a typical value of λ 3 ∼ 0.01, we expect the anharmonic contribution to start to dominate for q ≲ 0.5 keV.This kinematic regime does not strictly satisfy the conditions for the incoherent approximation which are assumed in this calculation.However, it is interesting to note here that the size of this anharmonic correction roughly matches onto the result for the 2-phonon structure factor in the long-wavelength limit (q ≪ 1/a) [25,37], where it was found that anharmonic interactions give up to an order of magnitude correction to the structure factor.At the edge of the Brillouin Zone q ∼ 2π/a ∼ O(keV), with the typical values used above, we find in the toy model an O(∼ 25%) correction at the boundary of the valid region for the incoherent approximation.
For k equal to 4, which corresponds to a quartic perturbation to the harmonic potential, the calculation proceeds similarly to the cubic case discussed above, except for some key differences.All the coefficients b i are either real or imaginary based on whether n is even or odd respectively, and hence the anharmonic corrections appear in all orders of λ 4 .We thus have corrections at O(λ 4 ).For even n, coefficients b i only have even powers of q, and thus cannot generate terms ∝ q 2 in the squared matrix element.The leading behavior for even n is thus ∝ q 4 .For odd n however, the leading behavior is ∝ q 2 , and the FIG. 4. Expansion of the structure factor in phonon number n, powers of q 2 /(2m d ω0), and powers of λ3 for a cubic perturbation (k = 3 in ( 36)).The right part shows the general behavior of the n-phonon term for n > 3, while the left part shows the expansion for n =1, 2, and 3. Shaded terms show the dominant contributions when q ≪ √ 2m d ω0, which comes from the anharmonic terms for n ≥ 2.Here we just illustrate the power counting; individual terms might not be present if there is a numerical cancellation in the coefficients.FIG. 5. Expansion of the structure factor in phonon number n, q 2 /(2m d ω0), and λ4 for a quartic perturbation (k = 4 in ( 36)).The right part shows the general behavior of the n-phonon term for n > 3, while the left part shows the expansion for n =1, 2, and 3. Shaded terms show the dominant contributions when q ≪ √ 2m d ω0, which comes from the anharmonic terms for n > 2. Similar to the above, individual terms might not be present if there is a numerical cancellation in the coefficients.minimum power of λ 4 is given by, For powers i greater than 1, the minimum power of λ 4 for any phonon number n is given by, We show the expansion of the structure factor in the powers of λ 4 and q 2 /(2m d ω 0 ) schematically in Fig. 5, where we drop the numerical coefficients for all the terms and only illustrate the behavior of the powers of λ 4 and q 2 /(2m d ω 0 ).Similar to Fig. 4, we are only illustrating the minimum allowed powers of λ k in perturbation theory for n > 3. Due to numerical cancellations, the leading λ k power can vanish in some cases.

Limitations of perturbation theory
Our analysis has focused on the regime q ≪ √ 2m d ω 0 because this corresponds to a low mean phonon number.For large enough n, perturbation theory will start to break down.Equivalently, for a given n, perturbation theory will only be valid for λ k sufficiently small.
For a particular phonon number n, if the energy correction in (38) is of the same order as the unperturbed energy eigenvalue (n + 1  2 )ω 0 , the perturbation can no longer be treated as small.Based on this, we set an up-per bound on |λ k | by requiring that At leading order, the correction for k equal to 3 (i.e. a cubic perturbation) is given by The equivalent result for k = 4 reads, Using the equations above, we get the critical value of λ 2 3 and λ 4 compatible with the perturbation theory expansion.These are shown in Fig. 6.With the analytic structures of the energy corrections shown above, we see that the perturbativity bound on λ 2 3 (λ 4 ) has a scaling ∝ 1/n 2 (∝ 1/n), where n is the phonon number.For typical values of λ 3 ∼ 0.01, we see that the perturbation theory is valid only up to n ∼ 6−7.Furthermore, perturbation theory is impractical for calculating corrections at small q and very high phonon number n, since these corrections will be a very high order in the anharmonicity parameter.
To deal with these limitations, we consider two different approaches in this paper.Since high n is associated with high ω and q, in the next section we will use the impulse approximation to account for anharmonic effects at high q.In Appendix D, we also study a special anharmonic potential, the Morse potential, where it is possible to obtain exact results.We use this as a case study to validate the perturbation theory and impulse approximation results.

C. Impulse Approximation for q ≫ √ 2mω0
As we have shown, perturbation theory quickly goes out of control beyond the first few number of phonons.Resumming the anharmonic interaction is usually needed for the structure factor when q or ω is large.Consider the following phase space Impulse regime: q ≫ √ 2mω 0 , It has previously been shown [25,50] in the harmonic case, that one can calculate the structure factor by using a saddle point approximation in the time-integral representation of the structure factor.This is called the "impulse approximation" since the steepest-descent contour is dominated by small times, which can be interpreted physically as an impulse.
We begin with the structure factor in Eq. ( 27), which can be decomposed as contributions from each atom d, S toy (q, ω) = d n d |f d | 2 S toy,d (q, ω).Then we rewrite the energy conservation delta function as a time integral S toy,d (q, ω) where in the second equality we use the fact that |Φ 0 ⟩ and |Φ f ⟩ are eigenfunctions of H, and in the third equality we use the completeness relation and the time-dependent position operator x(t) = e iHt xe −iHt .The final expression is the well-known structure factor in the time domain.
Using the above representation of the structure factor, We can further simplify this using the fact that e iqx acts as a translation operator on momentum p, e −iqx p e iqx = p + q.Applying the translation on the full Hamiltonian yields Here we generalize the impulse approximation to any 1D Hamiltonian, H(x, p) = p 2 2m + V (x), which satisfies One can also generalize impulse approximation to a generic potential V (x, p) as long as the above holds in the limit of large q. 2 In other words, we require that 2 In this case, the impluse regime in Eq. ( 52) needs to be replaced as ω ∼ q 2 2m + q m ⟨p⟩ and we impose Eq. ( 56) holds up to O ω 2 0 /q correction.
the Hamiltonian in the large momentum limit is dominated by the kinetic energy p 2 2m , not the potential.We can then obtain reliable theoretical predictions in the impulse regime even with large number of phonons.
Applying the above to Eq. ( 54), the structure function now reads where we translate the momentum in the first line and use Eq. ( 56) in the second line.Note that H = H(x, p) throughout and we drop the argument for brevity.The last line is exact for potentials that depend only on x.Now we can apply the saddle point approximation to evaluate the time integral.Defining H ′ ≡ H + pq m , we can write where In order to calculate this object, we can expand ln⟨e iH ′ t ⟩ in small t.The first few terms in this expansion are given by . . .
In the harmonic approximation, only the terms proportional to q 2 are nonzero.As a result, only the first few expansion terms are needed as long as t ≪ 1 ω0 since f (n+1) /f (n) is of order ω 0 .Then one can solve for the saddle point t I by solving f ′ (t I ) ≈ f ′ (0) + f ′′ (0)t I = 0, which gives where In the last equality we use the fact that ⟨p⟩ = 0 for a V (x) potential since ⟨p⟩ ∝ ⟨[x, H]⟩ = 0.Although t I is formally imaginary, its magnitude is small and close to the origin in the impulse regime.Since there is no pole around this saddle point, we can approximate the time integral by the saddle point and find For large energy depositions the Gaussian becomes narrowly peaked around ω = q 2 /2m, and this reproduces the nuclear recoil limit [25].
In the presence of anharmonic interactions, other powers of q will be present in the expansion of (60).In general, the f (n) term will have a q n term with coefficient of O(λ).In this case, f (n+1) /f (n) ∼ q ω 0 /m.Higher orders will then be important in the expansion of f (t) for sufficiently large q or t.For a given q, the higher order corrections become relevant for |t| ≳ m/ω 0 /q ∼ 1/ √ ωω 0 in the impulse regime.Including these corrections is difficult in general, but we can continue to use the second order expansion giving (63) as long as |t| ≲ 1/ √ ωω 0 .According to (61), this corresponds to a condition on how close ω is to q 2 /(2m).Since q 2 ∼ 2mω and σ 2 p ∼ mω 0 , this implies that We see that the distance of ω from q 2 2m sets the size of t I , which in turn tells us the regime for the validity for the approximation (63).The condition (64) is approximately the same condition that ω is within the Gaussian width in (63), and keeping terms in f (t) only up to f ′′ (0) is self-consistent near ω = q 2 2m .Therefore, in the presence of anharmonic interactions, the above structure factor result (63) remains valid in the impulse regime (52).The only modification is in σ 2 p .
Considering perturbations in V (x) up to x 4 and recalling that the expectation value is with respect to the full ground state, we find that at leading order in λ 3 , λ 4 .The nuclear recoil limit is again reproduced, with a small modification to the width of the Gaussian envelope due to anharmonic couplings.Note that in order to calculate the structure factor far from ω = q 2 2m , we must include additional orders in f (t) and t I .We do not perform these higher order calculations for the final results in this paper since they have a negligible effect on the integrated rates, but we provide the procedure for completeness in App. C.
Finally, we approximate the effect that introducing the full crystal lattice has on this single atom result.Up until the evaluation of various moments of H ′ , the impulse approximation is fully model-independent.We just have to make an adjustment to the final evaluation of ⟨p 2 ⟩.The states in the full crystal theory are smeared by the phonon density of states, so we calculate ⟨p 2 ⟩ via the following prescription where g(λ) is the anharmonic correction calculated in the single-atom potential.Essentially, we have used the average single phonon energy to calculate ⟨p 2 ⟩.In the harmonic limit, (63) then exactly matches the impulse result from [25].In summary, in this section we have demonstrated the general behavior of anharmonic effects with q and ω.We have shown that they are indeed negligible at high q and ω ∼ q 2 /2m d , consistent with the intuition that scattering can be described by elastic recoils of a free nucleus.The effects grow for q ≪ √ 2m d ω 0 and at low q they may dominate the structure factor.This roughly matches onto the results of Refs.[25,37], which found that for q < 2π/a anharmonic effects can have a large impact on the twophonon rate .

IV. NUMERICAL RESULTS FOR 1D ANHARMONIC OSCILLATOR
Having demonstrated the analytic behavior of the dynamic structure factor in the previous section, we now turn to obtaining numerical results using realistic potentials.We will perform concrete calculations for Si and Ge as representative materials while briefly commenting on others.As discussed in Sec.II B, we adopt an empirical model of interatomic interactions that encodes the anharmonicity in the potential.We use this empirical model to calculate a single atom potential, which we then use to evaluate the structure factor numerically.
As stated in Sec.II C, we start by fitting the single atom potential in a particular direction onto a 1D potential of the form, In the fit, ω 0 , λ 3 , λ 4 are free parameters but in order to reproduce the harmonic limit, we then make the replacement ω 0 = 1/ω −1 , which is calculated from the phonon density of states and gives a slightly different numerical value.This is motivated by the harmonic case discussed in Sec.II A. We do not consider anharmonic terms ∝ x k for k ≥ 5 as we observe that the anharmonic potential along any direction is dominated by the cubic and the quartic terms.We find that the maximum anharmonicity is typically along the nearest neighbor direction (x, y, z) = (1, 1, 1).For computing results, we will consider the potential along this direction, which represents maximum anharmonicity, as well as the potential in an orthogonal direction (x, y, z) = (1, −2, 1), which represents an intermediate value for the anharmonicity.Using the aforementioned interatomic models, we find anharmonicity strengths ranging from λ 3 ∼ 6 × 10 −3 to 10 −2 and λ 4 ∼ (2 − 3) × 10 −4 .For Si and Ge, the results are same for either atom in the unit cell.
Given the 1D potential in (67), we find exact solutions of the 1D eigenvalue and eigenvector problem using a simple finite difference method.We take a first order discretization of the Laplace operator and solve the discretized time-independent Schrödinger equation in a box.The box grid interval size must be small enough to resolve the maximum momentum scales of interest, which in this case depends on the highest excited state needed in the calculation.Also, the minimum box size required depends on the spatial extent of the highest excited state used.As seen in Sec.III C, the impulse approximation suffices for q > O(few) × √ 2m d ω 0 .Beyond this momentum, we no longer need to calculate excited states since the structure factor in the impulse limit is independent of the details of the highly excited states.The nth excited state is most relevant at momenta q ∼ √ n √ 2mω 0 .Therefore, to complete our calculation below the impulse limit, we include the first 10 excited states.The results for these eigenstates are converged above a box size of ∼ 10/ √ 2mω 0 and grid size of ∼ 0.1/ √ 2mω 0 .We now use these numerical eigenstates and energies to calculate the structure factor in Eq. ( 27).We apply a prescription for the energy-conserving delta function similar to that used in the harmonic 1D oscillator, Eq. ( 15).The final result at momenta below the impulse regime (q < 2 √ 2mω 0 ) is, where and f (n), |Φ 0 ⟩, |Φ f ⟩ are given by the numerically solved eigenenergies and eigenstates, respectively.D(ω) is the single phonon density of states calculated with DFT [51].
In this work we assume equal couplings of DM with all nucleons so that f d = A d , where A d is the atomic mass number.In the equations above, we have included a sum over all atoms in the unit cell d with density n d , and in general the atomic potentials and density states can also depend on d, although for Si and Ge we do not include this.
In the impulse regime (q > 2 √ 2mω 0 ), we have shown in Sec.III C that the structure factor for any positiondependent potential is approximated by a Gaussian envelope, where the the expectation values are all computed in the ground state and adjusted to the average single phonon energy via (66).Now we simply use the numerical ground state of the anharmonic potential (67) to calculate ⟨p 2 ⟩ and therefore obtain the structure factor.Note that the anharmonic contribution is essentially negligible in the impulse limit, since corrections to ⟨p 2 ⟩ are ∝ λ 2 3 , λ 4 .Fig. 7-8 shows numerical results on the structure factor for Si and Ge, taking the maximum anharmonicity in either case.In Fig. 7, the structure factor as a function of q is shown.As ω (and therefore minimum phonon q-dependence of structure factor: We the structure factor in the harmonic and anharmonic cases, where in the latter case the structure factor is calculated numerically with the maximal anharmonicity.The lines from top to bottom show the structure factor at different ω, corresponding to an increasing minimum phonon number n.There are large corrections for q ≪ √ 2m d ω0 when anharmonic interactions are included (dashed), and the corrections become more significant as the threshold is increased.For q ≫ √ 2m d ω0, both cases converge to the same result.For Si, we have √ 2m d ω0 ≈ 40 keV while for Ge, √ 2m d ω0 ≈ 50 keV.For other materials, this quantity is listed in Table I.The incoherent approximation momentum cutoff is qBZ < 2π/a ∼ 2.2 keV for both crystals.
Si, Multiphonon Structure Factor and Impulse Approximation FIG. 8. ω-dependence of structure factor: For different q values, we show the decomposition of the structure factor into individual n phonon terms, where the energy-conserving delta function has been smeared as in (68).Note that the maximum anharmonicity has been included in the numerical calculation, but the result is nearly identical to the harmonic result for these q values, as shown in Fig. 7.The dotted line shows the impulse approximation, which starts to become a good approximation as q increases above √ 2m d ω0.
number n) is increased, there is a larger anharmonic correction at small q.This can be understood by looking at the q scalings discussed in Sec.III B and illustrated in Fig. 4 and Fig. 5.At low q and thus DM mass, the contributions from the anharmonic structure factor can give smaller powers of q 2 2m d ω0 compared to the leading harmonic term , so the enhancement grows with n.At high q, results converge to the harmonic result, consistent with our discussion of the impulse regime in Sec.III C. We see this also in Fig. 8, which shows the structure factor at different q.The impulse approximation becomes better as q ≫ √ 2m d ω 0 , and is indistinguishable from the harmonic case.

A. Impact on DM scattering rates
We now use the numerical results for the structure factor to compute the DM scattering rates for a range of DM masses and experimental thresholds.Our results are summarized in Figs.9-10.We consider DM masses in the range ∼ 1 − 10 MeV.The lower end of the mass range is chosen such that the momentum transfers are large enough to satisfy the condition for the incoherent approximation (i.e.q > 2π/a), while at the upper end of masses it is expected that scattering is described by the impulse approximation [25].It is precisely this mass range where details of multiphonon production are important.We will also consider the two cases of scattering through heavy and light mediators.The goal will be to identify the region of parameter space where the anharmonic effects on the dynamic structure factor affect the scattering rates the most.
In the isotropic limit, the observed DM event rate per unit mass is given by [ where ρ χ is the DM energy density, ρ T is the mass density of the target material, m χ is the DM mass, µ χ is the DM-nucleon reduced mass, σ p is the DM-nucleon cross section, and f (v) is the DM velocity distribution.The structure factor S(q, ω) is given by our numerical results (68)-(72) and the integration bounds are determined by the kinematically allowed phase space where the energy threshold of the experiment is denoted by ω th .The q-dependence of the DM-nucleus interaction can be encapsulated in the DM form factor F (q), where F (q) = 1 indicates an interaction through a heavy mediator, and F (q) = q 2 0 /q 2 indicates an interaction through a light mediator for a reference momentum transfer of q 0 .Note that in general, the strength of the anharmonicity varies with the direction of the recoil of the nucleus, and the structure factor will depend on the direction of the momentum transfer.For simplicity, we are assuming that the anharmonicity strength is uniform in all directions.Our estimate with the maximum anharmonicity thus provides an upper bound on the anharmonic effects on DM scattering.
The DM mass sets the typical momentum-transfer scale q of the scattering, and the experimental energy threshold ω th sets the phonon number n.Hence, to identify the DM masses and experimental thresholds where anharmonic effects start to become important, we first need to understand the q-values where the anharmonic corrections are large for a particular phonon number n.We can estimate this using the perturbation theory results in Sec.III B. Note that in our numerical calculation, we find that λ 3 generally provides the larger anharmonic contribution, so we will focus on a purely cubic perturbation in this discussion.
For the analysis of a cubic perturbation discussed in Sec.III B, we showed that anharmonic effects introduced additional terms to the n-phonon structure factor of the form ∝ λ ν(n,i) 3 , see (39).Therefore when q is lower than the scale terms in the anharmonic structure factor can be of comparable size to the harmonic structure factor.In order to find the largest q-scale where the anharmonic contribution starts to become relevant, we can evaluate (76) for all positive i < n, and find the minimum possible exponent of λ 3 .For n = 2 or 3, the minimum exponent is achieved for i = 1, for which ν(n, 1) = 2.This gives a q-scaling of q ∼ √ 2m d ω 0 λ 1/(n−i) 3 . This tells us that for the 2-phonon case, the anharmonic contribution should begin to become important at q ∼ √ 2m d ω 0 λ 3 , while for the 3-phonon case, the anharmonic contribution becomes important at q ∼ √ 2m d ω 0 λ 1/2 3 .For a larger number of phonons, this scaling is approximately q ∼ √ 2m d ω 0 λ 1/3 3 .So we see that higher energy excitations have more significant anharmonic contributions at larger momentum transfers.Below the q-scale identified above, the anharmonic contributions are expected to increase substantially with decreasing q, as terms ∝ q 2i for i < n dominate the harmonic scaling ∝ q 2n .We now recast our analysis concretely in terms of DM mass and experimental energy thresholds as follows.For both massive and massless mediators, the event rate for n ≥ 2 phonons is always dominated by the large q portion of phase space and energy depositions near the threshold.Therefore the enhancement in the rate due to the anharmonicity roughly corresponds to the enhancement in structure factor evaluated at S(q = 2m χ v, ω = ω th ), where v is the DM velocity.Inserting q = 2m χ v into the condition in (76) gives a condition on the DM mass: where 10 −3 is the typical DM velocity.In order to determine the appropriate phonon number n for a given ω th we must take into account the subtlety that each excitation energy is smeared across a width, as discussed in Sec.II C and also given in (70).To solve for the smallest n that contributes appreciably above ω th , we solve the following equation: where σ is the single-phonon width as defined in (70) and we have for simplicity taken f (n) = n.Applying (77)-(78) to Si with ω 0 = 31 meV, σ = 18 meV, and m d = 26 GeV, we find the following results Below these masses, anharmonic corrections become large.The last line applies for thresholds above 160 meV which corresponds to n ≥ 4, and these n-phonon terms all give the same condition on DM mass.Note that this is only a heuristic, which does not include for example the combinatorial pre-factors or cancellations in the perturbation theory calculation.Nonetheless, we do see the same qualitative features in the complete numerical result which is given in Fig. 9.
In order to generalize (79) to other materials, we give the necessary energy scales in Tab.I. Despite large differences in ω 0 , the momentum scale √ 2m d ω 0 ends up being about the same in all crystals.Then the typical DM mass scale for anharmonic effects to become important is also about the same for a fixed phonon number n.However, the differences in ω 0 mean that the threshold corresponding to a given n can vary significantly.For a given threshold, GaAs and Ge have the largest phonon number.Since anharmonic corrections become more important with larger n, GaAs and Ge will therefore have larger anharmonic contributions compared to Diamond at the same threshold.
In Fig. 9, we present the ratio of scattering rates in the anharmonic case to the harmonic case in Si and Ge, taking two representative cases for the couplings.We also present the cross-sections corresponding to an observed rate of 3 events per kg-yr in Fig. 10.The bands depict the possible uncertainty that anharmonicity introduces to an experimental reach, with the solid line giving the harmonic result and the dot-dashed the result for maximal anharmonicity.We do not show the effects above the cross sections of σ n ≳ 10 −28 cm 2 as for these large interaction strengths, the DM is expected to lose a significant energy in 1 km of Earth's crust through scattering, thus rendering DM with such cross sections unobservable in underground direct detection experiments [52].
For m χ > 10 MeV, the typical q becomes similar or larger than √ 2m d ω 0 , where there is negligible difference in the anharmonic and harmonic structure factors.The rates will also start to be dominated by the impulse regime q ≫ √ 2m d ω 0 .In this case, the structure factor calculated with an anharmonic potential is nearly identical to that calculated in the harmonic case, as discussed in Sec.III C. We have also seen this behavior with numerical computations in Fig. 8.The anharmonic and harmonic scattering rates are also essentially identical for DM masses m χ > 10 MeV.
For DM masses m χ < 10 MeV (i.e.q < √ 2m d ω 0 ), the ratio of the anharmonic to harmonic rate begins to grow with decreasing DM mass.As the typical q decreases with decreasing DM mass, the leading anharmonic term ∝ q 2 2m d ω0 grows faster compared to the harmonic term ∝ q 2 2m d ω0 n for n ≥ 2. The effect is more pronounced for higher thresholds or equivalently higher n, since the harmonic term is even more suppressed.Therefore at larger thresholds, the anharmonic effects start becoming important already at larger masses and also grows much more quickly as the DM mass is decreased.For a given DM mass, this also implies that the spectrum of events will have larger anharmonic corrections on the high energy tail of events.However, the rates are also highly suppressed in this tail, and only observable for high scattering cross sections.
At DM masses m χ < 1 MeV, the slope of the ratio of the anharmonic rate to the harmonic rate starts to decrease slightly, which is an artifact of the Brillouin zone momentum cutoff that we apply across all rate calculations.The incoherent and subsequent approximations are not guaranteed to be justified in this regime, so this effect should not be treated as physical.For sub-MeV DM masses, the phonons again should be treated as collective excitations, similar to the calculation of Ref. [37].
Lastly, we note an interesting feature that the anharmonic scattering rate is strictly greater than the har- Comparison of the cross section corresponding to 3 events/kg-yr in the harmonic (solid) and anharmonic (dot-dashed) cases.The anharmonic result is shown for maximal anharmonicity, and so the shaded band represents our estimate of the theoretical uncertainty due to anharmonic effects.The effects are primarily important for high thresholds and low DM masses, corresponding to large σn, which is generally in tension with existing astrophysical or terrestrial constraints.monic rate in the entire parameter space that we probe.This is a consequence of the sign of the leading q-scaling term q 2 2m d ω0 .For the production of an excited state |Φ f ⟩ in the crystal, the term in the dynamic structure factor ∝ q 2 can only come from the term |⟨Φ f |iqx|Φ 0 ⟩| 2 , as the 2 |Φ 0 ⟩ * and its conjugate are zero from orthogonality.Thus, the sign of the term ∝ q 2 in the anharmonic structure factor is strictly positive for producing an excited state, whereas there is no corresponding term ∝ q 2 in the harmonic case for n ≥ 2 phonons.Since we are probing the q ≪ √ 2m d ω 0 regime, this leading term quickly dominates the structure factor.Thus, the anharmonic scattering rate exceeds the harmonic rate in this regime.A consequence of this is that we expect the harmonic crystal result gives a lower bound on the scattering.

V. CONCLUSIONS
Scattering of DM with nuclei in crystals necessarily goes through production of one or many phonons for DM masses smaller than ∼ 100 MeV.Previous work has focused on calculating the multiphonon scattering rates in a harmonic crystal under the incoherent approximation (i.e.q > q BZ or DM mass ≳ MeV).In this work, we have studied the effects of anharmonicities in the crystal on the scattering rates, while still working within the incoherent approximation.
In order to obtain a tractable calculation of anharmonic effects, we have simplified the problem into a toy model of a single atom in a 1D anharmonic potential.In this toy model, scattering into multiphonons can still be well-approximated by applying a smearing on the spec-trum of quantized states to account for the phonon spectrum of a lattice.We extract anharmonic couplings by modeling the interatomic potentials of Si and Ge, which give rise to realistic single atom potentials.This approach allows us to obtain an analytic understanding and first estimate of the impact of anharmonicity, although the numerical results should not be taken as a definitive rate calculation.
We find that the harmonic crystal results of Ref.
[25] can be safely assumed for DM masses down to ∼ 10 MeV.Below ∼ 10 MeV, this assumption cannot be taken for granted.In this regime, we find that anharmonic effects on the scattering rates increase with decreasing DM mass and increasing experimental thresholds.Anharmonic corrections up to two orders of magnitude are possible for DM masses ∼ a few MeV and for experimental thresholds ∼ a few times the typical single phonon energy of the crystal.These findings are consistent with Refs.[25,37], which studied two-phonon production from sub-MeV DM and found up to an order of magnitude larger rate from anharmonic couplings.
The size of the corrections is dependent on the material through the anharmonicity strength of that crystal and also, non-trivially, through the typical single phonon energies of the material.For a particular energy threshold, crystals with lower single phonon energies exhibit larger corrections since they require larger phonon numbers to be produced.For example, anharmonic effects in Ge can be larger by almost an order of magnitude than those in Si for similar DM parameter space and thresholds, even though the anharmonic couplings in the two crystals are similar.This is a consequence of the difference in q scaling of the harmonic and anharmonic contributions, which become more pronounced with larger phonon number.Materials with low single-phonon energies, such as GaAs and Ge, therefore have the largest anharmonic effects.The effects will be reduced in Diamond and Al 2 O 3 , which have even higher single phonon energies than Si.
The relevance of anharmonic effects to direct detection experiments depends on the DM cross section.The effects are largest for low DM masses and high thresholds, in other words on the tails of the recoil spectrum where the rates are small.For a typical benchmark exposure of 1 kg-yr, the anharmonic corrections become sizeable for DM-nucleon cross sections above ∼ 10 −34 cm 2 .Being agnostic about any terrestrial or astrophysical constraints on the DM model and only requiring the DM to be observable in underground direct detection experiments, the upper bound on the DM cross section is σ n ≲ 10 −28 cm 2 [52].This comes from considering an overburden of ∼ km.On the other hand, these very high DM-nucleon cross sections are typically excluded by terrestrial and astrophysical constraints for the simplest sub-GeV dark matter models [53,54].DM-nucleon cross sections σ n ≳ 10 −41 cm 2 (σ n ≳ 10 −31 cm 2 ) are constrained for typical models with a heavy mediator (light dark photon mediator) for a DM mass ∼ MeV.With these constraints, we see from Fig. 10 that the anharmonic effects can only impart corrections of at most an order of magnitude for experiments with kg-yr exposure.
Experiments with exposures above kg-yr could see larger anharmonic effects, since they would be more sensitive to the events at high phonon number for MeVscale DM.However, for solid-state direct detection experiments, achieving exposures significantly bigger than a kg-yr is challenging.Thus, for near-future crystal target experiments, we conclude that the anharmonic effects are only important up to O(1) factors at masses of ∼ a few MeV for the simplest DM models.
where the coefficients b j are given by, In general, the coefficient b j is schematically given by, 0 ⟩ + ... + ⟨ψ (1)  n |e iqx |ψ To study the powers of q appearing in b j , we first need to understand the structure of the matrix element ⟨n 1 |e iqx |n 2 ⟩ for general eigenstates |n 1 ⟩ and |n 2 ⟩ of the unperturbed harmonic oscillator.This matrix element is given by the following, We learn that the matrix element ⟨n Note again that the Debye-Waller factor e − q 2 4m d ω 0 is not included in this power counting since e − q 2 4m d ω 0 ≈ 1 in the regime of interest.
Combining this information with the structure of b j in (B8) and the structure of |ψ (j) n ⟩ in (B5), the powers of q in b j can be identified: Note that only those terms with powers of q larger or equal to 1 are present.Terms ∝ q 0 have to cancel as they otherwise lead to q 0 terms in the squared matrix element |⟨Φ n |e iqx |Φ 0 ⟩| 2 , which is forbidden due to orthogonality of eigenstates.
As the kinematic regime under consideration is of q ≪ √ 2m d ω 0 , we will focus on powers of q less than n, which corresponds to the harmonic case.We see from the equation above that the lowest powers of q decrease with increasing values of j.Thus, higher order corrections in λ k appear with lower powers in q.Eventually, at a sufficiently high power of λ k , we get a coefficient b j with the minimum power of q equal to 1.The squared matrix element can then be written in general as, where the first term on the right hand side ∝ q 2n is the harmonic term, and the anharmonic corrections are expanded in powers of q 2 which are denoted by i, with i ≥ 1.Every power i appears with a minimum allowed power ν(n, i) of λ k .
To study the behavior of ν(n, i), we first note that, for even k, the matrix element ⟨Φ n |e iqx |Φ 0 ⟩ is purely real or purely imaginary, depending on whether n is even or odd respectively.For instance, if n is even, then b 0 is purely real.Higher orders in λ k lead to insertions of (a + a † ) k and therefore matrix elements where the difference in the harmonic oscillator states is also even, so that all coefficients b j are real in this case.But for odd k, the b j coefficients will alternate in being real and imaginary.This changes the structure of the squared matrix element depending on k, as we will see below.
Odd k: We will first consider odd k.In this case, the squared matrix element can be written as, Thus we see that we get corrections at even orders in λ k , with the lowest non-zero power being λ 2 k .In general, at O(λ j k ) for an even j = 2j ′ , the lowest power of q 2 is n − (j ′ × k), and the highest power is n + (j ′ × k).Note that only terms with positive powers of q 2 are present.The term ∝ q 2 can also subtly cancel in some cases as there is no term ∝ q 0 in coefficients b j .We will deal with this case later below.But to get a power i > 1 of q 2 , the lowest non-zero j ′ is ⌈ |n−i| k ⌉, with the lowest j given by 2 × ⌈ |n−i| k ⌉.Thus, in the squared matrix element, the lowest non-zero power ν(n, i) required is given by, To get the lowest power i = 1 of q 2 i.e. the term ∝ q 2 , the only possible way is to get the term ∝ q 1 in the coefficient b j as there is no term ∝ q 0 .For odd n, the term ∝ q 1 in b j can only be generated at an even j, since that is the only way to satisfy n − jk = 1.For every even j = 2j ′ , the powers of q in b j range from n − (2k) × j ′ to n + (2k) × j ′ .The lowest j ′ to get a term ∝ q 1 is then given by ⌈ |n−1| 2k ⌉, with j given by 2×⌈ |n−1| 2k ⌉.For an even n, the term ∝ q 1 in b j can only be generated for an odd j.For every odd j = 2j ′ −1, the lowest power of q in b j is n + k − (2k) × j ′ .The lowest j ′ to get a term ∝ q 1 is then given by ⌈ |n+k−1| 2k ⌉, with j given by 2 × ⌈ |n+k−1| 2k ⌉ − 1.In the squared matrix element, the lowest non-zero power ν(n, 1) required is given by, 2m d ω 0 Thus we see that we get corrections at all orders in λ k , with the lowest non-zero power being λ k .In general, at O(λ j k ), the lowest power of q 2 is n − (j × k)/2, and the highest power is n + (j × k)/2.Following similar arguments to the case of odd k discussed earlier, ν(n, i) for i > 1 is given by, Another difference between the case of even k considered here and that of odd k is that we do not get an i = 1 term for even n, as all terms in the coefficients b j contain even powers of q.This means that the leading term will always go as q 4 , with a λ k power determined by (B20) for i = 2.For odd n, the lowest power of q in b j is n − k × j.Thus, in the squared matrix element, the lowest non-zero power ν(n, 1) required is given by, The calculations in this appendix up to this point consider the overall scaling behavior of the powers of q 2 and λ k in the squared matrix element.We have neglected combinatorial factors at several steps in the calculations that enter into the numerical coefficients a n,i in (B11).Sometimes, the numerical coefficients can also cancel with each other, and the naive leading behavior estimated in this section can vanish.In order to give concrete examples of the numerical coefficients, we perform explicit calculations of the squared matrix element using perturbation theory with k = 3 (i.e. a cubic perturbation), and phonon numbers n = 1, 2, 3, and 4. We perform this explicit calculation only up to O(λ 23 ).The results of various numerical coefficients are presented below.
For a single-phonon production (i.e.n = 1), the coefficients a n,i are given by, Note that the coefficients a 4,1 and a 4,2 amount to zero because of a numerical cancellation between the two terms in the b 1 coefficient in Eq.B7.The leading behavior of the terms proportional to q 2 and q 4 in the structure factor is instead q 2 λ 6 3 and q 4 λ 4 3 , respectively.As these numerical coefficients appear through combinations and interferences of several combinatorial factors at various steps of the calculation, it is hard to provide a general expression for them.By looking at the examples above however, we can make some general observations.Typically, we see that the coefficients follow a pyramid structure, with a n,i being the largest for i near n, and decreasing with i away from n.We also find that the coefficients can vary by orders of magnitude from each other.The terms with i near n receive contributions from several individual matrix elements, and in general seem to be larger.We expect to see this pattern continue for higher phonon numbers as well.The exact values of these coefficients play a role in determining where the anharmonic corrections dominate, and so our power counting approach only gives an O(1) estimate.

Appendix C: Impulse approximation
In Sec.III C, we calculated the structure factor via the saddle point approximation in the regime defined by (52).This regime corresponded to values of ω near q 2 2m and within the Gaussian width of (63).As discussed in the main text, in order to calculate the tail of the structure factor far from ω = q 2 2m , more expansion terms are needed in f .Here we discuss this extension of the impulse approximation.
First, in the special case of a harmonic potential, we can start from the full result in Eq. (35).After rewriting the energy conservation delta function as a time integral, we find that f (t) = −iωt + q 2 2m d ω 0 (e iω0t − 1) (C1) Solving f ′ (t) = 0 gives the exact result Using the saddle point approximation for ω ≫ ω 0 , we find S toy,d (q, ω) ∼ 1 √ ωω 0 e −2Wtoy(q) q 2 2mω ω ω 0 e ω ω 0 .(C3) The same result can also be derived by approximating the sum over phonon states as an integral in Eq. (35).
The saddle point approximation for the harmonic oscillator holds as long as ω ≫ ω 0 , and we no longer have a condition on how close ω is to q 2 2m .In the impulse regime, ω ∼ q 2 2m , one can check that it reduces to the previous result in Eq. ( 63).We see in this exact result that the tail at large ω is Poissonian instead of Gaussian.
For general potentials, this exact analytic result is no longer possible, but we can still calculate corrections to the tail.First, we start by giving the exact saddle point equation: which is valid at all orders.We begin by noticing that saddle point equation (C4) is satisfied exactly at ω = q 2 2m by t I = 0.Then, ω-derivatives of t I at ω = q 2 2m can be found by taking ω-derivatives of (C4) and solving for t ].This allows us to calculate t I [ω = q 2 2m ] in an iterative fashion.The first few terms are I [ where t (n) I denotes the nth ω-derivative of t I .In the harmonic case, this series resums to (C2).For general potentials, one can then use the expansions (60) and (C5) to calculate S toy,d (q, ω) ≈ 2π −f ′′ (t I ) e f (t I ) (C6) to a desired order.Comparison of analytic structure factor in the Morse potential and the numerical calculation for Si as described in Sec.IV.We find that the two methods give almost the same result due to the fact that the Morse potential well approximates the single-atom potential along the nearestneighbor direction.
FIG.1.(Left) Due to the computational challenges of obtaining the multiphonon scattering rate in crystals, analytic approximations are valuable.Here we show a classification of regimes in which a multiphonon calculation has been performed, as well as approximations made in each case.In this work, we show that anharmonic corrections can be significant for q ≲ √ 2mN ω0 (Sec.III B) but are negligible when q ≫ √ 2mN ω0 (Sec.III C).We obtain results for all q using numerical calculations (Sec.IV A). (Right) To estimate anharmonic effects, we take a toy model of dark matter scattering off an atom in a 1D anharmonic potential.We obtain the anharmonicity by fitting to empirical models of interatomic potentials.

4 FIG. 6 .
FIG.6.Perturbativity bound on λ2  3 and λ4 as a function of phonon number n.The bound is based on the criteria of (49) that the leading correction to the energy En is at most 10%.The dashed line shows the typical coupling sizes in Si and Ge crystals.

2 ] 4 ω 2 ]
FIG. 7.q-dependence of structure factor: We the structure factor in the harmonic and anharmonic cases, where in the latter case the structure factor is calculated numerically with the maximal anharmonicity.The lines from top to bottom show the structure factor at different ω, corresponding to an increasing minimum phonon number n.There are large corrections for q ≪ √ 2m d ω0 when anharmonic interactions are included (dashed), and the corrections become more significant as the threshold is increased.For q ≫ √ 2m d ω0, both cases converge to the same result.For Si, we have √ 2m d ω0 ≈ 40 keV while for Ge, √ 2m d ω0 ≈ 50 keV.For other materials, this quantity is listed in TableI.The incoherent approximation momentum cutoff is qBZ < 2π/a ∼ 2.2 keV for both crystals.

4 ωth 4 ωthFIG. 9 .
FIG. 9.Ratio of anharmonic to harmonic rate.For each material (Ge and Si) we consider two representative values of the anharmonic couplings.The larger set corresponds to a direction of maximal anharmonicity while the other set corresponds to an orthogonal direction of intermediate anharmonicity.Anharmonic effects become more important for DM masses near the MeV scale and for larger energy thresholds.
FIG.10.Cross section uncertainty.Comparison of the cross section corresponding to 3 events/kg-yr in the harmonic (solid) and anharmonic (dot-dashed) cases.The anharmonic result is shown for maximal anharmonicity, and so the shaded band represents our estimate of the theoretical uncertainty due to anharmonic effects.The effects are primarily important for high thresholds and low DM masses, corresponding to large σn, which is generally in tension with existing astrophysical or terrestrial constraints.

2 ] 2 ω 4 ω
FIG. 11.Comparison of analytic structure factor in the Morse potential and the numerical calculation for Si as described in Sec.IV.We find that the two methods give almost the same result due to the fact that the Morse potential well approximates the single-atom potential along the nearestneighbor direction.