Introduction

Organic molecular crystals (OMCs) possess complex crystal structures with an intricate interplay of electronic and structural degrees of freedom. Their electronic properties range from metallic to semiconducting or insulating, and they can host ferroelectricity1, magnetism2, and superconductivity3. OMCs are also versatile materials with broad applications, for example in electronics4, light-emitting diodes5,6, spintronics7,8, batteries9,10, and solar cells11,12,13. The charge carrier mobility is a key figure of merit for organic materials in these devices14,15,16,17. Yet, understanding charge transport and accurately predicting the mobility in OMCs remain open challenges. The mobility varies significantly across OMCs due to their wide range of carrier effective masses and electron-phonon (e-ph) interaction strengths16,18. Even in the same organic crystal, electron and hole carriers can exhibit different transport regimes, and the mobility can vary by orders of magnitude for different crystallographic directions.

Charge transport in OMCs is often classified into two limiting cases − the band transport and polaron hopping regimes, each entailing specific transport mechanisms19. In band transport, charge carriers are delocalized, the e-ph coupling is weak, and the mobility is correspondingly high (>10 cm2V−1s−1) and characterized by a power-law decrease with temperature. Band transport in OMCs is usually governed by scattering of carriers with low-energy acoustic and intermolecular phonons, and the corresponding e-ph interactions are often modeled with the Peierls Hamiltonian20. These weak e-ph interactions can be described with lowest-order perturbation theory, and the OMC mobility can be accurately predicted using the Boltzmann transport equation (BTE)18,21,22,23, including all phonon modes on the same footing.

In the polaron hopping regime, the charge carriers interact strongly with phonons, forming self-localized (small) polarons, which are often modeled with the Holstein Hamiltonian to describe intramolecular e-ph interactions24. The resulting charge transport is dominated by thermally activated polaron hopping and is often described with Marcus theory25,26. The mobility in the polaron hopping regime is relatively small, usually below 0.1 cm2V−1s−1, and is more challenging to predict from first principles.

Between these two limiting scenarios, OMCs also exhibit an intermediate transport regime, for which neither the band transport nor the polaron hopping pictures are fully adequate19,26. In the intermediate regime, the mobility exhibits a bandlike power-law temperature dependence26,27,28, but polarons can still be present and low mobility values (<1 cm2V−1s−1) are common19. A signature of intermediate transport is the violation of the Mott-Ioffe-Regel limit29, whereby the carrier mean-free paths become smaller than the intermolecular distance26, making the BTE description inadequate.

Various approaches have been employed to study intermediate transport in OMCs; they typically employ a Holstein Hamiltonian or a (Peierls-type) dynamical disorder Hamiltonian, or a combination of both, and obtain the mobility via linear-response theory30,31,32,33, diffusion simulation34,35, surface hopping method36, or transient localization calculation37,38. These methods are highly valuable for studies of OMCs, although they usually rely on simplifying assumptions such as including only specific phonon modes and e-ph interactions, or fitting model parameters to experiments. To date, first-principles approaches to predict charge transport in the intermediate regime with quantitative accuracy are scarce, especially within rigorous treatments based on many-body perturbation theory.

In this work, we develop rigorous calculations of the mobility in the intermediate charge transport regime in OMCs. Focusing on naphthalene crystal as a case study, we employ a finite-temperature cumulant approach27 to capture the strong e-ph interactions and polaron effects characteristic of the intermediate regime, and employ Green-Kubo theory to compute the electron mobility. All phonon modes are included and treated on equal footing. This cumulant plus Kubo (CK) approach is shown to predict the electron mobility in the intermediate regime with a high accuracy, within a factor of two of experiments between 100 and 300 K for crystallographic directions parallel to the naphthalene molecular planes. We additionally show the failure of the BTE to describe the mobility in the intermediate regime.

