Fully ab-initio all-electron calculation of dark matter–electron scattering in crystals with evaluation of systematic uncertainties

,


I. INTRODUCTION
There has been rapid progress in direct-detection searches of sub-GeV dark matter (DM) by looking for electron recoils from DM-electron scattering in noble liquids and crystals (see, e.g., ).A theoretical description of such DM-electron scattering processes requires a quantitative description of the electronic structure of the detector material.This can be achieved by utilizing methods from condensed-matter physics and quantum chemistry.In particular, density density functional theory (DFT) has been demonstrated to be a powerful tool for determining the ground-state electronic structure in a wide variety of materials from first principles, as well their response to various perturbations [27].Also, DFT is the necessary first step to performing calculations using more advanced methods to treat systems with, e.g., stronger electron-electron correlations [28] Several choices need to be made when calculating the electronic structure and DM-electron scattering rate from DFT.This includes choosing the type of basis functions to describe the electronic wavefunctions; whether to separate the core electrons in the material from the valence electrons by, e.g., pseudopotentials [27] or the projector-augmented wave (PAW) method [29]; and the exchange-correlation (XC) functional, which incorporates the many body effects of the electron-electron interactions [27,30].These choices often represent tradeoffs between accuracy and computational efficiency [31], which differ for different types of materials.Moreover, once these choices are made, the relevant properties must be tested for convergence with respect to the numerical parameters of the calculation.
DFT-based techniques have been applied to calculating DM-electron scattering in a variety of materials relevant for detectors including semiconductors [5,[32][33][34][35], semimetals [36], superconducting nanowire single-photon detectors [37,38], quantum dots [39], etc.For the most part, studies of electron recoils in semiconductors have taken the approach of a plane wave basis set, core electrons frozen in pseudopotentials of PAW potentials, and either local/semilocal XC functionals based on the local density approximation (LDA), the generalized gradient approximation (GGA), or hybrid functionals that include a fraction of exact electron-electron exchange interaction [40].DFT using plane wave basis sets are the one most commonly used for studying solids in condensedmatter physics and materials science.
However, there are motivations for choosing an alternative approach.First, it was recently shown that explicit treatment of the core electrons significantly affects the DM-electron scattering rates [34,35,41].While sub-GeV DM does not typically excite an electron from a core orbital to the conduction band, the inclusion of the rapidly oscillating part of the valence all-electron wavefunction near the atomic cores (where it must be orthogonal to the core electron wavefunctions [42]) is necessary to capture DM-electron scattering events with high momentum transfer.For a calculation with a plane wave basis to be tractable, these oscillations must be smoothed through the use of pseudopotentials, though the all-electron wavefunction can be reconstructed if PAWs are used [34,35].
A basis set made up of localized functions, e.g., atomcentered Gaussians, can treat core and valence electrons on the same footing without significant increase in computational cost.In addition, such basis sets can be used either with periodic boundary conditions for solids [43] or for finite systems, such as molecules or nanostructures [43], increasing the flexibility to explore different DM detector materials.Finally, using localized basis sets allows the use of quantum chemistry methods, which allows many-body correlations to be included when calculating the wave functions for atoms, molecules, liquids, and solids.
In this work, we develop the computational methodol-ogy to perform all-electron calculations of DM-electron scattering based on localized Gaussian basis sets.The resulting code, which we call Quantum Chemistry Dark (QCDark), is based on the python-based simulations of chemistry framework (PySCF) [43][44][45] package, which allows for DFT and quantum chemistry methods to be used on both finite and extended systems . 1 Via benchmark calculations on silicon and germanium, we show that the basis sets can be converged.PySCF has previously been employed in the context of DM-electron scattering in isolated atoms and molecules in [46].We compare our results for DM-electron scattering with previous work, and find good agreement with EXCEED-DM [34,35], which reconstructs all-electron effects with PAW reconstruction.We also quantify uncertainties related to the choice of XC functional and numerical convergence parameters.
The rest of the paper is organized as follows.In §II, we describe electronic structure calculations, include a comparison between plane-wave and atomic-centered bases, discuss how to include all-electron effects, and describe the DM-electron scattering rate calculations with quantum chemistry basis sets.In §III, we describe the results for Si and Ge, including the systematic uncertainties and the effects of the secondary ionization modeling; we also calculate the annual modulation rates, and compare our results with those of previous works.We conclude in §IV.Two appendices contain technical information, including the properties of Cartesian Gaussians (Appendix A) and a derivation of the scattering rate formulae (Appendix B).

II. CALCULATING DARK MATTER-ELECTRON SCATTERING RATES
In this section, we introduce the computational approach that underlies QCDark, comparing the all-electron localized basis set approach used in this work with previous implementations based on plane waves and pseudopotentials.

A. Electronic structure
The electronic structure of the material is described by Kohn-Sham (KS) [30] DFT, in which the equations that describe the system of many interacting electrons are mapped onto a set of single-particle equations in an effective potential constructed to reproduce the groundstate electron density and total energy.The energy functional of the density n is given by (Gaussian units are 1 The code is available at https://github.com/asingal14/QCDark.assumed throughout) where T S is the sum of the kinetic energies of the noninteracting orbitals; v ext is the external potential given by, e.g., the atomic nuclei in an all-electron calculation or the ions if pseudopotentials are used; and E xc is the so-called exchange-correlation energy, which accounts for the many-body and quantum effects that are neglected in T S and the Hartree electron-electron interaction (third term in Eq. ( 1)) [27].Writing the density as a sum over auxiliary single-particle orbitals, i.e., n(r) = i |ψ i (r)| 2 , and minimizing Eq. ( 1) with respect to variations in ψ results in the KS equations, where /δn(r), with the last term being the functional derivative of the exchangecorrelation energy with respect to the density, which results in the XC potential v xc (r).The KS equations in Eq. ( 2) must be solved self-consistently, since the effective potential depends on the density.
Several choices exist in the above calculation, which we discuss in subsequent sections.First, the KS wavefunctions ψ must be expressed in terms of some basis functions ( §II A 1); the choice of these basis functions has a significant impact on many other aspects of the calculation, including the possible choices for the boundary conditions and the treatment of core electrons ( §II A 2). Also, the exact form for the exchange-correlation potential is not known, and the choice of approximate v xc (r) can result in qualitatively different results for the electronic structure of the material ( §II A 3).

Basis sets
Previous works [5,[32][33][34][35] used plane waves as basis function (PW-basis), together with periodic boundary conditions.Then the KS wavefunctions for a given wavevector k in the first Brillouin Zone can be written as where K is a reciprocal lattice vector, e i(k+K)•r are the basis functions, u i (k + K) are the coefficients, and V is the volume of the crystal.The accuracy of the PW-basis for describing ψ is governed by the number of reciprocal lattice vectors included in the basis, which is usually specified as a kinetic energy cutoff of the plane waves.
The PW-basis is good at describing relatively delocalized states, e.g., the states near and above the Fermi level in most solids; however, as we will discuss in §II A 2, it is computationally expensive to capture the core electrons with a PW basis, since their localized nature requires the inclusion of very high energy plane waves and hence a large PW-basis size.
In this work, we use atom-centered Cartesian Gaussian basis sets.These basis sets are efficient at treating localized states, including the core electrons, and may be used for periodic or finite systems.The building blocks of this basis are primitive Gaussians, where A is the atomic position, ξ µ is an adjustable parameter, and the Cartesian exponents i, j, and k are all integers that satisfy the condition i + j + k = l, with l being the angular quantum number of the shell.Then, the basis functions are given by contracted Gaussians, i.e., weighted sums of N prim primitive Gaussians: where α = {κ, ijk} is a composite index that runs over all nuclei κ in the system (or unit cells for periodic boundary conditions) as well as the Cartesian exponents, N µ is the normalization of the primitive Gaussian, and c µ is the coefficient of the primitive Gaussian.Note that N µ and c µ are fixed throughout the DFT calculation.For periodic calculations, the atomic orbitals in the unit cell are then where R are the real-space lattice vectors.Upon selfconsistently solving the DFT hamiltonian, we obtain a coefficient matrix for each k, C iα (k), so that the Kohn-Sham wavefunctions (often referred to as molecular orbitals) are where N cell is the number of unit cells in the crystal.Compared to plane waves, Gaussian basis sets are significantly more complex.First of all, Gaussian basis sets are element-specific and the total basis set will be the sum of those from the individual atoms.For a given element, the construction of the basis set requires defining c µ , N µ , and ξ µ for each atomic shell n and angular momentum l = i + j + k.The choice of basis presents a source of systematic uncertainty in our calculation, which we estimate by varying the basis sets for each system.
The size of the basis set determines the number of orbitals, with multiple important parameters to consider.Moreover, there are several naming conventions that are widely used.In this paper, we use several basis sets, all taken from [47].The number of "zetas" in a basis set refers to the number of orbitals as a function of the number of occupied orbitals.An N -zeta (NZ) basis set would have N basis functions for each fully or partly occupied orbital in an atomic species.Si has an electronic configuration of 1s 2 2s 2 2p 6 3s 2 3p 2 (where the 3s 2 3p 2 are valence electrons), and so has 3 s-orbitals and 2 × 3 p-orbitals, and so a DZ (double zeta) basis set would have 6 s-and 4 × 3 p-orbitals.The treatment of l ≥ 2 orbitals differs between spherical and cartesian gaussians.This is because there exists a linear combination of, for example, three d-orbitals in Cartesian Gaussians (ijk) = (200), (020), and (002) orbitals, which mimics (in this case) an s-orbital.In addition, certain basis sets contain N zeta only for the valence orbitals, with one atomic orbital for core states.This is called a NZ-valence or NZV basis set.One can further add polarization functions on the valence orbitals, which adds new orbitals.For example, for Si, a DZP (double-zeta polarized) basis set would add a d-orbital to allow for the electrons to be polarized in the atoms.In general, addition of extra orbitals and polarization functions improves convergence, though diffuse components in basis sets may become pseudo-linearly dependent in periodic calculations [43,45].
It should be noted that the aforementioned details represent a very brief overview of the very complex field of quantum chemical basis sets [48], and different basis sets, especially with different naming conventions can add complexities to these considerations.
We employ the TZP and the def2-TZVP basis set for Si and Ge respectively, and show the uncertainties associated with the choice of basis set in §III C 1.

Treatment of core electrons in DFT
The treatment of core electrons, i.e., those tightly bound to the nuclei, is closely connected to the choice of basis set.For atom-centered Gaussians, it is straightforward to include all of the electrons in the DFT calculation, since localized basis functions can just as easily describe core electrons as those near and above the Fermi level.However, for plane waves, describing such localized wavefunctions would require a prohibitively large energy cutoff.In addition, including core orbitals would require the wavefunctions of the valence electrons to be orthogonalized to the wavefunctions of the core electrons; this would introduce rapid oscillations in the region around the atomic nuclei, which would also require too many plane waves to be computationally tractable.Therefore, plane-wave calculations usually freeze the core orbitals via pseudopotentials, effective core potentials (ECPs), or the projector-augmented wave (PAW) method.This has three main effects.First, the core orbitals do not participate in hybridization and bonding in the crystal; this is usually an excellent approximation, since such orbitals are so tightly localized around the atoms, so that there is negligible overlap between atoms.Second, since they are not explicitly included, DM-electron scattering transitions between core orbitals and the conduction band are neglected; this is also not usually an issue, since the energies of such transitions is often beyond the scope of light DM searches.Third, and most crucially, the "pseudowavefunctions" are smooth in the core region, since they no longer must be explicitly orthogonalized to the core orbitals.
As was shown in [34,41] (see also §4), all-electron effects are crucial for describing DM-electron scattering events with high momentum transfer, since such events couple high-frequency modes, which are only present in the all-electron valence and conduction bands due to the rapid oscillations near the core.In [34], the all-electron wavefunctions were recovered after a plane-wave calculation via the PAWs.In §III, we benchmark this methodology against the full all-electron calculation allowed by our Gaussian basis set.

Exchange-correlation functional
In practice, the exact form of the exchange and correlation energy in Eq. ( 1) is not known; however, there are many well-motivated approximations (see [49]).The choice of exchange and correlation functionals presents a major source of systematic uncertainty in the calculation.This systematic uncertainty can be estimated by calculating the electronic structure with different functionals, and comparing the results.The broad categories of functionals are 3. Meta-GGAs (mGGA) contain higher order derivatives of n, including terms like ∂ • ∂n and ∇ 2 n.
4. Hybrid functional are GGA and mGGA functionals with added exact Hartree-Fock exchange.
5. Double hybrid functionals add Møller-Plesset perturbation theory at second order ("MP2 level") to hybrid functionals in an effort to better model correlations.
We use the well tested PBE0 hybrid exchangecorrelation functional.We further test several GGAs, mGGAs, and other hybrid functionals to show the dependence of DM-electron scattering rates on the choice of exchange-correlation functionals.Moreover, the choice of E xc affects the calculated bandgap E gap of the material, and so we apply a scissor correction to match the experimental bandgap.

Theory
As described in §II A, we self-consistently solve the Kohn-Sham equations to obtain the coefficients C iα (k) in Eq. (7).The key quantity required from the DFT calculation for calculating DM-electron scattering is the crystal form factor (equivalent to Eq. (3.17) of [5]), where E ik is the energy of the i th orbital from ground state at k in the reciprocal cell and U = k + K − k.
Here, E e and q are the energy and momentum transferred from the DM particle to the electron, and θ v and φ v refer to the polar and azimuthal angles of v, respectively.Indices i and j run over occupied and unoccupied orbitals, respectively, and The DM-electron scattering rate for a DM particle of mass m χ and local density ρ χ is where is the DM-electron interaction cross section assuming a bosonic mediator with mass m V .In general, we assume two limits, m V q (called the heavy-mediator limit) and m V q (called the light-mediator limit).The factor f e /f 0 e 2 is a screening factor discussed next (see also [34]), while η(v min (q, E e )) is the average inverse speed of DM in the galaxy (for more details, see [5] or Appendix B) and We describe the calculation of the matrix elements, Eq. ( 9), in Appendix A.
If the DM-electron interaction is mediated by a dark photon or scalar, the interaction is screened due to the in-medium effects [36,50].Recent works [32,51] have emphasized the importance of this electrostatic screening, especially for recoils at low energy transfer.In the results shown below, we follow the prescription of [34,35], and multiply the crystal form factor, |f crystal | 2 by a fac- Here is the dielectric function, which is modelled as [52] (q, E e ) = 1+ where 0 ≡ (0, 0) is the empirically measured static dielectric constant, τ is a fitting parameter (we use the results from [52], which fit their dielectric function to the results of [53]), ω p is the plasma frequency, and q TF is the Thomas-Fermi momentum, with values listed in Table I.This equation ignores the tensorial nature of , since the crystals we consider here are cubic.This equation only estimates the real part of the dielectric function, and we do not account for the imaginary part in our screening.
Ref. [35] compares the differences in the expected DMelectron scattering rates using the analytical model in Eq. ( 11) for the dielectric function with an RPA calculation of the dielectric function.We expect an O(10%) correction for creating 3 or fewer electrons-hole pairs (Q ≤ 3e − ) with a numerically calculated dielectric function using a preliminary RPA dielectric function calculation, and O(1%) or less for creating more than 3 electronhole pairs.In practice, the numerical integration in k-space requires the replacement where N k is the number of k−points in the chosen k−grid.
Moreover, to discretize |f crystal (q, E e )| 2 in q and E e , we define bins with width ∆q and ∆E e and bin centers {q n } and {E m e }, respectively.Then, the discretization procedure follows as Consequently, the numerical crystal form factor is Here, as for Eq. ( 8), the recoil electron is assumed to be transferred from occupied molecular orbital |i, k to unoccupied orbital |j, k + K .
In principle, there should be n k = N T points modelled  14)) for silicon and germanium, respectively.The details of the calculation are described in §III A. The region beneath the black line is kinematically inaccessible for halo DM, as it would require vmin(q, Ee) > vEscape + v Earth (see Appendix B for more details).
in the calculation, where N T is the number of unit cells in this crystal, with N T O 10 23 for a 10 g crystal.However, such a dense mesh is impractical.Since we need to calculate rates for transition from each k to each k , the computational complexity scales as n 2 k , and the calculation quickly becomes unfeasible.
We generally model the grid in reciprocal space as a Γcentered Monkhorst-Pack grid, with n k,i points along the b i reciprocal lattice vector.We use the shorthand n k,1 × n k,2 × n k,3 to denote the mesh in the reciprocal space, which has n k = 3 i=1 n k,i k−points.Because both Si and Ge crystallize in FCC cells, it is reasonable to set n k,1 = n k,2 = n k,3 .We choose n k,i = 4 and 6 for Si and Ge, respectively, but check the dependence of DM-electron scattering rates on our choice of k−grid in §III C 3.
In addition, it becomes computationally expensive to include a large number of K vectors, and so in practice, we must limit the vectors to some q max .The choice of q max is another source of systematic error.We choose q max = 25αm e and q max = 20αm e for Si and Ge, respectively, and further show the dependence of DM-electron scattering rates on q max in §III C 5.

III. RESULTS & DISCUSSION
In this section, we describe the results from our calculation of the crystal form factors (Eq.( 14)) and the DM-electron scattering rates (Eq.( 10)) for Si and Ge performed with QCDark.We evaluate the various systematic uncertainties, and compare our results to those from other available codes.Finally, we look at the annual modulation rate as a function of the DM mass, m χ and the DM form factor, F χ .
Table II lists the values of crystal parameters used for our calculation of DM-electron scattering rates in both Si and Ge, including the experimental bandgap used for the scissor correction procedure.

A. Crystal Form Factor
The calculated crystal form factors using Eq. ( 14) are shown in Fig. 1, for silicon and germanium crystals in the left and right panels respectively.The region below the black lines are kinematically inaccessible for halo DM, i.e., the halo DM, irrespective of the mass of the fermion, is unable to transfer energy E e with a momentum transfer q < E e /(v Escape + v Earth ).
Both panels show an enhancement of the crystal form factor at q 4 αm e compared to Fig. 5 of [5], which is a consequence of including all-electron effects in our calculation.The enhancement in the crystal form factor for germanium crystals at E e 30 eV corresponds to the transitions from the semi-core 3d-shell to the conduction bands, which has an energy of −28.6 eV relative to the top of the valence bands in our calculation.

B. Dark matter-electron scattering rates
Fig. 2 shows the DM-electron scattering rates expected in silicon and germanium crystals for m χ = 10 MeV, 100 MeV, and 1 GeV for both heavy and light mediators, assuming the crystal form factors shown in Fig. 1.In this and subsequent figures, we plot where p(Q, E e ) is the probability that a transition with recoil energy E e excites Q electrons (for more details, see §III D).We use the ionization model at 100 K from [54] for Si.For Ge, we use an electron-hole-pair creation model, where Θ(x) is the Heaviside step function, E pp is the electron-hole-pair creation energy (E pp = 2.9 eV for Ge), and E gap is the bandgap of the material.
Panels (a) and (b) of Fig. 2 show DM-electron scattering rates in a silicon crystal mediated by a heavy and a light mediator, respectively.For Si, we use a TZP basis set with a PBE0 exchange-correlation functional, 4 × 4 × 4 k−grid and q max = 25αm e .Panels (c) and (d) of Fig. 2 show DM-electron scattering rates in a germanium crystal mediated by a heavy and a light mediator, respectively, calculated using a def2-TZVP basis set with a PBE0 exchange-correlation func- tional, 6 × 6 × 6 k−grid and q max = 20αm e .
Because the 3d-dominated bands in Ge are flat in k−space (i.e., they are highly localized in real space), we need a denser k−grid to reduce the numerical noise in the (unbinned) rate spectra dR/dE e .An ionization model for Ge akin to the model in [54] for Si (which includes a Fano factor) remains unavailable, which would smooth out the numerical noise while calculating ∆R Q (see §III D).
In principle, increasing the density of the k-grid would reduce the noise, at the expense of computation time.However the rates, barring ∼ 10% systematics coming from the numerical noise at high Q, are robust (see §III C 3 and Fig. 6 for more details), and indicate that for m χ O(100 MeV) the rates are higher for Q ≥ 11 e − than for lower Q (for lower masses, interactions with high q transfers are kinematically suppressed).This would imply that germanium-based detectors with relatively high thresholds can still probe significant regions of DM parameter space, assuming m χ 100 MeV and a heavy mediator.
The effects of all-electron modes are visible for heavy mediators, and not as much for light mediators.This is because of an effective |F χ | 2 ∝ 1/q 4 suppression in the DM-electron scattering cross-section in the case of light mediators.

C. Evaluation of Systematic Uncertainties
There are multiple sources of theoretical uncertainties as well as several convergence parameters (i.e., parameters that can be improved with more computational time) in our calculation of DM-electron scattering rates.The choice of the exchange-correlation functional, E xc [n], in Eq. ( 1) is a source of theoretical uncertainty.The realspace cut-off for constructing our Bloch atomic orbitals, the size of the k−grid, and the choice of q max are convergence parameters.The choice of the atomic centered Gaussian basis set is both a convergence parameter (since increasing the number of basis functions allows us to model the conduction states better) and a theoretical uncertainty (since different basis sets are optimized for different types of calculations, be it molecular or periodic boundary conditions).In this section, we go through each of these choices and determine their effects on our DM-electron scattering rate calculation.

Choice of basis set
There are many choices of atom-centered Gaussian basis sets available for use [47].However, most of these basis sets are optimized for molecular calculations, and we have to choose among the few optimized for periodic boundary conditions.In addition, the size of the basis sets determines the number of conduction bands.
Fig. 3 shows the DM-electron scattering rates calculated for various basis sets.While the DM-electron scattering rates are consistent across all the basis sets we test, the cc-pVQZ (correlation-consistent polarized valence quadruple zeta) and cc-pVTZ (correlation-consistent polarized valence triple zeta) are the best optimized basis sets that we test for Si and Ge, respectively (these are also computationally very expensive).For Si, DM-electron scattering rates calculated using TZP (shown in Fig. 2) differ by only ∼ 5% from those derived using cc-pVQZ.Similarly for Ge, DM-electron scattering rates calculated using def2-TZVP (shown in Fig. 2) differ by only ∼ 5% on average from those derived using cc-pVTZ.Both of For Si, we use here a TZP basis set, at 4 × 4 × 4 k-grid and PBE functional to calculate these results.For Ge, we use here a def2-TZVP basis set, with the same k-grid and PBE functional.Note that a good description of high momentum transfer q 3αme does not require a large Rcut for either element.the TZP and def2-TZVP basis sets provide a good balance of computational efficiency and accuracy, and we use these in further analyses.
Fig. 3 also shows that DM-electron scattering rates in Ge are heavily dependent on the choice of basis set, especially for large Q 9 e − .This is because basis sets like 3-21G and DZP are unable to capture conduction bands well, while the energy of the semi-core 3d electrons is highly dependent on accurate modelling of core shells.

Real space cutoff
Our atomic orbitals are Bloch sums in real space as in Eq. ( 6); in principle, one must sum over an infinite number of Gaussians displaced by real-space lattice vectors R to form each atomic orbital.In practice, however, Gaussians are rapidly decaying functions, and so it suffices to include a finite number of neighbors depending on the exponents in the contracted Gaussians.This is generally accomplished by setting a real space cut-off, which we call R cut .
While PySCF is capable of choosing a dynamic real space cut off for each orbital, which lowers the computa- tional cost to calculating the matrix elements (see Eq. 9), this complicates our analytical approach to calculating the matrix elements (see Appendix A for more details on an analytical calculation).Hence, we choose a constant R cut for all orbitals, the value of which is chosen via the following procedure.We first calculate the electron loss function, i.e., the imaginary part of the inverse dielectric function Im − (ω, q) −1 .For this, we assume that the real part of the dielectric function is modelled by Eq. ( 11), and calculate the imaginary part of the RPA dielectric function using Eq. ( 16) of [55].We then integrate ωIm − (ω, q) −1 over energy for a given magnitude of q; by the f -sum rule, the result should equal π/2ω 2 p , where ω p is the plasma frequency of the material.The f -sum rule is only achieved in limit of a complete basis set, however we have found that the convergence of this quantity is a useful diagnostic as to whether a given R cut is sufficient for the relevant range of q [51].We show the electron loss function integrated over E e up to 50 eV for Si (using the TZP basis set) and Ge (using the def2-TZVP basis set) in Fig. 4.
The high momentum transfer modes q 3αm e are q max = 8αm e q max = 10αm e q max = 15αm e q max = 20αm e q max = 25αm e 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Heavy Mediator q max = 8αm e q max = 10αm e q max = 15αm e q max = 20αm e q max = 8αm e q max = 10αm e q max = 15αm e q max = 20αm e q max = 25αm e captured well even by low R cut .One may understand this in the real space from the viewpoint of our KS wavefunctions -the high frequency modes of the molecular orbitals, resulting from orthogonalization with inner orbitals, are more localised near the nuclei, and hence have smaller overlaps in real space with counterparts from more distant atoms.The low q modes, on the other hand, correspond to long distance behavior of the matrix elements, and so necessitate the usage of larger R cut in order to satisfy the f -sum rule.The curves going to zero at very low q 0.3 αm e is an artifact of our choice of a sparse k−grid, and is not the true behavior of f sum rule.
Motivated by these results, our final calculations use a hybrid real space cut-off.For Si (Ge), we use R cut = 4 (3) cells for q ≤ 3αm e and R cut = 2 (1) cells otherwise.Fig. 5 shows that this is also cautious, as the low momenta deviation only occurs in a prohibited region of the parameter space (v min (q, E e ) > v Escape + v Earth ) for DM-electron scattering (see Appendix B for more details).

Convergence of k−mesh in reciprocal space
A potential source of systematic error in the rate calculation comes from the density of k−points in the first Brillouin Zone (1BZ).As discussed above, because the computational cost scales as the square of the number of k−points, it is infeasible to include N T k−points, where N T O 10 23 is the number of unit cells in the crystal.In this section, we discuss the effects of modelling the 1BZ with an N × N × N k−grid, with a total of N 3 k−points.
Fig. 6 shows the convergence of our calculations with k−grid for both Si and Ge.Note that Si is already converged at a 4 × 4 × 4 k−grid.One reason for this is that the numerical uncertainties are smoothed out from applying the ionization model (see §III D for more de-tails).For Ge, accurately describing transitions from the relatively dispersionless 3d-derived band at E e ∼ 29 eV below the Fermi level to the conduction bands requires a finer k−grid.We find that a 6 × 6 × 6 grid performs adequately, with errors of 15% compared to 8 × 8 × 8 in each bin.A probabilistic ionization modelling for Ge, akin to 54 for Si, will aid in smoothing out the recoil spectrum.

Exchange-correlation functional
Fig. 7 shows DM-electron scattering rates in Si and Ge crystals for various exchange-correlation functionals.We test the commonly used PBE GGA functional, along with SCAN and TPSS mGGAs.For Si, we test multiple hybrids -PBE0, SCAN0 and TPSS0, along with a hybrid semiempirical functional optimized for molecules rather than crystals (B3LYP).For Ge, we test PBE, SCAN, TPSS, and PBE0.
It is important to note that a scissor correction has been applied to the band gaps of Si and Ge, so DMelectron scattering calculations have the same gap regardless of functional.For materials where the experimental gap is not known, the differences in gaps predicted by different functionals is expected to lead to a significant source of variation in the scattering rates.A related issue observed for Ge (right panel of Fig. 7) is the dependence of the energy of the 3d shell.This results in significant differences in the DM-electron scattering rates in the 8-11 electron-hole-pair bins.It is apparent that PBE, TPSS, and SCAN underestimate the electron binding energy for the 3d-shell electrons, with values ∼25 eV from the top of the valence band.The PBE0 functional results in values between 28.6 and 29.0 eV, which are much closer to the experimental values of ∼29.5 eV of 3d-shell electrons, respectively [56].For Si, we scissor correct the bandgap, and core orbitals do not get involved until energies of ∼ 99.2 eV [56].

Maximum momentum transfer
The implementation of atom-centered basis sets without using an effective core potential has one direct effect-we are able to capture the high-momentum transfer regime of the crystal form factor.These high-q contributions come from orthogonalizing the valence and conduction bands against the core orbitals, which introduces high wavenumber modes to the valence and conduction wavefunctions.This allows the wavefunctions to be modelled to arbitrarily high wavenumbers, and allows us to fully capture the crystal form factor. Fig. 8 shows the impact of adding high q modes on the rates of 1 GeV DM particle interacting with Si (left panels) or Ge (right panels) via a heavy (top row) or light (bottom row) mediator.16), E pp = 3.8 eV Si m χ = 1 GeV, σe = 10 −39 cm 2 Heavy Mediator 4x4x4 k-grid 6x6x6 k-grid FIG. 9.The effects on the DM-electron scattering rates in silicon of using the secondary ionization modeling from [54] ("R&K") versus the step-function model from Eq. ( 16).We use the PBE exchange-correlation functional and qmax = 10 αme.
As expected, DM-electron scattering mediated by a light mediator is not significantly influenced by the high-q contributions in silicon.This is due to the |F χ (q)| 2 ∝ q −4 dependence of the rate in the integrand of Eq. (10).Ge, on the other hand, is sensitive to q max even for scattering through a light mediator, since the 3d-shell in germanium dominates the high q regime.
For interactions mediated by a heavy boson, however, there are important high-q contributions even for relatively small charge bins.Moreover, when including the high-q contributions, we see that DM with m χ 50 MeV and scattering through a heavy mediator (F χ ≈ 1), the rates from Q ≥ 11 e − bins dominate over the 1 e − ≤ Q ≤ 10 e − bins.Similarly, for Si with a heavy mediator, high q contributions are important, with, e.g., a ∼ 75% increase in rates if we go from q max = 8αm e to q max = 25αm e for the 10 e − bin.

