Understanding correlations in BaZrO3:Structure and dynamics on the nano-scale

Barium zirconate BaZrO3 is one of few perovskites that is claimed to retain an average cubic structure down to 0K at ambient pressure, while being energetically very close to a tetragonal phase obtained by condensation of a soft phonon mode at the R-point. Previous studies suggest, however, that the local structure of BaZrO3 may change at low temperature forming nanodomains or a glass-like phase. Here, we investigate the global and local structure of BaZrO3 as a function of temperature and pressure via molecular dynamics simulations using a machine-learned potential with near density functional theory (DFT) accuracy. We show that the softening of the octahedral tilt mode at the R-point gives rise to weak diffuse superlattice reflections at low temperatures and ambient pressure, which are also observed experimentally. However, we do not observe any static nanodomains but rather soft dynamic fluctuations of the ZrO6 octahedra with a correlation length of 2 to 3nm over time-scales of about 1ps. This soft dynamic behaviour is the precursor of a phase transition and explains the emergence of weak superlattice peaks in measurements. On the other hand, when increasing the pressure at 300K we find a phase transition from the cubic to the tetragonal phase at around 16GPa, also in agreement with experimental studies.


I. INTRODUCTION
Perovskite oxides constitute a prominent class of materials with a wide range of different properties, such as ferroelectricity, colossal magnetoresistance, electronic and/or ionic conductivity, piezoelectricity, superconductivity, metal-insulator transition, luminescence, and many more [1].
The prototypical oxide perovskite structure is cubic, with the general chemical formula ABO 3 , where the A and B sites can accommodate a wide variety of elements from the periodic table.Many perovskites are cubic at high temperatures but upon cooling most undergo one or several structural phase transitions, which depend sensitively on the choice of A and B. These phase transitions are often related to tilting of the BO 6 octahedra, typically referred to as antiferrodistortive transitions.Commonly they are out-of-phase and in-phase tilting phonon modes related to instabilities at the R and/or M-points of the Brillouin zone.
Barium zirconate BaZrO 3 is rather unique among the oxide perovskites.Neutron powder diffraction studies show that BaZrO 3 at ambient pressure maintains its high temperature cubic structure down to temperatures close to zero Kelvin [2][3][4].While the antiferrodistortive Rtilt mode softens substantially with decreasing temperature, its frequency remains positive as the temperature approaches zero Kelvin [5,6].
While the latter experiments have clearly established the long-range order, the short-range order of the cubic BaZrO 3 phase is more controversial.Raman spectra show pronounced peaks despite that first-order scattering is prohibited by symmetry reasons for cubic systems [7][8][9].This has been interpreted as evidence for distorted nanodomains with lower than cubic symmetry, giving * goran.wahnstrom@chalmers.se rise to first-order broad Raman spectra [7,10].A somewhat similar idea, an "inherent dynamical disorder", has been put forward to account for the apparent local deviation from the cubic structure identified by Raman spectroscopy [9].On the other hand, Raman studies of BaZrO 3 single crystals associated these spectral features to second-order events, but stated that it is likely that the overall scattering intensity finds its origin in some other type of local disorder [11].It has also been argued that a structural "glass state" may be formed upon cooling due to the extremely small energy differences between the phases allowed from condensation of the R mode [12].The structural order could then be distorted on the local scale but appear cubic in diffraction experiments.
Recent electron diffraction experiments by Levin et al. [13] suggest that BaZrO 3 undergoes a local structural change associated with correlated out-of-phase tilting of the ZrO 6 octahedra when the temperature is reduced below 80 K.They found weak, but clear, diffuse scattering intensity at the R-point (3/2, 1/2, 1/2), where the soft mode connecting the cubic to the tetragonal phase is located; yet their average structure remained cubic.The authors suggested that nanometer-sized domains ("nanodomains") with local tetragonal structure could explain the diffraction results.The size of these domains was estimated to be about 2 to 3 nm based on the full width half maximum (FWHM) of the diffraction peaks.They stated that the emergence of these relatively sharp superlattice reflections resembles a phase transition more than dynamic correlations, but their measurements could not conclusively discern between static and dynamic effects.
The pressure dependence of BaZrO 3 at room temperature has been investigated by several authors [7,[14][15][16].In a recent combined X-ray diffraction and Raman spectroscopy study [16], it was found that BaZrO 3 undergoes a single phase transition around 10 GPa from the cubic (P m3m) to the tetragonal (I4/mcm) phase and retains that structure up to 45.1 GPa.No second phase transition to an orthorhombic or any other tilted phase was arXiv:2310.05565v1[cond-mat.mtrl-sci]9 Oct 2023 observed.This confirms a previous high pressure X-ray diffraction study from 0 to 46.4 GPa [14], where also a single transition from the cubic to the tetragonal phase was obtained, but at the considerably higher pressure of 17.2 GPa.However, a recent study based on Raman spectroscopy [15] found two structural phase transitions, the first from the cubic to a rhombohedral (R3c) phase at 8.4 GPa and the second from the rhombohedral to the tetragonal phase at 11 GPa.
Here, we construct a machine-learned potential using the neuroevolution potential (NEP) approach trained with density functional theory (DFT) data to be able to simulate the system over long time-scales (100 ns) using large systems (15 million atoms) with near DFT accuracy.The phase diagram for BaZrO 3 is mapped out as a function of temperature and pressure and compared with experiments.The static and dynamic structure factors as a function of wavevector and frequency are computed and their dependence on temperature and pressure are investigated.Detailed and direct comparison is made with the electron diffraction data by Levin et al. [13] and the dynamics close to the R-point is clarified.Finally, the spatial and temporal correlations of the local tilt angles for each individual ZrO 6 octahedron are computed to elucidate the three-dimensional structure and dynamics of BaZrO 3 as a function of temperature and pressure.DFT calculations based on the CX functional yield a lattice parameter for cubic BaZrO 3 of 4.20 Å, for which the phonon dispersion curves show only a very weak instability at the R-point [4].When decreasing the lattice parameter the instability at the R-point increases and for 4.00 Å the dispersion curves also show an instability at the M-point (Fig. S5).

A. Instabilities and phase diagram
In Fig. 1 we show the static energy landscape along the R-tilt mode as function of the oxygen atom displacement.For 4.00 Å a clear double well energy landscape is obtained with depths equal to −29.8 meV/f.u. and located at ±0.25 Å.This corresponds to a tilt angle of 7.1 • .We note that the M-mode instability for 4.00 Å is barely visible on the same energy scale (Fig. S6).
Next, we consider the system at finite temperatures and pressures.MD simulations are carried out in the NPT ensemble where the length of the cell vectors are allowed to fluctuate but the angles between them are kept fixed at 90 • .The system is cooled at constant pressure from high temperature at a rate of 40 K/ns, which is sufficiently slow to avoid any noticeable hysteresis.We also note that it is due to the second-order nature of the phase transition that we can sample and observe it directly in MD simulations.
To monitor the dynamic evolution of the system we use the temperature dependence of the lattice parameters a i and the phonon mode coordinates Q λ .The latter are obtained by phonon mode projections [17,18].The atomic displacements at each time are scaled back to the original cubic supercell and these scaled displacements u(t) are then projected on a tilt phonon mode λ according to where e λ is the supercell eigenvector for mode λ.The mode eigenvectors are obtained using phonopy [19] and symmetrized such that each of the three degenerate modes corresponds to tilting around x, y, and z directions, respectively.A gliding time average with width 20 ps is applied along the trajectory of the cooling simulation allowing us to extract the lattice parameters as well as the phonon mode coordinates Q λ as a practically continuous function of temperature.In Fig. 2a the temperature dependence of the lattice parameters and R-mode coordinates are shown at 10 GPa when the system is cooled from 300 K.At around 160 K the lattice parameter in the x direction deviates from the other two forming a tetragonal structure at the same time as the R x mode condensates.This indicates a phase transition from the cubic a 0 a 0 a 0 (P m 3m) to the tetragonal a 0 a 0 c − (I4/mcm) phase.
Cooling runs are repeated for various pressures and the resulting phase diagram is determined and shown in Fig. 2b.At 300 K we find a phase transition to the tetragonal phase at 16.2 GPa.We do not see any condensation of the M-tilt modes (in-phase tilting) at any pressure or temperature.Furthermore, the phase transition only occurs to the tetragonal phase (I4/mcm), not to any orthorhombic or rhombohedral phases, except for a small region below 20 K and around 4 to 5 GPa, where the rhombohedral (R3c) structure becomes stable.However, for these low temperatures quantum fluctuations become important and we expect these to stabilize the tetragonal structure as discussed below.The observed lattice parameters as a function of temperature and pressure is shown in Fig. S8, and agrees well with experimental work [20].
Below about 100 K quantum effects have been shown to be important to correctly model the stability of the cubic phase [6].Therefore, the phase diagram obtained here using classical MD simulations becomes less accurate at low temperatures.This is indicated in Fig. 2 by the increased transparency of the color at low temperatures.We note here that while the classical MD simulations predict that the system becomes tetragonal at zero temperature and pressure, it is likely not the case if quantum fluctuations are included (see Fig. S9 and Ref. 6).
The phase transition to the tetragonal phase as function of pressure has been investigated experimentally by Raman spectroscopy [7,15,16] and XRD measurements [14,16].The experimental results are rather scattered.In Raman studies the phase transition to the tetragonal structure was observed at 11 GPa [7], 10 GPa [16], and 19.2 GPa [15] at room temperature, and in XRD measurements at 17.2 GPa [14] and 12 GPa [16].Our observed phase transition, from the cubic to the tetragonal phase, at 16.2 GPa falls approximately in the middle of experimentally observed range.In Ref. 15 a transition to a rhombohedral (R3c) structure was also obtained at 8.4 GPa.This type of transition was, however, neither observed in the other experimental studies nor does it appear in the present simulations.Dynamical structure factor, S(q, ω), for various temperatures at 0 GPa.Here, q = 2π a (3/2, 1/2, 1/2), corresponding to the R-point in the second Brillouin zone.
Next, we consider the temperature dependence of the structure factor at ambient (zero) pressure.The intermediate scattering function is defined as where r i (t) denotes the position of atom i at time t, N is the number of atoms, and ⟨. ..⟩ indicates a time average.The dynamic structure factor S(q, ω) is obtained by a temporal Fourier transform of F (q, t).
We calculate the intermediate scattering function F (q, t) from MD simulations in the NVE ensemble.The total simulation time for a run is 1 ns and F (q, t) is averaged over 100 independent such simulations.The corresponding dynamic structure factor S(q, ω) is shown in Fig. 3 at the R-point q = 2π a (3/2, 1/2, 1/2).The lower peaks, below 10 meV, correspond to the R tilt-mode and the peaks around 13 meV correspond to the acoustic mode.The R tilt-mode shows a strong temperature dependence, softening with decreasing temperature.This is in good agreement with previous experimental and theoretical modeling [5,6].
Next, we consider the static structure factor S(q), which is related to the intermediate scattering function and the dynamic structure factor via The partial static structure factors S αβ (q) are then evaluated as where α and β denote the atom types (α = Ba, Zr or O), the summation runs over all atoms of the given type and N α is the number of atoms of type α.The static structure factor is calculated from MD simulations in the NVT ensemble at 0 GPa and different temperatures.For each temperature, we average over 40 independent simulations that are each 100 ps long.We consider first the partial static structure factors S αβ at 100 K calculated along the Brillouin zone path q = 2π a (3/2, 1/2, x) with x : 0 → 1, shown in Fig. 4a, corresponding to a path M → R → M also used by Levin et al. [13] (see their Fig. 5).The oxygen-oxygen part gives rise to large intensity at the R-point (x = 1/2), in agreement with the soft oxygen tilt mode at R, as well as a background intensity.The other partial static structure factors only give rise to a very weak background intensity.
The temperature dependence of the static structure factor S(q) is shown in Fig. 4b at 40, 100, 200, and 300 K.For all temperatures there is a peak at the Rpoint (x=1/2).To further understand this, consider the static structure factor for a harmonic system [21] in the classical limit where the sum runs over all phonon modes for the given q-point and F ph λ (q) is the phonon structure factor containing the Debye-Waller factor and mode selection rules FIG. 4. Structure factors calculated from MD simulations at 0 GPa.Here, q = 2π a (3/2, 1/2, x), which corresponds to a horizontal 1D slice of the data given in Fig. 5, i.e., starting from an M-point (x = 0), passing an R-point (x = 1/2), and ending at another M-point (x = 1).a) Partial structure factors S αβ (q) at 100 K. b) Total structure factor S(q) at 40, 100, 200 and 300 K. c) Scattering intensity I(q) at 40, 100, 200 and 300 K.The dashed line for 40 K shows a Gaussian fit.[22].Therefore, we roughly expect the intensity to increase linearly with temperature T and to scale with frequency as 1/ω 2 .The peak height of S(q) is almost constant with temperature whereas the background increases linearly with temperature in accordance with a harmonic system.The constant peak height is due to that the tiltfrequency of the R-mode softens from 9 meV to about 3 meV between 300 K and 40 K, since 300/9 2 ≈ 40/3 2 .Thus, the structure factor at the R-point remains more or less constant with temperature.The clear peak at 40 K is therefore a result of the tilt-frequency of the R-mode softening with temperature.
The present MD simulations are based on classical mechanics.Quantum fluctuations of the atomic motion start to become important for the R-mode frequency below 100 K [6].We have tested the effect quantum fluctuations on the above peak height using a self-consistent phonon approach (see Fig. S10).The peak height at the R-point is slightly reduced by including the quantum effects, but a clear peak at 40 K is still present.
Lastly, to get a one-to-one comparison with the electron beam diffraction experiments carried out by Levin et al. [13], we determine the intensity, I(q), by weighting the partial structure factors with their corresponding electron atomic scattering factors according to Here, f α (q) are the q-dependent electronic scattering factors for the ions Ba 2+ , Zr 4+ , and O 2 -, with numerical data taken from Ref. 23 (see Fig. S11).The scattering factors are roughly proportional to the atomic number, reducing the oxygen contribution significantly.The peak at the R-point for the intensity I(q) is reduced in height (relative to the background) (Fig. 4c) and there is barely any visible peak above 100 K.This is in good agreement with the observation by Levin et al. [13] that a weak and diffuse, yet discrete spot appears at the R-point below about 80 K.For 40 K we also carry out a Gaussian fit of I(q) with a constant background to extract the FWHM of 0.23/ Å, which is in very good agreement with the value of 0.22/ Å, reported by Levin et al. [13].
FIG. 5. Intensity normalized by temperature, I(q)/T , calculated from MD simulations at a) 40 K, b) 100 K, and c) 300 K. Here, q = 2π a (3y, y, x).The lower left corner thus corresponds to the q-point q = (0, 0, 0), the upper left corner to q = 2π a (3, 1, 0), the lower right corner to q = 2π a (0, 0, 1), and the upper right corner to q = 2π a (3, 1, 1).The center of the heatmaps corresponds to the R-point q = 2π a (3/2, 1/2, 1/2).Note also that the color scale is set such that diffuse scattering is visible, but in practice the intensity at the Γ-points are orders of magnitude larger.
Next, we extend the calculation of the intensity I(q) to the same 2D space of q-points as highlighted by Levin et al. [13] in their Figure 4.Because the intensity increases almost linearly with temperature (Eq.( 5)), we plot I(q)/T to enable easier comparison between temperatures.These normalized intensities are shown as heatmaps in Fig. 5. Heatmaps for the partial intensities I αβ (q) at 100 K can be found in Fig. S12.Most of the intensity heatmaps in Fig. 5 look very similar for all three temperatures.The larger intensities in the corners (Γ points) corresponds to the Bragg peaks.The intensity between Bragg peaks, the diffuse scattering, arises due to thermal motion.The only real notable difference between the temperatures is the increased intensity in the middle of the heatmap (at the R-point q = 2π a (3/2, 1/2, 1/2)) for lower temperatures.At 300 K there is almost no peak visible at the R-point compared to the intensity level for the surrounding q-points, whereas for 40 K there is a very clear peak (as also seen in Fig. 4c).The notable intensities at x=1, y=1/3 and x=1, y=2/3 arise from the low frequency Ba and Zr modes between Γ and M (and close to X) in the phonon dispersion (cf .Fig. S13 and Fig. S12).

