Extended Calculation of Dark Matter-Electron Scattering in Crystal Targets

We extend the calculation of dark matter direct detection rates via electronic transitions in general dielectric crystal targets, combining state-of-the-art density functional theory calculations of electronic band structures and wave functions near the band gap, with semi-analytic approximations to include additional states farther away from the band gap. We show, in particular, the importance of all-electron reconstruction for recovering large momentum components of electronic wave functions, which, together with the inclusion of additional states, has a significant impact on direct detection rates, especially for heavy mediator models and at $\mathcal{O}(10\,\text{eV})$ and higher energy depositions. Applying our framework to silicon and germanium (that have been established already as sensitive dark matter detectors), we find that our extended calculations can appreciably change the detection prospects. Our calculational framework is implemented in an open-source program $\texttt{EXCEED-DM}$ (EXtended Calculation of Electronic Excitations for Direct detection of Dark Matter), to be released in an upcoming publication.

Reliable theoretical predictions of target-specific transition rates are important not only for current experiments, but also for planning the next generation of detectors. Compared to the DM-induced electron ionization rate in noble gases like xenon [1][2][3][4][5][6][7][8], calculations for the DM-electron scattering rate in a crystal are more complicated. Ionization rates for noble gases can be calculated by considering each noble gas atom as an individual target, where the calculation simplifies to finding the ionization rate from an isolated atom, for which the wave functions and energy levels are well tabulated [43]. However, for crystal targets the atoms are not isolated and more involved techniques are required to obtain an accurate characterization of DM-electron interactions in a many-body system.
There have been a variety of approaches taken to compute the DM-electron scattering rate in crystals. One of the first attempts, Ref. [2], computed the rate with semi-analytic approximations for the initial and final state wave functions, and used the density of states to incorporate the electronic band structure. Later, Ref. [3] continued in this direction and used improved semi-analytic approximations for the initial state wave functions. Meanwhile, a fully numerical approach was advanced in Refs. [1,10,11] where density functional theory (DFT) was employed to calculate the valence and conduction electronic band structures and wave functions. The latter approach, as implemented in the QEdark program and embedded in the Quantum ESPRESSO package [44][45][46], has become the standard for first-principles calculations of DM detection rates. Recently, in Refs. [15,16] we used a similar DFT approach as implemented in our own program for a study of DM-electron scattering in a variety of target materials. More recently there has been work utilizing the relation between the dielectric function and the spin-independent scattering rate [47][48][49], which also properly incorporates screening effects.
The goal of this work is to further extend the DM-electron scattering calculation in several key aspects, and present state-of-the-art predictions for Si and Ge detectors using a combination of DFT and semi-analytic methods. As we will elaborate on shortly, the time-and resource-consuming nature of DFT calculations presents an intrinsic difficulty that has limited the scope of previous work in this direction to a restricted region of phase space; typically only bands within a few tens of eV above and below the band gap were included and electronic wave functions were cut off at a finite momentum. We overcome this difficulty by implementing all-electron (AE) reconstruction (whose importance was previously emphasized in Ref. [50]) to recover higher momentum components of DFTcomputed wave functions, and by extending the calculation to bands farther away from the band gap using semi-analytic approximations along the lines of Refs. [2,3]. As we will see, the new contributions computed here have a significant impact on detection prospects in cases where higher energy and/or momentum regions of phase space dominate the rate, including scattering via a heavy mediator, and experiments with O(10 eV) or higher energy thresholds. We also stress that in contrast to the recent work emphasizing the relation between spin-independent DM-electron scattering rates and the dielectric function [47][48][49], our calculation can be straightforwardly extended to DM models beyond the standard spinindependent coupling. Furthermore, we do not make assumptions about isotropy for the majority of our calculation, and our framework is capable of treating anisotropic targets which exhibit smoking-gun daily modulation signatures [15,23,24] (see also Refs. [51,52] for discussions of daily modulation for phonon excitations).
Our calculation is implemented in an open-source program EXCEED-DM (EXtended Calculation of Electronic Excitations for Direct detection of Dark Matter), to be released in  Si (left) and Ge (right), divided into core, valence ("val"), conduction ("cond") and free. Shaded regions indicate the range of energies for each type of electronic states. In a scattering process, electrons transition from either core or valence states, below the Fermi surface at E = 0, to conduction or free states above the band gap E g . As outlined in Sec. 1.1 and explained in detail in Sec. 2, we compute the valence and conduction states numerically using DFT (including all-electron reconstruction), model the core states semi-analytically with RHF wave functions, and treat the free states as plane waves. an upcoming publication. Currently a beta version of the program is available here [53]. We also make available our DFT-computed wave functions [54] and the output of EXCEED-DM [55] for Si and Ge.

