Dark matter-electron scattering in dielectrics

A number of direct detection experiments are searching for electron excitations created by scattering of sub-GeV dark matter. We present an alternate formulation of dark matter-electron scattering in terms of the dielectric response of a material. For dark matter which couples to electrons, this approach automatically accounts for in-medium screening effects, which were not included in previous rate calculations for semiconductor targets. We show that the screening effects appear for both scalar and vector mediators. The result is a non-negligible reduction of reach for direct detection experiments which use dielectric materials as targets. We also explore different determinations of the dielectric response, including first-principles density functional theory (DFT) calculations and a data-driven analytic approximation using a Mermin oscillator model.


I. INTRODUCTION
An increasingly diverse set of underground direct detection experiments remains the most promising direct probe of the nature of the dark matter (DM). In a set of pioneering papers, Essig et. al. [1][2][3] generalized the search for DM beyond the traditional nuclear recoil paradigm, by showing that it is feasible to search for electron recoils from scattering of sub-GeV dark matter in noble liquids and semiconductors. There have since been numerous studies of DM scattering in these materials [4][5][6][7][8][9][10][11][12][13][14][15], as well as proposals for other targets that are sensitive to electron recoils [16][17][18][19][20][21][22][23][24][25][26][27][28]. Nowadays electron recoils are leveraged by every major experimental collaboration, and they are or will be a primary detection channel for experiments such SENSEI [29], DAMIC [30], Su-perCDMS [31] and LBECA [32]. For semiconductor targets in particular, the full calculation of the DM-electron scattering rate was first performed in [3] through an explicit calculation of the electronic wave functions with density functional theory (DFT) methods. Their calculation was recently extended to a broader range of semiconductors in [22,24], using a similar procedure.
In this paper, we formulate a new approach to calculate the DM-electron scattering rate in the broad class of dielectric materials, by expressing the rate in terms of the dielectric response (ω, k). As the dielectric response is dominated by the electron response for energies (ω) above the band gap, this gives an alternate way to understand DM-electron scattering and leads to quantitatively different scattering rates, as screening effects are automatically included. Furthermore, the dielectric function is extensively studied in both condensed matter theory and experiment. Rewriting DM scattering in this way thus provides a more direct translation between quantities of interest for condensed matter and dark matter physicists. * simon.knapen@cern.ch † jkozaczuk@physics.ucsd.edu ‡ tongyan@physics.ucsd.edu The dielectric response of a material determines the energy loss function (ELF), which is defined as the imaginary part of the inverse dielectric function This quantity describes the rate to lose momentum k and energy ω for a charged particle passing through the material. 1 The ELF is closely related to the dynamic structure factor S(ω, k), which describes the rate to create density fluctuations in the medium, independently on the nature of the external probe. We can therefore directly relate the ELF to the dark matter scattering rate. By writing the scattering rate in terms of the ELF or S(ω, k), we are moreover accounting for in-medium effects such as screening, as well possible collective excitations such as plasmons [33,34]. The same ELF also plays an important role in the Migdal effect, or inelastic DM-nucleus scattering, which we studied in a companion paper [35]. (See also [36].) For any scalar or vector mediator coupling to electrons, we can treat DM as an external source which couples to electron number density n(r, t). In linear response theory, perturbations of the electron number density in the medium can be determined by the susceptibility given here in Fourier space, with V the volume and k and ω respectively the momentum and energy of the perturbation. The expectation value in (2) includes the thermal average. The fluctuation-dissipation theorem relates the susceptibility to the dynamical structure factor, which parametrizes the rate at which excitations are emitted or absorbed by the system: Here the dynamic structure factor is defined as with β = 1/k B T and Z the partition function of the system. Eq. (4) should remind the reader of Fermi's golden rule, and S(ω, k) is directly proportional to the differential DM-electron scattering rate. Using the relationship between the susceptibility and dielectric response we can write the structure factor as This relation is well known in the condensed matter literature, see e.g. [37].
In the remainder of this paper we explore the consequences of this relationship for dark matter electron scattering. The main difference with previous works in the literature is essentially that, writing the ELF as Im( (ω, k))/| (ω, k)| 2 , we see that a screening factor of 1/| (ω, k)| 2 is included inside the dynamic structure factor. Previous works studying DM scattering in semiconductors [3,22,24] primarily considered the approximation | (ω, k)| 2 ≈ 1. Since the DM scattering rate is dominated by k > ∼ keV, this assumption is not unreasonable, but with detailed calculations we find that screening can affect the rate by a factor of a few in Si and Ge. In addition, while the importance of accounting for screening has been well understood for vector mediators, screening for scalar-mediated scattering was only pointed out more recently in Ref. [38] (see also Ref. [39] for discussion of in-medium effects for scalars). In this work, we put scalar and vector mediated scattering on the same footing and show how they lead to identical response functions. We also show how scattering form factors discussed in the literature relate to the dielectric response, and perform detailed calculations of the screening effect in semiconductor targets relevant for current low-threshold experiments.
In the following section, we show how the DM-electron scattering rate relates to the dynamic structure factor or ELF. In section III, we discuss different ways to determine the dielectric function and thus the ELF, including the details of our DFT calculations for semiconductors. In section IV we present the implications for DM scattering in semiconductors and superconductors. We conclude in section V.

