EXCEED-DM: Extended Calculation of Electronic Excitations for Direct Detection of Dark Matter

Direct detection experiments utilizing electronic excitations are spearheading the search for light, sub-GeV, dark matter (DM). It is thus crucial to have accurate predictions for any DM-electron interaction rate in this regime. EXCEED-DM (EXtended Calculation of Electronic Excitations for Direct detection of Dark Matter) computes DM-electron interaction rates with inputs from a variety of ab initio electronic structure calculations. The purpose of this manuscript is two-fold: to familiarize the user with the formalism and inputs of EXCEED-DM, and perform novel calculations to showcase what EXCEED-DM is capable of. We perform four calculations which extend previous results: the scattering rate in the dark photon model, screened with the numerically computed dielectric function, the scattering rate with an interaction potential dependent on the electron velocity, an extended absorption calculation for scalar, pseudoscalar, and vector DM, and the annual modulation of the scattering rate in the dark photon model.


I. INTRODUCTION
Discovering the nature of dark matter (DM) is one of the most important goals in all of physics.
Constituting 26.5% of the universe's energy density [1], DM is roughly five times more abundant than ordinary matter; yet any understanding of its particle content eludes us. One class of ongoing experiments, which search for DM-Standard Model interactions in a laboratory, are known as direct detection experiments. The canonical example of a direct detection experiment looks for nuclear recoils induced by a scattered DM particle. Many of these experiments have been performed, or are proposed, e.g., LUX [2], SuperCDMS [3], ANAIS [4], CRESST [5], SABRE [6], DAMA [7], PandaX [8], DAMIC [9], NEWS-G [10], LZ [2], XENON10/100/1T [11][12][13], KIMS [14], DM-Ice [15], DarkSide [16], and the absence of any signal has led to stringent bounds on a wide variety of DM models. However, all of these experiments share a common kinematics problem when looking for DM with mass lighter than a GeV; when the DM mass drops below the target nucleus mass, the energy deposited in a scattering event rapidly falls below the detection threshold. This is an unfortunate problem because many well-motivated DM models can possess a light, sub-GeV, DM candidate (see Refs. [17][18][19] for recent reviews). This problem has been remedied by an extraordinary influx of new detection ideas for covering this region of DM parameter space (see Refs. [19][20][21][22] for recent reviews). Currently, the preeminent method for searching for sub-GeV DM is via electronic excitations across band gaps in crystal targets . With band gaps of O(eV), DM with masses as light as O(MeV) can be searched for via a scattering event, and O(eV) masses can be searched for via absorption. Experiments such as SuperCDMS [62][63][64], DAMIC [65][66][67][68][69], EDELWEISS [70][71][72], SENSEI [73][74][75], and CDEX [76] are currently in operation searching for DM induced electronic excitations with a combination of Silicon (Si) and Germanium (Ge) targets.
The calculation of DM-electron interaction rates in these targets is more involved than the standard nuclear recoil calculation. This is because electrons in crystal targets cannot be treated as free states. The crystal creates a lattice potential, distorting the electronic wave functions.
Therefore, estimates of DM-electron interaction rates are crucially dependent on how the electronic wave functions, and band structure, are modeled. With experiments already being performed, and the nature of DM still unknown, it is crucial to have the tools to make predictions for the broadest set of DM models [39,78]. The earliest calculations used analytic approximations of the wave 6 not currently have routines for calculating DM-electron scattering rates for general operators, the primitive building blocks, discussed in Sec. II B, are available. Moreover, EXCEED-DM is structured such that calculations of more general DM-electron scattering rates need only be written, and coded, in terms of transition matrix elements, T ≡ F |O|I , in order to be added. This separates the pieces of the calculation dependent on the electronic state approximations, from the observables being computed. Adding new models is straightforward, and in Sec. IV B we discuss computing the DM-electron scattering rate via a scattering potential dependent on the electron velocity which arises from a very simple DM model. EXCEED-DM serves a different purpose than DarkELF, which can compute DM-electron interaction rates, for some DM models, starting from the dielectric function. EXCEED-DM is compatible with DarkELF since it can compute the (q, ω dependent) dielectric function, which can then be used as input to DarkELF. Additionally, EXCEED-DM can account for screening effects, the importance of which was emphasized in Refs. [40,43]. There are two ways to account for screening in EXCEED-DM, via an analytic model of the dielectric function, or using a pre-computed dielectric function. An analytic model for screening was used previously in Ref. [44] and in Sec. IV A we use the dielectric function computed with EXCEED-DM to screen the DM-electron interactions. However, this is not the only pre-computed dielectric function one can use; if available, a file containing the measured dielectric function can also be used to screen the interaction.
The goal of this paper is to familiarize a user with the types of calculations that can be performed with EXCEED-DM, first by introducing a formalism which encompasses all the calculations, and then detailing the inputs to EXCEED-DM and performing new calculations. While these new calculations will serve as examples, they are also new physics results in their own right, and are the most precise determination of DM-electron interaction rates in Si and Ge detectors to date. Our goal is not to derive all the relevant DM-electron interaction formulas from first principles (these are studied in detail in numerous other, referenced, works), but rather illustrate how EXCEED-DM can be used once these formulas are derived. This paper is organized as follows, in Sec. II we discuss the framework under which the EXCEED-DM calculations are performed. While the formalism for DM-electron interaction rates has been developed previously, it is often developed for different purposes, and therefore written in slightly different languages. Since EXCEED-DM can perform many different calculations it is important to have a common theoretical formalism for all the calculations. Moreover, this makes it clear the language in which future extensions, e.g., calculations for other DM models or different electronic state approximations, need to be written in for easy addition to EXCEED-DM. In Sec. II A we discuss the electronic configuration and the approximations regarding the initial and final states that can be used. In Sec. II B we discuss the transition matrix element, the key link between different electronic state approximations and calculations of observables. Calculations in EXCEED-DM are functions of the transition matrix elements, and the transition matrix element is only a function of the electronic configuration. In Sec. II C we overview the calculations that can be performed in v1.0.0 of EXCEED-DM. In Sec. III we discuss the input files needed to perform a calculation with EXCEED-DM. In Sec. IV we perform four new calculations with EXCEED-DM to illustrate what it can do, and how it can be used for publication quality results. Lastly, in Sec. V we discuss future extensions to EXCEED-DM.