Our analysis of the electronic spectral functions reveals the presence of a broad satellite next to the quasiparticle (QP) peak, explaining the breakdown of the BTE and the band transport picture. The broadening of the QP peak is mainly due to coupling with low-energy intermolecular phonons, while the polaron satellite peaks are due to strong coupling with intramolecular phonons. Therefore, the intermolecular modes are directly responsible for scattering the electrons and limiting the mobility, whereas the intramolecular modes limit the mobility indirectly, by transferring spectral weight from the QP peak to the polaron satellites. For charge transport normal to the molecular planes, we find that both the BTE and CK approaches cannot correctly predict the mobility, which experimentally is nearly temperature independent and governed by small-polaron hopping39,40,41. This finding restricts the applicability of the CK method to intermediate e-ph coupling strengths.

Taken together, our work provides an accurate first-principles method to study polaron transport in OMCs, and unravels the interplay of low- and high-energy phonon modes in the intermediate regime. Our results provide a blueprint for studying charge transport in a wide range of organic crystals.

Results

Computational approach

We compute the ground state electronic structure of naphthalene crystal using plane-wave density functional theory (DFT) calculations with the Quantum Espresso code42,43. We employ the generalized gradient approximation44 and norm-conserving pseudopotentials45 from Pseudo Dojo46. The DFT band structure is refined using GW calculations (with the Yambo code47,48) to better capture dynamical screening effects. Maximally localized Wannier functions49 are generated with the Wannier90 code50 following a procedure similar to Lee et al.18. We compute the e-ph interactions and charge transport separately at four temperatures (100, 160, 220, and 300 K), using different experimental lattice parameters at each temperature51 and relaxing the atomic positions with DFT. We obtain the lattice dynamics and e-ph perturbation potentials from density functional perturbation theory (DFPT)52, and compute the e-ph interactions with the Perturbo code53. Additional numerical details are provided in the Methods section.

Using the computed e-ph interactions, we study charge transport in the BTE and CK frameworks with the Perturbo code53. In the BTE, the mobility tensor μαβ is computed in the relaxation time approximation (RTA):

$$\begin{array}{l}{\mu }_{\alpha \beta }(T)=\frac{2e}{{n}_{\rm{c}}{V}_{\rm uc}}\int dE\left(-\frac{\partial f(E,T)}{\partial E}\right)\\ \qquad\qquad\,\,\,\,\,\times\, \mathop{\sum}\limits_{n{{{\bf{k}}}}}{\tau }_{n{{{\bf{k}}}}}(T){v}_{n{{{\bf{k}}}}}^{\alpha }{v}_{n{{{\bf{k}}}}}^{\beta }\delta (E-{\varepsilon }_{n{{{\bf{k}}}}}),\end{array}$$
(1)

where α and β are Cartesian directions parallel to the crystal principal axes, T is the temperature, e the electronic charge, nc the carrier concentration, Vuc the unit cell volume, f the electronic Fermi-Dirac distribution, and E is the electron energy. Here and below, n is the band index and k the crystal momentum of the electronic states. The BTE mobility depends on the electron band energies εnk, the corresponding band velocities vnk, and the state-dependent e-ph relaxation times τnk obtained within lowest-order perturbation theory53,54. As a sanity check, we compute the mobility at 220 K by solving the full linearized BTE with an iterative approach (ITA)53,55, and find that in naphthalene it gives results identical to the RTA, justifying our use of the RTA.

To properly treat strong e-ph interactions and include polaron effects in the mobility, we employ a finite-temperature cumulant approach in which the retarded electron Green’s function \({G}_{n{{{\bf{k}}}}}^{{{{\rm{R}}}}}\) is written using the exponential ansatz27,56,57,58,59,60

$$\begin{array}{r}{G}_{n{{{\bf{k}}}}}^{{{{\rm{R}}}}}(t,T)={G}_{n{{{\bf{k}}}}}^{{{{\rm{R,0}}}}}(t){{{{\rm{e}}}}}^{{C}_{n{{{\bf{k}}}}}(t,T)},\end{array}$$
(2)

where \({G}_{n{{{\bf{k}}}}}^{{{{\rm{R,0}}}}}\) is the non-interacting Green’s function and Cnk is the cumulant function, obtained here at finite temperatures from the lowest-order e-ph self-energy (see Methods section). The electron spectral function is obtained from the Green’s function at each electron energy E using

$${A}_{n{{{\bf{k}}}}}(E,T)=-{{{\rm{Im}}}}{G}_{n{{{\bf{k}}}}}^{{{{\rm{R}}}}}(E,T)/{{{\rm{\pi }}}}.$$
(3)

