Phonon-mediated Migdal effect in semiconductor detectors

The Migdal effect inside detectors provides a new possibility of probing the sub-GeV dark matter (DM) particles. While there has been well-established methods treating the Migdal effect in isolated atoms, a coherent and complete description of the valence electrons in semiconductor is still absent. The bremstrahlung-like approach is a promising attempt, but it turns invalid for DM masses below a few tens of MeV. In this paper, we lay out a framework where phonon is chosen as an effective degree of freedom to describe the Migdal effect in semiconductors. In this picture, a valence electron is excited to the conduction state via exchange of a virtual phonon, accompanied by a multi-phonon process triggered by an incident DM particle. Under the incoherent approximation, it turns out that this approach can effectively push the sensitivities of the semiconductor targets further down to the MeV DM mass region.

The Migdal effect has attracted wide interest recently because the study in Ref. [43] has shown that in theory the suddenly struck nucleus can produce ionized electrons more easily than anticipated for an incident sub-GeV DM particle, so exploring relevant parameter region is plausible for the present detection technologies. Although the Migdal effect has not been directly observed in a nuclear collision, attempts to make the first such measurement from neutron-nucleus scattering are underway [44][45][46]. After Ref. [43], there has emerged numerous theoretical proposals [18,43,[47][48][49][50][51][52][53][54][55][56] and experimental efforts dedicated to detecting the sub-GeV DM particles via the Migdal effect in liquids [57], and in condensed matter targets [58][59][60][61].
Compared with the typical ionization energy thresholds in atoms ε g ∼ O (10) eV, semiconductor targets have a much lower thresholds ε g ∼ O (1) eV, which makes them ideal materials for further exploiting the the Migdal effect in the probe of light DM particles. However, generalizing the boosting argument in isolated atoms proposed in Ref. [43] to the crystalline environments faces both conceptual and technical obstacles: while one keeps pace with the recoiling nucleus, the ion lattice background will move in opposite direction, which brings no substantial convenience in mitigating the original complexity. Thus the semiconductor target at rest is still a preferred frame of reference. In Ref. [49], we made a tentative effort to describe the Migdal effect in semiconductors using the tightbinding approximation, where a Galilean boost operator is imposed specifically onto the recoiled ion to account for the highly local impulsive effect caused by the collision with an incident DM particle, while the extensive nature of the electrons in solids is reflected in the hopping integrals. Refs. [50,51] managed to describe the Migdal effect in solids in a manner analogous to bremsstrahlung calculation, where the valence electron is excited to the conduction state via the bremsstrahlung photons emitted by the recoiling ion.
The bremsstrahlung-like approach is an effective description of the Migdal event rates for DM masses m χ ≥ 50 MeV [50]. However, below this mass, the picture of a recoiling ion in the solid begins to break down and the effects of phonons become important. In Refs. [51] we proposed that the Migdal effect in solids can alternatively be described by treating the phonon as the mediator for the Coulomb interaction in the lattice between the abruptly recoiling ion and itinerant electrons. Thus the objective of this work is to provide a complete and self-contained theoretical foundation for this idea. Within this framework, numerous phonons, rather than an on-shell ion, are produced from the DM-nucleus scattering, especially in the low energy regime, where the scattering is coherent over the whole crystal. In the large momentum transfer limit however, the recoiling on-shell ion is expected to reappear as a wave packet supported by a large number of phonons. Such an asymptotic behavior should self-consistently justify the impulse approximation adopted in the bremsstrahlung-like approach. While the multi-phonon process has been thoroughly discussed in literatures (e.g., Ref. [62] and references therein, and see Refs. [50,[63][64][65] for recent discussions related to Migdal effect and DM searches), the fresh idea in this paper is to incorporate the generation of phonons, and the excitation of the electron-hole pairs, as well as the medium effect in solids, into a common framework. By doing so, it is no longer necessary to match the bremsstrahlung-like calculation onto the phonon regime, and the inherent conflict between the picture of a recoiling ion and that of the scattered phonons can be resolved altogether.
For convenience, our discussions are carried out by using the machinery of the quantum field theory (QFT), a language more familiar to the particle physics community. This approach proves intuitive and effective. As an interesting example, we derive the Debye-Waller factor with the Feynman diagram method, circumventing the awkward techniques associated with the operator commutator algebra (see Appendix A 3). Based on the calculated Migdal excitation event rates using this phonon-mediated description, we are able to push the sensitivities of the semiconductor detectors down to the MeV DM mass range. This paper is organized as follows. We begin Sec. 2 by giving the QFT framework for the multi-phonon process induced by DM particles. Based on this discussion, we then generalize the formalism to the Migdal excitation process in Sec. 3. We conclude and make some comments on the methodology in Sec. 4. A short review on the electrons and phonons in the context of the QFT, as well as other supporting materials are provided in Appendix A.

MULTI-PHONON PROCESS
In this section we first derive the formula for the scattering cross section between a DM particle and the target material, and then discuss the asymptotic behavior of the phonon spectrum towards the large momentum transfer limit. For simplicity, here we only consider the case of the monatomic simple crystal at 0 K. Figure 1. The diagram of the process χ (pχ) + target → χ p ′ χ + target + (k1, α1) + (k2, α2) + · · · (kn, αn). See text for details.
We consider the scattering process where n phonons {k j , α j } , (j = 1, 2, · · · , n) are generated by an incident DM particle in the context of the QFT, where {k j } and {α j } represent the phonon wavevectors in the first Brillouin zone (1BZ), and phonon polarization branches, of the final states, respectively. The relevant diagram is shown in Fig. 1, where the initial (p χ ) and final (p ′ χ ) DM states are replaced with an external source. With such replacement it is convenient to switch the scattering theory at zero-temperature to the linear response theory at a finite temperature, where the interest is focused on the response of target material to external perturbations. A more complete treatment of the composite lattice at a finite temperature lies beyond the scope of this work, and will be pursued in further investigation. Using the Feynman rules summarized in Appendix. A 5, the amplitude is read as where q = p ′ χ −p χ is the momentum transferred to the DM particle, G's are reciprocal lattice vectors, N is the number of the unit cells in the crystal, which equals the number of the atoms in a monatomic simple crystal, V is the volume of the material, m N is the nucleus mass, ǫ kj,αj and ω kj,αj are the phonon eigenvector and the eigenfrequency of branch α j at wavevector k j , respectively; V χN (q) represents the DM-nucleus contact interaction V χN (x) in momentum space, which connects to the DM-nucleon cross section σ χn through with A being the atomic number of the target nucleus, and µ χn = m n m χ / (m n + m χ ) representing the reduced mass of the DM (χ)-nucleon (n) pair system. W (q) = k,α |q·ǫ k,α | 2 4N mN ω k,α is the Debye-Waller factor at zero-temperature. Since the lattice is not perfectly rigid, the Debye-Waller factor accounts for the effect of the quantum and thermal uncertainties of the positions of the nuclei in the scattering. At T = 0 K, only the zero-point fluctuation is relevant. Thus, the total cross section of the DM-target scattering is expressed as where ω p ′ p = p ′ χ 2 / (2m χ ) − |p χ | 2 / (2m χ ) is the energy transferred to the DM particle, and v is its incident velocity.
Note that the sum {kj,αj } runs over all possible phonon vibration modes as the final states. In the above expression, the integration of the out-going DM momentum p ′ χ is traded for that over the transferred momentum q. Since there are n identical phonons in a final state, the integration over momenta is divided by n!. A convenient correspondence , G is adopted in evaluating the amplitude squared. A detailed discussion on the quantization of vibrations in solids using the path integral approach is arranged in Appendix. A.
Moreover, note that the momentum q can be uniquely separated into certain reciprocal lattice G q , and a remainder part [q] within the 1BZ, such that q = G q + [q], and thus the summation over q can be equivalently expressed as the sum Gq [q]∈1BZ . The integration over [q] can always be integrated out from the sum G δ j kj +q, G for an arbitrary set of {k j } without noticeably affecting the values of other integrand functions (· · · ) q that are coarsely dependent on q. The variation of the integrand over the 1BZ is expected to be irrelevant as long as the momentum transfer q = |q| is much larger than the length of the 1BZ, i.e., q ≫ O (1) keV. In this case, one has the following incoherent scattering approximation, where a unique k 0 ∈ 1BZ satisfies G δ i ki+k0, G = 1. This approximation amounts to smoothing out q within the 1BZ as if one can only see a momentum transfer with a resolution comparable to the length of a reciprocal lattice. Next, we further approximate that the simple lattice is isotropic. In this case, eigenenergy ω k,α remains invariant under any rotational operation O R acting on wavevector k, while ǫ k,α also transforms as a vector under the same O R , and thus one has where E R (q) = q 2 / (2m N ). This result also holds for a monatomic cubic system [62]. In the right-hand-side of Eq. (2.5), we relabel the eigenmodes {k, α} with a single notation {i} for brevity, and Eq. (2.3) in the incoherent approximation is recast as where the scattering function S (q, −ω p ′ p ) is defined in the third line, n i represents the occupation number of the energy ω i , and ω = 3N i=1 ωi 3N is the phonon frequency averaged over the density of states (DoS). Note that is a combined Poisson distribution, so the key problem is to determine the probability density of the random variable ω = 3N i=1 n i ω i for this distribution. While it is difficult to derive an analytical expression on a general basis, one can prove that the factor S (q, ω) converges to a Gaussian form in the large q region, i.e., e − (ω−E R (q)) 2 2E R (q)ω / 2πE R (q) ω, by using an argument analogous to that used in the proof of the central limit theorem. Additionally, this Gaussian form converges to the delta function in terms of the energy conservation towards the large q limit, which validates the impulse approximation. In the large q regime, the scattering function S (q, ω) can be approximated with an asymptotic expansion with respect to parameter ω/E R (q) as follows (see Appendix A 8 for a detailed discussion), with To get some sense, taking silicon target for example, in the top row of Fig. 2 we show the non-dimensional function ωS (q, ω) for parameters E R (q) /ω = 5 and 10, respectively. It is evident that in the regime E R (q) /ω 5, the compound Poisson distribution in the third line of Eq. (2.6) already well resembles the asymptotic Gaussian form in shape, except for a minor displacement of the central value E R (q).
For a small transferred momentum q, the asymptotic expansion above is no longer valid. In this case, one can alternatively utilize the functions {T n (ω)} in the last line of Eq. (2.6) so as to calculate the scattering function S (q, ω) in a numerical fashion. Defined as . Bottom: Comparisons between the function ωS (q, ω) of silicon for the multi-phonon distribution (blue histograms) obtained by the numerical recursive method and the impulse Gaussian (red dashed curves) for the ratios ER (q) /ω = 1, 2, and 5, respectively. Similar discussion can be found in Ref. [50]. See text for details. where these {T n (ω)} can be determined by following an iterative procedure (see Appendix A 9 for further details): In the bottom row of Fig. 2 we present the non-dimensional function ωS (q, ω) of silicon target computed with recursive method for parameters E R (q) /ω = 1, 2, and 5, respectively. It illustrates the transition from the multiphonon spectrum into a Gaussion form with an increasing momentum transfer q. Taking typical semiconductors such as silicon for instance, where ω = 40.3 meV, the condition E R (q) /ω 1 translates to a momentum transfer q O (10) keV, which still guarantees the validity of the incoherent approximation. In the limit q → ∞, the width of the Gaussian becomes much smaller than the central value E R (q), and hence the Gaussian reduces to the δ-function δ(ω p ′ p + E R (q)). Then inserting Eq. (2.2), and taking the correspondence q ∼ V (2π) 3´d 3 q, the cross section in the limit q → ∞ becomes which, as expected, is equal to the sum of N incoherent DM-nucleus cross sections for monatomic simple crystal structure.
1 n − 1 n · · · |i |j Figure 3. The diagram for the Migdal effect, where an electron-hole pair is generated, i.e., an electron is elevated from a valence state |j to a conduction state |i , via the exchange of a virtual phonon, along with multiple on-shell phonons generated by the DM external field.

MIGDAL EFFECT AS A MULTI-PHONON PROCESS
The prospect of describing the Migdal effect in terms of phonons and electrons has been originally sketched out in Ref. [51]. Here we explore this approach in more details. The Migdal excitation process is illustrated in Fig. 3, where an electron-hole pair is excited through a virtual phonon, along with a bunch of on-shell ones produced from the collision with a DM particle. In essence, the electron-phonon interaction reflects the Coulomb forces between the distorted ion lattice and the itinerant electrons, for which we provide a short review in Appendix. A 4.
Following the Feynman rules summarized in Sec. A 5, one can read off the amplitude for the process illustrated in Fig. 3, where Z ion is the number of the valence electrons of the material atom, α e is the fine structure constant, and ε i (ε j ) denotes the energy of the conduction (valence) state |i (|j ). The n-phonon sector in the first line has been thoroughly discussed in the preceding section. In the second line lies a phonon mediator with its two ends linking the multiphonon blob and the bare phonon-electron vertex. For typical semiconductors, the band gaps ε g ∼ O (1) eV are much larger than the phonon eigenenergies ω k,α ∼ O 10 −2 eV, so the term in second line can be reduced to In the derivation, we use the contraction relation Recall that in Refs. [50,51] the same amplitude is obtained in the the soft limit, i.e., p N ·(k + G) /m N ≪ ε i −ε j and |k + G| ≪ |p N |, with p N being the momentum of the recoiling nucleus. In the low-energy regime however, neither the soft approximation nor the concept of a freely-recoiling nucleus still holds. In contrast, in the context of electrons and phonons, this particular form of amplitude naturally extends to the low-energy regime. The excitation rate (in the incoherent approximation) can then be written as where ρ χ represents the DM local density, f χ (v) is the DM velocity distribution. Note that the number of the nuclei N in solids, which equals the number of the primitive cells for a simple lattice structure, is explicitly represented with N T here. The factor 2 in the last line counts the two spin orientations for each valence state. For the present we have not taken into account the renormalization effect in our discussion, which can displace the locations of the phonon poles, and induce the screening of the Coulomb interaction. Since the band gaps are far larger than the phonon eigenenergies, only the screening effect that leads to a reduction of the scattering rate is relevant for our discussion.
Here we take the homogeneous electron gas (HEG) for a schematic illustration. As shown in Appendix. A 7 and explained in Ref. [50,51,66], the screening of the electron-phonon vertex adds an inverse dielectric function ǫ −1 (k, ω) to the amplitude analogous to Eq. (3.1) of a crystal structure, while the last line in Eq. (3.3) corresponds to Im [ǫ (k, ω)] at the random phase approximation (RPA) level. Therefore, the overall screening effect is encoded in the energy loss function (ELF) 3 q and k ∼ V (2π) 3´1 BZ d 3 k, the above event rate can be recast as where the nondimensional factor represents the averaged energy loss function, with Ω being the volume of the unit cell, and Im ǫ −1 G,G ′ (k, ω) (see Appendix A 6) being the EFL for the crystal structure. F (ω) has been calculated for diamond and silicon targets in Ref. [51]. Eq. (3.4) applies for the crystal targets that can be considered as isotropic, in which case only the one-dimensional DM speed distribution is relevant for the calculation of the excitation rate in Eq. (3.3), and thus an isotropic velocity distribution f χ (v) is assumed. In the derivation, we first integrate out the angular variable of velocity v with respect to q, which converts to the integral over variable E in Eq. (3.4), and restore the full integration over velocity distribution (by adding a factor 1/2) for convenience. Then we integrate out the solid angle of momentum transfer q using the Legendre addition theorem, which leads to the factor F (ω).
It is interesting to compare the event rate Eq. (3.4) with the one derived in the picture of the bremsstrahlung-like process proposed in Ref. [51], which is expressed as with p N and µ χN = m N m χ / (m N + m χ ) being momentum of the recoiled nucleus and the reduced mass of the DM-nucleus pair, respectively. Θ is the Heaviside step function. In the sub-GeV mass regime, µ χN ≈ m χ . Since the integrand functions {T n (E)} vanish if E < 0, we set the lower limit of the integral to be 0 in Eq.
This relation holds as long as the weight of S (q, E) lies below qv − q 2 2mχ − ω. In the right panel of Fig. 4 we present the expected 90% C.L. sensitivity of silicon target to cross section σ χn with 1 kg·yr of exposure, based on the phonon-mediated (solid ) and bremsstrahlung-like (dashed ) approaches, for a single-electron (blue) and a two-electron (orange) charge bins, respectively, under the zero background assumption.

SUMMARY AND DISCUSSION
In this paper we build up a phonon-mediated description of the Migdal effect in semiconductor targets in the context of the solid state QFT, in which the phonons, namely, the quantized collective vibrations of the ions, rather than the on-shell ions, are used to describe the Migdal excitation process.
In order to ease the discussion, three major simplifications of the problem are made: (1) we assume the solid target is a monatomic simple crystal, possessing an approximately rotational symmetry; (2) we use the zero-temperature QFT in formulating the multi-phonon scattering event rates, rather than a more general finite temperature QFT framework; (3) we take the incoherent approximation in calculating the multi-phonon process. While the isotropy approximation in (1) is valid for the diamond structure materials (e.g., silicon and germanium), it may result in uncertainty for some anisotropic materials (e.g., sapphire). The second assumption is sufficient for experiments operated at cryogenic temperatures. In fact, it is straightforward to generalize the zero-temperature formalism to the finite temperature one by simply replacing the propagators of the zero-temperature case in Eq. (3.1) with those in the finite temperature scenario. As for the third approximation, note that for DM masses around an MeV, the typical momentum transfer can be as small as q ∼ O (1) keV, which is comparable to the size of the 1BZ, and hence beyond the regime of validity for the incoherent approximation. Further study for m χ < 1 MeV is needed.
Based on the formalism, we numerically calculate the Migdal excitation event rates for the silicon semiconductor target. As expected, the multi-phonon energy spectra are found to well converge to the Gaussian form in the large q limit, justifying the impulse approximation used in the bremsstrahlung-like description. Although the behavior of the phonon scattering function S (q, ω) differ from that of a free nucleus in the low and intermediate scattering energy region, the Migdal excitation rates calculated from the impulse approximation are found to be well consistent with that obtained using the phonon-mediated approach throughout the relevant DM mass range. Finally, it is tempting to apply the phonon-mediated approach to the probe of sub-MeV DM particles through the Migdal effect in novel narrow-gap materials with band gaps of O (10) meV (e.g., Dirac materials), where the picture of the free-recoiling nucleus turns invalid altogether. We leave it for the future work.
Note added. After this work was published, Kim Berghaus suggested us that the contribution of T 0 was omitted in our original numerical implementation of Eq. (3.4), which as a consequence can remarkably suppress the calculated event rates for a small q (an upcoming paper [68] also discusses the Migdal effect in semiconductor detectors). After T 0 term is included, it is found that the Migdal event rates calculated from the phonon-mediated approach and the impulse approximation coincide quite well even in the low DM mass range. For easy reading, we provide an elementary introduction to relevant theoretical background to the main text in this appendix, including the treatment of the phonons and electrons, as well as their interactions in the context of the QFT. Part of the material can be found in Ref. [69].

Quantization of vibrations in solids
For simplicity here we only consider the case of the monatomic simple lattices. The dynamics of the crystal vibration is described with the following equation, where m N is the nucleus mass, u ℓ is the displacement of the nucleus at lattice site ℓ, and the force strength matrix elements are explicitly expressed as follows, with U being the potential between the nuclei at sites ℓ and ℓ ′ , and σ, σ ′ = {x, y, z} denoting the three space directions. The Fourier transform of the force strength matrix is called the dynamical matrix V (k), i.e., which is real symmetric matrix for the Bravais lattice at each wave-vector k in the 1BZ, and thus can be diagonalized with an orthonormal basis set of vectors {ǫ k,1 , ǫ k,2 , ǫ k,3 }, such that with ω (k, α) being corresponding eigenfrequency of the mode (k, α). Besides, since Φ (ℓ) = Φ (−ℓ) for Bravais lattice, it is straightforward to see that V (k) = V (−k), ǫ k,α = ǫ −k,α and ω k,α = ω −k,α . With these preparations, one first substitutes the displacements with eigen-vibration modes: where Q k,α encodes the vibration amplitude for the mode (k, α), and N is the number of the unit cells in the material, and then obtains the Lagrangian of vibration system, and the equations of motionQ One then follows the conventional quantization procedures to give pairs of the canonical position and momentum operatorsQ k,α (t) = 1 2 ω k,α â k,α e −iω k,α t +â † −k,α e iω k,α t , that satisfy the equal time commutation relation Q k, , and vice versa. In the context of the path integral quantization, the action is written aŝ from which one constructs the free phonon propagator

Quantization of electrons in solids
One can construct the path integral formalism for electrons in solids in a similar fashion, except for the anticommutation nature of the Grassmann algebra. Here we summarize some important results.
The action of the electron field can be drawn from the Schrödinger equation as the following: where −∇ 2 ψ e /2m e + V is Hamiltonian for a single electron, with {u i (x)} and {ε i } being its eigenwavefunctions and corresponding energies, respectively. Thus one obtains the electron propagator (ω p ′ p , q) Figure 5. The effects of the incident DM particle on the target material can be regarded as an external field.
The coupling term in the action between the DM particle and nuclei in solids can be directly written as where the interaction is instantaneous such that Since one is only interested in the target material that hosts the phonons and electrons, it is convenient to integrate out the DM component and regard it as an external field. To be specific, one draws DM external field from the amplitude of the scattering process χ (p χ ) +target→χ p ′ χ +target (excited) (see Fig. 5 for illustration) and obtains the effective Lagrangian · · · · · · · · · · · · + + = + · · · Figure 6. The diagram for the n-phonon scattering process containing full contributions of phonon loops, which generates the Debye-Waller factor denoted as the gray blob on the right-hand-side of the equation. See text for details.
is the Fourier transform of the DM-nucleus contact interaction V χN (x). It should be noted that in above discussion we adopt the discrete momentum convention.
Based on above discussion, one can derive the LSZ reduction formula for the multi-phonon scattering process. For example, the S -matrix for an n-phonon scattering process subject to a DM external field can be expressed as where {k i , α i } (i = 1, 2, · · · , n) label the n-phonon final states. Now one is able to investigate the amplitude for the full n-phonon scattering process shown in Fig. 6, which contains not only the tree level piece, but also higher orders of the phonon loop. Note that a self-closed propagator for mode (k, α) reads as as well as a companion factor −iq·ǫ k,α √ N mN −iq·ǫ k,α √ N mN . Thus every loop in diagram corresponds to a factor − k,α |q·ǫ k,α | 2 2N mN ω k,α . On the other hand, a specific external leg (k i , α i ) corresponds tô Besides, the effect of the symmetry factor should also be taken into account. For example, one considers the n-phonon process containing m self-interacting loops in Fig. 6. Determining the number of the contractions in Eq. (A.16) is equivalent to enumerating all possible ways the internal phonon lines interconnect among themselves, which is illustrated in Fig. 7. It is not difficult to verify that the overall constant that encodes the symmetry effect is equal to · · · · · · 1 n · · · · · · · · · 1 n n + 2m Figure 7. Illustration of possible ways in which n phonon external lines connect to the n internal ones while the left 2m internal members self-connect with each other to form m loops. See text for details.
counts which 2m lines self-connect among all n + 2m lines from the left in Fig. 7, (2m)! m! 2 m is the number of the ways these specific 2m lines connect with each other, n! describes the interchange of the external legs at the right hand in Fig. 7, and 1 (n+2m)! corresponds to the factor 1 s! in Eq. (A. 16). Putting all these pieces together, the sum of all diagrams on the left-hand-side in Fig. 6 can be expressed as where e − k,α | q·ǫ k,α| 2 4N m N ω k,α is no other than the Debye-Waller factor at the zero-temperature, which is represented with the gray blob on the right-hand-side of Fig. 6. In the derivation we interchange ǫ k,α and ǫ −k,α whenever necessary. From above discussion we can see one benefit of the path integral approach: one no longer has to resort to the cumbersome operator commutator algebra to obtain the Debye-Waller factor. Propagators do the job.

Electron-phonon interaction
The interaction between the ions and electrons can be directly written as where U e (x − ℓ) = −Z ion α e / |x − ℓ| is Coulomb potential between the ion located at ℓ and an electron at position x. Thus, the electron-phonon interaction Lagrangian term is written as with v (k + G) = 4πα e / |k + G| 2 .

Feynman rules
Based upon above preparation, the Feynman rules describing the processes involving multi-phonon and electronphonon interactions in momentum space can be summarized as follows.
• In discussion the vertex of the DM particle is replaced with an external field as shown in Fig. 5. Such an external source corresponds (−i) V χN (q) /V .
• Each blob also contributes the energy-momentum conservation condition presented as discrete delta functions • 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 ε i ), where p is the net momentum sinking into the vertex.
• A phonon-electron vertex is read as • A phonon internal line with one end connecting a phonon blob and the other connecting an electron vertex corresponds to the sum over propagators of all modes {k, α}, i.e.,

Random phase approximation
A short review of the random phase approximation (RPA) has been provided in Ref. [51], so here we only summarize some results relevant for our present discussion. Within the framework of the RPA, the Lindhard dielectric function for homogeneous electron gas (HEG) is expressed as where n i (n j ) denotes the occupation number of the state |i (|j ), with ε i (ε j ) being corresponding eigenenergy. In crystalline structure the translational symmetry for space-time reduces to that for the periodic crystal lattice. The momentum transfer q in Eq. (A.22) is expressed uniquely as the sum of a reciprocal lattice vector G, and corresponding reduced momentum k confined in the 1BZ, i.e., q = k + G. In this case, the microscopic dielectric matrix is used to describe the screening effect in solids. For more details, see Ref. [51].
= + + + · · · = + + · · · + Figure 8. Top: The renormalized electron-phonon vertex (black square box) presented at the RPA level, where the one-particleirreducible blob is represented by an electron-hole pair loop. Note that the gray lines represent the phonon propagators, while the black wiggly lines represent the Coulomb interaction. Bottom: The dressed phonon line (double wiggly line) presented at the RPA level. See text for details.

Electron-phonon interaction renormalization at the RPA level
Here we discuss how to renormalize both the bare phonon propagator and effective electron-phonon interaction at the RPA level. In RPA, the polarization bubble is approximated as a simple electron-hole pair bubble, and hence the dressed electron-phonon vertex shown in the top row of Fig. 8 can be expressed (when the reciprocal lattice vectors are suppressed) as the sum where (−i) Π corresponds to the electron-hole loop. It is evident that the effect of the normalization is suppressing the bare vertex with the dielectric function ǫ (k, ω). The electron-phonon interaction also bring about a correction to the position of the phonon poles (in the RPA), as shown in the bottom row of Fig. 8. Since the typical band gaps are much larger than the phonon eigenenergies, these small corrections are irrelevant for our purpose in this study. One can use the DarkELF package to take into account the screening effect in various materials [70].

Asymptotic behavior of the multi-phonon distribution
Here we investigate the asymptotic behavior of the combined distribution of n independent Poisson distributions encountered in Sec. 2, which justifies the validity of the impulse approximation in the large q regime. To make the discussion concise, the following parameters are introduced: , and then consider the distribution P (n 1 , n 2 , · · · , n 3N ) = p (n 1 , λµ 1 ) p (n 2 , λµ 2 ) · · · p (n 3N , λµ 3N ) where p (n i , λµ i ) is the Poisson distribution of variable n i with mean λµ i . We first try to obtain the distribution of a random variable in the form of z = 3N i=1 (ǫ i /λ) n i , and λ being a parameter that characterizes the scale of the problem. To achieve this goal, we calculate the relevant characteristic function ϕ Z (t) = {ni} P (n 1 , n 2 , · · · , n 3N ) e i zt λµ i e i(ǫi/λ)t − 1 , (A. 26) with which the distribution of the variable x is explicitly expressed and expanded in λ (λ ≫ 1) as follows, ω 2 i /ω 2 in the last line. The integral in the last step can be evaluated by integrating over its saddle point z 0 = −i √ λ (x − 1) along the path (z 0 − ∞, z 0 + ∞), on which z 3 term is suppressed by 1/ √ λ, as long as |x − 1| is not too much larger than the width 1/ √ λ.

Iterative calculation of multi-phonon process
Here we discuss how to calculate the multi-phonon spectrum at T = 0 K following a recursive procedure. First, we explicitly derive the scattering factor introduced in Eq. (2.6) as the following,  Then it is easy to verify the following recursive relation (see Ref. [62] for a general discussion on the case of a finite temperature) It is straightforward to see which means T 1 (ω) = 0, if ω < 0. This feature can be easily generalized to the case of an arbitrary p such that T p (ω) = 0, (ω < 0). In practice, we utilize Eq. (A.29) to obtain {T n (ω)} and to further calculate the spectrum of the multi-phonon process with Eq. (2.6). This recursive method requires only the phonon DoS for a solid target (e.g., for monatomic simple lattice) as follows, which is normalized such that´+ ∞ 0 g (ω) dω = 1. For illustration, in Fig. 9 we present the normalized DoS for the bulk silicon, which is calculated using PHONOPY code [71], while the force constants are computed using VASP package [72] based on the density functional theory [73,74] with Perdew, Burke, and Ernzerhof form [75] of the generalized gradient approximation on the exchange-correlation functional.