Multi-Channel Direct Detection of Light Dark Matter: Target Comparison

Direct detection experiments for light dark matter are making enormous leaps in reaching previously unexplored model space. Several recent proposals rely on collective excitations, where the experimental sensitivity is highly dependent on detailed properties of the target material, well beyond just nucleus mass numbers as in conventional searches. It is thus important to optimize the target choice when considering which experiment to build. We carry out a comparative study of target materials across several detection channels, focusing on electron transitions and single (acoustic or optical) phonon excitations in crystals, as well as the traditional nuclear recoils. We compare materials currently in use in nuclear recoil experiments (Si, Ge, NaI, CsI, CaWO$_4$), a few which have been proposed for light dark matter experiments (GaAs, Al$_2$O$_3$, diamond), as well as 16 other promising polar crystals across all detection channels. We find that target- and dark matter model-dependent reach is largely determined by a small number of material parameters: speed of sound, electronic band gap, mass number, Born effective charge, high frequency dielectric constant, and optical phonon energies. We showcase, for each of the two benchmark models, an exemplary material which has a better reach than in any currently proposed experiment.

As new experiments are being planned and detection technologies are being discussed and improved, it is important to identify the most promising targets in order to prioritize the experimental program. There are two questions in this respect: (i) what types of excitations can be utilized as efficient detection paths with current and developing technologies, and (ii) what materials have the strongest response to DM scattering?
It is the purpose of this paper to initiate a discussion on these questions, and provide theory input to the op-timization of experimental strategy. We consider several complementary detection channels: • nuclear recoils, sensitive to the heaviest DM masses, down to O(100 MeV) at best; • electron transitions across band gaps in crystals, covering DM masses down to O(100 keV); • single phonon excitations in crystals, reaching the lightest DM masses, down to O(keV).
The last two detection channels rely on collective properties of the target, which makes calculating the DM model reach more involved than the standard nuclear recoil calculation. While nuclear recoil was proposed long ago [57,58], electron transitions in semiconductors (proposed in Refs. [12,13,16]) and phonon production from sub-MeV DM in crystals (put forth in Refs. [19,46,47]) have a much shorter history. Now that all these ideas are available, we hope to find materials which have a strong response in all channels, in order to cover a broad range of DM masses. We begin in Sec. II with a brief review of each detection channel. A common framework to calculate the reach via all three channels is presented in a companion paper [59], which makes it clear that the detection rate factorizes into the particle-level scattering matrix element squared and a material specific dynamic structure factor that captures the target response. Here we summarize the main results of Ref. [59]. Our goal is to find materials with strong responses (a large dynamic structure factor) in each channel over the kinematically allowed mass region.
Toward this goal, in Secs. III and IV, we carry out a detailed comparison of target materials, focusing on two benchmark DM scenarios to illustrate how to optimize  Table I. Target materials studied in this work and their key parameters. The four blocks contain materials currently in use in nuclear recoil experiments, those considered for proposed near-future experiments, those with superior properties for some specific DM models discussed in this paper, and the remaining ones in alphabetical order, respectively. Sensitivity of electron transitions relies heavily on the band gap Eg, for which experimental values are shown (those with asterisks are measured at low temperature). Nuclear recoils and acoustic phonon excitations in the nucleon-coupling benchmark model are largely determined by the speed of sound of longitudinal acoustic phonons c LA s and atomic mass numbers Aj. For optical phonon excitations in the light dark photon mediated model, relevant parameters are the Born effective charges Z * , high-frequency dielectric constant ε∞, optical phonon energies ωO as well as Aj, all of which combine into a quality factor Q, defined in Eq. (27), which determines the reach at high mass. Barred quantities are properly averaged values; see Appendix A 3 for details.
target choice for the best sensitivity. Our study covers a total of 24 crystal materials, whose key properties that determine sensitivity to DM scattering are summarized in Table I

II. DETECTION CHANNELS
We begin by briefly reviewing the detection channels, which are discussed thoroughly in our companion paper [59]. Generally, for a DM particle χ, the event rate per unit target mass is given by where ρ T is the target mass density, ρ χ is the local DM energy density, m χ is the DM mass, and f χ (v) is the incoming DM's velocity distribution in the target rest frame. The event rate Γ(v) for an incoming DM particle with velocity v is usually normalized against a reference cross section, defined from the particle-level scattering matrix element M (in the nonrelativistic normalization) evaluated at a reference momentum transfer q 0 . Here we adopt the following definitions, for DM-nucleon and DM-electron interactions, respectively, where µ χn , µ χe are the reduced masses, and v 0 is the dispersion of the DM's velocity distribution. They coincide with the total particle-level scattering cross sections in the case of a heavy mediator. As we show in Ref. [59], for spin-independent (SI) scattering off a target material via tree-level exchange of a mediator, the matrix element factorizes into a DM component that is universal, and a target response component captured by a dynamic structure factor S(q, ω) that is target and excitation specific, such that Here σ, µ represent either σ n , µ χn or σ e , µ χe , q is the momentum transfer from the DM to the target, and is the corresponding energy deposition. The mediator form factor is given by 1 The dynamic structure factor, which captures the target's response to a general energy-momentum transfer ω, q, is given by where V is the total volume, |i , |f are the initial and final states of the target system, and F T is the quantum mechanical operator acting on the target Hilbert space that the DM couples to.
For an isotropic target, the dynamic structure factor depends only on the magnitude but not the direction of q, so the velocity integral can be evaluated independently, giving for which analytic expressions can be obtained assuming a boosted truncated Maxwell-Boltzmann (MB) distribution. On the other hand, for the more general case of anisotropic target response, the dynamic structure factor depends on the direction of q, and we can utilize the delta function in Eq. (7) to evaluate the velocity integral first, giving which can be computed analytically for the usually assumed boosted truncated MB distribution. In the following subsections, we consider each detection channel in turn, summarizing the formalism presented in Ref. [59] on the dynamic structure factors and detection rates, building on the discussion in previous works (particularly [16,47,78]).