D. Effects of the secondary ionization model for silicon
Our results for the DM-electron scattering rates in silicon are shown using the ionization modelling from [54].In Fig. 9, we compare these rates with those from a simple step function model from Eq. ( 16) for Si.For the latter, we use E gap = 1.1 eV and E pp = 3.8 eV (see, e.g., [57]).We see significant differences between the two ionization models for the 1 e − and 2 e − -bin, although the rates are similar for the bins with Q ≥ 3.As we observed for Ge, for which only a step-function model is available, the probabilistic model from [54] smoothes out the numerical fluctuations introduced by the sparse k−grid.10. modulation amplitude, f mod from Eq. ( 18), versus Q for Si (left) and Ge (middle), for mχ = 10 MeV and 1 GeV, for both heavy and light mediator-mediated scattering.The right panel shows f mod versus mχ, for Q = 1e − and Q = 5e − , for both Si and Ge.We use the PBE0 exchange-correlation functional, along with qmax = 25 αme (20 αme) and a 4×4×4 (6×6×6) k−grid for Si (Ge).11.This plot shows the reach for both, obtaining 2.3 events for 1 kg•yr exposure of our target material (solid line), as well as the threshold for a 5−σ discovery by annual modulation with the same exposure (dash-dotted line).Panel (a) shows the upper bounds on cross-section that can be placed for heavy mediators, while panel (b) shows the same for light mediators.We assume no background for either panel, and include 1 ≤ Q ≤ 10 e − for Si (black) and 1 ≤ Q ≤ 15 e − for Ge (blue).

