Crystal responses to general dark matter-electron interactions

We develop a formalism to describe the scattering of dark matter (DM) particles by electrons bound in crystals for a general form of the underlying DM-electron interaction. Such a description is relevant for direct-detection experiments of DM particles lighter than a nucleon, which might be observed in operating DM experiments via electron excitations in semiconductor crystal detectors. Our formalism is based on an effective theory approach to general non-relativistic DM-electron interactions, including the anapole, and magnetic and electric dipole couplings, combined with crystal response functions defined in terms of electron wave function overlap integrals. Our main finding is that, for the usual simplification of the velocity integral, the rate of DM-induced electronic transitions in a semiconductor material depends on at most five independent crystal response functions, four of which were not known previously. We identify these crystal responses, and evaluate them using density functional theory for crystalline silicon and germanium, which are used in operating DM direct detection experiments. Our calculations allow us to set 90% confidence level limits on the strength of DM-electron interactions from data reported by the SENSEI and EDELWEISS experiments. The novel crystal response functions discovered in this work encode properties of crystalline solids that do not interact with conventional experimental probes, suggesting the use of the DM wind as a probe to reveal new kinds of hidden order in materials.


I. INTRODUCTION
An invisible and unidentified mass component, dark matter (DM), is the leading form of matter in galaxies, galaxy clusters and the large scale structures we observe in the Cosmos [1]. It is responsible for the bending of light emitted by distant luminous sources and provides the seeds for the gravitational collapse modern cosmology predicts to be at the origin of all objects we see in the Universe [2]. And yet, the nature of this essential, but elusive cosmological component remains unidentified.
In the leading paradigm of astroparticle physics, DM is made of hypothetical, as yet undetected particles that the Standard Model of particle physics cannot account for [3]. While different methods have been proposed in the past decades to unveil the identity of DM, only a direct detection of the hypothetical particles forming the Milky Way DM component will likely prove the microscopic nature of DM conclusively [4].
DM direct detection experiments play a central role in this context [5,6]. Typically, they operate ultralow background detectors located deep underground to search for signals of interactions between DM particles from our galaxy and nuclei forming the detector material [7]. So far, this experimental technique has mostly focused on the search for nuclear recoils induced by Weakly Interacting Massive Particles (WIMPs) -a class of DM candidates with interactions at about the scale of weak interactions and mass ranging from a few GeV to a few hundreds of TeV [8]. The upper bound arises from the requirement that the WIMP annihilation cross section is unitary [9], whereas the lower bound guarantees that the predicted density of WIMPs in the present Universe matches that measured via cosmic microwave background (CMB) observations [10]. The possibility of explaining the present DM cosmological density in terms of particle masses and coupling constants is one of the defining features of the WIMP paradigm and is based upon the so-called WIMP chemical decoupling mechanisms, which occurs when the rate of WIMP annihilations in the early Universe equals the rate of expansion of the Universe [11].
In spite of about four decades of searches at DM direct detection experiments, WIMPs have so far escaped detection [12]. While WIMPs remain a central element of modern cosmology, the fact that they have not been detected has recently motivated a critical reconsideration of the assumptions underlying the WIMP paradigm, and in particular the DM direct detection technique [13]. One important example is the restricted mass range within which the search for DM particles has so far been performed. Indeed, DM particles of mass smaller than about 1 GeV would be too light to induce an observable nuclear recoil, and this might explain why experiments searching for WIMPs via nuclear recoils have so far not been able to report an unambiguous discovery. On the other hand, a DM particle of mass in the MeV -GeV range might have enough kinetic energy to induce an observable electronic transition in a target material -a possibility which has recently attracted a great deal of attention [14].
From the theoretical side, the search for sub-GeV DM particles involves the modelling of DM-electron interactions in detector materials, e.g. [13,21,39]. One generic feature of models for DM-electron interactions is that they involve a new mediator particle, in addition to the DM candidate [14]. A new mediator is needed to reconcile the chemical decoupling mechanism with current observations on the present DM density. In this context, the most extensively investigated model extends the Standard Model by an additional U (1) gauge group, under which only the DM particle is charged [40]. The associated gauge boson is referred to as the "dark photon". Interactions between the DM particle and the electrically charged particles in the Standard Model arise from a "kinetic mixing" between ordinary and dark photons. While this framework allows to interpret the results of present and future DM direct detection experiments, it is also rather restrictive, as it a priori excludes the possibility that the amplitude for DM-electron scattering depends on momenta different from the transferred momentum [41]. Scenarios in which the amplitude for DM-electron scattering depends on the initial electron momentum, in addition to the transferred momentum, include models where the DM-electron interaction is generated by an anapole moment or a magnetic and electric dipole [42,43].
In this work, we extend the formalism that we developed in Ref. [41] for the scattering of DM particles by electrons without making any restrictive assumption on the form of the underlying DM-electron interaction, to the case of electrons bound in crystals used in operating DM direct detection experiments. Our formalism is based on an effective theory approach to nonrelativistic DM-electron interactions, and a set of crystal response functions defined in terms of electron wave function overlap integrals. Effective theory methods have previously been used in the scattering of DM particles by nuclei [44,45], in modelling collective excitations in DM direct detection experiments [46] and in a study of the DM scattering by bound electrons in isolated atomic systems [41]. This latter work introduced the notion of "atomic response" to DM-electron interactions that we here extend to the case of semiconductor crystals. By applying our new formalism to the study of DMelectron scattering in crystals, we discover that, under standard assumptions for the local DM velocity distribution, the rate of DM-induced electronic transitions in a semiconductor material depends on at most five independent "crystal response" functions. We express these response functions in terms of electron wave function overlap integrals and evaluate them numerically using QEdark-EFT [47], our extension of the QEdark code [21], which relies on the integrated suite of open-source computer codes, QuantumEspresso [48]. Leveraging on this finding, within our effective theory framework we are able to set 90% confidence level exclusion limits on the strength with which DM can couple to electrons, for general models. As illustrative examples we explicitly give the limits for specific DM models that were not tractable before our work, including the anapole, magnetic and electric DM-electron interaction models. Four of the five crystal response functions computed in this work were not known previously, and have explicitly been identified for the first time here. From a practical point of view, these can be used to compute the rate of DM-induced electron excitations in virtually all DM models where the free amplitude for DM-electron interactions does not explicitly depend on the mass of the particle that mediates the underlying interaction. By promoting the coupling constant to a general function of the momentum transfer our formalism could in principle be applied to any interaction model. At a more speculative level, the novel crystal response functions discovered in this work encode properties of crystals that have so far remained hidden, and that could be revealed if it becomes possible to use the DM particles that form our Milky Way as a probe in scattering experiments with semiconductor targets.
This work is organised as follows. In Sec. II, we introduce our formalism to model the scattering of DM particles by electrons in semiconductor crystals. In Sec. III, we introduce and numerically evaluate the novel crystal response functions that our formalism predicts. In Sec. IV, we present the 90% confidence level exclusion limits on the strength with which DM can couple to electrons, both within our effective theory framework and within specific DM models. We summarise and conclude in Sec. V. Finally, details underlying our analytical and numerical calculations are presented in the appendices.