A. Nuclear Recoils
For each nucleus species, where m N is the nucleus mass, f n , f p and f N = f p Z + f n (A − Z) are the DM-neutron, DM-proton and DMnucleus couplings respectively, and F N (q) is the Helm form factor r n 1.14 A 1/3 n fm , s 0.9 fm , which approaches 1 in the q → 0 limit. The differential rate with respect to energy deposition, generalized to the case of multiple nucleus species, is where v min = q 2µ χN . The conventional nuclear recoil calculation is valid when each nucleus can be considered independent of the other nuclei. In a crystal target, this is true if the scattering happens at a timescale 1/ω much shorter than the inverse phonon frequencies 1/ω ph , i.e. if the energy deposition ω ω ph ∼ O(100 meV), or equivalently, q √ m N ω ph (note that this momentum cutoff is essentially the inverse of the spatial extent of nucleus wavefunctions in a harmonic potential). For lower energy depositions, the scattering event proceeds by direct production of (single or multiple) phonons. We discuss single phonon excitations in Sec. II C. We will see that single phonon excitation rates are suppressed by the Debye-Waller factor for q √ m N ω ph , which shows the complementarity between the two channels.

B. Electron Transitions
In solids, electrons form band structures with energy eigenstates labeled by a band index i and a wave vector k within the first Brillouin zone (1BZ). In an insulator or semiconductor, all electrons occupy the valence bands at low temperatures, and can be excited across the band gap to conduction bands. The dynamic structure factor encapsulates all such transitions from i 1 , k 1 to i 2 , k 2 : up to screening effects. Here G = n 1 b 1 + n 2 b 2 + n 3 b 3 , with n 1 , n 2 , n 3 ∈ Z and b 1,2,3 are reciprocal primitive vectors. The crystal form factor is defined by where u i (k+G) are Bloch wavefunction coefficients computed from DFT (see Appendix A 1). We neglect possible spin dependence of the electron band structures, and simply sum over contributions from the degenerate spin states. The total rate is given by where q = k 2 − k 1 + G and ω = E i2,k2 − E i1,k1 are assumed. Note that unlike in nuclear recoils, the dynamic structure factor for electron transitions is generally not isotropic in q for all energy-momentum depositions. When anisotropies are significant, the rate cannot be expressed in terms of η(v min ), and the g function in Eq. (10) should be used instead. The physical implication is that the rate depends on the direction of the DM wind and exhibits daily modulation. An example of this is discussed in Ref. [59].

C. Single Phonon Excitations
Phonons are quanta of lattice vibrations in crystals. For a three-dimensional crystal with n atoms/ions in the primitive cell, there are 3n phonon branches, with dispersions ω ν,k (ν = 1, . . . , 3n), where the wave vector k is in the 1BZ. The dynamic structure factor has the general form where Ω is the volume of the primitive cell, j = 1, . . . , n runs over the atoms/ions in the primitive cell, x 0 j are their equilibrium positions, and m j are their masses. Y j contains the DM-atom/ion couplings, whose general definition is given in Ref. [59]. We explicitly state the expression of Y j for each benchmark model below. ν,k,j are the phonon polarization vectors. k is the momentum within the 1BZ that satisfies q = k + G for some reciprocal lattice vector G -only those phonon modes that match the momentum transfer up to reciprocal lattice vectors can be excited, as a result of lattice momentum conservation. At large q, the dynamic structure factor is suppressed by the Debye-Waller factor, given by We obtain the total rate where m cell = ρ T Ω is the mass contained in a primitive cell. The phonon dispersions ω ν,k and polarization vectors ν,k,j that enter this equation are obtained from DFT calculations (see Appendix A 2).