Overview of the Calculation and Key Results
Before delving into the technical details, let us give a brief overview of the calculation and highlight some key results. We divide the electronic states in a (pure) crystal into four categories: core, valence, conduction and free, as illustrated in Fig. 1 for Si and Ge and discussed in more detail in Sec. 2. At zero temperature, electrons occupy states up to the Fermi energy, defined as the maximum of the valence bands and denoted by E = 0. The band gap E g , i.e. the energy gap between the occupied valence bands and unoccupied conduction bands, is typically O(eV) for semiconductors, e.g. 1.11 eV for Si and 0.67 eV for Ge; this sets a lower limit on the energy deposition needed for an electron transition to happen.
The electronic states near the band gap deviate significantly from atomic orbitals and need to be computed numerically. We apply DFT methods (including AE reconstruction) for this calculation, and refer to the DFT-computed states as valence and conduction. Specifically, for both Si and Ge, we take the first four bands below the gap to be valence, which span an energy range of −12 eV to 0 and −14 eV to 0, respectively, and take bands above the gap up to E dft = 60 eV to be conduction.
QEdark [11] c → c  Figure 2. Selection of results from Sec. 4, for DM-electron scattering via a heavy mediator in a Ge target. Left: Contribution from each of the four transition types, valence to conduction (v→c), valence to free (v→f), core to conduction (c→c), and core to free (c→f) to the scattering rate binned in energy deposition (with ∆ω = 1 eV) for a 1 GeV DM at a given reference cross section σ e = 10 −40 cm 2 . Right: 95% C.L. projected limit (3 events) on σ e assuming 1 kg-year exposure, for energy thresholds corresponding to 1 and 5 electron-hole pairs. We compare our results with QEdark calculations in Refs. [10,11] and the semi-analytic model of Lee et al [3]; see text for details.
With more computing power we can in principle include more states, both below and above the band gap, in the DFT calculation. However, since the states far from the band gap can be modeled semi-analytically with reasonable accuracy, computing them with DFT is inefficient. Below the valence bands, electrons are tightly bound to the atomic nuclei. We model them using semi-analytic atomic wave functions and refer to them as core states. These include the 1s, 2s, 2p states in Si and 1s, 2s, 2p, 3s, 3p, 3d states in Ge (the 3d states in Ge are sometimes referred to as semi-core, and we will compare the DFT and semi-analytic treatment of them in Secs. 2.2 and 3.3). Finally, above E dft = 60 eV, we model the states as free electrons as they are less perturbed by the crystal environment.
With the electronic states modeled this way, we compute the rate for valence to conduction (v→c), valence to free (v→f), core to conduction (c→c) and core to free (c→f) transitions induced by DM scattering, as discussed in detail in Sec. 3. The total rate is the sum of all four contributions. We then use our calculation to update the projected reach of direct detection experiments in Sec. 4, and compare our results with previous literature. Figure 2 gives a glimpse of some of our key results. Here we consider the case of DM scattering via a heavy mediator in a Ge target. The impact of core (3d) to conduction contributions is clearly visible from both the differential rate (left panel, for m χ = 1 GeV) and the projected reach (right panel). They dominate the total rate for m χ 10 MeV, and, as we can see from the right panel of Fig. 2, lead to significantly more optimistic reach compared to previous DFT calculations implemented in QEdark [10,11]; this is especially true for higher detector thresholds (corresponding to higher Q values). Note that while Refs. [10,11] included the 3d states in their DFT calculation, their contributions were significantly underestimated due to the absence of AE reconstruction. The importance of AE reconstruction is also seen from the valence to conduction differential rate in the left panel of Fig. 2, where our calculation predicts a much higher rate at ω 15 eV compared to the QEdark calculation in Ref. [11]. Meanwhile, accounting for in-medium screening (see Sec. 3.5) we find, consistent with Ref. [48], a lower rate at energy depositions just above the band gap, and as a result, weaker reach at low m χ , compared to Refs. [10,11]. On the other hand, our modeling of the core (3d) states is similar to the semi-analytic approach of Ref. [3], and indeed we find very similar reach at large m χ ; however, the approach of Ref. [3] overestimates the rate at smaller m χ due to reduced accuracy in the modeling of the valence and conduction states. We reserve a more detailed comparison with the literature for Sec. 4.1.

Electronic States
To compute the DM-electron scattering rate one must understand the electronic states of the target. In targets with a periodic potential, Bloch's theorem states that the energy eigenstates can be indexed by a momentum, k, which lies within the first Brillouin zone (1BZ). These Bloch states, ψ i,k , where i represents additional quantum numbers, are eigenstates of the discrete translation operator such that ψ i,k (x + r) = e ik·r ψ i,k (x), which means the electronic wave functions can be written as where u i,k (x + r) = u i,k (x) and V is the target volume. For every k there exists a tower of eigenstates (labeled by i) of the target Hamiltonian which constitutes the complete set of states in the target. Unfortunately this complete set is not known for a general material and therefore a combination of approximations must be used to calculate them. As discussed in Sec. 1.1 and illustrated in Fig. 1, we divide the states into core, valence, conduction and free, and use a combination of numerical calculations and semi-analytic modeling. In this section, we expand on the treatment of each type of electronic states. We first discuss the DFT calculation for valence and conduction states in Sec. 2.1, and then move on to explain the semi-analytic treatment of core states in Sec. 2.2. Our main results are contained in Fig. 3 where we compare the average magnitude of electronic wave functions, binned in momentum (see Eq. (2.7)), computed with and without AE reconstruction, discussed further in Sec. 2.1.1, and, for the highest energy core states (2p in Si and 3d in Ge), those computed using the core approximation discussed in Sec. 2.2. We find that the AE reconstruction includes a significant contribution from wave functions at large momentum as expected, and that for the core states, the semi-analytic approach reproduces the large momentum components of these AE reconstructed DFT wave functions. Lastly we will discuss the analytic treatment of the free states in Sec. 2.3.  , computed with DFT with (red, "AE") and without (blue, "no AE") AE reconstruction, and the semi-analytic core approximation of Eq. (2.9) (green, "core"). Shaded bands indicate the maximum and minimum values across all the bands belonging to the state type indicated in the upper right corner of each panel. AE reconstruction, discussed in Sec. 2.1.1, recovers the large momentum behavior of the electronic wave functions. Core electronic states, such as those shown in the right panels and discussed in Sec. 2.2, can be well modeled semi-analytically with atomic wave functions, as seen by the good agreement between the "core" and "AE" curves. When applicable, the semi-analytic parameterization is advantageous since the electronic wave functions are then known to arbitrarily large momentum.

