Describing the Migdal effect with a bremsstrahlung-like process and many-body effects

Recent theoretical studies have suggested that the suddenly recoiled atom struck by dark matter (DM) particle is much more likely to excite or lose its electrons than expected. Such Migdal effect provides a new avenue for exploring the sub-GeV DM particles. There have been various attempts to describe the Migdal effect in liquids and semiconductor targets. In this paper we incorporate the treatment of the bremsstrahlung process and the electronic many-body effects to give a full description of the Migdal effect in bulk semiconductor targets diamond and silicon. Compared with the results obtained with the atom-centered localized Wannier functions (WFs) under the framework of the tight-binding (TB) approximation, the method proposed in this study yields much larger event rates in the low energy regime, due to a $\omega^{-4}$ scaling. We also find that the effect of the bremsstrahlung photon mediating the Coulomb interaction between recoiled ion and the electron-hole pair is equivalent to that of the exchange of a single phonon.


I. INTRODUCTION
Although there has been overwhelming evidence for the existence of the dark matter (DM) from astrophysics and cosmology, its particle nature still remains mysterious. In the past decades, tremendous efforts have been invested into the search of the weakly interacting massive particles (WIMPs), one of the most promising DM candidate, but conclusive evidence has not yet emerged. Owing to the continuous developments in detector technologies in recent years, the frontier of the dark matter direct detection has been pushed to the mass range below the GeV scale, where the traditional detection methods based on the nuclear scattering are expected to lose sensitivity. So more and more theorists and experimentalists have begun to shift to other alternative proposals based on new detection channels and materials, such as with semiconductors [1][2][3][4], Dirac materials [5][6][7], superconductors [8,9], superfluid helium [10][11][12], and via phonon excitations [13][14][15] and bremsstrahlung photons [16,17], as well as other proposals and analyses .
The Migdal effect has aroused wide interest recently because it is theoretically clarified that the suddenly struck nucleus can produce ionized electrons more easily than anticipated for an incident sub-GeV DM particle [39], so exploring the sub-GeV parameter space is possible for the present detection technologies. There has been numerous theoretical proposals [40][41][42][43][44] and exper-imental efforts [45][46][47] dedicated to detecting sub-GeV DM particles via the Migdal effect. Ref. [48] first investigated the Migdal effect in semiconductor targets by exploring the connection between the DM-electron scattering in semiconductors and Migdal processes in isolated atoms. In our previous study [42], we proposed to describe the Migdal effect in semiconductors under the framework of the tight-binding (TB) approximation, where a Galilean boost operator is imposed on the recoiled ion to account for the highly local impulsive effect brought by the incident DM particle, while the extensive nature of the electrons in solids is reflected in the hopping integrals.
Meanwhile, nontrivial collective behavior in condensed matter system has also attracted attention and has been considered as a possible origin of observed signal features in various experiments [49]. Ref. [49] put forward two mechanisms based on the plasmon production induced by DM particles to explain the existing signal lineshapes that are at odds with the standard interpretations of the DM-electron scattering. As a typical many-body phenomenon in solids, the plasmon excitation cannot be understood in terms of standard two-body scattering, or non-interacting single-particle states. Therefore, the many-body physics may shed new light on the interpretation of the DM-detector interactions.
On the other hand, Ref. [50] examined the above postulates by concretely calculating the bremsstrahlung emission of plasmon induced indirectly by a DM particle. While the many-body effect, namely, the plasmon resonance can be well described in the same way as in the electron energy loss spectroscopy (EELS) analysis, the effect of a fast charged electron transversing the material is replaced with an abruptly recoiling ion, and this part of physics can be summarized with the bremsstrahlunglike process, in the context of classical electrodynamics or quantum mechanics.
In this study we integrate above two aspects, i.e., the electronic many-body effect and the bremsstrahlung-like description of the recoiled ion, into a coherent method to describe the Migdal effect in solids. In this picture the suddenly recoiled ion excites a primary electronhole pair via the bremsstrahlung photon, while the many-body physics is encoded in the dielectric function ǫ (k, ω), which is responsible for the screening of the pure Coulomb interaction and the emergence of collective plasma oscillations. It is tempting to investigate the consequence of the many-body physics on the DM semiconductor detector. In order to take into account the local field effect, the calculation of ǫ (k, ω) is extended to the case in crystalline environment using the random phase approximation (RPA). Based on this method, we concretely estimate the Migdal excitation event rates in bulk diamond and silicon semiconductors, respectively. Moreover, considering there have existed other methods describing the Migdal effect in solids, such as the aforementioned TB approximation, thus we also make comparisons with the event rates obtained from the TB approximation, or the atom-centered Wannier functions (WFs) [42].
This paper is organized as follows. We begin Sec. II by giving the quantum mechanical formalism and relevant scheme for calculating the Migdal excitation event rates induced by DM particles. Based on this approach, we then specifically calculate the event rates for diamond and silicon semiconductors, respectively in Sec. III. We conclude and make some comments on the methodology in Sec. IV. Some details of the derivation in the main text are provided in the Appendices.

