Ionic polarization and dielectric function play a fundamental role in the optoelectronic properties of hybrid perovskites, currently one of the most studied materials for next generation photovoltaics. The hybrid nature of the crystal, with molecular dipoles that can reorient within the inorganic lattice, gives rise to a complex dielectric response in the bulk material that has been largely studied and debated. Here, we investigate the nature and the relaxation properties of the dielectric polarization of hybrid perovskites at finite temperature by means of classical molecular dynamics. We provide evidence that a simple ionic model of classical interatomic forces is able to explain qualitatively the temperature and frequency dependence of the dielectric constant providing a picture that is fully consistent with experimental data. The constant dielectric function in the low-temperature phase is controlled by ionic displacements, while the temperature-dependent paraelectric behavior of the tetragonal phase is due to reorientation of dipoles that are responsible for the discontinuity at the orthorhombic-to-tetragonal transition. In the frequency domain, the molecular reorientations give rise to a broad band that is located in the 0.1 THz timescale at room temperature and that shifts down to the GHz timescale when cooling the system toward the tetragonal-to-orthorhombic phase transition. The relation between relaxation time and maximum absorption frequency is also clarified.

## I. INTRODUCTION

Hybrid perovskites of the form MAPbX_{3} with *X* = Cl, Br, I anions and MA = $CH3NH3+$ (methylammonium) cation have attracted great interest in recent years for their excellent photovoltaic properties^{1} and for potential applications also in optoelectronics,^{2} ionizing radiation imaging,^{3} nanoelectronics,^{4} and thermoelectricity.^{5} MAPbX_{3} is characterized by a high dielectric constant at room temperature (*ε*′ ∼ 60 for iodides and *ε*′ ∼ 50 for bromides)^{6} and by a sizable dependence on temperature and frequency.^{6–8} The bulk polarization and the dielectric constant have important implications in the screening of photogenerated carriers,^{9,10} the formation of polarons,^{11} and the exciton binding energy^{12,13} and represent an important issue for a deep fundamental understanding of the material. Although the magnitude of the dielectric constant is well-known from experiments, the nature of the polarization and the role of molecules are not easily accessible and have been largely discussed in the literature also in terms of hypothetical ferroelectric ordering of molecules.^{8}

Early studies on the dielectric constant of hybrid perovskites date back to the 1990s. Based on millimeter-wave spectroscopy,^{7} the dielectric constant was measured at frequencies of 50–150 GHz and temperatures of 100–300 K, indicating a discontinuity in the dielectric constant at the orthorhombic–tetragonal phase transition for the whole family of hybrid perovskites MAPbX_{3} (X = I, Br, Cl). The authors also provided an interpretation of data in terms of the Debye relaxation model according to which the dependence of the dielectric constant (*ε* = *ε*′ + *iε*″) on temperature and frequency is

where *T* is the temperature, *ω* the frequency of the applied electric field to measure the dielectric constant, and *τ*(*T*) is the relaxation time of the dipoles. They also proposed an Arrhenius *T*-dependence of the relaxation time $1\tau \u223ce\u22121T$ that was attributed to the fast reorientation of molecular dipoles (within the picosecond timescale) along the eight high symmetry orientation of the cubic (tetragonal) crystal. For measurements below kHz frequency, it is *ωτ* ≪ 1 and $\epsilon 0\u2032\u223cCT$ is nearly independent of the frequency, but decreases with *T*.

In 1992, Onoda-Yamamuro *et al.*^{6} extended measurements between 20 Hz and 1 MHz and at low temperatures confirming the sizable *T*-dependence, with sharp bend of the real part of the dielectric permittivity occurring at the orthorhombic–tetragonal phase transitions. The real part of the static dielectric constant was measured to be ∼36 in the orthorhombic phase (*T* < 162 K).