II. DM-ELECTRON SCATTERING AS DIELECTRIC RESPONSE
The most common models which predict dark matterelectron scattering involve a scalar or vector mediator which couple respectively to the electron number density and the electron current. In the nonrelativistic limit, the leading interactions of the mediator are the same for both cases: (7) since scattering via the 0 th component of the vector dominates. Here n χ and n are respectively DM and electron number densities. This makes is it manifest that in the non-relativistic limit the scalar and vector mediators ought to give identical rates, up to the rescaling of the coupling constants. Note that the vector here could represent a kinematically-mixed dark photon in the interaction basis, or another vector.
Given the similarity in these interactions, we can thus consider a general mediator with coupling to electrons g e and coupling to the DM g χ . We will write the mass of the mediator as m V , although it could also be a scalar. The coupling between the electron density perturbation n k and the external potential to the DM is then given by where the term in the parentheses represents the external and thus unscreened potential due to the DM (where x is DM position). In this basis, all in-medium corrections will be included in S(ω, k), as the propagator itself receives no corrections. In the particle physics literature the interaction term in (8) is often written in terms of the total potential felt by the electrons, especially so in the context of a kinetically mixed dark photon mediator. In this basis the propagator receives a multiplicative correction of the form 1/ (ω, k), and one defines a different structure factor, without the screening factor. The approaches are equivalent. However by working with the external rather than the total potential, the parallel between the scalar and the vector mediator in (7) is more manifest. Evaluating the Hamiltonian in (8) between initial and final DM states of momentum p i and p f , respectively, as well as initial and final electron fluid states |i , |f , we find the matrix element where in the continuum limit we can write the Kronecker delta function as a Dirac delta function, . We now use Fermi's Golden rule, and sum over initial states |i weighted by e −βEi /Z, as well as over final states. Inserting a factor of unity as dωδ(ω + E i − E f ), we obtain a DM scattering rate where ρ T is target density, µ χe is DM-electron reduced mass, and f χ (v) is the DM velocity distribution. Here we used the conventional definition of DM-electron scat-tering cross sectionσ e in terms of couplings [3]: and the DM-mediator form factor is defined as Plugging in (6), we arrive at our master formula for the scattering rate To compare this form of the rate with previous works in the literature, we use the Lindhard form for (ω, k). The Lindhard dielectric function, also known as the random phase approximation (RPA), is the leading-order polarization due to electron-hole excitations. It is given by [40,41] where we sum over states labeled by momentum p and band . There is also an implicit sum over spin states. The thermal occupation of the electron state with energy ω p, is f 0 (ω p, ) = 1/ [exp(β(ω p, − E F )) + 1], where E F is the Fermi energy. Using (14) in (6), we find Here we recognize a rate to create single-electron excitations, but with a screening factor of 1/| (ω, k)| 2 . Many works have considered this screening effect for vector mediators, by defining an effective coupling in the medium.
Here we show that it should apply to scalars too, and include it inside the dynamic structure factor.
Our formulation of DM-electron scattering in terms of a structure factor is then identical to that of Refs. [22,24] when the dielectric function is computed in RPA, with the only difference being the screening factor appearing in S(ω, k). Similarly, our results are equivalent to those of Ref. [3] when the RPA dielectric function is used and | (ω, k)| 2 → 1. More explicitly, we find that the crystal form factor of Ref. [3] is given by with V cell the volume of the unit cell. The relationships between the different conventions for DM-electron scattering are discussed in more detail in Appendix B.