In the CK method, the mobility tensor is computed directly from the spectral function using the linear-response Green-Kubo formula27,57,61:

$$\begin{array}{r}{\mu }_{\alpha \beta }(T)=\frac{1}{{n}_{\rm{c}}e}\int \,dE\,{{{{\Phi }}}}_{\alpha \beta }(E,T),\end{array}$$
(4)

where the integrand is the transport distribution function (TDF). Under the approximation of neglecting vertex corrections, the TDF reads61

$${{{{\Phi }}}}_{\alpha \beta }(E,T)=\frac{\pi \hslash {e}^{2}}{{V}_{{\rm uc}}}\mathop{\sum}\limits_{n{{{\bf{k}}}}}{v}_{n{{{\bf{k}}}}}^{\alpha }{v}_{n{{{\bf{k}}}}}^{\beta }| {A}_{n{{{\bf{k}}}}}(E,T){| }^{2}\left(-\frac{\partial f(E,T)}{\partial E}\right)\,,$$
(5)

where vnk are the unperturbed electron band velocities, the same as those used in Eq. (1)61. The CK mobility defined in Eq. (4) is obtained from the cumulant spectral function, therefore it takes into account the strong e-ph coupling and polaron effects. The CK calculations have been shown to provide results in close agreement with the BTE-RTA in the limit of weak e-ph interactions (see Zhou and Bernardi27 for a calculation on GaAs).

Electron mobility

The crystal structure of naphthalene consists of molecular planes in the a and b crystallographic directions, stacked along the plane-normal c* direction (Fig. 1a). We first discuss charge transport in the molecular planes. For hole carriers, we have previously shown that such in-plane transport is bandlike and well described by the BTE18. In this work, we focus on the mobility of the electron carriers, which due to their flatter electronic bands with greater effective masses compared to holes (Fig. 1b) are expected to exhibit lower mobilities and a range of transport regimes. Only the electronic bands formed by the lowest unoccupied molecular orbital (LUMO) and the next-higher-energy orbital (LUMO+1) contribute to electron transport in the 100−300 K temperature range studied here, so we consider only these two bands in our mobility calculations.

Fig. 1: Structure.
figure 1

a Monoclinic crystal structure of naphthalene, with molecular a-b planes stacked in the plane-normal c* direction. b Band structure of naphthalene, showing the two lowest-energy electron (LUMO and LUMO+1) and hole (HOMO and HOMO-1) bands.

Figure 2 shows the in-plane electron mobilities computed with the BTE and CK methods, and compares them with experimental data39. We fit each mobility curve with a T−n power-law temperature trend and give the exponent n next to each curve. The results show that the BTE predicts a much stronger temperature dependence of the mobility than in experiment, with errors in the computed exponents for transport along the a and b crystallographic directions (mobilities μa and μb in Fig. 2, respectively) of over 100% for μa and 270% for μb relative to the exponent n obtained by fitting the experimental results. Due to this error, the BTE greatly overestimates the mobility at low temperatures − for example, μa at 100 K from the BTE is an order of magnitude greater than the experimental value.

Fig. 2: In-plane electron mobility in naphthalene.
figure 2

a Electron mobility in direction a. b Electron mobility in direction b. Results obtained from BTE and CK calculations are compared with experimental data from Madelung et al.39.

These results are a strong evidence of the failure of the Boltzmann equation to describe electron transport in naphthalene; the physical origin of this failure is examined below. Note that the BTE failure is not a consequence of our use of the RTA, as the full solution of the BTE53 gives results nearly identical to the RTA (see the ITA points at 220 K in Fig. 2). The fact that the mobility has a power-law temperature dependence but is not correctly predicted by the BTE is a hallmark of the intermediate transport regime26,28.

The CK calculations give significantly improved results (Fig. 2). The CK mobilities are within a factor of 2 of experiment for μa and 1.3 for μb in the entire 100−300 K temperature range. The error in the T−n exponent is reduced to 20% for μa and 45% for μb relative to experiment, a five-fold improvement in accuracy over the BTE results. Achieving this level of accuracy for quantitative predictions of the mobility in OMCs has recently become possible in the band transport regime18 but has so far remained challenging in the intermediate regime. As we discuss below, by combining the cumulant and Green-Kubo frameworks, our CK approach can capture key polaron effects in the intermediate regime such as higher-order e-ph coupling and spectral weight transfer, resulting in improved mobility predictions.