Recently, Govinda *et al.*^{14} focusing on MAPbBr_{3} compared the strong temperature dependence of the near room temperature dielectric constant measured in the 1 kHz−1 MHz range with the nearly temperature-independent dielectric constant of CsPbBr_{3}. This comparison of the MA system with the corresponding Cs system indicated a strong temperature dependence only for the hybrid case, confirming the picture of dipoles rotating faster with respect to the time scale of measurement. The authors also concluded that for measurements up to kHz, there is very little frequency dependence. The occurrence of slight deviations at high-temperature was interpreted by the Maxwell–Wagner expression as the relaxation of grain boundaries and other extrinsic contributions. In contrast to the above reports, Anusca *et al.*^{8} recently reported a sizable frequency dependence in the GHz range associated with relaxation times in the timescale of 100–200 ps. However, this dependence was found in a range of temperatures (150–201 K) around the orthogonal-to-tetragonal transition. Overall, most experimental studies indicate that molecules contribute to the dielectric constant of hybrid perovskites and to its temperature dependence at the phase transitions.

The dielectric function of hybrid perovskites has been investigated theoretically using *ab initio* methods based on density functional theory (DFT).^{15} By using the random-phase-approximation,^{16} the electronic contribution to the static dielectric constant of tetragonal MAPI was calculated to be $\epsilon 0\u2032\u2009\u2243\u20096$. This value is small with respect to the measured static dielectric tensor indicating large ionic contributions, as expected for polar and relatively soft crystals as hybrid perovskites. The ionic contribution can be taken into account by density functional perturbation theory (DFPT), predicting static dielectric components in the range 18–37^{17} for cubic MAPI and ∼20 in MAPBr.^{11} A review of theoretical *ab initio* results can be found in Ref. 15. *Ab initio* calculations based on modern theory of polarization provide the best accuracy for the study of specific atomic configurations of hybrid perovskites such as perfect crystals (although successful examples have been reported also for the calculation of the dielectric function of nanostructures^{18}). However, temperature related phenomena such as cation rotations and configurational disorder require large configurational averages that are typically out-of-reach of *ab initio* methods.

In this work, we apply classical molecular dynamics to study the bulk polarization of hybrid perovskites at finite temperature, focusing on the MAPbI_{3} system. We show that a simple ionic model of interatomic forces (MYP) developed by Mattoni *et al.*^{19} and extensively applied to study hybrid perovskites^{20,21} is able to reproduce qualitatively the temperature and frequency dependence of the dielectric constant, making it possible to separate organic and inorganic contributions. In particular, we provide clear evidence that the dielectric discontinuity in temperature at the orthorhombic-to-tetragonal phase is mainly due to the onset of molecular reorientations with partial ordering under electric field; concerning polarization kinetics, we show that molecular contributions give rise to a broad absorption band in the 0.1 THz range at room temperature that shifts down to 10 GHz as *T* approaches the orthorhombic phase, and we discuss the relation between the absorption frequencies and rotational relaxation times. Overall, MD provides a microscopic picture that is consistent with the large body of experimental data supporting the scenario of a paraelectric material with a sizable polarization by orientation at room temperature and anti-ferroelectric ordering at low temperature.^{8,14}

## II. METHODOLOGY