III. DIELECTRIC RESPONSE
In this section, we discuss two approaches for determining (ω, k): a data-driven analytic approximation, and density functional theory (DFT) calculations. The various calculations of (ω, k) are compared with each other and with experimental data. Readers who are interested primarily in the DM scattering reach can proceed directly to Sec. IV.

A. Mermin oscillator model
A semi-analytic approximation to (ω, k) is valuable to quickly obtain results for many materials, as compared to numerically expensive DFT calculations. Direct measurements of (ω, k) at several ω, k points are moreover often available to anchor such a semi-analytic description. To fully make use of these measurements, a self-consistent interpolation is however needed which preserves the various sum rules and symmetries associated with (ω, k). For some materials the available data is also restricted to the optical (k = 0) limit, and a well-motivated extrapolation to finite k is therefore desirable.
One of the simplest, analytic models of (ω, k) is the Lindhard model for a homogeneous electron gas, for which the dielectric function can be characterized entirely by its Fermi velocity v F = k F /m e and plasma frequency ω p = 4πα em n e /m e , with n e the electron number density. The dielectric function is also isotropic in k in this case, and can be directly evaluated in (14) by inserting plane wave states. The result is [40] where we have explicitly separated the results into its real and imaginary parts. The main shortcoming of the Lindhard model is that the plasmon peak has zero width, which is certainly not the case in semiconductors such as Si and Ge. This problem is addressed in the Mermin model [42] Mer (k, ω; ω p , Γ) = 1 Lin(k,0)−1 (19) with Γ the width of the plasmon pole. By construction Mer (k, ω; ω p , 0) = Lin (k, ω; ω p ). Mer is moreover designed such that the various sum rules on the the dielectric function are explicitly satisfied.
Both the Mermin and Lindhard models however apply to a homogeneous electron gas, which is a far cry from a realistic material. This is often addressed in a phenomenological way by modeling the material as a superposition of many electron gas clouds with different densities. In other words, one describes ELF as a linear combination of Mermin dielectric functions. Here we follow the procedure outlined in [43,44] with where the ω p,i , Γ i , ω edge,i and A i (0) are fitted to experimental data. One can use as many Mermin oscillators as needed to describe the experimental data. The real part of 1/ (ω, k) can be obtained through a Kramers-Krönig transformation. This approach also makes it possible to include the semi-core electrons in a phenomenological manner [44], something which is computationally difficult to do in first principles DFT calculations. For our calculations we make use of the chapidif package 2 [45], with experimental inputs taken from [46,47], all obtained in the optical (k = 0) limit. In Fig. 1 we show ELF for Si and Ge, and compare with a DFT calculation with the GPAW code (see next section). For Si, we also compare with the finite-k data from Weissker et. al. [48], which is independent from the data used to fit to the Mermin model. We find good agreement between all three methods, except for the high-k, high-ω regime. Both the DFT and Mermin oscillator methods suffer from increased uncertainties in this regime: in the DFT calculation, higher values of ω require more bands to be included, which increases the computational complexity. With the Mermin oscillator method the uncertainties are expected to grow the further one deviates from the optical limit. This may be addressed by including finite-k data in the fit. For our numerical calculations of the DM scattering rate in semiconductors we will rely on the DFT method, and reserve a more detailed comparison to Appendix C.

