Ultrafast charge carrier dynamics of methylammonium lead iodide from first principles

Methylammonium lead iodide (MAPbI3) has been a major focus of photovoltaic research for the last decade. The unique interplay between the structural and electronic properties of this material contributes to its exciting optical properties especially under the action of an ultrafast laser pulse. First-principles methods like real-time time-dependent density functional theory (RT-TDDFT) enable performing corresponding simulations without the aid of empirical parameters: the gained knowledge can be applied to future studies of other complex materials. In this work, we investigate the ultrafast charge-carrier dynamics and the nonlinear optical response of MAPbI3 excited by a resonant pulse above the gap. First, we examine the electronic and optical properties in the static regime. Next, we impinge the system with a femtosecond field of varying intensity and follow the evolution of the photoexcited carrier density. A pronounced intensity-dependent response is observed, manifested by high-harmonic generation and nonlinear trends in the number of excited electrons and excitation energy. Our results provide relevant indications about the behavior of MAPbI3 under strong and coherent radiation and confirm that RT-TDDFT is a viable tool to simulate the photo-induced dynamics of complex materials from first principles.


Introduction
The transfer of energy from light sources to materials is a fundamental physical phenomenon of great technological relevance [1].While the interaction with weak and incoherent radiation such as solar light triggers the linear response, which can be probed by conventional techniques such as UV-vis absorption spectroscopy, intense sources give rise to more complex processes such as high-harmonic generation (HHG) [2].This effect [3,4] has attracted a lot of attention when detected in established semiconductors such as ZnO [4,5], SiO 2 [6,7], and GaSe [8].From a theoretical viewpoint, HHG can be described by semi-classical models [9] or from first principles in the framework of high-order perturbation theory [10].The latter approach, while very effective and predictive, is computationally very demanding and can hardly be applied to complex materials.A viable alternative is offered by real-time time-dependent density functional theory (RT-TDDFT) which provides a non-perturbative scheme to access the response of materials at all orders [11,12].In this formalism, it is possible to follow the dynamics of the electronic distribution under a time-dependent electromagnetic field of tunable intensity, shape, duration, and frequency, that evolves together with the excited system over a certain time window.This method has been successfully applied to disclose the response to coherent radiation of bulk crystals [13][14][15][16][17], molecules [18][19][20][21][22], as well as organic, inorganic, and hybrid interfaces [23][24][25][26][27][28].A recent development based on RT-TDDFT enables the fully ab initio simulation of two-dimensional electronic spectroscopy [29].
In the last decade, new classes of semiconductorls have emerged [30,31].Among them, hybrid halide perovskites have drawn particular attention thanks to their astonishing photovoltaic performance [32][33][34].These characteristics have stimulated extensive research on their interaction with visible light [35][36][37], with a particular focus on the regime of weak radiation that is relevant to photovoltaic applications.Further studies on the response of halide perovskites to intense light sources have revealed a plethora of intriguing phenomena occurring in these systems, ranging from phonon-driven Rabi splitting [38] to pulse generation and lasing [39,40].Yet, the complex structural and chemical properties of halide perovskites have so far inhibited systematic first-principle analysis exploring the response of these systems to intense, coherent radiation on a fundamental level.This knowledge is relevant to reveal which characteristics of the incoming electric fields give rise to optical nonlinearities.At the same time, applying RT-TDDFT to these materials represents an important test to assess its performance and limits, in view of future developments and applications to even more complex material classes.
We investigate the response of methylammonium lead iodide (MAPbI 3 ) to femtosecond (fs) pulses of varying intensities adopting the framework of RT-TDDFT.We start by analyzing the electronic structure of MAPbI 3 contrasting the results obtained in the local-density approximation (LDA) [41] with those provided by the more advanced r 2 SCAN functional [42].Supported by the qualitative agreement of the two methods, we proceed to the analysis of the optical properties using adiabatic LDA, first in the linear regime to identify the spectrum of excitations and subsequently by perturbing the system with a resonant pulse of increasing intensity, ranging from 6 GW cm −2 , for which the response of the material is still linear, up to 5 TW cm −2 driving nonlinear effects.We monitor the response of MAPbI 3 by inspecting the time-resolved current density and its Fourier transform corresponding to the HHG spectrum, and by considering the behavior in time of the excitation energy and the number of excited electrons during irradiation.The advantages and disadvantages of RT-TDDFT applied to MAPbI 3 are finally discussed.