DFT Wave Functions and Band Structures
In principle, DFT provides an exact solution to the many-electron Schrödinger equation by the Hohenberg-Kohn theorems that treat all properties of a quantum many-body system as unique functionals of the ground state density. They further show that the exact ground state density and energy can be found by minimizing the total energy of the system [56,57]. This becomes tractable by the Kohn-Sham (KS) equations that reduce the many-body problem to non-interacting electrons moving in an effective potential, V eff , where i is the orbital energy of the KS orbital ψ i [58]. The external potential V ext and Hartree potential V H , are known, while the exchange-correlation (xc) potential V xc , which contains the many-body interactions, must be approximated. Herein lies the deviation from the exact solution, and although various formulations of xc-energy functionals have been successful, the choice of xc-functional will affect the predicted electronic states and hence calculated transition rates. For Si, we use PBE [59], a type of generalized gradient approximation (GGA) xc-functional which is one of the most popular and low-cost choices. Local and semi-local based xc-functionals, such as PBE, suffer from a self-interaction error and band gap underestimation, which we modify with a "scissor correction" where the bands are shifted to match the experimentally determined values of band gap. For Ge, this underestimation results in zero band gap with PBE, therefore we instead use a hybrid functional, which mixes a parameterized amount of exact exchange into the xc-functional, correcting band gaps and band widths by error cancellation at the cost of increased computation time. We use the range-separated hybrid functional HSE06 [60,61], which applies a screened Coulomb potential to correct the long-range behavior of the xc-potential, giving high accuracy at a mid-level computational cost. Our computed band structures for Si and Ge are shown in Fig. 4.
The periodic Bloch wave functions, u i,k (x), Eq. (2.1), for band i and Bloch wave vector k are computed by finding the Fourier coefficients, u i,k,G (which satisfy the normalization condition, G | u i,k,G | 2 = 1): (2. 3) The number of reciprocal lattice vectors G kept in the sum is conventionally set by an energy cutoff, E cut , such that |k+G| 2 < 2m e E cut . These Bloch wave function coefficients u i,k,G for both Si and Ge are computed with the projector augmented wave (PAW) method [62,63] within vasp [64][65][66][67] up to E cut = 1 keV on a 10 × 10 × 10 uniform k mesh over the 1BZ. We then include AE reconstruction effects up to a higher energy cutoff, E AE = 2 keV, which recovers higher momentum components of the wave functions up to |k + G| 2 < 2m e E AE , as discussed in more detail in Sec. 2.1.1. The Bloch wave function coefficients, u i,k,G , for Si and Ge used for this work can be found here [54]. A final consideration of using DFT wave functions is that DFT is fundamentally a ground state method, and the KS conduction band states are only approximations to excited states. Excited state methodologies are much more computationally expensive than ground state KS-DFT. Furthermore, since excited state quasiparticles, such as excitons, have been argued to have a negligible effect on the calculation of DM scattering rates [11], they are neglected in our calculations.

All-electron Reconstruction
There are many different approaches to find the eigenstates of Eq. (2.2). The PAW method [62] is one such standard approach. The main idea of the PAW method is to split up the calculation of the eigenstates: near the ionic centers the wave functions resemble the eigenstates of an isolated atom, while further away they can be computed numerically with a pseudopotential. This greatly simplifies the numeric calculation since the focus is then on large distance (small momentum), and the small distance (high momentum) pieces can be self-consistently reintroduced after the main part of the DFT calculation. The large distance components of the wave function are known as "pseudo wave functions" (PS wave functions) and the total wave functions are known as the "all-electron wave functions" (AE wave functions), indicating that all of the wave function components are included. We will now give a brief overview of how the AE wave functions can be reconstructed from the PS wave functions, computed with PAW-based DFT codes, and refer the reader to Refs. [62,63,68,69] for more detailed information. 1 The AE wave functions, |Ψ AE are built from two components. Near the ionic core, or inside an "augmentation sphere," |Ψ AE is expanded in a set of basis functions, |φ AE j , which are simply taken to be the wave functions of an isolated atom, (2.4) Outside of the augmentation sphere, |Ψ AE = |Ψ PS . Near the ionic core the PS wave functions |Ψ PS are expanded in a set of basis functions |φ PS j that are computationally more convenient than the |φ AE j . Therefore, which simply replaces the components in |Ψ PS within the augmentation sphere with the AE wave function. To find the c coefficients we insert an identity, are projector functions, defined to satisfy this identity within the augmentation sphere. Therefore, This implies that all the PS states are related to AE states by this transformation T , such that c j = c j and the AE reconstruction can be written as In practice, we implement the AE reconstruction with pawpyseed [69], and the plane wave expansion cutoff of |Ψ AE , E AE , can be increased from the initial E cut . We use E AE = 2 keV.
To visualize the effect of AE reconstruction, we plot in Fig. 3 the average magnitude of the Bloch wave functions, binned in q, where u i,k,G are the Fourier components of the Bloch wave functions, defined in Eq. (2.3). Each bin in momentum space extends from q to q + ∆q with ∆q = 1 keV, and N q is a normalization factor equal to the number of points in a bin, N q = k G θ(q + ∆q − |k + G|) θ(|k + G| − q). We see that AE reconstruction recovers the high momentum components, which as we will see can significantly affect the DM-induced transition rate for processes which favor large momentum transfers (such as processes mediated by heavy particles), or processes limited to larger ω (e.g. higher experimental thresholds where large q processes are the only kinematically allowed transitions). Previous DFT calculations of DM-induced electron transition rates, with the exceptions of Ref. [15,16,50], used only the pseudo wave functions, |Ψ PS as opposed to the AE wave functions, |Ψ AE , and have therefore underestimated detection rates in several cases.

