Direct detection of sub-GeV dark matter with semiconductor targets

Dark matter in the sub-GeV mass range is a theoretically motivated but largely unexplored paradigm. Such light masses are out of reach for conventional nuclear recoil direct detection experiments, but may be detected through the small ionization signals caused by dark matter-electron scattering. Semiconductors are well-studied and are particularly promising target materials because their O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(1 eV) band gaps allow for ionization signals from dark matter particles as light as a few hundred keV. Current direct detection technologies are being adapted for dark matter-electron scattering. In this paper, we provide the theoretical calculations for dark matter-electron scattering rate in semiconductors, overcoming several complications that stem from the many-body nature of the problem. We use density functional theory to numerically calculate the rates for dark matter-electron scattering in silicon and germanium, and estimate the sensitivity for upcoming experiments such as DAMIC and SuperCDMS. We find that the reach for these upcoming experiments has the potential to be orders of magnitude beyond current direct detection constraints and that sub-GeV dark matter has a sizable modulation signal. We also give the first direct detection limits on sub-GeV dark matter from its scattering off electrons in a semiconductor target (silicon) based on published results from DAMIC. We make available publicly our code, QEdark, with which we calculate our results. Our results can be used by experimental collaborations to calculate their own sensitivities based on their specific setup. The searches we propose will probe vast new regions of unexplored dark matter model and parameter space.


The search for dark matter
There has been tremendous progress in the last three decades in the direct detection search for weak-scale dark matter (DM) using underground detectors. The original aim was to probe the scattering through Z-exchange of DM candidates with roughly weak-scale mass against nuclei [1]. Now, experiments searching for these DM-induced nuclear recoils [2][3][4] are sensitive to scattering cross sections many orders of magnitude below the Z-exchange cross section, for candidates in the O(10 GeV-10 TeV) mass range. The motivation behind this incredible experimental achievement has been the theoretically appealing, and dominant, Weakly Interacting Massive Particle (WIMP) paradigm: DM as a weak-scale thermal relic associated with new physics that solves the hierarchy problem. However, the era of this paradigm's preeminence appears to be ending due to both the lack of a DM discovery, which excludes significant regions of WIMP parameter space, and the absence of non-Standard Model (SM) physics at colliders, which has undermined the theoretical motivation behind it. More importantly, several other theoretically motivated candidates exist for resolving this great mystery of particle physics. Motivated particle-DM candidates have been proposed over a vast range of masses, from ultra-light bosonic fields such as a QCD axion [5][6][7], to non-thermal GUT-scale relics [8]. While these have inspired a diverse array of experimental searches, techniques for probing them are far less developed than the WIMP search program. One well-motivated candidate that has received increased attention recently and is the focus of this paper is light dark matter (LDM), with DM masses in the MeV to GeV range. LDM is often motivated by production mechanisms that go beyond the standard freeze-out and may be found in several frameworks in which the sub-GeV mass scale arises naturally. In addition, the origin for the DM relic density can be naturally addressed by several mechanisms that suggest that LDM interacts with SM particles via, for example, an exchange of a light "dark photon", an axion, or through an electromagnetic dipole moment. There is a large range of parameter space of such models that evades both laboratory and astrophysical bounds .
Investigating LDM is an important and natural direction to pursue in the DM search effort. An essential part of this pursuit is extending direct detection searches to this low mass range. Several possible ways to do this were described in [9]. Fortunately, much of JHEP05(2016)046 the impressive technology being developed for the Weak-scale direct detection program can be readily adapted to search for LDM. An example of this was described in [31], obtaining the first direct detection limits on DM with masses as low as a few MeV using published XENON10 data. In this work, we study in detail the even more promising possibility of semiconductor-based LDM searches, significantly expanding the preliminary work done in [9]. Other, complementary techniques to search for LDM have been discussed in .

Direct detection of sub-GeV dark matter
Current direct detection experiments are limited to probing DM masses above a few GeV due to the high energy thresholds required for detecting nuclear recoils. The challenge in probing lower DM masses is twofold: for lower masses, not only is the total kinetic energy of the DM particle decreased, but so is the fraction of energy that is transferred to the nucleus. As a result, the energy of the nuclear recoils is much lower and one must drastically reduce the threshold energies to detect it. This is an experimentally challenging task, although it may be possible to probe masses down to a few hundred MeV, see [62,63]. Instead, as discussed in [9], scattering channels other than elastic nuclear recoil are likely to be far more fruitful.
A very promising avenue is to search for the small ionization signals caused directly by DM-electron scattering. The lightness of the electron and the inelastic nature of the DM-electron scattering process allow DM particles to transfer a large fraction of their kinetic energy to the electron when they scatter, enabling DM as light as ∼1 MeV to cause an ionization signal. Furthermore, detecting small ionization signals is already a well-developed part of direct detection technology. In fact, the XENON10 experiment was already sensitive to the ionization of a single electron [64], and results of a short singleelectron-sensitive run [65] were used in [31] to place direct detection bounds on DM with masses as low as a few MeV. This serves as a proof-of-principle, motivating dedicated LDM searches in other dual-phase noble liquid experiments such as XENON100 and LUX. However, semiconductor targets have the potential to probe even smaller cross sections. In semiconductors such as silicon or germanium, the band gap (the threshold to "ionize" an electron by exciting it from a valence band to a conduction band) is ∼1 eV -a factor of 10 to 20 times lower than the ionization threshold in liquid xenon. The consequences of this lower energy threshold are significant. Not only could this allow sensitivity to DM down to masses below an MeV, but it would also mean a substantial increase in event rate for all DM masses [9,25]. The reason for this is that, given the characteristic velocities of DM particles and electrons, ∼1 eV recoil energies are typical, while recoil energies of ∼10 eV require velocities that are only found on the tails of the DM and electron velocity distributions. Moreover, although the background that causes an ionization signal at such low energies is still poorly understood, it is reasonable to expect that background event rates in semiconductors may be significantly lower than in xenon-based detectors [66] (especially since they may be operated cryogenically). 1 There is currently an active program in, for JHEP05(2016)046 example, both the SuperCDMS and DAMIC collaborations to develop germanium-and silicon-based detectors that are sensitive to single electron-hole pairs [66,67], enabling a leap forward in LDM detection.
This developing experimental program presents a new theoretical challenge: the calculation of the expected signal rate. Unlike for elastic nuclear recoils, this calculation is highly non-trivial. In this paper, we tackle this calculation head on and present detailed new results for germanium and silicon targets.

The challenge of calculating event rates
Several factors complicate the calculation of DM-electron scattering rates. Bound electrons in dense media have a) typical speeds of order α ≈ 1/137 or greater, much faster than DM particles (with v ∼ 10 −3 ), b) indefinite momentum, with even very large momenta having non-zero probability, and c) a complicated structure of energy-levels. This greatly modifies the scattering kinematics and breaks the simple link between momentum transfer and energy deposition. As we discuss in more detail below, event rates can be highly sensitive to the energy-level structure and the tails of the electrons' momentum distributions. In addition, the quantum nature of both the initial and final electron states is important, and they cannot be correctly treated classically. As a result, approximate calculations which do not fully account for these details may not give accurate results. This becomes even more important for the large energy depositions, well above O(eV), since these rely on the tails of the electron's momentum distribution. Once correctly calculated, the effect of all these complications can be completely encoded in an atomic form factor [9]. This function is different for each specific target material, but is independent of the DM model. Once it is known, event rates can be calculated relatively simply.
When the target is an isolated noble gas atom, the combination of spherical symmetry and previously-compiled bound-state electron wavefunctions makes calculation of the ionization form factor relatively straightforward. Refs. [9,31] used this as an approximation for the form factor of a liquid xenon target. However, calculating the form factor for a crystal target (such as a semiconductor) is far more challenging. A periodic crystal lattice is a complex multi-body system, with outer-shell (valence) electrons delocalized and occupying a complicated energy band-structure. Accurate wavefunctions of the valence electrons cannot be found analytically, but must be computed numerically with an expansion in a discrete set of plane waves. 2 Taking this approach, a first calculation was done in [9], assuming a single-electron threshold in a germanium target. A second approach was taken in [25], which succeeded in simplifying the calculation until it was analytically tractable. However, the approximations required for this were so extensive that the result might be considered only as an order-of-magnitude estimate. A third, semi-analytic approach was taken in [68] (see "Note added"), where numerical bound-state wavefunctions for free germanium and silicon atoms are used and the outgoing electrons are described by plane waves. The latter approach gives answers much closer to our full numerical calculation, but important differences remain. JHEP05(2016)046

Overview of the paper
In this paper, we present the results of a detailed numerical computation of DM-electron scattering rates in germanium and silicon targets as a function of the electron recoil energy. This significantly expands on the previous calculation in [9]. Higher recoil energies for the scattered electron allow a larger number of additional electron-hole pairs to be promoted via secondary scattering. Using a semi-empirical understanding of these secondary scattering processes, we convert our calculated differential event rate to an estimated event rate as a function of the number of observed electron-hole pairs. These results will allow several experimental collaborations, such as DAMIC and SuperCDMS, to calculate their projected sensitivity to the DM-electron scattering cross-section, given their specific experimental setups and thresholds. It will also allow them to derive limits on this cross section in the absence of a signal, or the preferred cross section value should there be a signal, in forthcoming data. Achieving low ionization thresholds could allow these experiments to probe large regions of LDM parameter space in the near future, as illustrated in figure 1.
In section 2, we briefly discuss the direct detection prospects for a few popular LDM models. We will see that the upcoming generation of experiments with semiconductor targets play an essential role in testing these models. In section 3, we outline how to calculate the rate for DM to scatter off bound electrons. We provide an intuitive understanding of the scattering kinematics. Our discussion is general and applicable to both electrons bound to (free) atoms as well as electrons in semiconductor targets. The details of this calculation as well as comprehensive formulas are contained in appendix A, significantly expanding on the information contained in [9,31]. We then focus on semiconductor targets, and describe the numerical computation of the scattering rates in section 4. We describe our code QEdark, which is an additional module to the publicly available code Quantum ESPRESSO [69]. The latter calculates the band structure and all electron wavefunctions using density functional theory (DFT) and pseudopotentials, two established condensed matter computational tools, to calculate the Bloch wavefunction coefficients for the initial and final state electrons. In QEdark, we use this information to calculate the crystal form factor for DM-electron scattering as well as the scattering rates. QEdark and the crystal form factors will be publicly available at this link. In section 5, we discuss the conversion from the energy of the primary scattered electron to the size of the final ionization signal. We present a conversion formula and discuss the uncertainty associated with it. In section 6, we present the results of our computation, showing the cross-section sensitivity as a function of detector threshold, as well as the potential discovery reach using annual modulation. We also provide detailed sensitivity estimates for two representative, near-term experiments that may soon reach the required sensitivity to detect LDM, namely DAMIC and SuperCDMS. We conclude in section 7. The appendices contain additional technical details: appendix A provides a detailed derivation of the formulae for the scattering rate and crystal form-factor, appendix B describes our choice of local DM velocity distribution, appendix C discusses the convergence of our numerical results, appendix D studies the effects of inner-shell electrons on the overall scattering rate, appendix E presents details of the systematic study of secondary interactions, and appendix F gives a brief review of DFT and pseudopotentials. We note that our main results are contained in figures 1, 2, 6, and 9 and described in section 2 and section 6.