Acoustic vs. Optical Phonons
It is useful to distinguish acoustic and optical phonons, as they are sensitive to different types of DM interactions. Among the 3n phonon branches, three are gapless with linear dispersions ω ν,k ∼ c s |k| near |k| = 0 (with c s the sound speed), as a result of spontaneous breaking of translation symmetries; these are acoustic phonons that, in the long wavelength limit, correspond to in-phase oscillations of atoms/ions in the same primitive cell. The remaining 3(n−1) branches are gapped "optical" phonons, corresponding to out-of-phase oscillations.
Due to the nature of in-phase oscillations, acoustic phonons can be efficiently excited if DM couples to different atoms/ions in a correlated way. An example is a DM particle coupling to nucleons via a scalar or vector mediator. In this case, Y j is proportional to a linear combination of A j and Z j , and can have the same sign and similar magnitudes for all j.
By contrast, the out-of-phase oscillations associated with gapped phonon modes have enhanced sensitivity to DM coupling to the atoms/ions in the same primitive cell differently. This is the case for dark-photonmediated DM scattering with polar materials. The dark photon mediator kinetically mixes with the SM photon, and as a result, Y j point in opposite directions for oppositely charged ions. We follow convention and call all gapped phonon modes "optical," though only in polar materials where there are both positively and negatively charged ions in the primitive cell (e.g. GaAs) do these modes couple strongly to the (dark) photon via the oscillating dipole. Diamond, Si and Ge, for example, all have gapped phonon modes, but none of these materials has a strong coupling to the dark photon as the primitive cell does not contain oppositely charged ions.

III. TARGET COMPARISON: KINETICALLY MIXED LIGHT DARK PHOTON MEDIATOR
A well motivated model of light dark matter involves interaction with the SM via a light dark photon A that kinetically mixes with the photon: where D µ = ∂ µ − ie A µ , and the DM χ can be either a complex scalar or a Dirac fermion. The gauge boson kinetic terms can be diagonalized by redefining A µ → A µ +κA µ , which gives J µ EM a charge under the dark U (1) of κe. The reference cross section, utilized in present results for this model, is given by The projected 95% condifence level (C.L.) exclusion reach on σ e assuming zero background (i.e., the cross section needed to obtain three events) from electron transitions and single phonon excitations are shown in Fig. 1, for m A → 0 and an exposure of one kg-yr. 2 In the rest of this section, we describe in detail the features in this plot, and also discuss nuclear recoils.

A. Single Phonon Excitations
Optical phonon excitation is the dominant detection mode for dark photon mediated scattering. As shown in Ref. [59], in the low q limit (which dominates the momentum integral for a light mediator since F 2 med ∝ q −4 ), the interaction is described via the Born effective charges of the ions, Z * j (which are generally 3 × 3 matrices), where ε ∞ is the high-frequency dielectric matrix. The total rate is given by Eq. (20). Only polar materials, or those which have differently charged ions in the primitive cell, can couple phonon modes to the dark photon, which explains the absence of phonon reach curves for Si and Ge in Fig. 1.
As explained in the previous section, optical phonon modes involve out-of-phase oscillations and are gapped. Because the optical modes are the dominant contribution to the rate, the properties of the optical modes determine the shape of the phonon excitation curves in Fig. 1: when there are sharp changes in the reach as a function of mass, it is because there is a transition in the dominance of a particular optical mode. For low momentum transfer, the dispersion of the gapped modes is approximately a constant, such that the lowest DM mass reachable is determined by setting the maximum kinetic energy of the incoming DM, m χ v 2 max /2, equal to the energy of the low- Projected reach from single phonon excitations (dashed) and electron transitions (solid) for DM scattering mediated by a kinetically mixed light dark photon (the smallest-gap target InSb suffers from slow convergence in the electronic transition calculation at mχ < 1 MeV, for which we show results of the two most accurate runs with solid and dotted curves, see Appendix A 1 for details). Nuclear recoils (not shown) can also probe this model, but the conclusion on which targets are superior is the same as for the light hadrophilic mediator model. A detector threshold of 1 meV is used for the phonon calculations, and all transitions with energy deposition greater than the band gaps are included in electron excitations. The freeze-in benchmark is taken from Refs. [12,80], corrected by including plasmon decay for sub-MeV DM [81]. Stellar constraints are from Ref. [82] and direct detection constraints are from DAMIC [61], DarkSide-50 [83], SENSEI [62], SuperCDMS [68], XENON10 [14,21], and XENON100 [83,84] Thus materials having low energy optical phonon modes are desirable to search for light dark matter; CsI, for example, has particularly low-lying optical phonon excitations, and its sensitivity to the lightest DM masses is seen in Fig. 1.
We can also see that at higher masses, single optical phonon production rates vary widely between materials. This can be understood analytically. Consider first the simplest case of a diatomic polar crystal (e.g. GaAs). The dominant contribution to the q integral in Eq. (20) is well within the 1BZ and therefore we can set G = 0, W j 0, and g(q, ω) ∝ q −1 . Approximating Z * j Z * j 1, and noting that Z * 1 = −Z * 2 ≡ Z * , we see that the rate is dominated by the longitudinal optical (LO) mode, for which one can show LO,k,1 and LO,k,2 are anti-parallel, and | LO,k,j | = µ 12 /m j in the limit k → 0, where µ 12 is the reduced mass of the two ions. Further approximating the phonon dispersion as constant and ε ∞ ε ∞ 1, the rate simplifies to We call Q a quality factor, since it is the combination of material-specific quantities that determines the direct detection rate. A higher-Q material has a better reach in the high mass regime. More concretely, we find Note that although we have focused on the special case of diatomic polar crystals in order to derive analytic estimates, similar considerations apply for more complicated crystals. For example, it is not surprising that larger Born effective charges and lighter ions are helpful. When comparing the targets, we adopt the following prescription for the quality factor, where n is the total number of ions in the primitive cell, and ω O is the directionally averaged optical phonon energy of the highest mode near k = 0, given in Table I. In our list of materials LiF has the largest quality factor, with SiO 2 second. We choose to highlight SiO 2 in Fig. 1 because LiF is a less desirable experimental target due to large backgrounds [85]. A further consideration for optimizing Q given a fixed chemistry (atomic species) is to maximize the Born effective charges. For example, cubic tungsten trioxide (WO 3 ) has been reported to have anomalously high Born effective charges of up to +12.5 and −9.1 on W and O respectively [86]. Materials with such high Born effective charges, a manifestation of highly covalent bonding character, provide a further route for maximizing Q. 4 We comment in passing that also in the case of a heavy dark photon mediator, the rate is largely determined by the quality factor defined in Eq. (27) for sub-MeV DM; for heavier DM, couplings to ions cannot be simply captured by the Born effective charges at high momentum transfer, and the total rate is more challenging to compute [59].