Atomic Wave Functions
If one could reconstruct the AE wave functions arbitrarily deep into the band structure, and to arbitrarily high momentum, one could calculate an accurate representation of the complete set of electronic states with a DFT calculation. In practice, however, this is neither feasible nor necessary. States deep in the band structure are more isolated from the influence of the crystal environment, and so an isolated atomic approximation becomes valid. We refer to these inner, tightly bound electrons as core electrons. In Si, we will show that the 2p states and below can be treated as core, while in Ge, the 3d states and below can, as alluded to in Fig. 1. The purpose of this subsection is to expand on the atomic approximation for core electrons and discuss its accuracy.
More precisely, the initial states of a transition should be taken as a linear combination of isolated atomic wave functions that is in Bloch form (known as Wannier states): where κ labels the atom in the primitive cell, n, l, m are the standard atomic quantum numbers, x κ is the equilibrium position of the κ th atom, r sums over all primitive cells in the lattice, and N is the total number of cells. In contrast to the valence and conduction states discussed in the previous subsection, the core states are labeled by (κnlm) rather than band index i. The corresponding periodic (dimensionless) u functions can be easily obtained via Eq. (2.1): where Ω = V /N is the primitive cell volume.
In general, the atomic wave functions ψ atom κnlm are not known analytically, but are expanded in a basis of well-motivated analytic functions. The basis coefficients are then fit by solving the isolated atomic Hamiltonian, giving a semi-analytic expression for ψ atom κnlm . We use a basis of Slater type orbital (STO) wave functions whose radial component is where a 0 = 0.53 Å = (3.7 keV) −1 is the Bohr radius, and Z is an effective charge of the ionic potential. Including the angular part, the atomic wave functions are then where C jln,κ , Z jl,κ , n jl,κ are tabulated in Ref. [43], and Y m l (x) are the spherical harmonics with the Condon-Shortley phase convention [70].
To assess the accuracy of the atomic wave function approximation, we temporarily push the DFT calculation beyond its default regime (valence and conduction), to the highest core states -2p states in Si and 3d states in Ge, where it is still computationally feasible -and compare the numerical wave functions to the semi-analytic ones discussed above. The results, in terms of the average magnitude of Bloch wave functions defined in Eq. (2.7), are shown in the right panels of Fig. 3. 2 We see that the atomic approximation accurately reproduces the numerical wave functions up to the momentum cutoff √ 2m e E AE 50 keV for E AE = 2 keV. These plots also show the limitation of DFT calculations. While AE reconstruction recovers higher-momentum components of electronic wave functions, it is not feasible to expand the plane wave basis set to arbitrarily high cutoff. However, having verified the atomic approximation for the highest core states, we can use it for all core states with confidence, allowing us to more easily include the high momentum components beyond the DFT cutoff.

Plane Wave Approximation
With the inclusion of the semi-analytic core states, all of the states below the band gap have been modeled. States above the band gap can also be computed with DFT methods, as described in Sec. 2.1. Similar to valence bands, there are practical limitations to how many conduction bands can be included. To remedy this in the simplest way possible, we model states far above the band gap as plane waves, where G is a reciprocal lattice vector, and plays the role of a band index. (To understand this, simply note that every momentum can be decomposed into a k vector inside the 1BZ and a reciprocal lattice vector. Integrating over the momentum of plane wave states amounts to a k integral within the 1BZ and a G sum.) From Eq. (2.1) we see that the corresponding periodic u functions are simply The plane wave approximation is often used in atomic ionization rate calculations, with the inclusion of a Fermi factor, F (ν), where E is the final state electron energy, and Z eff is an effective charge parameter, which enhances the transition rate at low E to account for the long range behavior of the Coulomb potential. See Refs. [1][2][3]6] for more details. In atomic ionization calculations one usually takes Z eff to be related to the binding energy of the initial state, E B , where n is the principal quantum number. Since the rate is proportional to the Fermi factor, Z eff = 1 is seen as the conservative choice. Later in Secs. 3.2 and 3.4 we quantify how much of an effect this has on the transition rate. This uncertainty is only important for very high experimental thresholds, and generally we find that Z eff = 1 leads to a smoother match (within an O(1) factor) to conduction band contributions from DFT calculations.

