Microscopic calculation of the optical properties and intrinsic losses in the methylammonium lead iodide perovskite system

Perovskite crystals containing organic cations have gained a lot of attention in recent years, especially for solar cell applications. While high efficiency has been reached, the relatively short life-time of perovskite based solar cells currently still constitutes a significant challenge.^1–3 The perovskites of interest here can be summarized by the chemical formula AMX₃ where A needs to be a large compound such as Cs or even better, to increase the stability, an organic cation. For M, typical atoms are Sn or Pb, and X can be one of the halides I, Cl, or Br as well as a mixture of these. Usually, the organic cation is methylammonium (MA) which is $CH3NH3+$ or formamidinium (FA) which is CH(NH₂)$2+$⁠, but there are also other candidates such as azetidinium or aziridinium.⁴ Most of the perovskite combinations that can be realized with the listed compounds are semiconductors with band gaps ranging from approximately 1.0 eV to 3.8 eV⁵ where the choice of the halide atom has the strongest impact on the actual band gap value. Depending on temperature and composition, perovskites can be in a cubic, tetragonal, orthorhombic, rhombohedral, or monoclinic structure. For changing system temperatures, interesting phase transitions have been observed.^6–9 

In the current study, we investigate the optical properties of the tetragonal phase of methylammonium lead iodide (MAPbI₃). This material system has a band gap of approximately 1.64 eV at room temperature¹⁰ and is one of the most intensively analyzed perovskites for solar cell applications. The band structure of MAPbI₃ was previously studied using density functional theory (DFT)^11–13 or the GW approach.^5,14–16 It was shown that DFT calculations which neglect the spin-orbit interaction usually yield band gaps close to the experimentally measured values. However, taking spin-orbit interaction into account severely reduces the band gap^11–15,17 and slightly shifts the conduction band minimum away from the Γ point of the Brillouin zone such that the optical transitions become somewhat indirect.¹⁵ More advanced DFT calculations utilizing the self-consistent PBE0¹¹ or LDA-1/2¹³ approach were shown to perform equally well as GW calculations when spin-orbit interaction is included.

There have been several studies investigating the optical properties of MAPbI₃ in order to understand the details behind the good energy conversion in this material. A significant range of values has been reported for the static and the high-frequency dielectric constants.^{5,12,14,15,18–21} Since the computed value for the exciton binding energy E_b depends on the choice of the dielectric constant, also different values for E_b have been obtained.^{5,19,20,22–24} Furthermore, there is an ongoing discussion to which degree phonons participate in the screening of the Coulomb interaction potential of the charge carriers. In particular, a large phonon contribution would result in a severe reduction in the exciton binding energy^5,19,22 and at room temperature, free carriers would predominate which is one possible explanation for the good energy conversion.^19,22 An alternative explanation is given by Bokdam et al.⁵ who predicted the formation of a polaronic state with a smaller band gap and reduced probability of exciton formation. The range of currently discussed exciton binding energies is summarized in Ref. 25. The fact that the perovskites under investigation have several active LO phonon modes whose eigen-frequencies are currently not unambiguously identified^18,26,27 further adds to the problem to identify a reliable value of the exciton binding energy.

Under many conditions, e.g., during laser operation (high carrier densities) or for solar cell operation (lower carrier densities), the performance of the opto-electronic device is limited by intrinsic carrier loss processes such as spontaneous emission and/or Auger recombination. Several experimental studies analyzed the strength of the spontaneous emission and Auger processes in MAPbI₃, yielding spontaneous emission coefficients in the range of 6 × 10⁻¹¹ cm³/s to 2 × 10⁻⁹ cm³/s^28–32 and Auger coefficients in the range of 10⁻²⁸ cm⁶/s to 10⁻²⁷ cm⁶/s.^29–31 A comparison between the bands involved in the Auger processes of hybrid perovskites and those of III-V semiconductors can be found in Ref. 12. As pointed out in Ref. 25, systematic theoretical investigations on the Auger processes in perovskite materials are deemed highly necessary.

