## Abstract

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.

## Introduction

Electron–phonon (*e*–ph) interactions are essential to understanding electrical transport, nonequilibrium dynamics and superconductivity. Using density functional theory (DFT) and related methods, it has become possible to compute *e*–ph interactions from first principles, and use them to predict the carrier scattering rates and mobilities, both in simple and in complex materials with up to tens of atoms in the unit cell^{1,2,3,4,5,6,7,8,9}. In the typical workflow, one takes into account only the leading-order *e*–ph scattering processes, which involve scattering of the carriers with one phonon. Nearly all ab initio work to date has relied on such leading-order perturbation theory, tacitly neglecting higher-order *e*–ph processes.

Yet many compounds, including polar semiconductors, oxides and organic crystals, exhibit *e*–ph interactions that cannot be treated with lowest-order perturbation theory. Polar materials are a paradigmatic case in which *e*–ph interactions due to polar optical modes—in particular, the longitudinal optical (LO) mode—are long-range^{10,11}, which can lead to higher-order *e*–ph scattering and polaron formation. In weakly polar materials, such as III-V and II-VI semiconductors and high-mobility oxides, one expects that perturbation theory is still valid, but that higher-order *e*–ph processes are significant and need to be included.

Studies of second- and even third-order *e*–ph processes in semiconductors exist^{12,13,14,15,16,17}, but they are limited to simplified models, restricted to the conduction band minimum, or only valid at zero temperature, and are therefore inadequate for quantitative predictions. In first-principles *e*–ph and charge transport calculations, next-to-leading-order effects have been typically assumed to be negligible, even though their contributions have been estimated to be important using simplified models^{13,14,15,18,19}. In materials with intermediate or strong *e*–ph coupling, such as oxides with low room-temperature mobility and ionic compounds (e.g., alkali halides), the *e*–ph interactions can lead to regimes beyond the reach of perturbation theory, including the formation of polarons^{20}. This coupling regime has been investigated with diagram-resummation techniques such as the cumulant method, both analytically^{21} and more recently ab initio^{22,23,24}.

Higher-order processes are generally important in quantum field theories of condensed matter. An example are light–matter interactions, where phonon-assisted^{25} and two-photon^{26,27} absorption have been studied extensively. Even calculations in quantum electrodynamics beyond the leading order can provide important lessons—for example, higher-order corrections are essential to accurately predict large-angle Bhabha scattering in electron-positron collisions^{28}, and calculations up to the tenth order have been carried out to compute the magnetic moment of the electron^{29}.

However, it is a daunting task to systematically go beyond the leading order due to the rapid increase in the number of Feynman diagrams and their computational complexity. For *e*–ph interactions, the next-to-leading order involves electron scattering events with two phonons, which requires computationally challenging Brillouin zone integrals over two crystal momenta. This computational complexity has hampered next-to-leading-order *e*–ph calculations for decades^{12}, and the relevant next-to-leading-order diagrams (see Fig. 1) have not yet been computed numerically from first principles.

Here we formulate and compute from first principles next-to-leading-order *e*–ph interactions, focusing on electron scattering processes involving two phonons (hereafter denoted as 2ph processes). We compute and analyze their contributions to the *e*–ph scattering rates, using as an ideal case study a weakly polar semiconductor, GaAs. We find that the 2ph scattering rates are surprisingly large, and comparable in magnitude to the lowest-order rates due to one-phonon (1ph) processes. Our analysis shows that the relative importance of the 2ph contributions is nearly temperature independent at 200–400 K, and rationalizes the peculiar energy dependence of the 2ph processes. The results are sensitive to the lifetimes of the intermediate states, which need to be included to avoid divergences due to resonance effects; we develop a self-consistent scheme to overcome this challenge. The 2ph processes are also shown to play an important role in accurately predicting the electron mobility in GaAs. We formulate and iteratively solve a linearized Boltzmann transport equation (BTE) that includes the 2ph processes, showing that this level of theory can correct the discrepancy with experiment of the mobility predicted with the BTE including only 1ph processes. Our work proposes an approach for systematically improving the accuracy of ab initio *e*–ph calculations beyond the leading order. While this method is broadly applicable, it is particularly well-suited for weakly polar (III-V and II-VI) semiconductors and high-mobility oxides, in which the weak polaron effects can be treated perturbatively.