II. FORMALISM
EXCEED-DM is built to compute DM-electron interaction rates, e.g., DM scattering and absorption rates. While the formalism of DM-electron interactions has been previously discussed many times [23, 24, 28, 34, 35, 38-40, 43, 46, 50, 51, 57, 59], rarely is there a unifying theme to these derivations; each is usually performed for a specific type of interaction, or approximation concerning the electronic states. However, since the purpose of EXCEED-DM is to compute any of these interaction rates, the building blocks must be abstracted away from a specific calculation. In this section we will introduce these building blocks abstractly, and then discuss the explicit realizations which are currently implemented in EXCEED-DM. This should provide a solid foundation for future extensions relying on the same abstract building blocks, but with different realizations.

A. Electronic State Configuration
The principal building block is the electronic state, represented as |I . These are unit normalized, I|I = 1, eigenstates of the target Hamiltonian,Ĥ 0 , A target is then nothing more than a collection of (orthonormal) electronic states {|I }. Within a target we can define two types of these electronic states: initial, those that are initially filled and final, those that are empty and may be transitioned to. An electronic state configuration, or electronic configuration, is simply a collection of initial and final electronic states, {{|I }, {|F }}.
Clearly, an electronic state configuration depends on the target material, e.g., a Si target will contain a different set of initial and final electronic states than a Ge target. An electronic configuration file is the main target dependent input to EXCEED-DM. 2 Moreover, this is where approximations about the target can enter; a different combination of electronic states constitutes a different approximation about the target. Both initial and final electronic states in the electronic configuration file may be chosen from the available electronic state approximations.
Currently, the available electronic state approximations are specific to crystal targets, i.e., those with a periodic potential. By Bloch's theorem these states can be indexed by a band number, i, and a momentum vector, k, which lies inside the first Brillouin zone (1BZ). We will refer to these states specifically as Bloch states throughout. The position space representation of these states is where s is the component of the spin in theẑ direction, and u i,k,s are the, periodic and dimensionless, Bloch wave functions, where r is a lattice vector. Therefore a Bloch state can be represented by a function which returns u i,k,s . Note that this allows for another layer of abstraction; to implement other approximations for Bloch states one only needs to extend the Bloch wave function implementation with methods for computing u i,k,s . There are three different types of Bloch states which can currently be used in an electronic configuration file, each corresponding to a different basis: bloch_PW_basis, bloch_STO_basis, and bloch_single_PW. Different coefficients must be supplied for each basis.
The bloch_single_PW are simply plane waves, which only requires specifying a momentum vector, p, for which the corresponding k ∈ 1BZ is determined. 3 In the bloch_PW_basis basis, the coefficients, u i,k,s,G must be specified to define the Bloch wave function, 2 See the documentation for more detailed information on explicitly constructing the electronic configuration file. 3 When a bloch_single_PW state is included as a final state, a Fermi factor, F , In the bloch_STO_basis the Bloch wave functions are a linear combination of Slater type orbitals (STOs) [86,87] of each atom in the unit cell. The bands are indexed by the principal quantum numbers, n, , m, as well at the index of the atom in the unit cell, κ, where Ω is the size of the unit cell, y r,κ ≡ x − r − x κ , x κ is the position of the atom in the unit cell. 4 Therefore to specify a bloch_STO_basis state, the central ingredients are the STO coefficients {C j, ,κ , Z j, ,κ , n j, ,κ }. A standard resource for these coefficients are tabulated for 2 ≤ Z ≤ 54 at Ref. [86]. They are also included in the utilities/RHF_wf_data/RHF_wf_data.hdf5 file in the EXCEED-DM repository for convenience.