II. ELECTRONIC EXCITATION IN THE QUANTUM FIELD THEORY DESCRIPTION
We begin this section with a short review of the description of the Migdal effect in crystal structures in the context of the quantum theory. The Migdal effect in solids refers to the process in which the incident DM particle collides with a nucleus in crystal, and the recoiled nucleus excites an electron across the band gap from valence state |j to a conduction state |i . As has been pointed out in our previous study [42] that the boost argument for the Migdal effect in an isolated atom used in Ref. [39] no longer applies for the crystalline solids, because one can not regard the whole crys-tal target as a big nucleus: it is the struck nucleus in the crystal that is recoiling, but not the whole crystal. A description fixed to the crystal is preferred. Inspired by the treatment of the plasmon production in Ref. [50], in this study we discuss the Migdal effect with similar description in the hope that sudden acceleration effect of the nucleus will be well captured with a bremsstrahlung-like process. This process can also be understood as the twobody scattering between the DM particle and a nucleus, with the scattered DM particle and the nucleus, as well as the electron-hole pair as the final states, i.e., the process χ (p χ )+N → χ p ′ χ +N (p N ) + electron-hole pair, where χ stands for the DM particle, N (p N ) (N (p ′ N )) is the nucleus before (after) the collision. The relevant Feynman diagrams are presented in Fig. 1, where e (i) and h (j) represent the electronic conduction and valence states, respectively, from which the amplitude is read as V is the volume of the crystal, m N is the nucleus mass, ε ij = ε i −ε j is the energy difference between the conduction states {|i } and valence states {|j }, V χN (q) represents the DM-nucleus contact interaction V χN (x) in momentum space, which is connected to the DM-nucleus cross section with the following relation, where A is the atomic number of the target nucleus, σ χn and µ χn = m n m χ / (m n + m χ ) represent the cross section and the reduced mass of the DM-nucleon pair, respectively. V N e (k) is the ion-electron Coulomb interaction propagator where Z ion is the number of valence electrons, and α is the fine structure constant. In the soft limit where p N · k/m N ≪ ε ij and |k| ≪ |p N |, the nucleus propagator in Eq. (2.1) can be simplified as [50] 1 Figure 1. The diagrams of the process χ (pχ) + N → χ p ′ χ + N (pN ) + h (j) + e (i).
At this stage, one may wonder, whether the bremsstrahlung-like process (where an electron is excited via exchanging a photon) amounts to a complete description of the actual physical process, considering that in discussion of the Migdal effect in an isolated atom, the excitation process is described with a full evolution of the Hamiltonian. So if one can show that bremsstrahlung-like description is equivalent to the boost argument for the Migdal effect in an isolated atom, the bremsstrahlung-like narrative can be convincingly generalized to the crystalline solids. To this end, we apply above discussion to an isolated atom by substituting the valence (|j ) and conduction (|i ) states, as well as the valence charge Z ion , with initial (|α ), final (|β ) atomic states, and the atomic effective charge Z eff , respectively, and the amplitude in Fig. 1 can be read off as the following, 5) where v N = p N /m N represents the velocity of the recoiled atom, and the correspondence 1 On the other hand, as demonstrated in Ref. [51] that the excitation rate derived with the boost argument in Ref. [39] is connected to above bremsstrahlung-like approach through the relation where m e is the electron mass, and E α (E β ) is the energy of the initial (final) state. Thus the bremsstrahlunglike framework can also describe the Migdal effect in an isolated atom. At first sight, one may conclude that this relation also holds in solids, and hence the atomic formulation can be directly generalized to the semiconductors [44,48]. However, this is not true because the derivation of Eq. (2.6) rests on a central Coulomb potential −Z eff e 2 /r, while in the crystalline environment the itinerant electrons are also subjected to the Coulomb interaction from all neighboring ions, and thus such generalization is invalid. In contrast, our rationale is that, since the bremsstrahlung-like narrative has been justified in the case of isolated atoms, we generalize this de-scription to the case in semiconductors, which leads to Eq. (2.1), rather than i|e imevN ·x |j .
Moreover, if one takes into consideration the screening effect by introducing the dielectric function ǫ −1 (k, ω) [49], the total event rate for the DM flux impinging on a semiconductor and then exciting a primary electron-hole pair can be written as * : where the bracket · · · denotes the average over the DM velocity distribution, ρ χ represents the DM local density, N T is number of the nuclei in target, f χ (v) is the DM velocity distribution, v ion = p N /m N is the velocity of the recoiled nucleus, ω denotes the relevant energy difference, and Θ is the Heaviside step function; the factor 1/3 stems from the average over all directions of k, given the isotropic nature of the electron gas; in the last line the factor 2 counts the two degenerate spin states. For the given p N and ω, function v min determines the minimum kinetically accessible velocity for the transition: being the reduced mass of the DM-nucleus pair.
In practice, we take ρ χ = 0.3 GeV/cm 3 , and the velocity distribution can be approximated as a truncated Maxwellian form in the Galactic rest frame, i.e., , with the Earth's velocity v e = 230 km/s, the dispersion velocity v 0 = 220 km/s and the Galactic escape velocity v esc = 544 km/s.