Electronic Transition Rates
We now present the DM-induced electron transition rate calculation. We begin with a general discussion and then in Secs. 3.1-3.4 consider the four different transition types in turn: valence to conduction (v → c), valence to free (v → f), core to conduction (c → c) and core to free (c → f). Finally, in Sec. 3.5 we discuss the treatment of in-medium screening.
The general derivation has been discussed previously (see e.g. Refs. [2,5,10,15,50]), and we repeat it here for completeness and clarity, as a variety of conventions have been used. Beginning with Fermi's Golden Rule, the transition rate between electronic states |i, s and |f, s due to scattering with an incoming non-relativistic DM particle, χ, with mass m χ , velocity v, and spin σ is given by where |p, σ; i, s = |p, σ ⊗ |i, s , q is the momentum deposited onto the target, p = m χ v, p = p − q, δĤ is the interaction Hamiltonian, V is total volume of the target, and ω q is the energy deposition: We assume that all quantum states are unit normalized. Modulo in-medium screening effects, discussed below in Sec. 3.5, we can write Eq.
We find where ψ i (k) = √ V k|i . We will limit our analysis to matrix elements which only depend on q, and assume that the electron energy levels are also spin independent, which allows the spin sums to be easily computed: where |M| 2 is the spin averaged matrix element squared and we have defined a crystal form factor f i→f , written in terms of both momentum and position space representations of the wave functions. The transition rate per target mass, R i→f , is then given by where ρ T is the target density, ρ χ = 0.4 GeV/cm 3 is the local DM density, and f χ is taken to be a boosted Maxwell-Boltzmann distribution. The total rate, R, is then simply the sum over all possible transitions from initial to final states. Since the only v dependence in Eq. (3.7) comes from the energy conserving delta function, we perform the v integral first and define g(q, ω) = 2π d 3 vf χ (v)δ(ω − ω q ). This integral can be evaluated analytically (see e.g. Refs. [15,23,52]): where ω = E f − E i is the deposited energy, and N 0 is a normalization factor such that d 3 vf χ (v) = 1. We take the DM velocity distribution parameters to be v 0 = 230 km/s , v esc = 600 km/s, and v e = 240 km/s. The total rate then becomes Here we focus on simple DM models, such as the kinetically mixed dark photon or leptophilic scalar mediator models. In these models M(q) can be factorized as M(q) = M(q 0 )F med (q 0 /q) f e /f 0 e , where F med (q 0 /q) = 1 for a heavy mediator and F med (q 0 /q) = (q 0 /q) 2 for a light mediator, and f e /f 0 e is a screening factor discussed in more detail in Sec. 3.5. As in previous works, we choose the reference momentum transfer to be q 0 = αm e . We can then finally write the rate in terms of a reference cross section, and find Another useful quantity is the binned rate (the rate for energy deposition between ω and ω + ∆ω), ∆R ω , defined as

Valence to Conduction
We begin with valence to conduction band transitions. The initial (final) states are indexed by band number, i(f ), and Bloch momentum, k i (k f ) inside the 1BZ. The wave functions in Eq. (2.1) can be substituted into the crystal form factor in Eq. (3.6), where the integral is over the primitive cell with volume Ω, and we have used the identity r e iq·r = N G δ q,G . The total rate in Eq. (3.12) is then is the number of valence (conduction) bands. This is identical to the rate formulae derived in [10,15,16] but written in terms of the periodic Bloch functions, u i,k (x), instead of their Fourier transformed components, u i,k,G , similar to Ref. [50]. Numerically the position space form is superior since the integral over the primitive cell can be computed by Fast Fourier Transform. This reduces the computational where N G is the number of G points, i.e. the number of Fourier components in the expansion of u i,k in Eq. (2.3).
In Fig. 5 we show the scattering rate from valence to conduction transitions binned in energy deposition, defined in Eq. (3.13), for a 1 GeV DM. The main difference between the calculation performed here and in previous works is the effect of the AE reconstruction, as discussed in Sec. 2.1.1. For the case of DM with a heavy mediator, the rate, even with experimental thresholds as low as ∼ 10 eV, is significantly enhanced relative to previous work. The AE reconstruction plays less of a role in the light mediator case since the transition rate is dominated by small momentum transfers. However, at high thresholds, where only larger momentum components can contribute, the AE reconstruction can still significantly boost the scattering rate by fully including the contributions neglected in the pseudo wave functions.
Since most earlier works computing DM-electron scattering include only valence to conduction transitions, it is useful to understand for which DM masses these are the only kinematically allowed transitions. If ω < E g − E core max , where E core max is the maximum energy of the core states, then the core states cannot contribute; if ω < E dft the free states are not available. Therefore if ω < min{E dft , E g − E core max } only the valence to conduction transitions are allowed, which can be related to a DM mass via ω max (m χ ) < min{E dft , E g − E core max }, where Requiring that ω max > E g , where E g is the band gap, sets a lower bound on the minimum detectable mass, m min χ , For Si (Ge), with a band gap of 1.11 (0.67) eV, m min χ is 0.28 (0.17) MeV. Lastly, we remark that for DM interactions characterized by higher-dimensional operators (not considered in this work), the scattering rate scales with higher powers of q and therefore is even more sensitive to AE reconstruction (and also c → c contributions discussed below in Sec. 3.3), which must be included in the analysis.

