Ab initio electron-two-phonon scattering in GaAs from next-to-leading order perturbation theory

Electron-phonon (e–ph) interactions are usually treated in the lowest order of perturbation theory. Here we derive next-to-leading order e–ph interactions, and compute from first principles the associated electron-two-phonon (2ph) scattering rates. The derivations involve Matsubara sums of two-loop Feynman diagrams, and the numerical calculations are challenging as they involve Brillouin zone integrals over two crystal momenta and depend critically on the intermediate state lifetimes. Using Monte Carlo integration together with a self-consistent update of the intermediate state lifetimes, we compute and converge the 2ph scattering rates, and analyze their energy and temperature dependence. We apply our method to GaAs, a weakly polar semiconductor with dominant optical-mode long-range e–ph interactions. We find that the 2ph scattering rates are as large as nearly half the value of the one-phonon rates, and that including the 2ph processes is necessary to accurately predict the electron mobility in GaAs from first principles.

(2ph) nk in GaAs at 300 K, computed using only the LO mode (curve labelled "2ph polar") and with all phonon modes included (curve labelled "2ph all modes").

Supplementary Figure 2
Supplementary Figure 2. Mobility calculation in BaSnO3. Electron mobility in BaSnO3 computed using the same methods as in Fig. 7 of the main text. It is seen that the iterative BTE with only 1ph scattering overestimates the mobility, while including 2ph scattering improves the agreement with experiment, similar to GaAs. Experimental data are shown with markers at room temperature and dashed lines at lower temperatures. The experimental data are from Figs. 4 and 5 in Ref. [1]. We used the single crystal mobility data for a 0.5 − 1.5 × 10 20 cm −3 carrier concentration range. Our calculation used a 0.8 × 10 20 cm −3 carrier concentration.  In this section, we calculate the contributions to the electron-phonon (e-ph) scattering rates from the next-toleading-order self-energy diagrams. We use the Matsubara technique to calculate the two-loop self-energy, whose imaginary part is related to the total scattering rates via the optical theorem (see Supplementary Figure 4). We focus on the scattering processes with two external phonons.
The Feynmann rules, which have been derived in Mahan [2], will be adapted here to our context. The starting point is the e-ph Hamiltonian where a nk and b νq are the annihilation operators for electrons and phonons with energies ε nk andhω νq , respectively, g mnν (k, q) is the e-ph coupling constant, and N Ω is the number of unit cells in the crystal. Comparing with Eqs.
[2], we have introduced the dependence on electron crystal momentum k for the e-ph couplings and the electron band indices m and n. Also note that for the Hamiltonian to be Hermitian, the e-ph couplings must satisfy g * mnν (k, q) = g nmν (k + q, −q).
3. With each vertex, associate the e-ph coupling constant g mnν (k, q). Beware of the direction of q.
4. Conserve momentum and complex frequency at each vertex and sum over the internal degrees of freedom.

Multiply the expression by
where F is the number of closed Fermion loops. The (2S +1) factor is a summation over spin degrees of freedom, and 2S + 1 = 2 for electrons. The integer L is the number of loops, and β = 1/k B T , where T is temperature.

Electron Self-Energy
We consider below the 1-loop diagram that gives the lowest-order (one-phonon) self-energy and the three relevant two-loop diagrams for the electron self-energy. Diagram IIc will not contribute to the two-phonon processes and thus will not be considered in the following.