B. Migdal effect with random phase approximation
As illustrated in general expression of the Migdal event rate Eq. (2.7), the dielectric function ǫ −1 (k, ω) plays the central role in the estimate, accounting for the manybody effects. In this paper we calculate the dielectric function with the random phase approximation (RPA), which corresponds to the Lindhard formula: with which the total event rate including the electronic many-body effect is recast as follows This formula applies for the textbook model for the homogeneous electron gas. In crystalline structure the translational symmetry for spacetime reduces to that for the periodic crystal lattice. In this case, above expression is rewritten as the following † : where the nondimensional factor represents the averaged energy loss function, with Ω being the volume of the unit cell. The transferred momentum k in Eq. (2.10) is expressed uniquely as the sum of a reciprocal lattice vector G, and corresponding reduced momentum [k] confined in the first BZ, i.e., k = [k] + G. F (ω) can be obtained from inverting the microscopic dielectric matrix where n i (n j ) denotes the occupation number of the state |i (|j ). The non-vanishing G-vectors in the microscopic dielectric matrix reflect the variation of the microscopic field over a unit cell. For simplicity, the crystal structure is still approximated as isotropic in discussion. If one presumes that the fine structure constant α well suppresses the matrix elements, and then takes only the terms up to the first order in the resolvent where the identity matrix I and M represent the first and the the second terms in Eq. (2.13), respectively, the inverse matrix can be approximated as and hence which exactly corresponds to the case without the screening effect. In this study the inverse matrix is calculated in a straightforward way. † For more details of the calculation, see Appendix A 2 c