## Results

### Two-phonon processes and their scattering rates

We use the Matsubara frequency sum technique^{30} to derive the electron self-energy due to *e*–ph interactions up to \({\mathcal{O}}\)(*g*^{4}) (see Fig. 1), where *g* is the *e*–ph coupling constant; our treatment focuses on computing the imaginary part of the self-energy and the related 2ph scattering rates. The derivations are lengthy and tedious, and are given in detail in the Supplementary Information.

All 2ph processes considered in this work consist of two consecutive 1ph scattering events^{12}, as we show in Fig. 2. They include processes in which the electron absorbs one phonon and emits another phonon, or vice versa (both processes are denoted as 1e1a), and processes in which the electron emits or absorbs two phonons (denoted as 2e and 2a, respectively). The intermediate electronic state accessed after the first scattering event is associated with a virtual electron that cannot be observed and whose energy can take any value. Processes are defined as on shell when the intermediate electron energy equals a band electron energy, and off shell when it does not (see Fig. 2). As we discuss below, on-shell processes lead to divergent integrals, which are regularized by including the intermediate electron lifetime.

The key result of our analysis is the scattering rate due to the 2ph processes, \({\Gamma }_{n{\bf{k}}}^{(2{\rm{ph}})}\), for an electronic state with band index *n* and crystal momentum **k**, which can be expressed as

where the process-resolved 2ph scattering rates \({\widetilde{\Gamma }}^{({i})}\) (*i* = 1e1a, 2e or 2a) depend on the two-phonon momenta **q** and **p** and their respective mode indexes *ν* and *μ*, and *n*_{2} is the band index of the final electronic state, whose momentum is fixed to \({{\bf{k}}}^{\prime}\equiv {\bf{k}}+{\bf{q}}+{\bf{p}}\) by momentum conservation; to correctly normalize the sum, we divide it by \({N}_{\Omega }^{2}\), the number of points sampled in the (**q**, **p**) space. The process-resolved 2ph scattering rates are defined as

where *ε* are electron energies relative to the chemical potential and *ω* are phonon energies; the delta function imposes energy conservation and the constants *α* for each process are defined as

The square amplitudes of the three processes are

where *n*_{1} is the band index, Σ is the self-energy of the intermediate electronic state and *η* is a positive infinitesimal. The prefactors *A*^{(i)} contain the thermal occupation numbers of electrons and phonons (denoted by *f* and *N*, respectively) and are defined as:

Let us examine the process amplitudes in Eq. (3). Since all three processes have similar expressions, we focus on a 1e1a process in which the electron absorbs one phonon and then emits another as an example, though our considerations are general. In this process, which is shown schematically in Fig. 3, an electron first absorbs a phonon with energy *ω*_{νq}, transitioning to an intermediate state \(\left|{n}_{1}{\bf{k}}+{\bf{q}} \right\rangle\) with energy *E* = *ε*_{nk} + *ω*_{νq}, and then emits a phonon with energy *ω*_{μp}, reaching the final state with energy \({\varepsilon }_{{n}_{2}{{\bf{k}}}^{\prime}}\). Note that the energy *E* of the intermediate (virtual) electron does not need to match its band energy \({\varepsilon }_{{n}_{1}{\bf{k}}+{\bf{q}}}\). The amplitude for this process, from Eq. (3), reads

where we omitted the *i**η* infinitesimal for clarity. This expression can be understood as follows: the two *e*–ph coupling constants are associated with each of the two consecutive 1ph scattering events, and the denominator is due to the propagator of the intermediate electron, which is proportional to \(1/(E-{\varepsilon }_{{n}_{1}{\bf{k}}+{\bf{q}}}-{\Sigma }_{{n}_{1}{\bf{k}}+{\bf{q}}})\). The amplitude of the 2ph process is thus inversely proportional to its off-shell extent, \(\Delta E=E-{\varepsilon }_{{n}_{1}{\bf{k}}+{\bf{q}}}\), implying that processes with large intermediate off-shell extents Δ*E* are unlikely to occur. In an on-shell 2ph process, *E* is equal (or very close) to the intermediate-state band energy \({\varepsilon }_{{n}_{1}{\bf{k}}+{\bf{q}}}\), resulting in a small denominator Δ*E* ≃ 0 and thus in a large amplitude. The inverse lifetime of the intermediate state, which is proportional to \({\rm{Im}}{\Sigma }_{{n}_{{\rm{1}}}{\bf{k}}+{\bf{q}}}\), prevents the on-shell process amplitude from diverging.