B. DFT calculations
As an alternative to the phenomenological approach of the previous subsection, it is also possible to determine the dielectric response of a material from first principles. In contrast to the case of a homogeneous electron gas, in a crystal, response functions are only invariant under lattice periodicity, so that in momentum space the full dielectric function is written as a matrix GG (q) where q is restricted to the first Brillouin zone (1BZ) and G, G are reciprocal lattice vectors. The dielectric function is then treated as a matrix in reciprocal lattice vectors. Here, we provide an overview of microphysical calculations of this dielectric response in the framework of time-dependent density functional theory (TDDFT), while additional details are reviewed in Appendix A.
In the TDDFT approach, one maps the system of interacting electrons in the presence of an external (timedependent) potential to a system of non-interacting electrons in the presence of an effective potential. The latter is known as the Kohn-Sham (KS) system [50] and is much simpler to work with, since one has only to deal with effective single electron wavefunctions. Quantities When a measurement is available, it is overlaid as well. For k = 0 the Si and Ge data are taken from respectively [49] and [46]. At finite k for Si the measured ELF is taken from the Weissker et al. dataset [48]. The Mermin oscillators were fit to optical (k = 0) data [46,47]; the Weissker et al. data for finite k values is independent and not included in this fit. The discrepancy at high ω is due to the fact the GPAW calculation only includes the lowest 70 bands in computing the ELF, and hence does not capture the dielectric response for ω > ∼ 70 eV. Note that for DM with maximum speed of ∼ 750 km/s, only the phase space with k > ∼ 4 keV × ω/(10 eV) contributes to DM-electron scattering. The sharp plasmon resonance in the first two columns therefore does not contribute to the scattering rate. Those panels are meant only as a validation of our methods.
such as the susceptibility χ GG (q, ω) and the polarizability P GG (q, ω) can be related to their counterparts computed in the simpler Kohn-Sham system by requiring the change in charge density in response to a small change in the external potential (in the full system) and the effective potential (in the KS system) to be the same.
We are ultimately interested in the microscopic dielectric function, which is related to the polarizability by [51] GG (q, ω) = δ GG − 4πα em |q + G| |q + G | P GG (q, ω). (22) In the random phase approximation, the polarizability is approximated with the KS susceptibility P GG (q, ω) ≈ χ KS G,G (q, ω) (see Appendix A), and thus which, neglecting the off-diagonal pieces, is simply the Lindhard dielectric function of Eq. (14) computed with KS wavefunctions and extended to momenta k = q + G outside the 1BZ (see also Eq. (B1)). By solving for the susceptibility in the relatively simple KS system, one arrives at an approximation for the full microscopic dielectric function. There exist several DFT tools to compute the KS susceptibility, and hence the RPA dielectric response. We use the public code GPAW [52,53] for this purpose and focus on Si and Ge semiconductors. First, the KS wavefunctions are computed. This is done at zero temperature, using a plane-wave basis with a cutoff of E cut = 500 eV, corresponding to |k| < ∼ 22 keV. The Brillouin zone is sampled using a gamma-centered Monkhorst-Pack grid with 8 × 8 × 8 k points for Si, while for Ge we use a 12 × 12 × 12 grid. The finer grid for Ge was chosen to improve convergence of the results with respect to the grid spacing. Seventy bands are included for each spin. The KS wavefunctions are computed using the TB09 exchange-correlation functional [54], and a scissor correction is applied to match the experimentally measured Si and Ge bandgaps at T = 0. Note that the 3d electrons in Ge are treated as part of the frozen core, in contrast to e.g. [3].
Next, the longitudinal dielectric matrix is computed in the RPA using (23) for all q ∈ 1BZ sampled by the Monkhorst-Pack grid. We will work in an approximation where we neglect the directional dependence of the response, as well as the off-diagonal components of the dielectric matrix. To this end, we define an angularaveraged dielectric function where N (k) ≡ q,G δ k,|q+G| , the q sum runs over all 1BZ points sampled by the Monkhorst-Pack grid and and the G sum over all reciprocal lattice vectors up to the plane wave cutoff momentum. This quantity can then be used as an approximation to the full dielectric function in the ELF, Im(−1/ (ω, k)) Im(−1/ (ω, k)). This approach neglects so-called "local field effects" (LFEs), since the off-diagonal components of the dielectric matrix are dropped altogether.
In practice, some information about the off-diagonal components of the dielectric matrix can be included by replacing GG → 1/( −1 GG ). This is known as the "inclusion of LFEs" in the literature. Using this quantity in the ELF results in a better fit to experiment (see e.g. [48] and Fig. 1), since at low momentum transfer this procedure amounts to averaging out the effects of the off-diagonal components of the dielectric matrix [51]. Approximating the loss function with (ω, k) with or without LFEs does not make a substantial difference in the experimental sensitivity to DM-electron scattering presented in the next section. We include local field effects except where stated otherwise, so that the loss function predicted by GPAW more closely matches experimental results.
The results computed by GPAW for the ELF in Si and Ge are illustrated in Fig. 1 for various values of k. We see that generally the DFT results agree well with both experimental results (where available) and the Mermin approach described in the previous subsection. The discrepancies at large ω are due to the fact that GPAW only includes the lowest 70 bands in computing the loss function, so does not yield reliable results above ∼ 70 eV for Si and Ge.

