Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

npj Computational Materials volume 8, Article number: 141 (2022 ) Cite this article

We introduce in the many-body GW scheme the modulation of the screened Coulomb interaction W arising from the macroscopic dielectric response in the infrared. We derive expressions for the polaron binding energies, the renormalization of the effective masses and for the electron and hole relaxation times. Electron and hole mobilities are then obtained from the incorporation of appropriate scattering rules. Zinc-blende GaN and orthorhombic MAPbI3 are used as test beds finding fair agreement with results from rigorous electron-phonon coupling approaches. Although limited to polar phonons, our method has a negligible computational cost.

Nowadays band structures are routinely evaluated from first-principles using density functional theory (DFT) or more accurate methods based on it1. The most popular of the latter is the GW method which takes its name from the approximation used for writing the electron self-energy Σ. This is obtained as the product of the one particle Green’s function G with the screened Coulomb’s interaction W2.

In the GW approach the one-particle Kohn–Sham equation is replaced by a Schrödiger like one in which the exchange and correlation potential has been replaced with an energy dependent self-energy operator. As this is intrinsically complex, the quasi-particle energies (i.e., the eigen-values of such equation) gain beside a real part also an imaginary one which gives the electron or the hole inverse relaxation time.

If only electron-electron interactions are taken into account in the evaluation of Σ, the electron and the hole can relax only through Auger decay. An hole (electron) added to the system would decay to a higher (lower) energy level if the decay can provide enough energy for promoting an electron from the valence to the conduction manifold. This process can happen only if the energy gained by the hole (electron) when jumping to an higher (lower) level is larger than the quasi-particle band gap. Hence, quasi-particle energies are expected to be purely real in vicinity of the band gap.

The inclusion of electron-phonon interactions in Σ leads to a variation of both the real and imaginary part of the quasi-particle energies3. Quasi-particle energies close to the band gap acquire an imaginary part so that the relaxation time indicates the average time an electron or an hole travels before being scattered by a phonon. At the same time the real part of hole and electron energy changes determining a modification of the gap. This is due to the binding of electrons and holes with the vibrating atomic lattice leading to the formation of electron and hole polarons. A consequence is the renormalization of effective masses which register a general increase.

The inclusion of electron-phonon coupling in the calculation of band structures from first-principles is relatively new because of the high-computational cost due to slowly converging integrals over the Brillouin’s zone. Luckily, the introduction of the Wannier’s functions representation4,5 permitted to significantly reduce the computational burden6 allowing application to real materials7. The approach has then been extended to the evaluation of key transport parameters as electron and hole mobilities8. We will refer to this approach as GW-EPW from the name of the software developed by its authors9.

Very recently, the accurate GW reproduction of band structures coupled with a thorough treatment of electron-phonon coupling permitted to successfully calculate from first-principles mobilities in hybrid organic inorganic perovskites10,11,12. These materials show an astonishing potential for the fabrication of optoelectronic devices. Stable perovskite solar cells reached efficiencies as high as 25%13. Part of this success is due to the high carrier mobility which in turn is determined by the coupling with polar, infrared (IR) active, phonons. Indeed, this has been shown3 to be one of the major channels for quasi-particle relaxation in polar semiconductors8. Notwithstanding, it is worth noting that adding electron-phonon calculations upon the GW evaluation of band structures still remains an heavy computational task so that electron-phonon calculations are still relatively rare.

Recently, we showed how electron-phonon interactions modulate the exciton binding energy in hybrid organic inorganic perovskites14. Indeed, the coupling with polar phonons determines a reduction up to ~50% of the exciton binding energies. This has been possible using a simple method, we introduced in Refs. 14,15, for including from first-principles the dynamical modulation of the screened Coulomb interaction W into the many-body BSE equation for the calculation of exciton energies and amplitudes16,17,18. We assumed for W at infrared frequencies the same behavior of the macroscopic complex dielectric function which in turn is calculated using density functional perturbation theory (DFPT)19.

Here, we develop the same ideas within the GW formalism. This leads to an additional complex term in the self-energy Σ. The real part accounts for the variation of the quasi-particle energies and the imaginary one gives the relaxation time for the scattering of charges with polar phonons. We then derive electron and hole mobilities from the relaxation times assuming scattering with virtual phonons. These are defined in order to reproduce the energy loss function. The method can also be applied directly to DFT band-structures by-passing the GW step. All of this is obtained at a negligible computational cost.

We validate our scheme addressing GaN in the zinc-blende form and the prototypical MAPbI3 perovskite. In these systems scattering from longitudinal phonons is dominant12,20. Zinc-blende GaN offers a simple crystalline structure while MAPbI3 is important for its technological relevance.

We start from from the G0W0 Schrödinger’s like equation1 for the quasi-particle energy E and amplitude ϕ :

where with an hat we indicate operators, \({\hat{H}}^{KS}\) is the Kohn-Sham (KS) Hamiltonian, \({\hat{V}}_{xc}\) is the exchange and correlation potential, \({\hat{V}}_{x}\) is the (exact) exchange one and \(\hat{{{\Sigma }}}\) is the correlation part of the self-energy operator. The self-energy operator is expressed as the convolution of the DFT Green’s function G0 with the screened Coulomb interaction W0. As the latter reads \({\hat{W}}_{0}=v{\hat{\chi }}_{T}v\) , with v the bare Coulomb interaction and \({\hat{\chi }}_{T}\) the time-ordered reducible susceptibility, we write:

where \({{{\bf{x}}}}=\left({{{\bf{r}}}},\sigma \right)\) is a coupled spatial and spin coordinate, η is a positive infinitesimal which selects the contour for integration in the Riemann’s plane. As the susceptibility \({\hat{\chi }}_{{{{\rm{T}}}}}\) accounts, in principles, also for electron-phonon interactions, we can isolate the contribution to the self-energy due to vibrations:

where \({\chi }_{{{{\rm{T,el}}}}}\left(\omega \right)\) is the time-ordered complex dielectric reducible susceptibility accounting only for electronic interactions.

The Green’s function G0 reads:

where the index v runs over the Nv KS valence states ψv of energy Ev and the index c runs over the Nc conduction states ψc of energy Ec, η is a positive infinitesimal so that the formula should be read as \(\lim \eta \to {0}^{+}\) .

We isolate now the contribution to \({\hat{{{\Sigma }}}}^{{{{\rm{ep}}}}}\) from a given ψv in the G0:

Proceeding in an analogous way as done in Refs. 14,15, we can replace the time-ordered operators with the retarded (or causal) ones which correspond to observable quantities. The retarded inverse dielectric tensor \({\hat{\epsilon }}_{{{{\rm{R}}}}}\) is written in terms of the retarded susceptibility \({\hat{\chi }}_{{{{\rm{R}}}}}\) as:

The relations between retarded and time-ordered responses lead us to write:

where we indicate with \({\hat{\epsilon }}_{{{{\rm{R,el}}}}}\) the inverse retarded dielectric tensor calculated including only electron-electron interactions:

The detailed derivation of these formulae is reported in the Supplementary Information (SI). Analogously, the contribution to \({\hat{{{\Sigma }}}}^{{{{\rm{ep}}}}}\) from a given ψc in the G0 is:

Moving to the retarded responses yields:

Until this point the formulae we derived are exact. For an easier evaluation, we now approximate the dielectric tensor in Eqs. (8) and (10) with the complex macroscopic dielectric function \({\epsilon }_{R}^{{{{\rm{M}}}}}(\omega)\) times the identity operator. To this end we write:

This is the main approximation entering our method. In this way, the electron-phonon coupling is introduced phenomenologically. As we are interesting at isolating the sole contribution to \(\hat{{{\Sigma }}}\) due to phonons we calculate \({\epsilon }_{{{{\rm{R}}}}}^{{{{\rm{M}}}}}(\omega)\) in the infrared frequencies region through DFPT theory19 and then we add a single oscillation at the electronic band gap energy ωgap in order to account for the electronic screening. The corresponding oscillator strength Fel is chosen in order to furnish the high-frequency dielectric constant ϵ∞ which in turn is calculated from DFPT21:

where the index n runs over the Nph phonon modes at center zone (Γ) of frequency \({\tilde{\omega }}_{n}\) The oscillator strength Fn reads22:

where the index I runs over the Nat atoms inside the primitive cell of volume V, j runs over the 3 Cartesian directions of the atomic displacements, i runs over the 3 Cartesian directions of the electric field polarization, MI is the atomic mass of the I-th atom, Z* is the Born effective charge tensor and un is the normalized phonon eigenmode of frequency ωn. The phononic ηph and electronic ηel inverse lifetimes enter our method only as broadening factors so that the final results do not depend on them. The use of a single peak approximation for the entire electronic response is employed with the sole purpose of isolating the phononic one. Indeed, we verified that using very different values for ωgap as those from DFT and GW determines only negligible (within ca. 2%) changes in the calculated properties. It is worth stressing that with our method we consider only the coupling with polar IR-active phonons.

We apply now the method developed in the previous section to the prototypical case of a direct band polar semiconductor with band extreme at the Γ-point. The KS states take the form of Bloch’s orbitals:

where the functions uk,i have the periodicity of the direct lattice. We focus on the variation of the valence band maximum (VBM), relative to the band v due to interaction with polar phonons. First we consider the case when we start directly from a DFT calculation. The renormalised VBM energy \({E}_{{{\Gamma }}v}^{ep}\) is found solving the self-consistent one-variable equation:

where \({E}_{{{\Gamma }}v}^{DFT}\) is the energy at the VBM from DFT. This equation can be solved numerically. However, it is often observed a linear behavior of the self-energy with respect to energy23,24. This bring us to write:

In the case we start from a GW calculations also the dependence of the electronic part of the self-energy operator \({\hat{{{\Sigma }}}}^{{{{\rm{el}}}}}=\hat{{{\Sigma }}}-{\hat{{{\Sigma }}}}^{{{{\rm{ep}}}}}\) should be considered. Assuming for the latter a linear behavior in energy yields to:

We now illustrate in details the evaluation of the expectation values of \({\hat{{{\Sigma }}}}^{{{{\rm{ep}}}}}\) . As a first simplification, we address contributions only from bands degenerate or close in energy with the band v:

where with the sum \(\tilde{{{\Sigma }}}\) we denote a sum limited to the bands which are degenerate or close in energy with v. For indicating products in real space of wavefunctions, we introduce the notation:

we now move to the representation in reciprocal (or G) space so that:

where BZ stays for Brillouin’s zone and we have introduced the functions:

In Eq. (22) the wavefunctions \({u}_{{{{\bf{k}}}},v^{\prime} }\) change smoothly on k while the Coulombian potential \(\frac{4\pi {e}^{2}}{{\left|{{{\bf{k}}}}+{{{\bf{G}}}}\right|}^{2}}\) and the factor \({\tilde{A}}_{v^{\prime} }\) change sharply on it. We note also that the factor \({\tilde{A}}_{{{{\bf{k}}}}v^{\prime} }\left(\omega \right)\) depends on k solely through the \({E}_{{{{\bf{k}}}}v^{\prime} }\) energy:

This permits us to simplify:

where the sum is over a coarse mesh of k points and the symbol ∫k indicates that the integral is restricted over a fraction of the reciprocal cell around k so that summing up the entire BZ is spanned. It is worth noting that inside the square brackets the wave functions depends on k and not on \({{{\bf{k}}}}^{\prime}\) .

For integrating, the accurate behavior of \({E}_{{{{\bf{k}}}}v}^{{{{\rm{DFT}}}}}\) with respect to \({{{\bf{k}}}}^{\prime}\) is required. Wannier’s4 or Shirley’s25 interpolation could be used. For simplicity, in this work, we assume parabolic bands defined by the GW or DFT effective masses. The integration is then performed numerically. At this point, we can solve the self-consistency problem of Eq. (17) and find \({E}_{{{\Gamma }}v}^{{{{\rm{ep}}}}}\) .

For the conduction band minimum (CBM), we can easily find in a analogous way the following approximation for the expectation value of the electron-phonon part of the self-energy:

where with the sum \(\tilde{{{\Sigma }}}\) we denote a sum limited to the bands which are degenerate with c and the function \({\tilde{B}}_{{{{\bf{k}}}}c}\) reads:

We have then verified in our calculations that the expectation values \(\left\langle {\psi }_{{{\Gamma }},v}\right|{{\Sigma }}\left(\omega \right)\left|{\psi }_{{{\Gamma }},v}\right\rangle\) and \(\left\langle {\psi }_{\Gamma ,{{{\rm{c}}}}}\right|{{\Sigma }}\left(\omega \right)\left|{\psi }_{\Gamma ,{{{\rm{c}}}}}\right\rangle\) show a linear behavior close to the (renormalised) band edge energies \({E}_{{{{\boldsymbol{\Gamma }}}}v}^{ep}\) and \({E}_{{{{\boldsymbol{\Gamma }}}}c}^{ep}\) as shown in Fig. 1. This allows for an easy evaluation of the renormalization of hole and electron effective masses due to the coupling with polar phonons. To this goal we introduce the further approximation of neglecting, close the Γ-point, the k-point dependence of the self-energy. For the hole mass, we write in the case we start directly from a DFT calculation:

with \({m}_{{{{\rm{h}}}}}^{{{{\rm{ep}}}}}\) and mh and the renormalised and the DFT effective mass, respectively. Instead, when starting from a GW calculation, the derivative of the expectation value of the entire self-energy operator should be used.

Real part of the expectation values of the phononic part of the self-energy Σep(ω) calculated at the Γ-point for the (lowest) electron band (blue), heavy hole bands (red) and light holes (green) for zinc-blende GaN. The zero of the energy axis corresponds to the VBM for the holes and to the CBM for the electrons. Only for the electrons the values on both axis have been multiplied by the factor (−1). The dashed lines are a guides for the eye. The solutions of Eq. (17) are found at the crosses with the black line Σep(ω) = ω.

Analogously, the renormalised effective electron mass \({m}_{{{{\rm{e}}}}}^{{{{\rm{ep}}}}}\) is written as:

with me the DFT one.