B. Electron Transitions
The typical band gaps between valence and conduction bands, E g , range from a fraction of an eV (InSb and Ge) to as high as 10 eV (e.g. SiO 2 ). This gap sets the lightest DM mass to which the experiment is sensitive, as kinematics requires that 4 Cubic WO 3 is dynamically unstable giving imaginary frequencies in the phonon band structure. Therefore we do not include it in phonon comparison plots, and leave a study of other stable isomorphs for future work.
Thus, small gap materials will generally have better reach. For example, InSb is superior to Si for m χ MeV, as seen in Fig. 1; in fact, the sub-eV band gap of InSb allows for a significant G = 0 contribution that is absent for larger gap materials, and this contribution dominates at m χ MeV, greatly extending the reach. However, note that Ge, which has a smaller band gap than Si, does not have a better reach. The difference here is due to a direct vs indirect band gap. 5 When depositing energy via a scattering process, there must be some momentum transfer, and therefore, strictly speaking, E g in Eq. (28) should be replaced by the minimum kinematically allowed energy difference. For direct gap materials this means that m χ,min will increase, as it does in Ge, which is why Ge has worse reach than Si. Note that there is a complementarity between single phonon excitations and electron transitions. In the phonon case, materials with the best sensitivity tend to be insulators, as they have small values of ε ∞ . However, for electron transitions, one prefers materials with smaller band gaps, which generally have larger values of ε ∞ . This is because loop corrections to the in-medium photon propagator are larger for a smaller band-gap: virtual electrons can be more easily created because of the smaller energy difference.
For higher masses an analytic comparison is not tractable. The wavefunction coefficients in Eq. (17) cannot be modeled well analytically, and hence the reach must be computed numerically. Note that for Si, Ge, NaI, CsI, GaAs, and diamond, our results are roughly consistent with previous calculations in Refs. [16,18,23], where the DFT calculation is implemented differently. However, we find discrepancies in the semi-core electron contributions, which are subdominant for our light mediator benchmark, but become important for a heavy mediator. We will investigate this issue in detail in an upcoming publication. Another improvement of the calculation that we plan to address is the treatment of inmedium screening effects (see Ref. [59] for further discussion), which we have neglected in the present calculation. Such effects are expected to be weak for materials with band gaps larger than about 1 eV. However, for sub-eV gap targets such as InSb, for masses below ∼ 1 MeV, the result here should be taken with caution, as the effects may not be negligible.