Electron spectral function

The electron spectral function is central to understanding polaron effects27 and intermediate charge transport. The spectral function can be viewed as the density of states of a single electronic state, and it integrates to one over energy due to a well-known sum rule57. In Fig. 3a, we show the spectral function at three temperatures, using results obtained with our cumulant method for the electronic state at the conduction band minimum (CBM) (Γ point in Fig. 1b). At 100 K, next to the main QP peak we find a broad spectral feature associated with the combined excitation of an electron QP plus one or two phonons. This broad satellite combines contributions from multiple satellite peaks, as shown by the arrows in Fig. 3a, and is a signature of polaron formation27. At higher temperatures, the QP and satellite peaks broaden and ultimately merge into a continuum at 300 K. The coexistence of a well-formed QP peak and broad satellites shows that large-polaron effects, characteristic of e-ph interactions with intermediate strength, are a key characteristic of the intermediate transport regime.

Fig. 3: Spectral function and e-ph coupling.
figure 3

a Spectral functions computed at three temperatures for the CBM electronic state. Results from the cumulant approach (\({A}_{{{{\rm{CBM}}}}}^{{{{\rm{C}}}}}\)) are compared to the Dyson-Migdal spectral function (\({A}_{{{{\rm{CBM}}}}}^{{{{\rm{DM}}}}}\)). The QP peak is chosen as the zero of the energy axis for each spectral function. The transport distribution function (TDF) in arbitrary units is also shown at each temperature. b Electron spectral function for the LUMO and LUMO+1 bands along a high-symmetry path, computed at 100 K using the cumulant method. The solid line is the GW band structure and the dashed line shows the renormalized cumulant band structure obtained by connecting the QP peaks of the spectral functions. c Gauge-invariant e-ph coupling strength as a function of phonon energy. The energies ω1 and ω2 of the two phonon modes with strongest e-ph coupling are shown with vertical dashed lines. d Atomic displacements for the two intramolecular modes with the strongest e-ph coupling.

The cumulant spectral functions for multiple electronic states in the LUMO and LUMO+1 bands can be combined to obtain a polaron band structure renormalized by the e-ph interactions. Figure 3b compares the band structures at 100 K computed with the GW method and with our cumulant calculations that use the GW band structure as input. The cumulant band structure, obtained by connecting the QP peaks of the cumulant spectral functions at neighboring k-points, captures polaron effects such as QP mass and weight renormalization. At 100 K, where the QP peaks are well-defined, we calculate the renormalized effective masses from the cumulant band structure, and find a moderate effective mass enhancement of 15−35% compared to the effective masses in the GW band structure in the in-plane directions. This shows that the cumulant approach can capture the bandwidth-narrowing due to strong e-ph interactions and polaron effects in OMCs62.

The physical origin of the polaron satellite in Fig. 3a is of key importance. In the prototypical case of a polar inorganic material with strong e-ph coupling with longitudinal optical (LO) phonons, the satellite peaks are located at the LO-mode energy ωLO (and its multiples) relative to the QP peak27,63. Here, due to the presence of a large number of phonon modes in OMCs (108 in napthalene), the satellites merge into a broad spectral feature resembling a long tail of the QP peak, with contributions from various phonon modes. To explain the origin of this broad satellite, in Fig. 3c we analyze the e-ph coupling strength for an electronic state near the CBM, as quantified by the absolute value of the gauge-invariant e-ph coupling, ∣g∣ (see Methods).