Models of Light Dark Matter
Theories of LDM have been receiving increased attention in recent years. Here we illustrate with just a few benchmark LDM models how the upcoming generation of experiments with semiconductor targets, including SuperCDMS and DAMIC, play an essential role in the search for LDM. Classes of models that are probed by LDM direct detection include DM that scatters through a dark-photon mediator or through a dipole moment interaction. We

JHEP05(2016)046
focus on DM coupled to a dark photon, leaving a discussion of dipole moment interactions [25,70], the SIMP [29,30], and other models that can be constrained by electron recoils to an upcoming publication [71].
For illustration, we consider models of LDM based on the vector-portal, in which the dark sector (and the DM particle, χ) communicates with the SM through a U(1) D gauge boson A . The A is kinetically mixed with the SM hypercharge U(1) Y via the interaction causing it to couple dominantly to electrically charged particles at low energies. Here is the kinetic mixing parameter, θ W is the Weinberg mixing angle, and F µν DM particles can scatter off electrons in direct-detection experiments through A exchange. In the notation of section 3.2 below, the DM-electron reference cross section is given by where µ χe is the DM-electron reduced mass and α D ≡ g 2 D /4π (with g D the U(1) D gauge coupling). We note that this expression is the same for DM that is a complex scalar or a fermion. The corresponding DM form factor is where q is the momentum transfer between the DM and electron.
In figure 2, we illustrate the parameter spaces of both the m A αm e and m A αm e regimes, taking the fermionic and complex-scalar cases separately for the former. We study three cases, which highlight different possible production mechanisms, and show the interplay between different experimental probes.
(i) Freeze-out via the vector portal: complex scalar LDM.
We consider the phenomenologically interesting and predictive region m A > 2m χ , corresponding to F DM (q) = 1. Annihilation to SM particles occurs via an off-shell A (χχ * → A * → SM). This process is p-wave suppressed, allowing the DM abundance to be set by thermal freeze-out while evading constraints from the cosmic microwave background (CMB), e.g. [72,73], and from gamma-rays in the Milky-Way halo [74]. We show the parameter space for this scenario in figure 2 (top left), taking m A = 3m χ for concreteness. The thick blue curve shows the cross-section for which the correct relic abundance is obtained from freeze-out [73] (this is largely insensitive to the specific choice of m A ). Above this line, an asymmetric DM component may complete the DM abundance. Below it, the abundance is naively too large, but this region may be viable with alternate hidden-sector freeze-out channels. We also show various constraints on this model. The black curve labelled "XENON10" shows JHEP05(2016)046 Figure 2. Prospects for benchmark models: selected 95% C.L. exclusion reach for the DAMIC (green curves) and SuperCDMS-silicon (dark red curves) experiments, compared with other constraints for the benchmark models discussed in section 2. White regions are unconstrained, while thick blue curves illustrate possible predictive mechanisms for generating the DM abundance. Top: DM interacting via a massive dark photon (F DM (q) = 1), for complex-scalar DM with freeze-out abundance (left), and Diracfermion DM with asymmetric abundance (right). Bottom: DM interacting via an ultralight dark photon (F DM (q) = (αm e /q) 2 ), with an abundance generated by freeze-in. The DAMIC and SuperCDMS projections assume 100 g-year and 10 kg-years background-free exposures, with 2and 1-electron thresholds, respectively, in a silicon target. See text for details.
the electron-recoil DM constraint set with XENON10 data [31]. The black curve labelled "Current NR Constraints" shows constraints from conventional nuclear-recoil searches from [3,75,76]. Some measurements only constrain as a function of m A . Among these, we only show the strongest constraints, which are a BaBar search for e + e − → γ + invisible [49,51,52] as well as electroweak precision tests (EWPT) [80,81]; however, to guide the eye, we also show the "favored" 2σ-region for which the A can explain the discrepancy between the measurement and SM prediction for the muon anomalous magnetic moment, a µ [82]. We translate these into the σ e versus m χ plane by using the constraint on α D from either perturbativity [83] JHEP05(2016)046 or χ self-interactions [22]. For these we require that α D is less than 1.0 and small enough so that σ self−int /m χ 1 cm 2 /g for clusters [84]. A second set of constraints bound some combination of , α D , and m A : the electron beam-dump E137 [57,85] and the proton beam-dump LSND [60,86,87]. We again use the constraint on α D from self-interactions and perturbativity to translate these into the σ e versus m χ plane. We also show a rough bound on N eff , see [53,54,73]; the presence of additional relativistic degrees of freedom could allow this bound to be evaded. For a complementary representation of this parameter space see [59].
(ii) Freeze-out via the vector portal: Dirac fermion LDM.
In figure 2 (top right), we consider the same scenario as in (i) but take χ to be a Dirac fermion. This also corresponds to F DM (q) = 1. The main difference between this scenario and (i) is that the annihilation cross section is now s-wave, so that constraints from the CMB preclude the abundance being set by freeze out. Instead, we assume the abundance to be asymmetric [88][89][90], and require the symmetric component to be small enough after freeze-out to avoid the CMB bounds [20]. This provides a lower bound on the annihilation cross-section and thus on σ e , shown with a black solid line. As before, this lower bound is model-dependent and can be evaded with additional annihilation channels. The other constraints are similar.
(iii) Freeze-in via the vector portal.
In figure 2 (bottom), we consider an ultra-light A mediator (m A αm e ), corresponding to F DM (q) = (αm e /q) 2 . Here the couplings are so small that the DM would never have thermalized with the SM sector. The χ abundance can receive an irreducible "freeze-in" [91] contribution from 2 → 2 annihilation of SM particles to χχ as well as Z-boson decays to χχ, computed in [9] (see also [24]). The parameters required for the abundance again uniquely constrain σ e versus m χ , as shown by the thick blue curve. In addition to the XENON10 electron-recoil constraint [31], we also show the bounds from conventional nuclear-recoil searches. The nuclear recoil cross-section, σ NR , can be related to the electron recoil cross-section by where the target nucleus has mass m N and atomic number Z, E NR is the nuclear recoil energy, and v min,NR = √ 2m N E NR /(2µ χN ), v is relative velocity of the DM, and η is the inverse mean speed defined in appendix B. Since this recoil spectrum is peaked towards low energies more than for a contact interaction, determining accurate DM constraints requires a careful analysis of the experimental data. We place approximate bounds from "CDMSLite" [92] and LUX [3] results, taking the former to have 6.2 kg-days germanium exposure, a 0.84 keV threshold, 100% signal efficiency and 10 observed events, and taking the latter to have 10 tonne-days xenon exposure, a 5 keV threshold, 50% signal efficiency and 0 observed events. Due to the smallness of the couplings, the other constraints seen in the previous scenarios are absent in this one.

JHEP05(2016)046
In each of the figures in figure 2 we show the prospects for DAMIC (100 g-years, silicon target, 2-electron threshold) and SuperCDMS (10 kg-years, silicon, 1-electron threshold), discussed in section 6.5. We note that a magnetic-dipole-moment interaction would also give F DM (q) = 1, while an electric-dipole-moment interaction would give F DM (q) = αm e /q. We see that these models above all have concrete predictions that the upcoming generation of direct detection experiments can test.

Direct detection of dark matter by electron scattering in semiconductors
In this section, we review the theory of DM scattering with bound electrons. We begin in section 3.1 by considering the simple kinematics of LDM scattering with both nucleons and electrons. This makes clear the motivation for using electron recoils to probe LDM. The discussion also shows that the DM-electron scattering rate is expected to be sensitive to the details of electron binding in the target, especially for higher energy/ionization thresholds.
A consequence of this is that to calculate accurate scattering rates, detailed modeling of the electronic structure of the target material is required, involving knowledge of the wavefunctions of all accessible occupied and unoccupied electron levels.
In section 3.2, we summarize how this scattering-rate calculation is formulated, with a focus on the case of semiconductor targets. The key results are eqs. (3.13) and (3.17). The former gives the differential scattering rate in terms of the DM model, the DM velocity profile, and crystal form factor. The latter gives the crystal form-factor, which encodes all the relevant electron binding effects for a given target material. This reviews and extends the discussion from ref. [9]. In appendix A, we provide a full derivation of all the results given here. For the interested reader, in appendix A.3 we present a derivation of the ionization rate in free atomic targets, as is relevant for xenon targets and which was used in refs. [9,31].

