Revealing the free energy landscape of halide perovskites: Metastability and transition characters in CsPbBr 3 and MAPbI 3

Halide perovskites have emerged as a promising class of materials for photovoltaic applications. A challenge in these applications is how to prevent the crystal structure from degradation to photo-voltaically inactive phases, which requires an understanding of the free energy landscape of these materials. Here, we uncover the free energy landscape of two prototypical halide perovskites, CsPbBr 3 and MAPbI 3 via atomic scale simulations using umbrella sampling and machine-learned potentials. For CsPbBr 3 we find very small free energy differences and barriers close to the transition temperatures for both the tetragonal-to-cubic and the orthorhombic-to-tetragonal transition. For MAPbI 3 , however, the situation is more intricate. In particular the orthorhombic-to-tetragonal transition exhibits a large free energy barrier and there are several competing tetragonal phases. Using large-scale molecular dynamics simulations we explore the character of these transition and observe latent heat and a discrete change in structural parameters for the tetragonal-to-cubic phase transition in both CsPbBr 3 and MAPbI 3 indicating first-order transitions. We find that in MAPbI 3 the orthorhombic phase has an extended metastability range and furthermore identify a second metastable tetragonal phase. Finally, we compile a phase diagram for MAPbI 3 that includes potential metastable phases.