In naphthalene, the 12 lowest-energy phonon modes are intermolecular, and the remaining 96 modes are intramolecular vibrations18. In Fig. 3a, the mode with the strongest e-ph coupling, an intramolecular phonon with energy ω1 ≈ 0.1 eV, generates a satellite peak in the spectral function at energy ω1 relative to the QP peak. The intramolecular phonon with the second strongest e-ph coupling, with energy ω2 ≈ 0.2 eV, gives a second contribution to the broad satellite, followed by a plateau at higher energy. Finally, the inflection point in the spectral function at energy ω1 + ω2 is due to higher-order e-ph coupling from the two modes with energies ω1 and ω2. The atomic displacements associated with these two intramolecular modes are shown in Fig. 3d. Both modes involve vibrations of the hydrogen atoms, in one case in the carbon ring plane and in the other case normal to the carbon rings. Our analysis demonstrates that these intramolecular phonons are responsible for the formation of polarons in naphthalene. This strong coupling with intramolecular phonons and the associated satellite peak in the spectral function are consistent with recent results from the Holstein–Peierls model32.

Interestingly, lowest-order theory is wholly inadequate to describe this polaronic regime with intermediate e-ph coupling strength. To show this point, we compute the Dyson-Migdal (DM) spectral function [see Equation (11)], which is obtained from the lowest-order e-ph self-energy and therefore does not include polaron effects. From the comparison of the cumulant and DM spectral functions in Fig. 3a it is clear that the DM spectral functions have a Lorentzian shape and lack any satellite structure. As a result, the subtle interplay between inter- and intramolecular phonons in the QP and satellite peaks cannot be captured in lowest-order theory. As we discuss below, this is the origin of the failure of the BTE to describe transport in the intermediate regime.

Failure of the Boltzmann equation

It is important to understand the microscopic origin of the failure of the BTE, and the success of the CK method, to describe transport in the intermediate regime. In the Green-Kubo framework, the mobility is given by an integral over electron energies [see Eq. (4)], which in principle combines contributions from all features of the spectral function. To quantify the contributions of the QP and satellite peaks to charge transport, we analyze the mobility integrand, the TDF in Eq. (5), and plot it together with the spectral functions in Fig. 3a. We find that the TDF decays rapidly outside the QP peak, within an energy ω1 of the QP peak at low temperature and ω2 at 300 K. Therefore any spectral function feature with energy greater than ω2 does not overlap with the TDF and cannot contribute to charge transport between 100 and 300 K. In this temperature range, although the mobility is mainly governed by the QP peak, polaron effects still contribute in important ways.

First, due to higher-order e-ph coupling with both inter- and intramolecular phonons, the linewidth of the QP peak in the cumulant spectral function is different than in the DM spectral function from lowest-order theory (see Fig. 4), whose linewidth is the scattering rate entering the BTE mobility calculation. This QP linewidth discrepancy is temperature and energy dependent (Fig. 4), which explains why the BTE cannot correctly predict the value and temperature dependence of the mobility in the intermediate regime, corroborating our results in Fig. 2. Second, the broad satellite in the cumulant spectral function limits the carrier mobility indirectly, by transferring spectral weight away from the QP peak (recall that the spectral function integrates to one over energy). The satellite peak at ω1 contributes directly to transport only above ~200 K, where the QP peak broadens, merging with the satellite, and overlapping with the TDF.

Fig. 4: Quasiparticle peak broadening.
figure 4

Full width at half maximum of the QP peaks, shown at two temperatures for both the cumulant and DM spectral functions. The energy zero is set to the conduction band minimum.

The picture that emerges is that electron transport in the naphthalene molecular planes is mainly governed by the scattering of QPs with renormalized weight, which couple directly (via the main QP peak) with low-energy intermolecular phonons and indirectly (via weight transfer to the satellites) with higher-energy intramolecular phonons. The latter can also contribute directly to charge transport as the temperatures increases above 200 K. The ability of our CK approach to address these subtle e-ph interactions enables accurate predictions of the mobility and its temperature dependence in the intermediate regime, where the BTE with lowest-order e-ph coupling fails to capture these essential polaron effects.

Comparison with hole mobilities

Finally, we present mobility results for hole carriers in naphthalene to contrast their behavior with electron carriers. We compute the hole mobility in naphthalene between 100 and 300 K using the CK approach, and compare the results to BTE calculations and experiments (see Fig. 5). The BTE calculations are a refinement of those we presented in Lee et al.18, obtained here using a more accurate Wigner-Seitz cell summation procedure, as implemented in PERTURBO and described in detail in Zhou et al.53. The revised BTE mobilities follow an identical temperature trend as in our previous results in Lee et al.18, but their value is now greater than experiment, a physically meaningful trend for phonon-limited mobilities.