### Numerical calculations and intermediate lifetime

We implement the expressions derived above in our in-house developed code, Perturbo^{31}, and carry out first-principles calculations on GaAs (see Methods). In our approach, we sum the 2ph scattering rates in Eq. (1) over all possible pairs of phonons, and thus include both on-shell and off-shell processes on the same footing. The treatment of the intermediate lifetimes in our implementation deserves a detailed discussion. When the intermediate state is on shell, the scattering process results in resonance effects, and as discussed above the intermediate state lifetime is crucial to prevent the 2ph scattering rate from diverging. Here and below, the intermediate state lifetime is defined as the inverse scattering rate of the intermediate state, which is obtained from the imaginary part of its self-energy as \(2\ {\rm{Im}}\Sigma /\hslash\).

The 2ph scattering rates depend on the intermediate state self-energy, Σ in the denominator of Eq. (3), whose value needs to be chosen carefully. We neglect its real part, which only corrects the band structure and barely affects the 2ph calculation. The imaginary part of Σ includes in principle scattering effects from all possible sources. In practice, we approximate \({\rm{Im}}\Sigma\) with the total *e*–ph scattering rate, including both the lowest-order (1ph) rates and the 2ph rates, using \(\left|{\rm{Im}}\Sigma \right|=\hslash /2[{\Gamma }^{(1{\rm{ph}})}+{\Gamma }^{(2{\rm{ph}})}]\), where Γ^{(1ph)} is the usual 1ph scattering rate^{3}. This approach makes Eq. (1) a self-consistent problem. We iterate Eq. (1) until the 2ph scattering rate equals the 2ph contribution to the intermediate state lifetime. In each iteration *m*, the lifetime is due to the sum of the lowest order plus the 2ph scattering rates at the previous iteration, namely, \(\left|{\rm{Im}}\Sigma (m)\right|=\hslash /2\ [{\Gamma }^{(1{\rm{ph}})}+{\Gamma }^{(2{\rm{ph}})}(m-1)]\). The initial Γ^{(2ph)} is set to zero, and the convergence process is performed separately at each temperature.

### Analysis of the two-phonon scattering rates

Figure 4 shows the first iteration and the converged result for the 2ph scattering rates in GaAs at 300 K, for states near the bottom of the conduction band (see Methods). In this energy range, the *e*–ph interactions in GaAs are dominated by the long-range Fröhlich interaction due to the LO mode^{3}, with nearly negligible contributions from all other phonon modes. The converged 2ph scattering rates are surprisingly large—they are smaller than the leading-order 1ph rate, thus justifying the perturbative approach, but they are nearly half the value of the 1ph rates at all energies. We find that the rainbow diagram (the first diagram in the second row of Fig. 1) provides the dominant contribution to the 2ph scattering rate, while the contribution from the vertex-correction diagram (the second diagram in the second row of Fig. 1) is negligible, in agreement with Migdal’s theorem. Therefore, lowest-order perturbation theory can capture only part of the dynamical processes due to the long-range *e*–ph interactions with the LO mode in GaAs. Though the results shown in Fig. 4 include only the dominant *e*–ph interaction with the LO mode, we have verified that including all phonon modes gives nearly unchanged 2ph scattering rates in GaAs (see Supplementary Fig. 1).

The 2ph scattering rate exhibits a trend as a function of energy with three plateaus near the bottom of the conduction band. This trend is a consequence of the dominant *e*–ph interactions in GaAs, which consist of absorption or emission events of LO phonons with energy *ω*_{LO} ≈ 35 meV. To rationalize the energy dependence of the 2ph scattering rate, we define three energy regions, denoted as I, II and III in Figs. 5 and 6, which correspond respectively to electron energies below *ω*_{LO}, between *ω*_{LO} and 2*ω*_{LO}, and greater than 2*ω*_{LO} referenced to the conduction band minimum, which is taken hereafter to be the electron energy zero. The band structure in all three energy regions is nearly parabolic. We plot the rates of each of the 1e1a, 2e and 2a process contributions in Fig. 5, and sketch the LO phonon emission and absorption processes in each energy region in Fig. 6.