II. DARK MATTER-INDUCED ELECTRONIC EXCITATIONS IN CRYSTALS
In this section, we extend the formalism of DMelectron scattering in crystals to general non-relativistic DM-electron interactions, including the anapole, magnetic and electric dipole couplings. We first introduce an expression for the rate of DM-induced electronic excitations that applies to arbitrary target materials and interactions, Sec. II A. Then, we narrow it down to the case of crystals, Sec. II B, and general non-relativistic DMelectron interactions, Sec. II C.

A. General rate of electronic excitations
For a generic target material and arbitrary DMelectron interactions, the rate of DM-induced electronic excitations from an initial electron bound state |e 1 to a final state |e 2 is given by [41] where the initial and final states, |e 1 , and |e 2 , are energy, not momentum eigenstates, and E i (E f ) is the initial (final) energy of the DM-electron system. In Eq. (1), the squared electron transition amplitude, |M 1→2 | 2 , is defined in terms of the initial and final momentum-space electron wave functions, ψ 1 and ψ 2 , and of the amplitude for DM scattering by free electrons, M, as where a bar denotes an average (sum) over initial (final) spin states and we integrate over the electron momenta, . Here, q = p − p is the momentum transferred, p is the outgoing DM particle momentum, and p = m χ v, where v is the initial DM particle velocity, while m χ and m e are the DM and electron mass, respectively. Eq. (1) also depends on the DM number density, n χ , and the local DM velocity distribution in the detector rest frame, f χ . For the latter, we assume where With this definition for N esc , f χ (v) is normalised to one. In all numerical applications, we set the most probable speed to the value of the local standard of rest, v 0 = 220 km sec −1 [49], the local galactic escape velocity to v esc = 544 km sec −1 [50], the speed of the Earth in the galactic reference frame (where the mean DM particle velocity is zero) to v ⊕ = 244 km sec −1 and the local DM number density to n χ = 0.4 GeV/cm 3 /m χ [51].