C. Tilt angle correlations: Temperature dependence
To obtain a more local picture we now consider the tilt angle of each individual ZrO 6 octahedron, and its static and dynamic correlations.We first extract the Euler angles for each octahedron from MD simulations.We employ the polyhedral template matching using ovito [24] as done in Ref. 25.This allows us to extract tilt angles around the α axis (α = x, y, z) for an octahedron located at (n x , n y , n z ) at time t, θ α (n x , n y , n z , t).Here, we follow a similar notation as in Refs.26, 27.
The distribution P (θ) over θ α (n x , n y , n z , t) averaged over α, (n x , n y , n z ), and t can now be determined (Fig. 6a).We notice that for all temperatures the distribution exhibits a Gaussian profile with zero mean and with a standard deviation that increases with temperature.The standard deviations are σ = 0.82 • , σ = 1.21 • , σ = 1.62 • , and σ = 1.90 • for 40 K, 100 K, 200 K, and 300 K, respectively.In a classical harmonic system we expect the variance σ 2 to increase linearly with temperature but here, due to the softening of the R-tilt mode, the distribution over tilt angles shows a weaker temperature dependence.
Next, we consider the static tilt-angle correlation function between θ α (n x , n y , n z , t) and its neighboring octahedra.Here, we only consider neighbors along the [100], [010] and [001] directions.The static correlation function in the [100] direction is calculated as where d corresponds to the number of neighbor distances between two octahedra in the x-direction and ⟨. ..⟩ denotes an average carried out over (n x , n y , n z ) and t.Sim- The result for the static tilt-angle correlation functions are shown in Fig. 6b.For both G ⊥ (d) and G ∥ (d) the correlation alternates between positive and negative values when increasing the neighbor distance, since the Rtilt mode dominates that motion [13,25,26].We thus only show |G ⊥ (d)| and |G ∥ (d)| in Fig. 6b.The alternation of the correlation is demonstrated in Fig. S14, where the joint probability distribution over two angles is shown.For the correlation perpendicular to the rotational axis, |G ⊥ (d)|, we find a strong correlation between nearest neighbor octahedra which decays towards zero after about 4 to 5 neighbor distances, corresponding to about 2 nm (Fig. 6b).In the direction parallel to the rotation axis, the correlation is weaker and decays faster.This is related to the soft phonon mode at the M-point corresponding to in-phase tilting of the octahedra that thus to some extent counteracts the out-of-phase tilting by the R-mode.For both G ⊥ (d) and G ∥ (d) we see that the correlation increases with decreasing temperature.This is connected to the softening of the R-mode frequency which causes the correlation length to increase.
Lastly, we consider the dynamic autocorrelation function for the tilt angles θ α (n x , n y , n z , t) defined as where ⟨. ..⟩ corresponds to an averaged carried out over (n x , n y , n z ) and t.The result for the correlation function, averaged over α, is shown in Fig. 6c.For all four temperatures the correlation function oscillates at short-time scales and then goes to zero after a few picoseconds.This is clear indication that there are no static or "frozen in" tilts for the temperatures considered, but rather the tilts are dynamically changing on a picosecond time-scale.

