Directional Detection of Light Dark Matter with Polar Materials

We consider the direct detection of dark matter (DM) with polar materials, where single production of optical or acoustic phonons gives excellent reach to scattering of sub-MeV DM for both scalar and vector mediators. Using Density Functional Theory (DFT), we calculate the material-specific matrix elements, focusing on GaAs and sapphire, and show that DM scattering in an anisotropic crystal such as sapphire features a strong directional dependence. For example, for a DM candidate with mass 40 keV and relic abundance set by freeze-in, the daily modulation in the interaction rate can be established at 90\% C.L. with a gram-year of exposure. Non-thermal dark photon DM in the meV - eV mass range can also be effectively absorbed in polar materials.

Sub-GeV dark matter (DM) has become an important direction in dark matter searches in recent years. While low mass DM models have long been recognized to be viable theoretically, only recently have they come within experimental reach for direct searches. The correct dark matter abundance can be obtained in a variety of models, such as Hidden Valleys [1], secluded DM [2,3], Asymmetric Dark Matter [4][5][6], freeze-in dark matter [7], supersymmetric hidden sectors [8][9][10], and strongly interacting massive particles [11]. All of these models can contain light mediators which couple the DM to the Standard Model, and in some cases such a mediator is crucial to set the relic abundance. The presence of light, weakly coupled mediators has made detection in high energy collisions challenging, while the sensitivity of direct detection experiments has been limited by the small kinetic energy of light DM in the Milky Way. Recent advances in low threshold techniques have put direct detection within reach, however.
There is now a wide range of upcoming experimental probes for DM with mass between 1 MeV and 1 GeV (see Ref. [12] for a comprehensive overview). For DM-electron couplings, currently the best sensitivity in this mass range is achieved by semiconductor target experiments such as Super-CDMS [13], SENSEI [14] and DAMIC [15], as well as noble liquid experiments such as DarkSide [16] and XENON10/100 [17]. A graphene target in a prototype of the PTOLEMY experiment also has directional sensitivity to sub-GeV DM with electron couplings [18]. For DM-nucleon couplings, the CRESST collaboration [19] has obtained sensitivity to masses as low as 0.5 GeV. Other experiments targeting nuclear recoils in this mass range include DAMIC [15], NEWS-G [20] and SuperCDMS SNO-LAB [21], while there are proposals to use liquid helium [22], molecules [23], or crystal defects [24,25] as detection targets.
Extending the reach to keV-MeV DM particles presents greater challenges, but also new opportunities. The main challenge is to detect the extremely small energy depositions in ultra pure targets.
Superconducting targets [26,27] (with a meV electronic band gap), and Dirac materials [28] (with an arbitrarily small gap) have been shown to be promising low threshold targets for both scattering and absorption of low mass dark matter, while molecular magnets [29] have been considered for absorption. Furthermore, a DM particle with mass less than ∼ 1 MeV has a deBroglie wavelength that is longer than the inter-particle spacing in typical materials, implying that the DM effectively couples to the collective excitations of atoms (phonons) in the target. Such DM-phonon scattering processes have different kinematics and allow for a greater amount of energy to be extracted from the DM than for scattering off a single free nucleus. This was the idea of Ref. [30], where it was shown that a DM collision can produce multiple phonon modes in superfluid helium, extending the reach of such a target by 2-3 orders of magnitude in DM mass compared to ordinary nuclear recoils. (Note that there is a phase space penalty for emitting multiple states in a restricted configuration, as discussed in detail in Ref. [31].) The emission of single or multiple phonon modes also allows for absorption of meV-eV mass bosonic DM in both superconductors [32] and semiconductors [33].
In this paper we consider the production of a single excitation -an optical phonon in a polar material -from the interaction of a sub-MeV DM particle. Optical phonons occur in all materials with more than one atom in their primitive cell, including e.g. germanium and diamond crystals.
Unlike acoustic phonons, where the atoms (or rather, the ions) in the primitive cell oscillate in phase, optical phonons arise when the inequivalent atoms in the primitive cell oscillate out of phase. The optical phonons are gapped at low momentum with typical energies of 10 -100 meV, which is wellmatched to the total kinetic energy of sub-MeV DM. In a polar material such as GaAs or sapphire (Al 2 O 3 ), the inequivalent atoms in the primitive cell also have different effective charges, such that the coherent, out-of-phase motion of the ions in the optical modes generates a strong oscillating dipole.
In Ref. [34], we argued that this dipole is particularly advantageous for DM interactions through a dark photon mediator, which can couple directly to the dipole. Polar materials moreover tend to be semiconductors or insulators, which means that the dark photon field does not experience the strong screening effects inherent to conductors (e.g. superconductors [26]). GaAs and sapphire are also wellunderstood materials, where the technology already exists to make ultra pure single-crystal targets.
This is in contrast to Dirac and Weyl materials, which show similar theoretical promise for sub-MeV DM [28] but are not yet feasible for large-scale high-quality synthesis.
In our previous analysis [34], we used several analytic approximations that allowed us to obtain the DM scattering rate in the isotropic limit for a relatively simple polar material, GaAs, which has only 6 phonon modes. In the present paper, we study a more complex material, sapphire, which we argue is better suited for direct detection. To this end we employ more advanced numerical condensed matter techniques, notably Density Functional Theory (DFT), which allows us to accurately compute the scattering rate in sapphire and to validate the analytic treatment of GaAs used in Ref. [34]. A key reason why sapphire is a more promising target is that its crystal structure is anisotropic, implying that the DM scattering rate will depend on the angle of the DM momentum with the primary crystal axis. This manifests itself as a modulation in the rate with the sidereal day, a smoking gun signature for DM that can be used to test the origin of any signal. 1 The outline of this paper is as follows. In Section II A, we begin by presenting the theoretical case for polar materials and elaborate on their benefits, focusing specifically on GaAs and sapphire. In the remainder of Section II, we describe in detail the crystal structures, the method for computing phonon band structures, and set up the formalism for calculating the direction-dependent DM scattering rate.
Next, we consider a few specific model benchmarks: in Section III, we present the reach and daily modulation for DM scattering via an ultralight dark photon mediator, or millicharged DM with photonmediated scattering. Section IV considers nucleon scattering mediated by a light scalar, where we find 1 Other proposals with directional sensitivity to low mass DM scattering include using graphene [18] or Dirac material targets [28] for DM-electron scattering, and taking advantage of direction-dependent thresholds for defect production in crystals [25] or nuclear recoils in semiconductors [35]. modulation rates as large 30% for m X ∼ 50 keV. We lastly consider absorption of dark photon DM with mass of 10-100 meV into optical phonons and multiphonons in Section V, and conclude in Section VI.
More details on the derivation of DM scattering rates are provided in Appendices A-C, while some results on scalar-mediated electron scattering are contained in Appendix D. Finally, some additional details on our method to estimate the statistical discrimination for the daily modulation signal are given in Appendix E. A detailed description of the experimental setup as well as the calculation of the most important backgrounds will be presented in an upcoming paper. In the meanwhile, for a brief description of the experimental setup and estimates of the backgrounds we refer the reader to Ref. [34].

II. POLAR MATERIALS
In this section we first lay out the qualitative features that make polar materials excellent targets for sub-MeV dark matter. For the purposes of this paper, we focus on the examples of GaAs and sapphire 2 since they are commonly used and well-characterized materials, including in some existing or proposed direct detection experiments 3 . We then review the crystallographic properties of GaAs and sapphire, setting up the theoretical framework necessary for performing DM scattering calculations.