Computational methods
The calculations performed in this work are based on density functional theory (DFT) [41,43], providing the ground-state electronic properties of the material, and on its time-dependent extension (TDDFT) [44], to follow in real-time the response of the system to ultrashort time-dependent electric fields of varying intensities.The LDA is adopted for the exchange-correlation potential.However, given its known limitations in describing band gaps, the ground-state electronic structure of MAPbI 3 is additionally calculated using r 2 SCAN [42].
In the framework of RT-TDDFT for periodic systems, the time-dependent Kohn-Sham (KS) equations are solved in the adiabatic LDA assuming the velocity-gauge formulation [45] with the explicit inclusion of an external ultrafast electric field, where A(t) is the vector potential driving the electronic density out of equilibrium.The response of the system is expressed by its time-dependent macroscopic current density, where Ω is the primitive cell volume and is the microscopic time-dependent current density, where ϕ nk are the time-dependent KS states.The expression reported in equation ( 3) is general and holds regardless of the specific TDDFT implementation.
Here, we adopt the code OCTOPUS [46] (see computational details below) which employs pseudopotentials to treat core electrons [47].In this case, the current density J(t) includes an additional contribution accounting for the nonlocal part of the pseudopotential, as extensively discussed in [47] (see equation (5) therein).Following the approach adopted in [27,48], the HHG spectrum, H(ω), is computed by taking the squared norm of the Fourier transform of the time derivative of the unit-cell averaged (i.e.macroscopic) current density (equation ( 2)): The excitation energy per atom and the number of excited electrons per atom at a certain time t are calculated as and respectively, where ε bk are the KS energy eigenvalues, N el and N at the total number of electrons and atoms in the system, respectively, and P b,k (t) = occ b ′ |⟨ϕ bk (0)|ϕ b ′ k (t)⟩| 2 with 0 ⩽ P b,k (t) ⩽ 2 is the (time-dependent) occupation of state ϕ bk (t) summed over spin .
All calculations are performed with the code OCTOPUS [46] adopting a spacing of 0.18 Å for the real-space grid and a 6×4×6 k-mesh to sample the Brillouin zone.We checked the robustness of this k-sampling in the nonlinear regime, when the material is subject to the strongest laser intensity.In all calculations, we adopt Hartwigsen-Goedecker-Hutter [49] pseudopotentials implemented in the velocity gauge accounting for non-local terms and preserving gauge invariance [45].The treatment of dynamical properties can be enhanced by incorporating the effects of core electrons [50,51].The band-structure calculations with the r 2 SCAN+rVV10 functional [52] are carried out using the plane-wave code VASP [53] on a 4 × 3 × 4 k-mesh, adopting an energy cutoff of 37 Ry and ultrasoft pseudopotentials.Spin-orbit coupling is included in the band-structure calculations.To compute the linear absorption spectrum, the propagation of the KS states is triggered by an instantaneous kick [54] with amplitude κ = 10 a.u.applied along each Cartesian direction.The time step is set to 4.84 as for a total duration of 20 fs reached with the approximated enforced time-reversal symmetry propagator [55].The time-dependent electric field is modeled by a Gaussian-enveloped pulse centered at t = 10 fs with a standard deviation σ = 4 fs and a carrier frequency of 3.12 eV.Five different intensities are considered (expressed in decreasing order): 5 × 10 12 W cm −2 , 1 × 10 12 W cm −2 , 1 × 10 11 W cm −2 , 1 × 10 10 W cm −2 , and 6 × 10 9 W cm −2 .