D. The structure factor -pressure dependency
Let us now consider the pressure dependence of the dynamic and static structure factors at 300 K.The system is studied from 0 GPa to 18 GPa and at 16.2 GPa it transforms from the cubic to the tetragonal phase.
The pressure dependence for the dynamic structure factor is shown in Fig. 7.For 0 GPa we see a clear peak at around 9 meV corresponding to the R-tilt mode.The peaks above 12 meV correspond to the acoustic mode.The frequency of the R-tilt mode decreases with increasing pressure and the magnitude of the dynamic structure factor increases substantially (Notice the logarithmic scale on the y-axis.)At the same time the damping of the mode increases and at around 15 GPa it becomes overdamped.
The corresponding static structure factor is shown in Fig. 8.The static structure factor has the shape of a Lorentzian peak on a log-scale.The width of peak decreases when approaching the phase transition, indicating that the correlation length increases.Further, the value of static structure factor increases exponentially at the R-point as one approaches the phase transition pressure.This can be understood from the fact that the frequency of the R-tilt mode approaches zero at the phase transition and thus S(q) diverges (cf.Eq. 5).Furthermore, the large values of S(q, ω) observed for higher pressures at low frequencies is directly related to the di- Frequency (meV) S(q, ω) The dynamical structure factor, S(q, ω), calculated from MD for various pressures at 300 K. Here, q = 2π a (3/2, 1/2, 1/2).
vergence of the static structure factor S(q) (Fig. 8), as can be seen from Eq. 3.