B. Rate of electronic excitations in crystals
In the case of crystalline materials, we label the initial (final) electron state by a band index i (i ) and a wavevector in the first Brillouin Zone (BZ) k (k ), i.e. in the notation of Eq. (1), 1 ≡ {ik} and 2 ≡ {i k }. Furthermore, we express the initial electron position space wave function at x in the Bloch form and similarly for ψ i k . Here, V = N cell V cell is the volume of the crystal, V cell is the volume of a single  [41,44,45]. Se (Sχ) is the electron (DM) spin, v ⊥ el = v − /me − q/(2µχe), where µχe is the DM-electron reduced mass, v ⊥ el is the transverse relative velocity and 1χe is the identity in the DM-electron spin space. (1) using wave functions ψ ik and ψ i k of the type in Eq. (5), we denote the corresponding rate of DM-induced electronic excitation by R ik→i k .
In the case of DM direct detection experiments using semiconducting crystals such as silicon and germanium as target materials, the observable quantity is the total rate of valence to conduction band electron excitations in the whole crystal, R crystal , which is given by The factor of 2 is a result of the spin degeneracy and consequent double occupation of each crystal orbital. In order to evaluate Eq. (6), we expand the free scattering amplitude, M, in the small electron momentum to electron mass ratio. As we will see in detail below, this nonrelativistic expansion allows us to identify the response of crystals to general DM-electron interactions in a modelindependent manner. Note that, in principle, one should obtain this total rate by summing Eq.
(1) over all unfilled conduction bands and all filled valence bands. However, in practical applications one has to truncate the number of conduction bands included in the sum to a manageable number, as we will see in Sec. III B.

C. Non-relativistic expansion
Assuming that both initial and final electron and DM particle move at a non-relativistic speed, the free scattering amplitude, M, can in general be expressed as a function of two momenta only [45]: the transferred momentum, q, and the transverse relative velocity v ⊥ el , i.e. the component of v that in the case of elastic DMelectron scattering is perpendicular to q. Namely, M = M(q, v ⊥ el ), where v ⊥ el = v − /m e − q/(2µ χe ) and µ χe is the DM-electron reduced mass. By expanding M in the electron momentum to electron mass ratio, | |/m e 1, we find While Eq. (7) applies to any model for DM-electron interactions as a first order expansion in | |/m e 1, it is an exact equation in the case of the so-called nonrelativistic effective theory of DM-electron interactions, where the free scattering amplitude is by construction expressed as a sum of interaction operators in the DM and electron spin space that are at most linear in v ⊥ el , where the interaction operators O i for spin 1/2 DM are defined in Tab. I, q ref is a reference momentum given by q ref ≡ αm e and α is the fine structure constant. Angle brackets in Eq. (8) denote matrix elements between the two-component spinors ξ λ and ξ s (ξ λ and ξ s ) associated with the initial (final) state electron and DM particle, respectively. For example, O 1 ≡ ξ s ξ s ξ λ ξ λ . Finally, the dimensionless coefficients c s i and c i in Eq. (8) are the coupling constants of the interaction operators in Tab. I. When c s i = 0 and c = 0, we refer to the interactions in Tab. I as of contact type; we refer to them as of long-range type when c s i = 0 and c = 0. In the non-relativistic limit, almost any model for DMelectron interactions in crystals can be matched onto the free scattering amplitude in Eq. (8). By substituting Eq. (8) into Eq. (2), we find where are electron wave function overlap integrals. The first one (Eq. (10)) arises within the standard treatment of DM-electron interactions in crystals [21]. The second one (Eq. (11)) is responsible for the novel crystal responses we define below in Sec. III. An explicit calculation performed in Appendix B shows that from each term in Eq. (9), we can collect a Dirac delta function and rewrite Eq. (9) as follows where ∆G ≡ G − G and Consequently, where |M ik→i k | 2 is defined through the comparison of Eq. (9) and Eq. (B9). In the non-relativistic limit, the initial and final DM-electron energies, E i and E f , respectively, are given by where v = |v| is the initial DM particle speed and E ik (E i k ) is the energy of the initial (final) electron bound state.
where q = |q| and θ is the angle between the momentum transfer q and the initial DM particle velocity v. By replacing Eq. (18), Eq. (15), Eq. (5) and Eq. (1) in Eq. (6), for the total rate we find We can now use the Dirac delta function in Eq. (19) to perform the integration over the polar angle θ in the velocity dependent part of the total excitation rate R crystal that we here denote by where We find where In the first step of Eq. (22) we follow Essig et al. [21] and use the simplification f χ (v) = f χ (v), performing the integration over θ by using the Dirac delta function. In the second step of Eq. (22), we introduce the azimuthalangle-averaged squared transition amplitude, defined as Again following [21], in the third step of Eq. (22) we restore the angular dependence of f χ (v) = f χ (v) and introduce the linear operator η(q, ∆E ik→i k ). Acting on Cg(v) where C is a constant and g(v) a function of the DM particle velocity in the detector rest frame, it gives In terms of η and M i,k→i ,k 2 θ=θ(q,∆E ik→i k ) , we can finally express the total excitation rate in a crystal as An interesting class of models in the search for DM via electronic excitations is that in which DM couples to electrons via higher order moments in the multipole expansion of the electromagnetic field [41,43]. If χ (ψ) is a Majorana (Dirac) spinor describing the DM particle, g a dimensionless coupling constant and Λ a mass scale, the anapole, magnetic dipole and electric dipole DM models are described by the following interaction Lagrangians [41], where F µν = ∂ µ A ν − ∂ ν A µ is the photon field strength tensor, and A ν the photon field. In the non-relativistic limit, the free electron scattering amplitudes associated with the Lagrangians in Eq. (29) are [41] M anapole = 4eg where g e 2 is the electron g-factor. From the free electron scattering amplitude in Eq. (30), we find that, in the case of anapole DM-electron interactions are the only coupling constants different from zero. From the amplitude in Eq. (31), we find that the only coupling constants different from zero in the case of magnetic dipole DM-electron interactions are Finally, from the amplitude in Eq. (32), we find that in the case of electric dipole DM-electron interactions, one coupling constant only is different from zero,