Molecular dynamics simulations are performed by using the LAMMPS code.^{22} Constant-temperature constant-pressure simulations are performed within the NPT ensemble with a triclinic fully flexible simulation cell. A time step as small as 0.5 fs has been used to properly integrate equations of motions. Relaxation times of 100 and 500 time steps have been used for the thermostat and barostat, respectively. The crystals were annealed for at least 0.5 ns at each temperature in the absence of electric fields. In the presence of periodic electric fields, the simulation times were chosen to be several times larger than the period of the electric field and long enough to get converged results at all temperatures, shapes, and magnitudes of electric field. Typical simulations times are in the range 200 ps–2 ns. The MYP0 force-field^{23} used in this work has been extensively applied to study several properties of hybrid perovskites. The model is able to reproduce satisfactorily the structural, elastic, thermal,^{21} and vibrational properties^{24} of the different crystal structures of hybrid perovskites. In particular, it reproduces the Pnma orthorhombic crystal, the correct tilting patterns of PbI octahedra, and the orthorhombic-to-tetragonal transition at 160 K.^{23} The model consists of a combination of Buckingham–Coulomb interactions for inorganic and GAFF^{25} for molecules. A comprehensive description of the model can be found in Refs. 19 and 23. The simulation using LAMMPS requires a hybrid style (lj/charmm/coul/long and buck/coul/long) and a cutoff as large as 10 Å. The long-range electrostatic interactions in triclinic cells require the Ewald method. A relative accuracy of 10^{−6} has been used. The typical systems^{23} used to perform molecular dynamics of orthorhombic or tetragonal crystals contain 256 formula units. Smaller systems containing 32 units have been used in this work to decrease the computational time. Since the orthorhombic-to-tetragonal phase transitions for such small crystals are size dependent, it is necessary to rescale the simulated temperatures. For the case of 32 units, we find that a factor of 0.47 makes it possible to match the experimental data. By this strategy, it was possible to largely reduce the computational cost and better control convergence of slow components of polarization. Boundary conditions^{26} can also, in principle, affect the calculated electrostatic energies, the polarization, and, in turn, the actual values of the dielectric functions. By comparing orthorhombic crystals1b of different sizes and shapes (from 32 to 500 units), we have found differences within 5%. Similar errors in the actual value of the dielectric constant are not critical for the present study since the effect of temperature and frequency is studied comparatively on systems with the same size. Time-varying electric fields $E(t)$ were applied during the dynamics by adding the instantaneous force $qE$ to internal forces. Depending on the type of calculations (see below), we applied periodic (harmonic, square-, and linear-waves) and non-periodic (e.g., Heaviside single step) electric fields by defining suitable instantaneous variables and averages within the LAMMPS input. The force field makes use of constant atomic charges that is an acceptable choice for crystals but also to study degradation phenomena.^{27} Under the effect of the electric field, the crystal center of mass is invariant (due to crystal neutrality), while positive and negative charges translate or rotate in opposite directions, hence inducing polarization. The (variation) of *i*th Cartesian component of polarization *P*_{i} is obtained during the dynamics by time-integrating the atomic currents according to the formula

where *i* is the Cartesian coordinate, *a* labels the atomic index, *q*^{a} is the (constant) atomic charge, and $via$ is the velocity in the same Cartesian component. Such a definition makes it possible to separate the contributions from different groups of atoms, and it also avoids the calculation of dipoles that are ill-defined in periodic systems. The electric susceptibility *χ* is obtained by comparing the induced polarization *P* and the macroscopic applied field $E$, according to relation $\chi =P/(\epsilon 0E)$. *ε*_{0} is the vacuum permittivity; the permittivity and relative permittivity are obtained from the calculated susceptibility *ε* = *ε*_{0}(*χ* + 1) and *ε*_{r} = *χ* + 1, respectively. *ε* and *ε*_{r} depend on temperature and frequency and are referred to as dielectric functions.

## III. TEMPERATURE DEPENDENCE OF THE BULK POLARIZATION

The equilibrium crystal structure of MAPbI_{3} at low temperature is orthorhombic. At 160 K, the crystal undergoes an orthorhombic-to-tetragonal phase transition. Such a thermodynamic evolution of the crystal is well reproduced by the present interatomic force-field.^{23} Molecular dynamics shows that the ionic polarization *P*(*t*) of the MAPbI_{3} crystal depends on temperature *T*, with a qualitative different behavior in the orthorhombic and tetragonal phases. In Fig. 1, panel (b) is reported as an example of the total polarization (black line) for the orthorhombic crystal induced by a periodic square-wave electric $E(t)$ in the *x* direction [panel (a)]. The amplitude and period are $E0=0.05$ MV/cm and 20 ps, respectively. The organic (MA) and inorganic (PbI_{3}) contributions are separated and shown as orange and blue lines, respectively. All the polarization components (thin lines) fluctuate rapidly. The averages (thick lines) follow the square profile of the applied electric field with very short transients (∼1 ps) at each inversion. The polarization is due to opposite displacements of the positive (Pb and MA ions) and negative (iodines) sublattices forming the crystal [panel (d) of Fig. 1, arrows] induced by the electric field. The major contribution to the polarization of the orthorhombic crystal is due to the polar distortions of the inorganic lattice [see red and black arrows in panel (d)].