Valence to Free
For valence to free transitions the initial states are identical to those from Sec. 3.1, labeled by band number i and Bloch momentum, k i . The final state wave functions are simple plane waves given by Eq. (2.12), labeled by a momentum k f in the 1BZ with the bands labeled by G. We can therefore directly substitute Eq. (2.13) into Eq. (3.14) derived in the previous subsection, and obtain the crystal form factor: where u i,k i ,G are the Fourier components of the Bloch wave functions defined in Eq. (2.3).
Incorporating the Fermi factor correction discussed in Sec. 2.3, we find the rate in Eq. (3.12) is given by With a change of variables, G = G f − G and defining k ≡ k f + G f (and then dropping the prime for simplicity), the rate becomes In Fig. 6 we compare the binned rate from the valence to conduction (v→c) calculation in the previous subsection to the valence to free (v→f) one performed here, again for a 1 GeV DM. We see that for large ω, where the v→c calculation is limited by the number of conduction bands included, the v→f calculation extrapolates the results to higher ω as expected. There is some uncertainty due to the choice of the effective charge parameters, which is why the results are shown in bands. The lower edge corresponds to the conservative choice of Z i,k i eff = 1 for all i, k i , and the upper edge corresponds to the value set by the binding energy, Eq. (2.15) with E B = −E i,k i . We find that the conservative choice Z i,k i eff = 1 is a better match to the edge for the v→c calculation, and will use this in our final projections in Sec. 4. Note that as the threshold increases, the effect of v→f transitions becomes more important, and for a heavy mediator non-negligible constraints can be placed even with O(100) eV energy thresholds.

Core to Conduction
We now turn to core to conduction transitions. The initial core states are indexed by κ, the atom in the primitive cell, the usual atomic quantum numbers, n, l, m, and the Bloch momentum, k i . The final states are the DFT computed conduction states. The crystal form factor is simply obtained from Eq. (3.14) by substituting u i,k i → u κnlm,k i :

23)
The total scattering rate is then  where N a is the number of atoms in the primitive cell, N κ p is the largest principal quantum number for atom κ, and ω = E f,k f − E κnl . The core wave functions u κnlm,k i (x) are given by Eq. (2.9), and involves a sum over primitive cells. Since the integral in Eq. (3.23) is just over one primitive cell, only the atoms in this and neighboring cells can have a significant contribution. In other words, the sum over r converges very quickly due to the localized nature of atomic wave functions. We therefore restrict r to be summed over only the 3×3×3 nearest cells.
The contribution of core to conduction (c→c) transitions to the binned rate, for m χ = 1 GeV, can be seen in Fig. 7. In most cases the v→c transitions are dominant compared to the c→c, but there are two main scenarios where this is not true. First, when the experimental threshold is raised; this excludes the v→c transitions and causes the c→c contribution to be dominant. heavy mediator (bottom left panel of Fig. 7). If the experimental threshold is ∼ 50 eV the c→c contribution from the 2p states in Si gives the dominant contribution. Second, for a Ge target, and a DM model with a heavy mediator, the 3d states dominate the rate even at the lowest experimental threshold. To understand this in more detail we present Fig. 8 which compares the binned rate taking different modeling approaches for the 3d states in Ge. We see that the large momentum components of the wave function, recovered only after AE reconstruction in the DFT calculation, dominate the rate, which explains why previous works have underestimated the importance of 3d electrons. Meanwhile, we see explicitly at the scattering rate level that the semi-analytic approach accurately reproduces the DFT calculation at low ω, and extends the latter beyond its cutoff at high ω, consistent with the observation at the wave function level in Fig. 3.

Core to Free
The last transition type we consider involves a core electron initial state and a free electron final state. The crystal form factor is most easily obtained by substituting Eqs. (2.8) and (2.12) into its definition, Eq. (3.6):  Figure 9. DM-electron scattering rate from core states to conduction bands (c→c) and to free states (c→f) binned in energy deposition (with ∆ω = 10 eV) for 1 GeV DM, light (top row) and heavy (bottom row) mediators, assuming σ e = 10 −40 cm 2 . As in the v→f calculation in Fig. 6, the upper edge of the shaded bands corresponds to Z eff from Eq. (2.15), and the lower edge corresponds to Z eff = 1.
where the Fourier transform of the RHF Slater type orbital (STO) core wave functions, given in Eq. (2.10), are known analytically [71]: The direct detection rate is then where q = k f − k i + G, and ν κnl = ν(Z κnl eff , ω + E κnl ). We can now shift the G f variable, G ≡ G f − G and define k = k i + G and k = k f + G. Therefore, q = k − k and which is the closest expression to the vacuum matrix element, with just the inclusion of the core wave functions acting as a form factor. In Fig. 9, we compare the binned rate from the core to conduction (c→c) calculation to the core to free (c→f) calculation and see a reasonable extrapolation to higher ω. As with the transition region between v→c and v→f shown in Fig. 6, we find Z eff = 1 gives a better match between c→c and c→f. While the total number of electrons from these transitions is expected to be much less than lower energy transitions, this is the best available calculation for thresholds up to the kinematic limit of ω max .