E. Tilt angle correlations -pressure dependency
Finally, we consider the tilt angle and its static and dynamic correlations as function of pressure, shown in Fig. 9.The highest pressure, 18 GPa, is located above 16.2GPa, the pressure where the system transform from the cubic to the tetragonal phase.The data for 18 GPa is therefore calculated using only the direction α for which the tetragonal structure is tilted around.For the three lower pressures the data are obtained by making an average of the three different α directions.
The distribution for the tilt angle P (θ) as function of pressure is shown in Fig. 9a.The distribution widens with increasing pressure.For 18 GPa, where the system has undergone the phase transition to the tetragonal phase, the distribution develops a symmetric double peak distribution.This can be fitted well with two Gaussians with mean values µ t = ±2.65 • and standard deviation σ t = 2.19 • .
The static tilt-angle correlation function as a function of neighbor distance is shown in Fig. 9b.The static correlations increase as function of pressure and the decay distance increases.Above the phase transition the correlations do not decay to zero and the correlation function approaches the constant value for 18 GPa, reflecting the (global) long-ranged tilting in the tetragonal phase.Similar behavior is also seen in the dynamic tilt-angle autocorrelation function C α (τ ) in Fig. 9c.For pressures below the phase transitions C α (τ ) decays to zero in the long-time limit, whereas for 18 GPa C α (τ ) approaches the same constant value as the static correlation function, i.e., C α (τ → ∞) = 0.59.It is interesting to note that just below the phase transition the decay time increases substantially.This has also been seen in similar simulation studies of halide perovskites [18,26].