The behavior is qualitatively different in the tetragonal crystal. The PbI_{2} lattice expands due to thermal fluctuations, and the molecules are able to reorient^{28} [see panel (e) of Fig. 1]. The Orientational disorder is well reproduced by the present force field^{19} with results that are in agreement with experiments based on electromagnetic^{6,7} and neutron spectroscopies.^{28} In the tetragonal crystal (*T* > 160 K), the polarization increases in magnitude [see Fig. 1, panel (c)] due to an increase in both organic and inorganic contributions. Remarkably, the ordering of organic/inorganic contributions is inverted when passing from orthorhombic to the tetragonal phase with the molecular polarization that becomes dominant. A second important qualitative difference is the occurrence of longer polarization transients, i.e., slower relaxation kinetics. The transients can be fitted by a periodic signal with single exponential decay,

where $g(t)=(\u22121)tt0$ inverts the sign at each half-period *t*_{0} and $s(t)=t\u2212t0tt0$ is the time relative to each half-period.

The proposed time evolution of $P(t)$ corresponds to the linear differential equation of a dissipative system with a single relaxation time $\tau $ under the effect of an external electric field $E(t)$,

where, for simplicity, $P(t)$ and $E(t)$ are normalized to the corresponding initial values. The response to an electric impulse at $t=t0$, i.e. $E(t)=\delta (t\u2212t0)$ is the susceptibility (i.e. Green function) that is $P(t)=\theta (t\u2212t0)e\u2212(t\u2212t0)/\tau $, where $\theta (t)$ is the Heavside step-function, i.e. $\theta (t)=0$ for at $t<0$ and $\theta (t)=1$ for $t>0$. The response to the a finite step $E(t\u2032)=[\theta (t\u2032\u2212t)\u2212\theta (t\u2032)]$ is $P(t)=\u222bdt\u2032G(t\u2212t\u2032)1\tau E(t\u2032)$ i.e. $P\u223c(1\u2212e\u2212t/\tau )$. This corresponds to the model of Eq. (2). The model is able to fit the atomistic data at all temperatures, providing information on relaxation time and polarization. If the step is long enough, the value $P0(E0)$ approaches the static polarization of the system $P0=\chi \u2002\epsilon 0\u2002E$ and gives the static electric susceptibility and dielectric constant *ε*_{r} = *χ* + 1.