In-medium Screening
DM-electron interactions mediated by a dark photon or scalar are screened due to the inmedium mixing between the mediator and the photon. The relevance of screening has been recently emphasized in Ref. [48]. The screening factor, f e /f 0 e , is related to the longitudinal dielectric, f e /f 0 e = (q · ·q) −1 , where is the dielectric tensor. It can be computed from in-medium loop diagrams or extracted from optical data. Here we model the dielectric of Si and Ge following Ref. [72]: (q, ω) = 1 + 1 and ij = (q, ω) δ ij . Here, 0 ≡ (0, 0) is the static dielectric, α is a fitting parameter, q TF is the Thomas-Fermi momentum, and ω p is the plasma frequency. The parameters used for Si and Ge are listed in Table 1, and we plot the dielectric as a function of q, ω in Fig. 10. Naively one might expect that the effect of the dielectric is to screen the rate by an O(100) factor due to the fact that the static dielectric, 0 , is O(10). However, this is only the value of the dielectric function at q = ω = 0, while as q → ∞ and ω → ∞ the dielectric approaches unity. Therefore, the effect of the dielectric crucially depends on the region of the kinematic phase space being probed. For a given energy deposition, ω, the momentum transfer is limited to q ω/v where v ∼ 10 −3 is the DM velocity. Therefore, the absolute minimum momentum transfer is q min ∼ E g /v ∼ O(keV), for O(eV) band gap targets. This is parametrically the same size as the Thomas-Fermi momentum q TF , so the dielectric is expected to slightly deviate from one, which causes only an O(1) shift to the scattering rate, as seen in Fig. 11.

Projected Sensitivity
We now compile the results from the previous sections to compute the projected sensitivity. We also compare the relative importance of each transition type, and discuss differences between our results and previous calculations in the literature. When there are large differences, it is typically because of the inclusion of AE reconstruction and core states in the calculation. Since AE reconstruction and core states contribute predominantly at higher momentum transfer and energy deposition, we will find the largest differences typically occur for a massive mediator and higher detector threshold, where the effects in some cases can be more than an order of magnitude (especially for Ge). For the case of a massless mediator and lower detection threshold, the differences with previous literature are much smaller and mostly due to the inclusion of in-medium effects.
In Fig. 12 we show the contribution to the binned rate from each of the four transition types, for a 1 GeV DM. We see that valence to conduction (v→c) has a higher peak than the other three transition types, except for the Ge, heavy mediator case, where core to conduction (c→c) has the highest peak. For comparison, Refs. [10,11] compute the valence to  Figure 11. Effect of screening on the binned rate (top row, for 1 GeV DM) and total rate (bottom row, as a function of m χ ) from v→c transitions for DM models with a light (red) and heavy (blue) mediator. The unscreened rate R no scr is obtained with = 1, and the screened rate R scr is obtained with the model of the dielectric function given in Eq. (3.31).
conduction rates with DFT, including also the 3d states in Ge, but without AE reconstruction. As expected, we find a lower rate at the lowest energy depositions due to the inclusion of in-medium screening, and a much higher rate at high ω due to AE reconstruction and inclusion of core states.
The impact of these observations on the reach depends on the energy threshold. Assuming charge readout (e.g. via a CCD), the relevant quantity is the number of electron-hole pairs, Q, produced in an event. For an energy deposition ω, this is given by  , valence to free (v→f), core to conduction (c→c), and core to free (c→f). We assume σ e = 10 −40 cm 2 , and take Z eff = 1 for all effective charges in the Fermi factor. Note that the c→c and c→f transitions involve semi-analytic treatment of 2p (3d) states and below in Si (Ge), which has been validated with DFT calculations including AE reconstruction; see Fig. 3. We also overlay the binned rate from Ref. [11] which computed the v→c contribution using QEdark (treating 3d states in Ge as valence, without including AE reconstruction effects).
where the values for ε are 3.6 eV and 2.9 eV for Si and Ge respectively. In Fig. 13, we show the total rate as a function of the DM mass, for Q ≥ 1, 5, 10. The threshold only affects the v→ c rate, as the other three transition types involve energy depositions corresponding to Q > 10, and are therefore always fully included. We see that for Q ≥ 1, the valence to conduction (v→c) contribution dominates the total rate with the exception of the Ge, heavy mediator scenario, where core to conduction (c→c) is dominant for m χ   Figure 13. DM-electron scattering rate as a function of the DM mass, for light (top row) and heavy (bottom row) mediators, from all four transition types: valence to conduction (v→c), valence to free (v→f), core to conduction (c→c), and core to free (c→f). We assume σ e = 10 −40 cm 2 , take Z eff = 1 for all effective charges in the Fermi factor, and show results for several threshold Q values which significantly impact the v→c contribution.
lower rate. We also see that v→f and c→f contributions are subdominant in all cases. Finally, we present the projected reach on the DM-electron reference cross section σ e in Figs. 14 and 15, for Q ≥ 1 and Q ≥ 10, respectively. Our new calculation yields several important differences compared to the previous literature, and we discuss them in detail in the following subsection.