We illustrate the method outlined in the previous sections considering zinc-blende GaN. We used the Quantum-Espresso DFT package which is based on the plane-waves pseudopotentials paradigm26,27. We coded an additional post-processing routine in order to perform the integrals of Eqs. (24) and (25) and we made it freely available (The code is available at https://gitlab.com/paoloumari/easy-eph). We choose the PBE exchange28 and correlation functional and the full-relativistic treatment of spin-orbit coupling (SOC). We adopt norm-conserving pseudo-potentials and choose a cutoff of 120 Ry for defining the plane-waves basis sets used for representing the KS wave-functions and an equally spaced mesh of 123 k-points for sampling the Brillouin’s zone. This leads to a calculated theoretical lattice constant ath = 8.60 Bohr in good agreement with the experimental (\({a}_{\exp }=8.54\) Bohr) value29. For calculating the vibrational part of the complex macroscopic dielectric function \({\epsilon }_{R}^{M}\) of Eq. (14) we use DFPT for obtaining the effective charges and the vibrational frequencies and displacements. The calculated high-frequency dielectric constant is ϵ∞ = 6.2 in agreement with the experimental value \({\epsilon }_{{{{\rm{\infty ,exp}}}}}=5.3\) 29. The static dielectric constant ϵ0,th = 10.6 is also in reasonable agreement with the experimental \({\epsilon }_{{{{\rm{0,exp}}}}}=9.7\) value29.

In Supplementary Fig. (1) we see that the function \({\tilde{A}}_{{{{\bf{k}}}}v}\) goes rapidly to zero for large Ev − ω, this corroborates the use of the approximation of Eq. (24).

The PBE-SOC band structure of zb-GaN features a direct gap at the Γ-point. Here, the lowest conduction band is doubly degenerate while the uppermost valence one is fourfold degenerate. As the degeneracy is removed moving away from the Γ-point the four upper valence bands split into two twofold degenerate heavy hole branches. A gap ΔSO = 13 meV due to spin-orbit coupling separates them from the twofold degenerate light hole branch. This figure is in good agreement with the value (19 meV) calculated in Ref. 30. In Table 1 we report the corresponding calculated effective masses along the Γ − X (100), Γ − K (110), and Γ − L (111) together with the average values. These have been obtained fitting band dispersions up to kmax = 0.1 Bohr−1. In the same table we also show OEPx(cLDA)+G0W0 results from Ref. 31. As these calculations neglect SOC, the branches we labeled hh1 and hh2 become degenerate. Indeed, the largest difference is found for the hh1 effective masses (2.374 vs. 1.335) while much smaller differences are found for the other valence bands. The conduction band effective masses are almost equal.

We illustrate how our approach works in practice considering the electron and the heavy and light hole branches. The expectation values of the phononic contribution to the self-energy at the Gamma-point are reported in Fig. 1 as a function of the energy difference with respect to the valance band maximum for the holes and the conduction band minimum for the electrons, respectively. For the case of electrons we multiply both energies displayed on the horizontal and vertical axis by the factor (−1) (i.e., both are negative). For the two heavy holes branches, which are degenerate at the Γ-point we average, for simplicity, the values of the self-energy. This procedure is justified by a found maximum difference in the calculated energies of ca. 5 meV. The shifts of the energy bands are then found solving Eq. (17). This is done graphically in Fig. 1 and the results are reported in Table 2.

The results are in agreement with those from quasi-particle self-consistent GW or DFT plus estimates of electron-phonon coupling of Refs. 32 and33, respectively, for a temperature T = 0 K. The agreement is also good with the accurate results obtained in Ref. 34 from a rigorous many-body perturbation theory treatment of the electron-phonon scattering, again at zero temperature, in this work only the change in the gap has been reported.

We now use Eqs. (27) and (28) for calculating the renormalization of the effective masses due to the interaction with longitudinal phonons. As reported in Table 1 we find an increase of ~20% for the hole branches and a smaller one of ~10% for the electrons.

The imaginary part of the self-energy is the key ingredient for the evaluation of electron and hole mobilities. Indeed, it is proportional to the inverse relaxation time. As the underlying scattering channels involve either the creation or the annihilation of phonons, we will add in our scheme simplified scattering rules which will depend on the Bose-Einstein statistics of phonons at a given temperature. We focus on the electrons and holes close to the Γ-point. We begin writing Eq. (24) for the valence band at the k-point k (close to Γ):

Now we can plug \({A}_{{{{\bf{k}}}}^{\prime} ,v}\left(\omega \right)\) from Eq. (8) and employ the diagonal approximation for the dielectric operator. We note that only the term in \({\epsilon }_{R}^{-1}\) contributes to the imaginary part of \({A}_{{{{\bf{k}}}}^{\prime} ,v}\left(\omega \right)\) :

For simplicity, we refer here to band extremes (\({E}_{{{\Gamma }}v}^{{{{\rm{DFT}}}}}\) and \({E}_{{{\Gamma }}c}^{{{{\rm{DFT}}}}}\) ) and effective masses (mh and me) taken from a starting DFT calculation. In case we would start from a GW one, they should be replaced by the corresponding GW quantities.

We suppose that the energy loss function \({{\epsilon }_{{{{\rm{R}}}}}^{{{{\rm{M}}}}}}^{-1}(\omega)\) can be written in terms of single peaks:

The peak positions ωj and intensities Lj are found by fitting the \({{{\rm{Im}}}}{\epsilon }_{R}^{{{{\rm{M}}}}}{(\omega)}^{-1}\) function calculated through DFPT. We then verified that small variations in the quality of the fit, even affecting the number of peaks, do not yield any significant change in the final results.

Now, according to Eq. (8) we have to evaluate the following integral for:

where the sign − 1 comes from the sign of the argument of \({\hat{\epsilon }}_{{{{\rm{R}}}}}^{-1}\) in Eq. (8) and the Dirac delta function assures that we have \(\left|{{{\bf{k}}}}\right| \,<\, \left|{{{\bf{k}}}}^{\prime} \right|\) , so we re-write the integral as:

where the factors M are defined as:

For simplicity we neglect self-consistency. Hence, for a state ψk,v we just calculate the imaginary part of the expectation value of the self-energy at the energy \({E}_{{{{\bf{k}}}}v}^{{{{\rm{DFT}}}}}={E}_{{{\Gamma }}v}^{{{{\rm{DFT}}}}}-\frac{{k}^{2}}{2{m}_{v}}\) . We set

With an analogous reasoning we find for the case of conduction bands:

The quantum nature of the lattice vibrations implies that an electron (or hole) of energy ϵ, with respect to the CBM (VBM), can be scattered through the creation of a phonon of energy ϵph only if ϵ > ϵph. In contrast a phonon of arbitrary energy can always be annihilated. The probability of the latter event is proportional to the corresponding occupancy factor given by the Bose-Einstein distribution \({\left({{{{\rm{e}}}}}^{\frac{{\epsilon }_{{{{\rm{ph}}}}}}{{k}_{{{{\rm{B}}}}}T}}-1\right)}^{-1}\) .

This is in contrast with Eqs. (39) and (40) which on one side do not depend directly on phonons energy but on the frequencies (labeled as ωj) of the peaks in the energy-loss function (\(-{{{\rm{Im}}}}{\epsilon }_{{{{\rm{R}}}}}^{-1}\) ) and on the other do not involve phonon statistics. To tackle this issue we first associate to each peak ωj of \({{{\rm{Im}}}}({\epsilon }_{{{{\rm{R}}}}}^{-1})\) a virtual phonon energy \({\omega }_{j}^{{{{\rm{ph}}}}}\) assuming the same behavior of an isolated vibration35:

It is worth noting that, in case of well energy-separated peaks in the energy loss function, virtual phonons frequencies match the physical ones. Then, we include the scattering rules at temperature T through:

where σ is the Heaviside step function. Analogously, for conduction electrons we write:

For evaluating the electron and hole mobility at a given temperature T we rely on the equilibrium Boltzmann transport equation35. We address only the the topmost valence and lowermost conduction bands and these are considered to be parabolic. This yields the following expressions for the effective valence \({{{{\mathcal{N}}}}}_{v}\left(T\right)\) and conduction \({{{{\mathcal{N}}}}}_{c}\left(T\right)\) densities of states:

Now, we can find the chemical potential \(\mu \left(T\right)\) :

where we suppose to have an intrinsic (not-doped) material. The chemical potential enters the definition of the intrinsic density of charge carriers:

for the valence band and

for the conduction one. For achieving the mobilities we adopt the Boltzmann transport equation in the relaxation time approximation. It is worth noting that it can be quite approximate with respect to iterative solutions36.

The diagonal terms of the hole mobility μv tensor are then easily achieved:

where α runs on the three Cartesian directions, \({{{{\bf{v}}}}}_{v}\left({{{\bf{k}}}}\right)\) is the velocity of the hole of band v and momentum k, \({\tau }_{v}\left({{{\bf{k}}}}\right)\) is its relaxation time which is given by the relation:

and fv,k is the hole Fermi–Dirac distribution:

within the parabolic band approximation, the velocity vector reads:

The diagonal terms electron mobility μc tensor are given by:

where the electron Fermi-Dirac distribution reads:

Finally, we average over the three Cartesian directions α.

For calculating charge carrier relaxation times and mobilities along Eqs. ((43), (43), (50), (54)) we need the high-frequency dielectric constant, the vibrational frequencies and infrared intensities at the sole Γ point, which can be easily evaluated through DFPT, band gap and effective masses, which can be evaluated either with DFT or with the GW method. Interestingly, using the DFT band gap instead of the GW one induces only negligible changes in the expectation values of the phononic part of the self-energy and in the relaxation times. Even the mobilities are within 2% those obtained with the proper GW band-gap. In contrast,the dependence on effective masses and high-frequency dielectric constant is much stronger.

For performing these calculations, we have developed a python code which takes as input the output files of the pw.x and ph.x codes of the Quantum-Espresso distribution. As a first step, our code computes the complex dielectric function and the energy loss function \(-{{{\rm{Im}}}}({{\epsilon }_{{{{\rm{R}}}}}^{{{{\rm{M}}}}}}^{-1}(\omega))\) . The latter is then fitted with Gaussian functions. Their centers correspond to the positions of the peaks of \(-{{{\rm{Im}}}}({{\epsilon }_{{{{\rm{R}}}}}^{{{{\rm{M}}}}}}^{-1}(\omega))\) higher than 1/100 of the highest peak. This allows for a straightforward and automatic calculation of the Lj and ωj parameters. We observed only negligible changes setting the threshold to 1/1000. Then, the factors M of Eq. (37) are found by interpolating from a grid of values previously evaluated numerically (We used the Wolfram Mathematica Package). We have made this code publicly available (The python code is available at: https://gitlab.com/paoloumari/easy-mobility).

For a significant test of our approach we opted for the perovskite MAPbI337 well known for its employment in solar cells38,39,40. We focused on the low temperature orthorhombic phase41. For this structure Giustino and co-workers performed full state-of-the-art electron-phonon calculations at the GW level of relaxation times11 and mobilities12,42 which we can use as a benchmark. Moreover, we make use, for consistency and better comparison, of the DFT-SOC effective masses (mh = 0.138, me = 0.118) in units of the electron mass when comparing with Ref. 11 and of the GW-SOC effective masses (mh = 0.182, me = 0.211) in units of the electron mass when comparing with Ref. 42. These values together with the band gap (1.38) and high frequency dielectric constant (5.9) are taken from Ref. 12. For a better comparison we also took IR frequencies and activities43 reported by the same group.

In Fig. 2, we display our calculated imaginary part of the quasi-particle energy for the hole branch as function of the hole energy (defined with respect to the VBM) at the temperature T = 1 K together with accurate GW-EPW results from Ref. 11. We note a satisfactory agreement with an overall discrepancy of ca. 25% despite an almost negligible computational cost (the calculation took ca. 2 mins on a single laptop core).

Imaginary part of the expectation value of the phononic part of the self-energy as a function of the hole energy (w.r.t to the energy of the VBM) for orthorhombic MAPbI3 at the temperature T = 1 K. Solid line: results with this method. Gray squares: results from Ref. 11.

The main difference between our approach and the standard Fan-Migdal self-energy coupled with a model Fröhlich vertex and a parabolic bands approximation11,44 stems from the diagrammatic treatment of the electron-phonon coupling of the latter to be contrasted with the mean field coupling of the first. Interestingly, some anharmonic effects could, in principles, be accounted by our method directly from the modulation of the macroscopic complex dielectric function, which could, in turn, be calculated through molecular dynamics45.

In Table 3 we report electron and hole mobilities at ambient temperature (T = 300 K). We see a reasonable correspondence with the GW-EPW figures (average error of 17.0 cm2 V−1 s−1) which must be taken as a point of reference. As a consequence we find lower values than those from Refs. 46,47 which neglect coupling with multiple phonons. It is worth noting that the experimental value reported in the same table and measured for a single crystal is somehow lower than the theoretical counterparts as additional relaxation channels are in action. These include scattering with defects and grain boundaries42. We verified that changing ηph and ηel of Eq. (14) from 0.1 meV to 0.01 meV results in a change of the mobilities within 0.3%.

We also calculated mobilities in MAPbI3 as a function of temperature (See Supplementary Figure (1)) We found a behavior ∝ T−1.1.

For the sake of completeness we calculated the zb-GaN mobility and the renormalization of MAPbI3 bands. For zb-GaN at T = 300 K we found μe = 1474 cm2 V−1 s−1 and μh = 75.7 cm2 V−1 s−1 using PBE-SOC effective masses. Using PBE-SOC-EPH ones we found μe = 1300 cm2 V−1 s−1 and μh = 55.4 cm2 V−1 s−1. Experimental values are available for epitaxially grown zb-GaN. In Ref. 48 an hole mobility of 39 cm2 V−1 s−1 is reported while in Ref. 49 an hole mobility of 215 cm2 V−1 s−1 is given both measured at room temperature. In free standing zb-Gan, at room temperature, an hole mobility of 47 ± 3 cm2 V−1 s−1 has been measured50.

We calculated band renormalization and effective mass renormalization for orthorombic MAPbI3 using the DFT band-structure in analogy with Ref. 11. In our method both quantities do not depend on temperature and should be compared with those at T = 0. We found mass enhancement factors of λh = 16% for the hole and of λe = 13% for the electrons in agreement with the value λ = 18% for both charge carriers from the same Reference. We also found that the VBM moves up of 8 meV and the CBM moves down of 7 meV. Interestingly, as our method accounts only for polar interactions, we register only reductions in the band gaps.

With our work we have introduced and validated a first-principles scheme for the rapid calculation of band renormalizations, relaxation times and mobilities. This involves the following assumptions and approximations: (1) the coupling with the phonons is limited to the modulation of the screened Coulomb interaction W which is approximated considering a diagonal dielectric matrix depending only on the macroscopic complex dielectric function (Eq. (13)). (2) In the Green’s function G entering the self-energy only the upper valence and lowest conduction bands are considered and their energies are evaluated with a quadratic approximation. Then, for evaluating the imaginary part of the self-energy and charge mobility, we further assume: (3) the energy loss functions can be written in terms of single peaks (Eq. (34)). (4) We assign to each peaks a virtual phonon.

In contrast with the standard Fan-Migdal self-energy3, which is built using the many-body perturbation formalism for both electrons and phonons, we have introduced a simple phenomenological treatment of the electron-phonon coupling through the modulation of W. This brings the advantage of having to calculate, or even to take from experiment, the sole complex macroscopic dielectric function. In principles, this could be calculated accounting for non-harmonic terms.

In a good number cases, being limited to longitudinal phonons is not an issue as they give the most relevant contributions.

However, it has been shown20,44,51,52 that in some crystals longitudinal polar modes play no role or only a marginal one. Some prototypical cases are Diamond, Silicon, c-BN and 3C-SiC20. Indeed, our method should be applied for systems where scattering from longitudinal optical phonons is dominant. From a look at Eqs. ((23), (26), (43), (44)) it can be seen that stronger scattering is expected for lower frequencies of the LO phonon modes at center zone and higher IR oscillator strengths. In Diamond and Silicon the IR oscillator strengths are vanishing because of symmetry while the lowest LO modes of c-BN and 3C-SiC are found at relatively high frequency (1280 cm−1 53 and 955 cm−1 54). We think that a good criterion for estimating the range of validity of our method is to look at the frequency of the lowest IR active LO mode. Largest such a frequency, less reliable results are expected. Nevertheless, for a better estimation the strength of the oscillator strengths should also be considered. Unfortunately, such an analysis cannot determine the relevance of other scattering channels.

The calculation can be initiated either directly from DFT or following a GW step. We have made the codes for evaluating bands renormalization and mobilities publicly available (The code is available at https://gitlab.com/paoloumari/easy-eph; The python code is available at: https://gitlab.com/paoloumari/easy-mobility) together with the data for reproducing the calculations. In contrast with other more accurate approaches, ours can be advantageous because of its very limited computational cost. Indeed, it could be easily used within automatized tools for materials discovery55. In a future work, we envisage to improve our scheme in order to account for anisotropies in the band structure and in the dielectric response.

DFT calculations were based on the plane-waves pseudo-potential paradigm. Wave-functions are represented using plane-waves basis set and the interaction between valence and core electrons is modeled with norm-conserving pseudo-potentials. We used local and semi-local approximations for the DFT exchange and correlation potential.

Input files for reproducing the results regarding band renormalizations are publicly available (The code is available at https://gitlab.com/paoloumari/easy-eph). Input files for reproducing the results regarding mobilities are publicly available (The python code is available at: https://gitlab.com/paoloumari/easy-mobility).

DFT calculations have been carried with the Quantum-Espresso package26,27 which is freely available. The code for the evaluation of bands and effective masses renormalization is freely available (The code is available at https://gitlab.com/paoloumari/easy-eph). The code for the calculation of charge mobility is freely available (The python code is available at: https://gitlab.com/paoloumari/easy-mobility).

Martin, R. M., D. M., C., Lucia, R. Interacting Electrons: Theory and Computational Approaches (Cambridge University Press, 2016).

Hedin, L. New method for calculating the one-particle green’s function with application to the electron-gas problem. Phys. Rev. 139, A796–A823 (1965).

Giustino, F. Electron-phonon interactions from first principles. Rev. Mod. Phys. 89, 015003 (2017).

Marzari, N. & Vanderbilt, D. Maximally localized generalized wannier functions for composite energy bands. Phys. Rev. B 56, 12847–12865 (1997).

Souza, I., Marzari, N. & Vanderbilt, D. Maximally localized wannier functions for entangled energy bands. Phys. Rev. B 65, 035109 (2001).

Giustino, F., Cohen, M. L. & Louie, S. G. Electron-phonon interaction using wannier functions. Phys. Rev. B 76, 165108 (2007).

Giustino, F., Cohen, M. L. & Louie, S. G. Small phonon contribution to the photoemission kink in the copper oxide superconductors. Nature 452, 975–978 (2008).

Poncé, S., Margine, E. R. & Giustino, F. Towards predictive many-body calculations of phonon-limited carrier mobilities in semiconductors. Phys. Rev. B 97, 121201 (2018).

Poncé, S., Margine, E., Verdi, C. & Giustino, F. Epw: electron-phonon coupling, transport and superconducting properties using maximally localized wannier functions. Comput. Phys. Commun. 209, 116–133 (2016).

Wright, A. D. et al. Electron–phonon coupling in hybrid lead halide perovskites. Nat. Commun. 7, 11755 (2016).

Schlipf, M., Poncé, S. & Giustino, F. Carrier lifetimes and polaronic mass enhancement in the hybrid halide perovskite ch3nh3pbi3 from multiphonon Fröhlich coupling. Phys. Rev. Lett. 121, 086402 (2018).

Poncé, S., Schlipf, M. & Giustino, F. Origin of low carrier mobilities in halide perovskites. ACS Energy Lett. 4, 456–463 (2019).

Jeong, M. et al. Stable perovskite solar cells with efficiency exceeding 24.8% and 0.3-v voltage loss. Science 369, 1615–1620 (2020).

Umari, P., Mosconi, E. & De Angelis, F. Infrared dielectric screening determines the low exciton binding energy of metal-halide perovskites. J. Phys. Chem. Lett. 9, 620–627 (2018).

Adamska, L. & Umari, P. Bethe-salpeter equation approach with electron-phonon coupling for exciton binding energies. Phys. Rev. B 103, 075201 (2021).

Rohlfing, M. & Louie, S. G. Electron-hole excitations in semiconductors and insulators. Phys. Rev. Lett. 81, 2312–2315 (1998).

Albrecht, S., Reining, L., Del Sole, R. & Onida, G. Ab initio calculation of excitonic effects in the optical spectra of semiconductors. Phys. Rev. Lett. 80, 4510–4513 (1998).

Strinati, G. Application of the green’s functions method to the study of the optical properties of semiconductors. Riv. del Nuovo Cim. 11, 1–86 (1988).

Baroni, S., de Gironcoli, S., Dal Corso, A. & Giannozzi, P. Phonons and related crystal properties from density-functional perturbation theory. Rev. Mod. Phys. 73, 515–562 (2001).

Poncé, S. et al. First-principles predictions of hall and drift mobilities in semiconductors. Phys. Rev. Res. 3, 043022 (2021).

Baroni, S., Giannozzi, P. & Testa, A. Green’s-function approach to linear response in solids. Phys. Rev. Lett. 58, 1861–1864 (1987).

Umari, P. & Pasquarello, A. Infrared and raman spectra of disordered materials from first principles. Diam. Relat. Mater. 14, 1255–1261 (2005).

Hybertsen, M. S. & Louie, S. G. Electron correlation in semiconductors and insulators: band gaps and quasiparticle energies. Phys. Rev. B 34, 5390–5413 (1986).

Kotani, T., van Schilfgaarde, M. & Faleev, S. V. Quasiparticle self-consistent GW method: a basis for the independent-particle approximation. Phys. Rev. B 76, 165106 (2007).

Prandini, G., Galante, M., Marzari, N. & Umari, P. Simple code: optical properties with optimal basis functions. Comput. Phys. Commun. 240, 106–119 (2019).

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).

Giannozzi, P. et al. Advanced capabilities for materials modelling with quantum ESPRESSO. J. Phys. Condens. Matter 29, 465901 (2017).

Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 77, 3865–3868 (1996).

Levinshtein, M. E., Rumyantsev, S. L. & Shur, M. S. Properties of Advanced Semiconductor Materials GaN, AlN, InN, BN, SiC, SiGe (Wiley, 2001).

Kim, K., Lambrecht, W. R. L., Segall, B. & van Schilfgaarde, M. Effective masses and valence-band splittings in GaN and AlN. Phys. Rev. B 56, 7363–7375 (1997).

Rinke, P. et al. Consistent set of band parameters for the group-III nitrides AlN, GaN, and InN. Phys. Rev. B 77, 075202 (2008).

Lambrecht, W. R. L., Bhandari, C. & van Schilfgaarde, M. Lattice polarization effects on the screened Coulomb interaction W of the GW approximation. Phys. Rev. Mater. 1, 043802 (2017).

Nery, J. P. & Allen, P. B. Influence of Fröhlich polaron coupling on renormalized electron bands in polar semiconductors: results for zinc-blende GaN. Phys. Rev. B 94, 115135 (2016).

Kawai, H., Yamashita, K., Cannuccia, E. & Marini, A. Electron-electron and electron-phonon correlation effects on the finite-temperature electronic and optical properties of zinc-blende GaN. Phys. Rev. B 89, 085202 (2014).

Grosso, G. & Parravicini, G. P. Solid State Physics (Elsevier, 2006).

Poncé, S., Li, W., Reichardt, S. & Giustino, F. First-principles calculations of charge carrier mobility and conductivity in bulk semiconductors and two-dimensional materials. Rep. Prog. Phys. 83, 036501 (2020).

Kojima, A., Teshima, K., Shirai, Y. & Miyasaka, T. Organometal halide perovskites as visible-light sensitizers for photovoltaic cells. J. Am. Chem. Soc. 131, 6050–6051 (2009).

Im, J.-H., Lee, C.-R., Lee, J.-W., Park, S.-W. & Park, N.-G. 6.5% efficient perovskite quantum-dot-sensitized solar cell. Nanoscale 3, 4088–4093 (2011).

Liu, M., Johnston, M. B. & Snaith, H. J. Efficient planar heterojunction perovskite solar cells by vapour deposition. Nature 501, 395–398 (2013).

Burschka, J. et al. Sequential deposition as a route to high-performance perovskite-sensitized solar cells. Nature 499, 316–319 (2013).

Baikie, T. et al. Synthesis and crystal chemistry of the hybrid perovskite (CH3NH3)PbI3 for solid-state sensitised solar cell applications. J. Mater. Chem. A 1, 5628–5641 (2013).

Xia, C. Q. et al. Limits to electrical mobility in lead-halide perovskite semiconductors. J. Phys. Chem. Lett. 12, 3607–3617 (2021).

Pérez-Osorio, M. A. et al. Vibrational properties of the organic–inorganic halide perovskite CH3NH3PbI3 from theory and experiment: factor group analysis, first-principles calculations, and low-temperature infrared spectra. J. Phys. Chem. C 119, 25703–25718 (2015).

Miglio, A. et al. Predominance of non-adiabatic effects in zero-point renormalization of the electronic band gap. npj Comput. Mater. 6, 167 (2020).

Dubois, V., Umari, P. & Pasquarello, A. Dielectric susceptibility of dipolar molecular liquids by ab initio molecular dynamics: application to liquid HCl. Chem. Phys. Lett. 390, 193–198 (2004).

Frost, J. M. Calculating polaron mobility in halide perovskites. Phys. Rev. B 96, 195202 (2017).

Lacroix, A., de Laissardière, G. T., Quémerais, P., Julien, J.-P. & Mayou, D. Modeling of electronic mobilities in halide perovskites: adiabatic quantum localization scenario. Phys. Rev. Lett. 124, 196601 (2020).

Lin, M. E., Xue, G., Zhou, G. L., Greene, J. E. & Morkoç, H. p-type zinc-blende GaN on GaAs substrates. Appl.Phys. Lett. 63, 932–933 (1993).

As, D. J. Electrical and optical properties of Mg doped MBE grown cubic GaN epilayers. Phys. Status Solidi B 210, 445–449 (1998).

Novikov, S. V. et al. Growth and characterization of free-standing zinc-blende (cubic) GaN layers and substrates. Semicond. Sci. Technol. 23, 015018 (2007).

Jhalani, V. A., Zhou, J.-J., Park, J., Dreyer, C. E. & Bernardi, M. Piezoelectric electron-phonon interaction from ab initio dynamical quadrupoles: impact on charge transport in wurtzite GaN. Phys. Rev. Lett. 125, 136602 (2020).

Brunin, G. et al. Electron-phonon beyond fröhlich: dynamical quadrupoles in polar and covalent solids. Phys. Rev. Lett. 125, 136601 (2020).

Ohba, N., Miwa, K., Nagasako, N. & Fukumoto, A. First-principles study on structural, dielectric, and dynamical properties for three BN polytypes. Phys. Rev. B 63, 115207 (2001).

Tong, Z., Liu, L., Li, L. & Bao, H. Temperature-dependent infrared optical properties of 3C-, 4H- and 6H-SiC. Physica B Condens. Matter 537, 194–201 (2018).

Pizzi, G., Cepellotti, A., Sabatini, R., Marzari, N. & Kozinsky, B. AiiDA: automated interactive infrastructure and database for computational science. Comput. Mater. Sci. 111, 218–230 (2016).

We acknowledge support from Fondazione Cariparo through grant Ricerca Scientifica di Eccellenza 2018 project SYNERGY.

Dipartimento di Fisica e Astronomia, Università di Padova, via Marzolo 8, I-35131, Padova, Italy

CNR-IOM Istituto Officina Dei Materiali, Consiglio Nazionale Delle Ricerche, Area Science Park, Strada Statale 14, km 163.5, I-34149, Trieste, Italy

You can also search for this author in PubMed Google Scholar

P.U. conceived the work, implemented the codes, carried the calculations and wrote the paper.

The author declares no competing interests.

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

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/.

Umari, P. Inclusion of infrared dielectric screening in the GW method from polaron energies to charge mobilities. npj Comput Mater 8, 141 (2022). https://doi.org/10.1038/s41524-022-00812-9

DOI: https://doi.org/10.1038/s41524-022-00812-9

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative

npj Computational Materials (npj Comput Mater) ISSN 2057-3960 (online)

Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.