In this work, we use an approach that combines DFT³³ calculations for the electronic properties of MAPbI₃ with a microscopic many-body scheme to compute optical properties such as carrier density dependent absorption and photoluminescence (PL) spectra as well as the temperature and density dependent spontaneous emission and Auger loss coefficients. This approach has previously been used successfully to study a wide range of III-V semiconductor material systems.^34,35 Additionally, we calculate the direction dependent exciton binding energies and discuss these results in the context of the current debate.

The DFT calculations were performed using the Vienna Ab initio Simulation Package (VASP)^36–39 with the Projector Augmented-Wave (PAW) pseudopotential method.^40,41 In all DFT calculations, the k-space was sampled with a 4 × 4 × 3 Monkhorst-Pack mesh.⁴² For the lattice optimization, the local density approximation (LDA) as parameterized by Perdew and Zunger⁴³ was utilized, converging the Hellmann-Feynman forces on the atoms below 10⁻² eV/Å and the energies to an accuracy of 10⁻⁶ eV. An energy-cutoff for the plane waves of 800 eV was chosen. Optimization of the atom positions, cell shape, and volume was performed from a perfect tetragonal crystal structure as a starting point. This results in a pseudo-tetragonal lattice with a = 8.524 Å, b = 8.529 Å, and c = 12.623 Å which are somewhat smaller than the experimental lattice constants,⁸ and the lattice vectors are slightly non-orthogonal. We omitted the inclusion of van-der-Waals corrections because a recent study on several perovskite systems showed that these yield no systematic improvement.⁴⁴ For the calculation of the band structures and wave functions, an energy-cutoff of 400 eV was used, the energies were converged to an accuracy of 10⁻⁴ eV, and spin-orbit coupling was enabled. In addition, the LDA-1/2 method⁴⁵ was used as this method was shown to yield band gaps that are close to the ones obtained with GW for perovskites.¹³ The advantage over using GW is that the LDA-1/2 method is computationally equally demanding to using the LDA. It is based on Slater’s half occupation scheme⁴⁶ and aims at an improved inclusion of the electron self-energies. For iodine, we find a cutoff $RCUTI=3.71$ a.u. that maximizes the band gap which is similar to the value Tao, Cao, and Bobbert¹³ reported (3.76 a.u.). For lead, we obtained $RCUTPb=2.00$ a.u., while Tao, Cao, and Bobbert¹³ found $RCUTPb=2.18$ a.u. Yin et al.¹⁷ found that electronic states, the MA⁺ cations contribute to, are far from the band edges. Therefore, and because it was found that applying the LDA-1/2 technique to the cation has a negligible effect,⁴⁷ we do not calculate corrected potentials for the C, N, and H atoms. For the non-self-consistent calculations, the k-path was sampled throughout the first Brillouin zone with 577 k-points in the ΓA direction ([111] direction), 521 k-points in the ΓM direction ([110] direction), 445 k-points in the ΓR direction ([011] direction), 368 k-points in the ΓX direction ([010] direction), and 249 k-points in the ΓZ direction ([001] direction). This corresponds to approximately one k-point for 10⁻³ Å⁻¹.

From the wave functions, the dipole and Coulomb matrix elements are calculated. Together with the band structure, these are needed as input for calculating the optical properties. To determine the absorption/gain spectra, we use the semiconductor Bloch equations.^48–50 These are derived on a microscopic level by expanding the macroscopic optically induced polarization P(t) in a Bloch basis,⁴⁹ 

