Understanding the defect chemistry of lead-halide perovskites and its effects on the hot-carrier lifetime is of significance for both fundamental understanding and applications as solar cell light absorbing materials. In this study, the mechanistic details of hot carrier decay in hybrid perovskites are investigated using a newly developed non-adiabatic molecular dynamics method. In this approach, the nuclear trajectory is based on Born–Oppenheimer ground state molecular dynamics, which is then followed by the evolution of carrier wave function including the detailed balance and decoherence effects. We found the longer decay time for hot electrons due to the incorporation of interstitial iodine in the hybrid lead-halide perovskites (MAPbI_{3}), while the hot hole decay time is not affected significantly by the interstitial iodine. The underlying mechanism for such modulation of hot carrier dynamics is attributed to the changes of carrier density of states and the electron–phonon coupling strength. Hence, iodine interstitial is the necessary condition to create long-lived hot electrons in perovskites, which is further demonstrated by the comparative analysis with the pure MAPbI_{3}.

## INTRODUCTION

One of the most intensely studied recent materials for next generation solar cells is the organic–inorganic hybrid perovskite, namely, methylammonium lead iodide (MAPbI_{3}).^{1} Besides possessing many desirable properties for solar cell applications, including ease of synthesis, strong light absorption,^{2} low exciton binding energy,^{3,4} low free carrier recombination rate,^{5} and long carrier lifetime, MAPbI_{3} exhibits unusual excited state hot carrier dynamics.^{6,7} Its hot carrier cooling time (>1 ps) is much longer than most organic semiconductors as well as purely inorganic compounds such as GaAs (*τ*_{e} = ∼0.1 ps).^{8–10} Not only is this intellectually intriguing, such long carrier lifetime could also have implications for practical applications, e.g., raising the possibility for hot carrier collection to realize the hot electron solar cell to go beyond the Shockley–Queisser limit.^{5,11}

Another intriguing property of MAPbI_{3} is its extreme defect tolerance. Being synthesized in a wet chemistry condition, the MAPbI_{3} crystal is believed to be filled with imperfections. Despite such high density of defects, the system still has very long carrier lifetime. In contrast to expectation, MAPbI_{3} exhibits low density carrier traps (∼10^{11} cm^{−3}).^{12,13} Most of the common defects do not become free carrier killer. Adinolfi *et al.*^{12} reported that the traps in single crystalline MAPbI_{3} are distributed within 0.2–0.3 eV above and below the conduction band (CB) and valence band (VB), respectively. It is argued that MAPbI_{3} is inherently defect tolerant due to the fact that dominant defects introduce only shallow traps in the material bandgap.^{14–16} Such native defects with shallow states can also form Schottky complexes leading to self-compensation effects.^{12,17} Among all these native defects, Meggiolaro *et al.*^{18} rationalized that the iodine defect chemistry determines the so-called “defect tolerance” in metal-halide perovskites due to intrinsic charge traps leading to retarded charge recombination. There are also some discrepancy in the literature regarding the exact positions of the defect states. For example, the PBE calculation indicates that a point ‘I’ interstitial defect will form a shallow acceptor near the valence band with a 0/− charged transition state,^{14} while a hybrid Density Functional Theory (DFT) calculation including spin–orbit coupling shows a mid-gap +/− transition state as the lower energy defect state.^{19}

Recent experimental observations on the behavior of MAPbI_{3} upon exposure to I_{2} vapor reveals that the charge and ion transport and trapping activity can be tuned upon exposing MAPbI_{3} to I_{2}.^{20} Recent endeavors have also established that I_{2} spontaneously dissociates on MAPbI_{3}(001) surfaces to form pairs of negative/positive interstitial iodine, and the reverse process of I_{2} release becomes thermodynamically favored upon trapping a hole/electron pair via a photoinduced trap-curing mechanism.^{20} Thus, I_{2}-interstitial (MAPbI_{3.5}) is a reversibly tunable defect, hence an efficient approach to tune the property without any permanent disruption in the perovskite framework. It is also shown that I_{2} vapor directly affects the PL quenching of a perovskite film and consequently influences the power conversion efficiency of hybrid perovskite photovoltaic devices.^{21}