E. Annual Modulation
The DM-electron scattering rates are dependent on the DM flux incident on the target material, which in turn depends on the velocity of the detector in the galactocentric frame.For table-top experiments, there are three major contributions to this velocity.First, there is the local circular velocity, which we take to be v 0 = 230 km s −1 .The second contribution comes from the Sun's peculiar velocity, v − v 0 .We use the recommended value, v = 250.2km s −1 [58].Finally, the earth revolves around the sun with an average speed v ⊕ = 29.8km s −1 .This revolu-tion causes an annual modulation in DM-electron scattering rates, as the total velocity, v Earth = |v + v ⊕ |, varies from 220.4 km s −1 on December 2 to 280 km s −1 on June 2.
We calculate the modulation amplitude following [5], where ∆R Q,0 = ∆R Q,Sept.2= ∆R Q,Mar.2 .Even in the presence of backgrounds, a measurement of f mod could allow for the detection of DM in such an experiment.We plot f mod as a function of the DM mass, the mediator form factor, and the charge ionized in the target material Q for both Si and Ge in Fig. 10.
Comparing to Fig. 8 of [5], the difference in rates from the inclusion of high wavenumber modes in the crystal form factor allows the electron to scatter into a larger parameter space, which reduces f mod , especially for the m χ = 1 GeV case.The same is visible for Ge, and is in fact even more pronounced for the 3d-shells, which dominate the rate.
A measurement of the annual modulation signal will be an important step in confirming a potential DM signal.We calculate the 5σ-sensitivity by requiring where ∆S = f mod S tot is the modulation amplitude, S tot is the total number of signal events, and B is the number of background events.Here f mod is calculated using Eq. ( 18), except we sum the rates over 1 e − ≤ Q ≤ 10 e − for Si, and 1 e − ≤ Q ≤ 15 e − for Ge.Assuming no background events, we show the 5σ discovery reach in Fig. 11 with dash-dotted lines for Si (black) and Ge (blue).The left panel shows the reach for heavy mediators, while the right panel corresponds to light mediators.