B. Transition Matrix Elements
Generally, DM-electron interaction rates cause an electron from an initial state to transition to a, higher energy, final state. Fundamentally this means that the most important quantity to compute are matrix elements between initial and final states, or transition matrix elements, where O is some quantum mechanical operator generating the transition between initial state |I and final state |F . These transition matrix elements are the basis of EXCEED-DM. Given an electronic configuration, EXCEED-DM computes the resulting transition matrix elements, in parallel, which are then used in calculations of a specific observable, e.g., DM-electron scattering. The connection between the transition matrix elements and available calculations is discussed in Sec. II C. This places the transition matrix elements directly between the approximations concerning the electronic states in the target, and those which concern the observable. Moreover, any observable which can be written in terms of the transition matrix elements can, in theory, be computed with EXCEED-DM.
Currently, there are six transition form factors which can be computed, where v ≡ p e /m e , p e is the electron momentum operator, σ are the Pauli matrices, and we have left the I, F implicit, for readability. Additionally, each of these operators can be computed in the q → 0 limit where the exponential factor, e iq·x , is dropped. Typically, the q → 0 limit operators are only used in absorption calculations.
For each electronic state approximation, new routines must be added to compute these T . For example, for the Bloch states discussed in Sec. II A, T 1 (q), between an initial state indexed by I = {i, k}, and a final state indexed by F = {f, k }, is given by, where Ω is the volume of the unit cell (UC). See App. C for the derivation and general form of T for the other operators in Eqs. (11) -(16).

C. Calculations
With the electronic state configuration and transition matrix elements defined, any DM-electron interaction rate written in terms of these quantities can, in principle, be added to EXCEED-DM and computed. Currently there are three calculation types supported, binned_scatter_rate, ab-sorption_rate, and dielectric which are discussed in detail below. While each of these calculations computes a specific quantity, other variables in the input file, see Sec. III, can be set to specify different parameters to be used in the calculation. For example, the DM-electron scattering rate computed inside binned_scatter_rate can account for a variety of scattering form factors, which is useful when studying a variety of DM models. The purpose of this section is to give the explicit, general formulas that each calculation is computing, to understand the scope of each calculation.
We will not give derivations for each of these formulae, and encourage the reader to follow the references for more details about individual calculations.

Binned Scatter Rate
We begin by discussing the DM-electron scattering rate calculation. A single formula for DMelectron scattering, per kg-year exposure, R, derived in a variety of other papers [23, 24, 38-40, 43, 50, 51], encompasses a wide number of physical processes. Assuming the scattering rate does not depend on the DM velocity, R is given by, whereσ e is a reference cross section, dependent on the couplings in the UV Lagrangian, m χ is the DM mass, µ χe is the DM-electron reduced mass, ρ χ is the DM density, ρ T is the target density, V is the target volume, I, F index the initial and final electronic states in the target and q is the momentum transferred to the target. We will now go through each of the terms in the integrand in detail. g(q, ω) is the kinematic function, encoding the DM velocity distribution dependence, where f χ (v) is the DM velocity distribution, in the lab frame. A standard choice for velocity distribution is the Standard Halo Model (SHM) [88,89] which has three parameters, the velocity spread, v 0 , the galactic escape velocity, v esc , and the Earth velocity relative to the galactic center, v e . The kinematic function for the SHM model is known analytically, where erf is the error function. 5 Annual and daily modulation effects, e.g., the daily modulation rate in an anisotropic target, considered in Ref. [23], enter the calculation through specific choices of v e (t).
f scr and F med in Eq. (18) parameterize the propagator dependence. f scr incorporates screening effects due to mixing between the propagating dark particle, e.g., a dark photon, and the photon. In the absence of screening effects, f scr = 1, but different DM models can introduce different screening factors. The primary example where screening is important is in a model with a dark photon coupling to the electron vector current, analogous to the Standard Model photon. In this case, the screening factor is dependent on the dielectric function, EXCEED-DM allows the user to set f scr to 1, use an analytic model for ε(q, ω) [90], or use a precomputed dielectric function in an external file. F med is the mediator propagator, not including screening effects, parameterized by β, where β = 0 for a heavy mediator, and β = 2 for a light mediator.
The last term to discuss in Eq. (18) is the scattering form factor, F IF . This is intrinsically DM model dependent, and will take different forms depending on the DM model.
In addition to the total scattering rate given in Eq. (18), EXCEED-DM computes the scattering rate binned in energy and momentum deposition, i.e., the binned scatter rate. This is a matrix, R ij , whose (i, j) component contains the total event rate for events with ω, q in where ∆ω, ∆q are the bin widths, E g is the band gap, and we assume the first index is one. R ij is defined such that R = ij R ij , for any choice of ∆ω, ∆q. 6