IV. IMPLICATIONS FOR DM-ELECTRON SCATTERING
To show the impact of screening, we now evaluate the scattering rate in example dielectric materials. We will consider the 'massless mediator' limit where m V αm e with F DM (k) = (αm e ) 2 /k 2 and the 'massive mediator' limit where m V αm e with F DM (k) = 1. As discussed before, the results here apply for both vector and scalar mediators.
Our main results focus on Si and Ge semiconductors, which are used in a number of direct detection experiments. We use (ω, k) computed in the DFT framework as described in the previous section, taking as our default the RPA dielectric function including local field effects. Again, there is only a small difference in rate whether local field effects are included or not, and we show an explicit comparison in Fig. 5 in Appendix B. The Mermin oscillator determination of (ω, k) also gives comparable results as long as we do not consider ω too close to the band gap, which is reasonable for backgroundlimited experiments. The results with the Mermin oscillator method are given in Appendix C. For the DM velocity distribution, we assume the Standard Halo Model with v esc = 500 km/s, velocity dispersion v 0 = 220 km/s, and Earth velocity v e = 240 km/s. Fig. 2 shows the impact of screening on the differential rate spectrum, for an example DM mass of 10 MeV. Here the unscreened rate (dashed lines) is obtained by writing the ELF as Im( (ω, k))/| (ω, k)| 2 and taking | (ω, k)| 2 → 1. The screening effects are most noticeable for lower energy deposition ω, since in that case there is a larger contribution from lower momentum transfers where the screening is largest. Scattering at large ω is dominated by large k, with negligible screening. Similarly, we see that the effect of screening is larger for the massless mediator   4. Ratio of the screened rate to the unscreened rate, for different thresholds corresponding to 1, 2 and 3 electrons. We use Q = 1 + (ω − Eg)/ε where for Si Eg = 1.11 eV, ε = 3.6 eV and for Ge Eg = 0.67 eV, ε = 2.9 eV, following Ref. [3].