In region I, the electrons possess an energy smaller than *ω*_{LO}, and thus cannot emit two LO phonons since this would require a final state in the band gap; the rate of the 2e process is accordingly negligible in region I. The 2a process in which the electron absorbs two LO phonons is active in region I, but it is thermally activated and weak at 300 K since *ω*_{LO} ≈ 35 meV. The 2a process is nearly energy independent, and weak throughout the three regions. The 1e1a process in which an electron absorbs and them emits an LO phonon is the dominant contribution in region I. Recall that the 2ph rate is inversely proportional to the square of the off-shell extent, Δ*E*^{2}. The 1e1a channel in which an LO phonon is first emitted and then absorbed is thus suppressed in region I since its intermediate state, resulting from one LO phonon emission, is always off shell, as shown in Fig. 6a.

In region II, the 1e1a scattering rate increases rapidly at energy *ω*_{LO} = 35 meV because the intermediate state accessed after emitting an LO phonon can be on-shell, as shown in Fig. 6b. The rate of the 2e processes remains negligible in region II up to an energy of 2*ω*_{LO}. In region III, the contribution from the 2e process increases substantially as the electrons can emit two LO phonons and transition to the bottom of the conduction band through an on-shell process, as shown in Fig. 6c. The 1e1a scattering rate drops from region II to region III due to subtle reasons related to the lifetime of the intermediate state. Recall that for on-shell processes, Δ*E* ≃ 0, and thus the amplitude is proportional to the intermediate state lifetime (see Eq. (5)). An electron in region II can emit a phonon, transition to an on-shell intermediate state in region I, and then absorb another phonon to transition to a final state, as shown in Fig. 6b. As the scattering rates for states in region I are considerably smaller than in the other two regions, the on-shell intermediate states have correspondingly longer lifetimes, thus giving a large amplitude to on-shell 1e1a processes in region II. On the other hand, the 1e1a processes for electrons in region III lead to intermediate states in region II or above (see Fig. 6d), which have much shorter lifetimes than in region I, resulting in a smaller 1e1a rate in region III compared to region II. These arguments allow us to rationalize the 2ph process rates in Fig. 5.

Let us discuss the temperature dependence of the 2ph scattering processes, focusing on the ratio Γ^{(2ph)}∕ Γ^{(1ph)} of the 2ph scattering rates to the leading-order 1ph rates. We provide a brief analysis here and a more extensive discussion in the Supplementary Information. The temperature dependence of the 2ph rates originates from the intermediate state lifetimes in the denominators of Eq. (3) and the thermal factors *A*^{(i)} in Eq. (4). The intermediate state *e*–ph scattering rate increases with temperature due to an increase in phonon number, so the intermediate state lifetimes become shorter, lowering the 2ph rate Γ^{(2ph)} for increasing temperatures. On the contrary, the thermal factors *A*^{(i)}, which contain factors proportional to the phonon number *N*, increase rapidly with temperature, making Γ^{(2ph)} greater at higher temperatures. In the 200–400 K temperature range investigated here, these two effects compensate in GaAs, resulting in nearly temperature independent Γ^{(2ph)} ∕ Γ^{(1ph)} ratios. We conclude that the 2ph processes are equally as important relative to the leading-order 1ph processes over a wide temperature range near room temperature.

### Contribution to the electron mobility

Since the 2ph contributions are significant, one expects that they play a role in charge transport. Figure 7 shows the electron mobility in GaAs obtained by solving the BTE both within the relaxation time approximation (RTA) and with a more accurate iterative approach (ITA)^{4,7}, which we have extended here to include 2ph processes (see Supplementary Information). Results are given for calculations that either neglect or include the 2ph contributions. The iterative solution with only 1ph scattering overestimates the electron mobility in GaAs by 40–80% at 200–400 K, consistent with the results in ref. ^{8} (and in ref. ^{4}, provided the polar correction is not included for the acoustic modes, which artificially increases the acoustic mode contribution to scattering in ref. ^{4} and lowers the computed mobility). This result is puzzling, since the BTE can accurately predict the mobility in nonpolar semiconductors; the discrepancy with experiment is too large to be due to small errors in the electron effective mass (in our calculation, the effective mass is 0.055*m*_{0} versus an experimental value of 0.067*m*_{0}, where *m*_{0} is the electron mass). Redoing the mobility calculation with a band structure with the experimental effective mass of 0.067*m*_{0} (see Supplementary Fig. 3) reduces the mobility overestimate to roughly 20% compared to experiment for the iterative BTE solution with 1ph scattering, and including 2ph processes still improves the agreement with experiment of the mobility and its temperature dependence.