III. COMPUTATIONAL DETAILS AND RESULTS
Now we are in a position to put into practice the estimate of the Migdal excitation event rate. With Quantum Espresso package [52] plus a norm-conserving pseudopotential [53], we perform the density functional theory (DFT) calculation to obtain the Bloch eigenfunctions and eigenvalues using the local-density approximation [54] for the exchange-correlation functional, on a uniform 6×6×6 (5×5×5) k-point mesh for diamond (silicon) via the Monkhorst-Pack [55] scheme. A core cutoff radius of 1.3 Bohr (1.8 Bohr) is adopted and the outermost four electrons are treated as valence for diamond (silicon). The energy cut ε cut is set to 200 Ry (70 Ry) and lattice constant 3.577 Å (5.429 Å) for diamond (silicon) obtained from experimental data is adopted. The band structure of diamond (silicon) crystal is presented in the upper left (right) panel of Fig. 2, where the scissor correction is conducted to match the experimental value of the band gap E g = 5.47 eV (1.12 eV).  The matrix ǫ −1 G,G ′ is calculated via directly inverting the matrix Eq. (2.13) with the YAMBO code [56,57], with a smaller matrix cutoff E cut of 50 Ry (20 Ry) for diamond (silicon). An energy bin width ∆ω = 0.05 eV is adopted within the range from 0 to 50 eV. In the YAMBO implementation, the output of the Bloch wavefunctions are formatted in the form of periodic wavefunctions {u ik (x)}, normalized within a unit cell, with which the matrix element in Eq. (2.13) is explicitly written as where the integral is performed over the unit cell. This term can be understood as the Fourier transformation of the squared term u * i ′ k ′ (x) u ik (x) within the unit cell, and hence the calculation is practically performed using the discrete fast Fourier transformation (FFT) technique. So the resolution in the position space is associated with the truncation radius of the reciprocal Gvectors in Eq. (2.12), G, which is determined from the energy cut ε cut through the relation G = √ 2m e ε cut .
In practical evaluation of dielectric matrix Eq. (2.13), a small broadening parameter η = 0.1 eV is adopted for both diamond and silicon, instead of an infinitesimal energy width 0 + . Theoretically, the smaller the parameter η, the more accurate the computation is, but on the other hand, a smaller η in turn requires a finer energy width ∆ω and a denser k-point mesh to smear the spectra. Thus such choice of parameter η is the result of a balance between accuracy, efficiency, and smoothness. The integrals of the continuous k-points in the first BZ are replaced by the summations over a uniform discrete mesh of representative k-points as follows: 2) with N k being the number of k-points sampled in the first BZ. As mentioned above, a homogeneous set of 6×6×6 (5×5×5) k-points for diamond (silicon) is used in this study. Given these k-point meshes we compute the averaged energy loss functions F (ω) defined in Eq. (2.12) for various choices of energy cuts for diamond (left) and silicon (right) in the middle row of Fig. 2, while in the third row of Fig. 2 we show the calculated spectra F (ω) for various k-point meshes for the given energy cuts. Following the convention in computational condensed matter physics, we quantify the uncertainties in our computation in terms of the tendency of convergence depending on these k-points and cut-energies E cut adopted in calculation in Fig. 2. In our computation, the differences between the event rates (after integrating over 0 eV to 50 eV) calculated from two sets of parameters are found to be convergent within 5% for both diamond and silicon. Therefore the computational parameters for the k-point meshes (6 × 6 × 6 for diamond and 5 × 5 × 5 for silicon) and the energy cuts (50 Ry for diamond and 20 Ry for silicon) are sufficient to give a quantitative description of the Migdal effect at the present stage.
In the upper panel of Fig. 3 shown are the velocityaveraged energy spectra of the Migdal excitation for DM mass m χ = 10 MeV (blue), 100 MeV(orange) and 1 GeV (green), respectively, for a benchmark cross section σ χn = 10 −38 cm 2 . Due to a long tail brought by the non-vanishing η, which is further amplified by the factor ∝ ω −4 , the spectra do not exactly terminate at the band gaps, so we truncate the spectra at values slightly higher than the lower edges of their conduction bands. Also due to the dependence of the ∝ ω −4 factor, the spectrum features a peak roughly at 13 eV (5 eV) for diamond (silicon), and suffers a suppression at the high energy end, exhibiting a drastic variation in the whole energy range. The energy spectra are calculated up to ω = 50 eV, a value falls short to span all the energy range relevant for ionization signal for DM masses larger than 100 MeV, considering that the relevant spectra extend out beyond the energy window. However, in this energy window, one can still determine or constrain the DM particle parameters from the single-and two-electrons bins via the Migdal scattering. To fulfill this purpose, we adopt the model [3] where the extra electron-hole pairs triggered by the primary pair are described with the mean energy per electron-hole pair ε in high energy recoils. In this picture, the ionization charge Q is then given by where ⌊x⌋ rounds x down to the nearest integer. Thus, from the energy spectra we estimate the sensitivity of a 1 kg-yr detector in the bottom panel of Fig. 3, assuming an average energy ε = 13 eV (3.6 eV) for producing one electron-hole pair for diamond [58] (silicon [3]). The 90% C.L. upper limits on DM-nucleon cross section for both a single-electron (blue) and a two-electron (orange) bins are presented with no background event assumed.
It should also be noted that our calculation of electronhole production rate also includes the contributions from the bremsstrahlung plasmon production, and we assume that the plasmons decay dominantly to the electron-hole pairs [50].