Absorption Rate
We now turn to the absorption rate calculation. Similar to the DM-electron scattering rate calculation, the formalism has been studied in detail previously [24,28,34,35,59]. The general absorption rate, R, per kg-year, of DM particle φ is given by, where ρ φ is the DM density, m φ is the DM mass, n is the number of degrees of freedom of φ, Π λ φφ is the self-energy of the φ-like state in medium, projected on to the λth polarization. The difference between φ and the φ-like state is due to mixing effects between φ and A. Π λ 1 λ 2 φ 1 φ 2 is the self-energy between the λ 1 polarization of φ 1 , and the λ 2 polarization of φ 2 (when φ 1 = φ 2 , we simplify the notation by only having one polarization index, λ 1 = λ 2 = λ).
The self-energies in Eq. (27) can be written in terms of the transition form factors in the q → 0 limit. To simplify this we introduceΠ O 1 ,O 2 which is a (polarization independent) selfenergy computed from a one loop diagram [24,34], where δ is the electron lifetime (width). 7 Table II lists the, currently supported, DM model dependent expressions for Πφφ, in Eq. (27), in terms of theΠ's, as derived in Refs. [24,34].

Dielectric Function
The last calculation included in EXCEED-DMv1.0.0 is of the dielectric function, ε(q, ω), which can be computed from the Lindhard formula [92], and written in terms of T 1 as, where G is defined in Eq. (30), e is the electromagnetic charge, and V is the target volume. We include this calculation in EXCEED-DM for two reasons. First, to provide a consistent calculation of DM-electron scattering rates, the screening factor in Eq. (23) should be computed under the same assumptions as the scattering form factor, F IF . This was not done in previous EXCEED-DM calculations for Si and Ge targets, e.g., in Ref. [44] which utilized an analytic dielectric function.
Consistent calculations, performed with v1.0.0 of EXCEED-DM, can be found in Sec. IV A. Second, it has been shown that for a class of simple DM models, DM-electron interaction rates can be written solely in terms of the dielectric function [40,43]. For these models, once a dielectric function has been computed, the package DarkELF can use these pre-computed dielectric functions and compute DM-electron interaction rates. Therefore adding the dielectric function as an output of EXCEED-DM is not only useful as a point of comparison to different calculations of the dielectric function, but also useful as an input to DarkELF.
Technically, the dielectric function in Eq. (31) is not what is computed. While the dielectric function in Eq. (31) is a function of q, due to lattice momentum conservation it is only non-zero is the Bloch momentum of the initial (final) state, and G is a reciprocal lattice vector. Therefore, the requested q point may not be kinematically satisfied with the users choice of initial and final states. For example, if the datasets only contain k I = k F = 0, and q is not a reciprocal lattice vector, the result would be zero. We try to avoid constraining the user to pick an appropriate set of q points. The easiest way around this would be to keep track of all the q that are included in the dataset, but an array of size N I × N F × N G is generally too large to store.
To avoid these issues we compute an averaged dielectric function, where V q bin is the volume of the bin around q. This has the benefits of including all q points in the dataset, and the output is smaller in size since the sum is performed before saving the data. We compare this averaged dielectric function to the analytic model used previously [44] in Sec. IV A.

III. INPUTS TO EXCEED-DM
Before discussing specific applications of EXCEED-DM in Sec. IV, we will discuss the input files needed to run EXCEED-DM. For every calculation there are two necessary files: an input file containing a list of parameters to use for the calculation, and the electronic configuration file which contains all the information regarding the electronic configuration, discussed in more detail in Sec. II A. The discussion here regarding both of these input files is meant to serve as an introduction, so the setup of the calculations in Sec. IV is clear. A complete list of input parameters, and their defaults, can be found in the documentation , as well as the full file specification for the electronic configuration file. A plethora of example input and electronic configuration files can be found in the examples/ folder in the EXCEED-DM repository .

A. Input File
The input file is just a text file where specific variables are specified in different groups. Each group contains a collection of variables relevant to that group. The group is specified inside brackets, '[]', and variables are specified below it. For example, the control group, contains the variables, • calculation -Specifying which calculation is performed.
• out_folder -Specifying the folder where the output file should be saved.
• run_description -A short description of the calculation, appended to the end of 'EXDM_out_' when creating the output filename.
The corresponding input file would contain: Currently, there are five main input groups which control a calculation: • control -Controls what calculations are performed and where the output files are saved.
• elec_config_input -Contains information about the electronic configuration input file.
• material -Specifies some target material parameters.
• astroph_model -Sets the astrophysical parameters of the DM velocity distribution.
For each calculation there is also a numerics_<calculation> group which specifies some numerical parameters specific to that calculation, e.g., energy bin width in a scattering rate calculation.
Having already discussed the control group above, we will now discuss, and give an example of each of these other groups.
The elec_config_input group is the simplest group, currently containing only a single variable, filename, which contains the path to the electronic configuration file, where mX are, again, the DM masses in eV, med_FF is the power of the mediator form factor, β in Eq. (24), (note that more than one may be specified at once), particle_type is the type of DM particle, FIF_id specifies the scattering form factor, and rho_X_GeV_per_cm3 is the DM density in GeV/cm 3 . Not all of these variables are used in every calculation. For example, med_FF and FIF_id are specific to the binned_scatter_rate calculation, and particle_type is specific to the absorption_rate calculation.
Lastly, an example astroph_model group is given by,