This open problem is solved here by including the 2ph processes in the ITA (see Fig. 7), which lowers the mobility due to the additional 2ph scattering and gives mobility values in excellent agreement with experiments^{32,33,34,35,36,37,38}. Both the absolute value of the mobility and its temperature dependence improve when the 2ph processes are included. This result implies that the agreement with experiment of the RTA with 1ph processes^{3} is due to error compensation. We conclude that the 2ph contributions are crucial to improving the computed electron mobility in GaAs.

Since the long-range LO mode coupling is dominant in many polar semiconductors and oxides, we expect the 2ph processes to be important in a wide range of polar materials. We apply our approach to BaSnO_{3}, a weakly polar oxide with dominant LO-mode *e*–ph interactions and a high room-temperature electron mobility. We compute the 2ph scattering rates and the mobility including 2ph processes in BaSnO_{3} (see Supplementary Fig. 2). Similar to GaAs, we find that the mobility computed with the iterative BTE including only 1ph processes is substantially higher than the experimental value, while adding 2ph scattering significantly improves the agreement with experiment. This result confirms that our approach is broadly applicable to weakly polar semiconductors and oxides, in which perturbation theory is valid and the mobility near room temperature is usually limited by scattering with the LO mode.

## Discussion

Our calculations show quantitatively that the 2ph contributions are substantial even in a weakly polar material such as GaAs in which higher-order *e*–ph interactions and polaron effects are typically neglected. The approach introduced in this work is general, and it is expected to give accurate *e*–ph scattering rates and carrier mobilities in weakly polar III-V and II-VI semiconductors and in polar oxides with high room-temperature mobility (and thus, weak polaron effects^{20}); the perturbative approach is valid in these materials, where the long-range *e*–ph interactions make next-to-leading-order effects substantial.

An open question left out for future work is whether the 2ph processes due to acoustic modes can be important in metals and nonpolar materials. A practical observation is that the mobility computed using only 1ph processes matches experiment closely in elemental metals^{2} and nonpolar semiconductors^{8}; therefore, 2ph acoustic processes are expected to be negligible, although a rigorous proof (or numerical evidence) of this point is more challenging. As discussed by Kocevar^{12}, there are two *e*–ph perturbations contributing to 2ph processes. The one considered in this work is due to the first derivative of the potential with respect to lattice vibrations, taken to second order in perturbation theory, which leads to 2ph scattering consisting of two consecutive 1ph scattering events. The second contribution, which is not included in this work, is due to the second derivatives of the potential with respect to lattice vibrations; it leads to a direct 2ph interaction associated with the so-called Debye-Waller (DW) vertex, which represents an electron interacting simultaneously with two phonons^{12}. The DW vertex leads to additional self-energy diagrams and 2ph processes, some of which are illustrated in Fig. 8.

Due to translational invariance, for acoustic phonons in the long-wavelength limit there is a strong cancellation between the scattering processes considered here and those due to the DW interaction^{12,39}. This cancellation has also been hypothesized to be valid for acoustic phonons beyond the long-wavelength limit, although this point has never been proven. This result, together with the numerical evidence that 1ph processes are sufficient to explain the mobility in nonpolar semiconductors and elemental metals, is suggestive of a negligible role of 2ph acoustic processes. Since the DW vertex cannot currently be computed ab initio for two general phonon momenta, numerical studies of 2ph processes due to the DW interaction are left out for future work. In the calculations on the two weakly polar materials considered here, GaAs and BaSnO_{3}, the LO mode interaction is dominant, and our results are nearly unchanged if acoustic phonons are included (see Supplementary Fig. 1).

In summary, our calculations of the 2ph scattering rates and their contribution to the mobility pave the way to studying higher-order *e*–ph interactions and charge transport in polar materials from first principles. Together with recently proposed methods to treat charge transport in materials with polarons and stronger *e*–ph interactions^{24}, it is clear that ab initio calculations are becoming able to investigate *e*–ph interactions and charge transport in a wide range of polar materials.