F. Comparison with other codes
Quantum Espresso was used in the first numerical calculation of the crystal form factor for Ge in [1].It was also used in [5], which presented a detailed calculation of the crystal form factor for both Ge and Si, and made the resulting code, QEDark, publicly available.In Fig. 12, the crystal form factor is recalculated with QEDark with improved computational parameters, including a higher energy cutoff for plane wave calculations, a denser k−grid for both Si and Ge, and the analytical screening described in §II B 1. We use the PBE functional, which, as shown above, underestimates the energy of the 3d-shell in Ge, and also excludes the effects of high frequency modes that are visible in both Si and Ge, especially for heavy mediators.
DarkELF [32] emphasized the need for better screening, especially for low energy excitations.Here we use GPAW RPA dielectric function with Local Field Effects (LFE) for both Si and Ge.However, it also does not include the high frequency modes, which dominate the rates at high energy, and the Ge 3d-shell is frozen in the pseudopotential, which otherwise dominate the DMelectron scattering rates at E e 30 eV.
EXCEED-DM [35] is able to reconstruct the highfrequency modes, and is also able to capture the dielectric screening with an RPA dielectric function.It also employs the well-tested and commonly employed HSE06 functional.However, it reconstructs the semicore and core orbitals after a pseudopotential calculation.We find good agreement between the rates calculated with EXCEED-DM and QCDark, showing that the PAW method is accurate in these materials.
QCDark implements ab-initio calculation of the crystal form factor, along with an analytical approximation to the dielectric function.However, as [35] recently showed, the analytic screening only approximates the true screening, and does not capture all the effects completely.On the other hand, QCDark allows for a much better handle on systematics by giving users control over the theory parameters, as discussed in §III C.

