Quantifying Dynamic Tilting in Halide Perovskites: Chemical Trends and Local Correlations

Halide perovskites have emerged as one of the most interesting materials for optoelectronic applications due to their favorable properties, such as defect-tolerance and long charge carrier lifetimes, which are attributed to their dynamic softness. However, this softness has led to apparent disagreements between the local instantaneous and global average structures of these materials. In this work, we assess the local tilt angles of octahedra in the perovskite structure through large-scale molecular dynamics simulations using machine learned potentials based on density functional theory. We compare structural properties given by different density functionals, namely PBE, PBEsol, SCAN, and vdW-DF-cx, and establish trends across a family of CsMX3 with M=Sn or Pb and X=Cl, Br or I perovskites. Notably, we demonstrate a strong short-range ordering that persists even in the cubic phase of halide perovskites. This ordering is reminiscent of the tetragonal phase and bridges the disordered local structure and the global cubic arrangement. Our results provide a deeper understanding of the structural properties of halide perovskites and their local distortions, which is crucial for further understanding their optoelectronic properties.


I. INTRODUCTION
Halide perovskites have gained significant attention as promising materials for various applications, including high-efficiency solar cells [1][2][3], lasers [4], light-emitting diodes [5], and more [6].Their exceptional performance is attributed to various factors, such as their defect tolerance [7][8][9] and long carrier lifetimes [10][11][12][13].These properties are linked to the high dielectric constants of these materials and effective screening of charges [14,15], which in turn can be traced back to their dynamic softness and ability to dynamically respond to the presence of excess charges.Despite their importance, understanding and quantifying the underlying dynamics of halide perovskites remains challenging.
Several previous studies have highlighted the crucial role played by the dynamic local structure and octahedral tilting in halide perovskites [16].For example, some investigations have revealed that ignoring the local structure can lead to inaccurate interpretation of various experimental techniques, including X-ray absorption near edge structure spectroscopy [17], high energy resolution inelastic X-ray scattering [18], and pair distribution function analysis [19][20][21].Moreover, electronic structure calculations performed for simulations cells corresponding to the average perovskite structure have been shown to result in poor estimates of band gaps [21][22][23][24], emphasizing the criticality of accounting for transient local distortions.Moreover, octahedral tilting has been demonstrated to be vital in stabilizing the perovskite phase in FAPbI 3 [25], further highlighting the significance of this property.These findings underscore the importance of understanding the dynamic nature of halide perovskites and its implications for a wide range of applications.
Octahedral tilting, being one of the key properties of halide perovskites, has been assessed in several previous computational studies.One of the simplest indicators of how much a materials will distort at finite temperatures can already be extracted from static calculations of potential energy surfaces or phonon dispersion curves at 0 K, as has been done for instance in Refs.26, 27.This way, one can also compare the performance of various functionals [28], but it is not possible to predict the exact dynamic properties of a material from calculations done at 0 K alone.Molecular dynamics (MD) simulations have also been used to study the dynamics of halide perovskites more directly, for instance to contrast different functionals [29,30], to verify the effect of cation mixing [31] or to understand the local disorder and its effect on experimental observations [17,32].Various dynamic properties have been used in these studies to quantify the disorder, such as evolution of a single M -X-M angle [31], distributions of M -X-M tilt angles [17,33] or distributions of halogen displacements [29,33].
In this work, we quantify the local octahedral tilting of halide perovskites through large-scale MD simulations and assess trends across functionals, temperatures, and material chemistries.We use neuroevolution potential (NEP) models to accelerate MD simulations based on density functional theory (DFT) calculations.In the case of α-CsPbI 3 , we find that different functionals lead to similar trends in tilt angle distributions within the same phase, but drastically different phase transition temperatures.Focusing further on the vander-Waals-density functional with consistent exchange (vdW-DF-cx) functional, we assess the trends across the CsM X family of inorganic perovskites with varying halogen atoms (X=Cl, Br, and I) and metal cation (M =Sn, and Pb).We finally create a bridge between the local environment of a single octahedron and the global perovskite structure through assessing the correlation be-tween neighboring units.

A. Machine-learned potential construction
We constructed third-generation (NEP3) NEP models [34,35] using a boot-strapping strategy, as described in detail in Ref. 30.In this process, we included the cubic (P m 3m), two tetragonal (I4/mcm, and P 4/mbm), and one orthorhombic (P nma) structure, as well as the so-called delta-phase (P nma).The model construction was carried out using the gpumd package [35] as well as the calorine package [36] and hiphive packages [37] for data preparation and analysis.The performance of the final models is summarized in Table S1 of the Supplemental Material (SM).

