Limits of the phonon quasi-particle picture at the cubic-to-tetragonal phase transition in halide perovskites

The soft modes associated with continuous-order phase transitions are associated with strong anharmonicity. This leads to the overdamped limit where the phonon quasi-particle picture can breakdown. However, this limit is commonly restricted to a narrow temperature range, making it difficult to observe its signature feature, namely the breakdown of the inverse relationship between the relaxation time and damping. Here we present a physically intuitive picture based on the relaxation times of the mode coordinate and its conjugate momentum, which at the instability approach infinity and the inverse damping factor, respectively. We demonstrate this behavior for the cubic-to-tetragonal phase transition of the inorganic halide perovskite CsPbBr$_3$ via molecular dynamics, and show that the overdamped region extends almost 200 K above the transition temperature. Further, we investigate how the dynamics of these soft phonon modes change when crossing the phase transition.


INTRODUCTION
The vibrational properties of solids are pivotal for many physical phenomena, including but not limited to phase stability and thermal conduction. In crystalline solids, the vibrational spectrum is commonly described in terms of phonons as quasi-particle representations of the lattice vibrations. The phonon frequency ω 0 is typically much larger than the damping Γ, and the phonon relaxation time τ = 2/Γ is thus much longer than the oscillation period, such that the quasi-particle picture is well motivated [1][2][3][4][5][6]. In this so-called underdamped limit, the relaxation time decreases as the damping Γ increases.
By comparison, there are far fewer cases when phonon modes become overdamped, i.e., ω 0 τ < 1 [7,8]. This can occur either due to large damping or for very soft modes, usually in the immediate vicinity of a phase transition, as for example in the case of body-centered cubic Ti [9][10][11], rotationally disordered 2D materials [12], in ferroelectrics such as BaTiO 3 [13][14][15][16][17] or in halide perovskites [18,19]. In the overdamped limit, the relaxation time increases with increasing damping Γ, which calls into question the picture of a well-defined phonon mode with a frequency and relaxation time. Overdamped phonon dynamics is, however, usually limited to a rather narrow temperature window and under these circumstances the inversion of the relation between relaxation time and damping cannot be readily observed. Here, we demonstrate that the soft phonons modes associated with the phase transitions in the prototypical halide perovskite CsPbBr 3 are, however, outstanding manifestations of this exact behavior as the overdamped region extends almost 200 K above the tetragonal-cubic phase transition.
Halide perovskites are promising materials for photovoltaic and optoelectronic applications. Specifically, CsPbBr 3 has received a lot of attention in recent years [20]. With increasing temperature it undergoes phase transitions from an orthorhombic (Pnma) to a tetragonal (P4/mbm) and eventually a cubic phase (Pm3m) [21][22][23][24][25]. These phase transitions are connected to specific phonon modes and arise due to tilting of the PbBr 6 octahedra, corresponding to phonon modes at the R and M points (Fig. 1a) [26][27][28][29][30]. Experimentally, these modes have been shown to exhibit overdamped characteristics in the vicinity of the phase transitions [18,19,31]. The phase transitions have also been studied from first-principles and via molecular dynamics (MD) simulations, see, e.g., Refs. 32-35. Here, we reveal the dynamics of the octahedral tilt modes in CsPbBr 3 over a wide temperature range via MD simulations based on a machine-learned potential (MLP) that achieves close to density functional theory (DFT) accuracy (Supplementary note S2, Fig. S1) [36,37]. To obtain access to mode specific dynamics we project the MD trajectories onto different normal modes that are associated with the phase transitions in this material. As shown below, this requires both large large systems (comprising at least several 10 000 atoms) and sufficiently long times scales (∼ 50 ns to 100 ns) in order to achieve converged results (see Supplementary note S7, Fig. S9 and Fig. S8). The DFT data and the MLP models are provided as a Zenodo dataset [38].