III. NOVEL CRYSTAL RESPONSES
We now focus on the novel crystal responses that arise from the electron wave function overlap integral in Eqs. (10) and (11).

A. Excitation rate and crystal response functions
As shown in the appendix, the azimuthal-angleaveraged squared transition amplitude, can be expressed as the sum of r products between a DM response function R l (q, v) and a crystal response function W l (q, ∆E), l = 1, . . . r, where ∆E is defined implicitly via Eq. (37) below. As a result, the total electron excitation rate in crystals can be written as where the DM response functions R l (q, v) given in App. C depend on the coupling constants c s i and c i , in addition to v, q and ξ, and Here, This factorisation into DM and crystal responses of the total electron excitation rate in crystals is analogous to the factorisation of the total DM-induced ionisation rate in atoms we found in [41]. In App. C 3, we show that within our simplified treatment of the velocity integral the two vectorial crystal responses, W 6 and W 7 , arising from the overlap integrals, are zero. For this reason, they are not included in Eq. (37). This is a good approximation for isotropic materials such as silicon and germanium, but will not suffice for anistropic materials such as graphene. For the rest of this paper we will focus on the 5 responses relevant to silicon and germanium. Because of the simplified treatment of the velocity integral discussed in the previous section [21], we can now perform the integral over the direction of the transfer momentum q, finding where We can then use the Dirac delta function in Eq. (37) to perform the integral over d 3 q , obtaining which is our final expression for the novel crystal responses to general DM-electron interactions that we implement in the following.
For our self-consistent calculations for silicon and germanium we use the Si.pbe-n-rrkjus psl.0.1.UPF and Ge.pbe-dn-rrkjus psl.0.2.2.UPF pseudopotential provided with the QuantumEspresso package, which include the 3s 2 , 3p 2 electrons for silicon and 4s 2 , 4p 2 and 3d 10 electrons for germanium in the valence configuration. The electron-electron exchange and correlation are treated using the PBE functional [55]. We sample reciprocal space using a 6×6×6 Monkhorst-Pack k-point grid, supplemented with additional k-points at and close to Γ and half way to the zone boundary to give a total of 243 k points; this grid was shown in Ref. [21] to give good convergence of the calculated scattering cross sections. We take an energy cutoff, E cut of 120 Ry (1.6 keV) for silicon and 100 Ry (1.4 keV) for germanium to expand the plane-wave basis, as explained below. We perform our calculations at the established experimental lattice constants of a Si = 10.3305 a.u. and a Ge = 10.8171 a.u.. With the exception of E cut , the above parameters are the same as those used previously in Ref. [21]. Since the spurious self-interaction within the PBE exchangecorrelation functional causes the Ge 3d bands to lie ∼ 5 eV higher than observed experimentally [56], we apply a Hubbard U correction with a value of U eff = 9.45 eV to the Ge 3d orbitals using the approach of Cococcioni and de Gironcoli [57]. This has the effect of shifting the narrow Ge 3d band rigidly down in energy by ∼ U eff 2 so that its position below the Fermi level (∼ 30 eV) is consistent with experimental observations. Our resulting densities of states (Fig. 1) show the usual DFT underestimation of the band gaps, which we correct in our response calculations using a scissors correction to set the band gap of silicon (germanium) to the experimental value of 1.2 eV (0.67 eV).
In order to numerically evaluate Eq. (47) we discretise it by introducing binning in q and ∆E, where q n and ∆E m are the central values of the n th qbin and m th ∆E-bin respectively. We use 2000 bins in q and ∆E, letting q take values between 0 and 40.8 keV (37.3 keV) for silicon (germanium), and ∆E take values between 0 and 85 eV. From Eq. (23) we see that the minimum velocity required to cause excitation with these highest values of ∆E and q is 625 km/s (684 km/s), well below the maximum dark matter velocity, v esc +v ⊕ . As in Ref. [21], we replace the k-integral with a discrete mesh and numerically evaluate where the k-sums go over the 243 k-point grid mentioned above, each with weight ω k , such that k ω k = 2. The lattice vectors ∆G satisfy the cutoff relation causing q to have a cutoff roughly at √ 2m e E cut , where E cut is our plane-wave energy cutoff described above. Our values of E cut = 1.6 keV (E cut = 1.4 keV) for silicon (germanium) correspond to cutoffs in q of ∼ 41 keV (∼ 37 keV). Note that these values are considerably larger than those used by [21], which is necessary for two reasons: First, some of the new interactions that we consider here, in which the DM response functions R l (q, v) depend on positive powers of q, means that the integrand in the crystal excitation rate formula, Eq. (36), is significantly different from zero for higher q-values than in the previously studied case of the dark photon model. And second, because we cover a larger range of deposited energies. The higher energy cutoffs require in turn inclusion of a larger number of bands in our DFT calculation; we include 102 unoccupied bands, which again is almost twice as many as in Ref. [21].
Note that, as in Ref. [21] our transition matrix elements are calculated between the Kohn-Sham pseudowavefunctions for the occupied and empty states. It will be an interesting future direction to evaluate the effect of the use of all-electron rather than pseudo wavefunctions on the calculated crystal excitation rates for our novel response functions, as was recently done in Ref. [58] for the previously known W 1 response.