IV. SUMMARY AND CONCLUSIONS
In this paper we make an alternative attempt to describe the Migdal effect in semiconductor targets, in which the description of the Migdal effect separates into two parts. Firstly, we incorporate the bremsstrahlunglike process to account for the drag force exerted on the electrons by the suddenly struck ion, where the Coulomb interaction is mediated by the bremsstrahlung photons. Secondly, many-body effects such as screening, collective behavior are also taken into consideration in the calculation of excitation rate of the electron-hole pairs. Based on this method, the Migdal excitation event rates for diamond and silicon semiconductors are calculated.
As mentioned in previous section, the spectra are modulated by the factor ω −4 , which results in a significant boost to the excitation event rate in the low energy regime and a suppression towards the high energy end. Such phenomenon is especially remarkable for a target possessing a small band gap, such as silicon presented in right panel of Fig. 3. It is interesting to make a comparison between the event rates calculated in this work and the ones obtained with the localized WFs in Ref. [42]. For this purpose we present the two spectra for diamond crystal in Fig. 4. Compared to the spectra calculated with the localized WFs in Ref. [42] (in dashed), the event rates calculated in this study (in solid) are found to be significantly larger in the low energy range near the band gap, and turn moderately higher towards the high energy end. Apparently, the ω −4 scaling behavior has not been reflected in the crystal form factor F (q, E e ) calculated in Ref. [42].
Finally, we try to interpret this method in the context of the quantization of the crystal vibration. We begin with the diagram in Fig. 5 that represents a process where p χ p ′ χ |j |i · · · 1 2 n − 1 n k Figure 5. The diagram for the multiphonon process where the DM particle excites an electron-hole pair via the exchange of a single phonon, while most of the transferred momentum of the DM particle is taken by a bunch of n phonons. See text for details. a pair of electron and hole is excited via the exchange of single phonon, along with a bunch of phonons created by the scattering of the DM particle and nucleus. It is noted that the phonon propagator can be approximated as iD k,α = i/ ε 2 ij − ω 2 k,α ≃ i/ε 2 ij , because the phonon eigenfrequncy of branch α at transferred momentum k, ω k,α ∼ O 10 −1 eV , is much smaller than the band gap. As a result of this approximation, the sum over the products of the two vertices relevant for the excitation can be contracted as the following, where ǫ k,α is the phonon eigenvector of momentum k. Taking into consideration also the interactions V χN and V N e from the two vertices, respectively, as well as other necessary factors, one observes that the processes in Fig. 1 and Fig. 5 are equivalent, except for a difference in interpreting the recoil effects in the crystal: in Fig. 1 the final state is the recoiled nucleus, whereas in Fig. 5 the final states are a large number of n phonons. In fact, it is implied through an isotropic harmonic oscillator toy model that in the limit q → ∞, the effects of the multiphonon final states can be well summarized with a recoiled nucleus, which is namely the impulse approximation [59]. In this regard, the two descriptions are essentially identical: the bremsstrahlung-like process in the soft limit is equivalent to a multiphonon process in the impulse approximation, where a valence electron is excited across the band gap via the exchange of a single phonon. However, in the low momentum transfer regime (or equivalently, in the low DM mass regime), the impulse approximation may no longer be valid, and the recoiling effect should be accounted for with a multiphonon process. A detailed analysis in Ref. [51] indicates that for m χ ≤ 50 MeV the impulse approximation ceases to be reliable for silicon and germanium targets, where phonons, rather than a free nucleus, are more appropriate for the description of the recoiling effect.
Note added In the final phase of preparing this paper, which was enlightened by Ref. [50], Ref. [51] appeared, which calculated the Migdal effect in semiconductor targets silicon and germanium with the similar techniques, and provided heuristic explanation for the treatment of the bremsstrahlung-like process, while in this work we calculate the event rates for diamond and silicon. Ref. [51] also pointed out the importance of the contribution of the off-diagonal terms of the dielectric matrix to the total Migdal event rate, which motivates us to further take into account the off-diagonal terms for a complete calculation of the Migdal event rates in diamond and silicon targets. In this Appendix, we provide some detailed derivation of the formulas in the main text.