## Methods

### Numerical calculations

The DFT and lowest-order *e*–ph calculations in GaAs follow our previous work^{3}. Briefly, we carry out DFT calculations in GaAs using the Quantum ESPRESSO code^{40} with a plane-wave basis set. We employ the local density approximation^{41} and norm-conserving pseudopotentials^{42}. A relaxed lattice constant of 5.55 Å, a kinetic energy cutoff of 72 Ry and 8 × 8 × 8 **k**-point grids are used in all DFT calculations. Phonon dispersions are computed with density functional perturbation theory (DFPT)^{43} on an 8 × 8 × 8 **q**-point grid. The *e*–ph coupling constants, *g*_{nmν}(**k**, **q**), are computed using these coarse **k**- and **q**-point grids^{3} using DFPT together with our in-house developed Perturbo code^{31}, and interpolated using Wannier functions^{44} generated with Wannier90^{45}.

To compute and converge the 2ph scattering rates, we use Monte Carlo integration by sampling up to 3 billion random (**q**, **p**) pairs of Brillouin zone points drawn from a Cauchy distribution^{3}. The delta function in Eq. (2) is approximated by a Gaussian with a small broadening of 5 meV. Since in GaAs the LO phonon dominates *e*–ph scattering for electrons within ~100 meV of the conduction band edge^{3}, in the 2ph scattering rate calculations we use only the long-range coupling to the LO modes^{11,46,47}, which greatly reduces the 2ph computational cost while giving 2ph scattering rates nearly identical to calculations that include all phonon modes (see Supplementary Fig. 1). We use a carrier concentration of 10^{16} cm^{−3} in all scattering rate and mobility calculations, and set the chemical potential accordingly at each temperature. The mobility calculations are carried out by integrating over energies of up to 250 meV above the conduction band minimum; we have verified that this energy window is sufficient to converge the mobility.

We highlight that the 2ph scattering rate calculations are very computationally expensive due to the large number of points sampled in the BZ integral over two-phonon momenta. The average cost for computing the 2ph scattering rate in GaAs is ~500 CPU core hours per electronic state for each iteration over the intermediate state lifetime; a typical calculation requires ten iterations to converge, for a total cost of 5000 CPU core hours per electronic state. By contrast, it takes only 0.08 CPU core hours (5 min on a single CPU core) per electronic state to converge the 1ph scattering rate using our code. The computational cost of the 2ph scattering rates is thus approximately 10^{4}−10^{5} higher than the 1ph scattering rates.

## Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## Code availability

The code employed in this work, Perturbo, is available for download at https://perturbo-code.github.io/. The 2ph routines developed in this work are available from the corresponding author upon reasonable request.

## References

- 1.
Bernardi, M. First-principles dynamics of electrons and phonons.

*Eur. Phys. J. B***89**, 239 (2016). - 2.
Mustafa, J. I., Bernardi, M., Neaton, J. B. & Louie, S. G. Ab initio electronic relaxation times and transport in noble metals.

*Phys. Rev. B***94**, 155105 (2016). - 3.
Zhou, J.-J. & Bernardi, M. Ab initio electron mobility and polar phonon scattering in GaAs.

*Phys. Rev. B***94**, 201201(R) (2016). - 4.
Liu, T.-H., Zhou, J., Liao, B., Singh, D. J. & Chen, G. First-principles mode-by-mode analysis for electron-phonon scattering channels and mean free path spectra in GaAs.

*Phys. Rev. B***95**, 075206 (2017). - 5.
Zhou, J.-J., Hellman, O. & Bernardi, M. Electron-phonon scattering in the presence of soft modes and electron mobility in SrTiO

_{3}perovskite from first principles.*Phys. Rev. Lett.***121**, 226603 (2018). - 6.
Lee, N.-E., Zhou, J.-J., Agapito, L. A. & Bernardi, M. Charge transport in organic molecular semiconductors from first principles: the bandlike hole mobility in a naphthalene crystal.

*Phys. Rev. B***97**, 115203 (2018). - 7.
Li, W. Electrical transport limited by electron-phonon coupling from Boltzmann transport equation: an ab initio study of Si, Al, and MoS