IV. RESULTS
In this section, we numerically evaluate the five response functions in Eq. (37) focusing on silicon and germanium crystals. We then use these response functions to compute the expected crystal excitation rates for the anapole, magnetic dipole and electric dipole DM interactions, as well as for specific (linear combinations of) interaction operators in Tab. I. Finally, we apply our responses to set 90% exclusion limits on the strength of such interactions from the null result reported by the EDELWEISS [60] and SENSEI [59] experimental collaborations.
A. Novel response functions for silicon and germanium crystals We display our crystal response functions in the (q, ∆E) plane, focusing on silicon in Fig. 2 and on ger-  Fig. 3. In each panel of the two figures, the value of the corresponding response function is given by the color bar, with darker colors corresponding to higher values. As described in Appendix D, our first response function, W 1 , is equal to 8∆Eαm 2 e /q 3 times the crystal form factor introduced by Essig et al. in [21] (see our Eq. (D2)). Taking into account this q and ∆E dependent pre-factor, our first crystal response functions, W 1 , (upper left panel in Figs. 2 and 3) are in broad agreement with those given in Fig. 5 of Ref. [21]. There are, however, some important differences between our calculations and those in [21]. First, the corrected position of our Ge 3d levels means that the abrupt increase in W 1 as a function of ∆E occurs at a ∆E value of 30, rather than 25, eV. Second, the higher energy cutoffs that we used allow us to calculate our response functions over a larger energy range.

manium in
We now focus on the remaining four response functions, W 2 , W 3 , W 4 and W 5 , which we compute here for the first time. Since W 2 is a complex response function, we consider its real and imaginary part separately. Focusing on silicon crystals, the upper right panel of Fig. 2 shows (W 2 ), while the central left panel reports the corresponding imaginary part, (W 2 ). In the same figure, the central right panel shows W 3 , the lower left panel displays W 4 , while the lower right panel reports W 5 . Fig. 2 shows the analogous crystal response functions for germanium crystals. In all panels of Fig. 2 and Fig. 3, the points below the dashed black line are kinematically inaccessible for DM particles that are gravitationally bound to our galaxy, given our assumptions for the earth and local escape velocity. Note the different color-bar scales in each case. In particular, W 2 and W 5 have linear color bars since they can take both positive and negative values. We find that, W 2 and W 5 are orders of magnitude smaller than the other responses, and are therefore expected to give a sub-leading contribution to the total crystal excitation rate. In order to understand this result, it is useful to compare our findings with the results of [41] on DM-electron scattering in atoms.