Here, V is the crystal volume, $dkλν$ are the dipole matrix elements, and the microscopic interband polarizations are denoted as $pkνλ=⟨aλ,k†aν,−k†⟩$⁠, where the operator $aλ,k†$ creates an electron in the conduction band λ at the wave vector k and $aν,−k†$ creates a hole in the valence band ν at the wave vector −k, respectively. The carrier occupation probabilities in the bands are given by $fkλ=⟨aλ,k†aλ,k⟩$ (electrons) and $fkν=⟨aν,−k†aν,−k⟩$ (holes). In order to derive the time dynamics, we use Heisenberg’s equation of motion,

Here, $H$ is the many-body Hamilton operator of the electron-hole system interacting with a classical light field [see Eqs. (1)–(4) of the supplementary material].

In the following, we focus the discussion on the dynamics of the interband polarization since the methodology to obtain the equations for the carrier occupations is very similar and described in Ref. 50. For the interband polarization, Eq. (2) yields Eqs. (2.1)–(2.4) of Ref. 51 which in our case can be written as

Here

are the electron (superscript λ) and hole (superscript ν) energies that are re-normalized due to the Coulomb interaction. In these equations, the energy dispersion and Coulomb matrix elements obtained from DFT calculations enter in $ε$_k and V_k−q, respectively. Furthermore, Ω_k is given by

where the first term describes the coupling between states in bands λ and ν by the laser field E(t) and the second term describes again a re-normalization caused by the Coulomb interaction V_k−q. Here, the dipole matrix elements d_k and the Coulomb matrix elements V_k−q enter. In Eq. (3), all contributions beyond the Hartree-Fock approximation are included in the scattering term that we address later. The differential equation for the interband polarization $pkνλ$ couples to differential equations for the electron and hole occupations. This system of differential equations is the semiconductor Bloch equations which are solved numerically. From the microscopic interband polarizations, the time dependent macroscopic polarization is obtained via Eq. (1). To obtain the absorption spectra α(ω), a Fourier transformation to the frequency space is applied to both P(t) and E(t) such that we can evaluate

Here, n₀ is the background refractive index and c₀ is the speed of light in vacuum.

The scattering term in Eq. (3) includes all interaction contributions beyond the Hartree-Fock level. An exact evaluation of these is fundamentally impossible since one has to deal with the infinite hierarchy of many-body correlations of ever increasing complexity. Hence, one has to resort to physically motivated approximations. As long as the system is only weakly excited, i.e., at low excitation densities, mainly the scattering of carriers with phonons is important, whereas at higher densities, the Coulombic scattering among the charge carriers becomes dominant. In the dynamic equation of the microscopic polarizations, all these contributions lead to damping terms which describe dephasing,⁵⁰ i.e., the decay of the interband coherences induced by the resonant optical excitation. While the detailed Boltzmann-like equation for the phonon-carrier scattering can in principle be included in the semiconductor Bloch equations, we do not follow this approach in the present case due to the strong phonon-carrier coupling present in the perovskite materials⁵ and the significant uncertainties in the relevant interaction matrix elements caused by the large spread of reported values for the dielectric constants and the LO phonon frequencies discussed above. Therefore, in our current analysis, the dephasing caused by carrier-phonon scattering is approximated by a phenomenological dephasing constant γ = 165 meV (equivalent to a dephasing of 25 fs) such that the corresponding scattering term takes the form

Additionally, we include carrier-carrier scattering in the second term of Eq. (8) at the level of the 2nd Born approximation where we keep the contributions up to the quadratic order in the Coulomb potential.⁴⁹ Formally, we evaluate the quantum Boltzmann type scattering integral derived in Ref. 49 [Eqs. (10) and (11) within]. Here, we take Coulomb-screening by the excited charge carriers into account by evaluating the longitudinal dielectric function at the level of the Lindhard formula.⁵⁰ 

For the calculation of the spontaneous emission, we numerically evaluate the semiconductor luminescence equations, i.e., Eqs. (31)–(33) of Ref. 52. From these, the total spontaneous emission can be calculated by integrating over the luminescence spectra [Eq. (4) of Ref. 53]. Finally, we evaluate the equations yielding the Auger rate in the form of Eqs. (9) and (10) of Ref. 53. In addition to a summation over the k-points, a fourfold sum has to be evaluated such that all possibilities for the initial and final states are considered. The required input for the energy dispersion of each band and the Coulomb matrix elements are taken from the DFT calculations.