Kinematics of dark matter-electron scattering
Conventional DM direct detection experiments assume that the DM particle scatters elastically off a target nucleus. This recoiling nucleus then collides with the surrounding matter within the detector, giving off energy in the form of heat, phonons, ionized electrons, scintillation photons, etc, depending on the detector material. However, if the DM particle is light, the momentum transfer, q, between the DM and the target nucleus is small and may not provide enough energy for the recoil of the nucleus to be detected. We can see this through the following argument. The energy of the recoiling nucleus in nuclear scattering is For the scaling in the last step of this equation, we have taken the typical DM speed to be 300 km/s ≈ 10 −3 c, and assumed m χ m N . For m χ = 30 GeV, we find E NR ∼ 2 keV. However, if we consider lighter DM masses, such as m χ = 100 MeV, the recoil energy drops to E NR ∼ eV, which is well below the detection thresholds of current direct detection JHEP05(2016)046 Figure 3. The scattering of a DM particle with a bound electron. The DM transfers momentum q to the target, exciting it from the ground state X to an excited state X * , which can be either a higher-energy bound state or an ionized state.
experiments (e.g. ∼ 840 eV NR for CDMSlite [92] and ∼ 4 keV NR for LUX [3]). Note that the energy of the recoiling nucleus is also not efficiently transferred to electrons, and so is not nearly large enough to ionize or excite even a single electron; it is also well below current phonon detection thresholds. As a result, DM masses below a few hundred MeV escape detection no matter how large their cross section. Now consider a DM particle colliding directly with a bound electron, exciting it to a higher energy level or an unbound state, as illustrated in figure 3. The kinematics are very different from those of a nuclear recoil. Firstly, being in a bound state, the electron does not have definite momentum -in fact it may have arbitrarily high momentum (albeit with low probability). This breaks the direct relation between recoil energy and momentum transfer given in eq. (3.1). The energy transferred to the electron, ∆E e , can still be related to the momentum lost by the DM, q, via energy conservation: Here the ∆E N term accounts for the fact that the whole atom also recoils. In practice this term is small, which also allows us to replace µ χN with m χ . We thus define as the energy transferred to the electron. 3 Since an arbitrary-size momentum transfer is now possible, the largest allowed energy transfer is found by maximizing ∆E e with respect to q, giving This shows that all the kinetic energy in the DM-atom collision is (in principle) available to excite the electron. For a semiconductor with an O(eV) bandgap, ionization can be caused by DM as light as O(MeV).
What is the likelihood of actually obtaining a large enough q to excite the electron? This brings us to the second major difference compared to DM-nuclear scattering: the JHEP05(2016)046 electron is both the lightest and fastest particle in the problem. The typical velocity of a bound electron is v e ∼ Z eff α, where Z eff is 1 for outer shell electrons and larger for inner shells. This is much greater than the typical DM velocity of v ∼ 10 −3 . The typical size of the momentum transfer is therefore set by the electron's momentum, where v rel is the relative velocity between the DM and electron. Returning to eq. (3.2), the first term on the right dominates as long as m χ is well above the bound in eq. (3.4). This gives a formula for the minimum momentum transfer required to obtain an energy ∆E e : This scaling suggests that the typical available momentum is enough to cause a transition of just a few eV, such as for an electron being excited just across the germanium or silicon bandgap. Exciting a more energetic transition will require a momentum out on the tail of the electron's momentum-space wavefunction (or probing the tail of the DM velocity distribution), and its probability will be correspondingly suppressed (as can be seen clearly in figure 5 below, which we will discuss in section 6.1). Ionization of a xenon atom, requiring ∼10 eV energy, falls into the second category, as do most possible transitions to the conduction band in germanium or silicon. From this argument we expect the rate of DM-electron scattering to be sensitive to the precise forms of the electron energy levels and wavefunctions in the target. The computation we present below is designed to address this sensitivity by modeling in detail the electronic structure in germanium and silicon crystals. A corollary of this argument is that, given the v-dependence in eq. (3.6), the rate should also be sensitive to the DM velocity profile. As this varies over the year, we expect a significant annual modulation in the signal size, a potentially crucial test of the DM origin of a signal. We discuss the expected annual modulation in section 6.4.

General formulation for dark matter-induced electron transitions
If a DM particle scatters with an electron in a stationary bound state such as in an atom, it can excite the electron from an initial energy level 1 to an excited energy level 2 by transferring energy ∆E 1→2 and momentum q. The cross section for this process takes quite a different form to the free elastic scattering cross section.
If M free ( q ) is the matrix element for free elastic scattering of a DM particle and an electron, then we parametrize the underlying DM-electron coupling using the following definitions [9]:

JHEP05(2016)046
where |M| 2 is the absolute square of M, averaged over initial and summed over final particle spins. The DM form factor, F DM (q), gives the momentum-transfer dependence of the interaction -for example, F DM (q) = 1 results from a point-like interaction induced by the exchange of a heavy vector mediator or magnetic dipole moment coupling, F DM (q) = (αm e /q) for an electric dipole moment coupling, and F DM (q) = (αm e /q) 2 for exchange of a massless or ultra-light vector mediator (see section 2). σ e parameterizes the strength of the interaction, and in the case of F DM (q) = 1 is equal to the cross section for free elastic scattering. All sensitivity estimates or constraints on LDM will be given for σ e , which plays the analogous role to σ χN , the DM-nucleon scattering cross section, in (WIMP) DM scattering with nuclei. With these definitions, the cross section for a DM particle to excite an electron from level 1 to level 2 can be written as (see appendix A.1) where f 1→2 ( q ) is the atomic form factor for the excitation. It is given by where ψ 1 and ψ 2 are the normalized wavefunctions of the initial and final electron levels.
We now apply this general result to the special case of electrons in a periodic crystal lattice, such as a semiconductor.

Excitation rate in a semiconductor crystal
The periodic lattice of a semiconductor crystal has a continuum of electron energy levels, forming a complicated band structure (see figure 4). A small energy gap separates the occupied valence bands from the unoccupied conduction bands; exciting electrons across this bandgap creates mobile electron-hole pairs, which can be manipulated and detected. In order to perform practical calculations for this system, the true multi-body electron wavefunction must be replaced with a product of single-particle wavefunctions (this is a well-understood procedure, which we discuss further in section 4). Once found, these singleparticle wavefunctions can be used in eqs. (3.9) and (3.10), giving the cross-section to excite an electron between specific energy levels. To find the total rate, these cross sections are integrated over initial and final electron levels, and over the DM velocity distribution.
DM halo dependence. Neither the electron band structure, nor the electron wavefunctions, nor the DM velocity distribution are spherically symmetric. As noted in [9], the excitation rate will therefore depend on the orientation of the crystal with respect to the galaxy, an effect which may be extremely useful in verifying the DM origin of a signal. Here, however, we sidestep this complication by approximating the DM velocity distribution as being a spherically symmetric function g χ (v). All the relevant information about the DM velocity profile can then be encoded in the function η(v min ) (see appendix A.2), defined as where Θ is the Heaviside step function. When calculating rates, we assume a Maxwell-Boltzmann distribution with a sharp cutoff (we describe this in more detail, and give analytic formulas for η(v min ), in appendix B). The requirement of energy conservation is captured by v min (q, E e ), the minimum speed a DM particle requires in order for the electron to gain an energy E e with momentum transfer q (note that E e was also denoted as ∆E e in section 3.1). This is given by

JHEP05(2016)046
Differential rate. As we show in appendix A.4, the differential electron scattering rate in a semiconductor target (with the approximation of a spherically symmetric DM velocity distribution) can be written as where ρ χ 0.4 GeV/cm 3 is the local DM density, E e is the total energy deposited, and N cell = M target /M cell is the number of unit cells in the crystal target. (M cell = 2 × m Ge = 145.28 amu = 135.33 GeV for germanium, and M cell = 2 × m Si = 56.18 amu = 52.33 GeV for silicon.) We have written this in such a way that the first line gives a rough estimate of the rate, about 29 (11) events/kg/day for silicon (germanium) for ρ χ = 0.4 GeV/cm 3 , m χ = 100 MeV, and σ e 3.6 × 10 −37 cm 2 (the current limit from XENON10 [31]), while every factor in the second line is a roughly O(1) number for the preferred values of q and E e .

JHEP05(2016)046
All the necessary details of the target's electronic structure are contained in the dimensionless crystal form factor, f crystal (q, E e ), which is a property purely of the target material and is independent of any DM physics. The computation of this form factor is one of the main results of this paper.
Crystal form factor. In the periodic lattice of a semiconductor crystal, each electron energy level is labelled by a continuous wavevector k in the first Brillouin Zone (BZ), and by a discrete band index i. The wavefunctions of these states can be written in Bloch form, where the G's are the reciprocal lattice vectors. Here V is the volume of the crystal, and the wavefunctions are taken to be unit-normalized, so that Using this form for the wavefunctions, we can define the form factor for excitation from valence level {i k} to conduction level {i k }, The crystal form factor required in eq. (3.13) is then given by (See appendix A.4 for the derivation.) The band index i is summed over the filled energy bands, while i is summed over unfilled bands, and the momentum integrals are over the 1st BZ. E i k is the energy of level {i k}, and V cell is the volume of the unit cell. The numerator in the first factor has units of energy, with value 2π 2 (αm 2 e V cell ) −1 = 1.8 eV for germanium and 2.0 eV for silicon. The crystal form factor can be computed numerically using established solid-state computational techniques. Once it is known, it can be used to find event rates for any DM model and halo profile, using eq. (3.13), along with eqs. (3.7), (3.8), (3.11), and (3.12). We now turn to our own numerical evaluation of the crystal form factor.

Numerical computation of the form factor
Our aim is to compute the crystal form factor, given by eq. (3.17), for silicon and germanium targets with low energy thresholds ( < ∼ 30 eV). Once these are found, it is possible to calculate scattering rates for any DM model. Calculating the form factor requires knowledge of the electron wavefunction coefficients u i ( k + G) for all energetically accessible electron levels. To calculate these coefficients, we utilize the "plane wave self-consistent field" (PWscf)

JHEP05(2016)046
code within the Quantum ESPRESSO [69] package, based on the formalism of DFT. We then input these into our own postprocessing code, QEdark, to calculate the form factors. In this section, we summarize the key conceptual and numerical details of our computation. We provide a review of DFT in appendix F, detail the approximations used in the computation of the wavefunctions, and lay out the numerical methods. In appendix C, we discuss the convergence of our computation.