B. Density functional theory calculations
To generate input data for model construction, we performed DFT calculations using the projector augmentedwave method (PAW) method [38,39] as implemented in the Vienna ab-initio simulation package [40,41].Our calculations use a Γ-centered k-point grid with density of 0.18 Å−1 and a Gaussian smearing width of 0.1 eV.Four different exchange-correlation functionals were used, namely PBE [42], PBSEsol [43], strongly constrained and appropriately normed (SCAN) [44], and the vdW-DF-cx method [45,46].The valence electron configurations used in the PAW datasets are provided in Table S2 in the SM.We used a plane-wave energy cutoff to 520 eV for all materials.
C. MD simulations MD simulations were carried out using the gpumd package [35,47].During the simulations the temperature was continuously increased from 20 to up to 620 K over a period of 100 ns, yielding to a heating rate of 6 K ns −1 .For some compounds with low transition temperatures we terminated the simulations already at 420 or 520 K.The simulations were run in the N pT ensemble, with temperature and pressure controlled using stochastic velocity [48] and cell rescaling [49].The time step is set to 5 ps, and the pressure to zero.The simulation cell comprised 16 × 16 × 12 repetitions of the primitive unit cell of the orthorhombic perovskite structure, corresponding to 61 440 atoms.

D. Analysis of tilt angles
To extract the tilt angle for each M X 6 octahedron we applied polyhedral template matching using ovito  [50].The orientation of the octahedra is then represented by Euler angles as shown in Figure 1.This is done at each step of the heating runs, giving us access to the angle distribution as a practically continuous function of temperature.

A. Comparison of functionals
We first compare the distributions of octahedral tilt angles computed for CsPbI 3 according to NEP models representing four different functionals: the widely used PBE functional [42], its variation adapted to modelling of solids, PBEsol [43], the SCAN functional [44], which has shown good performance for halide perovskites in comparison with calculations based on the random phase approximation [51], and the vdW-DF-cx method [45,46].The distributions for all temperatures are shown in Figure 1 as maps of tilting angles.The transition temperatures from the orthorhombic (P nma) to the tetragonal (P 4/mbm) phase at lower temperatures and from tetragonal to cubic (P m 3m)) at higher temperatures as obtained in Ref. 30 are also marked.Phase transition temperatures based on the models for all four functionals are given in Table I and compared to experimental values.
At low temperatures, the angle splitting are similar for all models (functionals).The θ = φ angle, characterizing the tilt in the z direction in the orthorhombic structure has a magnitude of just below 10°.The ψ angle, related to the octahedral rotation in the x − y plane is slightly larger than 10°.The temperature evolution of the angles is also qualitatively similar, with the main difference being related to phase transition temperatures (see Table I), with vdW-DF-cx giving the highest transition temperatures (the closest to experimental values) and PBE giving the lowest.We also note that the ψ angle in the orthorhombic phase has a stronger temperature dependence for the models based on the PBEsol and vdW-DF-cx functionals, than in the case of PBE and SCAN.
To enable a more straightforward comparison of angular distributions between the models representing the different functionals, we now analyze tilting at three selected temperatures, 150, 300, and 500 K (Figure 2).At 150 K, all models (functionals) correctly predict CsPbI 3 to exhibit the orthorhombic structure.While the θ = φ distributions are almost the same for all functionals, the separation of the maxima in the ψ distribution clearly increases in the order PBE−SCAN−PBEsol−vdW-DF-cx.
Next we focus on 300 K, corresponding to the most commonly value used in MD simulations for halide perovskites.At this temperature, CsPbI 3 should exhibit the orthorhombic structure, according to experiments [52].Strikingly, among the considered functionals, only the vdW-DF-cx-based model predicts this phase to be stable at 300 K, which is reflected in angle splitting for all three Euler angles.With the PBE, PBESol, and SCAN-based models we observe θ = φ distributions of similar broadening with one maximum only, indicating the tetragonal phase.For the ψ angle, all functionals yield a bimodal distribution.The splitting of the maxima is very similar for all models (functionals) except for PBE, in which it is significantly smaller.The differences between the distributions of θ = φ achieved with different functionals highlight the importance of choosing the right computational setting for constant-temperature MD simulations.We thus expect that room-temperature simulations of CsPbI 3 in the orthorhombic phase using PBE, PBESol or SCAN could lead to spurious distortions in the cell due to the drive towards the tetragonal phase.
At 500 K, all functionals predict the same cubic phase, leading to similar angular distributions.For all three Euler angles, we observe a broad Gaussian distribution centered around 0°.We note that such a broad distribution of octahedral tilting angles is the reason for a large difference in the local and global symmetry in the halide perovskites [17,22].While on average the structure at high temperatures exhibits a high-symmetry cubic structure, locally each octahedron is likely to be significantly tilted, with a standard deviation (see Figure S4 in SM) from the average position of about 7°in each direction.
To further analyze the data, we fit the angle distributions with a symmetric double Gaussian function of the form , where µ and σ are mean and standard deviation respectively.The fitting is done with a small regularization for the free parameters µ and σ to specifically reduce noise when fitting the distributions from the cubic phase.
The positions of the maxima in the distributions (µ) as a function of temperature are plotted in Figure 2b.The figure indicates that the splitting of the tilt angles exhibits a similar temperature dependence for all models (functionals).We also note that the abrupt change in the positions of the maxima in the distributions agrees very well with phase transition temperatures (indicated by dashed lines in the figure) extracted from heat capacities [30].We also notice a clear correlation between the magnitude of the splitting in the angular distribution at low temperatures and the transition temperatures, with higher splittings corresponding to higher transition temperatures.