C. Nuclear Recoils
The dark photon mediator coupling in a target system is momentum dependent. At very small momentum 5 The HSE06 exchange-correlation functional used in our DFT calculations slightly underestimates the direct band gap of Ge whilst being a close match to the indirect band gap [87]. This leads to the prediction of a direct band gap when optimized lattice parameters are used, contrary to experiment.
transfers q → 0, the coupling is negligible as the total target is assumed to have no net charge. For q r −1 ion , where r ion is the size of an atom without the binding electrons, ionic charges, if present, can be coupled to. As the momentum transfer increases further, outer-shell electrons will respond incoherently, possibly transitioning to conduction bands independent of proton and innershell electron responses. On the other hand, in a nuclear recoil event, q √ m N ω ph r −1 ion . In this regime, protons respond coherently as long as F N (q) 1, since they are bound in the nucleus, whereas electron couplings are irrelevant since even the core electron wavefunctions do not have such high momentum components. Therefore, nuclear recoils can happen in an overall neutral crystal via coupling to the proton number of each nucleus without any atomic form factor suppression.
In order to compare against phonon and electron excitations, we express the reach in terms of σ e instead of σ n . This corresponds to replacing (f N /f n ) 2 → Z 2 N for each nucleus species, and µ χn → µ χe , q 0 → αm e , and lastly, σ n → σ e in Eq. (14). While we discuss material comparison in this subsection, nuclear recoil reach curves have been omitted in Fig. 1 in order not to further complicate the plot; they can be approximately rescaled from the reach curves in Fig. 2 below, and are straightforward to compute from Eq. (14). 6 The low mass reach of nuclear recoils is material and threshold dependent, and can be understood from kinematics. The maximum momentum transfer is given by q max = 2µ χN v max , and therefore the maximum energy deposited is given by ω max = 2µ 2 χN v 2 max /m N . Requiring that this be larger than the threshold sets the minimum DM mass. For a threshold around 500 meV (which almost saturates the validity bound for some of the crystal targets as discussed in Sec. II A), and v max = 10 −3 , m χ m N , the minimum DM mass within reach is m χ,min ∼ 100 MeV ω min 500 meV Therefore, materials with lighter nuclei are more favorable for kinematic matching.
At higher masses, kinematics is not a limiting factor, and we can obtain an analytic approximation for the rate. Assuming a singular nuclear species, A N = A, Z N = Z simplifies Eq. (14) to and we see that the rate is dominated by small ω. At masses above a few hundred MeV and small ω, η(v min ) 6 Changes in the constraint projections via single phonon excitations in Figs approaches η(0). The total rate then becomes and is approximately material independent. Note that if the dark photon mediator is heavy, the factor A 2 ω 2 in the denominator of Eq. (30) would be absent, and heavier (larger Z) elements are advantageous.

IV. TARGET COMPARISON: HADROPHILIC SCALAR MEDIATOR
As a second benchmark model, we consider a real scalar mediator φ coupling to the proton and neutron, where the DM χ is taken to be either a real scalar or a Dirac fermion. In the absence of electron couplings, the relevant search channels are single phonon excitations and nuclear recoils. We will quote the reach in terms of σ n , given by The 95% C.L. exclusion reach on σ n for a light (effectively massless) and heavy mediator are shown in Figs. 2 and 3 respectively, assuming f p = f n , an exposure of one kg-yr, and zero background events. In the rest of this section we explain in detail the features in these plots.