Feynman rules
In order to describe the scattering processes in solids it is necessary to derive the Feynman rules at zero and finite temperatures. Here a brief review is arranged. We begin with the propagators in the position space.
• The DM-nucleus propagator between space-time coordinates x and x ′ is in which one neglects the retardation of the interaction in the non-relativistic limit.
• The nucleus propagator between x and x ′ is with ε p = |p| 2 / (2 m N ).
• The ion-electron Coulomb potential propagator between x and x ′ is In addition, the propagator of electrons in the crystal is similar to that of the nucleus. It is not shown here because it does not appear in the Migdal scattering process at the tree level.
On the other hand, we summarize the Feynman rules for the external legs as the following.
• The incoming and outgoing states of nucleus at x are represented with e ipN ·x−iεp N tx / √ V and x / √ V , respectively, with their corresponding energies ε pN = |p N | 2 / (2 m N ) and ε p ′ • The incoming and outgoing states of electron in crystal at x are represented with ψ j (x) e −iεj tx and ψ * i (x) e iεitx , respectively, with their corresponding energies ε j and ε i .
It is straightforward to translate above Feynman rules in the position space into those in the momentum space after integrating out the position 4-vector coordinates at each vertex. They are summarized as follows.
• Except for the case concerning the external legs of electrons in solids, each vertex contributes a factor (2π) V and an energy-momentum conservation condition presented as discrete delta functions δ ( i pi),0 δ ( ε i ).
• Vertex that contains both the incoming and outgoing states (|j and |i ) of the electrons in solids contributes a factor (2π) i|e ip·x |j and the energy conservation condition δ ( ε i ), where p is the net momentum that sinks into the vertex.
• Each DM particle or each nucleus external leg contributes a factor 1/ √ V .
• Each internal line corresponds to the propagator of its kind, as well as a factor 1/V . For example, the nucleus internal line is read as i/ (ω − ε p + i0 + ) multiplied by a factor 1/V . Applying these rules, and summing up all the discrete momenta at each vertex, the T -matrix can be expressed in terms of the amplitude as the following, Above usage of Feynman diagrams and rules in scattering theory can be transplanted to discussions in the linear response theory in a parallel fashion, where the imaginary time t → −iτ is introduced to account for the statistical field theory at a finite temperature. For example, a fermion propagator in the position space is written as where τ and τ ′ are the "temporal" coordinates, β = 1/kT is the inverse temperature, with Boltzmann constant k, ω n = n + 1 2 (2π/β) is the Matsubara fermion frequency, with n being an integer to ensure the anti-periodicity, and u i (x) is the normalized ith eigen wavefunction, with its eigen energy ε i , which can be free or bound state. Base on these propagators, Feynman rules for finite temperature filed theory can also be summarized and applied to calculation of physical quantities such as dielectric functions, and polarizabilities in the context of linear response theory.

Dielectric function and random phase approximation (RPA)
a. dielectric function and polarizability In a system where there is an external electromagnetic perturbation φ ext , the charge is redistributed and gets polarized. In order to describe the resulting total potential φ tot we invoke the inverse dielectric function ǫ −1 in the following way [60], where V Cou (x, x ′′ ) is the instantaneous Coulomb interaction, and the polarizability χ r ρρ is the density-density correlation function in the context of linear response theory, with · · · denoting the thermal equilibrium average. In a system with translation-invariance, the polarizability only depends on the differences of the space-time coordinates, i.e., χ r (x, t; x ′ , t ′ ) = χ r (x − x ′ , t − t ′ ; 0, 0), so it is convenient to discuss relevant problems in the momentum-energy space, where the inverse dielectric function can be expressed as a product with V Cou (q) being the reciprocal counterpart of the Coulomb potential in momentum space. Therefore, the calculation of the polarizability plays the central role in our investigation of the response of the solids induced by electromagnetic perturbation, and the DM particle. While the momentum component of χ r ρρ (q, ω) can be obtained by a direct Fourier transformation of Eq. (A7), deriving the frequency part at a finite temperature needs to be carried out with imaginary time Green's functions, or the so-called Matsubara Green's functions. The procedure for the density-density correlation is briefly summarized as follows. First, one introduces the Matsubara Green's function, which is defined as whereT τ is the imaginary time-ordering operator, β > τ > 0, and ν n = n (2π/β) is bosonic frequency sinceρ is a bosonic operator, with n being an integer to ensure the periodicity. Then one obtains χ r ρρ (iν n ) with the Fourier = + + + · · · transformation With the Lehmann representation, it is observed that χ r ρρ (q, ω) and χ ρρ (q, iν n ) are just special case of the same function defined in the entire complex plane except for a series of poles lying along the real axis. Thus, once χ ρρ (q, iν n ) is obtained from Eq. (A10), the retarded polarizability can be derived by performing the analytic continuation χ r ρρ (q, ω) = χ ρρ (q, iν n → ω + i0 + ).