III. DISCUSSION
Structural instabilities and phase transitions in perovskite oxides are important and have therefore been investigated extensively.Strontium titanate SrTiO 3 (STO) is generally considered to be a model perovskite for the study of soft mode-driven phase transitions [28,29] and it may be instructive to compare the behavior of STO with BZO.
At ambient conditions STO is cubic and its antiferrodistortive transition to the tetragonal (I4/mcm) phase can be induced by either decreasing the temperature or increasing the pressure [30].The pressure induced transition at room temperature occurs at 9.6 GPa for STO [31].The same type of transition also occurs in BZO but at a somewhat higher pressure [14,16].On the other hand, the temperature induced transition at ambient pressure only occurs in STO, not in BZO.In STO the R-tilt mode softens and at about 105 K [28,32] it approaches zero and the material undergoes a phase transition to the tetragonal structure.When approaching this phase transition from above the scattering intensity near the R-point increases dramatically and the scattering peak narrows substantially in q-space [32].Our results demonstrate that a similar mechanism is also at play in BZO and detected in the experiments by Levin et al. [13].Yet in contrast to STO, one only reaches the initial narrowing of the peak as the phase transition never occurs at ambient pressure.The R-tilt mode softens but its frequency remains finite when the temperature goes to zero [5,6].Levin et al. [13] find a diffuse peak at the R-point with a width of 0.22/ Å.This magnitude corresponds roughly to the corresponding peak for STO at about 160 K, that is 50 K above the transition to the tetragonal phase [32].The "nanodomains" observed by Levin et al. [13] are thus dynamic correlations at the onset of a phase transition that never occurs in BZO at ambient pressures.