IV. CONCLUSION
In this paper, we present DM-electron scattering rates in silicon and germanium crystals calculated using a new code, which we make public as QCDark.We use a novel approach that naturally includes all core electrons, and treats them on the same level as valence electrons of the crystal.This implies that all-electron effects are automatically included from the beginning.Moreover, we present a systematic treatment of the theoretical uncertainties associated with the calculation, including those associated with DFT (basis set, exchange-correlation functional, and k−grid), along with uncertainties associated with the transition matrix elements (real space cutoff and the maximum momentum transfer modelled, q max ).
The major sources of systematic error include the choice of basis set and exchange-correlation functional, even after we apply the scissor correction, and the choice of q max , though the rates converge quickly in the E e ∈ [0, 50 eV] range for both Si and Ge.The rates also converge quickly as finer k−grids are chosen, assuming the secondary ionization model in [54] for Si.The rates also converge quickly for small values of the real-space cutoff, R cut .
We find that modelling high momentum transfers by including all-electron effects is necessary for accurately modelling DM-electron scattering rates, especially at high recoil energies, in line with the findings in [34,35].This is especially important for Ge crystals, in which the transition rates from the 3d-shell (when kinematically accessible) dominate the rates if the high momentum transfer modes are modelled accurately.and summed over final spins.
For bound electron initial and final states, say ψ 1 and ψ , respectively, the cross section is modified as where Moreover, because there is only one electron final state being considered, we can make the replacement (B4) For non-relativistic scattering, (B7)