Massive mediator
case, since the DM form factor F DM (k) enhances the rate from lower k values.
We show the corresponding effect on the DM mass and cross section reach in Fig. 3. The solid lines show the reach for scalar and vector mediators, accounting for screening effects. We assume kg-year exposure, zero background, and 95% CL projected reach to match with the convention in the literature. The threshold is set by the electron band gap. For m χ > ∼ 10 MeV, there is roughly a factor of (1.4) 2.5 suppression in the total rate for (massive) massless mediators. The ratio becomes larger near threshold in m χ , since for those points the rate is restricted to ω near the band gap, where screening is more important. The screening effect is therefore reduced somewhat with higher thresholds in ω, as shown in Fig. 4. For instance, the threshold to detect 2 electron-hole pairs is roughly 4.7 eV (3.6 eV) in Si (Ge). Setting this as the threshold, we find a screening suppression instead of 2-2.1 for massless mediators and m χ > ∼ 10 MeV. For massive mediators the dependence on the energy threshold is smaller.
The O(1) screening effects we find for Si and Ge align with our expectations for semiconductors with eV-scale electron band gaps, and it is therefore interesting to compare with a lower gap material where the screening is much stronger. We also show in Fig. 3 the reach in a metal, taking Al as an example. Such targets have been proposed to be used in their superconducting phase as low-threshold dark matter detectors [16,17,26,55]. We thus consider sensitivity to electron recoils in the energy range 10 meV -1 eV, such that the material can still be approximated with the dielectric response of a metal.
Here we use the Mermin oscillator method with the Al data from [47]. For ω > 10 meV, we find that the rates are in good agreement with those obtained with the Lindhard dielectric function for a free electron gas in (17), taking ω p = 15 eV. For ω < 10 meV the agreement between the methods is not as good, for reasons to be understood further. Out of an abundance of caution we therefore impose a ω > 10 meV threshold in Fig. 3. For massive mediators, the screening strongly limits the sensitivity to sub-MeV dark matter despite the lower thresholds. For massless mediators, in the absence of screening there is enhanced scattering with low k and lower thresholds, and the unscreened reach is many orders of magnitude below what is shown on the plot. Accounting for screening, we find that there is still substantial reach to sub-MeV dark matter scattering via a massless mediator. Thus, even with the large screening, such a low gap target could be sensitive to cosmologically interesting sub-MeV dark matter models such as that of freeze-in through a kinetically-mixed dark photon [3,56,57].

V. CONCLUSIONS
By considering the linear response of a dielectric material, we have shown that the differential DM-electron scattering rate in a dielectric is proportional to the energy loss function Im[−1/ (ω, k)] (see (13)), which contains all relevant many-body effects associated with the target material. The ELF is moreover very well studied theoretically and experimentally in the materials science literature, and thus provides a convenient way of mapping the detailed properties of the target material onto sensitivity estimates or limits for DM direct detection experiments. In particular, we find that screening effects need to be accounted for, both for scalar and vector mediators, which reduces the reach of any direct detection experiment with a dielectric target. We computed the ELF for Si and Ge using a first principles DFT calculation and using a datadriven, phenomenological model. Both methods broadly agree within their regime of validity. Using these results, we can quantify the importance of the screening effect in Si and Ge (see Fig. 4). There are a number of possible future directions to pursue, such as accounting for angular dependence in the dielectric response for semiconductors, and applying our methodology to a broader range of materials, including others already proposed for the direct detection of electron recoils.

ACKNOWLEDGMENTS
We thank Diego Redigolo for collaboration in early stages of this work and for useful discussions. We also thank Yonit Hochberg, Yoni Kahn, and Noah Kurinsky for useful discussions. We thank Maarten Vos for providing us with a β-version of his chapidif package and for his assistance with its usage and the interpretation of the results. TL is supported by the Department of Energy under grant DE-SC0019195 and a UC Hellman fellowship. JK is supported by the Department of Energy under grants DE-SC0019195 and DE-SC0009919.