A. Advantages for direct detection
As briefly discussed in the introduction, polar materials have several features that make them attractive as targets for the scattering and absorption of light DM. These properties are: i) Even in the limit of low momentum transfer, a relatively large energy deposition is possible when scattering into optical phonons; ii) Sapphire has an anisotropic crystal structure, allowing for directional detection; iii) The optical response allows for dark photon interactions, but is still sufficiently weak that screening effects do not hinder the sensitivity.
We detail each of these features in turn.
2 Note that we will use sapphire or Al2O3 here interchangeably to mean crystalline aluminum oxide (Al2O3); this is also sometimes called corundum (α − Al2O3) in the literature. 3 Notably, the CRESST collaboration recently published results on ∼GeV DM with a sapphire target [36]. Here the deposited energy is measured in phonons, but the initial scattering is a nuclear recoil. In contrast, for sub-MeV mass scales, we are considering the process where single phonons are directly excited by the DM. The acoustic modes (labeled transverse (TA) and longitudinal (LA) acoustic) have the standard gapless, linear dispersion at |q| ≈ 0 ("Γ" point in Fig. 2.) The slope is given by the speed of sound c s = ω/|q| near |q| → 0, where the longitudinal sound speed is c s ∼ 4000 m/s in GaAs and c s ∼ 10 4 m/s in sapphire, though for sapphire the sound speed is somewhat direction dependent [39]. In the long wavelength limit, these modes carry no energy, as they correspond to translations of the lattice as whole. In this sense, the acoustic modes can be considered as the Goldstone modes 4 of the spontaneous breaking of the translation invariance by the crystal. The optical phonons are not protected by Goldstone's theorem, and at q ≈ 0 they can be thought of as a standing, non-propagating 4 There are no Goldstone modes associated with the spontaneous breaking of the rotation invariance: in the presence of a broken translational symmetry, rotations do not give rise to a linearly independent set of Goldstone modes (see e.g. [40,41]). wave which stores a finite amount of energy.
A priori, the dark matter can excite both the optical and acoustic modes, but the energy deposited in the acoustic modes is much smaller and is only detectable in the most optimistic circumstances.
Concretely, for m X MeV, the DM momentum m X v keV is sufficiently small that it is only possible to excite a phonon mode within the first Brillouin zone. Consider a DM scattering with momentum transfer q and energy deposition ω, which excites a single acoustic phonon; the phonon must absorb all of the energy and momentum transferred. This leads to the scaling with v ∼ 10 −3 the DM velocity and assuming the speed of sound for sapphire. The threshold for near future devices will be at best in the 10 − 100 meV range, which means that single acoustic phonon excitations from light DM will be difficult or impossible to detect, depending on m X . However, the scaling in (1) does not apply for the optical modes since they have an energy of ω ∼ 30 meV or more as |q| → 0, as is evident from Fig. 2.
The gapped dispersion of optical phonons is a particularly appealing feature, as it allows nearly the maximum amount of DM kinetic energy to be extracted in the scattering, even when the momentum transfer is much less than a keV. This is in contrast to recoils off free nuclei, where the energy deposited from light DM is much less than the initial DM kinetic energy. The presence of optical phonons is also advantageous compared to a material such as superfluid helium. Superfluid helium does have gapped quasiparticle excitations (rotons), but they only occur at high q and are much lower energy that the optical phonons in a solid. Since single phonon production in superfluid helium is undetectable in the foreseeable future, one must resort to multi-phonon production to break the relation in (1), as was demonstrated in Refs. [30,31]. However, the rate is suppressed since this is a higher order process relying on anharmonic phonon couplings. On the other hand, in Sec. IV we will see that for scalar-mediated DM the rate for producing optical phonons in GaAs and sapphire is also somewhat suppressed due to destructive interference between the different atoms in the primitive cell. The end result is that, for the scalar-mediated model, the reach for polar materials is comparable to that of superfluid helium.

Crystal Properties and Directional Detection
The crystal structure of polar materials can be anisotropic. This is the case for sapphire, which has a rhombohedral lattice structure and therefore a primary crystal axis. This anisotropy manifests itself both in the spectrum of the phonons (see Fig. 2) and in the strength of their coupling to the dark matter. The latter feature depends on the type of mediator, and will be discussed in detail in Secs. III and IV. Regardless of the type of mediator, this means that the scattering rate will modulate with the sidereal day as the angle between the primary crystal axis and the DM wind changes due to the rotation of the Earth. In other words, different regions of the Brillouin zone are sampled at different times of the day, which results in a changing rate due to the q-dependence of the phonon energies and DM-phonon coupling. We set up our notation and conventions for the directionally dependent rate in Sec. II D.

Optical Properties and Dark-Photon Mediated Scattering
The optical response is particularly important if the DM scattering is mediated by a dark photon that is kinetically mixed with the SM photon. Since the dark photon has a coupling to the electromagnetic current, naively the best target for direct detection would be a population of free charges, such as the conduction electrons in a metal. However, the same free charges screen the dark photon field at the low frequencies that are of interest. The ideal material therefore has few conduction electrons but a large polarizability. Superfluid helium fails the latter criterion, as the polarizability of a He atom is very small, rendering helium transparent to both SM and dark photon fields [31].
To compare the characteristics of polar materials with other candidate targets like superconductors and Dirac/Weyl materials, it is useful to express the electromagnetic response in terms of the longitudinal and transverse in-medium polarization tensor Π T,L , which is expressed in terms of the relative permittivity of the material, and where q 2 = ω 2 − |q| 2 .
Accounting for the in-medium effects, the matrix element for scattering can be written as where P µν T,L are the projection operators for transverse and longitudinal polarizations, κ is the vacuum mixing parameter between the dark and SM photon, m A is the dark photon mass and g X is the gauge coupling of DM with the dark photon. Thus, as demonstrated by previous work on dark photon mediated interactions [27,42], the in-medium polarization tensor is essential in determining the reach.
Superconductors were the first targets to be considered for direct detection of sub-MeV DM [27].
In the limit of |q| ω, the dielectric function in metals and superconductors displays Thomas-Fermi screening: with λ 2 TF = 3e 2 n e /(2E F ) (few keV) 2 . The maximum momentum transfer for sub-MeV DM-electron scattering is |q| = 2m X v X keV, such that metal tends to be very large. This screening severely limits the sensitivity of superconductors to dark photon mediated DM.
This screening can be lifted in Dirac and Weyl materials. These have an arbitrarily small gap for excitations of electrons to the conduction band, but a gauge symmetry in the material protects the photon from obtaining a large in-medium polarization [28]. This effect can be understood as a result of gauge invariance, where charge is renormalized but the photon obtains no mass. Calculating a simple one-loop polarization graph with the linear dispersion typical of Dirac materials near the Dirac point, one obtains the dielectric response where b is the background dielectric constant supplied by the nuclei and bound electrons, and v F is the Fermi velocity, which is typically 10 −3 − 10 −2 . The resulting dielectric constant is typically an O(1) number, such that excellent reach to dark photon mediated DM can be obtained. However, Dirac/Weyl materials are still the subject of intense research, and it is not yet known how to fabricate large, ultra pure samples needed for DM detection.
In polar materials, there is a gap for electronic excitations on the order of 1-10 eV, so there is no screening by conduction electrons. And while electron excitations are forbidden, dark photon mediated DM can instead couple to the dipole moment of the optical phonons. The interaction is only screened by the valance electrons, an effect which is encoded in the high frequency dielectric constant ( ∞ ).
Its value for GaAs and sapphire is in the O(1 -10) range, and we derive the screening of the dipole interaction in detail in Appendix B. Thus polar materials satisfy the criteria of large polarizability with little screening by free charges. Compared to Dirac/Weyl materials, there is some penalty to coupling though a dipole moment but the phase space for the scattering process is larger, such that the projected reach is comparable. Polar materials moreover have the advantage that the technology already exists to fabricate the bulk, ultra pure targets needed for detection.

B. Crystal properties
GaAs adopts a cubic zincblende crystal structure (space group F-43m) with two atoms (Ga and As) in the primitive unit cell (left panel of Fig. 3). These two atoms in the primitive cell give a total of six degrees of freedom corresponding to the six phonon modes as was shown in the left panel of Fig. 2. The covalent bond in GaAs is polar, like in other III-V compound semiconductors, owing to the moderate difference in electronegativity between the Ga and As ions. This results in the Ga and As carrying opposite net effective charges. The phonon branches corresponding to the out-of-phase motion of these net-charged ions will therefore couple to electric fields, hence the name "optical phonons".
Sapphire with the chemical formula Al 2 O 3 has a more complex rhombohedral crystal structure (space group R-3c) with four Al and six O atoms in the primitive unit cell (right panel Fig. 3). Each Al ion is six-coordinated with oxygen, which form a corner-sharing network to make up the crystal lattice. The corresponding phonon spectrum was shown in the right panel of Fig. 2. Because of the inequivalent in-plane and out-of-plane crystal directions, sapphire has a primary crystal axis, which implies that the scattering rate depends on the angle between the momentum transfer and the primary axis. It also means that the isotropic approximation used in Ref. [34] does not hold for sapphire, though we will see that it works well for a more symmetric crystal like GaAs. For sapphire, the scattering rate must be computed numerically, using first-principles methods that incorporate the crystalline and chemical specificity of the sapphire crystal.
Tab. I lists some useful quantities for both materials from the point of view of DM-phonon scattering.
While the quantities are all temperature-dependent, the differences between room-temperature and liquid helium temperatures are typically percent-level or less [43]. Notably, the dielectric constants are O(1) quantities, which is relevant for dark photon mediated DM, as discussed in the previous subsection. We list both the low frequency ( 0 ) and high frequency ( ∞ ) dielectric constants, where the high frequency case refers to ω above the optical phonon frequencies, but still well below the electronic band gap of the material. At high frequencies, the ions in the lattice have no time to respond to a rapidly changing electric field, and the dielectric function only receives contributions from the valence electrons. At low frequencies, the optical phonons contribute to the dielectric function as well, such that 0 > ∞ in a polar material. We also note that our first-principles calculations are carried out at zero temperature, providing a close reference for liquid He temperatures.

