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

For opto-electronic and photo-voltaic applications of perovskites, it is essential to know the optical properties and intrinsic losses of the used materials. A systematic microscopic analysis is presented for the example of methylammonium lead iodide where density functional theory is used to calculate the electronic band structure as well as the dipole and Coulomb matrix elements. These results serve as input for a many-body quantum approach used to compute the absorption, photoluminescence, and the optical and Auger losses for a wide range of application conditions. To illustrate the theory, the excitonic properties of the material system are investigated and numerical results are presented for typical photo-voltaic operation conditions and for the elevated carrier densities needed for laser operation.

For opto-electronic and photo-voltaic applications of perovskites, it is essential to know the optical properties and intrinsic losses of the used materials. A systematic microscopic analysis is presented for the example of methylammonium lead iodide where density functional theory is used to calculate the electronic band structure as well as the dipole and Coulomb matrix elements. These results serve as input for a many-body quantum approach used to compute the absorption, photoluminescence, and the optical and Auger losses for a wide range of application conditions. To illustrate the theory, the excitonic properties of the material system are investigated and numerical results are presented for typical photo-voltaic operation conditions and for the elevated carrier densities needed for laser operation.
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][2][3] . The perovskites of interest here can be summarized by the chemical formula AMX 3 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 CH 3 NH + 3 or formamidinium (FA) which is CH(NH 2 ) + 2 but there are also other candidates such as azetidinium or arizidinium 4 . 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 5 where the choice of the halide atom has the strongest impact on the actual band gap value. Depending on the temperature and the composition of the respective perovskite, the crystal structure is cubic, tetragonal or orthorhombic. For changing system temperatures, interesting phase transitions have been observed [6][7][8][9] .
In the current study, we investigate the optical properties of the tetragonal phase of methylammonium lead iodide (MAPbI 3 ). This material system has a band gap of approximately 1.64 eV at room temperature 10 and is one of the most intensively analyzed perovskites for solar cell applications. The band structure of MAPbI 3 was previously studied using density functional theory (DFT) [11][12][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][12][13][14][15]17 and slightly a) Electronic mail: lars.bannow@physik.uni-marburg.de shifts the conduction band minimum away from the Γ point of the Brillouin zone such that the optical transitions become somewhat indirect 15 . More advanced DFT calculations utilizing the self-consistent PBE0 11 or LDA-1/2 13 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 3 in order to understand the details behind the good energy conversion in these materials. A significant range of values has been reported for the static and the high-frequency dielectric constants 5,12,14,15,[18][19][20][21] . Consequently, 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][23][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 results in a severe reduction of the exciton binding energy 5,19,22 , whereas 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. 5 who predict 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, the operation of optoelectronic devices are 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 3 , yielding spontaneous emission coefficients in the range of 6·10 −11 cm 3 /s to 2·10 −9 cm 3 /s [28][29][30][31][32] and Auger coefficients in the range of 10 −28 cm 6 /s to 10 −27 cm 6 /s [29][30][31] . A comparison between the involved bands of the Auger processes in hybrid perovskites and those in 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 33 calculations for the electronic properties of MAPbI 3 with a microscopic many-body approach to compute optical properties such as carrier density dependent absorption and photoluminescence spectra as well as the temperature and density dependent spontaneous emission and Auger loss coefficients. This approach has previously been successfully used 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][37][38][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 42 . For the lattice optimization the local density approximation (LDA) as parametrized by Perdew and Zunger 43 was utilized, converging the Hellmann-Feynman forces on the atoms below 10 −3 eV/Å and the energies to an accuracy of 10 −6 eV. An energy-cutoff for the plane waves of 800 eV was chosen. Optimization of the atom positions, cell shape and volume were performed from a perfect tetragonal crystal structure as a starting point. This results in a pseudo-tetragonal lattice with a = 8.540Å , b = 8.546Å and c = 12.611Å which are somewhat smaller than the experimental lattice constants 8 and the lattice vectors are slightly non-orthogonal. We omitted the inclusion of van-der-Waals corrections based on the recently reported poor performance of Van-der-Waals functionals for the description of the perovskite lattice 44 . 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 −4 eV and spin-orbit coupling was enabled. In addition, the LDA-1/2 method 45 was used as this method was shown to yield band gaps that are close to the ones obtained with GW for perovskites 13 . The advantage is that the LDA-1/2 method is computational equally demanding to using the LDA. It is based on Slater's half occupation scheme 46 and aims at an improved inclusion of the electron self-energies. For iodine we find a cutoff R I CUT = 3.71 a.u. that maximizes the band gap which is similar to the value Tao et al. 13 report (3.76 a.u.). For lead, we obtained R Pb CUT = 2.00 a.u. while Tao et al. 13 found R Pb CUT = 2.18 a.u. Yin et al. 17 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 47 , 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 ΓA direction ([111] where f 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 k k and V k k k−q q q , respectively. Furthermore, Ω k k k is given by where the first term describes the interband dipole coupling to the laser field E(t) and the second term describes again a re-normalization caused by the Coulomb interaction V k k k−q q q . Here, the dipole matrix elements d k k k and the Coulomb matrix elements V k k k−q q q obtained from DFT calculations enter. In our current study, we investigate Coulomb effects. While phonon-carrier scattering in principle can be included in the scattering term of Eq. 1, this is challenging due to the strong phonon coupling present in the perovskite materials 5 and the multitude of values reported for the dielectric constants and LO phonon frequencies as discussed above. Therefore, the dephasing caused by scattering processes involving phonons and carriers is approximated by a phenomenological dephasing constant of γ = 165 meV (equivalent to a dephasing of 25 fs) and the scattering term takes the form −iγp k k k . Additionally, we include carrier-carrier scattering in 2nd Born approximation in the scattering term and Coulomb-screening is taken into account whenever Coulomb matrix elements enter the above equations by evaluation of the Lindhard formula.
For the calculation of the photoluminescence the semiconductor luminescence equations are evaluated which are derived in Ref. 50 [Eqs. (31)-(33) therein]. From these, the spontaneous emission can be calculated by integrating over the luminescence spectra [Eq. (4) of 51 ]. Finally, the equations yielding the Auger rate can be found in Ref. 51 [Eqs. (9) and (10) therein]. 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 Ref. 52 and references therein. Unless otherwise stated, we use ε ∞ = 6.5 from Ref. 20 for the calculation of the optical properties.
Our DFT calculations yield a Γ-point band gap of E g = 1.48 eV which is significantly smaller than the band gap found by Tao et al. 13 (E g = 1.84 eV) who also applied the LDA-1/2 method. The reason is that Tao et al. 13 re-scale their lattice constants to approximate experimental values. 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 bandgaps correspond to the experimental values, e.g. 1.64 eV at 300 K 10 .
The top of Fig. 1 shows the obtained band structure of the tetragonal phase of MAPbI 3 prior to shifting the conduction bands for the A-Γ-X-M-Γ path and relative to the Fermi level E f . Only the two lowest conduction and two highest valence bands that were taken into account in the absorption calculations are shown. 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 shown in the bottom part of Fig. 1 is identical with 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 direction. An especially strong absorption for the ΓM direction is evident whereas the absorption for the ΓA and ΓX direction are 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 −1 and 3.5 µm −1 which agrees with the results summarized by Shirayama et al. 53 (2 − 5 µm −1 ). 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.
From the onset of the free-particle absorption which corresponds to the band gap E g and the position of the exciton peak, the exciton binding energy can be com- puted. The position of the peak and its height depend on the choice of the 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,54 . This is most likely due to the reduced phonon scattering at lower temperatures. The position of the maximum can be obtained accurately by using a smaller dephasing constant. It is slightly different for all directions and it depends on the choice of the highfrequency dielectric constant ε ∞ . Several values for ε ∞ have been reported for MAPbI 3 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 14 .
In the case of CsSnX 3 (X = Cl, Br, I), Huang and Lambrecht 55 argue that phonons contribute to the screening and the static dielectric constant ε 0 instead of the highfrequency dielectric constant is to be used. Whether this is the case for MAPbI 3 or not is still under debate 25 . While Bokdam et al. 5 argue that screening due to phonons can be neglected because the energies of the longitudinal phonon modes in MAPbI 3 are by about a factor of 5 smaller than their calculated E b = 45 meV, the authors of Ref. 19 claim that screening due to phonons is relevant. The reason for their argument is the large exciton radius of 204Å reported by Frost et al. 56 which, however, is calculated from a large dielectric constant (25.7 15 ). 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 ∝ ε −2 applies, our results obtained from using the high-frequency dielectric constant, need to be downscaled by a factor of (ε 0 /ε ∞ ) 2 . 55 The actual magnitude of this factor mostly depends on the choice of ε 0 for which different values are reported in the literature, e.g. 70 19 , 33.5 18 and 25.7 15 . Here, we use the values reported in Ref. 18 and find an exciton binding energy of about 1 meV. This is close to the value reported by Lin et al. 19 and even closer to the result of Frost et al. 56 (0.7 meV). Figure 2 shows the numerical results for different carrier densities n c using the DFT results in ΓM direction. We see, that with increasing carrier densities the absorption is shifted to higher energies and around n c ≈ 10 · 10 18 cm −3 , 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 blueshift is clearly visible. By integrating over the spectra, the rate of carrier loss due to 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 1 = 170 K (dotted  14,18,20 . line), T 2 = 235 K (dashed line) and T 3 = 300 K (solid line). 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 bandgap. 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 (orange) while it is smallest for the ΓZ direction (grey). For increasing densities, B decreases due to phase space fill- ing. Except for very high carrier densities, our values for B are all in the range of those obtained by experimental measurements [28][29][30][31][32] where the reported carrier densities are all well below 10 19 cm −3 . In our calculations for C, twelve conduction bands and 36 valence bands were considered which corresponds 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 having any other effects on the results. 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 1 = 170 K, T 2 = 235 K and T 3 = 300 K and plotted as a function of the carrier density n c in Fig. 5.
With the exception of the ΓZ direction (grey line), in the low density limit the Auger coefficients for all directions are very similar. Generally, C in ΓX direction (green line) is somewhat larger than compared to the other directions. For increasing densities, more and more states are filled and eventually resonances in the Auger coefficient are visible. These resonances occur whenever carriers can make a transition to higher bands without 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 (dotted line → dashed line → solid line; shown for the ΓA direction, i.e. in blue), the resonances are shifted to even higher densities (not shown). While at low carrier densities the Auger coefficient increases with temperature, it is almost equal from roughly n c ≈ 5 · 10 17 / cm 3 on to higher densities. All in all, these results are in good agreement with experimental findings that reported Auger coefficients in the range of 10 −28 cm 6 /s to 10 −27 cm 6 /s. [29][30][31] At carrier densities that are common for solar cell applications on the earth's surface 25 , we find an Auger coefficient in the range of 1 − 5 · 10 −27 cm 6 /s. Independent of temperature and carrier density, we find that the loss current due to spontaneous emission is usually larger by about one order of magnitude than that due to Auger recombination. An exception are 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 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 3 . We find a strong direction dependence for the absorption which is largest for the ΓM direction. Our calculations yield exciton binding energies (21 meV to 65 meV) well in the range of reported values that depend on the choice of the highfrequency dielectric constant. Additionally, we compute the reduction of the exciton binding energy in the case when screening due to phonons needs to be accounted for 19 . 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. Lastly, 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 15 − 10 19 cm −3 ) we note that the loss current due to spontaneous emission is by about one order of magnitude larger than the loss current due to Auger recombination, independent of temperature or direction.
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 numbers FA9550-16-1-0088 and FA9550-17-1-0246.