More details on the underlying microscopic many-body approach can be found in the supplementary material and in Ref. 50 as well as the references therein. Unless otherwise stated, we use $ε$_∞ = 6.5 from Ref. 20 and $n0≈ε∞=2.55$ which is in good agreement with Ref. 54 for the calculation of the optical properties.

Our DFT calculations yield a Γ-point band gap of E_g = 1.48 eV. This value is significantly smaller than the value E_g = 1.84 eV found by Tao, Cao, and Bobbert¹³ who also applied the LDA-1/2 method. To some degree, we attribute this difference to the fact that Tao, Cao, and Bobbert¹³ obtained larger lattice constants and additionally further increased these to approximate experimental values prior to their band structure calculation. To test this hypothesis, we performed calculations with the lattice constant of Tao, Cao, and Bobbert,¹³ finding that this only accounts for part of the difference between their and our results. An additional contribution could come from the different cutoff parameters obtained in the LDA-1/2 approach, especially for Pb.

To account for the fact that the DFT calculations are at 0 K, we extrapolate the temperature dependent band gaps reported in Ref. 10 for the tetragonal phase to 0 K finding E_g(0 K) = 1.55 eV which is in good agreement with our result. For the microscopic many-body calculations at different temperatures, we shift the conduction bands such that the resulting temperature dependent band gaps correspond to the experimental values, e.g., 1.64 eV at 300 K.¹⁰ 

The top of Fig. 1 presents the computed band structure of the tetragonal phase of MAPbI₃ prior to shifting the conduction bands for the A-Γ-X-M-Γ path and relative to the Fermi level E_f. For clarity, we only depict the two lowest conduction and two highest valence bands that were taken into account in the absorption calculations. The bottom part of the figure shows the material absorption calculated using the data from DFT calculations along the ΓA (blue), ΓX (red), and ΓM (green) directions from top to bottom, respectively. The color of the absorption spectra in the bottom part of Fig. 1 is identical to the color of the band structure for the corresponding direction in the top part of the figure. Both, the absorption spectra including Coulomb interaction and the free-particle absorption spectra have been calculated. For better visibility, a y-offset is added to the absorption for the ΓX and ΓA directions. An especially strong absorption for the ΓM direction is evident, whereas the absorption for the ΓA and ΓX directions is of about equal strength. The actual strength of the absorption depends on the orientation of the crystal in the laser field. According to our results, the absorption at the 1s exciton peak is approximately between 1 µm⁻¹ and 3.5 µm⁻¹ which agrees with the results summarized by Shirayama et al.⁵⁵ (2–5 µm⁻¹). Taking more bands into account does not change the absorption in the energy range shown here (E_g ± 164 meV) as further transitions are at higher energies.

The position of the exciton 1s peak in Fig. 1 and its height depend on the choice of the carrier-phonon induced dephasing constant. For Fig. 1, the dephasing was chosen such that the peak is barely visible. In fact, in experimental measurements around 300 K, the exciton peak is absent, while it emerges for lower temperatures.^10,56 This is most likely due to the reduced phonon scattering of carriers at lower temperatures.