By fitting *χ* on atomistic data, the dielectric constant *ε*(*T*) is obtained as a function of *T*, which is reported in Fig. 2. By decomposing *P* into organic and inorganic components $P=PMA+PPbI3$, the inorganic (blue) $PPbI3(\epsilon 0E0)\u22121$ and molecular (orange) $PMA(\epsilon 0E0)\u22121$ components of the static susceptibility can be obtained. The experimental data are reported in Fig. 2 showing that MD reproduces qualitatively the temperature dependence of the dielectric function. For a better comparison, the experimental and atomistic data are normalized to the corresponding value of the dielectric constant in the orthorhombic phase. It is found that the dielectric coefficient in the orthorhombic phase is almost *T* independent. At the orthorhombic-to-tetragonal transition, an abrupt increase in *ε* is observed, giving rise to a maximum at *T* ∼ 160 K where the dielectric coefficient is three times larger than that of the orthorhombic phase [i.e., *ε*_{r}(*T*) ∼ 3*ε*_{r}(0)]. By further increasing *T*, the dielectric constant decreases as the inverse of temperature. The MD data have a slower *T* decrease. At room temperature, *ε*_{r}(300 K)/*ε*_{r}(0 K) is ∼1.5 for experiments and ∼2.5 for MD. It is remarkable that the MD data are able to reproduce qualitatively the *T* dependence across the phase transition.

The actual values of the relative dielectric constant obtained by MD for orthorhombic (0 K) and tetragonal phases (300 K) are *ε*_{ort} ∼ 4.5 and *ε*_{tet} ∼ 12, respectively, and are quite underestimated with respect to the experimental values (∼30 and ∼60). The discrepancy is not surprising considering that the MYP0 model was not fitted on dielectric properties; and it tends to be stiffer (i.e., smaller atomic displacements and, in turn, smaller dipoles induced by external forces). Furthermore, the model does not include electronic polarization (*ε*^{elec} ∼ 6 according to DFT calculations^{16}). Quantitative improvements of the force-field are expected by including core–shell degrees of freedom or refitting parameters using the dielectric function in the database. However, the absolute value of the dielectric constant is not critical for the present study that investigates the nature of the polarization and temperature/frequency dependence.

### A. Polarization of the tetragonal crystal due to molecular dipole reorientation

Having shown that the molecules give the major contribution to the dielectric constant in the tetragonal phase (*T* > 160 K), it is interesting to characterize the nature of this polarization. To this aim, we consider the MAPbI_{3} crystal at *T* = 326 K and apply electric fields in opposite directions, as indicated by arrows in Fig. 3. The distribution of molecular dipoles over the unit sphere (see right panels of Fig. 3) shows that the largest fraction of molecules are not oriented along the electric field but distributed quasi-randomly over the unit sphere.^{19} A more quantitative analysis reveals, however, that there is a sizable polarization by orientation in the range of electric fields here investigated, i.e., 0.001–0.05 MV/cm. To quantify the polarization, we calculate the angular probability [i.e., the number of molecules $P(\theta ,\varphi )=n(\theta ,\varphi )/N$ along the direction *n*_{θ,ϕ} = (cos *ϕ* sin *θ*, sin *ϕ* sin *θ*, cos *θ*)] and we reconstruct the spherical potential *F*(*θ*, *ϕ*) through the formula *F* = −*k*_{B}*T* log *P*(*θ*, *ϕ*). The reconstructed potential is represented by a color-map (see Fig. 3, top left). In the absence of electric field, the spherical potential exhibits the tetragonal symmetry with six high probability dark regions (at the ⟨100⟩-equivalent directions of the PbI cubic-like cage) and eight red–yellow low probability regions (at ⟨111⟩-equivalent directions of the PbI cubic-like cage). The electric field breaks the inversion symmetry inducing dipolar ordering (i.e., polarization by orientation) (see the bottom panel where the blue area increases and the yellow one decreases). The opposite situation is found by reversing the electric field (compare middle and bottom panels at *E* > 0 and *E* < 0 in Fig. 3). In addition to the angular distribution, we calculate the average molecular orientation in the direction of the electric field, i.e., ⟨cos *θ*⟩, where the polar direction is that of the electric field (⟨100⟩). We verify that ⟨cos *θ*⟩ depends almost linearly on $E$ for values up to 0.05 MV/cm.

Present results on molecular susceptibility and angular distribution can be analyzed within the theory of polar liquids.^{29,30} The static susceptibility of a polar liquid can be related to the microscopic dipole moments *p*,^{30}

where *N*/*V* is the volume per dipole and *g* is a factor introduced by Kirkwood that includes the correlation between the orientations of neighboring molecules. The temperature favors the molecular disorder through the Boltzmann factor $1kT$. The average molecular orientation along the field $E$ can be calculated through the relation $\u27e8cos\u2061\theta \u27e9=PVNp$. By using the definition $P=\epsilon 0\chi E$, the average orientation in terms of microscopic dipole is obtained,

The above equations [Eqs. (4) and (5)] for susceptibility *χ* and orientation cos *θ* can be compared to the MD results for tetragonal MAPbI_{3} at room temperature, once the permanent dipole *p* of MA is calculated. For the present atomistic model, we find (after subtracting the net charge of the molecule) *p* = 0.35 *e* Å (oriented from CH_{3} to $NH3+$) that is within the range of previously reported theoretical^{15} and experimental^{6} values. Concerning the molecular component of susceptibility, *χ*_{MA}, we find that Eq. (4) for polar liquids reproduces the MD value at room temperature, *χ*_{MA} ∼ 7, provided that the correlation parameter is set to a reasonable *g* ∼ 4. This result suggests that the dielectric response of hybrid perovskites is consistent with the properties of polar liquids. We also find a qualitative agreement with Eq. (5) since the ⟨cos *θ*⟩ obtained from the analysis of the dipole orientation increases linearly with $E$. However, the actual MD average overestimates the prediction based on Eq. (5) by a factor of ∼20. The actual molecular ordering is larger with respect to the polar liquid theory but also with respect to the polarization calculated by MD. This discrepancy is likely due to the screening effect of the inorganic lattice and the competition between orientation and displacements of charged ionic dipoles. A more systematic analysis that is beyond the scope of present investigation is necessary to fully clarify the quantitative applicability of polar liquid theory to hybrid perovskites.