B. Comparison to atomic responses
Each of the response functions identified here has in principle an atomic analogue. Specifically, we find that our crystal response functions can be mapped on to the atomic response functions found in Ref. [41] if one replaces the crystal Bloch wave functions in Eq. (5) with the wave functions of electrons in atoms, and the summation/integration over lattice/crystal momenta and the summation over band indices in Eq. (37) with a summation over the atomic azimuthal and magnetic quantum numbers (compare Eq. (37) with Eq. (41) in [41]). This comparison shows that the atomic response functions associated with W 2 and W 5 through this mapping are exactly zero. This is the result of A ≡ mm f 1→2 f 1→2 being real and (anti) parallel to q, where the sum is over the initial and final electron magnetic quantum numbers while "1" and "2" label the initial and final state of the atomic electron, see [41]. Indeed, in the case of atoms A is expected to be proportional to q, q being the only preferred direction in an otherwise spherically symmetric system. Specifically, the electron wave function overlap integral f 1→2 (q) (the analogue of f i,k→i ,k ) is axially symmetric around the direction of q. Similarly, f 1→2 (q), the analogue of f i,k→i ,k , is spherically symmetric. In contrast, the spherical symmetry of f i,k→i ,k and the axial symmetry of f i,k→i ,k are only approximate in the case of crystals. This explains why, although W 2 is approximately real and W 5 is subleading, for semiconductor crystals W 2 and W 5 are not exactly zero. We expect more dramatic departures from axial symmetry in the case of anisotropic materials, such as 2D materials, or anisotropic 3D Dirac materials. For these systems, we therefore expect larger values for the responses W 2 and W 5 .

C. Predicted electron excitation rates
Once the crystal responses have been calculated we can insert them in Eq. (45) to obtain the expected excitation rates. In Fig. 4 we show the expected excitation rates for the 14 non-relativistic operators of Tab. I for DM masses of 0.5 MeV, 5 MeV and 50 MeV for long-and short-range interactions. A notable difference between the predicted excitation rates for silicon and germanium is that silicon has zero excitations for m χ = 0.5 MeV, which is due to the band gap of silicon being too large for galaxy-bound DM of this mass to cause excitations. We also see that for most operators the excitation rate is maximum for m χ = 5 MeV showing that both silicon and germanium are well suited to probe DM masses at the MeV-scale.
Here, we focus on the predicted excitation rates for the SENSEI [59] and EDELWEISS [60] experiments, which employ targets made of silicon and germanium crystals, respectively. Indeed, such experiments do not measure The band gap for germanium is considerably lower than that for silicon, which allows germanium target experiments to probe lower masses, as we discussed above.
Having specified the values of E gap and ε used in our analysis, we can now calculate the expected crystal excitation rate corresponding to a given number of electronhole pairs Q = 1, 2, etc... and to a given target material by setting the boundaries of the ∆E integration in Eq. (36) to the values required by Eq. (51). In order to illustrate the contributions to the excitation rate from the individual response functions, we find it convenient to rewrite Eq. (45) as where We start by computing the expected excitation rates for the anapole, electric dipole and magnetic dipole interactions using the relation between couplings and interaction scale Λ that we derived in Eqs. (33) to (35). For these interactions, Fig. 5 shows the expected excitation rates as a function of Q for silicon and germanium crystals, and for m χ = 10 MeV. The black dashed line gives the total crystal excitation rate while the solid colored lines give contributions from the individual responses. The light blue, light green, dark orange and yellow lines correspond to R 1 , (R 2 ), R 3 and R 4 respectively. In order to illustrate the dependence of our results on the DM particle mass, in Fig. 6 we plot the total excitation rate for 1 MeV, 10 MeV and 100 MeV. Focusing on events producing few electron-hole pairs, i.e. Q ∼ 1, we see that the excitation rate is maximum at around 10 MeV, and falls off both for higher and for lower masses. We also see a clear rise in the excitation rate for germanium at Q = 10 for m χ = 100 MeV caused by the 3d electrons. Without the Hubbard-U correction implemented in our DFT calculations, this peak would be located at Q = 9 (see Appendix A). Importantly, we find that the crystal excitation rate is not necessarily dominated by the crystal response function W 1 (the crystal form factor of Essig et al. [21]). This is apparent in the case of, e.g. magnetic dipole and anapole DM-electron interactions, but it is also true for other combinations of nonrelativistic effective operators in Tab. I.
In Fig. 5, we can see the expected excitation rates for different interaction types plotted against the number of electron-hole pairs produced. As expected, we see that these rates drop significantly for larger ionization signals since those require larger energy transfers. We can see that different responses are dominant for different inter-action types. Namely that for the case of an anapole interaction, a combination of R 1 and R 2 dominates, for magnetic dipole interaction at low electron-hole-pair production R 3 dominates whereas R 1 dominates elsewhere, and that the electric dipole interaction is dominated by the R 1 interaction. Figs. 7 and 8 show the crystal excitation rate as a function of the number of excited electron-hole pairs Q for two selected non-relativistic effective operators with coupling constants c s 7 = 1 or c 7 = 1 and c s 15 = 1 or c 15 = 1, respectively, and with a DM mass of m χ = 10 MeV. In both figures, the rate for silicon (germanium) is given in the left (right) panels and the results for short-(long-) range interactions are given in the top (bottom) panels. In Fig. 7 all the nonzero response functions give positive contributions of a similar magnitude, whereas in Fig. 8 the contribution from R 4 almost exactly cancels with that from R 3 , leading to a total crystal excitation rate that is orders of magnitude smaller than the excitation rate produced by R 3 alone. In the case of contact and long-range interactions of type O 15 in germanium crystals, we also find a peak at Q = 10 in the R 3 and R 4 contributions to the total excitation rate (see Fig. 8). This peak is due to the DM-induced excitation of 3d electrons. Interestingly, the previously mentioned cancellation between R 3 and R 4 washes out this peak, so it is not present in the total excitation rate. Even in this case, however, the 3d electrons have an impact on the total crystal excitation rate: they slow down the decrease of the crystal excitation rate for Q above about 10 electron-hole pairs. In the case of germanium crystals, we find a similar effect also for longand short-range interactions of type O 7 (see Tab. I and Fig. 7).