A. Single Phonon Excitations
We first consider DM creating a single phonon via the nucleon coupling. As shown in Ref. [59], where f j = f p Z j + f n (A j − Z j ) for the nucleus at site j in a primitive cell, and F Nj (q) is the nuclear form factor given by Eq. (13). As before, the total rate is calculated from Eq. (20). However, a major difference compared to the dark photon mediator model is that, if f p and f n have the same sign, the rate is dominated by acoustic and not optical phonons, assuming the energy threshold is low enough to access the acoustic phonons. This is because Y j points in the same direction for all j, resulting in stronger in-phase oscillations as discussed in Sec. II C 1. We first discuss Fig. 2, for the light mediator case,   . Single phonon and nuclear recoil reach for a light (m φ = 1 eV) hadrophilic scalar mediator. 1, 20, and 100 meV thresholds are shown for the single phonon reach (solid, dashed, and dotted respectively), and 500 meV threshold is assumed for the nuclear recoil reach (dot-dashed). For m φ = 1 eV the dominant constraint on fn is from fifth force experiments [88]. If mχ makes up all the DM then the dominant constraint on yχ is from DM self-interactions (SIDM) [88]. If mχ is only a subcomponent, we only require perturbativity yχ < 1 (Pert.); in this case the reach curves can be easily rescaled. 6 when the energy threshold, ω min , is 1 meV. While such a low threshold is experimentally challenging, the curves are easier to understand conceptually compared to the higher ω min curves. In fact, over most of the mass range, for most materials, the rate is dominated by single longitudinal acoustic (LA) phonon production. At the high mass end, the reach is material-independent, understood analytically as follows. The mediator form factor F med ∝ 1/q 2 , and therefore the rate is dominated by the lowest detectable momentum transfer. In this case, we can set G = 0 (or equivalently, q = k in Eq. (20)), W j 0, ω LA = c LA s q, F Nj 1, and g(q, ω) ∝ q −1 . Lastly, in this limit q · LA,j,k q m j /m cell . Thus the rate For f p = f n , we have f j ∝ A j ∝ m cell /m n , and the dependence on the target properties drops out. The reference cross section σ n corresponding to a given event rate R scales with mass as µ 2 χn /m 3 χ , as we see in Fig. 2. Note that as we go to higher m χ the reach on the couplings f 2 n y 2 χ gets worse as µ 2 χn m χ ; the apparent better reach at higher mass in Fig. 2 is due to the definition of σ n ∝ m −4 χ . For DM masses below ∼ 0.1 MeV in Fig. 2, kinematics causes the reach to diminish: the maximum momentum transfer, 2m χ v max , must be large enough to reach the minimum momentum transfer set by the detector threshold, ω min /c LA s . This sets the minimum reachable DM mass m χ,min ∼ 20 keV ω min meV To reach the lightest dark matter particle at low thresholds, an ideal material is then diamond, as it has the highest speed of sound. AlN and SiO 2 are the next best candidates from our search. As we move on to the curves with higher energy thresholds, ω min = 20 meV and 100 meV, the materials with lower sound speed lose reach altogether. (The ω min = 500 meV curves are derived from nuclear recoil; this is discussed in the next subsection.) The reason is that acoustic phonons are accessible only when ω min c LA s /a, where a is the lattice spacing. For materials with lower sound speed, the energy threshold may simply never be low enough to have any reach with an acoustic phonon. In addition, one can see where optical phonons start to play a role, as the slope of the reach curve changes at  There are no stellar constraints for m φ 400 MeV [88]. Currently, the best experimental nuclear recoil constraints in this region of parameter space are from DarkSide-50 [89] (assuming binomial fluctuations), and XENON1T (combined limits from [90,91]). We also show the constraint from CRESST-II [77], which is stronger than the DarkSide-50 constraint at low masses assuming no fluctuation in energy quenching. A more complete collection of nuclear recoil constraints can be found in Refs. [89,91,92]. The neutrino floor is taken from Ref. [93]. 6 lower masses, e.g. Si with an energy threshold of 20 meV. This feature will be present for all materials if the lowest kinematically reachable DM mass from optical phonon excitations, given in Eq. (24), is smaller than the lowest kinematically reachable DM mass from acoustic phonon excitations, given in Eq. (36). Next we turn our attention to Fig. 3, for the same hadrophilic scalar mediator benchmark, but with a heavy mediator. Again, we first focus on the case of a 1 meV threshold, as here the acoustic phonon contributions dominate and analytic simplifications can be made since the integrals are dominated by the high momentum behavior. There are four distinct regions in mass and we now discuss the mass and material parameters dependence of each of them.
In the lowest mass regime, m χ 10 −1 MeV, the reach ends when the acoustic modes are no longer kinematically available, just as in the massless mediator case, with minimum reachable mass again set by Eq. (36). Between 10 −1 and 1 MeV, the reach curves flatten and the order of the curves reverses: materials with a higher speed of sound have worse reach, which can be understood analytically starting with Eq. (20). For m χ 1 MeV the momentum transfer is within the 1BZ, so we can take q = k, W j 0, ω = c s q and g(q, ω) ∝ 1/q as in the light mediator case. For simplicity we ignore angular dependence, assume the ions are the same, A j ≡ A, m j ≡ m, set f n = f p , and consider only the longitudinal mode so that q · ε ∝ q. Then we have where the upper cutoff is due to kinematics and manifests in the g function, which goes to zero as k reaches the maximum allowed momentum transfer.
A similar derivation applies to the mass dependence in the next two regimes. For 1 MeV m X 10 MeV, the dominant momentum transfer is outside of the 1BZ, which means that ω can no longer be approximated by c s q. In fact, since ω is only a function of the phonon momentum in the 1BZ, it will vary rapidly as q increases. We therefore exchange ω with a q independent quantity, roughly thought of as the average of ω over the whole 1BZ, ω . The rate becomes Since the rate scales inversely with ω , materials with lower energy phonon modes are preferred. As ω is usually correlated with c s , the ordering of the curves is the same as in the previous regime. We have neglected the Debye-Waller factor in the analytic estimates above, because the momentum transfer is on the order of m χ v, and is less than the Debye-Waller cut-off around m N ω . However, for the last mass regime, above ∼ 10 MeV, this is no longer the case, and the momentum integral is cutoff by the Debye-Waller factor, Therefore, materials with heavier elements and higher phonon energies are preferred. In our search, CaWO 4 has the highest factor of A ω , with PbTe following, which is the reason we choose to highlight PbTe in Fig. 3. For higher thresholds, the optical phonon modes contribute to a greater degree, so the scaling arguments given above for the first two mass regimes no longer hold, but for the last two they do, which is why the curves are almost parallel.

B. Nuclear Recoils
For DM heavier than O(100 MeV), nuclear recoils offer a complementary detection channel to phonon excitations. The low mass behavior of the reach curves is understood in the same way as in Sec. III C (see Eq. (29)), and lighter elements are advantageous. At higher masses, the σ n reach depends on the mediator mass. To show this analytically we again consider a single nucleus species, A N = A, and f n = f p . In the case of a light mediator the differential rate in Eq. (14) becomes For DM heavier than a few hundred MeV, the m N dependence via η(v min ) is weak, as in the dark photon mediator case. The rate is then which is material independent. This is why all the reach curves coincide for large DM masses. We also see that as in the case of acoustic phonons, achieving lower energy thresholds is crucial for improving the reach. If the mediator is heavy, we have where for simplicity we take the η function to decrease sharply at the kinematic bound. We reach the conclusion that heavier nuclei are preferred, similar to the case of single phonon excitations with a heavy mediator. Note also that there is no threshold dependence for larger masses. Therefore a lower threshold only helps to reach lower DM masses, as opposed to the case of the light mediator.