The present atomistic model provides a qualitative yet coherent picture of the nature of polarization and dielectric constant, which explains several experimental facts: (i) the polarization of the orthorhombic phase is almost temperature independent and controlled by internal PbI displacements with a minor contribution from molecules; (ii) the polarization of the tetragonal phase is instead dominated by molecular dipole ordering, and it has its maximum at the phase transition and it decreases with temperature, in qualitative agreement with the theory of polar liquids; (iii) the fraction of ordered dipoles overestimates the standard theory of polar liquids, likely due to the interaction with the inorganic sublattice; and (iv) molecular disorder (induced by temperature and interaction with PbI) persists in the tetragonal crystal even at very high electric fields, indicating paraelectric (and not ferroelectric) behavior.

## IV. FREQUENCY DEPENDENCE OF DIELECTRIC CONSTANT

The dielectric response of hybrid perovskites is known to depend on the frequency. Based on general arguments on the dielectric function (i.e., Kramers–Kronig relation^{31}), it is known that the dispersion (i.e., change in real part of the dielectric function) is associated with absorption bands. No absorption and constant dielectric function is expected in the very low frequency range Hz–MHz. However, measurements are difficult due to current transients and hysteresis,^{32} ferroelastic domains,^{8} and extrinsic factors such as ionic defects and grain-boundaries.^{33,34} As a consequence, measurements of the static dielectric constant are typically performed in the kHz–MHz range.

On the other hand, beyond THz, the dielectric function is controlled by polar optical phonons (30–150 cm^{−1}) and IR-active internal modes of the methylammonium molecule (900–3300 cm^{−1}). The MAPbI_{3} absorption in this frequency region has been studied by IR^{35} or Raman^{36} spectroscopy. In Ref. 24, classical MD with the same force-field adopted here has been used to calculate the IR absorption through the charge-weighted velocity-correlation function.

Here, to obtain information on the frequency dependence of the dielectric function, we calculate the Fourier transform of the polarization *P*(*t*) induced on the MAPbI_{3} crystal by a time-dependent step-like electric field $E(t)=E0\theta (t)$. The Fourier transform of the step function is $\theta (t)\u2192F\pi \delta (\omega )+1i\omega $, and, since the polarization is the convolution of *χ* and $E$ [i.e., $P(t)=\u222b\u2212\u221e\u221e\u2002dt\u2032\u2002\chi (t\u2032)E(t\u2212t\u2032)$], it follows that the complex susceptibility at *ω* > 0 can be obtained from the calculated polarization as

The complex dielectric function is obtained from *ε* = *χ* + 1 with the absolute value $(\chi \u2032+1)2+(\chi \u2033)2$. The calculated absolute value of susceptibility is reported in Fig. 4 at four different temperatures across the orthorhombic-to-tetragonal phase transition. Polar optical phonon peaks are located at 1–5 THz (30–150 cm^{−1}). Two major peaks are present at low temperature and evolve into a broad band beyond the phase transition. In addition, the internal modes of methylammonium molecules are responsible for the sharp absorption IR bands at *ν* > 30 THz (900 cm^{−1}). The data of Fig. 4 (derived from polarization) are consistent with IR theoretical and experimental analysis reported in the literature.^{24}