IV. CONCLUSIONS
We have performed large scale MD simulations of barium zirconate, an oxide perovskite, using machinelearned potentials based on DFT calculations.Both the temperature and pressure dependence of the local and global structure and the dynamics were investigated, and compared with available electron diffraction results.
At ambient pressure it is now well established that BZO remains cubic down to zero Kelvin, although the R-tilt mode softens substantially [6].Our MD simulations predict a softening from 9 meV at 300 K to 3 meV at 40 K.We find that this mode softening gives rise to a clear oxygen related peak in the static structure factor at the R-point, which explains the superlattice reflection observed by Levin et al. [13] using electron diffraction.
Levin et al. [13] state that the peak is only visible be-low about 80 K.However, we show that it does exist also at higher temperatures, albeit with weaker intensity.The present study strongly suggests that the disappearance of the peak in the electron diffraction study at higher temperatures is due to the large background intensity from scattering of Ba and Zr at those temperatures.The oxygen related peak is the result of strongly correlated and dynamic tilting between neighboring ZrO 6 octahedra.By investigating the tilt angle correlations we find that the spatial extent of the correlated motion at 40 K is about 2 to 3 nm and with a short relaxation time of about 1 ps.We therefore conclude that the oxygen peak observed at the R-point is purely of dynamic origin.
The pressure dependence at room temperature was also investigated.It is known that BZO undergoes a phase transition from the cubic to the tetragonal phase.Here, we obtain this transition at about 16 GPa in the middle of the experimentally observed range.When approaching the phase transition from lower pressures, the frequency of the R-tilt mode approaches zero and close to the phase transition the motion becomes overdamped.At the same time the static structure factor at the Rpoint increases dramatically by several orders of magnitude.The dynamic tilt-angle autocorrelation function shows a rapid decay on the order of 1 ps, but close to the transition, the correlation function also develops a component with a considerably slower decay.At the phase transition this decay goes over to a constant finite value.The static tilt-angle correlation function shows a similar behavior: The decay rate becomes longer and longer and the correlation function approaches a constant value at the phase transition.
The present study shows that large scale MD simu-lations based on a machine-learned potentials with near DFT accuracy can provide immensely detailed and accurate atomic scale information on the local structure and complex dynamics close to phase transitions.
V. METHODS