Electronic structure
In this study, we consider bulk MAPbI 3 (chemical formula: CH 3 NH 3 PbI 3 ) in its orthorhombic phase (space group Pnma).This material is characterized by an inorganic lead-iodide scaffold intercalated by methylammonium cations (figure 1(a)) and modeled using the experimental lattice parameters obtained from single-crystal x-ray diffraction at 100 K [56]: a = 8.84 Å, b = 12.58 Å, and c = 8.55 Å.
In figures 1(b) and (c), we compare the band structures of MAPbI 3 obtained with the LDA and r 2 SCAN+rVV10 functionals, respectively.Qualitatively, the two results are very similar, both featuring a direct band gap at Γ.With LDA, it is equal to 1.47 eV in agreement with earlier predictions at the same level of theory [57,58] and also, remarkably, with experimental findings [59][60][61], due to error cancellation between spin-orbit coupling and many-body interactions [62][63][64].On the other hand, as expected, the r 2 SCAN+rVV10 functional yields a larger band gap of 2.1 eV, which is much closer to the result obtained from GW calculations [59].In the following, we compute the optical properties of MAPbI 3 on top of the electronic structure from LDA.The qualitative agreement with the results obtained with r 2 SCAN+rVV10 confirms the reliability of this approximation for the subsequent RT-TDDFT simulations of pulse-driven electronic dynamics.