Input File
[astroph_model] annual modulation signals to be calculated, see Sec. IV D for a more detailed discussion of annual modulation, and Refs. [23,30,45,56] for a more detailed discussion of daily modulation signals.

B. Electronic Configuration File
The electronic configuration file is stored in a hierarchical data format (HDF5). The pros of this approach are that the data is compressed, and easily (quickly) read in to EXCEED-DM. The cons are that the data is not easy to edit or view without external programs. We recommend the use of HDFView and the python program h5py for viewing and manipulating these files. In this section we will discuss the structure of these files, and leave specific details to the documentation . This should make clear where future extensions may add their components, and aid in understanding the example electronic configuration files, elec_config.hdf5, in each example in the examples/ folder on the Github page . In the directory trees below, folders end with a backslash, and datasets do not. Extra folders/datasets may be included, but they will not be used in the calculations (without modification to EXCEED-DM). This is useful mainly to store explanatory metadata with the file.
At the highest level there is simply an elec_states folder with two similar subfolders, init, and fin, For example, consider of set of Bloch states, I ∈ {i, k}, sampled uniformly in the 1BZ with N k k points, In this case the Jacobian was simple, however it need not be and will depend on the electronic states specified. Consider a set of bloch/single_PW, {p} states, sampled logarithmically in E = p 2 /2m e between E min and E max , and uniformly on the sphere in (θ p , φ p ) space, Uniformly sampling on the sphere involves changing variables to α = φ/2π, β = (1/2)(1 + cos θ), and sampling uniformly between [0, 1] in α, β. Sampling logarithmically in E involves changing variables to x = ln E, and sampling uniformly in x between [ln E min , ln E max ]. Therefore Eq. (37) becomes Ω (2π) 3 dp p 2 dθ sin θ dφ → where I ∈ {α, β, x}. We see that the Jacobian is no longer as simple. Since the Jacobian is inherently tied to the combination of states in the electronic configuration, we leave it to the user to take care of defining these properly.

IV. APPLICATIONS
Practically, the best way to understand what EXCEED-DM can do is to see its applications.
The example calculations in the examples folder offer an introduction to these calculations with small datasets. Additionally, EXCEED-DM has been used previously in a number of papers, see Refs. [23,24,34,44,52,76,93]. However, these papers have not gone in to detail on constructing the inputs for EXCEED-DM. Therefore the purpose of this section is two-fold: showcase the new additions to v1.0.0 of EXCEED-DM by performing new calculations, and explicitly detailing the input files so it is clear how to perform these calculations for other target materials, DM models, etc.
Our focus will be on Si and Ge targets, used in SuperCDMS [62][63][64], DAMIC [65][66][67][68][69], EDEL-WEISS [70][71][72], SENSEI [73][74][75], and CDEX [76]. The electronic configuration files are publicly available here  . The initial states contain a combination of states in the bloch_STO_basis and bloch_PW_basis, and the final states contain a combination of states in the bloch_PW_basis and bloch_single_PW basis. In the language of Ref. [44], these correspond to a set of "core" and "valence" initial states, and "conduction" and "free" final states. More specific details about the electronic states included in these files can be found in the documentation . All of the calculations shown below can be done with these electronic configuration files combined with the input groups shown in each section. The raw output files can be downloaded here . To avoid repetition, the material group variables in the input file for Si and Ge are, respectively. These should be added to the calculation specific input groups discussed in each section below. We will also omit filename input variables which will be computer specific. Specifically, the out_folder, and run_description variables in the control group, and filename in the elec_config_input group.
This section is organized as follows, in Sec. IV A we will perform two new calculations. First, we will compute the average dielectric function,ε(q, ω), with EXCEED-DM and compare it with the analytic dielectric function model used in Ref. [91] to screen the DM-electron interaction. Second, we will use this numerically computed dielectric function to screen the DM-electron scattering rate in the kinetically mixed dark photon model, and compare to the scattering rate with no screening, as well as the scattering rate from Ref. [44] which used the analytic dielectric function to screen the interactions. In Sec. IV B we study DM models beyond the simple dark photon model, and compute the DM-electron scattering rate in an example model which has a scattering potential dependent on the electron velocity. In Sec. IV C we compute scalar, pseudoscalar, and vector DM absorption rates up to DM masses of a keV. This extends the calculation done in Ref. [34] which was previously limited to only the lowest energy, valence → conduction transitions. Lastly, in Sec. IV D we demonstrate EXCEED-DM's ability to compute modulation effects, and compute the annual modulation of the DM-electron scattering rate in the dark photon model.