B. Chemical trends
We now focus on the vdW-DF-cx functional, which exhibits the closest agreement with experiment with respect to transition temperatures and thermal expansion [30], and compare different compounds.We first analyze the effect of varying the halogen atom within the family of compounds CsPbX 3 with X =Cl, Br, and I. Tilting angle maps as well as transition temperatures can be found in the SM, while in Figure 3 we show the distributions at three temperatures (150, 300, and 400 K) as well as the positions of maxima in the angular distributions.At low temperatures, all three materials exhibit the orthorhombic P nma phase consistent with bimodal distributions of all three Euler angles.The splitting between the maxima of both the θ = φ and ψ angles increases with the size of the halogen ion, following the sequence Cl-Br-I.As can be seen in Figure 3b, this ordering is retained at all temperatures.
At 300 K, CsPbCl 3 undergoes a phase transition from the tetragonal to the cubic phase, which explains the mixed nature of the ψ angle distribution.CsPbBr 3 is stable in the tetragonal phase, with a unimodal distribution of θ = φ and a bimodal distribution of ψ.CsPbI 3 is still orthorhombic at 300 K and exhibits bimodal distribution of all Euler angles.We note that for CsPbCl 3 , 300 K is exactly the phase transition temperature according to the vdW-DF-cx-based model, which explains the mixed tetragonal/cubic distribution of the ψ tilting angle.At 400 K, CsPbI 3 is tetragonal, while both CsPbBr 3 and CsPbCl 3 are cubic.This allows for a direct comparison between the later two.While the plots in Figure 3a are rather similar, the distribution for CsPbCl 3 is slightly more narrow, suggesting that the octahedral tilting is reduced for smaller halogen ions.
We now evaluate the effect of the metal cation, by making a comparison between CsPbI 3 and CsSnI 3 (Figure 4).At 150 K, both CsPbI 3 and CsSnI 3 adopt the orthorhombic phase.The tilt angle splitting is much more pronounced in CsPbI 3 than in CsSnBr 3 .The smaller magnitude of tilt angles at low temperatures for CsSnI 3 is correlated with much lower phase transition temperatures as compared to CsPbI 3 (see Table S3 in the SM).