One-Loop Diagram I
As a warm up exercise, we first derive the one-loop self-energy diagram, labelled as I in figure above. Since L = 1 and F = 0, the Feynman rules give To apply the Matsubara frequency summation method, we define the bosonic weighting function as in [2]: whose poles are at iω κ = i2κπ/β, with residues 1/β and integer κ values. The weighting function for fermions is n F (z) = 1 e βz + 1 , whose poles are at ik λ = i(2λ + 1)π/β, with residues −1/β and integer λ values. For this diagram, we will do the contour integral for f (z)n B (z) at the complex infinity. Since f (z)n B (z) decays faster than 1/z, we can apply the Cachy residue theorem, which gives (here and below, z are the relevant poles): Using this result, we get: The poles of f (z) are at z 1 = ω νq , z 2 = −ω νq and z 3 = −ik λ + ξ mk+q . Their residues are We also know that where N and f are the thermal occupation numbers for phonons and electrons, respectively. We also used the fact that ik λ = i(2λ + 1)π/β. Substituting this result in the self-energy expression, we get Employing the analytic continuation ik λ → E + iη, we obtain the off-shell lowest-order e-ph self-energy: We will be mainly interested in the scattering rate at the electron energy ξ nk and therefore we will set E = ξ nk to obtain the on-shell self-energy for the state with band n and crystal momentum k. Using the identity and Eq. (7.304) in Ref. [2], which states that the scattering rate Γ is obtained as Γ = −(2/h)ImΣ, we get withh placed back into the expression. This is the well-known lowest-order scattering rate commonly used in firstprinciples calculations.

Two-Loop Diagram IIa
Now we compute the first two-loop diagram, which is shown in the figure below. Since L = 2 and F = 0, the Feynman rules give where in A κ we collect all terms independent of iω σ . Let us sum over iω σ first. Performing the contour integral for Note that f (iω κ , z) has three poles, whose residues and bosonic weight factors are computed as: where we defined two functions, h (−) (iω κ ) and h (+) (iω κ ), as the first and second terms in the expression above. We then sum over iω κ , using again A subtle point is that the cases with n 3 = n 1 and n 3 = n 1 have different pole structures, and need to be discussed separately (see the figure above for diagram IIa; n 3 and n 1 label two intermediate electronic states in the self-energy diagram). Luckily, the two cases give the same expression for the two-phonon scattering processes, as we show explicitly below. Before carrying out the calculation, let us introduce some useful abbreviations. We will use in the Let us focus on h (−) for the case n 3 = n 1 first. In this case, h (−) is defined as It has five poles, which are given here together with their residues and bosonic weight factors: The terms related to the two-phonon processes are those containing 1/(ik − ξ 2 ± ω p ± ω q ). After repeating this procedure for h (+) , collecting terms and ignoring terms that are irrelevant to the two-phonon processes, we get The rates of the two-phonon processes emerge after analytically continuing ik λ to E + iη and taking the imaginary part of 1/(E − ξ 2 ± ω p ± ω q + iη). We also use the delta functions to set E = ξ 2 ∓ ω p ∓ ω q in some of the denominators. After carrying out these calculations, we obtain In this case, h (−) is defined as The function h(z)n B (z) has a pole of order 2 at z 1 = −ik λ + ξ 1 . By employing where n is the order of the pole, we get After substituting z 1 = −ik λ + ξ 1 , n B (z 1 ) = −f 1 and n B (z 1 ) = βf 1 (1 − f 1 ), we get The other three poles are simple poles and can be treated in the usual way. Repeating this procedure for h (+) and adding all the contributions, we get We perform the analytic continuation ik λ → E + iη, take the imaginary part of 1/(E − ξ 2 ± ω p ± ω q + iη), and get The second two-loop diagram, called here IIb, is shown in the figure below. Since this diagram also has L = 2 and F = 0, the Feynman rules give Using a notation we introduced above, A κ collects all terms independent of ω σ , and we use again abbreviations introduced in the previous section, such as ξ 3p ≡ ξ n3k+p , etc. Summing over iω σ first, we get We then sum over iω κ , and collect the relevant terms for two-phonon scattering processes, which are given below: After performing the analytic continuation and taking the imaginary part, we get

Two-Phonon Scattering Rates
Collecting the contributions from diagrams IIa and IIb, using Γ = −(2/h)ImΣ, and setting E to the band energy ξ nk , the scattering rate of the two-phonon processes becomes where we introduce the process amplitudes The scattering rate is proportional to the absolute square of the scattering amplitude: Therefore, we will insert the infinitesimals in a way that allows us to express the scattering rates in an absolute square form. To achieve this, first note that quantities such as q, p, ν and µ are dummy variables that are summed over or integrated, so we can rename them at will. Let us denote (νq ↔ µp) the term with its dummy variables swapped in the way indicated by the arrows, ν ↔ µ, q ↔ p, etc. Let us consider the processes with one phonon absorption and one phonon emission first, that is, the sum of terms (i) and (iv): Since the expression is already in square form, we have inserted the iη terms in a way that makes the expression become an absolute square. Using a similar approach for the process in which the electron emits two phonons, and for the process in which the electron absorbs two phonons, We thus get:

Resonance
The last expression is very close to the final result given in Eqs.
(1)−(4) of the main text. The last problem we need to solve is that the sum giving Γ (2ph) nk in the expression above diverges when the intermediate electron state is on shell, in which case the denominator in the γ terms given above vanishes, resulting in a divergent scattering rate. This phenomenon is called resonance. The problem is that the intermediate state will eventually transition into a different state, but using free propagators 1/(E − ξ + iη) for the on shell intermediate states implies an infinite intermediate state lifetime. The common practice in this situation, which also arises in other quantum field theories, is to consider the full electron propagator 1/(E − ξ + iη − Σ) as shown in the figure below, which introduces a finite lifetime for the intermediate electronic state. Diagramatically, this approach is equivalent to performing a resummation of diagrams to all orders, as is done in the well-known GW self-energy {see Eq. (5.54) in Ref.
[2]}. For our 2ph scattering rate expression, we simply add the intermediate state self-energy in the denominators of all the γ terms above, which removes the divergences.
The square amplitudes γ (i) for the different processes, i = 1e1a, 2e and 2a, are defined as where we have taken into account the resonance by adding the intermediate state self-energy in the denominators. The factors of A (i) contain the thermal occupation numbers of electrons and phonons, and are defined as These are our final expressions for the two-phonon scattering rates, which are given in Eqs.
Supplementary Figure 5. Temperature dependence of the 2ph scattering rates. Temperature dependence of the ratios of the 2ph scattering processes to the leading order e-ph scattering rate. From left to right, the panels are for electronic states with energies of 20, 45 and 90 meV above the conduction band minimum, and thus respectively in region I, II and III defined in the main text. and III defined in the main text. The ratios are give for both the total 2ph scattering rate and (for completeness) for the individual 2ph processes, 1e1a, 2e and 2a. In the 200−500 K temperature range considered in this work, the ratio of the total 2ph scattering rate to the 1ph scattering rate is nearly temperature independent in all three energy regions.

Supplementary Note 3: Boltzmann Transport Equation with Two-Phonon Contributions
Let us briefly summarize the formulation of the linearized Boltzmann transport equation (BTE) incorporating the 2ph scattering processes. A more extensive explanation of its formulation and derivation will be presented in a separate work. Defining the total e-ph scattering rate as where the relaxation time τ nk is the inverse of the total scattering rate, namely τ nk = 1/Γ nk , and the first and second terms in brackets are due to the lowest-order (1ph) and 2ph scattering processes, respectively. The function F nk is the unknown in the equation, and F 0 nk = τ nk v nk , with v nk the band velocity. After solving the equation, the electrical mobility in direction i can be obtained using where V uc is the unit cell volume and n c the charge carrier concentration. A common approach to computing F nk is the relaxation time approximation (RTA), which neglects the second term on the right hand side of Eq. (3) and approximates F nk to F 0 nk . A more accurate solution can be obtained through an iterative approach (ITA) [3], in which, starting from the RTA solution F nk = F 0 nk , one iteratively substitutes F nk in the right hand side of Eq. (3) to obtain an updated F nk , until a converged F nk is reached.
Our work presents calculations, both within the RTA and ITA, in which the 2ph processes are either included or neglected; when only 1ph processes are included, τ nk is set to the inverse of Γ (1ph) nk , and the second term in brackets in Eq. (3) is neglected. The ITA with 2ph contributions included is the most accurate level of theory, and the one that agrees best with experiment, while the ITA with only 1ph processes overestimates the experimental result. Note that interference between 2ph processes is neglected in the mobility calculation. This is a good approximation if for most of the (q, p) pairs there is only one process in Eq. (3) in the main text dominates, as in the case of GaAs that we have verified. [3] W. Li, Phys. Rev. B 92, 075405 (2015).