Deducing the exact position of the exciton 1s peak bears the problem that the peak is located on the rising shoulder of the absorption. By choosing an unrealistically low dephasing constant γ = 0.4 meV, the exciton peak becomes much larger and the position of the maximum can be accurately determined. We find that the energy of the maximum is not only slightly different for the different directions but it also depends on the choice of the high-frequency dielectric constant $ε$_∞. Several values for $ε$_∞ have been reported for MAPbI₃ in the literature ranging from 5.0 to 7.1.^{5,12,14,18,20,21} For the ΓM direction, we compute an exciton binding energy of E_b = 65 meV for $ε$_∞ = 5 taken from Ref. 18, E_b = 31 meV for $ε$_∞ = 6.5 taken from Ref. 20, and E_b = 25 meV for $ε$_∞ = 7.1 taken from Ref. 14. Table I gives a summary of the obtained exciton binding energies for five different directions in reciprocal space and three different values of $ε$_∞. With the exception of the highest value, all are well in the range of the binding energies summarized in Ref. 25. The direction dependence of the exciton binding energy is due to anisotropic effective masses.¹⁴ 

In the case of CsSnX₃ (X = Cl, Br, I), Huang and Lambrecht⁵⁷ argued that phonons contribute to the screening and the static dielectric constant $ε$₀ instead of the high-frequency dielectric constant should be used. Whether this is really the correct choice for MAPbI₃ is still under debate.²⁵ While Bokdam et al.⁵ argued that screening due to phonons can be neglected because the energies of the longitudinal phonon modes in MAPbI₃ are by at least a factor of 5 smaller than their calculated E_b = 45 meV, the authors of Ref. 19 claimed that screening due to phonons is relevant. The reason for their argument is the large exciton radius of 204 Å reported by Frost et al.,⁵⁸ which, however, is calculated from a large dielectric constant (25.7¹⁵). In Ref. 19, the authors find E_b ≈ 2 meV.

Here, we estimate how our computed exciton binding energies are modified by the inclusion of significant phonon screening. Since for the exciton binding energy E_b ∝ $ε$⁻² applies, our results obtained from using the high-frequency dielectric constant need to be down-scaled by a factor of $(ε0/ε∞)2$⁠.⁵⁷ The actual magnitude of this factor mostly depends on the choice of $ε$₀ for which different values are reported in the literature, e.g., 70,¹⁹ 33.5,¹⁸ and 25.7.¹⁵ Here, we use the values reported in Ref. 18 (⁠$ε$_∞ = 5; $ε$₀ = 33.5) and find an exciton binding energy of about 1 meV. This is close to the value reported by Lin et al.¹⁹ and even closer to the result of Frost et al.⁵⁸ (0.7 meV).

Figure 2 shows numerical results of the absorption for different carrier densities n_c using the DFT results in the ΓM direction. We see that with increasing carrier densities, the absorption is shifted to higher energies and around n_c ≈ 10 × 10¹⁸ cm⁻³, we obtain a certain amount of negative absorption, i.e., optical gain.

In Fig. 3, we show the computed PL for different carrier densities where the values of the spectra were divided by the square of the respective carrier density for better visibility. For very low carrier densities, we observe that the PL maximum is close to the exciton 1s peak position and, with increasing carrier density, a blue shift is clearly visible. By integrating over the spectra, the rate of carrier loss due to the spontaneous emission is obtained and from this, the coefficient B for the spontaneous emission losses is deduced.

The obtained B values are plotted in Fig. 4 as a function of the carrier density. For the ΓA direction, B is shown in blue for the temperatures T₁ = 170 K, in red for T₂ = 235 K, and in yellow for T₃ = 300 K (all solid lines and filled areas). We clearly see that B increases with decreasing temperature, which is in accordance with Ref. 32. The reason for this increase can be attributed to the temperature dependence of the occupation probabilities at a given density for states near the band gap. The computed values for B are also shown in Fig. 4 for the other directions at 300 K where it is largest for the ΓM direction (dotted line), while it is smallest for the ΓZ direction (dashed-triple-dotted). This direction dependence of B is caused by the anisotropic effective masses of MAPbI₃.¹⁴ For increasing densities, B decreases due to phase space filling. Except for very high carrier densities, our values for B are all in the range of those obtained by experimental measurements^28–32 where the reported carrier densities are all well below 10¹⁹ cm⁻³.