Computational framework
It is impossible in practice to obtain the exact many-electron wavefunctions that describe interacting electrons in a many-body system such as a crystal. Instead, several methods exist to obtain excellent numerical approximations to these wavefunctions. We use DFT, which reformulates the interacting quantum many-body problem in terms of functionals of the particle density n( r). For the case of electrons, the Hohenberg-Kohn theorems [93] imply that all properties of the interacting system are determined once the ground-state electron density is known. In order to obtain the ground-state density, we use the Kohn-Sham method [94] to map the system of interacting electrons into a system of independent electrons under the presence of an auxiliary potential that produces the same ground-state density. After this mapping, one has to solve the much simpler non-interacting electron system in order to obtain the ground-state energy and electron density.
The mapping from an interacting to a non-interacting many-body system comes at the expense of having to use an approximate auxiliary potential. Typically this potential is split into the mean-field Hartree potential and an exchange-correlation potential. The latter captures the quantum mechanical effect of having identical electrons and also attempts to capture the correlation energy among the interacting electrons. The exchange-correlation potential is not known exactly and needs to be approximated. We use the Perdew-Burke-Ernzerhof (PBE) functional [95], which belongs to the class of the Generalized Gradient Approximations (GGA). We discuss this further in appendix F.
Both silicon and germanium have a diamond lattice structure that contains two atoms in the unit cell. There are two s-shell and two p-shell valence electrons per atom (3s and 3p (2s and 2p) for germanium (silicon)), which makes a total of 8 electrons per cell. This translates to 4 valence bands, since each band is doubly degenerate in electron spin. In silicon, the core electrons have binding energies of ∼100 eV and above, and so are irrelevant for the energies we consider here. One must take more care with germanium, since the 3d electrons have binding energies of ∼30 eV, and so can be relevant for the higher energy thresholds we consider here. 4 In the computation, energetically-inaccessible core electrons can be replaced with a pseudopotential, which increases the computational efficiency by reducing the number of initial states required, and by reducing the resolution needed to describe the wavefunctions. We use ultrasoft pseudopotentials [96] in place of all but the outer two s-shell and two p-shell valence electrons. For germanium, we also use a pseudopotential that allows us to treat the 3d electrons as valence states. As a result, the computational cost for germanium is slightly higher than that of silicon. We use an JHEP05(2016)046 empirical "scissor correction" approach [97,98] to set the band gap to 1.11 eV for silicon and 0.67 eV for germanium [99]. 5

Discretization procedure and cutoff choices
In order to obtain the crystal form factor with a finite computation, several modifications must be made to eq. (3.17): • Binning in q and E e . The form factor must be evaluated for finite grid of q-and E e -values. We do this by averaging over bins of equal width in q and E e : Here q n is the central value of the n th q bin, and E m is the central value of the m th energy bin, and ∆q and ∆E are the widths of the bins. We use 500 E e -bins with ∆E = 0.1 eV and 900 q-bins with ∆q = 0.02 αm e .
• Discretization in k. The continuum of k-values in each energy band must be replaced with a discrete mesh of representative k-points. The k-integrals in eq. (3.17) are then replaced with finite sums: Here V BZ is the volume of the Brillouin Zone, V cell is the volume of the crystal's unit cell, and w k are the weightings of the k-points, with w k = 2 (following the convention of Quantum ESPRESSO). We use a uniform 243 k-point mesh.
• Cutoff in G, G . The wavefunctions are expanded in a finite size plane-wave basis whose reciprocal lattice vectors satisfy the "kinetic energy" cutoff (really a cutoff in the space of G-vectors) Note that since q = | k − k + G |, and since | G max | | k| and | k |, the momentum transfer q essentially has a cutoff of √ 2m e E cut . We choose a value of E cut = 70 Ry, which allows us to sample a large enough q space to obtain O(1%) accuracy for our rate calculations.
• Energy bands. As discussed above, we consider initial electron states in the 4 valence bands for silicon and the 4 valence bands + 10 outer core bands (corresponding to the 3d-shell electrons) for germanium. We include final-state energy bands up to the 52 nd conduction band in both germanium and silicon. The lowest conduction states not included are about 57 eV above the band gap, while the highest energy core states not included are more than 60 eV below the band gap. Our choice of bands therefore fully covers any energy transition below ∼57 eV.

JHEP05(2016)046
We can now write the form factor in the form that is implemented in our numerical code: Note that the first line here represents summing over bands, k-points, and reciprocal lattice vectors, and calculating the contribution to the form factor from each. The sums are all over finite ranges as discussed above. The second line represents adding each contribution to the appropriate {q, E e } bin. We present the results of our computation, including prospects for upcoming experiments, in section 6. In appendix C we discuss convergence with respect to the choice of k-point mesh and E cut .

Conversion from energy to ionization
The calculation described in the previous two sections gives the DM-electron scattering rate in a semiconductor crystals as a function of the total energy deposited by the dark matter, E e . However, experiments will not directly measure the deposited energy itself, but rather the ionization signal Q -i.e., the number of electron-hole pairs produced in an event. Linking the two is a complicated chain of secondary scattering processes, which rapidly redistribute the energy deposited in the initial scattering.
A realistic treatment of the conversion from energy to ionization is a crucial step in calculating the sensitivity of experiments. Unfortunately, exact modeling of the secondary scattering processes is extremely challenging and is beyond the scope of this paper. Instead, we assume a linear response, which we believe does a reasonable job of capturing the true behavior. Specifically we assume that, in addition to the primary electron-hole pair produced by the initial scattering, one extra electron-hole pair is produced for every extra ε of energy deposited above the band-gap energy. Here ε is the mean energy per electron-hole pair as measured in high-energy recoils. The ionization Q is then given by .

(5.2)
We devote section 5.1 and appendix E to a discussion motivating this simple treatment. We emphasize that, while our treatment is approximate, it (a) is quite separate from the systematic, first-principles calculation of dR/dE e described in section 3 and section 4, and does not affect that calculation's accuracy; (b) is probably conservative, since it does not JHEP05(2016)046 account for fluctuations that could push a low-energy event above the ionization threshold; and (c) should be possible to improve upon in the future, both with better theoretical modeling and with experimental calibration.

Understanding the secondary scattering processes
It is experimentally well-established that for high energy electron recoils ( > ∼ keV), the ionization signal is directly proportional to the deposited energy, with a constant average energy ε deposited per electron-hole pair created, ε is several times the bandgap energy, accounting for the fact that only a fraction of the energy deposited goes directly into pair production. Fluctuations around the average ionization are quite small, with the Fano factor (defined as the ratio of the variance to the mean) measured to be [101,102] F At the low energies we are interested in, O(1-50 eV), the energy-ionization relationship has not been directly measured. Fortunately, there is reason to expect that the highenergy response can be extrapolated to lower energies. It has long been understood (see e.g. [100,103]) that following a high-energy electron recoil, an electronic cascade occurs that rapidly redistributes the energy between many low-energy electrons and holes. Roughly speaking, any electron or hole is expected to re-scatter and create an additional electronhole pair, so long as it has sufficient energy to do so. This repeats, distributing the energy over an exponentially increasing number of electron-hole pairs, until all electrons and holes have energy below the pair-creation threshold. Note that this threshold is larger than the band gap energy due to the constraints of momentum conservation [100]. The excess energy carried by the electrons and holes after the cascade is slowly lost to phonons, as is a fraction of the energy during the cascade. As a result of the cascade, the vast majority of secondary scatterings that occur after the initial electron recoil are low energy scatterings. This means that, for example, a single 10 keV electron recoil is approximately equivalent to 100 recoils with 100 eV each, or 1000 recoils with 10 eV each. This justifies the extrapolation of the high-energy behavior to low energies.
The linear response described by eq. (5.1) is not the only tractable approach. Other, less simplistic approaches can be taken without resorting to a full first-principles treatment. For comparison, in appendix E we construct a phenomenological Monte Carlo model of the secondary scattering cascade, following [103]. The model is intended to capture the general features of the cascade, without knowledge of the specific microscopic structure of the target material. The model reproduces the known high-energy behavior well with only a single tunable parameter, and can be used instead of eq. detectability. For typical masses, however, we find that the two approaches agree to within a few 10's of percent (see figure 17). We conclude that the linear treatment of eq. (5.1) is a reasonably realistic approximation, and postpone a more careful treatment to future work.

Results
In this section, we present the results of our calculation of the DM-electron scattering rates in silicon and germanium detectors. We show the potential reach for single-electronsensitive experiments, as well as the effect of higher experimental thresholds. We also give the full recoil spectra and the annual modulation fraction, which may be crucial for discriminating a possible signal from background. Lastly we discuss near-term prospects, focussing on upcoming searches expected from the SuperCDMS and DAMIC collaborations.
Experimental thresholds are set in terms of the ionization signal Q (the number of electron-hole pairs produced in an event) rather than the deposited energy E e . In the following results, we model the conversion of deposited energy to ionization with the linear treatment described in section 5. We take the DM halo to have a local density of ρ DM = 0.4 GeV/cm 3 [104,105], and a Maxwell-Boltzmann velocity distribution with a mean velocity v 0 = 230 km/s and escape velocity v esc = 600 km/s, and we take the average Earth velocity to be v E = 240 km/s (see appendix B for explicit formulae). In appendix C we discuss the numerical convergence of our results.
Event rates as a function of Q, for an extensive range of DM masses, are available online at this link. The crystal form-factor, as a function of q and E e , is also available there. Using eq. (3.13), the information online can be used to re-derive rates using a different DM form-factor or velocity profile, or using a different treatment of the energy-to-ionization conversion.