C. Theoretical description of phonons
With advances in first-principles modeling of materials and in high-performance computing, it is routine to calculate the electronic and vibrational properties of crystals from first principles. For DM direct detection in particular we need the phonon spectrum over the whole Brillouin zone, since this is an input for the DM scattering (or absorption) rate calculation. Here we briefly discuss how these calculations are performed and establish the notation for the remainder of the paper. Readers familiar with the subject or only interested in the results can choose to skip the remainder of this section.
The positions of the atoms (or ions) in the crystal are denoted by u j,l + r 0 j + l, where u j,l is displacement of the atom relative to its equilibrium position, r 0 j is the equilibrium position of the atom relative to the origin of the primitive cell and l is a lattice vector labeling the primitive cell. The index j therefore runs over the atoms in the primitive cell. In what follows, a boldface symbol always refers to a tensor or vector in position or reciprocal space. The potential energy V is a function of the displacements and can be expanded as where the V l,l ,j,j etc are the force constants. The force constants can be calculated from ab initio density functional theory methods. For this work we use density functional theory (DFT) as implemented in the VASP package [44] with full calculation details given in Appendix A. Firstly, the equilibrium crystal lattice and atomic positions are found by minimizing the forces on the atoms and stresses in the crystal cell. From this optimized crystal structure, the force constants can be calculated using two different methods. The first displaces atoms in the cell in symmetry-inequivalent directions, calculates the resulting forces on the atoms in the unit cell, and from these builds up the force constant matrix.
The second method, the linear response method, uses density functional perturbation theory (DFPT) to calculate the forces. In this work, we will use the former method, known as the 'frozen-phonon' method, to calculate the full phonon spectrum as it is computationally less expensive. For the Born effective charges, we will use DFPT.
In the harmonic approximation, one only considers the leading non-vanishing order, keeping only l,j,l ,j . The displacement operator is then quantized in terms of phonon modes: where theâ ν,q (â † ν,q ) are the creation (annihilation) operators of a phonon mode in branch ν with momentum q. The total number of branches is 3n, where n is the number of atoms in the primitive cell. ω ν,q is the energy of phonon branch ν with momentum q, and e ν,j,q is the unit vector indicating the direction of oscillation of atom j for phonon branch ν. Finally, m j is the mass of atom j, and N is the number of cells in the lattice. The q form an N -point discretization of the Brillouin zone, such that the variance of the displacement u j,l · u j ,l is an intrinsic quantity under N → ∞.
In Fourier space, the equations of motion for the displacements then reduce to a standard eigenvalue problem for a given momentum vector, j D q,j,j · e ν,j ,q = ω 2 ν,q e ν,j,q , where the eigenvectors are normalized such that j e * ν,j,q · e ν,j,q = 1. The dynamical matrix D q,j,j is given by For a rigorous derivation, see e.g. [45]. Once the force constants are known from a DFT calculation, the eigenvalue problem can be solved for ω ν,q and e ν,j,q using the post-processing software package phonopy [38]. From these the phonon-derived properties such as the phonon band structures shown in Fig. 2, can be calculated.
The dynamical matrix receives an additional non-analytic contribution from the Born effective charges, which modifies the frequencies of the LO phonons. The Born effective charge is the electric charge that effectively contributes to the polarization induced during an atomic displacement, and is used to quantify the coupling between optical phonons and electric fields. Formally, the Born effective charge tensor Z * is defined as the change in polarization P resulting from a displacement u: where Ω is the unit cell volume, and e is the electric charge. It can also be written in terms of the change in the force F in a direction i with respect to a homogeneous electric field E in direction j. The Born effective charge tensor Z * can be calculated using DFPT. This uses density functional theory to calculate the response of the system to a finite electric field as detailed in [46,47]. The Born effective charges are hence computed from the change in the Hellmann-Feynman forces and mechanical stress tensor due to the changes in the wavefunction.
The calculated Born effective charges for GaAs and Al 2 O 3 are and We note that due to the high symmetry of the GaAs crystal, the Born effective charges are a scalar quantity. For Al 2 O 3 , owing to the different anisotropic chemical environment surrounding Al and O atoms, the Born effective charges tensors are in general tensors that differ for inequivalent atoms, as listed in Appendix A 2. For our numerical calculations we use the diagonal, cell-averaged values for Al and O given above.
The LO phonon modes correspond to ions with opposite effective charges moving in opposing directions alongq, causing long-range macroscopic electric fields in a polar crystal. In contrast, TO phonons correspond to oppositely-charged ions moving in adjacent planes parallel to each other, resulting in no long-range Coulomb interaction (see Fig. 1). The additional force created by the electric field interaction with the LO phonon modes results in a frequency change in the LO mode as q → 0.
The lifting of the degeneracy between the LO and TO phonon modes at the Brillouin zone centerthe so-called LO-TO splitting -can be calculated by including the non-analytic contribution to the dynamical matrix, given by in Lorentz-Heaviside units. Hence calculating the Z * allows one to determine the corrected LO modes.
Here we use diagonal ∞ tensors, as given in Table I. The ω ν,q and e ν,j,q obtained from phonopy will be the most important inputs for the DM scattering rate calculations. The next missing ingredient is the effective coupling of the dark matter to the displacement operator in (7). This coupling is model dependent and we treat it separately for dark photon and scalar mediator cases in Secs. III and IV respectively.

D. Crystal alignment relative to dark matter flux
Before turning to the scattering rate computation, we first establish our assumptions and conventions regarding the DM velocity distribution and the orientation of the DM wind in the frame of the crystal, which will determine the directional signal. The incoming DM velocity in the lab frame is modeled in a standard way, with a boosted Maxwell-Boltzmann distribution: with v 0 = 220 km/s, and truncated by the escape velocity v esc = 500 km/s. The velocity of the Earth with respect to the DM wind is indicated with v e , with |v e | ≈ 240 km/s on average.
The orientation of v e relative to the crystal changes as the Earth rotates around its axis. Combined with the anisotropic crystal structure, this sources a daily modulation of the scattering rate. We neglect the yearly modulation due to the Earth's orbit around the Sun. The orientation is illustrated in Fig. 4, where θ e ≈ 42 • is the angle between the Earth's rotation axis and the direction of its velocity and θ lab gives the latitude at which the experiment is constructed. We choose the crystal orientation and coordinate system such that the z-axis in the crystal frame is aligned with the Earth's velocity at t = 0. For GaAs, we choose the z-axis in the crystal frame along one for the faces of the cubic lattice. For sapphire, the z-axis is taken to be aligned with the primary crystal axis, which is the axis defined by the Al atoms in Fig. 3. This implies that at t = 1/2 day, the primary axis of the sapphire crystal is at an angle of roughly 90 • with the DM wind. While we have not explicitly optimized for the crystal configuration, we expect that the choice here should (nearly) maximize the amplitude of the daily modulation since the biggest anisotropies in sapphire are that between the crystal axis and the axis perpendicular to it.
Since we explicitly orient the crystal relative to the dark matter wind, there is no dependence of the DM flux or scattering rate on the latitude at which the experiment is located. As a function of time, the unit vector of v e in the crystal coordinate frame iŝ sin θ e cos θ e (cos φ − 1) with φ = 2π × t/24h the angle parametrizing the rotation of the Earth around its axis.