A. Reference calculations
The energy, forces, and virials are obtained for the training structures via DFT calculations as implemented in the Vienna ab-initio simulation package [33][34][35] using the projector-augmented wave [36,37] setups in version 5.4.4 with a plane wave energy cutoff of 510 eV.The considered valence configurations for Ba, Zr and O are 5s 2 5p 6 6s 2 , 4s 2 4p 6 4d 2 5s 2 and 2s 2 2p 2 , respectively.The Brillouin zone is sampled with a Monkhorst-Pack grid, with the maximum distance between two points being 0.19/ Å along the reciprocal lattice vectors.This leads to a 8 × 8 × 8 k-point grid for the primitive cell with a lattice parameter of 4.20 Å.
For the exchange-correlation functional we employ the van-der-Waals-density functional with consistent exchange (vdW-DF-cx) [38,39], here abbreviated CX.This functional is a version of the vdW-DF method [40], with the aim of accurately capturing competing interactions in soft and hard materials [41,42].It has been applied to BaZrO 3 before [4,43] and been found to give a very good account of its structural and vibrational properties.In particular, the anharmonicity of the R-tilt mode at ambient pressure is well described compared to recent experiments on BaZrO 3 [6], and so is the thermal expansion [43] (for more details see Fig. S8 and Fig. S7).