The crystal form factor
Much of the behavior of the scattering rates can be understood from the behavior of the crystal form factor, |f crystal (q, E e )| 2 , in eq. (3.17). We show the crystal form factor in figure 5 as a function of q and E e , for both silicon and germanium. The rapid fall-off as q increases is clearly visible. The solid line in the figure corresponds to v min = v esc + v E from eq. (3.12) as m χ → ∞. The region below this line is kinematically inaccessible for any DM mass. The dashed line uses the velocity of a typical DM particle in the halo, i.e. v min = 300 km/s. We see that larger recoil energies require larger q, for which the crystal form factor is suppressed. The implication of this is that the DM-electron scattering rates increase dramatically for smaller recoil energies, resulting in a dramatic increase in sensitivity as detector thresholds are lowered.  .17)). In the region below the solid line, v min > v esc + v E for any DM mass, and electron scattering is thus kinematically inaccessible. The dashed line corresponds to v min = 300 km/s (a typical DM halo velocity) in the heavy DM limit; the region below this line is only kinematically accessible to DM particles with velocities larger than the average velocity. For energies above ∼ 10 eV, the scattering rate is suppressed by both the form factor and DM velocity distribution. We see that the 3d electrons in germanium give a sizable contribution to |f crystal (q, E e )| 2 for E e > 25 eV. thresholds, Q th , of 1, 5, and 10 detected electron-hole pairs, respectively, which correspond to deposited energies, E e , of 0.67, 12.3, and 26.8 eV in germanium, and 1.1, 15.5, and 33.5 eV in silicon (to get the corresponding ionization energy thresholds, subtract 0.67 eV for germanium and 1.1 eV for silicon from these numbers, see eq. (5.1)). The three plots show results for different DM form factors, corresponding to different classes of DM models: F DM (q) = 1 (top left), F DM (q) = αm e /q (top right), and F DM (q) = (αm e /q) 2 (bottom), see section 2 for details.
As expected, the reach dramatically improves when the threshold is lowered, since the crystal form factor strongly suppresses the rate for high electron recoil energies. This improvement is most pronounced for F DM (q) = (αm e /q) 2 , since lower q tends to correspond to lower recoil energies. With a single-electron threshold, the difference in sensitivity for silicon and germanium targets can be accounted for by the fact that germanium is 2.6 times heavier, and so has correspondingly fewer valence electrons per kg. However, germanium targets are sensitive to slightly lower DM masses due to their lower band-gap. In addition, germanium targets become comparably more sensitive than silicon targets for ionization thresholds of Q th 9 due to the additional contribution from the 3d-shell electrons (see below). Figure 7 shows the spectrum of events as a function of the ionization signal Q, for different DM form-factors and masses, in silicon (left) and germanium (right) targets. The fast fall-off with Q shows the large gain to be made from lowering experimental thresholds JHEP05(2016)046 Figure 6. Dark matter-electron cross-section sensitivity: the 95% C.L. exclusion reach in the DM-electron scattering cross-section, σ e , of an experiment with 1 kg-year exposure and zero background events, for different experimental thresholds. Solid (dashed) lines show the reach for silicon (germanium) targets. Ionization thresholds of 1, 5, and 10 electron-hole pairs are shown with blue, green, and red lines, respectively. The corresponding energy thresholds are 0, 11.6, and 26.1 eV in germanium, and 0, 14.4, and 32.4 eV in silicon. The gray shaded region shows the existing constraint from XENON10 data [31]. The three plots assume different DM form factors, F DM (q) = 1, αm e /q, (αm e /q) 2 , corresponding to different DM models. towards a 1-or 2-electron threshold, especially for the steeper DM form-factors and for lower DM masses. The shape of the spectra may be useful in discriminating a signal from background.
In germanium, the 3d-shell electrons dominate the rate for E e 24 eV, corresponding to Q 9 electron-hole pairs (the 3d-shell electrons lie about ∼15 eV below the bottom of the valence band and ∼24 eV below the bottom of the conduction band). The intuitive reason for the 3d-shell electrons dominating over the valence electrons at large E e can be seen from eqs. (3.5) and (3.6). The typical velocity of the 3d-shell electrons, and hence the typical momentum transferred from the DM, is larger than for the valence electrons, so that the 3d-shell electrons can dominate if they are kinematically accessible. 6 We show the effect of neglecting the 3d-shell electrons in figure  cross-section reach and the recoil spectrum generated by DM scattering with and without the inclusion of the 3d-shell electrons. The effect is significant for ionization thresholds above Q th ≈ 7 or 8, but not important at lower thresholds.
We note that there are some differences between our results and those in [9,25]. In [9], only the case Q th = 1 was considered and we find that the new computation predicts a somewhat lower rate. We find that the shape of the recoil spectra in [25] is noticeably different from ours, which gives rise to several differences in the expected limits for the different Q thresholds. Furthermore, for germanium, we also include the 3d-shell electrons, which can be important as discussed above.

Comparison with existing XENON10 limit and discussion of background
We see from figure 6 that to surpass the existing limits obtained with XENON10 data [31], a germanium-or silicon-based experiment with an ionization threshold of 10 electrons would require a background-free exposure of less than 1 kg-year. However, with a single-electron threshold, such an experiment would surpass the XENON10 limit at all masses with a background-free exposure of around 1 kg-day for F DM (q) = 1, 10 g-days for F DM (q) = αm e /q, or just a 1 g-day for F DM (q) = (αm e /q) 2 . In addition, with any exposure such an experiment would place the first ever bounds in the ∼1-5 MeV mass range, below the threshold of the XENON10 search. The XENON10 detector had a threshold of one electron with an O(1) detection efficiency, but to obtain one electron required an energy of at least JHEP05(2016)046 12.4 eV to overcome the binding energy of an electron in the outer shell. 7 Moreover, the background in the XENON10 data was much larger than conventional nuclear-recoil background, so that the number of DM events leading to single (two, and three) electrons was only limited, at 90% C.L., to be less than 8,550 (1,550, and 330) counts/kg/year, respectively.
While the single-electron background in the XENON10 data was rather large, its origin is likely specific to its dual-phase detector setup. Many of the single electron backgrounds likely had one, or a combination, of the following origins [31]: (i) electrons, trapped in the potential barrier at the liquid-gas interface, were randomly drawn into the gas phase (these transiently trapped electrons likely originated from other background events that caused xenon atoms to ionize); (ii) photo-dissociation of a negatively charged O 2 -ion, which received its negative charge from a drifting electron that originated from another event; (iii) field emission from the cathode. XENON100 and LUX may face similar challenges, although an analysis is still in progress.
The semi-conductor targets will not suffer from these same detector-specific backgrounds. They will, of course, have their own unique experimental challenges to deal with, including detector noise and dark current, as we will discuss in more detail in section 6.5 for DAMIC and SuperCDMS. These will likely be the limiting instrumental factors in setting the threshold for a particular experiment. Once these challenges are overcome, one needs to deal with the physics backgrounds. As argued in [9], neutrinos are not an important source of background even for the largest exposures considered in this paper (O(20 kg-years) for SuperCDMS, see section 6.5). Compton scattering or other events that produce recoiling electrons will usually lead to a much larger energy deposition and most of them could thus be vetoed, although some backgrounds will persist to the lowest energies. The size of this background will depend on the shielding and purity of the materials around the detector; for SuperCDMS at SNOLAB, the Compton background is estimated to be O(6 × 10 −3 ) events/kg/day/keV [106]. Assuming that it is flat at low energies, this translates into O(0.04) events/eV for 20 kg-years, which would be negligible. We do not expect there to be backgrounds from neutrons, and (cosmogenic) x-ray lines will lie well above our energies of interest. Surface events and other, unknown, backgrounds may be present at low energies. As experiments reach the required sensitivity to probe the few-electron events expected from LDM scattering off electrons, a better understanding of all backgrounds will emerge allowing for an attempt to mitigate them if necessary. It is possible that a spectral analysis of a signal will further allow for the removal of some background events. Our assumption of zero background for the plots should be taken as the best-case scenario.

Annual modulation
Even with a significant background event rate, it may be possible to distinguish a signal from background using annual modulation, as long as the background is stable on year timescales. Annual modulation is a distinguishing feature of a LDM scattering signal [9,107] occurring due to the change in the earth's velocity through the DM halo as it rotates around the Sun. For a standard smooth and isotropic DM velocity distribution, the modulation is approximately sinusoidal with year period and a peak around June 2nd (the presence of DM streams or non-trivial DM structure may complicate this, as may gravitational focusing by the Sun [68,108], which we do not include). The modulation fraction, f mod , is defined to be the ratio of the amplitude of the modulating signal to the median signal rate, which for our assumed halo profile (see appendix B) corresponds to where R i represents the rate at time of year i.
For DM scattering off electrons, the modulation fraction can be significantly larger than for the usual elastic scattering of (heavy) WIMPs off nuclei. As we saw in figure 5 (see discussion in section 3.1), DM-electron scattering relies on the tail of the DM velocity distribution, especially for energies above ∼ 5-10 eV. We plot the modulation fraction in figure 8. The left and center plots show f mod as a function of ionization Q for different masses and DM form-factors for the two elements. f mod rises from a few percent for singleelectron events to above 10% for events with more than ∼ 3 electrons. Comparing with the spectrum in figure 7, we see that there is a trade-off between modulation fraction and event rate. Events with several electron-hole pairs provide large modulation without sacrificing too much in the rate, and may give the best prospects for annual modulation searches depending on the background. The modulation fraction also rises near the mass threshold, as we show on the right of figure 8 for ionization thresholds of Q th = 1 and 5. Note that the high-mass value of f mod for the single-electron threshold, at 4-6%, is larger than the JHEP05(2016)046 Figure 9. Annual modulation: discovery reach: the 5σ discovery reach in the m χ -σ e plane of an annual modulation search for DM-electron scattering for an experiment with a 1 kg-year background-free exposure. Solid (dashed) lines show the reach for silicon (germanium) targets. Ionization thresholds of 1, 5, and 10 electronhole pairs are shown with blue, green, and red lines, respectively. The corresponding energy thresholds are 0, 12.3, and 26.8 eV in germanium, and 0, 15.5, and 33.5 eV in silicon. The gray shaded region shows the existing constraint from XENON10 data [31]. The three plots assume different DM form factors, as indicated, corresponding to different DM models. Thin lines are from figure 6, showing the exclusion reach of a search with the same exposure seeing no events. values in the Q = 1 bin of the left and center plots, because the total rate is not dominated by the single-electron events.
Once a signal is found in an electron scattering search, increasing the exposure of the experiment until the annual modulation can be tested will be a crucial step in claiming a DM discovery. In figure 9 we show the 5σ discovery reach of an annual modulation search in the mass-cross-section plane. We calculate this cross section 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. The thick curves in figure 9 show the discovery reach for different thresholds and DM form-factors, assuming a background-free exposure of 1 kg-year. (A non-zero background will of course weaken the reach, following the equation above.) This figure mirrors figure 6, which shows the exclusion reach obtained using a simple counting search instead of a modulation search, but otherwise with the same assumptions. The curves of figure 6 are replotted as the thin curves in figure 9, for comparison. We see that a substantial discovery reach is possible with a 1 kg-year exposure: at low masses for F DM (q) = 1, and at all masses for F DM (q) = αm e /q or (αm e /q) 2 . Finally, we comment that taking into account the directional (sub-daily) modulation, which is expected in crystalline detectors, will further allow for an improved sensitivity to a DM signal. As discussed in section 3.2.2, we have averaged-out such directional effects in this work, and we postpone their study to future work.