Average rate in a DM halo
The rate of the specific transitions induced by DM hitting a target electron is then where n χ and g χ (v) are the DM number density and velocity distribution, respectively.In this work, we use the parameters recommended by [58].Note that the velocity distribution of DM in the standard halo model implies that the speed of the DM wind we observe must follow v χ < v Escape + v + v ⊕ = v Escape +v Earth , where v Escape is the escape velocity at the Sun's location in the galactic gravitational potential well, v is the sun's galactocentric speed, v ⊕ is the earth's heliocentric speed, and v Earth is the Earth's galactocentric speed.
In this paper, we assume both DM velocity distribution and electron wavefunctions to be spherically symmetric, which is not true in general.We then use the integral over DM velocity to eliminate the δ-function in Eq. (B7), obtaining (B9) Here v min is the minimum velocity of the DM particle required for an energy-momentum transfer of (E e , q) to be feasible, v min (q, E e ) = E e q + q 2m χ .(B10) We define η (v min (q, E e )) ≡ and obtain R 1→2 = n χ σe 8πµ 2 χe d 3 q 1 q η (v min (q, E e )) × |F χ (q)| 2 |f 1→2 (q)| 2 . (B12)

Excitation rates in crystals
So far, we have not discussed the initial and final states of the electron being acted upon.In a crystal system, an electron may transition from an occupied orbital (core or valence) to an unoccupied state (conduction or free).Since we treat both core and valence shells equivalently, we simply call these occupied orbitals.Using the setup discussed in §II A, we describe now the transition form factors, f 1→2 (q).The electron is excited from an occupied state |ψ ik to a conduction state |ψ i k , and so f 1→2 (q) −→ f ik→i k , with where we have maintained the PySCF normalization, and Einstein summation over α and β indices is implied.We shall invoke the orthogonality relation The matrix element is then given by where we have used the definition of f [i k ,ik] (q) from Eq. ( 9).We can now plug this into Eq.(B12) to obtain To calculate the total event rate, we must sum over oc-cupied orbitals i and unoccupied orbitals i , and integrate over both k and k .Moreover, we must also consider the spin of the electrons in the occupied bands, giving (B17) Expanding this, and inserting the relevant δ−distributions in q and E e , To simplify this form, we take two steps.First, we define U ≡ k + G − k, and so δ 3 (k + q − k − G) = δ 3 (q − U), which can be further expanded as δ 3 (q − U) = 1 q 2 sin θ q δ(q − U )δ(θ q − θ U )δ(φ q − φ U ) .