III. DARK PHOTON MEDIATED SCATTERING
We begin by considering Dirac fermion DM that interacts with the SM via a kinetically mixed dark photon. The model is defined by the vacuum Lagrangian where κ 1 is the kinetic mixing parameter, and g X and m A are respectively the gauge coupling and mass of the A µ field. For finite m A , one can perform a diagonalization to the mass basis, where the electron picks up a small charge (in vacuum) of κe under the dark photon. On the other hand, in the limit where m A → 0, we can perform a field redefinition A → A − κA to write the Lagrangian as where the dark matter X has a small charge e ≡ κg X under the photon. In either of these cases, there is a coupling of the DM current to the electromagnetic current that is proportional to κg X e.
For sub-MeV dark matter, the relic abundance can be explained by freeze-in [48][49][50] via the outof-equilibrium process e + e − → XX. Since this production rate is proportional to the coupling combination κ 2 g 2 X , requiring that X is 100% of the dark matter predicts also a compelling target for DM scattering off SM particles. Requiring that the DM-A coupling satisfies self-interaction bounds and that the kinetic mixing κ is consistent with dark photon constraints, one finds that m A 10 −11 eV [51]. Since this mass is much smaller than a typical in-medium effective photon mass, we can take the massless A limit. We are then in the situation given by Eq. (18) above, where we can treat the DM as a millicharged particle for the purposes of our calculations. One could also consider the model above with m A = 0 as a specific model of millicharged DM 5 .
For interactions mediated by an ultralight dark photon, the long-range coupling of DM with a phonon in the crystal is then similar to that of electrons with phonons, but with an amplitude suppressed by e /e. Here, we specifically mean the interaction associated with a 1/r 2 Coulomb field.
(There are also short-range electron-phonon interactions in a material, but in a polar material these interactions are typically much smaller.) The long-range interaction between electrons and phonons in semiconductors and insulators is described by the Fröhlich Hamiltonian [54,55]. Physically, an electron injected into the crystal sources a small electric field, which induces an oscillation of the ions via the Born effective charges. This oscillation can then be identified with a phonon mode.
Since the DM scattering in polar materials behaves similarly to millicharged dark matter, we can directly use the Fröhlich Hamiltonian as a description for this process. The main difference with the electron case is that the DM is a free particle, while for electrons the appropriate in-medium wavefunctions must be included. On a practical level, this is a simplification of the computation since the plane wave approximation is sufficient to describe the DM. In the following section we will summarize the most important formulas and numerical results, and provide the relevant derivations in Appendix B.

A. Fröhlich interaction
Adapted for the DM case, the Fröhlich matrix element 6 for a periodic lattice is given by [56] M q+G,ν = iee This result is derived in Appendix B. Here ν, j and G are the phonon branches, the atoms in the primitive cell and the reciprocal lattice vectors, respectively. The momentum transfer is given by q + G, while q is the momentum of the phonon restricted to the first Brillouin zone. e is the electron charge 7 and Ω is the volume of the primitive unit cell. In general ∞ is a tensor, though it is well approximated by the scalar quantity ∞ times the unit tensor. The phonon eigenvectors e * j,ν,q , energies ω ν,q , and the Born effective charges Z * j are all computed from first principles, as described in Sec. II C. A similar formulation is often used in the literature to describe the coupling of electrons with optical phonons, albeit with the inclusion of the electron wavefunctions in the medium.
The expression above can be understood more intuitively by taking the isotropic and longwavelength approximation. In this limit, and assuming a single phonon branch ν contributes, the matrix element simplifies to where we have dropped the reciprocal lattice vector (since the result is dominated by G = 0). For a given phonon eigenmode, the physical displacements of atom j are proportional to e * j,ν,q / 2m j ω LO ; weighted by eZ * j , this is simply the dipole moment corresponding to the lattice displacements. The eigenvector is dotted into q, selecting for the longitudinal mode, while the overall 1/|q| scaling is that expected for a dipole-charge coupling. Finally, the field generated by the dipole is screened by the valence electrons, which accounts for the 1/ ∞ factor.
The above expression can be further simplified for certain crystals. In Ref. [34], we considered GaAs, which has a single LO phonon branch. As discussed in the previous section, GaAs has a cubic symmetry with Z * Ga = −Z * As and m Ga ≈ m As . We can then make the additional approximation in the long-wavelength limit, 6 Note that (19) differs from the expression in [56] with a phase factor, as we have used a different convention for the phase of e * j,ν,q . 7 We adopt Lorentz-Heaviside units where the vacuum permittivity ε0 = 1 and e = √ 4παem, while [56] uses SI units.
with µ the reduced mass of the Ga and As atoms. In the second equality, we expressed the Born effective charge Z * in terms of the measured quantities ω LO (the frequency of the LO phonon as q → 0), ∞ , and 0 . The derivation for this identity is given in Appendix B. None of these simplifications apply for sapphire, however, and there we must numerically sum over all phonon eigenmodes.

B. Reach
The scattering rate for an incoming DM particle with velocity v i is obtained from Fermi's golden rule, where the momentum integral is over the first Brillouin zone. For scattering for sub-MeV dark matter, we simplify the matrix element in Eq. (19) by observing that the momentum transfer q ∼ v X m X is small compared to the size of the Brillouin zone, except for m X approaching 1 MeV. In addition, the matrix element is proportional to 1/|q| and therefore the rate is dominated by those phonon modes well within the first Brillouin zone. We can therefore neglect Umklapp processes where phonons are created with wavevectors outside the first Brillouin zone; this amounts to setting G = 0 in Eq. (19).
(We expect that the reach does extend to higher masses via such processes, but will not consider this regime further here.) The integral above can be performed analytically if we take the isotropic limit for the matrix element in Eq. (22): where v i is the initial velocity of the DM and the Heaviside Θ-function enforces energy conservation.
For the full numerical result as well as for the analytic approximation, the scattering rate for the target is obtained by integrating over the initial DM velocities: with f (v i ) the dark matter velocity distribution in Eq. (14) and ρ T the mass density of the target material.
To estimate the reach, we compute the expected 90% exclusion on e assuming zero events observed with no expected background. 8 To compare with existing constraints and other proposed experiments, 8 Backgrounds from coherent photon and coherent neutrino scattering are estimated to be no higher than a few events / kg year exposure [34].
which corresponds to the typical cross section of dark matter with a bound electron, e.g. in a semiconductor-based experiment. µ Xe is the DM-electron reduced mass, α is the fine-structure constant and m e is the electron mass. The result is shown in Fig. 5 for both GaAs and sapphire. For GaAs, we compare the isotropic limit with the numerical result including phonon eigenmodes and find excellent agreement. Also shown are existing stellar cooling [58], BBN [59] and Xenon10 [60] constraints, as well as the projected reach of other experimental proposals [12,26,28,50]. Interestingly, we find that as little as a gram-month exposure would suffice to reach the freeze-in benchmark. In the sub-MeV range, an experiment based on a Dirac material [28] is currently the only other proposal which could compete with polar materials. Given that Dirac materials have not yet been fabricated in the quantities needed for a dark matter detector, we expect that the polar material concept could be realized on a substantially shorter timescale. Also shown in Fig. 5 (dashed blue) is the expected sensitivity for sapphire if one requires a daily modulation signal at 2σ. We elaborate on the daily modulation in the next section.  6. Modulation of the scattering rate for sapphire over a sidereal day, assuming a 25 meV threshold.