In our calculations for C, twelve conduction bands and 36 valence bands were considered which correspond to E_CBM + 1.82 eV and E_VBM − 1.83 eV at the Γ-point, where E_CBM is the conduction band minimum and E_VBM is the valence band maximum. Including more valence bands only slightly increases the Auger coefficient without causing any other modifications of the results (see Fig. 1 of the supplementary material). The next higher conduction bands are energetically too far away to contribute. The complete k-path from Γ to the edge of the first Brillouin zone is included in the computation of the Auger coefficient. C is obtained for three temperatures, namely, T₁ = 170 K, T₂ = 235 K, and T₃ = 300 K, and plotted as a function of the carrier density n_c in Fig. 5.

With the exception of the ΓZ direction (gray line), the Auger coefficients for all directions are very similar in the low density limit. Generally, C in the ΓX direction (green line) is somewhat larger as compared to the other directions. For increasing densities, more and more states are filled and eventually resonances in the Auger coefficient become visible. Due to charge carrier screening of the Coulomb matrix elements, the energetic positions of the bands are shifted with increasing densities [see Eq. (4)]. At those densities where resonances in the Auger coefficient are visible, the shift is such that it allows carriers to make a transition to higher bands without a momentum transfer in k-space. Since for the five directions considered here, the effective masses differ, these resonances show up at different carrier densities. With increasing temperatures, the resonances are shifted to even higher densities (see Fig. 1 of the supplementary material). While at low carrier densities, the Auger coefficient increases with temperature, it is almost equal for the three temperatures shown from roughly n_c ≈ 5 × 10¹⁷ cm⁻³ onto higher densities. All in all, these results are in good agreement with experimental findings that reported Auger coefficients in the range of 10⁻²⁸ cm⁶/s to 10⁻²⁷ cm⁶/s.^29–31 At carrier densities that are common for solar cell applications on the earth’s surface,²⁵ we find an Auger coefficient in the range of 1–5 × 10⁻²⁷ cm⁶/s.

Independent of temperature and carrier density, we find that the loss current due to the spontaneous emission is usually larger by at least one order of magnitude than that due to the Auger recombination (see Fig. 2 of the supplementary material). An exception is the conditions for which some of the resonances in the Auger rates occur which can enhance the corresponding loss current at that particular carrier density such that it becomes larger than the loss current due to the spontaneous emission.

In conclusion, we used an ab initio based approach combining density functional theory calculations for the electronic properties and microscopic many-body computations to obtain optical properties of MAPbI₃. We find a strong direction dependence for the absorption which is largest for the ΓM direction. Our calculations yield exciton binding energies (21 meV–65 meV) well in the range of reported values that depend on the choice of the high-frequency dielectric constant. Additionally, we compute the reduction in the exciton binding energy in the case when screening due to phonons needs to be accounted for.¹⁹ Furthermore, we analyze the optical absorption for varying carrier densities. We find a blue shift of the absorption with increasing density and identify the onset of the gain region. This blue shift is also evident in the results of our photoluminescence calculations. Finally, we investigated the magnitude of the intrinsic loss processes for different temperatures by calculating the spontaneous and Auger recombination coefficients, finding a reasonably good agreement with existing experimental data. We extend the experimental results by also calculating the dependence on the carrier density. This is especially strong in the case of the Auger coefficient C. For the complete range of densities used (10¹⁵–10¹⁹ cm⁻³), we note that the loss current due to the spontaneous emission is by at least one order of magnitude larger than the loss current due to the Auger recombination, independent of temperature or direction.

See supplementary material for further details on the quantum many-body theory and additional results for the Auger losses.

The Marburg work was funded by the DFG via the GRK 1782 “Functionalization of Semiconductors”; computing time from HRZ Marburg, CSC Frankfurt, and on the Lichtenberg high performance computer of the TU Darmstadt is acknowledged. The Arizona work was supported by the Air Force Office of Scientific Research under Award No. FA9550-17-1-0246.