In view of the pertinent role of I interstitial defect in the operation of MAPbI_{3} and its unique property of hot carrier cooling, it is thus of great interest to study the effects of the I interstitial defects to hot carrier relaxation in such a system. Owing to the way to controllable synthesis of such defects, we will call the system MAPbI_{3.5} although it does not mean that the interstitial iodine will still be in a molecular form inside the MAPbI_{3} framework. Despite many experimental and theoretical studies aiming at the understanding of hot carrier cooling in the pure MAPbI_{3},^{7,22–24} there is still a lack of studies of the effects of defects on hot carrier cooling, which is the goal of this work.

The MAPbI_{3} system affords three distinct types of chemical bonding. The first is the interlinked and covalently bonded [PbI_{6}] octahedra, the second is the bond between [PbI_{6}] octahedra and the I-chain, and the third is the weak hydrogen bond involving N–H–I.^{25} The interplay among these three kinds of chemical bonds gives rise to a “soft” fluctuating structure.^{26} Our previous work showed that the carrier mobility of MAPbI_{3} is driven by the dynamic disorder due to the random orientations of the MA cations that causes fluctuation of the electrostatic potential.^{27,28} Consequently, the electron and hole wave functions were localized and spatially separated. The localized carriers are able to conduct random walks, which constitute their carrier mobility at the ground state. However, the mechanism of hot carrier relaxation dynamics and the effect of wave function localization to the carrier relaxation time are still far from being understood, which are impeded by a lack of adequate computational approach to carry out simulations for the non-adiabatic process for sufficiently large system-size and longer timescale.

In this study, we apply a newly developed non-adiabatic molecular dynamics (NAMD) method and identify the role of iodine defects toward the hot-carrier relaxation process. Through our simulation, we reveal that a combined effect of the distorted PbI_{x} polyhedra (*x* ≥ 6) enhanced hydrogen bonding interaction, which in combination with the change of the density of states (DOS) leads to a prolonged hot electron lifetime in MAPbI_{3.5} compared to MAPbI_{3}, while the hot hole lifetime is not changed dramatically.

## METHOD

### First principles simulation details

The calculations are separated into two parts. The first is to investigate the structure and formation energies for MAPbI_{3} added with extra I (MAPbI_{3.5}). The second is to simulate the carrier dynamics in such a system using NAMD. For the formation energy calculations, taking advantage of the compound energies already calculated in the Material Project database, which used the Vienna Ab initio Simulation Package (VASP) code,^{29,30} we have also used the VASP calculation with the same pseudopotential and energy cutoff. More specifically, the projector augmented-wave^{31} pseudopotential is used with a cutoff energy of 520 eV. k-point number multiplied with the number of atoms is set to be greater than 1000. We have examined different configurations of the I_{2} within the MAPbI_{3} framework and have estimated their phase stability to determine the stable structure. To calculate the stability for a given configuration, we have used grand canonical linear programming (GCLP), using the hundreds of known compounds from the materials project.^{32} Basically, we decompose a given system I_{2}:MAPbI_{3} into compounds such as NH_{4}I, CH_{4}, PbI_{2}, and C. We have used enumeration to generate twenty one structural configurations of MAPbI_{3.5} while considering the cutoff distance as 1.39 Å (covalent radius for iodine) to initially construct the structure with interstitial iodine. Ten most stable DFT-optimized MAPbI_{3.5} structures and their relative stabilities are shown in Fig. 1 of the supplementary material. We selected the one with the lowest formation energy and used that for the following NAMD simulation.

We have used the PWmat package^{33,34} for NAMD simulation. The generalized gradient approximation (GGA/PBE) exchange correlation functional^{35} is used, and the atomic structure is relaxed prior to molecular dynamics (MD) simulations. The PWM pseudopotentials^{36} is used with a 60 Ryd plane wave kinetic energy cutoff. The total energy vs energy-cutoff plot to check the convergence is provided in Fig. 3 of the supplementary material. This pseudopotential yields similar results as that of the SG15 pseudopotential^{37} for both the atomic structure and electronic structure. The comparison of the PWM pseudopotential and the SG15 pseudopotential is provided in Fig. 5 of the supplementary material.

Carrier (electron/hole)–phonon coupling constant (CP) is calculated for a given pair of input electronic states for all the phonon modes as (*i*_{1} and *i*_{2} are the wave function for state index i, k-point k)