A. DM-Electron Scattering Rate Screened with Numeric Dielectric Function
Our first application is a new calculation of the DM-electron scattering rate in the canonical, kinetically mixed, dark photon model [94]. As a quick review, this model contains a new, potentially broken, U (1) symmetry whose gauge field is the dark photon, A µ . The DM candidate is a new fermion, χ, which is charged under this new U (1) symmetry. The symmetries of the Lagrangian allow a kinetic mixing between the dark photon and photon, with coupling parameter κ. This interaction can be rotated to the mass basis by shifting the photon field, A µ → A µ + κA µ , generating small "dark electric" charges for the electrons, The DM-electron scattering rate can be shown to be in the form of Eq. (18) (see Ref. [44] and App. B) with, where we have assumed that the targets are isotropic enough to approximateq · ε ·q ≈ ε, and F med takes it standard form for the heavy/light mediator scenarios.
While this rate has been computed for Si and Ge targets many times, [38,39,44,[50][51][52]61], only recently were screening effects included numerically. The most recent calculation using EXCEED-DM uses an analytic model for a dielectric function. While the calculation in Refs. [42,43] includes the screening effect numerically, the all-electron reconstruction effects [36,44], which alter the electronic wave functions near the Fermi surface, were ignored, along with more deeply bound "core" states, and high energy "free" states, in the language of Ref. [44]. The calculation here includes these effects by computing the dielectric function with EXCEED-DM, with the same electronic configuration as used previously (which includes "core" and "free" states, in addition to the all-electron reconstruction corrected valence and conduction states).
The input for the dielectric function calculation is, Since they are only used for scattering rate calculations, the only relevant differences between them are for kinematically allowed q, ω, sectioned off by the black lines. The dashed black line corresponds to the band gap, E g , and the solid black line bounds ω < q/v max χ , where v max χ ∼ 10 −3 . There are O(10%) differences at small q, ω, which will propagate to slightly different DM-electron scattering rates for models that are screened. Both dielectric functions asymptote to one for large q, ω. dielectric function is compared to the analytic model in Fig. 1, defined in Ref. [90] with 0 = 11.3 (14), α = 1.563, ω = 16.6 (15.2) eV, q T F = 4.13 (3.99) keV for a Si (Ge) target. We find that the screening factor is dominated by low q, ω transitions, i.e., the dielectric function contribution from the "core" and "free" states is small, and in the kinematic regions where these transitions are dominant, the dielectric function can be approximated as one.
Moreover, we only find O(10%) changes between the analytic dielectric function model, and the numerically computed dielectric function. This will translate to an O(10%) discrepancy in the DM-electron scattering rate when screened with the analytic model and EXCEED-DM computed dielectric function.
With the dielectric function computed, we now turn to computing the scattering rate with no screening, screening with the analytic dielectric function in Eq. (44), and screening with the numerically computed dielectric function in Eq. (31). Additionally, we will study the dependence on a heavy and light mediator. The screening independent input groups are, for the analytic screening, and numerically computed screening scenarios, respectively (no screening is the default, so no input groups have to be added). Entries in red are specific to Si targets, and those in blue specific to Ge targets. dielectric_filename is the path to the previously computed dielectric function, E_bin_width and q_bin_width are the specific binning parameters used when computing the dielectric function (note that they match the numerics_dielectric input previously discussed). width_id specifies which width parameterization to use (from widths input for dielectric function calculation).
With the inputs discussed, we now discuss the results of the DM-electron scattering rate calcu- solid lines correspond to screening with the numerically computed dielectric function. We see that while the inclusion of screening is important, as noted in Refs. [42,43], the difference between the analytic and numerically computed dielectric function, when used as the screening factor, is small.
lation. In Fig. 2 we compare the binned scattering rate in Si and Ge targets, assuming m χ = 1 GeV, σ e = 10 −40 cm 2 , and different models for the screening factor, f scr . We find that while screening is an important effect, especially for small ω, the difference in using the analytic and numerically computed dielectric functions is small, O(10%), with the difference being smaller in Ge than Si targets. In Fig. 3 we compare the cross section constraints and find similar results for the other DM masses considered. Additionally, we show that the effects of screening for different electron-hole pair number thresholds, Q, defined as, where E g is the target band gap, taken to be 1.11 eV for Si and 0.67 eV for Ge, ω is the energy deposited in a single electron transition, and ε is a material dependent quantity (not the dielectric function) which is taken to be 3.6 eV and 2.9 eV for Si and Ge targets respectively.  (48). We assume a kg-year exposure of Si (left) and Ge (right) targets, assuming the DM model given in Eq. (46), and V µ is light. The red, yellow, and green curves correspond to Q ≥ 1, 2, 3, respectively.
Our next calculation involves a new DM model, different than the kinetically mixed dark photon from Sec. IV A. This DM model was chosen for two reasons: its simplicity, and the fact that the DM-electron scattering potential depends on the electron velocity. This generates a scattering rate dependent on the T v transition matrix element in Eq. (12), and therefore only calculable with a complete spectrum of electronic states with EXCEED-DMv1.0.0. 9 The Lagrangian is given by, 9 While QEDark-EFT [38,39] also supports computing this DM-electron scattering rate, the calculation does not include all electron reconstructed valence and conduction states, and ignores the effects of 'core' and 'free' states, in the language of Ref. [44]. See Sec. I for more details on the discrepancies between the codes.
where V µ is a, light (m V keV), dark photon, and χ is the DM candidate. The scattering form factor, derived in Appendix B, is given by, The reference cross section is,σ and there is no screening, f scr = 1. Fig. 4
In these examples an incoming DM particle deposits its mass energy to the target, driving a, nearly, vertical transition. For some DM models, this absorption rate can be related to the measured, long wavelength dielectric function, ε(0, ω). Specifically, in isotropic targets, e.g., Si and Ge, the absorption rates for axion-like particle and vector DM, e.g., a dark photon, models can be related to ε(0, ω). However, even a simple model of scalar DM cannot be related, see Ref. [34] for more details. Therefore the first principle calculations provided by EXCEED-DM serves three purposes: 1) in DM models which cannot be related to ε(0, ω), provide the only available calculation, 2) for DM models which can be related, the calculation provides a data driven check on the input electronic configuration, and 3) for novel targets, which do not have a rigorously measured ε(0, ω), EXCEED-DM can facilitate a comparison with previously considered targets.
In this section we extend the constraints on DM absorption via electronic excitations to a larger mass range, up to m φ = 1 keV, in the three DM models, scalar, pseudoscalar, and vector DM, considered in Ref. [24,34]. The interaction Lagrangians of each model are given in Table II. The first principles calculations provided by EXCEED-DM for Si and Ge targets in Ref. [34] were limited to m φ < ∼ 60 eV due to a more general implementation not being available in previous versions of EXCEED-DM. Here we use an expanded electronic configuration which includes deeply bound core states, and high energy free states, in the language of Ref. [44]. 10 The inputs for the calculation are, where the particle_type input option specifies which DM model is considered in Table II, and the entries in numerics_absorption_rate follow the same conventions as those for numerics_dielectric in The constraints on the couplings for the three DM models of interest, in Si and Ge targets assuming a kg-year exposure, are given in Fig. 5. Note that the couplings have been conventionally The results for m φ < ∼ 60 eV are in agreement with previous results from Ref. [34], and constraints for where g e is defined in the Lagrangian in Table II and M Pl is the Planck mass.