Fig. 5: Hole carrier mobilities in naphthalene.
figure 5

(μh) shown for the three crystallographic directions, a direction a, b direction b, and c direction c. The mobilities are computed both with the BTE in the relaxation time approximation (gray solid line) and with the CK approach (blue, green, and red). The computational settings are the same as in the electron mobility calculations. The BTE results presented in Lee et al.18, which were obtained using a previous implementation of e-ph interactions with less accurate Wigner-Seitz cell summations, are shown as dashed lines for comparison. The exponent n from a T−n power-law fit of the temperature dependence is given next to each curve. The experimental data (black) are from Madelung et al.39.

For hole carriers, both the BTE and CK methods give accurate predictions of the hole mobility, within a factor of 2−3 of experiment at all temperatures. The temperature dependence is nearly identical for the CK and the BTE mobilities, as shown by fitting the mobility curves with a T−n power-law and giving the exponent n next to each curve in Fig. 5. These findings demonstrate that the band transport picture of the BTE, which is inadequate for electron carriers, is sufficient to describe transport for hole carriers due to their more dispersive bands (see Fig. 1) and overall weaker e-ph coupling (see Supplementary Fig. 1).

Discussion

In naphthalene, measurements of the mobility in the direction normal to the molecular planes (c* direction in Fig. 1) point to a transport regime different from the in-plane directions. In experiments, the mobility along c* is <1 cm2V−1s−1 and is nearly temperature independent between 100 and 300 K39,40,41. These trends suggest that charge transport normal to the molecular planes may occur in the small-polaron hopping regime, where the carriers are strongly localized and the e-ph interaction is so strong that a diagram-resummation technique such as the cumulant method is not expected to give accurate results.

We calculate the plane-normal mobility using both the BTE and CK methods, and compare the results with experiments in Fig. 6. The computed mobility decreases with temperature in both the CK and BTE approaches, deviating substantially from the nearly temperature independent mobility found in experiment. It is encouraging that the CK mobility agrees well with experiment at 100 K and its temperature dependence is weaker than in the BTE − fitting the temperature dependence with a T−n power law gives an exponent n = 1.72 in the CK and n = 3.78 in the BTE method, versus n = 0.04 in experiment. However, although the CK provides a significant improvement over the BTE, it is clear that neither method can accurately describe charge transport normal to the molecular planes.

Fig. 6: Electron mobility in the plane-normal direction in naphthalene.
figure 6

The plot compares CK and BTE calculations with experimental data from Madelung et al.39.

The electron bands in the GW band structure are nearly flat in the plane-normal c*-direction (Γ − Z direction in Fig. 3b), with large effective masses of order 15 me for the GW calculation done on the 100 K structure, and greater at higher temperatures. Combined with the absence of a power-law temperature trend in the experimental mobility, this relatively flat band suggests that electrons are nearly localized to a single molecular plane and that transport in the plane-normal direction occurs via small-polaron hopping. The failure of the CK approach in this regime highlights the need for predictive first-principles approaches to study charge transport in the small-polaron hopping regime in OMCs.

In summary, we studied the electron mobility in naphthalene crystal as a paradigmatic case of intermediate charge transport in OMCs. Combining a finite-temperature cumulant method and Green-Kubo transport calculations, we demonstrated accurate predictions of the electron mobility and its temperature dependence in the intermediate regime. Our results reveal a subtle interplay between inter- and intramolecular phonons: the low-energy intermolecular phonons determine the broadening of the QP peak, while the intramolecular phonons are responsible for forming polarons and the associated satellite peaks. Both types of phonons contribute to limit the mobility. The broad satellite removes spectral weight from the QP peak, modifying the mobility and its temperature dependence. By capturing these subtle polaron effects, our CK approach addresses the shortcomings of the BTE for modeling the intermediate transport regime. We also highlighted the limitations of the CK approach to describe polaron-hopping between molecular planes. Taken together, our work advances microscopic understanding of the intermediate transport regime and paves the way for accurate first-principles calculations of the carrier mobility in OMCs.

Methods