Appendix A: Dielectric response in a crystal
In our analysis we considered the scalar longitudinal dielectric function. The more general quantity in a crystal is a dielectric tensor that is a matrix both in spatial indices and in reciprocal lattice vectors. For completeness, here we introduce the dielectric tensor and detail the approximations made in the main text.
The dielectric tensor describes the electrical response of a system to an external electric field, E ext . In terms of microscopic quantities, the relationship between the external and total electric fields is which serves as the definition of −1 ij (ω, r, r ). The Latin subscripts correspond to spatial indices. The fields can be written in terms of their Fourier components and similarly for E ext . In the above expression, q lies in the first Brillouin Zone (1BZ) and G, G are reciprocal lattice vectors. The dielectric tensor can be Fourier transformed as where we have used the fact the microscopic electronic response of a crystal is invariant under translations by a lattice vector R, so ij (ω, r, r ) = ij (ω, r + R, r + R).
In Fourier space (A1) is then where again q ∈ 1BZ. In considering dark matter-electron scattering we are primarily interested in the longitudinal response. The longitudinal field is defined as E L (ω, k) ≡ k·E(ω, k)/ |k|, where k is a general momentum vector and not necessarily confined to the 1BZ. The transverse field is E T ≡ E − (k · E(ω, k)) k/ |k| 2 . The external field and dielectric tensor can be similarly decomposed. The scalar dielectric function is obtained by projecting both the total and external electric fields onto their longitudinal components: In (A5) we have defined the (symmetrized 3 ) longitudinal dielectric function −1 which describes the longitudinal response to a longitudinal external field. One can also define the matrices −1 LT, T L, T T to describe the other components of the response, but for nearly isotropic crystals and the energies of interest for dark matter scattering, the purely longitudinal contribution dominates both E ext and the response. For compactness, we denote −1 The microscopic dielectric function can be computed from the density response function χ GG (q) in density functional theory, as described in Sec. III B. One can show (see e.g. [51]) that the susceptibility in the full interacting system is related to the KS susceptibility χ KS via a Dyson equation where f xc is a so-called "exchange correlation kernel" which is defined such that the charge density of the KS system exactly matches that of the full system. Exact knowledge of f xc would require solving for the wavefunctions of the full interacting system, however in TDDFT calculations one typically approximates this term using simple physically-motivated models such as the "adiabatic local density approximation" (ALDA), or dropping it altogether. Setting f xc → 0 in the expression above corresponds to the "random phase approximation" (RPA) which we use throughout this study. Furthermore, the susceptibility and polarizability are 3 Note that there exists an alternative definition of the longitudinal dielectric function in the literature. If a factor of |q + G| / |q + G | is included on the right hand side of (A5), one arrives at the "unsymmetrized longitudinal dielectric function". The G = G elements of the symmetrized and unsymmetrized quantities are the same, but the off-diagonal components are not, so one should be careful to use a consistent definition in considering local field effects.
related by a separate Dyson equation Comparing (A8) to (A7) with f xc → 0, we see that in the RPA the polarizability of the full system is governed by the same Dyson equation as the KS susceptibility, motivating the approximation P G,G (q, ω) ≈ χ KS G,G (q, ω). This gives the RPA dielectric function of Eq. (B1) (see below), and reduces to that of Eq. (23) if the off-diagonal components are neglected and one restricts the momentum transfer to lie within the 1BZ.