C. Daily modulation
The anisotropy in the crystal structure induces a dependence of the scattering rate on the crystal orientation, which translates to a modulation over the sidereal day. Here there are two effects that lead to modulation: the directional dependence in the phonon couplings to the DM model, and in the phonon energies. For GaAs, which has a high degree of symmetry, we find that this modulation is negligible. Instead, for sapphire, there is a sizable anisotropy in the DM scattering rate. In the rest of this section, we discuss the dominant effects and present results on the modulation.
The daily modulation in sapphire for several DM masses is shown in Fig. 6. Here we assumed a threshold of 25 meV, well below the energies of the optical phonons. Since any backgrounds are expected to be either flat in time, or at least out of phase with the sidereal day over many periods, this can be used as an additional indicator of a DM signal. Assuming a kg-year exposure, we can estimate the cross section needed to reject the null hypothesis of non-modulating scattering at the 2σ level. This is given by the dashed blue line in Fig. 5, which requires that in 50% of the simulated signal datasets, the null hypothesis can be rejected. The shaded band indicates the ±1σ band around the mean: specifically, the cross section needed if we instead require this to be true in 16% or 84% of the simulated datasets (assuming only statistical fluctuations). We refer to Appendix E for details on our statistical treatment.
To understand the origin of the modulation in Fig. 6, it is useful to deconstruct the total scattering rate in terms of the rate from individual phonon branches. In Fig. 7a  a single mode appears to dominate the matrix element. This is also the most energetic mode in the spectrum, mode 30, with ω ≈ 104 meV. (We label the 30 phonon modes according to their energy in the vicinity of the origin of the Brillouin zone, from least energetic "mode 1" to most energetic "mode 30".) Fig. 7a also highlights the directional-dependence and the q-dependence of the phonon couplings, which enters directly into the scattering rate. As can be seen from Fig. 4, we have assumed that the crystal axis is aligned with the DM wind at t = 0, so that the scattering is preferentially along the crystal axis (the degree to which this is true depends on the DM mass, of course). Meanwhile, the crystal axis is nearly perpendicular to the DM wind at t = 0.5 day, with the dominant scattering into those modes which have large dipole along the q ⊥ directions. Fig. 7a thus suggests that the highest rate occurs along the q direction, corresponding with t = 0, which is consistent with the location of the maximum in Fig. 6 for m X 50 keV.
For m X 50 keV, mode 30 is kinematically forbidden, and the modulation pattern therefore changes. In this case, the lower-lying mode 16 (ω ≈ 60 meV) takes over, as shown in Fig. 7b. For such low m X , threshold effects from the directionally-dependent phonon energies dominate the modulation, . Adobe Acrobat reader is required to view these animations.
amplitude along the crystal axis as compared to mode 30. As such it is subdominant, unless mode 30 is kinematically inaccessible. It is worthwhile to inspect the modulation patterns of the contributions from mode 16 and mode 30 separately, which we present in Fig. 9 for several DM masses. One can see that mode 16 gives rise to a much larger amplitude, and that its phase is shifted with respect to that of mode 30, especially at low mass. This explains the dramatic change in the modulation pattern in Fig. 6 for m X = 25 keV, for which mode 30 is forbidden.
The amplitude of the modulation decreases for higher m X , where larger q values are sampled in the Brillouin zone. We expect that in this limit, scattering starts to transition towards scattering with a single nucleus, which is isotropic. In other words, at high momentum transfer the DM is blind to the long-range crystal structure. In practice, this effect manifests itself in a gradual randomization of the eigenvectors as q is increased on a particular phonon branch. To illustrate this effect, the animation in Fig. 8c shows mode 30 for a point near the edge of Brillouin zone with |q| ∼ 1 keV, which displays less coherent oscillations within the unit cell.
Finally, we comment on the theoretical uncertainties of our first principles calculation. As a proxy for the uncertainty, we have also calculated total scattering rates and modulation patterns using inequivalent rather than averaged Born effective charges for each of the different Al and O atoms in the primitive cell (see discussion in App. A 2). We find small differences in the total rate at the level of a few percent for m X 75 keV, though the difference grows for lighter DM, up to a factor of ∼ 4 for m X ≈ 25 keV. The modulation amplitude differs by roughly a factor of two between the two assumptions for m X ≈ 25 keV, but this difference reduces to roughly 50% for m X 50 keV. We find that the overall modulation pattern remains unchanged.

IV. DM-NUCLEON SCATTERING
In this section, we consider the benchmark where DM couples primarily to nuclei through a light scalar mediator. The underlying interactions are where we assume a scalar DM particle, X, and an identical coupling y n to both neutrons and protons.
We further take m φ small compared to the typical momentum transfer; the so-called light-mediator regime. Similarly to the dark photon mediator, this model is already subject to a number of astrophysical and terrestrial constraints, and we refer the reader to Ref. [51] for a detailed discussion. 9

A. DM-phonon form factor
The scattering of DM off the nuclei of a lattice is similar to the scattering of cold neutrons, except for an additional form factor associated with the light mediator. There is extensive literature on the scattering of cold neutrons (for a review, see for example Ref. [45]), as this process is important to experimentally measure phonon dispersion relations. In the cold neutron case, the cross section for a 9 In particular, this model is most motivated for sub-component DM. In our figures we will however assume that X is 100% of the DM, as the reach is easily rescaled to a particular subcomponent fraction.
neutron to scatter off a single nucleus N is written as σ = 4πb 2 N , where b N is the average neutronnucleus scattering length. Accounting for the lattice structure requires summing over nuclei, weighted appropriately by the phonon wavefunctions. For DM scattering via a light scalar, the techniques for cold neutron scattering in the lattice can then be directly applied.
For a nearly massless mediator, the differential cross section diverges as 1/|q| 4 , though the divergence is cut off by the experimental threshold. Since this threshold varies for different experiments, it is conventional to introduce an effective DM-nucleon cross section, where q 0 ≡ v 0 m X is a reference momentum and µ Xn is the DM-nucleon reduced mass. (The choice for q 0 is merely a convention, and drops out in the scattering rate.) We can similarly define an effective DM-nucleon scattering length from the relation σ n = 4πb 2 X . The primary quantity for describing the response of a crystal to an incident neutron or DM particle is the dynamic structure factor S(q, ω). Here we provide only the final expressions for S(q, ω); we summarize their derivation in Appendix C. In particular, if the momentum transfer is below the size of the Brillouin zone, S(q, ω) can be written as: where the sum runs over the phonon modes (ν). The phonon form factor F ν (q) is given by where the sum runs over the atoms j in the primitive cell, and A j is the atomic mass number. Here we have used that for our benchmark model, the DM-nucleus scattering length is given by b j = A j b X , since we have a coherent sum over all nucleons in the long-wavelength limit. The Debye-Waller function, W j , measures the average motions of the atoms 10 in a phonon excitation, and is given by Note the quantity is finite, since the 1/N factor is compensated by the sum over all phonon modes k.
For all our results, it is a good approximation to take W j ≈ 0, as the spread on the motions of the atoms is small compared to the inverse momentum transfer. Taking m j > 16 GeV since the lightest nucleus is O and ω ν,k > meV for the most optimistic experimental threshold, we still find √ m j ω ν,k = 4 keV, which is larger than the typical momentum transfer for scattering of sub-MeV DM. 10 Here we use 'atom' and ''nucleus' interchangeably to refer to the scattering center.
As derived in Appendix C, the integrated scattering rate per unit of target mass is given in terms of the dynamical structure factor, where ρ T is the mass density of the target and Ω is the volume of the primitive unit cell. The (q 0 /|q|) 4 form factor is the result of the light mediator. The expressions for the massive mediator limit can be obtained by dropping this form factor and substituting q 4 0 with m 4 φ in equation (28).

B. Reach
Contrary to the case with a dark photon mediator, all atoms in the primitive cell contribute with the same sign to the form factor in (30). The modes which couple most strongly to the dark matter are those where all atoms move in the same direction, and thus interfere constructively in (30). In addition, the q · e ν,j,q factor indicates that only the longitudinal modes with motion of the atoms parallel to the momentum q contribute, to leading order in the small q expansion (see Appendix C). Thus, the DM coupling to the longitudinal acoustic mode will be the largest. The optical modes also contribute, but since at least some atoms move in opposite directions, there are inevitably cancellations (destructive interference) between the contributions of various atoms. These effects can be seen most easily for GaAs, where the form factor can be approximated by where the + sign applies for the LA mode and the − sign for the LO mode. We have included relative phases for the motion of the Ga and As, to account for the fact that the motion will not be perfectly in phase away from the long-wavelength limit (see also (9), where the phases appear explicitly in the dynamical matrix).
Since Ga and As have similar mass numbers, we see from the equation above that there is destructive interference for the optical mode, which leads to a suppression of the rate by several orders of magnitude compared to the acoustic mode. For sapphire, the mass hierarchy between the two elements is slightly larger, but since there six O atoms as compared to four Al atoms in the primitive cell, both elements end up contributing a similar amount to the scattering rate. To fully remove the suppression due to the destructive interference, it would be interesting to consider a polar material with a large mass difference between the elements, such as PbS.
Here we use the numerically computed phonon eigenmodes to calculate the scattering rate, while we previously applied the analytic approximation of (33) to GaAs in Ref. [34]. As before, we estimate the reach by computing the projected 90% CL limit under the assumption of no backgrounds and no events observed. The result for both GaAs and sapphire is shown in Figure 10 for a kg-year exposure. The  (33). Also shown is a projection for a superfluid helium target that is sensitive to multiphonon production from DM, with kg-year exposure and meV threshold [31].
analytic approximation for GaAs matches the numerical result very well for the acoustic branch (dark purple), and for the optical branch (light purple) it reproduces the numerical result to within a factor of ∼ 3. As expected, the reach dramatically improves if the threshold is low enough to pick up the acoustic modes, and in this case substantially outperforms a superfluid helium detector in multiphonon mode [31].
If only the optical modes are accessible, the reach is comparable or somewhat weaker than that of superfluid helium. In this case only one LO mode contributes for GaAs, and it is imperative that the threshold is lower than 30 meV. For sapphire, there are several modes in the spectrum which contribute comparably to the total rate. The cross section and the reach therefore differ for different experimental thresholds in Fig. 10, as more phonon modes can be accessed for lower thresholds. This is to be contrasted with the case of the dark photon mediator, where mode 30 alone was responsible for almost all of the rate, provided that it is kinematically accessible. The threshold dependence of the rate is thus not present for the dark photon mediator, and could be a discriminating variable between the models, should a signal be observed. As for GaAs, the sapphire reach would increase substantially if the acoustic phonons could be accessed. In particular, the improved reach for the 25 meV threshold and m X > 200 keV in sapphire is due to one of the acoustic modes: at this point, the momentum transfer becomes just large enough to access a portion of the acoustic branches (see Fig. 2). This substantially enhances the rate, giving rise to the feature in Fig. 10.