_{2}.*Phys. Rev. B***92**, 075405 (2015). - 8.
Ma, J., Nissimagoudar, A. S. & Li, W. First-principles study of electron and hole mobilities of Si and GaAs.

*Phys. Rev. B***97**, 045201 (2018). - 9.
Sohier, T., Campi, D., Marzari, N. & Gibertini, M. Mobility of two-dimensional materials from first principles in an accurate and automated framework.

*Phys. Rev. Mater.***2**, 114010 (2018). - 10.
Fröhlich, H. Electrons in lattice fields.

*Adv. Phys.***3**, 325–361 (1954). - 11.
Vogl, P. Microscopic theory of electron-phonon interaction in insulators or semiconductors.

*Phys. Rev. B***13**, 694–704 (1976). - 12.
Kocevar, P.

*Multiphonon Scattering*167–174 (Springer US, 1980). - 13.
Alldredge, G. P. & Blatt, F. J. On the role of two-phonon processes in the energy relaxation of a heated-electron distribution.

*Ann. Phys.***45**, 191–231 (1967). - 14.
Sher, A. & Thornber, K. K. Resonant electron-phonon scattering in polar semiconductors.

*Appl. Phys. Lett.***11**, 3–5 (1967). - 15.
Thorbergsson, G. I. & Sak, J. Mobility of an acoustic polaron at very low temperatures.

*Phys. Lett. A***116**, 325–328 (1986). - 16.
Mora-Ramos, M. E., Rodríguez, F. J. & Quiroga, L. Polaron properties of III-V nitride compounds: second-order effects.

*J. Phys. Condens. Matter***11**, 8223–8238 (1999). - 17.
Smondyrev, M. A. Diagrams in the polaron model.

*Theor. Math. Phys.***68**, 653–664 (1986). - 18.
Wang, Z. & Mahan, G. D. Low-temperature resistivity from electron-dual-phonon processes for alkali metals.

*Phys. Rev. B***39**, 10753–10760 (1989). - 19.
Woods, L. M. & Mahan, G. D. Nonlinear electron-phonon heat exchange.

*Phys. Rev. B***57**, 7679–7685 (1998). - 20.
Emin, D.

*Polarons*(Cambridge University Press, 2012). - 21.
Gunnarsson, O., Meden, V. & Schönhammer, K. Corrections to Migdal’s theorem for spectral functions: a cumulant treatment of the time-dependent Green’s function.

*Phys. Rev. B***50**, 10462–10473 (1994). - 22.
Story, S. M., Kas, J. J., Vila, F. D., Verstraete, M. J. & Rehr, J. J. Cumulant expansion for phonon contributions to the electron spectral function.

*Phys. Rev. B***90**, 195135 (2014). - 23.
Nery, J. P. et al. Quasiparticles and phonon satellites in spectral functions of semiconductors and insulators: cumulants applied to the full first-principles theory and the Fröhlich polaron.

*Phys. Rev. B***97**, 115145 (2018). - 24.
Zhou, J.-J. & Bernardi, M. Predicting charge transport in the presence of polarons: the beyond-quasiparticle regime in SrTiO

_{3}.*Phys. Rev. Res.***1**, 033138 (2019). - 25.
Noffsinger, J., Kioupakis, E., Van de Walle, C. G., Louie, S. G. & Cohen, M. L. Phonon-assisted optical absorption in silicon from first principles.

*Phys. Rev. Lett.***108**, 167402 (2012). - 26.
Murayama, M. & Nakayama, T. Ab initio calculations of two-photon absorption spectra in semiconductors.

*Phys. Rev. B***52**, 4986–4997 (1995). - 27.
Hayat, A., Ginzburg, P. & Orenstein, M. Observation of two-photon emission from semiconductors.

*Nat. Photonics***2**, 238 (2008). - 28.
Bern, Z., Dixon, L. & Ghinculov, A. Two-loop correction to Bhabha scattering.

*Phys. Rev. D***63**, 053007 (2001). - 29.
Aoyama, T., Kinoshita, T. & Nio, M. Revised and improved value of the QED tenth-order electron anomalous magnetic moment.

*Phys. Rev. D***97**, 036001 (2018). - 30.
Mahan, G. D.

*Many-Particle Physics*, 3rd edn (Springer, 2000). - 31.
Zhou, J.-J. et al. Perturbo: a software package for