Prospects for near-term experiment
In this subsection, we discuss the near-term prospects for electron-scattering searches with the DAMIC and SuperCDMS experiments.

DAMIC
DAMIC [78,109,110] uses thick, fully-depleted silicon CCDs for their target material. These CCDs are ten times more massive than conventional CCDs, allowing them to be competitive targets for DM direct detection. In [78], DAMIC used one 0.5 g CCD to perform an engineering run, obtaining an exposure of 107 g-days. They were able to constrain DM-nuclear scattering for DM masses almost as low as 1 GeV. Work is ongoing to increase the total mass of the detector (by using more CCDs) as well as the detector's sensitivity to low threshold energies (by using so-called "Skipper CCDs") [67].
The first direct detection limit using a semi-conductor target. Here we investigate the (albeit weak) constraints on DM-electron scattering from their existing result, and give reach estimates based on their projected detector improvements. For the engineering run [78], DAMIC used a single 0.5 g CCD, for an exposure of 107 g-days. They obtained the following values for the read-out noise and the dark current: (i) A readout noise of below 2 electrons/pixel, corresponding to 2 × 3.6 eV = 7.2 eV of r.m.s. readout noise. The CCD has about 4.2 million pixels, so that one requires a threshold of ∼ 13 electrons (∼ 47 eV) for the noise to produce a signal above threshold in less than one pixel. DAMIC chose a threshold of ∼ 40 eV (∼ 11 electrons) for the search for DM-nuclear scattering; we expect ∼ 35 pixels to reach this threshold. In our recast of their data for DM-electron scattering below, we will use the same 40 eV threshold.
(ii) A dark current of ∼ 1 electron/CCD/day (at the chosen 120 K operating temperature). Since the exposure of the CCD is a few hours, before being read-out within a few minutes, the threshold is limited by the read-out noise, and not the dark current.
We can use the result in [78] to constrain DM-electron scattering. We will assume that the efficiency to select electron recoil events is the same as selecting nuclear recoils, i.e. 7×10 −3 . Figure 12 in [78] shows the data that was recorded by DAMIC, and we see that zero events

JHEP05(2016)046
were recorded in the first bin above the threshold (40 eV to 100 eV). This may be as a result of the efficiency being very low in this bin; nevertheless, it could also be a sign that backgrounds may be small at such low energies, boding well for future runs with even lower energy thresholds. In any case, this information is sufficient to derive the current DAMIC limit on LDM, which we show with a green shaded region (bounded by a green line) in figure 1. We see that with the current threshold, the form-factor suppression is too large for this constraint to compete with the existing XENON10-based limit [31]. Nevertheless, this is the first direct detection limit for sub-GeV DM using a semi-conductor target.
Projections for future DAMIC runs with improved "Skipper" CCDs. There are two main challenges that need to be overcome by DAMIC (and similar experiments) to allow them to push to low thresholds [67]: (i) reduce the noise in reading out the ionization deposited in the detector, and (ii) reduce the dark current. The read-out noise can be reduced substantially by taking more time to read the CCD, while the dark current (i.e. genuine electron-hole pairs produced by thermal excitations in the silicon substrate) can be reduced by lowering the temperature and improving the quality of the silicon. The contribution from the dark current will increase with the readout time, so it will take some optimization to find a way to reduce the readout noise while keeping the contribution from the dark current manageable. Lowering the temperature also reduces the electron mobility in the substrate, requiring a careful trade-off.
Here we project what future data runs can achieve with the improved DAMIC Skipper CCDs. The DAMIC Collaboration has been working on so-called "Skipper CCDs", which will reduce the r.m.s. read-out noise down to 0.2 electrons/pixel/day, with the possibility of going down to 0.1 electrons/pixel/day [67,78]. This is done with a new output circuit that enables multiple read-outs. The size of the CCDs can be anything up to 4 × 4 = 16 million pixels (Mpix) [111], but we will assume 8 Mpix for the projections below. The 0.2 (0.1) electrons/pixel/day correspond to read-out noise of 0.72 (0.36) eV; assuming gaussian noise, the 0.1 electron/pixel/day will allow for sensitivity down to single electrons, although the 0.2 electrons/pixel/day may require a threshold of two electrons to avoid the noise faking a DM signal. However, for such low read-out noise, the dark current becomes the limiting factor in determining the energy threshold. Significant non-gaussian tails could change this conclusion.
The dark current has been measured currently at 5×10 −3 electrons/pixel/day [111]. As mentioned above, lowering the operating temperature and improving the silicon substrate quality will decrease this, and it is reasonable to expect further improvement; the theoretical lower limit is O(10 −7 ) [111]. Below we provide projections under both assumptions.
The effect of the dark current on the threshold depends on the number of pixels and the exposure length of the CCD. The CCD is read pixel-by-pixel, and can be read continuously from one side to the other, before cycling back again to the beginning. We will assume that the 8 million pixels of the CCD are all read in one hour, so that its exposure is one hour for the purposes of calculating the dark current. We will consider the following two scenarios for the number of CCDs, the mass, and exposure (we assume an efficiency of 1 for making our projections below):  Table 1. Expected number of events with at least Q th electron-hole pairs under different assumptions for the dark current and exposure (5 × 10 −3 electrons/pixel/day and 10 −7 electrons/pixel/day) and assuming either (i) 4 CCDs (10 g) and an exposure of 1 kg-day; or (ii) 40 CCDs (100 g) and an exposure of 100 g-years. In both cases, we assume that it takes one hour to read the entire CCD, and that it is read continuously. Projected exclusion reach based on a simple counting experiment are given in figure 1 for the entries marked with a single star ( ). A projected discovery reach based on seeing the annual modulation of the signal, with negligible background, is also given in figure 1 for the entry marked with a dagger ( †) (we find that the prospects from annual modulation for the entry marked with a diamond ( ) are very similar). See text for details.

JHEP05(2016)046
(i) There are currently four prototype skipper CCDs, each with a mass of 2.5 g, which were produced as part of an R&D project (these will be deployed at the MINOS near site this year). For our first set of projections, we will assume that data is taken over 100 days (livetime), for a total exposure of 1kg-day.
(ii) If the testing of the skipper CCDs at MINOS goes well, one can expect that several more of them will be deployed to search for DM. Thus, for our second set of projections, we will assume that 40 CCDs are deployed (for a total mass of 100 g) and that data is taken again over 365 days (livetime) for a total exposure of 100 g-years. Table 1 gives the expected number of events with at least Q electron-hole pairs for the two scenarios assuming Poisson statistics. We see that DAMIC could have a threshold of 2 electrons if the dark current can be reduced to below 10 −7 electrons/pixel/day, but a 3-electron threshold is required for the present dark current rate of 5 × 10 −3 electrons/pixel/day. In figure 1, we show solid green lines that indicate the 95% C.L. prospects for the entries marked with a star ( ), i.e. we show the cross section to obtain 3.6 signal events, assuming zero background events. We also show the reach of an annual modulation search for the entry marked with a dagger ( †), i.e. a 2-electron threshold with no background in a 100 g-year exposure. This is shown by the dashed green line in figure 1. We checked that the prospects from annual modulation for the entry marked with diamond ( ) are very similar. We see from these projections that DAMIC can significantly improve upon the current XENON10 limit, especially at the lowest DM masses.

(Super)CDMS
The CDMS experiment uses cryogenic solid state detectors operated at temperatures below 100 mK. In the WIMP search, they distinguish electron from nuclear recoils by measuring the ratio of the ionization versus phonon energy deposited into the crystal. This ratio will be smaller for nuclear recoils than for electron recoils. Here, we are interested in their ability to detect electron recoils.

JHEP05(2016)046
The signal from a low-energy recoiling electron can be dramatically enhanced by applying a relatively large bias voltage, V b ∼ O(50-100 V), across the target material. The work done in drifting an electron-hole pair out of the crystal, eV b , is emitted as Luke-Neganov [112][113][114] phonons, which will be picked up by the phonon sensors. This was done for "CDMSlite", which has yielded electron recoil thresholds with an O(1) detection efficiency of 170 eV (i.e. O(50) electrons) [92]. Here we discuss the prospects of future versions of SuperCDMS.
In figure 1, we show 3 projections for SuperCDMS, two for silicon (with an exposure ∼10 kg-years) assuming an electron-hole pair threshold, Q th , of either 4 or 1, and a signal detection efficiency of 0.7, and one for an annual modulation search assuming Q th = 2 with the same exposure and efficiency. The cross-section reach for germanium is very similar to those of silicon.
The electron-hole-pair thresholds are based on the following assumptions. The Q th = 4 threshold is based on the numbers used by SuperCDMS for Snowmass [63], while the Q th = 1 threshold is based on an ambitious but achievable best-case scenario. For Snowmass, a phonon energy resolution of 50 eV was assumed. As there may be non-gaussian tails, a 7σ threshold was assumed, corresponding to a threshold of E e = 350 eV. Taking the bias voltage to be 100 V, this translates into Q th = 3.5, which we round up to 4. For the second set of projections, we assume that further R&D can push the noise threshold down to better than ∼14 eV, which is ambitious but achievable in principle [115]. A 7σ threshold corresponds to ∼100 eV, so that a bias voltage of 100 V is sufficient to achieve sensitivity to Q th = 1. In practice, the bias voltage can be optimized as well. A larger bias voltage would allow for a reduced threshold in terms of the number of electron-hole pairs needed to pass the phonon energy threshold, but could lead to breakdowns of the substrate. However, it has been demonstrated that the bias voltages needed for sensitivity to Q th = 1 are achievable in both silicon and germanium [66].
As can be seen from figure 1, SuperCDMS has the potential to improve drastically upon the existing XENON10 limit, especially at low DM masses.

Conclusions
Direct detection experiments have so far primarily focused on searching for WIMPs, and as a result of an intense research effort, the path forward in this direction is rather well-defined. Within the next decade, WIMPs will either be found or become significantly less motivated. However, other theoretically motivated candidates exist that could constitute the DM in our Universe. In this work, we focused on a class of DM candidates that have a mass between a few-hundred keV to a GeV. We showed that tremendous progress can be made in exploring the direct-detection parameter space of these candidates over the next few years, by searching for DM-induced electron recoils in experiments with targets that consist of semiconductor materials. The technology currently used in WIMP searches can be adapted for such light-DM searches by improving the ionization sensitivity, and this is being actively pursued. The backgrounds are expected to be quite different in nature to those in WIMP searches, and there is reason to believe that they will be small and controllable.