The absorption observed at intermediate frequencies (i.e., in the range of GHz–THz)^{8} is more subtle and requires a deeper analysis. These frequencies are below the lowest polar optical phonons and far from internal molecular motions. Anusca *et al.*^{8} reported absorption bands in the GHz range with periods as long as 100–200 ps at temperatures in the range of 150–210 K. They proposed that the bands are related to the molecular relaxation times. Here, by MD we will confirm this picture. To this aim, we start by analyzing the timescale of exponential transients of polarization of Fig. 1. The fitted relaxation times *τ*(*T*) of the polarization are reported in Fig. 5. The highest curve (orange) is associated with molecules. The T-dependence of relaxation times *τ* is similar to that of the susceptibility with an abrupt increase in the orthorhombic-to-tetragonal transition. The relaxation times *τ*(*T*) obtained here from polarization can be compared with the correlation times $\tau \u0303$ derived from angular correlation in the work of Mattoni *et al.*^{19} We find that within numerical errors, the two quantities are similar $\tau \u0303(T)\u223c\tau (T)$.

This is not accidental. The fluctuation–dissipation theorem states that the susceptibility $\chi =1\epsilon 0dPdE$ (i.e., response of a system to the application of an external electric field) is proportional to *P* fluctuations $\chi \u223c1\epsilon 012\u2002kBT\u27e8P(0)\u22c5P(t)\u27e90$ at zero field. Accordingly, the orientational susceptibility $\chi \u0303$ is related to dipole *n*(*t*) fluctuations ⟨cos *θ*(*t*)⟩ = ⟨*n*(*t*′) · *n*(*t* + *t*′)⟩. On the one hand, by analyzing polarization, we get *χ* ∼ *e*^{−t/τ}; on the other hand, by analyzing the angular correlation, the reorientation time $\tau \u0303$ is obtained from $\u27e8cos\u2061\theta (t)\u27e9\u223ce\u2212t/\tau \u0303$. From fluctuation–dissipation, we conclude that $e\u2212t/\tau \u223ce\u2212t/\tau \u0303$ and the time evolution should be similar $\tau \u223c\tau \u0303$.

Both times have a maximum at the orthorhombic-to-tetragonal transition temperature and decreases at higher *T* ($\u223c\u2009\u2009e1T$), while the reorientation kinetics become faster. At 300 K, both values are ∼1.5 ps. Having clarified the timescale of molecular reorientations (i.e., relaxations), we aim to identify the associated features in the dielectric function. The above analysis based on the Fourier transform (see Fig. 4) is not suitable as it does not provide sufficient resolution for the dielectric function at low frequencies. As an alternative method, we apply a sinusoidal electric field of frequency $\nu $ to the crystal and we calculate the corresponding component of (complex) susceptibility *χ*(*ν*),

The complex spectral dielectric function is easily derived *ε* = *χ*′ + 1 + *iχ*″. An example of *P*(*t*) response to sinusoids are reported in Fig. 6 and compared to the case of square wave. Calculations were performed by using 32 f.u. crystals and varying the electric field frequency within 40 GHz–10 THz (periods 240 ps–0.1 ps).

The results for imaginary and real parts of the dielectric function are reported in Figs. 7 and 8, respectively. The imaginary part *ε*″ (i.e., related to optical absorption) at low temperature (blue triangles) is almost flat but for the high frequency region 1–10 THz where the peaks of polar optical phonons are found, consistent with Fig. 4. The flat region of the orthorhombic complex dielectric function at lower frequency is consistent with the fact that molecules cannot reorient in the orthorhombic phase. At 326 K (purple squares), the peaks of polar optical phonons of the orthorhombic crystal are still present but sizably broadened, as expected.^{24} The major qualitative difference between *ε*″ of the tetragonal (purple) and orthorhombic (blue) phase is the presence of a band between 0.1–1 THz (1–10 ps), which is absent in the orthorhombic crystal. As *T* decreases toward the orthogonal-to-tetragonal transition, the new band in *ε*″ shifts toward lower frequencies (0.01–0.1) and broadens. We also observe that the shoulders in *ε*′ move to lower frequencies. We can explain qualitatively these features within the Debye relaxation model.^{7} The susceptibility of a dissipative system [Eq. (3)] with relaxation time *τ* is