D. Annual Modulation of DM-Electron Scattering Rate
In our last example, we will study a different aspect of DM-electron scattering in the kinetically mixed dark photon model than discussed in Sec. IV A. Here we will study the annual modulation of the DM-electron scattering rate signal [26,50,61,[101][102][103]. Annual modulation arises due to the Earths motion around the Sun; the Earths velocity relative to the DM wind changes as it revolves around the Sun, leading to changes in the DM velocity distribution. A signal which modulates annually would be a smoking gun signature of DM, since other backgrounds are not expected to have this feature.
the calculation here is an improvement due to more accurate modelling of the electronic structure in Si and Ge, the largest discrepancies are in DM models where these high momentum, high energy differences in the electronic structure models are important. Moreover, the calculation here uses the specific velocity distribution parameters which have been recommended to be used in direct detection calculations [89]. One could also reproduce all the scattering rate results in Sec. IV A with the output data , and compare the dependence on the velocity distribution parameters.
Our focus here will be on changing, v e , while keeping v 0 , v esc set to the recommended values [89].
The input file is the same as in Sec. IV A, with an updated astroph_model group, In Fig. 6 we study f mod , binned in energy deposition, for m χ = 1 GeV for Si and Ge targets. We find that the annual modulation fraction in both targets is generally greater than 10%. While the overall size of the modulation is consistent with previous studies [50,61]. Interestingly, the behavior and Ge (bottom row) targets, for both a light (left column) and heavy (right column) dark photon mediator, modulation fractions can be greater than 10%, and stay relatively large in higher ω transitions.
extends to larger energy depositions, meaning that high threshold experiments will see the same amount of annual modulation as low threshold experiments (albeit at lower overall event rate).
In Fig. 7 we study the annual modulation along a different axis; specifically, f mod versus DM mass for different electron-hole pair thresholds, Q. For low DM masses f mod becomes asymptotically large which is simply due to the shift in kinematically allowed DM masses. The lightest DM mass which can scatter is given by m χ > ∼ 2E g /v 2 max , where v max = v e + v esc . Therefore for different values of v e , the lowest detectable DM mass changes. Higher values of v e correspond to smaller DM masses, and therefore R 0 becomes kinematically inaccessible before R + as the DM mass decreases.
This causes f mod to diverge. At higher masses, f mod generally levels off, except when deeply bound states get involved, e.g., the 3d shell in Ge, whose effects can be seen at m χ ≈ 20 MeV in the bottom right panel of Fig. 7. The dependence on Q is also easily understood from Fig. 6; increasing the threshold picks out the pieces of the binned scattering rate which have larger f mod .

V. SUMMARY AND FUTURE DEVELOPMENT
EXCEED-DM has already been shown to be useful for a wide range of DM-electron interaction calculations: from target comparison studies [23] of DM-electron scattering, proof-of-principle daily modulation signals in anisotropic targets [23], extending the DM-electron scattering rate calculation Ge (bottom row) targets. The different colors correspond to different Q thresholds; the red, yellow, and green curves correspond to Q ≥ 1, 2, 3, respectively.
to include higher energy and momentum contributions [44], performing first principles DM-electron absorption rate calculations [34], and DM-electron interaction rate calculations in spin orbit coupled target [24]. More recently the CDEX collaboration has used it to place constraints with their Ge detector [76], and EXCEED-DM has also been used in studies of non-standard neutrino interactions [93]. Clearly, EXCEED-DM has proved to be a useful tool for studying new physics, even before the v1.0.0 release.  For further installation instructions see the documentation . Our starting point is the NR effective Lagrangian of a DM fermion, χ, coupling to the electron, ψ, where ψ, χ are the two component NR electron and DM fields, respectively, and O is the interaction operator. χ and ψ each have units of eV 3/2 , and O has units of eV −2 . Assume an incoming DM particle with momentum p and spin s scatters off an electronic state |I , creating to an outgoing DM particle with momentum, p = p − q (q is the momentum transferred to the target) and a final electronic state |F . Fermi's Golden Rule tells us that the transition rate due to the Lagrangian in Eq. (B1) is given by, where E I , E F are the energies of the electronic states, ω q = q · v − q 2 2mχ , v is the incoming DM velocity, and m χ is the DM mass. The volume, V , and exponential factor come from the incoming and outgoing DM states, χ|p, s = V −1/2 e ip·x ξ s , and ξ ↑ = where the overline notation indicates (1/2) ss . Assuming the DM velocity distribution is given by f χ (v) and g(q, ω) ≡ 2π d 3 vf χ (v)δ(ω − ω q ), the average transition rate between two electronic states becomes, The transition rate per kg-year, R, is then simply this interaction rate multiplied by the number of DM particles in the target, divided by the target mass, and summed over initial and final states, where ρ χ is the DM density, and ρ T is the target density. Eq. (B5) is most useful if one is interested in finding the rate as a function of the Lagrangian parameters in the NR EFT since it is directly related to the O appearing in Eq. (B1). However, it is useful, and common practice, to write this rate in terms of a reference cross section, which is simply a function of the underlying Lagrangian parameters. In terms of O, the reference cross section is given by, where µ χe is the DM, electron reduced mass. This function is evaluated at reference momentum scales, e.g., q = q 0 = αm e , and ignores screening effects, ε = 1. We can now write Eq. (B5) in terms of the reference cross section. At the risk of abusing notation, we define, and the rate can be written as, From here, the connection to Eq. (18) is a simple rearrangement, and we see that the scattering form factor is simply | F |e iq·xÔ ss |I | 2 with the propagator dependence factored out, With Eq. (B9) we can compute the scattering form factor for the dark photon and 'VA' DM model. We begin with the dark photon model which has the UV Lagrangian, where V µ is the dark photon and it has standard kinetic terms. If the coupling g e is due to kinetic mixing, then g e = eκ, where κ is the kinetic mixing parameter. It can be shown that the NR EFT Lagrangian and O, are, O ss ,σσ = g e g χ ε(q, ω) q 2 + m 2 V δ ss δ σσ , for a material with an isotropic dielectric function, where m V is the dark photon mass, and ω is the energy transferred to the target. We will find the scattering form factor in two limits, for a heavy and light dark photon. We begin by computingÔ ss , O ss = δ ss 1 ε(q, ω) In both cases, f scr = 1/ε, and F med takes its standard form for heavy and light mediators. The scattering form factor is the same in both the light and heavy dark photon models, Moving on to the VA DM model we follow the same steps, first by writing down the UV Lagrangian, L ⊃ g χ V µχ γ µ χ + g e V µē γ µ γ 5 e .
Following Ref. [23,34] after taking the NR limit of the currents and integrating out the dark photon, the NR effective Lagrangian is, where K e = 2k e + q, k e is the electron momentum, K χ = 2k χ − q, k χ is the DM momentum, and σ χ (σ e ) are the Pauli matrices acting in DM (electron) space. Contrary to the kinetically mixed dark photon model, there is no screening. Since the dark photon couples to electrons via the spin operator, σ e , there is no mixing between the dark photon-photon to generate screening terms.
The NR effective Lagrangian in Eq. (B16) can be further simplified since the last two terms are of order O(v χ ), whereas the first term is O(v e ).
Since v e ∼ α ∼ 10 −2 the first term will dominate and we will approximate the interaction with just the first term. Therefore, where K 0 e is some reference momentum, which we take to be αm e . Therefore the scattering form factor is, F IF = 1 (αm e ) 2 4m 2 e |T v·σ | 2 + 2m e T v·σ (q · T σ ) * + 2m e T * v·σ (q · T σ ) + |q · T σ | 2 .
where I, F are now just the band labels.