First-principles calculations

We carry out first-principles density functional theory (DFT) calculations of the ground state and electronic structure of naphthalene using the Quantum Espresso code42,43. Thermal expansion of the lattice is taken into account by employing lattice constants51 taken from experiments at four different temperatures of 100, 160, 220, and 300 K. All calculations are carried out separately at these four temperatures. The initial atomic positions are also taken from experiment64,65. We use a kinetic energy cutoff of 90 Ry together with 2 × 4 × 2 and 4 × 4 × 4 k-point grids for self-consistent and non-self-consistent calculations, respectively. The Grimme van der Waals correction66,67 is included during structural relaxation of the atomic positions. To improve the description of dynamical electronic correlations, we correct the DFT electronic band structure with G0W0 calculations, which include 500 bands in the polarization function and a cutoff of 10 Ry in the dielectric screening using the Yambo code48. The Wannier90 code50 is employed to obtain Wannier functions and the corresponding transformation matrices, using the selected-columns-of-the-density-matrix method68. The lattice dynamics and e-ph perturbation potential are computed with density functional perturbation theory (DFPT)52 calculations on a 2 × 4 × 2 q-point grid (here and below, q is the phonon wavevector). Using our Perturbo code53, the electron and phonon data are combined to form the e-ph coupling matrix elements53:

$${g}_{mn\nu }({{{\bf{k}}}},{{{\bf{q}}}})=\sqrt{\frac{\hslash }{2{\omega }_{\nu {{{\bf{q}}}}}}}\mathop{\sum}\limits_{\kappa \alpha }\frac{{{{{\bf{e}}}}}_{\nu {{{\bf{q}}}}}^{\kappa \alpha }}{\sqrt{{M}_{\kappa }}}\left\langle m{{{\bf{k}}}}+{{{\bf{q}}}}\right|{\partial }_{{{{\bf{q}}}}\kappa \alpha }V\left|n{{{\bf{k}}}}\right\rangle ,$$
(6)

where \(\left|n{{{\bf{k}}}}\right\rangle\) are electronic Bloch states, ωνq are phonon energies, ∂qκαV are e-ph perturbation potentials, \({{{{\bf{e}}}}}_{\nu {{{\bf{q}}}}}^{\kappa \alpha }\) are phonon displacement vectors, and Mκ is the mass of atom κ. The absolute value of the gauge-invariant e-ph coupling strength shown in Fig. 3b is computed for each phonon mode ν and phonon wavevector q as53

$$| {g}_{\nu}({{{\bf{k}}}}=0,{{{\bf{q}}}})| =\sqrt{\mathop{\sum}\limits_{mn}| {g}_{mn\nu}({{{\bf{k}}}}=0,{{{\bf{q}}}}){|}^{2}/{N}_{{{{\rm{b}}}}}},$$
(7)

where Nb is the number of selected bands. The mobility calculations use a fine k-grid of 60 × 60 × 60 for the BTE and 30 × 30 × 30 for the CK method. Both methods use between 105 and 106 randomly selected q-points.

Electron-phonon scattering rate

The relaxation time τnk used in the BTE is computed as the inverse of the scattering rate, defined as54

$$\begin{array}{l}{{{{\Gamma }}}}_{n{{{\bf{k}}}}}(T)=\frac{2\pi }{\hslash }\mathop{\sum}\limits_{m\nu {{{\bf{q}}}}}| {g}_{nm\nu }({{{\bf{k}}}},{{{\bf{q}}}}){| }^{2}\\ \qquad\qquad\quad\times\, \left[({N}_{\nu {{{\bf{q}}}}}+1-{f}_{m{{{\bf{k+q}}}}})\delta ({\varepsilon }_{n{{{\bf{k}}}}}-{\varepsilon }_{n{{{\bf{k+q}}}}}-{\omega }_{\nu {{{\bf{q}}}}})\right.\\ \qquad\qquad\quad\left.+\,({N}_{\nu {{{\bf{q}}}}}+{f}_{m{{{\bf{k+q}}}}})\delta ({\varepsilon }_{n{{{\bf{k}}}}}-{\varepsilon }_{n{{{\bf{k+q}}}}}+{\omega }_{\nu {{{\bf{q}}}}})\right],\end{array}$$
(8)