*i*_{1} is the CBM (VBM) state to compute the electron (hole)–phonon coupling. *i*_{2} is the near band-edge 20 adiabatic conduction band (valence band) states. CP_{i1−i2} represents a perturbed atomic force, which can also be decomposed into the contribution from each atomic site. Each point in Fig. 5 represents the averaged CP for 20 adiabatic states (*i*_{2}) near the band edge,

### NAMD method details

There are two different ways to theoretically study hot carrier cooling. One is to explicitly calculate the electron–phonon coupling and then use Boltzmann equation and Fermi golden rule to calculate the carrier relaxation.^{38,39} This approach is good for problems well described by the harmonic approximation for phonon oscillation and there is no multiple phonon process. For the MAPbI_{3} system, however, the “soft” fluctuation and rotation of the MA molecule might be well beyond the harmonic approximation. In such cases, one can use the second method to study hot carrier cooling, which is a direct NAMD simulation. It can handle the anharmonicity and multi-phonon process, although classical ways are used to describe the phonon movement. One challenge is the finite size of the simulated system. As we will test later, for the MAPbI_{3} system, a few hundred atom system seems good enough.

In our NAMD simulation, we first perform an *ab initio* molecular dynamics (BOMD) simulation. The plots in Fig. 4 of the supplementary material show the energy and temperature evolution during the adiabatic molecular dynamics simulations. The nuclear trajectory from BOMD is then used to perform the time-dependent Schrödinger’s equation (TDSE)

for the time evolution of the wave function. In doing so, we have ignored the back reaction from the electron movement to the nuclear movement. This so called classical path approximation (CPA) is good for large systems where there is no strong polaronic effect,^{40} and we care more about the electronic behavior rather than the nuclear behavior. In contrast with the existing NAMD methods, our time-dependent Schrödinger’s equation is reformulated using a density matrix formalism (called the P-matrix method).^{41}

The hot-carrier wave function [*ψ*_{l}(*t*)] is expanded in terms of the adiabatic basis *ϕ*_{i}(*t*) as follows: *ψ*_{l}(*t*) = ∑_{i}*c*_{il}*ϕ*_{i}(*t*), where the index “*l*” stands for the individual cases in a statistical assembling. The density matrix of the system, *D*_{ij}(*t*), under the basis of *ϕ*_{i}(*t*) is then defined as

where $wl$ is the statistical weight of *ψ*_{l}.

The charge density *ρ*(*r*) of the hot hole/electron is calculated by *ρ*(*r*) = ∑_{ij}*D*_{ij}*ϕ*_{i}(*r*)*ϕ*_{j}(*r*)^{*}. When *t* → *∞*, one can expect that the average population of hot carrier will be diminished to zero at the high energy state and it will be one at the band edge, CBM (VBM). To analyze the charge density redistribution while hot carrier with energy 0.5 eV decays down to the band edge, we have plotted *ρ*(*r*) at the initial and final NAMD steps, respectively.

Applying the TDSE, the equation of motion for the density matrix can be expressed as

and

where *V*_{ij} implicitly includes the effect of the electron–phonon coupling. The second term in Eq. (5) contains the decoherence time *τ*_{ij}(*t*) between state i and j. To overcome the problem of classical description of the phonon movement, the P-matrix formalism includes the detailed balance via splitting the density matrix *D*: *D* = *P* + *P*^{T}. *P*_{ij} represents the transition of electronic state population from state i to j, where $Pij\u2260Pji*$. The diagonal elements of the density matrix (*Dii* = 2*P*_{ii}) evolves as^{41}

The last two terms in Eq. (7) introduce the detailed balance, where Δ*ϵ*_{ij} = *ϵ*_{i} − *ϵ*_{j}. When the hot electron is considered as carrier, *f*_{ij} = 1(0) for Δ*ϵ*_{ij} > 0 (Δ*ϵ*_{ij} < 0). On the other hand, for the hot hole, *f*_{ij} = 0(1) for Δ*ϵ*_{ij} > 0 (Δ*ϵ*_{ij} < 0).

The off-diagonal element of the P matrix evolves as

The last term in Eq. (8) introduces the decoherence effect.