C. Daily modulation
Similar to the case of dark photon mediated scattering, the rate modulates with sidereal day due to anisotropies in the phonon spectrum and the phonon form factor. Here the directional dependence of the form factor is encoded in the eigenvectors e ν,j,k in Eq. (30). The modulation rates for different DM masses and possible experimental thresholds are shown in Fig. 11 for sapphire; similar to before, we find much smaller modulation rates for GaAs, with sub-percent modulation except for m X 30 keV. As for the dark photon mediated scattering, the modulation decreases for larger DM masses, as the e ν,j,k tend to be more randomized for higher q. However, the modulation amplitude drops more slowly compared to dark photon mediated scattering, and even for m X ≈ 200 keV the modulation can still be as large as ∼ 20%.
To understand the dependence of the modulation on the threshold, we first observe the lack of a substantial modulation for the lowest (1 meV) threshold. The reason is that the acoustic modes dominate in this case. Since all atoms move in phase on the acoustic branches, the primary modulation comes from the anisotropy of the sound speed, which is fairly small. For a higher threshold (>25 meV), we instead rely primarily on the optical modes. As explained in the previous section, in the q → 0 limit, the contributions from the different atoms tend to destructively interfere for the optical branches. The effect of finite q corrections is then to partially remove these cancellations; this effect will vary along different crystal directions, leading to a sizable directional dependence of the scattering rate. (One way of seeing this is to consider the effect of the phase factors in (33).) As the threshold is further increased, fewer optical modes can contribute to the rate. Since each mode has a unique modulation pattern, this means that the total modulation pattern depends on the threshold. In addition, different DM masses sample different regions in the Brillouin zone, which means that the relative weight of the phonon modes shifts as the DM mass is varied. This too has an effect on modulation pattern, as can be seen most clearly by comparing the curves for m X = 50 keV and m X = 200 keV benchmarks in Fig. 11. Both features may help with characterizing the DM mass, should a signal be observed.