b. RPA
The linear response theory can be discussed with our familiar language of path integrals, as well as necessary modifications accounting for the imaginary-time argument. The dielectric properties can be described in terms of diagrams in the top panel of Fig. 6, where the Schwinger-Dyson equation for the screened Coulomb interaction is represented by a bare interaction wiggly line attached to a geometric series over polarization bubbles. If we use W and Π to represent the screened Coulomb potential, and the self-energy corresponding to the 1-particle-irreducible (1PI) blob, respectively, the infinite sum is read as and hence one has On the other hand, adopting the random phase approximation (RPA) means that, only the electron-hole bubble is retained among all 1PI diagrams in calculation of the dielectric functions, which is illustrated in the bottom panel of Fig. 6. Thus within the framework of RPA, starting from Eq. (A9) and Eq. (A10), one first obtains the 1PI blob where n i (n j ) denotes the occupation number of the state |i (|j ), and then performing the analytic extension iν n → ω + i0 + and inserting the polarizability into Eq. (A8), one finally arrives at the Lindhard dielectric function and its inverse In above discussion we assume that the electronic system possesses a translational symmetry, which is true for the case such as homogeneous electron gas (HEG). However, for the case of crystal structure where the translational symmetry for continuous space reduces to that for the crystal lattice, the polarizability can no longer be expressed as differences of the space-time coordinates. In this case, the double periodic position function such as χ (x, x ′ ; ω) can be expressed in the reciprocal space as the following, where χ G,G ′ (q; ω) is the reciprocal matrix with G and G ′ being reciprocal lattice vectors and q is restricted to the first BZ, which can be determined with the Fourier transformation As a consequence, for an arbitrary momentum transfer Q, which can be split into a reduced momentum confined in the 1BZ, and a reciprocal one, i.e., Q = q + G, above Lindhard formula relevant for the excitation is generalized to the following microscopic dielectric matrix: where V Cou G,G ′ (q) = V Cou (q + G) δ G,G ′ = 4πα δ G,G ′ / |q + G| 2 is obtained from Eq. (A17), and the inverse dielectric function can be determined via matrix inversion.

c. Cross section and event rate
Here we first derive a general formula dedicated to the calculation of two-body scattering cross section for the case of discrete momenta. By modifying the derivation in textbook, we obtain the cross section between particles A and B, where v A − v B is the relative velocity between the two particles, p f is the momentum of the f th outgoing particle, and M (k A , k B → {p f }) is the amplitude in the discrete momentum space. Following the aforementioned Feynman rules, one can read and approximate the T -matrix from Fig. 1 as follows p ′ χ p N ; electron hole pair|iT |p χ 0 N = i where in the first line we separate out four-momentum conservation condition from other terms dependent on transferred momentum k of the Coulomb potential, and in the second line we assume that k is much softer than that of the recoiling nucleus p N , which is consistent with the soft limit condition and explains the summation over k in Eq. (2.1). Thus, one obtains the total cross section of an incident DM particle exciting an electron across the Fermi surface for the HEG via a recoiling nucleus, where v is the velocity of the DM particle in the frame of laboratory, v ion = p N /m N is the velocity of the recoiled nucleus, and −Im ǫ −1 (k, ω) is called the loss function. Since the amplitude i|e ik·x |j uniquely pins down the momentum k (through the momentum conservation), only the diagonal terms of momenta k's are relevant in above calculation for the HEG. Following Refs. [49,50], we adopt a Coulomb potential between the recoiled ion and electrons screened by dielectric function ǫ (k, ω) so as to obtain a general formula that includes self-interaction effects of the Coulomb propagator. To this end, we take the correspondence and use the relation |ǫ (k, ω)| 2 = Im −1 ǫ (k, ω) .
Thus, from Eq. (A21) one gives the Migdal exciting event rate where v min is defined in Eq. (2.8). In contrast to the case in HEG, in crystalline environment both the diagonal and the off-diagonal terms contribute to the total event rate ‡ , because in this case i|e ik·x |j only constrains momentum k up to a reciprocal lattice vector G, and thus the off-diagonal terms (G, G ′ ) survive when generalizing the above ‡ See also Ref. [51].