Overall, in this formalism, the detailed balance and decoherence effect are accounted. The decoherence influences the cooling rate, and the detailed balance describes the cooling process. Due to the use of density matrix, one single time evolution can account for an statistical assemble, including many branching (collapsing) of the wave function trajectory. This makes the simulation efficient without the need for repeated stochastic simulations. The parameter *τ*_{ij} in Eq. (8) is estimated explicitly by the fluctuation of the Δ*ϵ*_{ij} as described in Ref. 41 and subsequently being used to get the hot carrier energy decay time. Nevertheless, Fig. 6 of the supplementary material shows the hot carrier energy evolution (same initial condition) vs time for the MAPbI_{3.5} system while considering different decoherence time (dct). If one uses a high decoherence time (dct = 1000), hot carrier cooling will be slower. On the other hand, if one uses an extremely small decoherence time (dct = 1), hot carrier cooling will be faster.

To quantify the cooling time of hot carriers, we fitted the energy decay curves by a function consisting of an exponential and a Gaussian function as used in Refs. 24 and 42,

Here, *dE*(*t*) is the average energy of one assemble simulation above (or below) the band edge. Then, the effective cooling time for hot carriers is estimated as decay time constant (*τ*) as

Typically, “a” is close to 0.9. The fitting parameters are tabulated in Fig. 7 of the supplementary material.

## RESULTS & DISCUSSION

We employ a tetragonal MAPbI_{3} structure [see Fig. 1(a–i)] as the model framework with 4% I_{2}-interstitial (MAPbI_{3.5}) to study the effect of iodine-interstitial defect physics toward the hot carrier relaxation process. After testing many initial configurations followed by atomic relaxations, the stable configuration of MAPbI_{3.5} exhibiting the smallest Δ*E*_{decomposition} is determined from the convex hull analysis and shown in Fig. 1(a–ii). Such a configuration has two characteristics: (1) PbI_{x} (x > 6) polyhedra are formed along with the regular PbI_{6} octahedra; and (2) isolated I-interstitial in the framework forming linear I_{4} chains, combining two ‘I’ from PbI_{x} octahedra. From the structural overview [see Fig. 1(a–ii)], the MAPbI_{3.5} crystal structure exhibits three distinct types of chemical bonding: (1) covalently bonded [PbI_{x}] polyhedra (x ≥ 6); (2) the secondary bonds between the [PbI_{x}]-chains and I_{2} bridging molecules, forming the I_{4} chain; (3) N–H–I hydrogen bonds. Similar to MAPbI_{3}, MAPbI_{3.5} also exhibits phase-instability with respect to its decomposition into CH_{4}, NH_{4}I, NH_{4}I_{3}, PbI_{2}, and C with a decomposition energy (Δ*E*_{decomposition}) of 55 meV/atom. This is comparable to the decomposition energy of 40 meV/atom for MAPbI_{3} itself (which does not prevent the experimental synthesis of MAPbI_{3}). Nevertheless, favorable bonding interaction between the excess I_{2} and MAPbI_{3} leads to the negative formation energy (−0.47 eV/atom) of MAPbI_{3.5}.

We note that the above structural configuration of MAPbI_{3.5} with a finite percentage of ‘I’ doping is different from the isolated ‘I’ interstitial defect. The partial density of states (PDOS) of MAPbI_{3.5} is shown in Fig. 1(b), in comparison with the PDOS of MAPbI_{3} in Fig. 1(a). The electronic structure calculation shows that there is no obvious defect state in the middle of the bandgap (see Fig. 1). Figure 2 of the supplementary material shows the resemblance of the band-edge states of HSE and PBE for our defect structure of MAPbI_{3.5}, while conduction band minima wave functions are concentrated near the iodine centers in a linear fashion. However, since there is no isolated iodine defect, the wave function is not highly localized as in the isolated iodine defect case. For the conduction band, the peak in MAPbI_{3.5} DOS is reduced near the band edge and contributed by the I_{4}-species, reducing the bandgap of MAPbI_{3.5} compared to MAPbI_{3}. On the other hand, the valence band-edge is mostly consisting of the I-p state for both MAPbI_{3} and MAPbI_{3.5}. Thus, the overall elemental characteristics for the valence band, especially the I-p/Pb-p contribution, are not changed dramatically from MAPbI_{3} to MAPbI_{3.5}. In sharp contrast, the conduction band-edge is mostly consisting of I_{4}-species for MAPbI_{3.5}, while the Pb-p state predominantly constitutes the same for MAPbI_{3}. The NAMD wave functions, as will be discussed subsequently, show similar localization for the conduction band-edge states where the localization is centered around the I_{4} chains with strong I-p contributions and the typical excited states (in both conduction and valence bands). Note that the PDOS shown in Fig. 1(b) are calculated based on a single snapshot during the BOMD for both MAPbI_{3} and MAPbI_{3.5}; thus, there could be some fluctuations depending on the particular time used to calculate the PDOS.