B. Neuroevolution potential
We construct an NEP model for the potential energy surface using the iterative strategy outlined in Ref. 44.Training structures include cubic, tetragonal and rhombohedral primitive cells at different volumes and cellshapes, MD structures in a 4 × 4 × 4 (320 atoms) supercell at temperatures up to 500 K and pressures up to 40 GPa, structures with various tilt-modes imposed, cubic-tetragonal and tetragonal-tetragonal interface structures, and structures found by simulated annealing at different pressures.The MD structures are generated using an initial NEP model and are selected based on their uncertainty, which is estimated from the predictions of an ensemble of models [44].The final NEP model used in the production runs is trained on all the available training data (see Fig. S1).The NEP model accurately reproduces the energy-volume curves for the different phases, the phonon dispersions for the cubic phase as well as the static energy landscape of the tilt modes (R and M).More details pertaining to the validation of the NEP model including parity plots are provided in the Supporting Information.

C. Molecular dynamics
All MD simulations are run with gpumd [45].In all simulations we employ a timestep of 1 fs and equilibration time of 50 ps.For most simulations we use supercells comprising 24 × 24 × 24 cubic primitive cells (∼ 70 thousand atoms).However, the static structure factor S(q) is calculated from MD simulations with 144 × 144 × 144 cubic primitive cells (∼ 15 million atoms) in order to achieve an adequate q-point resolution.The static and dynamic structure factors are calculated from MD trajectories using the dynasor package [46].
The phase diagram is obtained from simulations in the NPT ensemble, static properties from simulations in the NVT ensemble, and dynamic properties from simulations in the NVE ensemble.NVT and NVE simulations are carried out with lattice parameters obtained from NPT simulations (see Fig. S8).

FIG. 1 .
FIG.1.The potential energy landscape for the R-tilt mode obtained with the NEP model as function of the oxygen atom displacement.The inset shows the atoms in a tilted structure, where red is oxygen, green zirconium and the blue faces show the ZrO6 octahedra.

FIG. 2
FIG. 2. a) Temperature dependence of lattice parameters and R-mode coordinates at 10 GPa from a cooling MD run in the NPT ensemble.The phase transition from the cubic to the tetragonal phase is observed at about 160 K.The effective lattice parameter for the tetragonal phase is given by V 1/3 , where V is the volume per formula unit.b) Phase diagram from cooling MD runs.Here, the region below 100 K is shown with increasing transparency to reflect the uncertainty due to the classical sampling in MD.Experimental Raman spectroscopy data from Refs. 7, 15, 16 and XRD data from Refs.14, 16.

FIG. 6 .
FIG. 6. Distributions of the tilt angle θ α and its correlations from MD simulations at 0 GPa and different temperatures.a) Tilt angle distribution, averaged over α = x, y, z.Solid lines represent Gaussian fits with zero mean.b,c) Static tilt-angle correlation functions G ∥ (d) and G ⊥ (d) as a function of neighbor distance d.Solid lines are guides to the eye.d) Dynamic tilt-angle correlation function C(τ ) as defined in Eq. (7) and averaged over α = x, y, z.

FIG. 9 .
FIG.9.Properties of the tilt-angle θ α and its static and dynamic correlations from MD simulations at 300 K and for four different pressures.The tilt-angle distribution and its correlations are averaged over α = x, y, z for 0, 8, and 16 GPa.For 18 GPa they are only calculated over the direction α for which the tetragonal structure is tilted around.a) Distribution over observed tilt-angles.Solid lines corresponds to Gaussian fits with zero mean, or symmetric Gaussians with nonzero mean for 18 GPa.The dashed lines correspond to the two symmetric Gaussians for 18 GPa.b) and c) The static tilt-angle correlation functions, G ∥ (d) and G ⊥ (d), as a function of neighbor distance d.Solid lines are guides to the eye.c) The dynamic tilt-angle correlation function C(τ ) as defined in Eq. (7).