Optical properties
We start the analysis of the optical properties of MAPbI 3 by examining its linear absorption spectrum.Due to the anisotropy of the crystal, we obtain three different signals by setting the laser polarization along each Cartesian axis (see figure 2).The absorption onset is around 1.5 eV, consistent with the LDA band gap of 1.47 eV.The error cancellation mentioned above in the analysis of the band structure manifests itself also in the optical spectrum, which overall matches optical measurements [59][60][61]65].This agreement is explained by a compensation between the quasi-particle correction enhancing the gap and the exciton binding energy as well as the spin-orbit coupling reducing it [62,[65][66][67].From the onset, characterized by three weak excitations polarized along each Cartesian direction, the absorption spectrum grows steeply and exhibits intense maxima around 2.0 eV, 2.5 eV, and 3.4 eV.All these features are in agreement with experimental results [68][69][70].At higher energies, the absorption becomes more isotropic with similar contributions from the three polarizations (figure 2).
With this knowledge, we examine the optical response of MAPbI 3 to an external time-dependent field.To this end, we excite the system with an ultrafast pulse of 20 fs duration, with bandwidth of 0.21 eV (see figure 2) and polarized along the x direction according to the reference system in figure 1(a), with carrier frequency ω 0 =3.12 eV, and varying intensity sweeping a range from 6 GW cm −2 up to 5 TW cm −2 .We did not explore other pulse widths or shapes as in our experience, these parameters affect only the details but not  the main outcomes of the calculations.The applied fields are visualized in black in figure 3, overlaid by the induced current densities J(t) (equation ( 2)) in color.Under a weak field intensity (6 GW cm −2 , figure 3(a)), the response of the system is linear: the current oscillates with the laser frequency and it is in quadrature with respect to the pulse oscillations (including a π/2 phase shift).As soon as the pulse is switched off, the amplitude of the current density drops, and the oscillations due to the residual coherence between groundand excited-state in the absence of dissipation channels persist until the end of the simulation.When the intensity increases to 10 GW cm −2 , the current density is enhanced accordingly during irradiation; after the pulse, weak oscillations remain, exhibiting a hint of amplitude modulation (figure 3(b)).This effect becomes more prominent when the laser intensity is further increased by one order of magnitude (100 GW cm −2 , see figure 3(c)): while the amplitude of the post-pulse residual oscillations is enlarged compared to the results obtained under weaker lasers, the response of the system remains linear.
The scenario changes dramatically when the field has a peak intensity of 1 TW cm −2 (see figure 3(d)).In this case, the amplitude of J(t) changes abruptly after approximately 10 fs, when the external field reaches its maximum.In the remaining 10 fs of laser activity, the current density keeps oscillating with about half the amplitude reached in the initial phase of irradiation.As soon as the pulse is turned off, J(t) decreases and, as in the previous cases, weaker oscillations coming from the residual electronic coherence induced by the laser manifest themselves.The abrupt change of the current around 10 fs can be ascribed to two effects: (i) the activation of third-order optical non-linear response (the second-order response is zero in a centrosymmetric material like MAPbI 3 ) that interferes with the linear one, as seen in the HHG spectrum in figure 4; (ii) the current density has an explicit dependence on the field and evolves proportionally and in-phase with the laser amplitude during irradiation.We can assume that upon large field strengths, there is interference with the linear-response signal that attenuates the total current density.
By further increasing the intensity up to 5 TW cm −2 (figure 3(e)), the non-linear response of the material becomes stronger and the amplitude of the current density is no longer regularly modulated by the envelope of the field.Moreover, the post-pulse signal exhibits an irregular evolution which is an indicator of highly nonlinear response [47,71].It is worth noting that the current densities displayed in figure 3 have compatible values with analogous results obtained for diamond on the same theoretical footing [16,72].In that material, pulse intensities of the order of 100 TW cm −2 with durations of 16 fs are predicted to induce optical breakdown [73].In a softer compound like MAPbI 3 , this threshold is reasonably reduced by about two orders of magnitude, as hinted by our findings.
To better characterize the (nonlinear) response of MAPbI 3 under laser irradiation, we compute the HHG spectra H(ω), see equation ( 4), at different laser intensities.The results shown in figures 4(a)-(c) indicate that for pulses with an intensity up to 100 GW cm −2 , only the first (fundamental) harmonic appears at ω 0 =3.12 eV, corresponding to the pulse frequency.This behavior confirms the qualitative observation reported above that the response of MAPbI 3 stays linear under pulses with peak intensities below I = 100 GW cm −2 and is in line with other works [71,74].As the intensity reaches 1 TW cm −2 , three peaks are visible in H(ω), see figure 4(d): the most prominent one at ω 0 =3.12 eV is accompanied by the third and fifth higher harmonics at ω 3 = 3ω 0 = 9.21 eV and ω 5 = 5ω 0 = 15.70 eV, respectively.These peaks are approximately one-half and one-fourth the intensity of the fundamental one.Finally, under the strongest pulse with I =5 TW cm −2 , the HHG spectrum of MAPbI 3 changes dramatically (see figure 4(e)).In this case, the third and fifth harmonics dominate over the fundamental one which is hardly visible despite the magnified scale.No further signals appear in our calculations above 20 eV.The threshold intensity in the TW cm −2 range for the appearance of the third and higher harmonics, as well as their relative ratio as a function of the laser intensity is in line with previous calculations [74].
We complete the characterization of the nonlinear optical response of MAPbI 3 by examining the time evolution of the number of excited electrons per atom (equation ( 6)) and of the excitation energy, calculated as the difference between the total energy at time t and the ground state energy divided by the total number of atoms [75]: ∆E ex (t) = (E(t) − E(t = 0)) /N at .When positive, ∆E ex (t) indicates an energy uptake of the system from the laser.This quantity is monitored over the entire 40 fs window considered in these simulations.The energy uptake is negligible for intensities below 100 GW cm −2 : the corresponding curves are hardly visible in figure 5(a).For I =100 GW cm −2 , ∆E ex (t) starts to grow between 10 and 15 fs when it reaches a maximum and forms a plateau for the remaining propagation time.This behavior, which is due to the absence of dissipation channels in the simulation, was reported for insulators [73], semiconductors [2], graphite [47], and dielectrics [76], investigated at the same level of theory.When the intensity increases by one order of magnitude (I =1 TW cm −2 ), the response of MAPbI 3 becomes nonlinear and ∆E ex rises shortly after 5 fs, reaches a maximum of about 0.25 eV atom −1 between 10 and 15 fs, and then slightly decreases until 20 fs, maintaining a constant value around 0.22 eV atom −1 until the end of the simulation (figure 5(a)).The result obtained for I =5 TW cm −2 indicates a huge non-linear increase in the energy uptake by the system from the external field.In this case, ∆E ex starts increasing immediately after the pulse is switched on and grows steadily, although some oscillations are present until t = 15 fs when it reaches its maximum value of 1.2 eV atom −1 .Subsequently, the excitation energy decreases slightly, forming a plateau at 1.15 eV atom −1 that persists until the end of the simulation.The energy absorbed per unit time and volume by the system is given by U [47] and the excitation energy per atom in the unit cell can be calculated as From these equations, it is evident that the amount of absorbed energy at each time step crucially depends on the amplitude and the relative phase between the induced current and the external field.When the current and the field are in phase (anti-phase) the system dissipates (absorbs) energy from the laser.For this reason, the overshoots are associated with the time frame in which the system absorbs the maximum amount of energy from the external field and the latter reaches its minimum energy.The maximum happens when the product between the induced current and electric field is minimum, i.e. when they are in anti-phase.
For a clearer visualization of the behavior of ∆E ex as a function of the pulse intensity, we plot this quantity at the end of the propagation time, t = 40 fs, see figure 5(b).The obtained trend confirms our findings discussed above: for intensities up to 10 GW cm −2 , the response of the system is linear and there is almost no energy uptake from the laser.From 100 GW cm −2 and, in particular, above 1 TW cm −2 , the nonlinear response of the material dominates, as indicated by the steep increase of the energy uptake.It is worth mentioning that the value of 1 TW cm −2 identified here as the threshold fornonlinear behavior in MAPbI 3 under 20 fs laser irradiation is compatible with the laser used in Z-scan experiments to access multiphoton absorption coefficients in this material [77][78][79].
In the final step of this analysis, we consider the time evolution of the number of excited electrons per atom, N ex .The trend displayed in figure 5(c) is qualitatively similar to the one for ∆E ex shown in figure 5(a): this is expected as the two quantities are intrinsically connected.For laser intensities below 100 GW cm −2 , the amount of excited electrons is negligible and the corresponding curves are hardly visible in the plot.With I = 100 GW cm −2 , only 0.03 electrons/atom are excited after 40 fs (see figure 5(d)): N ex (t) rises between 7 fs and 15 fs when it reaches its maximum and remains constant until the end of the simulation.Intensities equal to 1 TW cm −2 or exceeding this value trigger optical nonlinearities in the material and, consequently, a nonlinear increase in the number of excited electrons, mirroring the behavior of the excitation energy.N ex (t) grows steeply for the first 12 fs of the simulation when it reaches its maximum of 0.08 electrons/atom (see figure 5(c); afterward, it decreases until t = 20 fs when the plateau starts and lasts til the end of the simulation.When the pulse intensity reaches 5 TW cm −2 , the number of excited electrons grows up to 0.20 electrons/atom at t = 12 fs; subsequently, N ex (t) decreases and reaches the constant value of 0.15 electrons/atom (see figure 5(c)).In both ∆E ex (t) and N ex (t) (figures 5(a) and (c)), we notice the appearance of intensity-dependent fast oscillations on top of the slower build-up dynamics during the transient phase.A careful inspection of these features reveals that, while their amplitude depends on the field strength, their oscillations have always twice the frequency of the resonant laser.By using time-dependent perturbation theory, one can show that these rapid oscillations originate from the anti-resonant terms of the transition probabilities or populations P nk (t) in equations ( 5) and ( 6) during the initial transient phase and they disappear after the laser is turned off.At first order in the field strength, these oscillations of the transition probability are proportional to the laser intensity and inversely proportional to the square of the resonance frequency [10,80].
The number of excited electrons per atom as a function of the pulse intensity is summarized in figure 5(d).Qualitatively, the trend is similar to the one of the energy uptake (compare figure 5(b)), although the latter quantity grows more steeply for the strongest considered pulses.At t = 40 fs, about 0.05 electrons/atom are excited by the pulse with I = 1 TW cm −2 , while for I =5 TW cm −2 , N ex increases to 0.15 electrons/atom (see figure 5(d).We checked that the obtained results are qualitatively independent of laser polarization by performing an additional calculation with the pulse polarized in the (1,1,1) direction.

Summary and conclusions
In summary, we investigated the (nonlinear) optical properties of MAPbI 3 under (intense) laser irradiation using RT-TDDFT.As a preliminary step, we characterized the electronic structure of the material benchmarking the results obtained with LDA with those provided by the more advanced r 2 SCAN+rVV10 functional: the good qualitative agreement between them supports the use of adiabatic LDA in the subsequent time-dependent simulations.The linear absorption spectrum, consistent with experimental results from the literature [81], is anisotropic due to the orthorhombic symmetry of the considered MAPbI 3 crystal.The laser-driven charge-carrier dynamics is computed on an ultrashort time-window of 40 fs by applying a pulse in resonance with a strong absorption maximum at 3.12 eV, with an intensity ranging from 6 GW cm −2 up to 5 TW cm −2 .The response of the system is monitored through the time evolution of the macroscopic current density and its Fourier transform, yielding the HHG spectrum.From both quantities, the nonlinear character of the response is evident for pulse intensities beyond 1 TW cm −2 .In particular, under the strongest considered laser with I = 5 TW cm −2 , the fundamental peak at ω 0 = 3.12 eV is hardly visible in the HHG spectrum while the overtone signals dominate.These trends are confirmed by the energy uptake and the number of excited electrons monitored over the entire simulation window.
In conclusion, our study indicates I = 1 TW cm −2 as the threshold for the nonlinear optical response of MAPbI 3 to ultrafast pulses in resonance with its above-edge absorption maximum at 3.12 eV.The sole inclusion of the electronic degrees of freedom in the simulations inhibits a direct comparison between our findings and experimental data, given the prominent role of vibronic couplings in the probed halide perovskite samples [38].The results reported in this work provide a useful point of comparison to assess the contribution of the electronic dynamics in (intense) irradiation processes of halide perovskites.From a methodological perspective, this work represents an important case study for RT-TDDFT: under (intense) excitation well above the absorption onset, where excitonic effects are negligible, this method has proven able to describe the electronic dynamics of MAPbI 3 without the aid of empirical parameters and at affordable computational costs.
In future work, the role of the nuclear degrees of freedom needs to be addressed in detail by means of RT-TDDFT simulations in conjunction with Ehrenfest dynamics, possibly in a multi-trajectory flavor [82][83][84] to gain insight into the complex behavior of nuclear motion [85][86][87] and to ultimately achieve an all-around understanding of the nonlinear response of MAPbI 3 to ultrafast radiation.Further research directions include the analysis of scattering mechanisms induced by HHG and other optical nonlinearities in this material.To this end, more complex modeling of the sample including impurities [88], grain boundaries [89], or surface roughness [74] is needed.This analysis is expected to provide a better understanding of the nonlinear optical response of MAPbI 3 , carrying important implications for its photovoltaic efficiency and optoelectronic performance.

Figure 1 .
Figure 1.(a) Ball-and-stick representation of the unit cell of MAPbI3 in its orthorhombic phase (space group Pnma).Band structure of MAPbI3 computed with (b) LDA and (c) r 2 SCAN+rVV10: the Fermi energy is set to zero at the valence-band maximum.

Figure 2 .
Figure 2. Linear absorption spectrum of MAPbI3 computed as the imaginary part of the macroscopic dielectric function ϵ along different polarization directions x, y and z, parallel to the crystal lattice vectors a, b and c, respectively.The red, shaded area indicates the spectral region excited by the pump polarized along x.The arrow marks the carrier frequency of the laser, ω0 = 3.12 eV.

Figure 5 .
Figure 5. (a) Excitation energy per atom as a function of time, ∆Eex(t), computed at varying laser intensities; (b) ∆Eex(t = 40 fs) at increasing pulse intensity.(c) Number of excited electrons per atom, Nex(t), as a function of time, computed at varying laser intensities; (d) Nex(t = 40 fs) at increasing pulse intensity.