We have used the above-mentioned MAPbI_{3.5} (MAPbI_{3}) system with 400 (384) atoms (2 × 2 × 2 supercell) for further NAMD simulations. We find that the hot carriers in MAPbI_{3} with an excess energy of ∼0.5 eV above the band-edge decays has a ∼0.4 (∼0.6) picosecond lifetime (*τ*) for the hole (electron), as shown in Fig. 2. This timescale for decay agrees with the transient absorption spectroscopy measurement by Xing *et al.*,^{43} which shows that hot-hole cooling in MAPbI_{3} at 480 nm laser excitation has a lifetime of ∼0.4 ps. We note that the hot carrier excitation of MAPbI_{3} has already been simulated with similar non-adiabatic methods,^{24,44} and similar hot carrier cooling time was found. Here, to test whether our supercell is large enough, we have increased the supercell size to 4 × 2 × 2 consisting of 768 atoms compared to the original 2 × 2 × 2 supercell of the MAPbI_{3} system with 384 atoms and repeated the NAMD simulation using the PWM pseudo-potential. A comparative plot for the time evolution of hot electron energy is shown in Fig. 5 of the supplementary material. We found similar decay time in both the cases, indicating the size-consistency of the system in our simulation. This might be due to the relatively large atomic distortion of the system caused by the soft phonon modes. Such large distortion makes the system closer to a random system, avoiding the requirement of momentum conservation in a typical crystal, and usually makes the system-size convergence much faster.

For the MAPbI_{3.5} system, our NAMD results show a significantly prolonged lifetime for electron decay (*τ*_{e} = 1.5 ps), while only a modest increase for the hole decay (*τ*_{h} = 0.6 ps) compared to the pure MAPbI_{3} (see Fig. 2). To test the effect of MA fluctuation, we have artificially slowed down the movement of MA by increasing the masses of H, N, and C atoms by a factor of two. Notably, slowing down the MA-ion motion in pure MAPbI_{3} leads to an increase in the hot-carrier lifetime (see Fig. 2) for both electrons and holes, showing that the movement of MA plays an important role for the carrier cooling.

We have also studied the dependence of the hot carrier decay on the initial hot carrier energy away from the band edge. The results are shown in Fig. 3. In Figs. 3(a) and 3(b), the average hot carrier energies as functions of times are shown as the red curve, while individual instantaneous adiabatic states are shown in blue thin lines. Using Eq. (9), we obtained the hot carrier decay time, and the results are shown in Fig. 3(c). As we can see, the decay time varies significantly for the hot electron depending on the initial excitation energy. The higher the excited energy, the faster the initial decay for hot-electrons. This is probably due to the higher density of states for the electron at higher energy. On the other hand, the decay time for the hole is more or less unchanged as a function of the excitation energy.

To understand the microscopic origin of these two findings, i.e., (1) prolonged carrier lifetime and (2) difference between hot-electron and hot-hole characteristics, we have analyzed the charge density [*ρ*(*r*)] redistribution during the hot-carrier decay process (see Fig. 4), which reveals that there is a significant contribution from the iodine states for the CB-wave function for MAPbI_{3.5}, unlike the predominant Pb-contribution in the CB-wave function for pure MAPbI_{3}. As consistent with the PDOS plot, it is clear from the individual wave function plots. Thus, there are some major changes for the CB wave function characteristics, especially for the states near the conduction band minimum (CBM). This is partially responsible for the dramatic increase in the electron relaxation time from MAPbI_{3} to MAPbI_{3.5}. Another difference is the localization of the CB wave functions. Due to the existence of I_{4} chains, there is strong localization of the CB states localized in such chains [see Fig. 4(a–iv)]. This is probably is another reason why the hot electron relaxation in MAPbI_{3.5} becomes much slower. As for evidence, such slow down becomes more prominent when the carrier is closer to the conduction band edge, since there the wave function becomes more localized.