RESULTS AND DISCUSSION
Tilt-modes and phase transitions. The phase transitions in CsPbBr 3 , and similarly in many other perovskites, are driven by modes that correspond to the tilting of the PbBr 6 octahedra. These modes are located at the M (in-phase tilting) and R-points (out-of-phase tilting) in the phonon dispersion for the cubic structure ( Fig. 1a). They are three-fold degenerate corresponding to tilting around the three Cartesian directions. These tilt-modes exhibit a double-well potential energy surface (PES), which the MLP reproduces perfectly compared to DFT (Fig. 1b).
The MLP predicts temperatures of 300 K and 275 K for the cubic-tetragonal, T c↔t , and tetragonal-orthorhombic, T t↔o , transitions respectively (Fig. 1c). This is lower than the experimental values of 400 K and 360 K [20,22,23,25], a discrepancy that can be primarily attributed to the underlying exchange-correlation functional [54].
The mode coordinates of the tilt modes are useful order parameters for analyzing the phase transitions (Fig. 1d,e). At 300 K the system transitions from the cubic to the tetragonal phase as seen in both the lattice parameters and in the freezing in of one of the three M-tilt modes (M z ). For the tetragonal phase two Rmodes (R x and R y ) start to show larger fluctuations, and at 265 K the system transitions to the orthorhombic phase. Here, we also note the slight difference in character between these two phase transitions. For the cubictetragonal transition the order parameter (Q M ) and lattice parameter change sharply at the transition temperature T c↔t (closer in character to a first-order transition), whereas for the tetragonal-orthorhombic transition the order parameter and lattice parameter change more gradually around T t↔o (exhibiting a continuous character) in agreement with experimental observations of the transition character [20,21]. We note here that the mode coordinate is a global order parameter for the system. In the cubic phase even though the mode coordinate is on average zero, there still exists a strong local correlation between the neighboring octahedra. This connects to previous work on perovskites regarding the local atomic structure deviating from the cubic structure while globally still appearing cubic [28,29,[55][56][57][58][59].
Mode coordinate dynamics. The mode coordinates exhibit interesting dynamical behavior already in the cubic phase far above the transition to the tetragonal phase, which can be conveniently observed in the time domain, Fig. 2a, b). At 500 K regular (phonon) oscillator behavior is observed, whereas for 350 K (closer but still above T C ) a slower dynamic component becomes evident. Finally, at 280 K and thus below the phase transition, one ob-serves the common oscillatory motion superimposed on a long timescale hopping motion between the two minima, corresponding to the (degenerate) tetragonal phase (Fig. 1b). We note here that the hopping frequency depends strongly on system size, and is thus not a good observable on its own.
The mode coordinate can be analyzed by fitting the respective autocorrelation functions (ACFs) to a damped harmonic oscillator (DHO) model (Fig. 3). The ACF for a regular (underdamped) mode shows a clear oscillatory pattern as illustrated here by the highest optical mode at the R-point with a typical relaxation time of about 0.37 ps, which is longer than the mode period of about 0.2 ps (Fig. 3a). The M-tilt mode at 500 K has a similar damping but is much softer (yet still underdamped), and the ACF decays with a relaxation time of about 0.58 ps (Fig. 3b). At 350 K (Fig. 3c), however, the same mode is overdamped and in this case the DHO model becomes the sum of two exponential decays, see Eq. 3, with relaxation times τ L =5.22 ps and τ S =0.31 ps. It is interesting to note that the decay time of the ACF at 350 K is about ten-times longer than at 500 K. The DHO fits still match the data very well for both the underdamped and overdamped cases (see Fig. S11 for how the two exponential decays behave for Q and P in the overdamped case).
When Γ/ω 0 increases and the system becomes overdamped the dynamics of the modes is moving towards the diffusive Brownian motion regime. For overdamped modes the relaxation time of the ACF increases as Γ/ω 0 increases, opposite to the underdamped behavior. While this is a well known feature of a simple one-dimensional DHO, here one observes this behavior for phonon modes in a complex atomistic system. This phenomenon arises due to the free energy landscape being very flat close to the transition, resembling a bathtub. As a result of the high friction and weak restoring force, it therefore takes a long time for the DHO to move back and forth around zero ( Fig. 2c; also see Fig. S12 for the power spectra) [60].
Frequency and relaxation times vs temperature. The frequencies and relaxation times of the M-tilt and R-tilt modes are summarized as a function of temperature in Fig. 4. The frequency ω 0 softens significantly with decreasing temperature for both modes, whereas the relaxation time τ is more or less constant. The softening of the frequency thus drives the modes to the overdamped limit with decreasing temperature. The M-tilt and Rtilt modes only become underdamped above 480 K and 410 K, respectively, well above the transition temperature to the tetragonal phase at 300 K. This indicates that we expect the phonon quasi-particle for these modes to work better at high temperatures, which interestingly is the opposite behavior compared to most phonon modes which become more damped and anharmonic with increasing temperature. At the cross-over from the underdamped to the overdamped regime, the two time scales τ S and τ L emerge. When approaching T c↔t we see that τ L increases exponentially, whereas τ S → τ /2.
Self-consistent phonons and effective harmonic models. Next, we analyze the representation of these strongly anharmonic modes by commonly used phonon renormalization techniques, specifically different SCP schemes and EHMs (Fig. 5) (see Supplementary note S5 and Supplementary note S6 for a more detailed description of the methods). We construct 4:th order force constant potentials (FCPs) at each temperature which is used as input to all SCPs methods, see Supplementary note S4 and Fig. S6 for more details. There are several SCPs variants [45]. In SCP-alamode the Green's function approach is employed as implemented in the alamode package [51]. In stochastic self-consistent harmonic approximation (SSCHA) the harmonic free energy is minimized using gradient methods, as implemented in the sscha package [52]. In SCP-hiphive second order force-constants are obtained by iterative fitting to forces from displacements sampled from the harmonic model and forces obtained from the MLP as implemented in the hiphive package [50]. Here, we employ the "bare" SCP implementations in alamode and ss-cha. We note, however, that there are computationally more demanding corrections for both methods [35,61], the analysis of which is, however, beyond the scope of the present work. The EHMs (in this field also referred to as temperature-dependent potentials) are constructed from fitting second-order force constants to displacement and force data obtained from MD simulations with the MLP (see Supplementary note S6 for details).
Here, we find very similar behavior for both M-tilt and R-tilt modes. The three SCP methods (SCP-hiphive, SS-CHA, SCP-alamode) employed here are in great agreement with each other given the differences in theory and implementation between them. The SCP frequencies systematically overestimate the frequency ω 0 obtained from the ACFs by about 1 meV (see Supplementary note S5 for a more detailed description of the SCP methods). The EHMs constructed by fitting the forces from MD trajectories, on the other hand, show good agreement with the mode projection results. We note here that the trend for SCPs and EHMs to over and underestimate frequencies, respectively, appears to hold for all modes in the system, which is in line with previous studies [35,[62][63][64]. However, while EHMs from MD yield a better frequency for Phonon mode coordinates. Mode coordinate Q(t) for the M-tilt mode, a) at 500 K (well above Tc↔t), b) at 350 K (close to Tc↔t), and c) at 280 K (below above Tc↔t). The M-tilt mode is three-fold degenerate (x, y, z) but here only Mz show, and note that for 280 K the system switches the tilt axis at irregular intervals.  the tilt modes compared to SCP, this is not in general true (see Fig. S7 for details).
Behavior near phase transitions. Next, we look at how these modes behave as the system goes through the transition from the cubic to the tetragonal phase. While in the cubic phase the three M (and R) modes (denoted with subscripts x, y, z to indicate the Cartesian direction) are degenerate, the degeneracy is broken in the tetragonal phase and the z-direction becomes symmetrically distinct from the other two (Fig. 1c,d). Therefore, in order to distinguish these modes we will denote them by M xy and M z (analogously for R-modes) in the tetragonal phase. In the tetragonal phase there exist multiple global minima for the M-mode coordinate, as can be seen in Fig. 2c where the system jumps between these minima. To avoid capturing this (system-size dependent) hopping time in the ACFs we employ very large system sizes of up to 400 000 atoms, for which the system remains in the same tetragonal orientation throughout the entire simulation. Furthermore, we extrapolate the frequencies to the infinite system-size limit (Fig. S10).
The resulting frequencies are shown in Fig. 6. For the cubic-to-tetragonal transition the M-frequency does not go to zero at the transition temperature, which is in agreement with the character of the transition being first order, as observed experimentally [20,21,23]. For the tetragonal-to-orthorhombic phase transition, on the other hand, the R xy frequency does go to zero at the transition temperature, in agreement with a continuous transition as observed experimentally [21]. This leads to the long timescale in the DHO trending to infinity, τ L → ∞, as the temperature approaches T C . Additionally, the R xy mode exhibits a strong size dependence close to the transition temperature (see Fig. S10).
The Curie-Weiss law, ω 0 (T ) ∝ (T −T C ) p , provides very good fits for the temperature dependence of the modes driving the phase transitions. For the tetragonal-toorthorhombic transition the fitted critical temperature, T C =273 K, agrees very well with the observed transition temperature of 274(1) K, which is consistent with this transition being a continuous transition [21]. Furthermore, the fitted critical exponent of 0.55 is very close to the value of 1/2 suggested by Landau theory observed in many continuous phase transitions driven by soft modes [15,[65][66][67]. The cubic-to-tetragonal transition has firstorder character as evident from the finite frequency of the M mode at the transition temperature. As a result, fitting both the critical temperature and the critical exponent is ambiguous (due to the absence of data at temperatures for which the frequency goes to zero). We therefore fix the critical exponent to 1/2, which yields a critical temperature of 295 K, about 7 K lower than the transition temperature. Here, the critical temperature corresponds to the temperature at which the cubic phase becomes dynamically unstable, i.e., the point at which the free energy barrier between the two phases disappears. The parameter τ remains fairly constant in the tetragonal phase across its entire temperature range for all four modes (Fig. S13). Interestingly, once the M z mode freezes in (and the tetragonal phase is formed) both the M z and R z modes stiffen significantly with temperature. This results in the R z mode becoming underdamped again with decreasing temperature at around 290 K and both M-modes approaching the underdamped limit as the system approaches the orthorhombic transition.

CONCLUSIONS
We have carried out a detailed computational analysis of the dynamics in CsPbBr 3 , focusing in particular on the tilt modes. We observe overdamped modes for the cubic phase almost 200 K above the cubic-to-tetragonal transition temperature. These overdamped tilt-modes exhibit correlation on very long time scales (τ L ) compared to the typical relaxation time (τ ) or period (1/ω 0 ) of the mode. This is in line with the dynamics of the modes transi-tioning toward Brownian motion due to the frequency approaching zero. What we find here is that these modes can, however, still be mathematically well described as DHOs, which allows one to formally obtain a phonon frequency and relaxation time compliant with a quasiparticle picture. A physically more intuitive description is, however, obtained if the DHO model is described by two relaxation times, which can be approximately associated with mode coordinate and momentum, respectively. As a result of the soft character of these modes, the respective amplitudes can be large already at moderate temperatures. This implies that even for relatively modest electron-phonon coupling strengths, these modes should have a notable impact on the optoelectronic properties of these materials [57,[68][69][70][71]. A systematic investigation of these effects on a per-mode basis would be an interesting topic of further study.
In addition, we demonstrated that commonly used computational phonon renormalization methods agree very well with each other but without extensive correction schemes exhibit systematic errors in describing the frequencies of the anharmonic tilt modes considered here. Understanding the single-point frequencies obtained from such methods and their relation to the full dynamical spectra is thus very important when, e.g., comparing to experimental measurements.

METHODS
To analyze phonon modes directly from MD simulations we employ phonon mode projection [3,72,73]. The MD simulations are carried out using the gpumd package [37,74]. For more details on the MD simulations see Supplementary note S3. The atomic displacements u(t) and velocities v(t) can be projected on a mode λ, with the supercell eigenvector e λ via Q λ (t) = u(t) · e λ and P λ (t) = v(t) · e λ .
Here, the phonon supercell eigenvector of the tilt modes are obtained with phonopy [75], and symmetrized such that each of the three degenerate modes corresponds to tilting around the x, y and z direction respectively. The ACFs of Q and P are calculated in order to analyze the dynamics of the modes of interest as which can be modeled as the ACF of a DHO. The DHO is driven by a stochastic force and has a natural frequency ω 0 and a damping Γ. The ACF of the DHO splits into an underdamped regime (ω 0 > Γ/2) and an overdamped regime (ω 0 < Γ/2). In the underdamped regime the solution of the DHO is C DHO Q (t) = Ae −t/τ cos ω e t + Γ 2ω e sin ω e t , (2) where ω e = ω 2 0 − Γ 2 4 , the relaxation time is τ = 2/Γ, and A is the amplitude [11]. In the overdamped limit the solution becomes the sum of two exponential decays as where τ S,L = τ Here, τ S and τ L denote the short and long timescales, respectively. If the natural frequency approaches zero (e.g., for continuous phase transitions driven by a soft mode) we thus expect τ L → ∞ and τ S → τ /2. In this limit the resulting ACF, C DHO Q (t), would only consist of a single exponential decay, with a decay time approaching infinity, which corresponds to the behavior seen in Brownian motion.
Similar expressions are obtained for the ACF of the phonon velocity, which is C DHO . For the overdamped case it becomes The ACFs for Q and P are fitted simultaneously to the DHO model in order to extract ω 0 and Γ.

DATA AVAILABILITY
The DFT data and the MLP models are provided in a Zenodo dataset [38].