(B19)
Here θ U and φ U have the usual definitions as the inclination and azimuthal angles, respectively.From here, we can integrate over Ω q .Second, we differentiate the rate equation with respect to ln E e , and obtain Eqs. ( 8) and (10), dR crystal d ln E e = n χ N cell σe α m 2 e µ 2 χe d ln q E e q η (v min (q, E e )) |F χ (q)| 2 |f crystal (q, E e )| 2 , where (B20)

2 FIG. 1 .
FIG.1.Panels (a) and (b) show the calculated crystal form factor |f crystal (q, Ee)| 2 (see Eq. (14)) for silicon and germanium, respectively.The details of the calculation are described in §III A. The region beneath the black line is kinematically inaccessible for halo DM, as it would require vmin(q, Ee) > vEscape + v Earth (see Appendix B for more details).

1 ∆R Q [ 1
FIG.2.Panels (a) and (b) show the DM-electron scattering rates in a silicon crystal for heavy and light mediators respectively.These values are calculated using a TZP basis set with PBE0 exchange and correlation functionals, 4 × 4 × 4 k−grid, and qmax = 25αme.Panels (c) and (d) show the DM-electron scattering rates in a germanium crystal with interaction mediated by a heavy and a light mediator respectively calculated using a def2-TZVP basis set with PBE0 exchange-correlation functional, 6 × 6 × 6 k−grid, and qmax = 20αme.

FIG. 3 .
FIG.3.Panels (a) and (b) show the DM-electron scattering rates for Si and Ge respectively, calculated with a 4 × 4 × 4 k−grid, with DM parameters noted in the figures.The plots use qmax = 10 αme and PBE exchange-correlation. cc-pVQZ is the most accurate basis set we test for Si, while cc-pVTZ is the same for Ge.Moving further, we use TZP and def2-TZVP for Si and Ge respectively, which mimic the most accurate bases at lower computational costs.

R
FIG.4.Panels (a) and (b) show the electron loss functions integrated over electron recoil energy, Ee, for Si and Ge, respectively, with the dash-dotted line showing the theoretical upper bound from the f -sum rule.For Si, we use here a TZP basis set, at 4 × 4 × 4 k-grid and PBE functional to calculate these results.For Ge, we use here a def2-TZVP basis set, with the same k-grid and PBE functional.Note that a good description of high momentum transfer q 3αme does not require a large Rcut for either element.

m χ = 1
FIG. 5. Panels (a) and (b) show the DM-electron scattering rates in Si and Ge, respectively, for various values of the real space cut-off Rcut, assuming an exposure of 1 kg-year.We use a 4 × 4 × 4 k-grid, set qmax = 6αme, and use the PBE functional.

FIG. 7 .
FIG.6.Panels (a) and (b) show DM-electron scattering rates for Si and Ge, respectively, for various k−grid densities, assuming an exposure of 1 kg-year.We use the PBE exchange-correlation functional and set qmax = 10αme.It is evident that the Si calculation is converged, even at the sparse 4 × 4 × 4 k−grid level, while Ge only converges at higher k−grid densities.

m χ = 1 FIG. 8 .
FIG.8.Panels (a) and (b) shows DM-electron scattering rates for Si with DM-electron interaction mediated by a heavy and a light mediator, respectively, with different qmax cutoffs, assuming an exposure of 1 kg-year.Panels (c) and (d) show the same for a Ge crystal.We use a 4 × 4 × 4 k−grid for Si and a 6 × 6 × 6 k−grid for Ge, and the PBE0 exchange correlation funtional.

F χ = 1 F χ ∝ q − 2
FIG.11.This plot shows the reach for both, obtaining 2.3 events for 1 kg•yr exposure of our target material (solid line), as well as the threshold for a 5−σ discovery by annual modulation with the same exposure (dash-dotted line).Panel (a) shows the upper bounds on cross-section that can be placed for heavy mediators, while panel (b) shows the same for light mediators.We assume no background for either panel, and include 1 ≤ Q ≤ 10 e − for Si (black) and 1 ≤ Q ≤ 15 e − for Ge (blue).

FIG. 12 .
FIG. 12. Panels (a) and (b) show comparison among DM-electron scattering rates in Si assuming a heavy and a light mediator, respectively, calculated using different codes.This work and QEDark both implement the screening described in §II B 1, while EXCEED-DM and DarkELF use their numerically calculated dielectric function.Panels (c) and (d) show the comparative plots for Ge assuming a heavy and a light mediator respectively.

TABLE I .
[52]meters used in the dielectric function calculation in Eq. (11) for Si and Ge from[52].

TABLE II .
Parameters used for DFT calculation of electronic structure of Si and Ge crystals.Egap refers to the scissor corrected bandgap we employ for the materials.