JHEP05(2016)046
The calculation of the DM-electron scattering rate and the subsequent electron recoil spectrum in semiconductor targets is much more challenging than for DM-nuclei scattering. We have provided detailed formulae for the scattering rate and recoil spectrum, expressed in terms of simple DM properties and a target-dependent "crystal form factor", which encodes the quantum structure of the target electrons. We numerically calculated the crystal form factor for germanium and silicon with our code QEdark, which is based on the software package Quantum Espresso that calculates the crystal wave functions and energy levels.
Convergence tests indicate that our results are accurate at the few percent level. QEdark will be publicly available at http://ddldm.physics.sunysb.edu, together with the crystal form factors. With these, upcoming experiments can derive their own sensitivities or limits.
The crystal form factor is a steeply falling function of the electron recoil energy. Consequently, even a small improvement in an experiment's detector threshold translates into a significant increase in the sensitivity to DM-electron scattering. We have provided the projected sensitivity for a variety of experimental thresholds, showing that upcoming experiments including DAMIC and SuperCDMS can probe orders of magnitude of unexplored DM parameter space in the next few years. In addition to setting limits, sub-GeV dark matter can be discovered via its expected modulation signal. We showed that in the case of electron-scattering, annual modulation is sizable and could provide the necessary signal for discovery. Additional sub-daily modulation is expected due to the orientation-dependent nature of scattering in crystalline detectors. We have ignored directionality in this work, deferring it to future study.
Calculating the experimentally observable signal requires a conversion from energy deposition to the ultimate ionization signal. This conversion requires a detailed knowledge of the secondary scattering processes in crystals, at energies below the existing experimental sensitivity. We therefore used a phenomenological model for secondary interactions, and studied its possible systematic uncertainties using a Monte Carlo model. We find that our predictions suffer from systematic uncertainties of order a few tens of percent, and is likely conservative. Further theoretical and experimental study of secondary interactions would be useful to improve the modeling of this conversion.
To summarize, our work provides the necessary tools for experiments which use semiconductor targets to search for sub-GeV dark matter to derive accurate limits. Technologies adapted from WIMP searches and currently under development can be employed in searches for sub-GeV dark matter. This highly-motivated direction in dark matter searches is a natural progression from the WIMP program, and we expect that it will take a leading role in the search for dark matter.
Note added. While this work was being completed, ref. [68] appeared, which also deals with DM-electron scattering in germanium. Ref. [68] is complementary to our work, its main point being the effect of "gravitational focusing" on the modulation signal of DMelectron scattering. Ref. [68] derives scattering rates using a semi-analytic approach, which builds on the method of ref. [25], but is significantly less detailed than the method we have presented here. We find that their results are comparable to ours within a factor of a few, but with some notable differences. In particular, ref. [68] finds increasingly higher rates JHEP05(2016)046 than us at increasingly higher recoil energies. Most strikingly, ref. [68] finds that scattering of the 3d shell electrons dominates the total rate when it is kinematically accessible (we find the the 3d shells cause a bump in the spectrum, but with a rate subdominant to lower energy events). We attribute these differences to the inherent sensitivity of the calculation to the tails of the electron wavefunctions for energies above O(10 eV), as we discussed in section 3.1.

A Derivation of scattering rate formulae
A.1 General formula for dark matter-induced electronic transitions If a DM particle scatters with an electron in a stationary bound state such as in an atom, it can excite the electron from an initial energy level 1 to an excited energy level 2, by transferring to it energy ∆E 1→2 and momentum q. The cross section for this process can be derived in a standard way using non-relativistic quantum mechanics, but here we derive it starting from the usual formula for the cross-section in field-theory, in order to make easier connection with the underlying particle physics. We treat the electron as being bound in a static background potential -in other words we approximate the atoms as being infinitely heavy objects which can absorb momentum without recoiling. This is an excellent approximation (< 1% error), since the momentum-transfers we will be interested are typically of order keV.

JHEP05(2016)046
The cross section for free 2 → 2 scattering is given by where M free is the usual field-theory matrix element, and |M| 2 represents its absolute square averaged over initial spins and summed over final spins. If the electron were unbound, the non-relativistic scattering amplitude would be given by where |χ p , e k is plane-wave state for a DM particle of momentum p and an electron of momentum k, H int is the interaction Hamiltonian, and C is an unimportant coefficient. However, because the electron is bound it is instead given by where ψ 1 , ψ 2 are the (unit normalized) momentum-space wavefunctions of the initial and final electron levels. We have purposefully used plane-wave normalization for both the free and bound electron states: e k |e k = e 1 |e 1 = (2π) 3 δ 3 ( 0) ≡ V , where V is the volume of space (which always cancels in the end). To find the cross section for this excitation process, we can use the free 2 → 2 scattering cross section formula but with two replacements: one to account for the modified scattering amplitude, and the other to account for the different final-state phase space. Squaring eqs. (A.2), (A.3), we see that the bound-state scattering amplitude is accounted for by making the replacement where f 1→2 ( q ) is the atomic form factor, Fourier transforming eq. (A.5) gives the definition given in eq. (3.10).
Since there is only one final electron state being considered, we also need to remove the usual final-state phase space integral: Combining eqs. (A.1), (A.4), (A.6), we can write the formula for the cross-section for a DM particle to excite an electron from level 1 to level 2:

JHEP05(2016)046
Since we are in the non-relativistic regime, the energies are given by Using the following definitions [9] to parametrize the underlying DM-electron coupling the cross section simplifies to A.2 Average rate in a dark matter halo The actual rate of excitation events, for a given transition and a given target electron, is found by multiplying eq. (A.12) by the DM density and averaging over the DM velocity distribution g χ ( v), In general, both the electron wavefunctions and the DM velocity distribution will not be spherically symmetric. As noted in [9], the rate will then depend on the orientation of the target with respect to the galaxy. Here we ignore this interesting complication, and approximate the velocity distribution as being spherically symmetric. We can then use the d 3 v integral to eliminate the δ-function in eq. (A.12), giving Here η(v min ) has its usual definition, and v min is a function of q and the energy transfer given by

A.3 Ionizing an isolated atom
For the purposes of connecting with previous work [9], in this subsection we consider ionization of electrons bound in isolated atomic potentials. We derive the ionization rate of such a system, assuming a spherical atomic potential and filled shells. This approximation was used in [9] to model a liquid xenon target material, and the results below reproduce eqs. (5) and (6) of that paper. The full calculation of event rates in liquid xenon would require knowledge of electron wavefunctions in the dense, disordered xenon liquid. This is a more challenging calculation than for a semiconductor crystal, but can in principle be performed with similar methods -we leave this for future work. The corrections, however, can be argued to be small, lowering the ionization threshold and increasing the event rate. An electron ionized from an atom can be treated as being in one of a continuum of positive-energy bound states. These states are affected by the potential well of the atom, but can be approximated as free particle states at asymptotically large radii. We denote their wavefunctions as ψ k l m ( x), where l and m are angular quantum numbers, and k is the momentum at asymptotically large radius. The energy of such a state is therefore The ionization rate for such an atom is found by taking eq. (A.14), summing over occupied electron shells, and integrating over the phase space of all possible ionized states. Since these are asymptotically free spherical-wave states, the phase space is when the wavefunction normalization is, as in [9], taken to be Plugging this in, the ionization rate is given by where E Bi is the binding energy of occupied state i.
Since the potential is assumed to be spherically symmetric, and we are ionizing a full atomic shell, we can sum |f 1→k l m ( q )| 2 over initial and final angular momentum variables and the result cannot depend on the direction of q. This means we can define the dimensionless ionization form factor,

JHEP05(2016)046
After applying this definition to the previous equation, we can replace the d 3 q integral with 4πq 2 dq, giving This reproduces the formulae given in [9].

A.4 Excitations in a semiconductor crystal
In the periodic lattice of a semiconductor crystal, the electron energy levels form a complicated band structure, with an energy gap separating the filled valence bands and the unoccupied conduction bands (figure 4). Each possible electron level is labelled by a band index i and a wavevector k in the first Brillouin Zone (BZ). Due to the periodicity of the potential, the wavefunctions of these states are in Bloch form, where the G's are the reciprocal lattice vectors. Here V is the volume of the crystal, and the wavefunctions are taken to be unit-normalized, so that With this form for the wavefunctions, the form factor eq. (3.10) to excite from valence level {i k} to conduction level {i k } becomes We define the term in the absolute square in eq. (A.25) to be f [i k,i k , G ] : Inserting this into eq. (A.14), we can use the δ-function to eliminate the d 3 q integral, giving (A.27) The total excitation rate for an electron in level {i k} is found by summing eq. (A.27) over all unfilled final energy levels i ,

JHEP05(2016)046
Note that we do not sum over final electron spins here as that sum has already been included in the definition of σ e . The total rate of excitation events in the crystal, R crystal , is given by summing eq. (A.28) over all filled initial levels i, Here the extra factor of 2 is the sum over the two degenerate spin states of the filled valence bands.
Putting this together gives the total excitation rate in a crystal, (A.30) where again q = | k + G − k|. Note that this is the total event rate for the whole crystal, and so it is appropriate that it is proportional to the volume V of the whole crystal. Since the dependence on the DM velocity distribution and interaction type are entirely encoded in η and F DM , which are functions only of the momentum transfer q and energy deposited E e , it is useful to insert delta-functions into the above expression as follows: Using V = N cell V cell , where V cell is the volume of the crystal's unit cell and N cell is the number of cells, the differential rate can then be written in the form of eq. (3.13), (A.32) where the crystal form-factor is defined as in eq. (3.17), In this section, we will derive analytic expressions for η(v min ). For simplicity we assume that the DM velocity distribution, g χ ( v χ ), takes the form of a Maxwell-Boltzmann distribution in the galactic rest frame, with a hard cutoff at the galactic escape velocity. In the Earth's frame the velocity distribution then takes the form