V. ABSORPTION OF DARK PHOTONS
The presence of optical phonons in polar materials also makes it an excellent target for absorption of dark photon DM. In the sub-keV regime, dark photons are a viable DM candidate, and can be detected by an optical absorption signal if there is a small mixing with the SM photon. Similar to Section III, we consider the Lagrangian We begin with a review of the absorption of dark photons in optically isotropic materials, such as GaAs. The mixing present in Eq. (34) is modified in the presence of an in-medium polarization for the photon, which can be written as [27,61] Π γγ (q, ω) = ω 2 (1 −n 2 ).
Note that the above result holds for both longitudinal and transverse polarizations, and we have taken the limit of |q| → 0, appropriate for absorption processes, such that we can write Π γγ (q, ω) ≡ Π(ω). n = n + ik is the frequency-dependent complex index of refraction, and is related to the permittivitŷ and to the optical conductivityσ of the material: Note that the real part ofσ, σ 1 , appears in the imaginary part of the polarization tensor Π(ω). It can thus be seen that σ 1 is the absorption rate of SM photons. For energies near the LO and TO phonon frequencies, the permittivity of a polar material can be described analytically as [62] ( where we have included a product over all optical branches ν. Each branch is split into longitudinal and transverse modes with energies ω TO,LO , while γ TO,LO are the damping parameters. ∞ is the contribution of the electrons for energies below the electronic band gap. It is at the LO phonon frequencies where (ω) becomes suppressed. For GaAs, there is one active branch and data on the parameters at low temperatures can be found in Ref. [63]. However, note that the permittivity above does not include multiphonon absorption, and where possible we will supplement the above result with the measured index of refraction.
Including the in-medium polarization from Eq. (35) in the Lagrangian and diagonalizing, we obtain a coupling of the dark photon with the EM current given by κ eff J EM , where the effective in-medium kinetic mixing parameter is where we took ω = m A in the second step. The dark photon absorption rate per unit target mass is then determined in terms of the photon absorption, and can be written as where ρ T is the target density. In Ref. [34], we applied the above result to GaAs. The phonon absorption is temperature dependent, so we have selected low-temperature results whenever available.
For the absorption into phonons (m A <eV), we used calculations of the zero-temperature absorption coefficient α into single and multiple phonons from Ref. [64], where α = σ 1 /n, and we use Eq. (37) to determine n. It can be seen in the right panel of Fig. 12 that using only Eq. (37) misses a large portion of the absorption, due to multiphonons. Fig. 12 also shows that the peak of the absorption is actually at ω LO , even though the photon absorption is peaked at ω TO . This is due to the relatively suppressed κ eff at ω TO . For eV and greater masses, we used room-temperature data onn from Ref. [65].
The absorption of dark photons in sapphire differs from that of GaAs because sapphire is a birefringent material, meaning that the complex index of refraction depends on the polarization of the vector field relative to the optical axis (the crystal axis or c-axis in sapphire). In the optical phonon regime, this should not be too surprising: as discussed in the previous sections, there is significant anisotropy in the phonon dipole moments and energies for modes parallel or perpendicular to the c-axis. Data on the index of refraction is typically quoted separately for ordinary rays ( E ⊥ c-axis) and for extraordinary Effective rate of absorption GaAs, data GaAs, analytic Al2O3, data Al2O3, analytic For dark photons as the DM, we expect the field to have a random polarization with respect to its k-vector and to the orientation of the c-axis. In particular, the coherence time for the dark photon field is ∼ 1/(m A v 2 ) 1µs for the masses considered here and with v ∼ 10 −3 , and so the polarization will change randomly on a time scale much faster than the rotation of the crystal, for instance. As such, we simply take the average of the absorption rate for polarizations perpendicular and parallel to the c-axis, where the subscripts indicate the ordinary (o) and extraordinary (e) directions, respectively.
In the left panel of Fig. 12, we show the effective absorption rate κ 2 eff σ 1 /κ 2 for both polarizations in sapphire, as well as the weighted average we use in computing the sensitivity. The data is obtained from Ref. [66], which compiled measurements at room temperature. Similarly to GaAs, while the strongest absorption into photons is at the TO frequencies, we actually find strong dark photon absorption peaks at the LO frequencies due to the in-medium κ eff . In particular, we find the strongest absorption at the mode with ω LO ≈ 110 meV, which we identified earlier as having the largest dipole moment.
In the right panel of Fig. 12, we compare the room temperature data with the result using Eq. (37) and best fit parameters measured at 77 Kelvin from Ref. [67]. The parameters we used are reported in Frequencies are from Ref. [68] with error bars of less than 0.5%. Widths are from Ref. [67] at 77K, and only reported for the o-ray case; for the e-rays, we adopt the same values as in the o-ray case for similar phonon frequency. All values are quoted in units of 1/cm.
Tab. II. It can be seen that the bulk of the absorption is described by the broad resonances in single optical phonon production, and there is good agreement in the two approximations. Ideally, one would obtain data at even lower temperatures, but we did not find any in the literature. We expect that reducing the temperature further would lead to reduced phonon widths by an O(1) factor, and thus somewhat narrower peaks. Depending on the details of the eventual experimental setup, it may be possible to measure the low temperature absorption rate during a calibration run. Fig. 13 shows the resulting sensitivity to dark photon DM, parameterized in terms of the vacuum kinetic mixing κ. We again assume kg-year exposure and zero background, and find that polar materials provide an excellent broadband target in the mass range of few meV up to 0.1 eV via the multiphonon signal.

VI. CONCLUSION
Except for the simplest of crystals, most materials have gapped lattice vibrations (optical phonons) with energies between 10 meV and 100 meV. This matches the typical kinetic energy of DM in the Galaxy for masses between the ∼ 10 keV warm DM limit and up to 1 MeV, allowing for single optical phonons to be excited in DM collisions with the crystal. We used Density Functional Theory (DFT) methods to compute the rate for DM to create an optical phonon in GaAs or sapphire in the zero temperature limit. Both crystals are examples of polar materials, where the optical phonon modes give rise to long range electric fields in the crystal. This implies a coupling to any DM candidate that scatters through an ultralight dark photon mediator, which is a challenging scenario for other direct detection proposals targeting sub-MeV DM such as superconductors [27] or superfluid helium [34].
In previous work [34], we studied the example of GaAs with an analytic treatment. Here we go significantly beyond the earlier work in several ways: we validated the analytic treatment for GaAs using DFT methods, we extended the calculations to the more complex but potentially more promising example of sapphire, and we studied the directional dependence of the scattering rate in sapphire. In particular, sapphire has higher energy optical phonon modes that can be more readily accessed in an For m A > eV, the dark matter is absorbed into electron excitations. Also shown are existing direct detection constraints from DAMIC [69], SuperCDMS [13], Xenon10 [61,70], and Xenon100 [33,70] (shaded blue) and constraints on emitting dark photons in the Sun [42,71]. The dotted lines are projections from Al superconductors [32], Ge and Si semiconductors [33], Dirac materials [28] and molecules [72]. See Ref. [70] for absorption on GaAs for m A > eV. Molecular magnets [29] have a reach in the κ ∼ 10 −17 − 10 −15 range for 10 −2 eV m A 10 eV. experiment, and the crystal anisotropy leads to a sizable directional dependence, which is manifest as a modulation in rate over a sidereal day. This directional dependence is much smaller in GaAs due to the more isotropic nature of the crystal. The dependence of the modulation pattern and amplitude on the target material suggests that if a signal were to be observed, one could employ a number of different polar material targets to extract details on the DM model and further confirm its cosmic origin.
We analyzed sub-MeV DM scattering via both a dark photon (vector) mediator and a scalar mediator. For the dark photon mediator, the scattering occurs dominantly into optical phonon modes, and the resulting reach and daily modulation are shown in Fig. 5 and Fig. 6 respectively. In the case of the scalar mediator, the best sensitivity can be obtained if acoustic phonon modes are accessible; the reach and modulation are shown in Fig. 10 and Fig. 11. The modulation pattern and which phonons are excited thus depend strongly on the DM model, and a definitive observation of the modulation could in principle be used to infer the DM mass and mediator spin. For the scalar mediator, we studied the example where the mediator couples to nuclei, but our analysis also applies to the scenario where the dark matter, and their sensitivity to various models. "Electron" and "Nucleon" refer to a scalar coupling to electrons and nuclei respectively. ( ) refers to cases where sensitivity is expected, but no calculation has been performed at this time. In addition, molecular magnets [29] have been shown to have good reach to dark photon absorption, and may be sensitive to scattering and/or scalar absorption processes.
scalar mediator couples to electrons. This is because the scattering into phonons is really a scattering off of the nucleus plus the inner-shell electrons rather than just the nucleus, so that one can estimate the rate by substituting the atomic mass numbers in Eq. (30) with the number of bound electrons in each atom. The results of this procedure are summarized in Appendix D.
Polar materials are also sensitive to the scenario where the DM is a boson with mass below ∼ eV, where the DM could be absorbed into single or multi-phonon excitations. Fig. 13 shows the reach for dark photon DM, for which the absorption rate can be related to the measured optical conductivity of the material. We expect that polar materials could also be sensitive to the absorption of scalar/pseudoscalar DM with a coupling to nucleons and/or electrons, but the absorption rate on optical phonons should be subject to the same destructive interference that we found for scalar-mediated scattering (Sec. IV). This implies that the multiphonon absorption could increase in importance, as compared to the case where a dark photon is absorbed. We reserve this computation for future work as it requires knowledge of anharmonic phonon interactions.
In Tab. III we provide a summary of the target materials that so far have been proposed for sub-MeV dark matter scattering and sub-eV dark matter absorption, and their sensitivity to different models.
Single element semiconductors such as Ge and Si have a similar phonon spectrum as polar materials and therefore have sensitivity to the same models, with the exception of dark photon mediated scattering.
The optical phonons in Ge and Si crystals do not give rise to long-range dipole fields, so dark photon mediated interactions cannot excite a single optical phonon in the long-wavelength limit (multiphonon excitations are still possible, and have been considered for dark photon DM absorption [33]). We conclude that polar materials can test a wide range of models for sub-MeV dark matter, with the added advantage of a directional dependence in the scattering rate for certain materials. That polar materials are readily available and well-understood crystals also makes them an exciting prospect for experimental realization. To calculate the phonon eigenmodes, ν, for a particular crystal, we require solutions to the eigenvalue equation: j D q,j,j · e ν,j ,q = ω 2 ν,q e ν,j,q with the dynamical matrix D q,j,j given by and V (2) l,j,l ,j are the force constants to be calculated, as shown in Section III. In this work we use the frozen-phonon method to calculate the force constants and the corresponding dynamical matrix. This method displaces each atom in the unit cell and calculates the resulting forces on the other atoms using DFT. From a combination of symmetry-inequivalent displacements, the full force-constant matrix can be built up using DFT calculations. A post-processing software package, phonopy [38], is then used to solve the eigenvalue problem for ω ν,q and e ν,j,q .

Computational Details for DFT
Our density functional theory calculations were performed with the projector augmented-wave (PAW) method [76] as implemented in the VASP code [44]. All calculations were performed using the Perdew-Becke-Ernzerhof (PBE) parametrization of the generalized gradient approximation (GGA) [77]. The wavefunctions were expanded using plane waves with an energy cutoff of 600 eV, and used a Monkhorst-Pack [78] k-point sampling mesh of 12x12x4 for 10-atom calculations, 6x6x4 for 30-atom calculations and 4x4x4 for 90-atom (tripled unit cell) calculations. We performed a full relaxation of the lattice constants and internal coordinates of the structure until the forces were converged to 0.01eV/Å. The phonon calculations and modulations of the phonon modes were performed using the frozen-phonon method as implemented in the phonopy [38] software.
The incoming and outgoing DM states can be modeled by plane waves p i − (q + G)| and |p i , such that the transition matrix element p i − q − G|H|p i is The expression for electronic transitions is identical, except that the appropriate in-medium wave functions must be used instead of plane waves. Since N is formally infinite, the matrix element in (B7) appears to go to zero. However, in Fermi's golden rule the squared matrix element is always be evaluated as a sum over q, which diverges as well for N → ∞. We can thus go to continuum limit by where the integral runs over the Brillouin zone. Since we work in the continuum limit for the calculations in Sec. III, it is convenient to absorb the √ N Ω factor directly into the matrix element, which then becomes which manifestly independent of the number of cells in the lattice.
In the isotropic, long-wavelength limit we can drop the dependence on the reciprocal lattice vectors G. Taking a 2-atom unit cell such as for GaAs, the expression reduces to where ω LO is the frequency of the optical phonon and µ ≡ (1/m 1 + 1/m 2 ) −1 is the reduced mass. Here we used that the eigenvectors are normalized within the unit cell (see condition above Eq. 9), so |e j | = 1/ √ 2 for a 2-atom unit cell, and that 1 With the identity (B10) then reproduces Eq. (22).
It now only remains to derive (B11). Following Ref. [55], we consider a harmonic oscillator with reduced mass µ, charge Z * , and natural oscillation frequency ω T O (this will be identified as the frequency of the TO modes, hence the notation). When the oscillator is driven by an electric field with amplitude E 0 and frequency ω, the amplitude of the oscillations is given by The macroscopic polarization vector is P = eZ * u 0 /Ω, where the 1/Ω is merely the number density of the oscillators. The displacement vector of the system is then with the frequency dependent dielectric function: The ∞ term is again the contribution from the valence electrons, while the second term is the contribution from the oscillators. At high frequencies the ions are too slow to respond and only the electron contribution remains. Gauss' law demands that k · D = 0, which is trivially satisfied for the transverse modes. For the longitudinal mode k E, this requires that (ω) = 0, which is satisfied at the frequency ω = ω LO with One may interpret the additional term as the self-energy correction to the LO mode from the backreaction of its induced electric field. Combining (B14) and (B15) yields the Lyddane-Sachs-Teller with 0 ≡ (0). Combining this with (B15) results in (B11).
Appendix C: Nucleon-scattering structure factor In this Appendix, we present the derivation of the dynamic structure factor for DM scattering in a lattice at zero temperature. We follow closely the discussion presented in Ref. [45], which reviews scattering of cold neutrons in a lattice. To compute the structure factor for hard sphere scattering, we treat the crystal as a regular, periodic lattice with N cells and n atoms in a unit cell, for a total of N × n atoms in the lattice. Summing the potential of the individual scattering centers gives the total potential where J sums over all the atoms in the lattice, b X is the DM-nucleon scattering length, and A J is the mass number of the nucleus J. In Fourier space, the potential is We then define the structure function by with λ i,f the initial and final states, and p(λ i ) is the thermal distribution over the initial states. Since we envision a very cold target, we only consider the ground state in the sum of the initial states, setting λ i = λ 0 . We have normalized S(q, ω) such that it is an intrinsic quantity under N → ∞. With this definition, the rate from Fermi's golden rule is where we treated the incoming and outgoing DM particle as plane waves. The integral is over the Brillouin zone and Ω is the volume of the primitive unit cell.
To compute the structure function, first we note that the squared matrix element in (C3) can be rewritten as a single correlation function as follows: Since λ 0,f are eigenstates of the Hamiltonian, we can replace the E λ 0 ,λ f with the operator H: In the last step we used the quantum evolution operator on the phase factor, and made the time dependence of the r J explicit. By applying the sum over the final states from (C3), we can use the completeness of the |λ f states to obtain In what follows we will drop the |λ 0 to facilitate the notation and all expectation values are understood to be with respect to the ground state.
To compute this two-point correlation function, we write the position vectors in terms of the atomic displacements u relative to their equilibrium position, where now l labels the lattice vector for a given primitive cell, and r 0 j are the equilibrium positions of the atoms relative to the origin of the primitive cell. We thus replace the sum over all atoms in lattice (labelled by J) with a sum over all lattice vectors l and atoms in a single primitive cell (labelled by j).
Since the mass numbers A J are identical within each cell, we can also take A J → A j . Inserting this in the correlation function, We wish to expand this in the displacements, and keep only the leading correlation function. We can do so by applying the Baker-Campbell-Hausdorff identity and truncate at leading order. Concretely, for two operators A = iq · u j,l and B = −iq · u j ,l we have Since we are in the small displacement (harmonic) approximation, the operators u j,l can be written as a linear combination of creation and annihilation operators. The commutator in (C13) is therefore proportional to the identity operator and we can pull it outside of the expectation value: In the last step we brought commutator back into the expectation value, again using that it is proportional to the identity as long as A and B are linear combinations of creation and annihilation operators.
where in the second line we expand the exponential to leading order and drop the constant piece that does not contribute to scattering. The two exponentials in front are the Debye-Waller factors, defined by W j (q) ≡ 1 2 (q · u j ) 2 . (C20) where we dropped the l index due to translation invariance over the lattice vectors. From (C20), we see that the Debye-Waller factor measures the average motion of atom j relative to the momentum transfer.
Putting the above results together, the structure function is then S(q, ω) = 1 2πN j,j ,l,l A j A j e iq·(r 0 j −r 0 j ) e iq·(l−l ) e −W j (q) e −W j (q) +∞ −∞ dt q · u j ,l (0)q · u j,l (t) e −iωt . (C21) To further simplify the sums, one can use the invariance of the two point function under lattice translations, which permits the replacement l,l e iq·(l−l ) → N l e iq·l : S(q, ω) = 1 2π j,j ,l A j A j e iq·(r 0 j −r 0 j ) e iq·l e −W j (q) e −W j (q) +∞ −∞ dt q · u j ,0 (0)q · u j,l (t) e −iωt .
It remains to compute the correlation function and the Debye-Waller functions in terms of the phonon eigenvectors and dispersion relations. To this end, we decompose the displacement operators in creation and annihilation operators, as in (7) u j,l (t) = where the index ν runs over all 3n phonon modes and k is a regular, N -point discretization of the first Brillouin zone. The 1/ √ N factor implies that S(q, ω) is an intrinsic quantity, as mentioned below (C3). Inserting this in the two-point function, we can trivially perform the Wick contractions, at least in the zero temperature limit. (For the finite temperature result we refer to section 9.12 of [45].) This results in (q · u j ,0 )(0)(q · u j,l )(t) = 1 2N √ m j m j ν,k 1 ω ν,k (q · e ν,j ,k )(q · e * ν,j,k )e iω ν,k t e −ik·l e ik·(r 0 j −r 0 j ) . (C24) The Debye-Waller function is just the special case where j = j , l = 0 and t = 0: W j (q) = 1 2 (q · u j,0 ) 2 = 1 4N m j ν,k 1 ω ν,k |q · e ν,j,k | 2 (C25) Putting everything back together, we find S(q, ω) = 1 2N j,j ,l A j A j e iq·(r 0 j −r 0 j ) e iq·l e −W j (q) e −W j (q) 1 √ m j m j (C26) × 3n ν k 1 ω ν,k (q · e ν,j ,k )(q · e * ν,j,k )e −ik·l e ik·(r 0 j −r 0 j ) δ(ω ν,k − ω) where we used +∞ −∞ dt e i(ω ν,k −ω)t = 2πδ(ω ν,k − ω). With the identity l e iq·l = N G δ q,G with G the reciprocal lattice vectors, this finally reduces to S(q, ω) = 1 2 G,k,ν 1 ω ν,k |F ν (q, k)| 2 δ k−q,G δ(ω ν,k − ω) (C28) with the phonon form factor F ν (q, k) ≡ j A j √ m j e −W j (q) q · e ν,j,k e i(q−k)·r 0 j (C29) Both energy and crystal momentum conservation are now manifest in these expressions. For scattering with sub-MeV dark matter, the momentum transfer is typically smaller than the size of the Brillouin zone, such that we can neglect the sum over the reciprocal lattice vectors and set G = 0. In this limit, the structure factor further simplifies to Note that this expression differs from the one in [34] by a phase factor, since a different convention was used for the phonon eigenvectors.

Appendix D: DM-electron scattering
In this appendix we comment on the reach for models where the DM couples to electrons through a scalar mediator. Such models tend to be extremely constrained by stellar cooling and Big Bang Nucleosynthesis bounds for m X 1 MeV [51,80], and at the moment we are not aware of models which can achieveσ e 10 −45 cm 2 . While it is likely difficult for near future experiments to access such low cross sections, we briefly discuss the reach for the sake of completeness.
Because the displacements involved in a phonon excitation correspond to displacements of the nucleus and tightly-bound inner shell electrons, a DM-electron coupling also results in an effective DM-phonon coupling, analogous to the discussion in Sec. IV for DM-nucleon couplings. The difference is that we must replace the mass number of the atom in the form factor (30) with the number of core electrons for each atom. Ga and As both have 28 core electrons, while O and Al have 2 and 10 core electrons, respectively. 11 Note that the form factors for coherently scattering off the electrons in Also shown are projections for Dirac materials [28] and superconductors [27].
the atom are constant for |q| 1 keV [81], and we can neglect their effect for the DM mass range of interest. For higher DM masses, these form factors are expected to suppress the rate.
The results are shown in Fig. 15 where y e (y X ) is the electron-mediator (DM-mediator) coupling, µ eX is the DM-electron reduced mass, and α is the fine structure constant. If only the optical branches are accessible, we find a reach that is competitive with that of Dirac material targets, in which the DM can create an electron excitation with ∼ meV threshold. In the optimistic case where the acoustic modes can also be resolved, polar materials could have a reach approaching that of a superconducting target.

Appendix E: Statistical power of daily modulation signal
To estimate the discriminating power of the daily modulation, we calculate how many events are needed to distinguish the scenarios where (i) all observed events are due to a hypothetical, nonmodulating background and (ii) all events are due to a modulating signal, as predicted in Sec. III.
For a given mass point m X and an expected number of events N ev , we generate simulated datasets for both scenarios above. The number of events in each dataset is Poisson distributed with average N ev , and for the modulating sample the probability distribution in t is given by the computations in Sec. III. For each dataset, we then perform a fit to the modulation, allowing for both a constant component and a modulating component with amplitude A, fixing the template for that m X . Denote the modulation amplitude as A non-mod for the datasets that are purely background, and A mod for the datasets that are purely signal. Repeating this procedure for many datasets, we generate the expected distribution in the modulation amplitude, shown in Fig. 16 for an example set of parameters. By construction, A mod = 1, and A nod-mod = 0.
For each N ev , we then compute the 2σ upper value (95% quantile) on A for the non-modulating data (A non-mod 2σ , indicated by the green arrow in Fig. 16). Interpolating in N ev , we can then find the number of events needed such that A non-mod 2σ is below the expected amplitude for the modulating sample, in other words A non-mod 2σ < 1. This gives the number of events needed so that in 50% of the purely signal datasets, we can reject the background hypothesis at 2σ.
Similarly, we can obtain ±σ quantiles about the mean expectation for the modulating signal (A mod ±σ , indicated by the blue arrows in Fig. 16). The ±σ bands are then obtained by demanding that A non-mod 2σ < A mod ±σ . The results of this procedure are shown in Fig. 17 as a function of m X , and translated in terms of cross section in Fig. 5 (blue shaded band).