where fnk and Nνq are electron Fermi-Dirac and phonon Bose-Einstein occupations in thermal equilibrium, respectively.

Cumulant method

The cumulant ansatz assumes that the retarded electron Green’s function in the time domain takes the exponential form in Eq. (2), where the cumulant function Cnk is defined as27,56

$${C}_{n{{{\bf{k}}}}}(t,T)=\int\nolimits_{-\infty }^{\infty }dE\frac{| {{{\rm{Im}}}}{{{\Sigma }}}_{n{{{\bf{k}}}}}(E+{\varepsilon }_{n{{{\bf{k}}}}},T)| }{{{{\rm{\pi }}}}{E}^{2}}({{{{\rm{e}}}}}^{-{{{\rm{i}}}}Et}+{{{\rm{i}}}}Et-1).$$
(9)

Here, εnk is the electron band energy, E the electron energy, and Σnk is the lowest-order (Fan-Migdal) e-ph self-energy57:

$$\begin{array}{lll}{{{\Sigma }}}_{n{{{\bf{k}}}}}(E,T)=\mathop{\sum}\limits_{m\nu {{{\bf{q}}}}}| {g}_{mn\nu }({{{\bf{k}}}},{{{\bf{q}}}}){| }^{2}\\ \qquad\qquad\quad\,\,\,\,\,\,\times\, \left[\frac{{N}_{\nu {{{\bf{q}}}}}+{f}_{m{{{\bf{k+q}}}}}}{E-{\varepsilon }_{m{{{\bf{k+q}}}}}+{\omega }_{\nu {{{\bf{q}}}}}+{{{\rm{i}}}}\eta }+\frac{{N}_{\nu {{{\bf{q}}}}}+1+{f}_{m{{{\bf{k+q}}}}}}{E-{\varepsilon }_{m{{{\bf{k+q}}}}}-{\omega }_{\nu {{{\bf{q}}}}}+{{{\rm{i}}}}\eta }\right]\end{array}$$
(10)

whose temperature dependence is due to the occupation factors Nνq and fnk. Due to its exponential form, our cumulant Green’s function includes contributions from higher-order e-ph interactions: it includes Feynman diagrams of all orders and sums over all the improper diagrams in which, at order n, the Fan-Migdal e-ph self-energy is repeated n times and weighted by a 1/n! factor57. After Fourier-transforming the retarded Green’s function in Eq. (2) to the energy domain, we obtain the electron spectral function using Eq. (3). We compute \({{{\rm{Im}}}}{{{\Sigma }}}_{n{{{\bf{k}}}}}(E)\) off-shell, using a fine energy E grid, and \({{{\rm{Re}}}}{{{\Sigma }}}_{n{{{\bf{k}}}}}\) on-shell at the band energy εnk, and use them as input to obtain the spectral function Ank27 as a function of electron energy E.

Dyson-Migdal spectral function

The Dyson-Migdal (DM) spectral function is given by

$${A}_{n{{{\bf{k}}}}}^{{{{\rm{DM}}}}}(E,T)=\frac{-{{{\rm{Im}}}}{{{\Sigma }}}_{n{{{\bf{k}}}}}(T)}{{[E-{\varepsilon }_{n{{{\bf{k}}}}}-{{{\rm{Re}}}}{{{\Sigma }}}_{n{{{\bf{k}}}}}(T)]}^{2}+{[{{{\rm{Im}}}}{{{\Sigma }}}_{n{{{\bf{k}}}}}(T)]}^{2}},$$
(11)

where Σnk(T) is the lowest-order e-ph self-energy [Eq. (10)] computed on-shell at the band energy εnk. The DM spectral function has a Lorentzian shape as a function of energy, with a linewidth of \(2\,{{{\rm{Im}}}}{{{\Sigma }}}_{n{{{\bf{k}}}}}(T)\) which is proportional to the e-ph scattering rate in Eq. (8), \({{{{\Gamma }}}}_{n{{{\bf{k}}}}}(T)=2\,{{{\rm{Im}}}}{{{\Sigma }}}_{n{{{\bf{k}}}}}(T)/\hslash\)54.