JHEP05(2016)046
where v χ is the DM velocity in the Earth frame, and v E is the Earth's velocity in the galactic rest frame. We take v 0 = 230 km/s for the typical velocity, and v esc = 600 km/s for the escape velocity. We take v E = 240 km/s for the average Earth velocity relative to the DM halo, adding (subtracting) 15 km/s for the Earth velocity in June (December). The normalization factor K is determined by requiring d 3 vg χ ( v) = 1, giving Using these values, we obtain K = 6.75 × 10 22 [cm/s] 3 or 2.50 × 10 −9 in natural units. We then define the function η(v min ), where c θ = cos θ is the angle between the velocity and the velocity of the Earth. We can explicitly solve eq. (B.3), but need to consider two cases: This gives us where the subscript corresponds to the case number. Note that the two cases converge to the same value for v min = v esc − v E .

C Convergence of the numerical results
In this section, we investigate the dependence of our calculation on the kinetic-energy cutoff, E cut , (see eq. (4.3)) and on the k-point mesh. The choice of E cut determines the maximum allowed three-momentum transfer q, which impacts the maximum E e that we sample. Truncating the range of q can have more of an effect for high DM masses and high electron thresholds, since these two regimes depend on high values of q. On the other hand, the k-points included in the mesh determine the computationally allowed values of q. A higher-resolution mesh is particularly important for low E e transitions to the bottom of the conduction band, and is therefore especially important for low DM masses and low electron thresholds.  Figure 10. The minimum velocity, v min as a function of q for m χ = 1 MeV, 10 MeV, and 1 GeV and E e = 1, 10, and 40 eV. The grey shaded region indicates where v min > v E + v esc and is thus prohibited. The vertical grey line indicates the maximum q value for a given E cut = 70 Ry.  Figure 11. Left: cross-section sensitivities for ionization thresholds of Q th = 5 and Q th = 11 electrons in silicon for E cut = 30, 50, 70, 90, and 110 Ry (we take a mesh consisting of 27 k-points). Note that most lines are on top of each other, demonstrating the weak dependence of σ e on E cut . Right: difference in the rate, R, for a given E cut , to that at E cut = 70 Ry, R 0 , in silicon for m χ = 1 GeV. We see that the error in choosing E cut = 70 Ry is < O(1%) for the thresholds considered in this paper.
In figure 10, we show the dependence of v min in eq. (3.12) on q for different m χ and E e . The choice of E cut (top axis) determines the range of q (bottom axis). In figure 11, we show the ratio of the rate for different values of E cut to the rate at E cut = 70 Ry for m χ = 1 GeV for silicon targets. We see that the error in the rate with our choice of E cut = 70 Ry is < O(1%). The E cut convergence in germanium is slightly worse due to the fact that we are solving for the 3d electrons instead of including them in the pseudopotential. This effect is greatest near the 3d shell energies (a few percent uncertainty), as seen in figure 12. In left plots of both figure 11 and 12, we see a step-like transition between 30 Ry and the other curves for 11e because 30 Ry is not a high enough energy cutoff to accurately calculate the JHEP05(2016)046  Figure 12. Left: cross-section sensitivities for ionization thresholds of Q th = 5 and Q th = 11 electrons in germanium for E cut = 30, 50, 70, 90, and 110 Ry (we take a mesh consisting of 27 kpoints). Note that most lines are on top of each other, demonstrating the weak dependence of σ e on E cut . Right: difference in the rate, R, for a given E cut , to that at E cut = 70 Ry, R 0 , in germanium for m χ = 1 GeV. The structure of the distributions arise from the effect of the 3d electrons. rate for 11e. The irregular behavior in the distributions on the right side of figure 12 are from the semicore electrons in Ge. We do not see the same behavior in figure 11, which considers Si and no semicore electrons. We investigate the effects of our choice of k-point mesh on our results in two ways. 8 First, we vary the number of k-points in our mesh and find that there is sensitivity to our choice at low masses and thresholds, see figure 13. Second, we perturb each point on the  Figure 14. The effects of our choice of k-point mesh on our germanium results by perturbing the mesh with random shakes of amplitudes up to half the lattice-spacing. On the left, we look at the standard deviation of the shaken runs over the mean value as a function of DM mass. On the right, we look at the standard deviation of the shaken runs over the mean value as a function of E e for m χ = 5 MeV. We see that our choice of k-mesh spacing is accurate to a few tens of percent for masses above 1 MeV.
mesh with a random shake of amplitude up to half the lattice-spacing so as to cover the entire k-space. We use an energy cutoff of E cut = 70 Ry and 243 k-points. The amplitude of our perturbations is ∆k = 0.08 a.u. as the lattice spacing for 243 k-points is 0.17 a.u. We run 10 independent simulations and plot the results in figure 14. We find that our choice of k-point mesh does not appreciably affect our results for masses above ∼ 1 MeV.

D The importance of the 3d-shell in germanium
The importance of the 3d-shell electrons in germanium are illustrated in figure 15. We see that they dominate the rate at high recoil energies and thus for high thresholds. We discuss this in more detail in section 6.2.

E A Monte Carlo model of secondary scattering
In the main results of this paper, we modeled the ionization response of a target crystal with the linear treatment described in section 5. For comparison, here we attempt to mock-up the secondary scattering with a Monte Carlo model, following [103]. The deposited energy E e is randomly split between an initial electron and hole. In each following step, each electron or hole with energy above a threshold E ion then generates an extra electron-hole pair, with the energy being randomly split between the three particles. This is iterated until all particles have energy less than E ion . The random energy splittings follows a distribution that weights all phase space volume equally, with the density of states assumed to grow as √ E above and below the bandgap, as in a simple 2-band free electron/hole system. Explicitly, for the initial 1 → 2 splitting the probability distribution for energy E 0 to split into energies E 1 and E 2 has the form while for the subsequent 1 → 3 splittings it has the form where electon/hole energies are measured above/below the upper/lower edge of the band gap. We ignore phonon losses during the cascade -these are understood to be quantitatively fairly small, and should not affect the qualitative conclusions. The output of the Monte Carlo model is a probability distribution P (Q|E e ) to get ionization Q given a deposited energy E e . Given the band-gap energy of E gap = 0.67 eV (1.11 eV) in germanium (silicon), we find that E ion = 2.67 eV (3.1 eV) reproduces the measured values of ε for high energy recoils (see eq. (5.2)). The distributions for both elements have Fano factors of F ≈ 0.1 for all energies above ∼ 10 eV, consistent with measurements. We illustrate the probability distributions in figure 16. In figure 17, we show the effect on the event rate of using this model rather than the naive linear model of section 5. For thresholds of 2 to 4 electron-hole pairs, downward fluctuations reduce the rate compared to the naive estimate. For higher thresholds, occasional upward fluctuations combined with the steeply falling recoil spectrum lead to an increase in the rate. However, the two models are consistent within a few tens of percent.

F Review of Density Functional Theory and pseudopotentials
In this appendix, we review the formalism of density functional theory (DFT), explain in more detail the approximations used in the computation of the wavefunctions, and further explain the numerical methods.   where α, β = 1, 2, · · · , N label electrons, I = 1, 2, · · · , M labels nuclei, and Z I is the atomic number of nuclei I. The first term in the Hamiltonian is the electron kinetic energy T , the second term is the Coulomb electron-nucleus attraction V ext and the third term is the electron-electron Coulomb repulsion V ee . The constant nuclei-nuclei term has been omitted. Even though the question is well-posed, obtaining the many-electron wavefunctions Ψ i ( r 1 , . . . , r N ) computationally is an extremely difficult task because of the exponential JHEP05(2016)046 to a discontinuity in the DFT exchange-correlation potential δE xc /δn( r) when electrons are added above the gap [116,117]. There are methods based on many-body perturbation theory to improve the DFT band gap and band shapes, such as the GW method [118]. However, since the largest contribution to the scattering rate comes from the low energy excitations, we choose to follow an empirical "scissor correction" approach [97,98]. In this approach a rigid shift is imposed on the conduction bands with respect to the valence bands in order to set the band gap to the experimental values of 1.11 eV for silicon and 0.67 eV for germanium [99]. It is worth noting that the semiconductor band gap features a temperature variation of around 10 meV [96], but we are choosing the room temperature band gap values for our calculation.

F.2 Energy density functionals
In order to be able to use DFT, a choice for the exchange-correlation functional E xc [n] is required. The Local Density Approximation (LDA) [94] has been remarkably successful because of its simplicity and transferability. In this method the exchange-correlation energy functional is based only on physical considerations and is approximated locally by the energy of a homogeneous electron gas (HEG) with the following density:

xc
[n] = d 3 r n( r) HEG x (n( r)) + HEG c (n( r)) . (F.5) The HEG exchange [119] and correlation [120] energy functionals are well established. There are some faults in the LDA which are known to be most dramatic where the electrons are highly localized and exchange repulsions are significant. In order to correct for that, the Generalized Gradient Approximations (GGA) introduce a dependence on the density gradient in the exchange-correlation energy density E GGA

xc
[n] = d 3 r n( r) GGA xc (n( r), |∇n( r)|). (F. 6) In this work we choose the well-established PBE functional [95] which is known to produce a broad set of properties of materials to accuracies of order a few percent [121]. Since LDA functionals tend to underestimate the energies of excited states compared to GGA functionals, we find a difference in cross-section sensitivity of around 10-20%, with a larger difference at higher thresholds.

F.3 Pseudopotentials
The valence electrons are responsible for the formation of interatomic bonds and their wavefunctions are in general delocalized, spanning over interatomic distances. The core electron wavefunctions, however, are very localized around the nucleus and they barely change from the isolated atom to the condensed matter phase. This fact allows to use the atomic core electron wavefunctions in the condensed matter phase by replacing the bare positive nuclear Coulomb potential and the negative Coulomb potential generated by the core wavefunctions with a pseudopotential in the Kohn-Sham problem eq. (F.2). The advantage is two-fold: first, the number of electrons in the problem is reduced to the JHEP05(2016)046 number of valence electrons and second, the only wavefunctions to be calculated are valence wavefunctions which, since they are rather smooth, do not require as fine a grid to represent them as a core electron wavefunction would, thus improving the computational efficiency. In this work we use Vanderbilt-type ultrasoft pseudopotentials [96]. The pseudopotential for Si includes the 3s and 3p electrons in the valence, while in the case of germanium, we use a pseudopotential which includes the 3d, 4s and 4p electrons in the valence.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.