Multiphonon excitations from dark matter scattering in crystals

For direct detection of sub-MeV dark matter, a promising strategy is to search for individual phonon excitations in a crystal. We perform an analytic calculation of the rate for light dark matter (keV $<m_{DM}<$ MeV) to produce two acoustic phonons through scattering in cubic crystals such as GaAs, Ge, Si and diamond. The multiphonon rate is always smaller than the rate to produce a single optical phonon, whenever the latter is kinematically accessible. In Si and diamond there is a dark matter mass range for which multiphonon production can be the most promising process, depending on the experimental threshold.

The quest to directly detect dark matter (DM) in a laboratory experiment has, in recent years, significantly diversified [1]; both theoretical and experimental developments have driven the search beyond the WIMP paradigm. A steady decrease in energy thresholds has enabled sensitivity to particle-like dark matter with a mass well below that of a typical WIMP, with masses as low as an MeV currently being probed. Next generation detectors aim to push down to the lower limit of particle-like dark matter, probing the keV-MeV mass range. The energy scales for excitations created in current and proposed detectors coincide with the energy scales typical of many-body excitations in condensed matter or atomic systems. For DM heavier than 1 MeV, electronic excitations in atoms or semiconductors with energy gaps in the ∼ eV range are well suited, if the DM couples to electrons [2][3][4][5][6][7][8][9]. For light DM with nucleon couplings on the other hand, one must rely on chemical bond breaking [10], nuclear de-excitations [11], crystal defects [12,13] or soft nuclear recoils, where the latter in particular require very low thresholds [14][15][16][17][18][19][20]. For DM lighter than 1 MeV, vibrational modes in crystals [9,[21][22][23], molecular systems [24,25] or superfluid helium [26][27][28][29] naturally have energy spectra in the required 1 -100 meV range. Possible alternative detection strategies in this mass range are electronic systems with ultra-low bandgaps [30][31][32][33][34], magnon excitations [35] and avalanche gains in molecular magnets [36].
On the theoretical side, it is necessary to understand dark matter interactions with vibrational modes rather than with individual nuclei, so as to reliably estimate sensitivity and, in the event of any signal, extract dark matter properties. The reason is that for DM lighter than ∼ MeV, its de Broglie wavelength exceeds the interparticle spacing in typical materials, and it becomes necessary to transition to a different effective theory by integrating out the nuclei and electron clouds. One can therefore expect new and interesting features in these interactions, as they are subject to different kinematics and symmetry principles than those which govern the interactions in conventional dark matter experiments.
In this work, we focus on theoretical calculations for DM to excite vibrational modes (phonons) in crystals. In a crystal, phonons can be characterised as either acoustic or optical. The acoustic phonons are the Nambu-Goldstone modes associated with the breaking of translation symmetry by the crystal lattice; they must, at low energies, obey a linear dispersion relation. This feature in particular poses an experimental challenge for the detection of DM with mass below ∼100 keV: the momentum transfer in this regime is comparatively low, and the linear dispersion relation of the acoustic branch with typical slope ∼ 10 −5 implies a very low energy transfer (∼ meV). The optical phonons, on the other hand, are gapped and have 10 meV of energy at arbitrarily low momentum transfer. This makes them experimentally much more attractive.
There are also important theoretical differences between acoustic and optical phonons, as their couplings to DM depend very strongly on the DM model [21,22,37]. Concretely, if the DM has a coupling proportional to electric charge, there is a dipole interaction with the optical branches while the coupling to the acoustic branches is strongly suppressed. If instead the DM has a coupling p i k 2 k 1 p f < l a t e x i t s h a 1 _ b a s e 6 4 = " K F K P c P s Y v L D J z L c / j j p g C R 7 G J / U = " > A A A D Z n i c b Z L f i 9 N A E M c 3 V 3 + c U c + e I j 7 4 E m w P T g k l S U F 9 O T j 0 x c c T 7 N 1 B W 8 p m M 0 m W b n b j 7 q Z Y Q / 5 F 3 / 0 L 9 A / w V X E 3 j X r t u R C Y f O Y 7 8 x 2 G i U t G l Q 6 C r 8 5 e 7 8 b N W 7 f 3 7 7 h 3 7 9 0 / e N A / f H i u R C U J T I h g Q l 7 G W A G j H C a a a g a X p Q R c x A w u 4 u V b m 7 9 Y g V R U 8 A 9 6 X c K 8 w B m n K S V Y G 7 Q 4 d L J Z D B n l d V q k K W X Q 1 C m s e S Z x m T e u 9 y / X k h f N c R g E / s v g u T s z j E G q a z r 2 a e T T 0 K o N k z T L d S 3 G v o h 8 8 Q e a n r I w d k 1 N Q 3 8 V t o m W l z n m W h S N 7 Y K r T 7 4 Y 7 y Y s X X X y V A J 8 B v e o V S R Y 5 a B 8 h m N g J 8 O P Q 5 + N F E 3 g p B 2 g q Y 3 L K u q a b a R N b U a y 5 p 1 H W 1 k P y w U d 2 r l 2 Y G q g 2 I b L R W h h t A M j C 0 1 P d w Y 8 u b K q v / / t W t 1 F f x C M g v Z 5 1 4 O w C w a o e 2 e L / r d Z I k h V A N e E Y a W m Y V D q e Y 2 l p s Q 2 n F U K S k y W O I O p C T k u Q M 3 r 9 i I a 7 8 i Q x E u F N B / X X k u v V t S 4 U G p d x E Z Z Y J 2 r 3 Z y F / 8 t N K 5 2 + n t e U l 5 U G T j Z G a c U 8 L T x 7 X l 5 C J R D N 1 i b A R F I z q 0 d y L D H R 5 g i 3 X O y d F a X Z m 9 l M u L u H 6 8 F 5 N A r H o + h 9 N D h 9 0 + 1 o H z 1 F z 9 A x C t E r d I r e o T M 0 Q c T 5 4 v x w f j q / 9 r 7 3 D n q P e 0 8 2 0 j 2 n q 3 m E t l 7 P + w 0 k H h b n < / l a t e x i t > q p i k 2 k 1 p f < l a t e x i t s h a 1 _ b a s e 6 4 = " 4 n f M U g l v K j X a a L k + m Y b 5 U H x Z Z U Q = " > A A A D Z X i c b Z L P b t N A E M b X D X 9 K o J A C 4 s I B i w S p V F Z k O x J w q V T B h W O R S F s p i a L 1 e j Z e Z b 1 r d t c R Y e V H 5 A F 4 g j 4 A V 5 B Y O 6 E k K S N Z G v / m m 5 n V p 0 k K z r Q J w x / e X u v W 7 T t 3 9 + + 1 7 z 8 4 e P i o c / j 4 X M t S E R g S y a W 6 T L A G z g Q M D T M c L g s F O E 8 4 X C T z D 3 X 9 Y g F K M y k + m 2 U B k x z P B K O M Y O P Q 9 N C j 4 w R m T F i a U 8 o 4 V J b C U s w U L r K q 7 f + r N e S 4 O o r C M H g T v m 6 P H e N A j W W D g M U B i 2 q 1 Y 4 r N M m P l I J B x I P 9 C N 1 P l b l 1 l W R Q s o q b Q 8 C L D w s i 8 q q f g 8 m s g B 7 u F m i 7 W c q o A v s F K k G K d g Q 4 4 T o C f 9 L 7 0 A t 7 X L I W T Z n 9 l 3 Z J F X G 1 K K + t e t I i v V z S d t l d M W a 9 + 1 g 6 k D s p t O J 9 G N Y x 3 Y F x D N 7 M 9 B p F u O H X 9 3 7 j a n n a 6 Y T 9 s w r + Z R O u k i 9 Z x N u 1 c j V N J y h y E I R x r P Y r C w k w s V o a R e u C 4 1 F B g M s c z G L l U 4 B z 0 x D Y H U f m v H E l 9 K p X 7 h P E b u t l h c a 7 1 M k + c M s c m 0 7 u 1 G v 6 v N i o N f T e x T B S l A U F W i 2 j J f S P 9 + r r 8 l C k g h i 9 d g o l i 7 q 0 + y b D C x L g b 3 N p S n 1 l e O N + c M 9 G u D z e T 8 7 g f D f r x p 7 h 7 + n 7 t 0 T 5 6 j l 6 i I x S h t + g U f U R n a I i I 9 9 3 7 6 f 3 y f u 9 d t Q 5 a T 1 v P V t I 9 b 9 3 z B G 1 F 6 8 U f a g 4 W u Q = = < / l a t e x i t > proportional to the atomic mass, the coupling to the optical branches is strongly suppressed. At the same time, detecting single acoustic phonons is expected to be extremely challenging experimentally. This motivates the study of processes where multiple phonons are produced, which can have larger energy transfer. This was studied already in the context of superfluid helium [26][27][28][29], where the sound speed is particularly low and multiphonons were found to extend the reach for sub-MeV DM.
The purpose of this paper is to compute multiphonon processes for cubic crystals such as Ge, Si, GaAs and diamond in the isotropic approximation. Such materials are either already being used or considered for direct detection experiments, and it was found previously that the isotropic limit matches the numerical result well for single phonon excitations in GaAs [22]. More complicated, strongly anisotropic crystals, such as sapphire, are left for future work. We focus on DM that couples proportional to atomic mass number of the target nuclei, as it is in this scenario where multiphonon corrections are the most important. We focus on two acoustic phonons in the final state, for which there is a well-known effective theory, and briefly comment on multiphonon excitations with optical phonons.

B. Summary of results
The main object we are computing is the structure factor S(q, ω), which parametrizes the scattering rate of an external probe to the crystal for a momentum transfer q and energy transfer ω (see Sec. II for a precise definition). There are two distinct contributions to S(q, ω) from the production of two phonons, represented by the diagrams in Fig. 1. The left-hand diagram relies on a contact interaction between the DM and two phonons, which originates from the matching between the low energy effective phonon theory and UV theory of nuclei and electrons. There are analogous operators with three, four or more phonons, for which each additional phonon comes with a factor of q/ √ m N ω, with m N the nucleus mass. For m DM < MeV, q/ √ m N ω is a good expansion parameter, rendering the ≥ 3 phonon contributions negligible. For higher DM masses, the breakdown of this expansion signals the transition to the regular nuclear recoil regime. A resummation procedure is needed in this transition regime, which we do not attempt in this paper.
The right-hand diagram in Fig. 1 only relies on the DM coupling to a single phonon, and is therefore lower order in q/ √ m N ω. This process instead occurs via an off-shell phonon and the phonon self-interactions, which arise from the anharmonicity of the crystal potential. We will see in Sec. III that in the low-momentum regime, the self-interactions of acoustic phonons are governed by multiple dimensionful parameters that are related to the elastic constants of the crystal. The sense in which the self-interactions are "small" can be most easily seen from the fact that the typical width of the longitudinal acoustic phonon, Γ, is very small compared to its energy, 1 in other words Γ/ω 1. In the m DM 1 MeV regime, it is instructive to further expand S(q, ω) in the low q limit, as this allows for a qualitative comparison between different channels and materials. The resulting scaling is represented schematically in Table I. For a single acoustic mode in the final state, S(q, ω) scales linearly with q and is by far the most favorable in terms of rate, but requires a very low threshold. For both the single optical mode 2 and the two-phonon contact interaction, S(q, ω) scales as ∼ q 4 , while for the anharmonic contribution it scales more favorably as ∼ q 2 . However, the latter requires an insertion of the phonon self-interaction, which also provides a suppression. Quantitatively, we find that the rates of both two-phonon contributions are smaller than the rate for the single optical mode, whenever the latter is kinematically accessible. Finally, it is interesting to compare crystals with superfluid helium, where the rate scales as ∼ q 4 . (See Sec. V B.) The phonon self-couplings in helium are however much stronger than in most crystals, and the reach ends up being competitive with or better than the cubic crystals we considered. Our final, quantitative results are shown in Fig. 4.
The paper is structured in the following way. In Section II, we first give general expressions for the crystal structure factors that determine the rate into one and two phonons (which could be evaluated with numerical phonon eigenmodes and couplings in full generality, without the following approximations); we further present formulae for the case of scattering in the long wavelength limit, relevant for light dark matter. In Section III we introduce the isotropic approximation, and detail the elasticity theory used to determine the necessary sound speeds and anharmonic parameters for channel low-q scaling typical threshold needed Ref.
TABLE I. Leading scaling of the structure factor S(q, ω) in the low q (low m DM ) limit for different channels, and required approximate thresholds to observe them. It is assumed that the DM couples proportional to the mass of the atoms. The # indicates that this channel vanishes in the limit where the (material dependent) phonon self-couplings are taken to zero. 1 Note this is different from superfluid He, where the phonon-roton self-interactions are much larger, but where the phonon decay is kinematically forbidden for part of the dispersion curve. 2 The single optical mode scales as q 4 for dark matter that couples proportional to mass [37], which is the situation considered here; otherwise, it scales as q 2 . the acoustic phonons. In Section IV we present sensitivity curves for various crystals and compare rates for single and two acoustic phonon final states. We provide an estimate of rates for two-phonon final states with one or two optical phonons in Section V, finding that they are subdominant to the single optical channel. We also briefly comment on the qualitative differences with superfluid helium. Section VI concludes with a short discussion. We include three appendices which detail more lengthy aspects of the calculations, where we aim for our results to be self-contained within this paper, and reproducible.

II. SCATTERING FORMALISM
The scattering rate for an incident DM particle to excite phonons in a crystal is given by the dynamic structure factor, or simply structure factor. In this section, we establish our notation and provide a derivation of the structure factors for single-and two-phonon excitations. In both cases, we obtain approximate formulae for the structure factors when the final states consist of long-wavelength acoustic phonons. In this limit, the acoustic modes have nearly linear dispersion and the structure factor can be expressed in terms of bulk properties such as sound speeds, target density, and so on.
We begin with the most general form of the potential seen by an incident DM particle of mass m DM : where the index J sums over all scattering centers (ions) with crystal position coordinate r J , b J is a factor that depends on the DM coupling with atom J, and the tilde indicates the Fouriertransformed function. Note that we assume from the start a crystal lattice containing N primitive unit cells and n ions per unit cell. Boldface symbols indicate 3-vectors in position or momentum space, while non-boldface symbols indicate scalar quantities (e.g. q ≡ |q|).
Two specific cases of Eq. (1) are of particular interest: a contact interaction between DM and nuclei for whichF (q) = 1, and scattering via a massless mediator withF (q) ∝ 1/q 2 . The DM wavelength is always much larger than the radii of the nuclei, so we set the nuclear form factors to 1 everywhere. We also assume a coupling proportional to atomic mass number, A J . We then have b J = 2πb n A J /m DM , where b n is the DM-nucleon scattering length and σ n ≡ 4πb 2 n is the DM-nucleon scattering cross section.
Pulling out the overall factor of 2πb nF (q)/m DM , we focus on characterizing the expectation value of the sum over scattering centers and define a dynamical structure factor given by where ω and q are respectively the energy and momentum transferred from the DM to the crystal. Φ f | represents the collection of final states, indexed by f and having energy E f . We have assumed that the system is in its ground state |0 before the collision; this is an excellent approximation since any dark matter experiment would necessarily be operating at very low temperatures, with negligibly small numbers of thermal phonons present. Each term then represents a scattering probability to excite a given final state. The differential cross section is moreover closely related to the structure factor; for example, taking the isotropic limit for a material, where v is the initial DM speed in the lab frame.
To evaluate (2) for final states with a specific number of phonons, we must expand the position vectors r J in terms of equilibrium positions and displacement vectors. For a crystal with repeating primitive cells, the sum over atoms J can be broken up into a sum over the lattice vectors for the primitive cells, indexed by , and the atoms in the primitive cell, indexed by d. The position operator can then be written as r J = + r 0 d + u ,d , where r 0 d is the equilibrium location of atom d relative to the origin of the primitive cell, and u ,d is the displacement of that atom relative to its equilibrium position.
We quantize the displacement vector in the harmonic approximation, following the convention in Ref. [22], here adapted to the Schrödinger picture operator: where there are 3n phonon branches, indexed by ν, for a primitive cell containing n atoms. Here m d is the mass of atom d,â † ν,k andâ ν,k are the creation and annihilation operators for the phonons, ω ν,k is the energy of phonon branch ν at momentum k, and e ν,d,k is the phonon eigenvector (normalized within a unit cell) for atom d.
Using (4), the structure factor can then be expressed as where W d (0) is the zero-temperature Debye-Waller factor for atom d, and M f,q,d is the matrix element associated with final state f , This expression represents the matrix element for scattering into the crystal final state labeled by f , at leading order in V (q); however, it is not yet practical for concrete calculations. As explained in the introduction, q/ √ m d ω ν,k < 1 for the DM mass range of interest, which means we can consistently expand the exponential factor. This amounts to an expansion in the number of phonons coupling to the DM, where the quadratic (two phonon) contribution is represented by the left-hand diagram in Fig. 1. Once crystal anharmonicity is included, we also expand in the phonon self-interactions. The leading contribution in terms of the phonon self-couplings is shown in the right-hand diagram in Fig. 1. In summary, the calculation amounts to a double expansion in the momentum transfer q and the phonon self-couplings.

A. Single-phonon structure factor
For the final state consisting of a single phonon with polarization ν and momentum k, the leading result for the matrix element is where G are the reciprocal lattice vectors, which satisfy e i ·(q−k) = N G δ G,q−k , with the Kronecker-δ enforcing momentum conservation in the crystal. Here we also used that phonon observables such as ω ν,q are invariant under q → q+G. While there can be anharmonic corrections to the above matrix element, they are negligible in the low q limit.
Summing over all possible single-phonon final states, this gives a structure factor identical to the result in Ref. [22]. For sub-MeV DM scattering, where q is typically well within the first Brillouin zone, it is a good approximation to neglect the sum over G as well as the Debye-Waller factors. Then the result simplifies to In the long wavelength (low q) limit, we can moreover approximate the acoustic modes as having real eigenvectors with magnitudes given by |e ν,d,q | ≈ m d /( n d m d ) and with polarization vector independent of d. It is therefore convenient to introduce "long-wavelength polarization vectors" with unit length by defining the (real) vector The difference between the two objects should be clear from the presence/absence of the index d labelling atoms in the primitive unit cell. We can then simplify sums over atoms, for example by In the last step, we have made the approximation that the bound atom masses are given by m d ≈ A d m p , with m p the proton mass. As can be seen from (7), transverse polarizations cannot contribute to the single phonon rate. Considering the one-phonon structure factor for longitudinal acoustic (LA) phonons, we can take e ν,q = q/|q| with the result The LA dispersion in the long-wavelength limit is linear with slope given by the sound speed associated with the LA mode, c LA (q), which in general can depend on the phonon propagation direction. Note that the factor of d A d will drop out in the expression of the rate per unit target mass, so that the rate to excite a single acoustic phonon depends only on the sound speed. 3

B. Two-phonon structure factor
For the case with two phonons in the final state, there are two pieces which contribute to the matrix element: a contact term from expanding the exponential in (6) to second order, and a piece resulting from anharmonic phonon interactions in the material. We define δH as the leading order anharmonic phonon interaction Hamiltonian; its precise definition we defer to Sec. III A. At leading order, the 2-phonon matrix element is then with where the factor s 1,2 ≡ (δ ν 1 ,ν 2 δ k 1 ,k 2 + 1) −1/2 accounts for Bose statistics. The contributions in (12) and (13) were shown diagrammatically in Fig. 1, and we refer to them as the contact term and the anharmonic term, respectively. Anharmonic phonon interactions also lead to a non-zero phonon width, Γ ν,k . This has been resummed in the phonon propagator in (13) and becomes relevant when the intermediate phonon goes on-shell. Details regarding the derivation of the above matrix elements are given in Appendix A.
In the long wavelength limit, we can again consider only the G = 0 contribution to the matrix elements and drop the Debye-Waller factors. It will then be convenient to express the three-phonon matrix element as where V is the volume of the crystal. As we show in Sec. III A, in the long wavelength limit M(q, k i , ν i ) is a function only of the momenta, long-wavelength polarization tensors and elastic constants of the material. In addition, eigenvectors are real in this limit, such that the matrix element M(q, k i , ν i ) is real as well. The two terms in (11) therefore do not interfere to leading order in the small q expansion, when neglecting terms higher order in Γ ν,k . Using the longwavelength polarization vectors defined in (9), the two-phonon structure factor can be simplified where we took the continuum limit by substituting is the mass density of the material. Similar to the single-phonon structure factor, the overall factor of d A d will drop out in the expression of the rate per unit target mass, so that the rate to excite two phonons depends only on bulk properties such as sound speeds, density, and elastic constants.

III. EVALUATION OF STRUCTURE FACTORS
In this section we provide explicit results and analytic formulae for the two contributions to the two-phonon structure factor, Eqs. (16)- (17). Even in the long-wavelength limit, the dispersions and anharmonic couplings are in general direction-dependent, substantially complicating the calculations. For cubic crystals, the isotropic limit is however known to be in excellent agreement with the general result for scattering to single phonons [22]. In this work we will therefore restrict ourselves to cubic crystals such as GaAs, Ge, Si and diamond, and approximate them as isotropic. We leave a fully general calculation of the multiphonon rate with Density Functional Theory (DFT) for future work, but we do not expect that accounting for anisotropy would qualitatively change our conclusions.
In the isotropic limit, both transverse acoustic polarizations are degenerate and the dispersion relations are simply ω LA,q = c LA q and ω T A,q = c T A q, with c LA and c T A the average sound speeds associated with the longitudinal acoustic (LA) and transverse acoustic (TA) modes, respectively. The structure factor for scattering to a single acoustic phonon then simplifies to For the multiphonon contribution, a description of the phonon self-interactions is needed, and this is where the isotropic approximation is most advantageous: as we will see in Sec. III A, the effective Hamiltonian is relatively simple in the isotropic and long-wavelength limit, containing 5 independent operators (this number grows to 9 if instead cubic symmetry is assumed). The coefficients of these operators can moreover be extracted from the elastic properties of the material. Each coefficient maps directly to a linear combination of the second order elastic constants (related to the bulk modulus and Young's modulus) and third order elastic constants; these quantities can either be measured or computed with ab initio methods.

A. Anharmonic term
To compute the anharmonic contribution, we use a low-momentum effective description of the phonon self-interactions. As for any effective theory, we first constrain the form of the Hamiltonian using the symmetries of the theory and subsequently fix the Wilson coefficients from measured observables, or by matching on to the full UV theory. It is hereby convenient to introduce a "longwavelength displacement operator", in analogy to the long wavelength polarization tensors defined in (9). Replacing the polarization tensors with their long-wavelength versions and averaging over the atoms in a unit cell, we can define where now we only sum over acoustic polarizations ν and we have replaced the individual atomic position vectors + r 0 d with the continuous position vector r. Once again, the long-wavelength displacement operators u can be distinguished from their more general counterparts u ,d by the index labels.
Assuming isotropy, there are only 5 independent operators to third order in the effective Hamiltonian [38,39]: with u ij ≡ ∂ i u j and the i, j running over the three spatial coordinates. Repeated indices are summed over. The coefficients α, β, γ, λ and µ can be determined from the measured or calculated elastic constants of the crystal. In particular, the parameters µ and λ are the Lamé parameters of the crystal and related to the bulk and Young's moduli. The parameters α, β and γ can be calculated from the third order elastic constants, as described in Appendix B. All five parameters have units of pressure and are reported in units of Giga-Pascal (GPa) in Tab. II for the crystals we consider.
Using (19)-(20), the anharmonic three-phonon matrix element can be written in the form of  (14), where the function M is given by: and we introduced the shorthand notation e = e ν,q , e 1 = e ν 1 ,k 1 etc. From (13) it follows that only the longitudinal polarization of the off-shell, intermediate phonon contributes. Depending on the polarizations of the outgoing phonons, different terms in (21) contribute. Concretely, there are four distinct combinations for which the matrix element is non-zero: • LA-LA • TA-TA with both phonons polarized in the plane spanned by the momenta • TA-TA with both phonons polarized orthogonal to the plane spanned by the momenta • LA-TA with the TA phonon polarized in the plane spanned by the momenta.
In the isotropic limit, the structure factor in (17) reduces to The anharmonic matrix element given in (21) can also be used to compute Γ LA,q , which we provide explicitly in Appendix C 1.

B. Contact term
With the definition of the long-wavelength polarization tensors in (9), the structure factor for the contact term in (16) reduces to which can also be evaluated analytically. Concretely, there are three final-state polarization configurations (LA-LA, TA-TA and LA-TA) which can contribute, where the TA modes must be polarized in the plane spanned by the momenta: with g (cont) where we again used δ = c LA c T A . Note that the expansion in (37) assumes x 1 δ . The exact expressions for g (cont) T AT A and g (cont) LAT A were used for all our numerical results (see Appendix C 2). Note that the O(q 4 ) scaling discussed in the Introduction is manifest in (32), (33) and (34).

C. Numerical comparison
The left-hand panel of Fig. 2 shows the different contributions to S(q, ω) for the example of GaAs, where we summed the different TATA contributions. We show the full kinematic range, where the most striking feature is the resonance at x = 1/δ for the anharmonic contributions, indicating that the intermediate LA phonon goes on-shell. Whenever the resonance is kinematically accessible, it dominates the rate to the extent that the off-shell multiphonon contribution is completely negligible. The LALA channel also cuts off for x > 1/δ, since in this regime it is not possible to simultaneously conserve energy and momentum. Except for the region near the resonance, all contributions scale as ω 4 with respect to our 10 meV reference value. The inset zooms in on the low momentum region and shows the ∼ q 2 and ∼ q 4 scaling of the anharmonic and contact contributions, respectively.
The long wavelength approximation necessarily breaks down at momenta approaching the edge of the Brillouin zone for two reasons: the dispersion relations of the acoustic phonons cease to be linear, and the description of the phonon self-couplings in terms of the elasticity parameters (Sec. III A) starts to break down. We show the dispersions in the right-hand panel of Fig. 2 for the example of GaAs, where the full dispersion relations [40] are compared with those in the long wavelength, isotropic approximation. To ensure that the calculation is not extrapolated beyond its regime of validity, we impose a maximum momentum cutoff of q cut = 0.7 keV for GaAs and Ge, q cut = 0.8 keV for Si, and q cut = 1.2 keV for diamond. This corresponds roughly to q cut ≈ q BZ /3 where q BZ ≡ 2π/a is the approximate boundary of the first Brillouin zone and a is the lattice spacing. The cut is indicated by the light gray shading in Fig. 2, and below this value the dispersions of the acoustic phonons in all four materials is close to linear. We also enforce this  Due to the relatively low sound speeds in GaAs and Ge, the phase space is substantially restricted by these consistency conditions. In this sense, our calculations should be viewed as a conservative estimate. The choice of q cut is to some degree arbitrary, and therefore we also compute all rates with a q cut that is twice the values reported in Tab. III. This provides a measure of the sensitivity of our results to q cut . We expect that the long wavelength formulas overestimate the structure factor when extrapolated beyond their regime of validity because of the strong growth in q and ω, and because the isotropic linear dispersions overestimate the mode energies at large momenta. In this sense we anticipate that the true answer is bracketed by the two cutoff choices. Numerically, we find that integrating fully out to the edge of the Brillouin zone does not change the rates appreciably in comparison to our upper choice of ∼ 2q BZ /3.

IV. RESULTS
Folding in the DM velocity distribution, the total rate per unit exposure is given by whereF (q) indicates a form factor whose functional form is determined by the properties of the particle mediating the DM-nucleon scattering process. The most common, limiting cases areF (q) = 1 if the mediator is heavier than the DM, andF (q) = v 2 0 m 2 DM /q 2 for a mediator which can be treated as massless in the scattering process. In addition, ω − is the energy threshold of the experiment, and where p i ≡ m DM v i and p f ≡ m DM v 2 i − 2ω/m DM are the magnitudes of the initial and final DM momenta respectively. The cuts involving q cut and ω cut ensure that the integral is not evaluated GaAs, m DM = 10 keV, σ n = 10 −30 cm 2 in a regime where the long wavelength approximation is invalid, as discussed in Sec. III C. For the DM velocity distribution f (v i ) we use the standard truncated Maxwellian distribution in the Earth's frame: and we take v 0 = 220 km/s, v esc = 500 km/s, and the Earth's average velocity to be v e = 240 km/s. Fig. 3 shows the differential scattering rate as a function of the deposited energy, assuming a massless mediator. All curves are cut off when the momenta of the final state phonons are outside the first Brillouin zone. The dotted vertical lines indicate values of ω cut , above which we expect that the long wavelength approximation starts to break down. Integrating the rate beyond ω cut to the edge of the Brillouin zone is likely to overshoot the true answer. For the m DM = 10 keV benchmark (left-hand panel), the single phonon resonance occurs for ω < 1 meV, while its enormous contribution to the scattering rate is visible for ω < 5 meV for the 50 keV benchmark (right-hand panel). The anharmonic terms always dominate over the contact terms. Fig. 4 shows the cross sections needed to obtain 3 events with a kg-year exposure, again assuming a massless mediator. 4 The most striking feature in Fig. 4 is the enormous enhancement of the reach once the single acoustic phonon becomes accessible. In this regime, integrating the multiphonon structure factor matches onto the single phonon scattering rate (see Appendix A) and we can simply use the single phonon structure factor. For a given experimental threshold ω − , the mass m * DM at which the single-phonon resonance appears may be analytically derived by requiring that the maximum momentum transfer supplied by the DM suffices to create an on-shell LA phonon above the threshold, or in other words: 2m * Using the sound speed for GaAs and a 1 meV threshold as an example, the single-phonon resonance will appear at m * DM ≈ 12 keV, as can be seen in Fig. 4. The very high sound speed of diamond then explains why this material maintains sensitivity to the single acoustic mode for most of the mass range, even for a threshold as high as ∼ 10 meV. (See [20] for a detailed study of diamond as a dark matter detector.) No backgrounds or experimental efficiencies have been included in Fig. 4, which is meant to both illustrate the most optimistic reach possible, as well as the relative importance of the various channels, rather than provide an accurate projection of the absolute reach. The single optical phonon channel is computed using an analytic approximation given in Sec. V A, with the (dispersionless) optical mode energy given in the figure labels. We see that the multiphonon channel is always subleading to the single optical, except at low m DM for Si and diamond. The reason is that the longitudinal optical mode in both these materials is relatively high energy, 60 meV and 140 meV respectively, and is not kinematically accessible in the low m DM region. For comparison, multiphonon production in superfluid helium is also shown in Fig. 4; in an idealized setting where all experimental effects aside from the threshold are neglected, it always outperforms both the single optical and multiphonon modes in crystals.
The shaded bands in Fig. 4 indicate the estimated uncertainty from taking the long wavelength approximation, by displaying the calculated rates using two choices for the momentum cutoff, as explained in Sec. III C. Concretely, the upper edge of the band corresponds to the values reported in Tab. III, whereas the lower curves assume twice these values. This source of uncertainty is negligible once the single LA channel is accessible, as this contribution is peaked at low ω (see right-hand panel of Fig. 3), and moreover it does not rely on the validity of Eq. (20). The size of the band is larger in GaAs and Ge because of the lower sound speeds and ω cut (Tab. III). This source of uncertainty is also more severe as the experimental threshold is increased, since this reduces the available phase space in Fig. 3, which leads to greater dependence on ω cut . For a 10 meV threshold, the lower value of ω cut severely restricts the phase space for the TATA channel, especially for GaAs and Ge. Meanwhile, for diamond ω cut has no effect on the rate, since it is always larger than the initial DM kinetic energy when m DM < m * DM . We therefore expect the long wavelength limit be an excellent approximation in this case.
Other sources of uncertainty are the values for the elasticity parameters, as to the best of our knowledge they have not yet all been measured at ultra low temperatures. As explained in Appendix B, we instead rely on ab initio calculations of these parameters, which in some cases carry O(1) uncertainties. This propagates to an O(1) uncertainty on the overall multiphonon rate, regardless of the DM mass. In addition, we expect corrections to the isotropic approximation once the detailed crystal structure is accounted for. These uncertainties are not included in the band in Fig. 4. Given the current experimental unknowns, we consider the uncertainties acceptable at this stage, especially given that multiphonon processes typically have a much lower rate than the single optical mode.
To conclude, we briefly comment on stellar cooling constraints, warm dark matter bounds and the material overburden. For millicharged particles with mass 10 keV, there are strong constraints from the cooling of white dwarfs, red giants and horizontal branch stars [43,44]. To our knowledge, the analogous computation has not yet been performed for light DM with a coupling to nuclei, but we expect that similar constraints should apply for m DM 10 keV. In this mass range, the DM is also generally considered as warm and there are constraints from structure formation, although these are alleviated if this candidate doesn't provide the entire DM abundance. The likely existence of both bounds is suggested by the gray shading in Fig. 4. Finally, for sufficiently large σ n , the DM is likely to scatter in the Earth's crust before reaching an underground detector.
To determine roughly where this occurs, we estimated the rate for DM to lose at least 1% of its initial kinetic energy by scattering in a silicon crust. The dotted line in Fig. 4 indicates where the corresponding mean free path is 1 km.

A. Multiphonons involving optical branches
As discussed in the introduction, the rate for scattering that excites a single optical phonon is suppressed when the DM coupling is proportional to the mass of the atom. Nevertheless, as seen in the previous section, processes involving optical phonons are still important, particularly for higher experimental thresholds. In this section, we briefly review the single longitudinal optical (LO) phonon calculation, before discussing two-phonon processes involving optical phonons.
To obtain an estimate of the rate to excite a single LO phonon, we use an approximation for the eigenmode in a cubic lattice with diamond or zincblende structure, valid at low q: where r 2 = (a/4, a/4, a/4) is the position of the second atom in the primitive cell and a is the lattice constant. Note that without the phase factor the structure factor would be exactly zero; we have included it to account for the subleading behavior [37]. Using eq. (8) and averaging over angles such that (q · r 2 ) 2 ≈ q 2 a 2 /16, we obtain where we have approximated the LO phonon's dispersion relation as flat, ω LO (q) = ω LO . This approximation reproduces the full numerical result for the DM reach in GaAs (see Ref. [22]) to within an O(1) factor.
There are two kinds of two-phonon processes involving optical phonons to consider: opticalacoustic, and optical-optical. We begin with the former, since they are the most relevant for light DM. Optical-acoustic scattering also has both contact and anharmonic contributions. For all of the materials we consider, there is a suppression of the contact contribution at low q when DM couples proportional to atomic mass. This can be seen from the expressions for the structure factor and matrix element in eqs. (5) and (12). When q = 0, momentum conservation requires k 1 = −k 2 and the sum over the unit cell in (5) vanishes due to the orthogonality of the eigenvectors. Using the low-q approximation for the LO eigenvector (43), one can explicitly see that the leading term in the small q expansion of the structure factor vanishes; the contact term then scales as q 6 and is negligibly small. Note that this result does not hold for general lattices, since with more complicated unit cells there can be mixed longitudinal-transverse optical modes which may only be orthogonal to the acoustic modes after also contracting the Lorentz indices of the eigenvectors.
The anharmonic contribution is more difficult to reliably calculate. It could be obtained from a first principles calculation of the anharmonic corrections to the lattice potential using Density Functional Theory, however this goes beyond the scope of the present paper. Here, we adopt a simpler method in order to obtain an estimate of the size of this contribution. We follow an approach that has been used in the literature to calculate the lifetime of LO phonons and describe the anharmonic three-phonon interactions via the Hamiltonian [45], where γ G ≈ 1 is the mode-averaged Grüneisen constant andc is the average of the LA and TA sound speeds. The above Hamiltonian can be obtained starting from eqs. (19)- (21) and then averaging over phonon modes and angles (see Ref. [45]). Since this model treats the lattice as an isotropic continuum it does not actually contain optical modes; nevertheless, eq. (45) has been used in the calculation of optical phonon lifetimes (e.g. [46,47]).
The dominant anharmonic contribution is that mediated by an off-shell LA phonon, since the LO mediated process has the same suppression as single optical scattering. Using the Hamiltonian (45) in eq. (13) we obtain the structure factors, where we have again assumed a flat dispersion relation for the optical mode. The expressions for the TO-LA and TO-TA processes can be obtained by the substitution ω LO → ω T O and multiplying by a factor of two. Integrating the structure factor to obtain the total rate we find that, for all the materials we consider, the LO-LA scattering rate is one-two orders of magnitude smaller than the single optical rate, depending on the q cut used for the acoustic phonons. The LO-LA process becomes further suppressed relative to the single optical with increasing DM mass. The LO-TA process is enhanced by the smaller TA sound speed, but is still suppressed compared to the single optical. A similar conclusion holds for optical-acoustic scattering involving TO phonons, although these processes could be relevant in a narrow range of DM masses that are above the threshold to excite a TO phonon but below the LO threshold. While eqs. (46) and (47) should only be considered as an estimate of the two phonon optical-acoustic rate, we do not expect a detailed DFT calculation to change the qualitative conclusion that it is sub-leading compared to single optical scattering. Next, we briefly discuss scattering into two optical phonons. This process only becomes kinematically accessible for heavier DM masses due to the higher energy threshold to excite two optical phonons. Unlike optical-acoustic scattering, there is no additional suppression of the contact contribution for DM that couples proportional to atomic mass. The LO-LO structure factor is then proportional to q 4 /(m p ω LO ) 2 . On the other hand, the single optical structure factor scales as q 4 a 2 µ/(m 2 p ω LO ), where µ is the reduced mass of the primitive cell. The two optical phonon contact contribution is then expected to be significantly smaller than the single optical. The anharmonic contribution is again challenging to reliably estimate; however, based on our above estimate for optical-acoustic scattering, where it was found to be sub-leading, we do not expect it to give a significant contribution. In summary, two phonon scattering processes involving optical phonons are expected to give only a sub-leading contribution to the total scattering rate.

B. Multiphonons in superfluid helium
Given that our results qualitatively differ from similar calculations in superfluid helium, it is worthwhile to compare the symmetries of both systems in a bit more detail. Crystals spontaneously break both translation and rotation invariance, but since the rotation operators are linearly dependent on the translation operators, there are only 3, rather than 6, Goldstone modes [48]. These are the 1 LA and 2 TA modes we have encountered throughout our discussion. Since translations are broken spontaneously, all amplitudes must vanish in the limit where one of the external (spatial) momenta go to zero. This symmetry principle explains the form of the amplitude in (21) and its scaling in the low q limit. Schematically, the matrix element for 3 acoustic phonons in a crystal is therefore always of the form Superfluid helium on the other hand does not break translation and rotation invariance, though the Bose-Einstein condensate breaks boost invariance as well as a linear combination of the time translation and particle number operators. All four broken operators are linearly dependent, such that there only exists a single Goldstone mode [48], which is the phonon-roton branch. In more mundane terms: phonons in superfluid helium do not have polarization tensors, which further restricts the form of the amplitude compared to (48). The Ward identity associated with the U (1) particle number symmetry moreover enforces that the amplitude vanishes in the q → 0 limit [49]. Bose symmetry on the final state momenta then implies that in the low q limit, the amplitude must be proportional to where the second ∼ follows from momentum conservation (q = k 1 +k 2 ). The key point here is that the O(q 2 ) contribution is absent, unlike in (48) for crystals. This behavior is a direct consequence of the broken U (1) symmetry and therefore a feature of the superfluid nature of the system. In concrete calculations it manifests itself as a delicate cancellation between seemingly unrelated terms in the amplitude. This cancellation has been observed both in a modern effective field theory treatment [28,29], as well as in an older, quantum hydrodynamic treatment (see Appendix A of [27]). In the latter case the cancellation only occurs if the interacting ground state of the theory is accounted for. Despite this suppression factor, the multiphonon rate in helium exceeds that in the crystals we considered (see Fig. 4), due to the stronger phonon self-couplings in helium.

VI. CONCLUSIONS AND OUTLOOK
In this work, we evaluated the rate for production of two acoustic phonons in crystals from scattering of sub-MeV DM. We considered cubic crystals such as GaAs, Ge, Si and diamond and worked in the isotropic and long wavelength approximations. In addition, we focused on DM which couples proportional to atomic mass, since in this case the rate for single optical phonon excitations is suppressed and multiphonon production is most relevant. However, for all four crystals, we found that the multiphonon rate is smaller than the single optical phonon rate whenever the optical mode is kinematically accessible. Similarly, the rate to excite a single acoustic phonon dominates whenever that mode is kinematically accessible. In diamond and Si there is, however, a range of DM masses between 10 keV and 100 keV for which the multiphonon process could be the only detectable channel, depending on the experimental threshold. We have also estimated the multiphonon rate with optical phonons and expect it to be sub-leading. In idealized experimental conditions, the multiphonon rate in superfluid helium exceeds that in all the crystals we have considered, despite the less favorable scaling of the structure factor in a low momentum expansion.
For GaAs and Ge, our approach here in taking the long wavelength approximation has a limited regime of validity, leading to appreciable uncertainties in the scattering rate. A more precise evaluation with Density Functional Theory methods would be desirable for these materials. Such a DFT treatment would also allow one to study anisotropic materials such as sapphire, which are expected to exhibit a sizable daily modulation in the multiphonon signal. where M (1−ph) is defined in Eq. (7), and V is a volume factor from the normalisation of the DM momentum eigenstates. The transition rate is directly related to the structure factor in Eq. (5) up to overall factors.
Note that for unstable final states (Γ ν,k = 0) w i→f vanishes. In this case the transition probability in Eq. (A3) does not grow with time (it is constant when → 0), since due to the exponential decay of the state only the last ∆t ∼ Γ −1 ν,k contributes significantly.

Two phonon
Next, consider scattering into the two phonon state |p f ; ν 1 , k 1 ; ν 2 , k 2 (with Γ ν 1 ,k 1 = Γ ν 2 ,k 2 = 0). In this case anharmonic effects enter at second (mixed) order in perturbation theory and can have a significant impact on the scattering rate. The transition rate is (A7) The contact and anharmonic contributions are shown diagrammatically in Fig. 1, with the matrix elements given in eqs. (12) & (13), and the δH matrix element discussed in Sec. III A. In the narrow width limit (Γ ν,k /ω ν,k → 0), and neglecting the interference terms, the anharmonic contribution reduces to the single phonon rate times the branching ratio to |ν 1 , k 1 ; ν 2 , k 2 . Similarly, while eq. (A6) is strictly only valid for scattering into stable final states, the narrow width approximation applied to multiphonon scattering justifies its use for final states with non-zero width.
which measures how a material responds under stress. Since x i = u i + a i by definition, we can use to rewrite the strain tensor as with u ij ≡ ∂ i u j . Note that η ij is manifestly symmetric.
The generalization of Hooke's law is [51] σ with C ijk the elastic constants and σ ij the stress tensor. This relation can be written in Hamiltonian form where the stress tensor σ ij acts as a source for the η ij . (B6) is then just the equation of motion of η ij given by this Hamiltonian. Dropping the source term, the Hamiltonian in (B7) can be further generalized to include the cubic response where the C ijk mn are the third order elasticity constants. C ijk is invariant under i ↔ j, k ↔ and (ij) ↔ (k ), C ijk mn is invariant under i ↔ j, k ↔ , m ↔ n and the permutations of the (ij), (k ) and (mn) pairs. In the most general case, C ijk and C ijk mn have therefore respectively 21 and 56 independent components.

The isotropic approximation
The cubic crystals we consider in this work are not completely isotropic but instead are only invariant under permutations of the x, y and z axes and parity transformations such as x → −x etc. The latter imply that all components of C (cub) ijk and C (cub) ijk mn for which a value of an index occurs an odd number of times must vanish (e.g. C (cub) 1222 = 0 etc). One can show that imposing these symmetries reduces the general elasticity tensors to 3 independent second order elastic constants, and 6 independent third order elastic constants. In order to express the 5 isotropic elasticity parameters µ, λ, α, β and γ in terms of these 9 measured elasticity parameters for the cubic crystals of interest an averaging procedure is needed.
Given that a 6-tensor such C (cub) ijk mn can be rather unwieldy, much of the literature has chosen to adhere to the Voigt convention, where each pair of double indices (ij), (k ) and (mn) is replaced with a single index running from 1 to 6 through the mapping This maps C ijk ) respectively, where we have used lowercase c for components of the elasticity tensors in Voigt notation. The independent elasticity parameters for a cubic crystal, as typically reported in the literature, are c 11 , c 12 and c 44 for the second order elastic tensor and c 111 , c 112 , c 123 , c 144 , c 166 and c 456 for the third order elastic tensor 5 , where we have dropped the (cub) superscript going forward. All other components either vanish or can be obtained by applying one of the symmetries listed above. An explicit representation of c ij and c ijk can be found in e.g. [52].