Appendix B: Comparison with previous works
The dynamic structure factor can be directly related to the DM scattering form factors appearing elsewhere in the literature when the Lindhard dielectric function (or random phase approximation) is used. In this section, we provide some additional formulae to help translate the presentation here in terms of (ω, k) to the results appearing in several previous studies. We also provide some plots comparing our results with those in Essig et. al. [3] and Griffin et. al. [22].
In the main text, the Lindhard dielectric function given in (14) was only valid for k within the 1BZ for a crystal.
Here we generalize to account for reciprocal lattice vectors and split k = q + G where q lies in the 1BZ. The Lindhard dielectric function can be written as [41] RPA where the matrix elements above are defined by with p within the 1BZ and G, G reciprocal lattice vectors. Here we have assumed that the real-space Bloch wavefunction for an electron in band can be written as where u p, (r) is periodic under under r → r+R, with R a lattice vector. In the second equality above, we have written the wavefunction in terms of the momentum-space coefficients u (p + G). The matrix element above can equivalently be written in momentum space as  The solid and dashed blue lines correspond to screened and unscreened rate obtained with our calculations of the dielectric function using GPAW, including local field effects. The dotted blue line shows the result if we compute the unscreened rate with the RPA dielectric function from GPAW without local field effects. This case is the one that corresponds closely to previous calculations, and numerically we find very good agreement with Griffin et al. [22] and with QEDark [3].
where in the last line we make contact with the notation of Refs. [3,22].
To compute the ELF, the dielectric matrix must be treated as a matrix in reciprocal lattice vectors and inverted to obtain the inverse dielectric function. In order to compare with results in the literature, we work in the approximation that the off-diagonal elements can be neglected, and restrict to G = G in (B1). Then taking the imaginary part of the dielectric function above gives We use (6), take the continuum limit, and now explicitly include a factor of 2 for the spin sum. (This was implicit in equations in the main text.) We find that the structure factor can be written as where we introduce the sum over G to select out the piece of the incident DM momentum k that brings it to the first BZ. As noted before, this agrees with the definition of the structure factor in Ref. [22] except for the 1/| (ω, k)| 2 screening factor appearing here.
To connect with the definitions in Ref. [3], which also averages over all directions in calculating the rate, we replace the 3-dimensional momentum delta function with a delta function averaged over the sphere, δ 3 (p − p + k − G) → δ(k − |p − p + G|)/(4πk 2 ). From this, we can immediately compare with the definition of the isotropic crystal form factor appearing there, and obtain (16) by neglecting the factor of (1−e −βω ) in the low temperature limit.
In Fig. 5 we show a comparison of various calculations of the cross section reach, taking here v esc = 600 km/s, v 0 = 230 km/s, and v e = 240 km/s for direct comparison. We show the screened and unscreened reach using our default calculations of the dielectric function, which account for local field effects, as discussed in Sec. III B.
There is a small difference in the rate if we use the RPA dielectric function without local field effects (dotted light blue). This unscreened rate corresponds to the calculation of Refs. [3,22], and with which our results agree very well. For Ge and scattering via massless mediators, there are somewhat larger differences in the unscreened rate across different calculations, which may be due to differences in the various DFT calculations (choice of exchange-correlation functionals, lattice constants, etc).

Appendix C: Mermin oscillator results
In this appendix we briefly present some results obtained with the Mermin oscillator method (see Sec. III A), and study how they compare with those obtained with the DFT calculation (see Sec. III B). With both methodologies, accessing the high-k regime is challenging, for different reasons. In the DFT calculation, an increasingly large basis set of wave functions is needed, which increases the computational cost of the calculation. In the Mermin oscillator method, the high k regime corresponds to a substantial extrapolation from the experi-mental data, which was taken in the optical limit (k = 0). For k > ∼ 12 keV our numerical results with the Mermin oscillator method in particular cease to be stable and we therefore impose a cut of k < 12 keV on the phase space in both calculations. We verified that the contribution of the omitted part of the phase space is negligible in the integrated rate, but it slightly affects the shape of dR/dω for ω > ∼ 15 eV..
With this assumption, Fig. 6 shows the differential scattering rate obtained with both methods. We find overall good agreement, except for low and high ω. Poor agreement at low ω is anticipated, since the Mermin oscillator method models the semiconductor as a linear combination of free electron gas systems, and is therefore expected to be less reliable for ω near the band gap of the material. Once we impose the 2e − threshold (dashed line), the agreement between both methods is largely satisfactory. The substantial deviations in the high ω regime are also straightforward to understand. For kinematical reasons, this regime corresponds to the higher k part of the phase space, which is challenging for both methods as discussed above. Further studies are needed to bring down the uncertainty in this region. On the other hand, for experiments with a 2e − threshold, this region provides a subdominant contribution to the rate, and is likely only relevant in the event of a discovery.
The integrated rate above the 2e − threshold is shown in Fig. 7. For Si, both methods agree to within 10%, and for Ge the agreement is within roughly 30% in most of the mass range. The uncertainties increase for low masses, due to the challenge of modeling the ELF accurately for ω close to the band gap.