where the spectral susceptibility is complex *χ* = *χ*′ + *iχ*″. The frequency $\nu \xaf\u2009=\u20091/\tau \xaf\u2009=\u20091/(2\pi \tau )$ gives the maximum of the imaginary part of the susceptibility (i.e., optical absorption) $|\chi \u2033|\u2009=\u2009(\nu /\nu \xaf)[1+(\nu /\nu \xaf)2]\u22121$ and the Half Width at Half Maximum (HWHM) of *χ*′. We point out that the period corresponding to maximum absorption is larger by a factor 2*π* with respect to the relaxation time. The correspondence between the relaxation time *τ* and the period $\tau \xaf$ is shown in Fig. 6 for *τ* = 6.3 ps and $\tau \xaf=2\pi \tau =40$ ps. The sinusoidal *P*(*t*) corresponding to the 40 fs-period electric field approximates the exponential *P*(*t*) transient with relaxation time 40/(2*π*) = 6.3 ps. The relaxation model makes it possible to interpret qualitatively the atomistic data. Relaxation times due to methylammonium rotations depend on temperature *τ*(*T*) and varies in the range of 1–5 ps. The corresponding absorption is expected at $\nu \xaf=1/(2\pi \tau (T))$, which corresponds to 0.03–0.1 THz (i.e., periods $1/\nu \xaf$ of 6–38 ps). This is consistent with the calculated *ε*″ (see Fig. 7). The real part *ε*′ of the susceptibility is also consistent with the Debye relaxation model. At low frequencies, *ε*′ is almost flat corresponding to the static dielectric constant; at higher frequencies, *ε*′ decreases with a shoulder at $\nu \xaf$. This calculated data reproduces qualitatively the expected behavior (see Fig. 8): the shoulder is within 0.03–0.1 THz for *T* = 186 K and 209 K and within 0.1–1 THz for T = 326 K. These values roughly correspond to the expectation based on relaxation times of methylammonium rotations $\nu \xaf=1/(2\pi \tau (T)$. We also observe that the static dielectric constant obtained in Fig. 8 is consistent with data reported in Fig. 5. At 326 K, for example, the ratio between tetragonal and orthorhombic is 3.5.

Overall, we conclude that the dielectric features of the tetragonal phase of MAPbI_{3} below THz frequencies correspond to dissipative relaxation events related to molecular rotations. The present analysis is consistent with experimental data discussed in Ref. 8.

## V. CONCLUSIONS

By using the simple ionic interatomic MYP force-field, we calculate the dielectric function of MAPbI_{3} hybrid perovskites. Although MD underestimates the absolute experimental value of the dielectric function, it is nevertheless possible to reproduce the qualitative features of the dielectric function clarifying the role of molecular orientations. The model reproduces the discontinuity at the onset of the tetragonal phase and the temperature dependence. Furthermore, the model makes it possible to describe the dependence on frequency of both real and imaginary parts of the dielectric function. It is possible in this way to explain the broad absorption bands associated with molecular relaxations in the GHz–THz region, well below the polar optical phonons and the internal IR modes of the molecules. We also clarify the relation between the frequencies of maximum absorption and the relaxation times obtained from polarization or the angular correlation function. The model provides a comprehensive description of polarization and the dielectric function consistent with available experimental and theoretical data.

## ACKNOWLEDGMENTS

We acknowledge MIUR for funding the project (Grant No. PON04a2 00490), M2M Netergit and PRACE for awarding access to Marconi KNL at CINECA, Italy, for projects UNWRAP (Grant No. 2016153664), DECONVOLVES (Grant No. 2018184466), and ISCRA B (MITOMASC).

## REFERENCES

_{3}, NH

_{3}, PbX

_{3}, (X = Cl, Br, I)

_{3}(X = Br and I)

_{3}

_{3}NH

_{3}PbI

_{3}missolution my liquid water

_{3}NH

_{3}PbI

_{3}perovskite

_{3}NH

_{3}PbI

_{3}from theory and experiment: Factor group Analysis, first-principles calculations, and low-temperature infrared spectra

_{3}NH

_{3}PbI

_{3}hybrid perovskite: Interplay of theory and experiment