C. Local tilt angle correlations
From analyzing the tilt angle distributions, it is clear that at high temperatures halide perovskites are cubic on a global scale corresponding to µ = 0.However, locally each octahedron can still be significantly tilted, with tilt angles reaching the values seen in the tetragonal and orthorhombic phases.In the tetragonal and orthorhombic phases long-range ordering occurs with respect to the octahedral tilt angles.In the cubic phase there is no longrange order, this does, however, not imply the absence of short-range order (correlation) between neighboring octahedra.In the context of local regions in the cubic phase appearing tetragonal or orthorhombic [17], it is thus of interest to analyze and quantify the short-range order, i.e., assess the correlation between tilt angles of neighboring octahedra, in the cubic phase.Strikingly, a recent study by Doherty et al. found that phases that in X-ray diffraction appear to be cubic, locally feature regions that are tetragonal or orthorhombic [25].Here, we investigate the joint distribution of the tilt angle of neighboring octahedra, P (ψ 1 , ψ n ), where n refers to the n-th neighbor shell in the [100] direction (perpendicular to the rotational axis of ψ) and ψ the tilt-angle around the z-axis (Figure 5a).Note that the same correlation occurs for the symmetrically equivalent [010] direction.This is analysis is carried out for CsPbI 3 , but the behavior is qualitatively the same for all materials considered here (Figure S5 in the SM).There is a clear correlation between ψ 1 and ψ n for small n, which can be quantified with the Pearson correlation p n defined as The correlation naturally occurs due to the nearest neighbor octahedra sharing one halogen ion (X=I, Br, Cl) atom (Figure 1a).If a given octahedron rotates in one direction its neighbors will thus favorably rotate in the opposite direction.This leads to the sign of the correlation alternating between odd and even neighbor shells along the [100] direction.The absolute value of the correlation |p n | decays exponentially with a typical length-scale of about 3 to 5 shells (corresponding to about 20 Å) with the number of the neighbor shell (Figure 5b).We note that at 500 K (and thus above the tetragonal-to-cubic transition temperature), there is a significant correlation between octahedra lying multiple shells apart, stretching as far as the seventh neighbor shell.Although the correlation weakens at temperatures between 600 and 700 Kelvin, which is far into the cubic stability region, the orientation of a specific octahedron still has a notable effect on its second and third neighbors.
The full temperature dependence of the correlation p n (Figure 5c) shows that the correlation increases as the temperature decreases towards the transition to the tetragonal phase.While the tilt angle around the z-axis (ψ) correlates strongly between neighboring octahedra along the [100] direction, there is almost no correlation in the cubic phase between neighboring octahedra along the [001] direction.[53]This is likely due to the fact that the shared X cation for these neighbors does not move when rotating around the z-axis, and thus the octahedra can rotate around the z-axis more independently of each other.In the tetragonal phase the correlation between all neighboring shells approaches the respective limiting values of ±1, indicating that long-ranged ordering is obtained.
The strong short-range ordering observed in the cubic phase may have implication for refinement of experimental diffraction data.Since some of the tetragonal ordering is still present in the cubic phase, the structure determination can be difficult in these materials.Specifically, our findings put into question the frequently adopted approach of modelling thermal displacement parameters through Debye-Waller factors in this class of materials.The latter (as commonly used) incorporate namely the assumption that the local energy landscape can be approximated as harmonic.Halide perovskites, however, exhibit extreme anharmonicity [33,52], in particular for the zone boundary modes that are connected to the octahedral tilting [32].

IV. CONCLUSIONS
In conclusion, in this study we have explored the local structural properties of inorganic halide perovskites through large-scale MD simulations using machinelearned potentials based on DFT calculations.By assessing the phase transitions and local tilt angles of octahedra in the perovskite structure, we first analyzed the performance of different functionals as described by NEP models, including PBE, PBEsol, SCAN, and vdW-DFcx.We find that all models (functionals) underestimate the phase transition temperatures, with the vdW-DFcx-based model yielding the values that are closest to experiment.Focusing subsequently on the vdW-DF-cx method, we then established trends across a family of CsM X 3 perovskites.Finally, we demonstrated the presence of strong short-range ordering, reminiscent of the tetragonal phase, even in the cubic phase of halide perovskites for above the cubic-to-tetragonal transition temperature.This ordering connects the transiently strongly distorted local structure to the global cubic arrangement and sheds light on the complex interplay between local and global structural properties in these materials.This study provides direct insight into the local disorder in inorganic halide perovskites across a wide temperature range and contributes to the more complete understanding of their structural properties.

FIG. 1
FIG. 1.(a) Representation of the CsPbI 3 perovskite structure with PbI 6 octahedra indicating the definition of the three Euler angles θ, φ, and ψ describing the orientation of the octahedron.(b) Maps of tilt angles as a function of temperature in PBE, PBEsol, SCAN, and vdW-DF-cx functionals.Dashed vertical lines indicate the orthorhombic-to-tetragonal and tetragonal-tocubic phase transitions.
FIG. 2. (a) Probability distribution of octahedral tilts in CsPbI 3 as described by the three Euler angles θ, φ, and ψ.Data obtained with NEP models based on four different functionals (PBE, PBEsol, SCAN, and vdW-DF-cx) are shown at three different temperatures.Solid lines are double Gaussian fits to the data.(b) Position of the maximum in the tilt angle distribution as extracted from double Gaussian fits.Vertical dashed lines indicate phase transitions.

FIG. 3 .FIG. 4 .
FIG. 3. (a) Probability distribution of octahedral tilts at 150, 300, and 400 K in CsPbCl 3 , CsPbBr 3 , and CsPbI 3 , as described by the three Euler angles θ, φ, and ψ.(b) Position of the maximum in the tilt angle distribution as extracted from double Gaussian fits.Vertical dashed lines indicate phase transitions.

TABLE I .
[52]e transition temperatures in K found using different functionals for CsPbI 3 .Values were extracted from the temperature dependence of the heat capacity[30].Experimental values are given for comparison[52].