Comparison With Previous Results
We begin by comparing to our previous work, Ref. [16], shown in brown in Fig. 14. We previously restricted our analysis to the light mediator scenario, and Q ≥ 1, which is relatively  unaffected by AE reconstruction effects since the rate is peaked at small energy/momentum transfers, as seen in Fig. 5. The main reason the reach here is weaker is the inclusion of in-medium screening discussed in Sec. 3.5.
Next we compare to Ref. [10], shown in red in Figs. 14 and 15. Those results were computed solely from valence to conduction (v→c) transitions. The largest discrepancy is in the high m χ regime scattering off a Ge target via a heavy mediator. This is due to high momentum contributions to the 3d wave functions in Ge. Ref. [10] computed the 3d states with DFT without AE reconstruction, which as we saw in Fig. 3 is crucial for recovering the dominant part of the 3d wave functions at high momentum. As discussed in Sec. 2.2, our modeling of 3d electrons in Ge as core states reproduces their DFT-computed wave functions up to the AE reconstruction cutoff, and provides a robust parameterization of higher momentum components. Since the valence states in Ge also contribute an appre-  ciable amount, the Q ≥ 1 results in Fig. 14 only differ by about an order of magnitude. However, the difference is more stark when going to higher Q thresholds in Fig. 15, which essentially isolates the 3d electrons' contribution. In the low mass regime the difference is less significant, and primarily due to the inclusion of screening effects. Another difference that is important here is sampling of the 1BZ. Ref. [10] used a uniform 6 × 6 × 6 mesh with extra 27 points chosen by hand close to the center of the 1BZ, whereas here (as well as in Ref. [16]) we use a uniform 10 × 10 × 10 grid. While checking convergence we found our (unscreened) results using a 6 × 6 × 6 uniform mesh were a closer match to Ref. [10]; generally, increasing the number of k points reduces the rate toward convergence, i.e. R 10×10×10 < R 9×9×9 < R 8×8×8 . This can be seen more directly in the difference between the brown and red lines in the light mediator scenario (as both are computed without screening), and it affects Ge more than Si, as is expected due to the smaller band gap and greater dispersions of nearby bands requiring denser k point sampling for convergence.
Ref. [3] also computed DM-electron scattering rates in semiconductors, focusing on Ge. The approach taken in that paper was to semi-analytically model the Ge wave functions with the core wave functions (with the same set of RHF STO wave function coefficients tabulated in Ref. [43]) and treat the final states as free with a Fermi factor, analogous to the core to free calculation performed here. As we can see from Fig. 14, while for most of the mass range and mediators the estimates are too optimistic due to incorrect modeling of the valence and conduction states, in the high mass region with a heavy mediator (bottomright panel), where 3d states dominate, their estimates are in good agreement with ours presented here, as expected.
Finally, we discuss the comparison with the most recent work, Ref. [48], which was limited to valence to conduction transitions. To show the effect of screening, we show their projected reach with (purple) and without (green) screening in Fig. 14. Again the largest discrepancy is in the heavy mediator scenario with a Ge target, primarily due to the neglect of the 3d states in Ref. [48]. When these are not important, i.e. the low mass regime or a light mediator, we generally find good agreement, with our reach being a bit stronger. Notably this does not seem due to a mis-model of the dielectric, since the effect of screening relative to our previous results, Ref. [16], is consistent with their result. We also find that screening has a smaller effect at high masses in the heavy mediator scenario for Si. These small differences are harder to disentangle since they could be due to: 1) different xc-functionals used (PBE and HSE vs. TB09); 2) local field effects which are only partially included here since we assume the screening factor is isotropic; 3) the plane wave expansion parameter, E cut , taken to be 500 eV without AE reconstruction in Ref. [48], vs. 1 keV, AE corrected to 2 keV taken here; 4) DM velocity distribution parameters, studied in detail in Ref. [73], for which Ref. [48] assumed v esc = 500 km/s as opposed to v esc = 600 km/s chosen here; and 5) Ref. [48] took a directionally averaged dielectric, whereas here we only assume isotropy in the screening factor but not the matrix element itself.

Conclusions
Dark matter-electron scattering in dielectric crystal targets, especially semiconductors like Si and Ge, are at the forefront of DM direct detection experiments. It is therefore imperative to have accurate theoretical predictions for the excitation rates. In this work, we extended the scattering rate calculation in several key aspects. Much of the focus of previous calculations has been on transitions from valence to conduction bands just across the band gap, which will be accessible to near-future experiments. We performed state-of-the-art DFT calculations for these states, and highlighted the importance of all-electron reconstruction which has been neglected in most previous works. Along with this, we extended the transition rate calculation by explicitly including the contributions from core electrons and additional states more than 60 eV above the band gap using analytic approximations.
We updated the projected reach with our new calculation and found important differences compared to previous results. In particular, we found that in the heavy mediator scenario, 3d electrons in Ge give a dominant contribution to the detection rate for DM heavier than about 30 MeV. Also, the rate can be significantly higher than predicted previ-ously for higher experimental thresholds. This is exciting because new DM parameter space will be within reach even before detectors reach the single electron ionization threshold.
We also release a beta version of EXCEED-DM (available here [53]) that implements our DM-electron scattering calculation for general crystal targets, and make the electronic wave function data for Si and Ge [54], as well as the EXCEED-DM output [55], publicly available so our present analysis can be reproduced. We have previously used EXCEED-DM for a target comparison study [16], and to study the daily modulation signals that can arise in anisotropic materials [15]. The generality of EXCEED-DM means that potential applications are vast. It can be used to compute detection rates for other target materials (assuming DFT calculations of valence and conduction states are available), and can also be adapted to include additional DM interactions such as in an effective field theory framework similar to the study of atomic ionizations in Ref. [5] (see Ref. [74] for a recent effort in this direction). For momentum-suppressed effective operators, a full calculation in our framework is even more important, as the effects of all-electron reconstruction and core states (overlooked in Ref. [74]) are generally amplified. Moreover, the differential information that can be obtained from our program facilitates further studies including realistic backgrounds. Details of EXCEED-DM and additional example calculations will be presented in an upcoming publication.