*ab initio*electron-phonon interactions, charge transport and ultrafast dynamics. Preprint at https://arxiv.org/abs/2002.02045 (2020). - 32.
Rode, D. L. Electron mobility in direct-gap polar semiconductors.

*Phys. Rev. B***2**, 1012–1024 (1970). - 33.
Rode, D. L. & Knight, S. Electron transport in GaAs.

*Phys. Rev. B***3**, 2534–2541 (1971). - 34.
Blakemore, J. S. Semiconducting and other major properties of gallium arsenide.

*J. Appl. Phys.***53**, R123–R181 (1982). - 35.
Hicks, H. G. B. & Manley, D. F. High purity GaAs by liquid phase epitaxy.

*Solid State Commun.***7**, 1463–1465 (1969). - 36.
Wolfe, C. M., Stillman, G. E. & Lindley, W. T. Electron mobility in high-purity GaAs.

*J. Appl. Phys.***41**, 3088–3091 (1970). - 37.
Blood, P. Electrical properties of

*n*-type epitaxial GaAs at high temperatures.*Phys. Rev. B***6**, 2257–2261 (1972). - 38.
Nichols, K. H., Yee, C. M. L. & Wolfe, C. M. High-temperature carrier transport in n-type epitaxial GaAs.

*Solid-State Electron.***23**, 109–116 (1980). - 39.
Holstein, T. Theory of ultrasonic absorption in metals: the collision-drag effect.

*Phys. Rev.***113**, 479–496 (1959). - 40.
Giannozzi, P. et al. Quantum ESPRESSO: a modular and open-source software project for quantum simulations of materials.

*J. Phys. Condens. Matter***21**, 395502 (2009). - 41.
Perdew, J. P. & Zunger, A. Self-interaction correction to density-functional approximations for many-electron systems.

*Phys. Rev. B***23**, 5048–5079 (1981). - 42.
Troullier, N. & Martins, J. L. Efficient pseudopotentials for plane-wave calculations.

*Phys. Rev. B***43**, 1993–2006 (1991). - 43.
Baroni, S., de Gironcoli, S., DalCorso, A. & Giannozzi, P. Phonons and related crystal properties from density-functional perturbation theory.

*Rev. Mod. Phys.***73**, 515–562 (2001). - 44.
Giustino, F., Cohen, M. L. & Louie, S. G. Electron-phonon interaction using Wannier functions.

*Phys. Rev. B***76**, 165108 (2007). - 45.
Mostofi, A. A. et al. An updated version of wannier90: A tool for obtaining maximally-localised Wannier functions.

*Comput. Phys. Commun.***185**, 2309–2310 (2014). - 46.
Sjakste, J., Vast, N., Calandra, M. & Mauri, F. Wannier interpolation of the electron-phonon matrix elements in polar semiconductors: polar-optical coupling in GaAs.

*Phys. Rev. B***92**, 054307 (2015). - 47.
Verdi, C. & Giustino, F. Fröhlich electron-phonon vertex from first principles.

*Phys. Rev. Lett.***115**, 176401 (2015).

## Acknowledgements

This work was supported by the Air Force Office of Scientific Research through the Young Investigator Program Grant FA9550-18-1-0280. J.-J.Z. and H.-Y.C. were supported by the National Science Foundation under Grant No. ACI- 1642443, which provided for code development. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

## Author information

### Affiliations

### Contributions

M.B. conceived the research. N.-E.L. derived the analytic expressions, implemented the code and carried out the numerical calculations. J.-J.Z. helped with code development. H.-Y.C. helped deriving the equations. N.-E.L. and M.B. wrote the manuscript. All authors contributed to analyzing the results and editing the manuscript.

### Corresponding author

## Ethics declarations

### Competing interests

The authors declare no competing interests.

## Additional information

**Peer review information** *Nature Communications* thanks the anonymous reviewer(s) for their contribution to the peer review of this work.

**Publisher’s note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Supplementary information

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

## About this article

### Cite this article

Lee, NE., Zhou, JJ., Chen, HY. *et al.* Ab initio electron-two-phonon scattering in GaAs from next-to-leading order perturbation theory.
*Nat Commun* **11, **1607 (2020). https://doi.org/10.1038/s41467-020-15339-0

Received:

Accepted:

Published:

## Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.