D. Exclusion limits
Once the expected excitation rates have been computed we can compare them to the number of events measured experimentally to place limits on the mass and couplings of the DM particle. As anticipated, here we focus on the results reported by the SENSEI@MINOS [59] and EDELWEISS [60] experiments. The former operates silicon semiconductor detectors, the latter employs germanium targets. We compare our theoretical predictions with the experimental data by assuming that all events reported by the two experimental collaborations are caused by DM-electron interactions in the detector. By "event" we refer to the production of electron-hole pairs as the result of DM-induced electron excitations in the given semiconductor target. In order to compute the expected number of events associated with a given value of Q, we multiply the crystal excitation rate in Eq. (36) by the effective experimental exposure corresponding to that Q value. For the SENSEI@MINOS experiment, we restrict our analysis to Q = (1, 2, 3, 4), and we take the different effective exposures for the four Q values given by the experiment, namely: (1.38, 2.09, 9.03, 9.10) g-day. The observed number of events in each of the four Q bins is: (758, 5, 0, 0), as reported in Tab. 1 and explained at the end of page 4 in [59]. For EDELWEISS, we also consider Q = (1, 2, 3, 4), but assume the effective exposures (0.04, 0.22, 1, 1) × 33.4 g×58h, respectively. The observed number of events in each Q bin is in this case (5814, 44706, 2718, 227), as one can see by digitising Fig. 2 of [60].
For each experiment, we calculate 90% confidence level (C.L.) exclusion limits on the strength of a given DMelectron interaction by requiring that, in each of the four Q bins we consider, the cumulative distribution function of a Poisson probability density function of mean equal to the predicted number of events in that bin is at least equal to 0.1 if evaluated at the observed number of signal events. Fig. 9 shows our 90% C.L. exclusion limits on the strength of an anapole, magnetic and electric dipole DMelectron interaction as a function of the DM particle mass. For comparison, in the same figure we also dis-play the 90% C.L. exclusion limits on these coupling constants from the null results of XENON10 [16,61], XENON1T [18], and DarkSide-50 [15] as obtained in [41]. We find that SENSEI@MINOS provides the strongest constraints on DM masses below 5 MeV, with EDEL-WEISS taking over for sub-MeV DM masses.
Using the same procedure we also calculate the 90% C.L. exclusion limits on the strength of the interactions O 7 and O 15 that we used as benchmarks in the previous subsection. The resulting constraints on c s 7 and c 7 , as well as on c s 15 and c 15 are given in Fig. 10. The dashed (solid) lines correspond to short-(long-) range interactions. The dark green line gives the 90% C.L. exclusion limit from SENSEI@MINOS, while the orange line gives the 90% C.L. exclusion limit from EDELWEISS. Due to its larger exposure, SENSEI@MINOS generically produces the strongest constraints above about 1 MeV. Below this threshold, however, the lower band-gap of germanium implies that the strongest constraints arise from EDELWEISS.