In spite of the very large body of literature available, halide perovskites continue to surprise.For example, a previously not experimentally recognized crystal structure of the the widely studied perovskite CsPbBr 3 , was detected in a combination of convergent beam electron diffraction and electron ptychography [11].Additionally, several recent studies showed that short-range ordering in octahedral tilt angles characteristic of low-temperature phases persists even in the cubic phase of halide perovskites [6,7,12].Furthermore, the widely accepted low-temperature Pnma structure of MAPbI 3 has been put into question by a combination of inelastic neutron scattering and first-principles calculations [13].As an alternative the authors proposed that structures with lower local symmetry are needed to reproduce experimental data [13] and in fact several halide perovskite structures have been shown to be locally rather similar [14].The structural models differ, however, in properties that are crucial for photovoltaic applications, most importantly the band gap [15].Therefore, a more detailed understanding of the local atomic structure and dynamics as well as their coupling to phase transitions is needed for further optimizing these materials.
Similarly, the character of the phase transitions in halide perovskite has been debated.Both the tetragonalto-cubic and orthorhombic-to-tetragonal phase transitions are driven by soft phonon modes, specifically shear modes at the M and R-points on the Brillouin zone boundary of the cubic structure, which in the following we simply refer to as M and R modes for brevity.For a purely continuous transition the frequency of the soft modes should go smoothly to zero at the transition temperature.In CsPbBr 3 the tetragonal-to-cubic phase transition has, however, been observed to have first-order character through experimental studies [16][17][18], whereas the transition between the tetragonal and orthorhombic phases has been characterized as purely continuous [16].
In MAPbI 3 the character of the tetragonal-to-cubic transition has been analyzed in numerous experimental studies with different conclusions concerning the order of the transition [19][20][21][22][23], although most studies point towards it being first order.The orthorhombic-to-tetragonal transition in MAPbI 3 has been consistently observed to have first-order character.
One should also note that mixing both on the cation and anion sites involving two or more different species has been demonstrated as a promising route for improving the thermodynamic stability of these materials [24][25][26][27].This renders the configuration (and phase) space of these materials even more complex as it adds further compositional degrees of freedom.
From a computational standpoint, the very dynamic crystal structure of halide perovskites [28,29] presents a challenge.Specifically, vibrational properties play a key role in determining the free energy landscape and ground state density functional theory (DFT) calculations provide only very limited information in this regard.To capture the dynamic properties, ab-initio molecular dynamics (MD) simulations are a popular approach, and have shone light on, among other things, the band gaps of halide perovskites [30][31][32][33][34], local disorder [14,35], and defect energy levels [36].Yet, the computational cost of ab-initio MD simulations prevents extensive sampling and thereby limits access to structural and dynamic correlations that require careful convergence.Since machinelearned potentials (MLPs) have the ability to significantly speed up sampling with a negligible loss in accuracy, there has been a number of applications of these techniques to halide perovskites [4,7,[37][38][39][40][41][42][43][44][45][46].For example, Jinnouchi et al. [4] constructed a MLP for MAPbI 3 to study the phase transitions [4].However, the previous studies focused on specific, predetermined phases or transition between those.Therefore, our understanding of the free energy landscape of halide perovskites is still evolving.
Here, we uncover the free energy landscape of two halide perovskites, CsPbBr 3 and MAPbI 3 with respect to the tilting degrees of freedom, by constructing MLP that are both accurate and computationally very efficient on graphics processing units (GPUs) [46,47].The structural phase space is explored via MD simulations and umbrella sampling [48].We demonstrate that the free energy minima in CsPbBr 3 for the various phases are typically wide and soft, with small free energy differences between phases.For MAPbI 3 the free energy landscape indicates a large barrier between the orthorhombic and tetragonal phases as well as several competing tetragonal phases.We relate these results to the available tilt modes by analyzing the space of accessible tilt patterns using the concept of "Glazer space".This provides a geometrical view of the possible transitions that is transferable beyond the specific compounds and chemistries considered here.Our analysis enables us to clearly resolve the character of the different phase transitions in these systems in agreement with the majority of experimental assignments and eventually to construct a phase diagram for MAPbI 3 that includes potential metastable phases.

A. Density functional theory calculations
DFT calculations were performed using the projector augmented-wave method [49] as implemented in the Vienna ab-initio simulation package [50,51].The exchange-correlation contribution was represented using the strongly constrained and appropriately normed (SCAN) density functional [52], which has been previously established for the study of these systems [4,53].The Brillouin zone was sampled using a Γ-centered grid with a k-point density of 0.25 Å−1 and Gaussian smearing with a width of 0.1 eV.

B. Neuroevolution potentials
For CsPbBr 3 we employed the neuroevolution potential (NEP) model described in our recent work [45], and for MAPbI 3 we constructed a model by similar means using the gpumd package [47,48,54] in conjunction with the calorine [55] and ase packages [56].The models are based on a neural network for which local atomic environments are described by radial and angular components [54].In the NEP3 version [47] employed here the radial part of the atomic environment descriptor is constructed from linear combinations of Chebyshev basis functions while the three-body angular part is similarly build from Legendre polynomials.
Here, we used radial and angular cutoffs of 8 Å and 4 Å, respectively, and terminated the radial and angular expansions at orders 12 and 6, respectively.The neural network consisted of one hidden layer of 40 neurons with a hyperbolic tangent activation function and was trained using the natural evolution strategy [57] implemented in gpumd, using 200 000 and 500 000 generations for training the CsPbBr 3 and MAPbI 3 models, respectively.Both ℓ 1 and the ℓ 2 regularization terms were included, using λ 1 = λ 2 = 0.02 for CsPbBr 3 and λ 1 = λ 2 = 0.05 for MAPbI 3 (with λ 1 and λ 2 as defined in Ref. 47).Structures for training and validation were generated by means of a bootstrapping approach, i.e., first preliminary models were fitted, new structures were extracted from MD simulations based on these models, and then used to fit the next model generation, similar to the approach outlined in Ref. 46.The final model for CsPbBr 3 was fitted to forces, energies, and virials from DFT calculations for 642 atomic structures with a total of 176 920 atoms, and validated against a set of 72 structures with a total of 19 800 atoms.In this holdout set, the model achieved an root mean square (RMS) error of 0.6 meV atom −1 for energies, 15 meV atom −1 for virials, and 43 meV Å−1 for forces.The MAPbI 3 model was fitted to forces, energies, and virials from DFT calculation of 511 atomic structures with a total of 172 092 atoms, and validated against a set of 57 structures with a total of 18 804 atoms.It achieved a RMS error of 0.6 meV atom −1 for energies, 4.5 meV atom −1 for virials, and 38 meV Å−1 for forces (see Supporting Information (SI) for details).

C. Molecular dynamics and umbrella sampling
MD simulations were carried out using the gpumd package (version 3.2) [48] using a custom patch for umbrella sampling.In these simulations, time was propagated in steps of 0.5 fs in systems containing hydrogen, otherwise a time step of 1 fs was employed.The Berendsen thermostat and barostat [59] were used to control temperature and pressure.Systems were initialized from 4 × 4 × 3 supercells of the ideal primitive orthorhombic structure, for a total of 192 formula units.The system size was chosen to be large enough to avoid finite-size effects but small enough to exhibit the statistical fluctuations required for free energy integration with weighted histogram analysis method (WHAM) (see below).Furthermore, the global constraints in the free energy sam- pling become weaker for larger system sizes and the system may not stay in the desired crystal structure but instead form interfaces or domains.The cell was constrained to remain orthorhombic, i.e., the volume was free to fluctuate but the angles between cell vector were fixed to 90 • , throughout the simulations to simplify the analysis.
Umbrella sampling [60] was used to extract the free energy in "Glazer space" (Fig. 1a and Sect.III A).To this end, we biased the simulations using six collective variables: projections ξ of atomic displacements on the three M modes and the three R modes as done in Ref. 45.Each phonon mode projection was biased with a quadratic energy term, where ξ mode is the projection on the respective phonon mode displacement (with the full displacement vector normalized to 1), ξmode is a specific target projection, and the spring constant k = 0.1 eV −1 .We performed one MD simulation of 1 ns for each combination of values of ξmode in order to span the Glazer space along the diagonal x = y.If the target projection of any of the modes M i was non-zero, the corresponding R i target projection was set to zero.Free energies were calculated with WHAM [61] using the software provided by Grossfield [62] but extended to handle the three-dimensional case.At low temperatures (≤ 200 K for CsPbBr 3 , ≤ 50 K for MAPbI 3 ) we resorted to umbrella integration since fluctuations were too small to achieve sufficiently overlapping histograms from individual simulations.
Mapping the free energy across the Glazer space thus allows one to identify the phases corresponding to the free energy minima, their respective tilt angles, and gain insight regarding the barriers between phases [63].

A. The Glazer space
The perovskite crystal structure has the general formula ABX 3 .The B atoms reside in corner-sharing octahedra formed by the X atoms, which in the case of halide perovskites are halide species.These octahedra tilt according to patterns that are distinctive for each phase.Typically, there are at least three phases that are stable in different temperature intervals: an orthorhombic ground state, a tetragonal phase, and a cubic phase.The tilt patterns are associated with phonon instabilities.In this regard the three M modes and the three R modes of the cubic structure are particularly important [45].Each M mode corresponds to the rotation of an octahedron around one of the three major axes, with consecutive octahedra along that axis rotated in-phase.The R modes have a similar effect but the octahedra rotate out-of-phase.The different perovskite crystal structures are conveniently described by the projections on these modes, for which Glazer introduced a compact notation [64].For example, the commonly occurring Pnma phase can be written as a − a − c + , which corresponds to a structure with an out-of-phase rotation of a degrees about both the x and y axes (the R x and R y modes) and an in-phase rotation of c degrees about the z axis (M z ).The cubic phase, on the other hand, corresponds to a 0 a 0 a 0 (no average rotation about any axis).We will refer to the three-dimensional space spanned by in-phase and out-of-phase rotation about the three axes as the Glazer space (see Fig. 1a for a schematic).All known, stable perovskite phases of CsPbBr 3 and MAPbI 3 correspond to distinct regions in Glazer space.

B. CsPbBr3
Behavior on heating/cooling.First, we run heating and cooling MD simulations in the N P T ensemble for CsPbBr 3 with 61 440 atoms and a cooling/heating rate of about 10 K ns −1 (Fig. 2a,b).Two phase transitions are observed, at around 260 K and 300 K, corresponding to the orthorhombic-to-tetragonal transition and the tetragonal-to-cubic transition, respectively.The transitions are clearly visible as changes in the lattice parameters, but can also be determined from peaks and kinks in the heat capacity [46].For the tetragonal-tocubic transition some hysteresis is observed suggesting first-order character, whereas for the orthorhombic-totetragonal transition no hysteresis is observed, fully in line with a continuous phase transition.
Free energy landscape in Glazer space.It is now instructive to consider the temperature dependence of the free energy landscape in Glazer space (Fig. 1).At 350 K the free energy minimum is found at the origin of the Glazer space, i.e., the cubic a 0 a 0 a 0 structure (Pm 3m).At 300 K the free energy is lowest for the tetragonal a 0 a 0 c + structure (P4/mbm), although the free energy landscape is very flat along the z-tilt axis towards the cubic phase.The tilt angle of the tetragonal structure is about 8 • .At 250 K the free energy minimum is found for the orthorhombic a − a − c + structure (Pnma).Finally, at 200 K a distinct narrow free energy minimum is obtained for the a − a − c + phase, which continues to be the stable phase down to 0 K.The tilt angle about the z-axis is the same as in the tetragonal structure, while the tilt angle about x and y is about 4 • .Again, we note that the free energy landscape between the orthorhombic (a − a − c + ) and the tetragonal (a 0 a 0 c + ) phase is very flat.The local free energy minimum closer to the center of Glazer space corresponds to the a 0 a 0 c + phase but with multiple domains separated by anti-phase boundaries, which arises due to the global constraint on the phonon mode coordinate.These can for example occur in a 6 × 6 × 6 (cubic) supercell with four layers of primitive cells tilting with the same phase, while the other two tilt with the opposite phase.The resulting structure has a global order parameter (mode coordinate) equal to 1/3 of the a 0 a 0 c + ground-state.
The calculated free energy in Glazer space contains the correct phase transitions as observed in experiments [16][17][18][65][66][67][68].They occur, however, at slightly lower temperatures, which can be attributed to the exchangecorrelation functional used in the reference DFT calculations [46].
Character of the transitions.Next, we address the character of the transitions in CsPbBr 3 , which are analyzed via long (20 ns), large-scale (compared to abinitio calculations; 61 440 atoms) MD simulations in the N P T ensemble in the vicinity of the transition temperature.Energy, lattice parameters, and mode projections are recorded over 10 independent simulations for each temperature.The energy referred to here is the potential energy after subtracting the zero Kelvin energy as well as the Dulong-Petit term 3k B T /2.
In the vicinity of the upper transition between the tetragonal and cubic phases, the histograms over the observed energy and lattice parameter values are multimodal (see Fig. 3c,d for examples; see Fig. S3 in the SI for the complete data).This behavior is due to the superposition of the distributions of the cubic and tetragonal phases, and the total histograms can be expressed as a sum of normal distributions (as shown in Fig. 3c,d).We emphasize that the two phases are not observed in the MD simulations simultaneously; the system is rather switching between the two phases over time.
As a result of the multimodal character in the vicinity of the phase transition, the mean of the distribution (circles in Fig. 3a,b), which also corresponds to the time average of the respective quantity, changes gradually across the transition.[69] The mode of the distribution, i.e., the position of the maximum (or most likely value), however, exhibits a discrete jump at the transition temperature (squares in Fig. 3c,d), which is clear evidence for the first-order character of the transition.The magnitude of the change in the potential energy yields a very small but finite latent heat of about 0.6 meV atom −1 .This magnitude is consistent with a first-order transition with a (very) small free energy barrier separating the phases, allowing the system to jump between the two phases at the size and time scale of our MD simulations.This barrier is also clearly visible in the mode coordinate of the M mode (Fig. S2).
By contrast, for the lower (orthorhombic-tetragonal) Lattice parameters and heat capacity from heating and cooling MD simulations in the N P T ensemble for (a,b) CsPbBr3 and (c,d) MAPbI3.The heat capacity is reported per degree of freedom in the system.transition, the distributions are unimodal and mean and mode coincide throughout the transition, which implies continuous character (Fig. 3).
Our simulations thus show that the tetragonal-to-cubic transition is of first order while the orthorhombic-totetragonal is continuous, in agreement with experiments [16,18].This is also consistent with our recent work on the phonon dynamics in this material, where the M-mode frequency in the cubic phase was observed to be non-zero at the tetragonal-to-cubic transition temperature [45].

C. MAPbI3
Behavior on heating/cooling.In the case of MAPbI 3 we observe a transition during cooling at about 300 K from the high-temperature cubic phase (a 0 a 0 a 0 ) to a tetragonal structure (a 0 a 0 c − ; Fig. 2c,d).There are no further transitions as the system stays in the tetragonal phase down to 0 K.
On heating, starting from the orthorhombic groundstate structure (a − a − c + ), a transition occurs at around 150 K to a tetragonal (referring to a = b ̸ = c) structure [70] (a − a − c 0 ) followed by another transition at 300 K to the cubic phase.
The cubic (a 0 a 0 a 0 ) and orthorhombic (a − a − c + ) phases as well as the tetragonal phase obtained on cooling (a 0 a 0 c − ) are observed experimentally, whereas the tetragonal phase observed on heating (a − a − c 0 ) is not.This is a reflection of the stochastic nature of phase transitions.Their outcome depends on the free energy landscape at the transition temperature and is thus sensitive to both the free energy differences between possible struc-tures and the barriers that separate them.To understand the outcome of these simulations and the nature of the phase transitions (and pathways) in MAPbI 3 it is therefore instructive to inspect the free energy landscape as a function of temperature.
Free energy landscape in Glazer space.At 350 K the free energy minimum is located at the center of the Glazer space, corresponding to the cubic phase (a 0 a 0 a 0 , Fig. 1c).At 300 K the minimum has shifted to the tetragonal structure with anti-phase tilting (a 0 a 0 c − ).Similarly to the case of CsPbBr 3 (Fig. 1b), the free energy landscape close to the tetragonal-to-cubic transition is very wide and flat.
Interestingly, we also observe a local minimum in the free energy for a tetragonal a − a − c 0 structure at 100 K. Below 100 K the free energy minimum shifts to the orthorhombic a − a − c + structure.The experimental phase transition to the orthorhombic phase is around 160 K, which is slightly higher than the transition temperature we find here.
Character of the tetragonal-to-cubic transition.To gain further insight, we analyze the character of the tetragonal-to-cubic transition in MAPbI 3 in the same manner as for CsPbBr 3 (Fig. 3).We observe a discrete change in both energy (latent heat) and lattice parameters which indicates a first-order phase transitions.The full distributions and corresponding fits are shown in Fig. S4 of the SI.The very small free energy barrier between these two phases makes the transition appear more continuous in MD simulations.This is in agreement with the near tricritical character of this transition described by Whitfield et al., who based on experimental data reported coexistence between the cubic and tetrag- onal phase in the temperature range between 300 and 330 K [23].Here, we observe both phases in MD simulations in a narrow temperature interval between 305 and 310 K.
Molecular order and rotations.Replacing the Cs atom with the organic MA molecule gives rise to additional degrees of freedom and renders MAPI significantly more complex than CsPbBr 3 .Crucially, the MA units have an orientation (C-N bond-vector), which is associated with ordering, in particular in the orthorhombic phase but also in the tetragonal phases [22,53,[71][72][73][74][75].These orientational degrees of freedom also make it much more challenging to sample the free energy landscape at low temperatures as they lead to energetic barriers between phases, which in turn allows metastable phases to remain dynamically stable down to 0 K.In order to understand the energetics of the MAPbI 3 phases better we therefore conducted an extensive sampling of possible structures (Fig. 4). 10 000 different initial structures were constructed based on 2 × 2 × 2 supercells of the primitive cubic unit cell, by applying a random set of tilt modes (either M or R) for each Cartesian direction with randomized amplitudes, and combined with randomized MA orientations.Each structure was relaxed (atomic positions and cell shape) until the largest force fell below 1 × 10 −4 eV Å−1 .The relaxed structures were then classified into the Glazer structures by projection onto the M and R modes using a numerical tolerance of about 0.5 • .
The first thing to note is that this approach correctly captures the a − a − c + structure as the ground state.Furthermore, it is clear that the group of tetragonal a 0 a 0 c + structures is significantly higher in energy compared to a 0 a 0 c − .This is consistent with the fact that far fewer structures end up in the a 0 a 0 c − structure after relaxation.Yet the majority of structures ends up at slightly higher energies (>20 meV atom −1 ) which illustrates the rather rough energy landscape of this material.
At finite temperature, the orientations of the MA molecules fluctuate and can even flip at a significant rate even at modest temperatures [53,[76][77][78].This gives rise to a sizeable entropic contribution and a smoothing of the free energy landscape with increasing temperature.To analyze this aspect using the present approach, we ran MD simulations at several different temperatures starting from the lowest energy structure of each phase (Fig. 4) in supercells containing 96 000 atoms.Following an equilibration for 1 ns in the N V T ensemble (at the volume obtained previously from N P T runs), then the C-N bond vector r CN (t) of each MA unit was sampled in the N V E ensemble (Fig. S6).We then analyzed the rotational dynamics of MA in terms of the orientational autocorrelation function (ACF) where r i CN (t) is the C-N bond vector at time t for the i-th MA unit.For the orthorhombic structure (a − a − c + ), one observes very few MA reorientation or rotation events on the time scales sampled here, and hence they are not included in the following analysis.The ACF C(τ ) exhibits oscillations on a shorter time scale (about 100 fs) corresponding to atomic vibrations, and a slow decay on a longer time scale (about 10 ps to 100 ps) corresponding to MA reorientation (Fig. 5a,b).The latter slow rotational behavior can be modeled with one exponential decay, C(τ ) ∝ e −τ /τrot , where τ rot corresponds to a typical rotational time (see, e.g., Fig. S5) [53,76].One can model the temperature dependence of the rotation time with an Arrhenius expression, 1/τ rot (T ) ∝ e −E A /kBT , which fits the data very well (Fig. 5c).This yields effective activation barriers E A for these rotational events of 120 meV and 103 meV for the tetragonal and cubic phases, respectively.These activation barriers are in good agreement with experimental measurements [77,78].
Phase diagram and metastable phases.Finally, we can combine the information from the heating/cooling simulations (Fig. 2 and Fig. 3), the free energy landscape (Fig. 1c), and the zero Kelvin energy distributions (Fig. 4) to construct a free energy diagram that shows the stability ranges of the different phases observed in this work (Fig. 6).
The tetragonal-to-cubic transition occurs around 309 K according to the cooling simulations.While it has first-order character the latent heat is very small, whence it can be readily sampled in MD simulations (Fig. 3).This is consistent with the free energy landscape in Glazer space, which shows the distinction between the tetragonal a 0 a 0 c − and the cubic a 0 a 0 a 0 phase to vanish slightly above 300 K ( Fig. 1c).
The transition from the orthorhombic a − a − c + to the tetragonal a 0 a 0 c − phase involves switching from in-phase to out-of-phase tilts relative to the c-axis and is unsurprisingly associated with a free energy barrier (Fig. 1c).This prevents us from sampling this transition directly with MD simulations.Based on the free energy differences we can, however, locate this transition at approximately 90 K.As a result of this free energy barrier, the orthorhombic a − a − c + remains metastable up to about Free energy difference (meV/f.u.) Phase diagram of MAPbI3 constructed based on MD simulations, free energy calculations, and energy minimization carried out in this work.Symbols represent local minima in the free energy landscape obtained from umbrella sampling (Fig. 1).Lines are drawn as guide for the eyes.Free energy differences are given relative to the orthorhombic phase at each temperatures (suitably extrapolated at high temperatures).Grey horizontal lines indicate phase transition temperatures obtained from the free energy differences for the orthorhombic-to-tetragonal (a − a − c + to a 0 a 0 c − ) and from cooling MD simulations for the tetragonal-to-cubic (a 0 a 0 c − to a 0 a 0 a 0 ) transition.The dashed line indicates the transition from the orthorhombic a − a − c + to the tetragonal a − a − c 0 phase as observed in heating MD simulations (Fig. 2).150 to 150 K (Fig. 2d).Above this temperature it becomes unstable with respect to the alternative tetragonal a − a − c 0 phase, which itself is metastable with respect to the tetragonal a 0 a 0 c − phase that is also the one observed experimentally.The metastable a − a − c 0 phase (alternatively written as a 0 b − b − ) is identifiable as a local free energy minimum up to about 250 K beyond which point the free energy differences become numerically too small to be resolvable.
While the present modeling does not yield the same transition temperature as experiments (caused by the approximations implicit to any current exchangecorrelation functional [46]), we expect our results to be semi-quantitatively correct.Hence, it is suggested that the metastability of the orthorhombic a − a − c + is a genuine feature of MAPbI 3 , and the result of the orientational degrees of freedom associated with the molecular MA cation.This might explain some of the uncertainty in the experimental analyses of this material.Moreover it indicates the local structural variability that can be present as well as the complexity of the underlying dynamics.

IV. CONCLUSIONS
For CsPbBr 3 our simulations provide strong support for the cubic-to-tetragonal phase transition to have first order character, in agreement with the majority of experimental studies.The latent heat is only on the order of a few 0.1 meV atom −1 , which allows one to observe the transition directly by MD simulations, even in rather large systems.We find the orthorhombic-totetragonal transition in CsPbBr 3 to be completely continuous, which overall gives rise to a rather simple free energy landscape and phase diagram.
As a result of the additional orientational degrees of freedom associated with the molecular MA cation, the situation is more complex for MAPbI 3 .While the cubicto-tetragonal transition has first-order character with a similarly small latent heat as in the case of CsPbBr 3 , the low temperature of the phase diagram is notably more intricate.The orthorhombic-to-tetragonal transition has first-order character and is associated with a rather large free energy barrier that cannot be readily overcome on the length and time scales of MD simulations.In our simulations the orthorhombic a − a − c + phase is found to be metastable up to about 40 K above the orthorhombic-to-tetragonal phase transition.We also identified a metastable tetragonal a − a − c 0 phase (alternatively written as a 0 b − b − ).
As noted before, the increased complexity of MAPbI 3 compared to CsPbBr 3 can be attributed to orientational degrees of freedom associated with the molecular cation.This leads to an energy (and by extension free energy) landscape with many, usually rather shallow, local minima.Most notably it causes the preferred tetragonal structure to switch from in-phase (a 0 a 0 c + ), which is preferred by many inorganic halide perovskites including CsPbBr 3 , to out-of-phase (a 0 a 0 c − ) tilting.This suggests that the trivial transition path from the tetragonal to the orthorhombic phase involves switching from out-of-phase to in-phase tilting, suggesting an intuitive understanding for the larger free energy barrier of this transition compared to transitions in inorganic halide perovskites such as CsPbBr 3 .The increased complexity of the (free) energy landscape moreover gives rise to metastable phases, which can in principle appear (possibly transiently) in experimental samples as well.This might explain some of the uncertainty in the experimental analyses of these materials [13].
Orientational degrees of freedom are also present in other hybrid (organic-inorganic) perovskites, most notably FAPbI 3 , for which one can expect similar features.It will therefore be particularly interesting to discriminate these effects in mixed systems, involving but not limited to MAPbI 3 , FAPbI 3 , and CsPbI 3 .We hope that the present work provides a strong basis for future explorations in these directions.
FIG. 1.(a) Commonly occurring perovskite phases and their positions in Glazer space.The insets on the left-hand side show the tilting patterns along the x, y, and z directions, where red and blue indicate positive and negative tilt angles, respectively, whereas gray represents the absence of tilting.In addition to the four experimentally observed phases shown on the left, the figure also indicates the position of the a − a − c 0 structure (in purple), which can be observed on heating in MAPbI3 if the first-order transition from a − a − c + to a 0 a 0 c − is kinetically suppressed (see Fig. 2c).The atomic structures were visualized using ovito [58].(b,c) Free energy landscape of (b) CsPbBr3 and (c) MAPbI3 in Glazer space a • b • c • with a = b.The abscissa in the free energy heat maps corresponds to the tilt angle in degrees about both x and y-axes, where negative and positive values indicate out-of-phase and in-phase tilting, respectively.The ordinate is the tilt angle about the z-axis.

FIG. 3 .
FIG. 3. (a) Energy and (b) lattice parameters from MD simulations in the N P T ensemble around the cubic to tetragonal phase transition for CsPbBr3.Here, mean denotes the average over all MD data while mode refers to the point at which the probability is maximized.The solid lines correspond to the mean µ of the fits to (S1).The dotted line indicates the transition temperature between the orthorhombic and tetragonal phases.(c,d) Full distribution at 302 K, indicated by the dashed line in a) and b), and solid lines indicate fits to (S1).(e-h) The same data as in (a-d) but for MAPbI3, where the histograms are shown at 309 K.

FIG. 4 .
FIG. 4. Energy distribution of fully relaxed MAPbI3 structures obtained by generating and relaxing 10 000 possible tilted structures with randomized MA orientations in 2×2×2 supercells of the primitive cubic cell.

FIG. 5 .
FIG. 5. (a, b) Autocorrelation function C(τ ) for the orientation of the MA units in each phase.The spacing between the lines is 20 K. (c) Rotation rate 1/τrot obtained from (a-b) as a function of temperature.The solid lines correspond to Arrhenius fits 1/τrot(T ) ∝ e −E A /k B T , where EA is the activation barrier of the rotation process.