In contrast with the electron, the hot hole relaxation in MAPbI_{3} and MAPbI_{3.5} is similar. This can be understood from the PDOS of these two systems. As shown in Fig. 1(b), there is a more prominent peak in the PDOS of MAPbI_{3.5}, compared to that of MAPbI_{3} near the VB edge. However, if we focus on the band edge, the density is state is similar. A more detailed wave function plot in Fig. 4(b) also shows that the typical wave functions are similar for these two systems. There is no strong spatial localization for the hole states, although the I-p is the main component of the hole state.

To further understand hot carrier cooling, we have calculated the electron–phonon coupling for all the phonon modes and for the coupling between the CBM (VBM) states and the higher (lower) energy states. The results are shown in Fig. 5. As we can see, although the hole–phonon coupling is similar for both MAPbI_{3} and MAPbI_{3.5}, the electron–phonon coupling between the CB states is much smaller in MAPbI_{3.5} than in MAPbI_{3}. This, together with the change of DOS, is directly responsible for the slowing down the carrier relaxation process in the conduction band. The reduction in the electron–phonon coupling is probably a consequence of the spatial localization of electron wave functions in MAPbI_{3.5}.

## CONCLUSION

In summary, it is interesting to note that the incorporation of molecular iodine (I_{2}) and its partial dissociation within the pure hybrid perovskite framework introduces linear I_{4} moieties within the perovskite cage framework. Interstitial iodine also leads to hyper-coordination for otherwise perfect Pb-octahedra, leading to distortion in the inorganic cage of MAPbI_{3.5}. Such geometric factors and associated electronic structures enhance the lifetime (*τ*) for hot electrons (1.5–2.5 ps) compared with the pure MAPbI_{3} (0.6 ps). On the one hand, the hot-hole relaxation time is only modestly increased from the original 0.4 ps to 0.6 ps. Upon incorporation of I_{2}, the change in the conduction band is mainly due to the change in the electron state localization and the increase in I-p characteristics. All these lead to a smaller electron–phonon coupling constant, hence the increase in hot electron relaxation time. On the other hand, the extent of modulation for the hole states from MAPbI_{3} to MAPbI_{3.5} is relatively small. The wave functions are similar without strong localization, and the hole–phonon coupling is also similar. As a result, the hot hole relaxation time does not change significantly. The exceptionally long-lived excess electronic energy may make hot-carrier harvesting possible in MAPbI_{3.5}.

## SUPPLEMENTARY MATERIAL

See the supplementary material for a comparative analysis of PBE and HSE methods, total energy vs energy cutoff plot, time evolution of total energy and temperature of MAPbI_{3.5} during the adiabatic molecular dynamics simulations, convergence test for supercell size, time evolution of the energies of hot electrons and hot holes with different E_{ex} for MAPbI_{3.5}, and relaxation time (*τ*) and fitting parameters.

## ACKNOWLEDGMENTS

S.B. and X.Z. were supported by the National Key R&D Program of China (Grant No. 2016YFB0700700) and the National Natural Science Foundation of China (Grant Nos. 11774239 and 61827815). L.W.W. was supported by the Director, Office of Science (SC), Basic Energy Science (BES)/Materials Science and Engineering Division (MSED) of the U.S. Department of Energy (DOE) under Contract No. DE-AC02-05CH11231 through the Theory of Material project. We acknowledge the computational resources of the National Energy Research Scientific Computing Center (NERSC), Office of Science of the U.S. Department of Energy.

## REFERENCES

_{3}NH

_{3}PbI

_{3}perovskite solar cell absorber

_{3}

_{3}NH

_{3}PbI

_{3}: Effects of spin–orbit coupling and self-interaction error

_{3}: A probe of lead-halide perovskites defect chemistry

_{2}molecules and weak interactions in supramolecular assembling of pseudo-three-dimensional hybrid bismuth polyiodides: Synthesis, structure, and optical properties of phenylenediammonium polyiodobismuthate(III)

_{3}NH

_{3}PbI

_{3}

_{3}NH

_{3}PbI

_{3}

^{2+}-doped ZnO quantum dots

_{3}NH

_{3}PbI

_{3}