V. CONCLUSIONS
We considered spin independent DM direct detection through three channels -single phonon excitations, electron transitions, and nuclear recoils -in a wide variety of crystal target materials, and two well motivated DM models. Many of these materials are already being discussed for DM detection, but we have presented some new targets for consideration.
For each type of interaction, we specified the target material parameters which should be optimized in order to maximize the reach, and we found complementarity between targets depending on (i) the experimental threshold, (ii) the mass range, and (iii) the model. The experimental threshold dictates which modes are available: at higher recoil energies, only electron transitions and nuclear recoils are possible; as the threshold drops, optical and acoustic phonons become accessible. The phonon modes in materials with high sound speed become kinematically available at higher thresholds than in materials with lower sound speeds. Also, for a given threshold, materials with higher sound speeds have reach to lighter dark matter. Regarding the mass range, the smallest detectable masses are always set by a kinematic constraint, and the dependence on material parameters, and detection threshold, can be found in Eqs. (24), (28), (36), (29) for optical phonon, electron, acoustic phonon excitations, and nuclear recoils respectively. As for the model, we defined a quality factor (in Eq. (27)) for single optical phonon excitations from dark photon mediated scattering to indicate which targets will have the best sensitivity. On the other hand, for a hadrophilic mediator, target optimization for acoustic phonon excitations depends on the mediator and DM masses. We summarize our results in Table II.
An attractive feature of phonon and electron excitations is the possible daily modulation of event rates, as the dynamic structure factors in Eqs. (15) and (18) are generically anisotropic. In the context of phonon excitations, Al 2 O 3 has been considered in Ref. [47], and in our Light dark photon mediator (Sec. III, Fig. 1  We obtain the materials-specific responses using first principles calculations based on density functional theory (DFT) [94]. DFT is a standard method for obtaining solutions to the many-electron interaction problem, and can accurately predict materials properties ab initio ranging from electronic and magnetic to mechanical and vibrational properties. For this work, we used DFT to calculate the full electronic and phonon spectra for a range of materials, with the calculation details given below. However, since DFT is a ground-state method, it suffers from the famous 'band gap' problem where excited-state properties, including band gaps, are not accurately treated using standard DFT methods. We correct for this in two ways: (i) we performed beyond-DFT calculations (hybrid functional calculations) for several of the compounds where standard DFT gave a zero band gap, and (ii) we adjusted the band gaps to experimentally-reported values for all compounds. We note that the convergence parameters used for the electronic and phonon calculations are different owing to the different physical properties being calculated.
The list of materials calculated with their corresponding space groups and space group numbers is given in Table III, with the crystal structures depicted in Fig. 4. For compounds where several structural isomorphs exist, we considered the reported low-temperature ground state structure. The Brillouin zones for the crystal structures considered in this work are depicted in Fig. 5 with the high-symmetry points labelled. Both the electronic and phonon band structure plots take paths through these high-symmetry points.
All DFT calculations were performed using the Vienna Ab initio Simulation Package (VASP) [95][96][97] with projector augmented wave (PAW) pseudopotentials [98,99] using the Perdew-Becke-Ernzerhof (PBE) exchangecorrelation functional [100]. In the PAW scheme, we treated s and p electrons as valence for Li, C, N, O, F, Na, Al, Si, S, Cl, Ca, I, Cs and W, p electrons as valence for Mg and d electrons as valence for Zn, Ga, Ge, As, In and Sb. Below we summarize the convergence criteria used for the (i) electronic structure and wavefunctions, and (ii) phonon calculations.

Calculation details for electronic band structures and wavefunctions
For structural optimizations, a plane wave cutoffenergy of 950 eV was used with a 12 × 12 × 12 Γ-centered k-point grid. The energy and force convergence criteria were 1 × 10 −8 eV and 1 meV A −1 respectively.
All-electron wavefunction coefficients were extracted from PAW calculations using a modification of the pawpyseed code [133]. This enabled recovery of the full wavefunctions as normalized single-particle Kohn-Sham states from the pseudo-wavefunctions obtained by the PAWmethod. Initial PAW wavefunctions were calculated with a plane wave energy-cutoff of 1000 eV, from which the allelectron wavefunctions were constructed with a minimum energy-cutoff of 450 eV. Calculations were performed using Γ-centered Monkhorst-Pack grids, with a k-point density of at least 0.27 A −1 . Energy bands were included up to 60 eV above and below the valence band maximum. However, since there is no pseudopotential containing the low-lying 4d -states for indium in VASP, these bands are neglected from the calculations. In NaI and CsI the I 4d states are positioned at approximately 43 eV and 42 eV below the valence band maxima respectively. A scissor operator was applied to match the experimental band gaps given in Table I. For Ge, InSb and GaSb, the PBE functional gave partially occupied bands due to underestimation of the band gap. In these cases the HSE06 hybrid functional [134] was applied in a static calculation to introduce a band gap before applying the scissor correction. Electronic band structures were computed on a discrete k-mesh along the high-symmetry directions. PbS, PbSe and PbTe were excluded from the electron calculations because spin-orbit interactions are required to capture important features of the band structures and spin-orbit coupling is not yet implemented within the pawpyseed code.
Multiple k-point densities, energy-cutoffs, and energy bands included were tested for all materials to ensure convergence of the scattering rates to less than 2% at m χ = 10 GeV, less than 3% at m χ = 10 MeV, and less than 28%, 18%, 18% and 10% for GaAs, GaSb, Ge, and Si at m χ = 1 MeV respectively. InSb was tested with a 12 × 12 × 12 and 14 × 14 × 14 k-point grid, plotted as dotted and solid curves in Fig. 1 respectively. At m χ = 1 MeV the rate convergence is 5%, and decreases for larger masses. However at smaller masses the G = 0 contribution, from momentum transfers within the 1BZ, and energy depositions below ∼ 1 eV, dominate the rate. The slow convergence here is due to the fact that InSb has rapidly changing band structure near the Γ point, and more k-points are needed for better convergence. These uncertainties are plotted as shaded bands in Fig. 1, and accompanying figures in Appendix B, although most are invisible due to the plots being log-log.

Calculation details for phonon spectra
We obtained the phonon dispersions from Phonopy [135] using the 'frozen-phonon' method by diagonalizing the force matrix using VASP as the force calculator. For the VASP calculations, the electronic wavefunctions were expanded in a plane-wave basis with a kinetic-energy cutoff of 600 eV. The Brillouin zone sampling was no less than 0.8 A −1 in each direction of the unit cell with Monkhorst-Pack grids, and was correspondingly scaled for phonon supercell calculations. Born effective charges were calculated for polar materials using density-functional perturbation theory as implemented in VASP. Table I The experimental electronic band gaps, E g , are taken from references cited in Table III. The speed of sound, c LA s , was calculated by averaging ω LA /q over a uniform 20 × 20 grid on the surface of a sphere in reciprocal space with radius q ≈ 10 eV centered at the Γ point. The same averaging procedure was used in calculating the range of optical modes, ω O , when a range exists. The average Born effective charge, Z * , is defined as Tr Z * + /3, where Z * + is the Born effective charge of the positive ion(s) (the other charges can be found by requiring that the primitive cell is neutral). The average high frequency dielectric constant, ε ∞ , is defined as Tr [ε ∞ ] /3.

Appendix B: Additional Target Comparison Plots
In this Appendix, we provide plots for the remainder of the materials in Table I not presented in the main text. For concreteness, in all figures we take the local DM density to be ρ χ = 0.4 GeV/cm 3 , and assume a Maxwell-Boltzmann distribution with velocity dispersion v 0 = 230 km/s, truncated at the escape velocity v esc = 600 km/s, and boosted to the target rest frame by the Earth velocity in the galactic rest frame v E = 240 km/s. We take the direction of the Earth's velocity to be in thê z direction with respect to the crystal coordinates when computing the reach (for most of the target materials we consider here we expect modulation effects from the Earth's motion to be small). The constraints on σ correspond to a 95% confidence level (C.L.) assuming Poisson distributed counts and no events are seen (equivalently, the constraint corresponds to the cross section needed to obtain three events). We chose the 95% C.L. for easier comparison with previous literature. Because it is also standard to compute the 90% C.L. exclusion reach, we note that one simply has to mulitply the 95% C.L. exclusion reach by 2.3/3 (as the 90% C.L. constraint corresponds to the cross section needed to obtain 2.3 events). We also assume an exposure of one kg-yr.
(a) Diamond: diamond-C, Si, Ge. Two interpenetrating face centered cubic lattices, one offset by 1/4 along the cubic diagonal. Each atom has four nearest neighbors, forming cornersharing tetrahedra.
(b) Zincblende: ZnS, GaAs, InSb, GaSb. Same arrangement as diamond cubic, but with two atom types, each occupying one of the face centered cubic lattices.
(c) Rock salt: NaCl, MgO, LiF, NaF, NaI, PbS, PbSe, PbTe. The two atom types each form a face centered cubic lattice, offset by 1/2 along the cubic axis. One atom type is octahedrally coordinated to the other atom type and vice versa. (e) CsI. The two atom types form interpenetrating primitive cubic lattices, with an atom of one type at the center of each cube of the other type.
(g) Corundum: Al 2 O 3 . Each Al ion is bonded to six O ions, forming octahedra with a mixture of corner, edge and face-sharing connectivities.
(h) Rutile: MgF 2 . Each Mg ion is bonded to six F ions, forming octahedra with a mixture of corner and edge-sharing connectivities.
(i) Wurtzite: GaN, AlN, ZnO. One atom type is tetrahedrally bonded to the other atom type and vice versa.
The tetrahedra are corner-sharing and the structure is a member of the hexagonal crystal system.
(j) CaWO 4 . Each Ca ion is bonded to eight O ions, and each W ion is bonded to four O ions, forming cornersharing octahedra.  Table I.  Table I.  Table I.  Figure 8. Phonon dispersions calculated with VASP and phonopy [135] including non-analytic corrections. The path through the high symmetry points is found using SeeK-path [136].  Figure 9. Phonon dispersions calculated with VASP and phonopy [135] including non-analytic corrections. The path through the high symmetry points is found using SeeK-path [136].