{'fraction_non_alphanumeric': 0.0455435847208619, 'fraction_numerical': 0.021547502448579822, 'mean_word_length': 4.400822803408757, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Revisiting controversies in the prototypical antiferroelectric PbZrO3

The prototypical antiferroelectric PbZrO3 has several unsettled questions, such as the nature of the antiferroelectric transition, possible intermediate phase and the microscopic origin of the Pbam ground state. Using first-principles, we show that no phonon becomes truly soft at the cubic-to-Pbam transition temperature, and the order-disorder character of this transition is clearly demonstrated based on molecular dynamics simulations and potential energy surfaces. The out-of-phase octahedral tilting is an important degree of freedom, which can collaborate with other phonon distortions and form a complex energy landscape with multiple minima. Candidates of the possible intermediate phase are suggested based on the calculated kinetic barriers between energy minima, and the development of a first-principles-based effective Hamiltonian. The use of this latter scheme further reveals that specific bi-linear interactions between local dipoles and octahedral tiltings play a major role in the formation of the Pbam ground state.

Antiferroelectrics (AFEs) form a special class of materials that possess anti-polar displacements, which was first conceptualized by Kittel [1]. AFEs have received much attention in recent years because they hold great potential to reach high energy density for high-power energy storage [2][3][4]. Lead zirconate PbZrO 3 (PZO) is the first discovered AFE material, and crystallizes in a Pbam ground state below ∼505 K [5][6][7]. Besides the anti-polar distortions, this particular structure also exhibits strong long-ranged order tiltings of oxygen octahedra, which is often termed as antiferrodistortive (AFD) distortions.
Fully understanding the AFE ground state of PZO requires an explanation of the origin of the different types of distortions in the Pbam phase. Fthenakis et al. predicted that both the AFE and AFD modes completely soften in the paraelectric cubic phase at the AFE transition temperature (T A ) [8], but from the experimental side no soft phonon has been observed [9]. There is also no consensus regarding what is the most important degree of freedom leading to the ground state. In addition to the soft AFE and AFD modes mentioned in Ref. [8], Tagantsev et al. proposed that the AFE transition in PZO is driven by a non-critical softening of the transverse acoustic mode via flexoelectric coupling [9], whereasÍñiguez et al. believe that the AFD mode is playing the crucial role [10].
Furthermore, the AFE transition in PZO is accompanied by a strong dielectric anomaly at T A [11], which indicates a softening ferroelectric (FE) mode in the hightemperature phase. It means that there should exist a FE phase that is energetically close to the AFE state. This FE phase is likely to be the rhombohedral R3c phase, as it has been reported that the R3c phase has an energy merely a few meV/f.u. higher than Pbam at 0 K, according to ab initio calculations [10,12], and it has been identified as the field-induced FE phase [13]. Therefore, it would be interesting to understand why the ground state of PZO is Pbam rather than R3c (see Refs. [10,14,15] and references therein). The situation is even more complex when considering the possibility of the existence of an intermediate phase in the vicinity of T A [5,11,[16][17][18]. Such intermediate state can only occur over a very narrow temperature range of 3-5 K between the AFE phase and the cubic phase, and its structure has not been fully determined. In this Letter, we aim to address the aforementioned questions by using firstprinciple and first-principle-based techniques.
Phonons in the cubic phase at high temperatures are investigated by combining ab initio molecular dynamics (AIMD) simulations implemented in vasp [19] with the temperature dependent effective potential (TDEP) method [20][21][22]. The potential energy surface (PES) is calculated by density functional theory (DFT) [19], and the kinetic barriers are computed with the climbing image nudged elastic band (CI-NEB) method [23]. In particular, PES allows us to shed light on the existence and possible structure of the intermediate phase. We also develop an effective Hamiltonian for PZO to explain the microscopic origin of its complex ground state and to suggest other and more complex candidates for the intermediate state. Further details of the methods are given in the Supplemental Material (SM) [24].
Let us first check how phonon softens in the cubic phase as temperature decreases. The dispersion for different temperatures above 450 K are shown in Fig. 1a, as calculated by the TDEP method based on ab initio MD simulations with a 4×4×4 supercell (320 atoms). In con-trast to the 0K phonon computed by density functional perturbation theory (DFPT) that exhibits soft branches [25][26][27][28] (also shown in Fig. 1a via dashed lines), all the phonon frequencies are positive for T > 450 K, implying that the cubic structure is dynamically stable at high temperature. We also numerically find via AIMD simulations that, below 450 K, the structure is not cubic [29]. Interestingly, this latter temperature is rather close to the experimental AFE transition temperature T A of ∼505 K.
Moreover, the 0K phonon exhibits the largest instabilities at the Γ, R, and M points of the Brillouin zone, that correspond to q Γ = (2π/a)(0, 0, 0), q R = (2π/a)(1/2, 1/2, 1/2), and q M = (2π/a)(1/2, 1/2, 0), respectively (with a being the 5-atom lattice constant) and that are associated with the FE distortion, outof-phase (oop) tilting, and in-phase (ip) tilting, respectively. Furthermore, according to the calculated dispersion curves, when cooling from high temperature, these phonon modes show noticeable softening, as depicted in Fig. 1b. Note that the calculated zone-center softening agrees reasonably well with the measured data for the Γ-point mode, as evidenced in Fig. 1b [30]. Interestingly, no phonon becomes truly soft (that is none of them drops to zero frequency) when the temperature approaches 450K from above, as consistent with experimental data in the cubic phase [9,30]. In particular, Ref. [9] pointed out that the Σ-and R-point modes are not soft at T A , where q Σ = (2π/a)(1/4, 1/4, 0), which is consistent with our prediction in Fig. 1b, but contrasts with the theoretical work of Ref. [8] that predicted completely softened Σ and R modes.
The lack of soft phonon suggests that the paraelectricto-AFE transition (without counting the intermediate phase) is of first order, as also consistent with experiments [31][32][33]. Moreover, in order to check if this tran-sition has characters of an order-disorder type, as suggested by Ref. [30] based on the measured central mode, we analyze the MD simulation at 600 K, for which the cubic phase is stable. Figure 2a reports the pair correlation function g(r), which is related to the probability of finding an ion at a given distance from another ion. For comparison, we also plot g(r) from cubic stochastic sampling in Fig. 2b, in which the atomic displacements are centered harmonically around the cubic positions, as indicated by the Gaussian distribution of each pair. The only obvious difference between the two sampling schemes lies in the Pb-O pair, which has a single peak with cubic stochastic sampling but splits into two distances with MD sampling. This splitting is due to the finite tilting angle as long-life-time positions in the cubic phase, as evidenced from the MD simulation (see Figs. S2 and S3 in the SM [24] for more information). As the AIMD simulation is considered to be much more realistic than cubic stochastic sampling [34], it indicates that the cubic phase should be regarded as a dynamical average structure, in which thermal excitations cause the system to hop among local minima with finite tilting angles, with low-energy paths that bypass the ideal cubic positions (see Fig. S4 in the SM [24]). Such splitting also exists at even higher temperature, e.g., up to 1000 K, except that g(r) gradually becomes closer to the cubic stochastic case when increasing T . Further understanding of the structures and nature of the transition can be obtained from the PES. With the cubic structure as reference, we plot the double well potential with respect to a single mode distortion (of normalized magnitude denoted as Q), as shown in Fig. 3a. The energy can be lowered with finite FE distortion (Γ − 4 ), ip tilting (a 0 a 0 c + or M + 2 ), or oop tilting (a − a − c 0 or R − 5 ), which are the most significant instabilities in cubic PZO (see DFPT phonon in Fig. 1a). The largest energy gain is equal to 222 meV/f.u. and is provided by the R − 5 mode. Furthermore, the M + 2 or Γ − 4 mode causes less energy change. Such result implies that the oop tilting is likely to be the most important degree of freedom in PZO. This is in line with the work ofÍñiguez et al., which suggests that the same mode plays a crucial role in stabilizing the AFE ground state [10]. It is also worth to mention that all the low-energy structures reported in Ref. [10] (e.g., Pbam, R3c, Pnma, etc.) have the a − a − c 0 tilting as a component, further reflecting its role leading to low-energy states.
Next, we combine a − a − c 0 tilting with other distortions that can yield three competing polymorphs having low energy, viz., the Pbam, R3c, and Pnma phases. Using the crystallographic tool amplimode [35,36], the Pbam structure is mainly composed of three distortions: R − 5 , Σ 2 and S 2 . The R − 5 component captures 58.8% of the total distortion, which is the largest, while it is 37.1% and only 3.8% for Σ 2 and S 2 , respectively. To simplify the picture, we combine the Σ 2 and S 2 modes and compute the 2-dimensional PES by changing R − 5 and Σ 2 +S 2 . The resulting contour plot is shown in the left panel of Fig. 3c, and exhibits four minima that correspond to the Pbam structure [37] -which is a manifestation of collaborative couplings between the R − 5 , Σ 2 and S 2 modes. The energy of these minima relative to the cubic structure (center) is 311 meV/f.u., which is significantly larger than the barrier heights between these minima, which are found to be 89 or 113 meV/f.u. by reversing R − 5 or Σ 2 +S 2 , respectively. Note that using the CI-NEB method allows these barriers to be further lowered via the minimum energy pathways to 51 or 58 meV/f.u., respectively. Considering the energy scale of thermal energy (e.g., k B T = 86 meV at 1000 K), it appears unlikely that the ideal cubic structure is visited even at high temperature, therefore further implying an order-disorder nature of the paraelectric-to-AFE transition (intermediate phase not counted).
The middle (respectively, right) panel of Fig. 3c depicts the PES associated with the changes of R − 5 and a 0 a 0 c − +Γ − 4 (respectively, of R − 5 and a 0 a 0 c + +X − 5 ), and therefore involves the R3c (respectively, Pnma) phase. The coupling between R − 5 and a 0 a 0 c − +Γ − 4 (or a 0 a 0 c + +X − 5 ) is also collaborative, and there exist lowbarrier energetic paths without passing through the cubic structure too, which confirms that the high-temperature cubic phase should be taken as a dynamical average structure, with the system hopping between different states with finite AFD tiltings. This finding is in full accordance with the pair correlation function obtained from MD simulations (Fig. 2a).
Combining a − a − c 0 tilting with other phonon modes can therefore lead to these three low-energy structures, the ground state being Pbam (set as ∆E = 0), the R3c phase with an energy that is only 1 meV/f.u. higher, and ∆E = 32 meV/f.u. for the Pnma phase. For the most important distortion, a − a − c 0 , it is interesting to note that these three structures favor different oop tiling angles (ω R ) about the [100] or [010] axis, i.e., 7.5 • for Pbam, 6.3 • for R3c, and 8.6 • for Pnma. The three panels of Fig. 3d show the contour plots of PES by varying a 0 a 0 c − +Γ − 4 and Σ 2 +S 2 , on top of a − a − c 0 reference structures with oop tilting angles of 6.3 • , 6.9 • and 7.5 • , respectively. These two latter sets of distortions are repulsive to each other, with energy minima thus located on the x or y axis. The minima on the x (or y) axis correspond to the R3c (or Pbam) phase. By decreasing the angle from 7.5 • to 6.3 • , the lowest-energy structure changes from Pbam to R3c. On the other hand, as the angle increases from 7.5 • to 8.6 • (Fig. 3e), the Pnma phase gradually becomes more favorable than the Pbam phase. It is therefore interesting to note that the energy landscape and hierarchy of the minima sensitively depend on the oop tilting angle in PZO.
At very low temperature, the Pbam structure is the ground state with an oop tilting angle of 7.5 • . With increasing T , this average tilting angle will slightly decrease and oscillate around the average value with an increasing amplitude due to thermal excitation. MD simulations at temperatures below T A indicate that, although the stable structure is still Pbam, the deviation of the tilting from the average angle (close to 7.5 • ) can easily be a few degrees even at 100 K (see SM [24]). The reason that the structure remains Pbam is likely due to the kinetic barriers that prevent any transition. The 0K barrier from Pbam to R3c calculated by the CI-NEB method is 57 meV/f.u., and it is 53 mev/f.u. from Pbam to Pnma. This is schematically drawn in Fig. 3b. Furthermore, it is interesting to note that the barriers between equivalent Pbam structures is of similar magnitudes, i.e., 51 (or 58) meV/f.u. for reversing R − 5 (or Σ 2 +S 2 ). This brings to different possible scenario regarding the possible intermediate structure near the AFE transition temperature T A . For instance, the R3c state is energetically very close to Pbam at 0K, and may even have a free energy near T A that is slightly lower than Pbam, which would suggest that the R3c state can be the intermediate structure -which is in line with the ferroelectric intermediate phase with rhombohedral symmetry that was observed in some experimental studies [16-18, 32, 38]. Moreover, the hopping between the minima of Pbam of opposite oop tiltings has rather lower barrier (51 meV/f.u., see left panel of Fig. 3c), which can imply that the intermediate state is an average state of Pbam symmetry too but that does not possess any time-averaged finite a − a − c 0 tilting. However, as no orthorhombic intermediate phase has been reported in experiment, this latter Pbam phase probably does not exist. Another possibility for the intermediate state is a more complex case consisting of an average between R3c and some (or all) Pbam equivalent minima, reflecting the possible hoppings between all these states. Note that the difficulty of determining the symmetry of the intermediate phase in experiments may reflect such picture, i.e., PZO dynamically samples multiple local minima near T A . On the other hand, we do not expect the Pnma state to be (significantely) involved in the intermediate state because its 0K energy is 32 meV/f.u. higher than Pbam.
Furthermore, we want to add some insight about why the Pbam phase is the ground state of PZO, instead of R3c or Pnma, and suggest other possible intermediate states that are too complex to be treated by direct first- principles techniques. To this end, we developed an effective Hamiltonian for PZO, which includes energy terms involving local modes (proportional to the local electric dipole), local oxygen octahedral tiltings, strain, and the coupling among them (see SM for more details [24]); In particular, the three energy terms related to the couplings between the local dipoles and AFD distortions are given by: where the sum runs over all the sites (perovskite cell), with i and j being first nearest neighbors of each other. The three terms are, respectively, the bi-quadratic (only D xyxy here [24]), tri-linear, and bi-linear coupling between the local mode vector u and the pseudo vector characterizing tilting ω. We then carried out Monte Carlo simulations based on the effective Hamiltonian, up to 40,000 MC sweeps at 1 K. We typically adopted a 4×4×2 supercell in terms of 5-atom perovskite cells, which can thus accommodate the Pbam, R3c and Pnma structures. With the DFT determined coefficients, as shown by the vertical lines in Fig. 4, the Pbam phase is indeed found to be the ground state, while the R3c and Pnma phases are 4 and 27 meV/f.u. higher in energy, respectively. This is in excellent agreement with direct DFT calculations. We also constructed, relaxed and stabilized a complex structure with a 16×16×2 supercell, which has a long modulation corresponding to q LM = (2π/a)(3/16, 3/16, 0) [39]. This complex structure is merely 11 meV/f.u. above the Pbam phase, indicating that modulated structures may occur as intermediate phase, which is reminiscent of the incommensurate structure having a similar k-point and experimentally found in PZO under pressure [40]. To see the effect of each type of coupling on structural stability, we individually vary the coefficients D (related to u 2 ω 2type of coupling), D (corresponding to the uω 2 type of couplings), and D (associated with the uω type of couplings), and calculate the total energies (Fig. 4). It is important to note that each type of coupling has significant influence on one specific phase. For instance, larger magnitude of D (or D ) favors the R3c (or Pnma) phase, as it is known before in other perovskites [41,42]. PZO is especially interesting that the strength of the bi-linear coupling uω is found to play an important role to render the Pbam phase ground state, as recently proposed by Patel et al. [15]. These findings also reveal that the uω 2 tri-linear coupling that was proposed in Refs. [10,14,43] to stabilize the Pbam structure may, in fact, not be the main contribution for PZO, as it lowers the energy of Pbam only slightly (Fig. 4b) [44].
Our calculations therefore, e.g., (i) emphasize the order-disorder character (due to dynamical hopping); (ii) point out the importance of bi-linear couplings between the local dipoles and octahedral tiltings; and (iii) suggest possible complex candidates or phenomena for the intermediate phase of PbZrO 3 . We thus hope that our findings provide a better understanding of antiferroelectric perovskites.