V. SUMMARY
We performed a comprehensive and modelindependent study of interactions between DM particles of our galactic halo and electrons bound in crystals. By modelling the DM-electron interactions in crystals using an effective theory approach, we identified the most general amplitude for DM-electron scatterings and the general responses by the crystal to these interactions. Our effective approach allows predictions of scattering rates in crystals for virtually all DM models, it applies e.g. to the anapole, magnetic dipole and electric dipole DM-electron interaction models.
In particular our study focuses on short-and long-range DM-electron interactions with silicon and germanium crystals which are currently used as targets in DM direct detection experiments.
This led us to discover that there at most five ways a crystal can respond to an external probe (not necessarily a DM particle) in the limit of non-relativistic short-and long-range interactions of the most general type. We identified these five independent crystal responses and expressed them in terms of electron wave function overlap integrals. By performing state-of-the-art DFT calculations, we evaluated the five responses, focusing on silicon and germanium crystals, as these are used in operating DM direct detection experiments such as SENSEI and EDELWEISS.
We applied our crystal response functions to predict the rate of DM-induced electron excitations in silicon and germanium crystals. We performed this calculation for a set of 14 non-relativistic interaction operators (of shortand long-range type) and for specific models, such as the anapole, magnetic dipole and electric dipole models.
As a second important application of our novel crystal response functions, we computed the 90% C.L. exclusion limits on the strength with which DM can couple to electrons in crystals by comparing our predicted crystal excitations rates with data collected at the SENSEI@MINOS and EDELWEISS experiments. We performed this cal-  culation for the already mentioned set of non-relativistic DM-electron interactions, as well as for the anapole, magnetic dipole and electric dipole DM models. We compared these limits to constraints arising from different DM direct detection experiments and identified the range of masses where silicon or germanium detectors are expected to set the most stringent bounds on DM-electron interactions.
Within the field of DM direct detection, our novel crystal responses will enable the scientific community to perform predictions for DM particle models that were not tractable before. This is for example the case for the anapole and magnetic dipole interaction models, as they generate crystal responses that were not known previously. Furthermore, our crystal response functions will allow the community to calculate at which statistical significance a given DM-electron coupling is excluded by the null result of a DM experiment using germanium or silicon semiconductor detectors. Finally, they will allow us to interpret a discovery in a next generation DM direct detection experiment within a range of DM models that can not be covered by the current most widely used approach to data analysis, which is based on the use of a single crystal form factor.
On a more speculative level, our work paves the way for a new line of research lying at the interface of astroparticle and condensed matter physics: the study of yet hidden material properties that are encoded in the response of materials to external probes sharing the same interaction properties as the elusive galactic DM component.  [63], SciPy [64], Web-PlotDigitizer [65], Wolfram Mathematica [66]. The computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputer Centre (NSC) and the Chalmers Centre for Computational Science and Engineering (C3SE).

ACKNOWLEDGMENTS
Appendix A: Hubbard U correction for the germanium 3d-states In this appendix we illustrate the effect of the Hubbard U correction. In Fig. 11 we see how W 4 is affected by the Hubbard U correction. In particular, we see the high value region produced by the 3d electrons is shifted from 25 eV to 30 eV. In Fig. 12 we see the impact this has on the expected excitation rate for the anapole, electric dipole and magnetic dipole interactions. For Q-values of 9 and 10 the the excitation rate differs with orders of magnitude. In this appendix we cover the gap between Eqs. (9)-(11) and Eqs. (12)- (14). The electron wave-function in a crystal is given on Bloch form as We see that, in the absence of the correction, the lower binding energy, by ∼ 5 eV, of the 3d electrons leads to excitation rates that are orders of magnitude larger for 9 and 10 electron-hole pairs.
where V = N cell V cell is the volume of the crystal. Inserting this in Eqs. (10) and (11) gives and respectively, where u i,k ≡ u i (k + G) and u i ,k ≡ u i (k + G ). We now compute |f i,k→i ,k | 2 , f i,k→i ,k f * i,k→i ,k and (f i,k→i ,k · w) f * i,k→i ,k · w with w and w being arbitrary 3-vectors. For |f i,k→i ,k | 2 we find, where ∆G ≡ G − G and ∆F ≡ F − F. From the delta functions we see that the terms in the double sum is only non-zero when ∆G = ∆F. The sums over F and G are identical except for the labeling, i.e.
The calculation for f i,k→i ,k f * i,k→i ,k is analogous, and produces: where Finally, Inserting Eqs. (B5), (B6) and (B8) in Eq. (9), and setting w = (∇ M) =0 and w = (∇ M * ) =0 , we obtain our Eq. (15), (B9) 1. The first dark matter response The first dark matter response is produced by the first term in Eq. (12) and is The second term of Eq. (12) requires more treatment. In [41] we expanded the second term and found where A ≡ f i,k→i ,k f i,k→i ,k * is a complex 3-vector. In the case of atoms we found that A is both real and antiparallel to q. Neither of these simplifications apply to the crystal. In order to separate R(q, v, c i ) from W (q, ∆E), we rewrite the above equation by expressing it in terms of (A) and (A), so we use that From this we can define the second scalar overlap integral as In addition to this we also find two complex vectorial overlap integrals, which we can call For the scope of this paper we have treated the velocity distribution as isotropic, and we can use this to further simplify Eq. (C4). We start by decomposing v and A in a component parallel to q and a component perpendicular wheren ⊥ v andn ⊥ A are unit vectors in the plane perpendicular to q. Inserting this decomposition in the terms proportional to v ⊥ el · A in Eq. (24)  This result follows from our treatment of the velocity distribution as isotropic. Intuitively one can understand our treatment as not considering directionality. We expect directional effects such as a daily modulation to be small in silicon and germanium crystals justifying our treatment of the velocity distribution. In less isotropic materials such as 2D materials we expect these directional effects to be important and we will not in these cases be able to absorb B 6 and B 7 (fully) in the other responses.