Order-disorder transition in the prototypical antiferroelectric PbZrO$_3$

TL;DR

This study uses first-principles simulations to clarify the order-disorder nature of the antiferroelectric transition in PbZrO₃, highlighting the role of octahedral tilting and local interactions in its ground state.

Contribution

It introduces a first-principles-based effective Hamiltonian revealing the importance of bi-linear interactions in the Pbam ground state of PbZrO₃.

Findings

No phonon becomes truly soft at the transition.

Order-disorder transition is confirmed.

Bi-linear interactions are crucial for the ground state.

Abstract

The prototypical antiferroelectric PbZrO$_3$ 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, which contrasts with most of the previous explanations.

Equations

$$\Delta E^{\text{AFD}}=$$

$$+$$

Full text

Antiferroelectrics (AFEs) form a special class of materials that possess anti-polar displacements, which was first conceptualized by Kittel Kittel (1951). AFEs have received much attention in recent years because they hold great potential to reach high energy density for high-power energy storage Hao (2013); Rabe (2013); Xu et al. (2017). Lead zirconate PbZrO3 (PZO) is the first discovered AFE material, and crystallizes in a Pbam ground state below $\sim$505 K Sawaguchi et al. (1951); Shirane et al. (1951); Shirane (1952). 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.

Despite the efforts devoted to understanding the origin of the AFE ground state in PZO Tagantsev et al. (2013); Íñiguez et al. (2014); Hlinka et al. (2014), how distortions condense from the high-symmetry cubic structure remains an open issue. Unlike proper ferroelectrics where the ferroelectric (FE) distortion results from soft polar mode in the cubic phase, it is presently unclear if PZO is a proper AFE. As a matter of fact, on one hand, completely soft AFE and AFD modes have been predicted at the AFE transition temperature ($T_{A}$) by Fthenakis et al. based on an effective Hamiltonian method Fthenakis and Ponomareva (2017). On the other hand, no truly soft AFE or AFD modes was observed in Ref. Tagantsev et al. (2013). In this regards, 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 Tagantsev et al. (2013); however, this mechanism is rather exotic and has recently been ruled out Vales-Castro et al. (2018). Moreover, from an energetic point of view, Íñiguez et al.. believe that the AFD mode plays a crucial role and PZO may be an improper AFE Íñiguez et al. (2014).

In fact, the above considerations in Refs. Fthenakis and Ponomareva (2017); Íñiguez et al. (2014) rely on the existence of a displacive transition, for which the distortions in the antiferroelecrtic state are derived from sudden deviations from the ideal cubic structure at $T_{A}$. However, such assumption may be invalid for PZO, as the AFE transition have been advocated to possess a strong order-disorder character, as evidenced by (i) the existence of central mode Ostapchuk et al. (2001); Ko et al. (2013) and (ii) the fact that disordered Pb displacement and oxygen tiltings were found to occur even in the cubic phase Paściak et al. (2015); Zhang et al. (2015). Interestingly, this picture of dynamic disorder in the cubic phase is similar to the recent studies of the AFD transitions in CaMnO3 and inorganic halide perovskites Klarbring and Simak (2018); Klarbring (2019), where low-energy paths between local minima are found to be critical for understanding the finite-temperature transitions.

Furthermore, there exists a FE R3c phase that is energetically very close to the AFE state, merely a few meV/f.u. higher than Pbam at 0 K, according to ab initio calculations Reyes-Lillo and Rabe (2013); Íñiguez et al. (2014). Hence, it is interesting to understand what microscopic effect is mostly responsible for the stabilization of Pbam rather than R3c as the ground state. In particular, is it the commonly believed trilinear energy coupling Hlinka et al. (2014); Íñiguez et al. (2014), or rather something else Patel et al. (2016)? Another unsettled issue in PZO concerns the possible existence of an intermediate phase in the vicinity of $T_{A}$ Sawaguchi et al. (1951); Tennery (1965, 1966); Goulpeau (1967); Fujishita (1992); Liu (2018): such intermediate state can only occur over a very narrow temperature range of 3–5 K between the AFE and the cubic phase, and its structure remains a mystery.

In this article, we use first-principles methods to demonstrate that the order-disorder picture may hold the key to understand the AFE transition and the debated intermediate phase Sawaguchi et al. (1951); Tennery (1965, 1966); Goulpeau (1967); Fujishita (1992) in PZO. We also reveal that a very specific bi-linear (and therefore not trilinear) energetic coupling plays a crucial role in the stabilization of Pbam and can also yield other more complex structures, such as with long-periods and incommensurable ones.

Phonons in the cubic phase at high temperatures are investigated by combining ab initio molecular dynamics (AIMD) simulations implemented in vasp Kresse and Joubert (1999) with the temperature dependent effective potential (TDEP) method Hellman et al. (2011); Hellman and Abrikosov (2013); Hellman et al. (2013). The potential energy surface (PES) is calculated by density functional theory (DFT) Kresse and Joubert (1999), and the kinetic barriers are computed with the climbing image nudged elastic band (CI-NEB) method Henkelman et al. (2000). We also develop an effective Hamiltonian for PZO to explain the microscopic origin of its complex ground state and to suggest other more complex candidates for the intermediate state. Further details of the methods are given in the Supplemental Material (SM) sm .

A purely displacive picture requires truly soft phonon mode. 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 AIMD simulations with a 4$\times$4$\times$4 supercell (320 atoms). In contrast to the 0K phonon computed by density functional perturbation theory (DFPT) that exhibits soft branches Ghosez et al. (1999); Cockayne and Rabe (2000); Leung et al. (2002); Leung (2003) (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 not (a). Interestingly, this latter temperature is rather close to the experimental AFE transition temperature $T_{A}$ of $\sim$505 K.

Moreover, the 0K phonon exhibits the largest instabilities at the $\Gamma$, $R$, and $M$ points of the Brillouin zone, that correspond to $\bm{q}_{\Gamma}=(2\pi/a)(0,0,0)$, $\bm{q}_{R}=(2\pi/a)(1/2,1/2,1/2)$, and $\bm{q}_{M}=(2\pi/a)(1/2,1/2,0)$, respectively (with $a$ being the 5-atom lattice constant) and that are associated with the FE distortion, out-of-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 $\Gamma$-point mode Ostapchuk et al. (2001).

Interestingly, no phonon becomes truly soft (that is none of them drops to zero frequency) when the temperature approaches $\simeq$ 450K from above, as consistent with experimental data in the cubic phase Tagantsev et al. (2013); Ostapchuk et al. (2001). In particular, Ref. Tagantsev et al. (2013) pointed out that the $\Sigma$- and $R$-point modes are not soft at $T_{A}$, where $\bm{q}_{\Sigma}=(2\pi/a)(1/4,1/4,0)$, which is consistent with our prediction in Fig. 1b, but contrasts with the theoretical work of Ref. Fthenakis and Ponomareva (2017) that predicted completely softened $\Sigma$ and $R$ modes.

The lack of soft phonon suggests that the paraelectric-to-AFE transition (without counting the intermediate phase) is of first order, as also consistent with experiments Samara (1970); Whatmore and Glazer (1979); Fujishita et al. (2003). Moreover, in order to check if this transition has characters of an order-disorder type, as suggested by Ref. Ostapchuk et al. (2001) 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 sm for more information). As the AIMD simulation is considered to be much more realistic than cubic stochastic sampling not (b), 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 sm ). 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 ($\Gamma^{-}_{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 of 222 meV/f.u. is provided by the $R^{-}_{5}$ mode, whereas the $M^{+}_{2}$ and $\Gamma^{-}_{4}$ modes cause less energy change, i.e., 167 and 57 meV/f.u., respectively. 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 this mode plays a crucial role in stabilizing the AFE ground state Íñiguez et al. (2014). It is also worth to mention that all the low-energy structures reported in Ref. Íñiguez et al. (2014) (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. Nevertheless, the large barrier heights of reversing $R^{-}_{5}$ and $M^{+}_{2}$ imply that the transition is not likely to be displacive, as thermal energy at $T_{A}$ is not sufficient to bypass the ideal cubic positions, even for condensing a single tilting mode. Such fact is related to the finite tilting angles in the cubic phase from MD simulations.

Next, we combine $a^{-}a^{-}c^{0}$ tilting with other distortions to yield three competing polymorphs having low energies, viz., the Pbam (ground state), R3c, and Pnma phases. Using the crystallographic tool amplimode Orobengoa et al. (2009); Perez-Mato et al. (2010), the Pbam structure is mainly composed of three distortions: $R^{-}_{5}$, $\Sigma_{2}$ and S2. 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 $\Sigma_{2}$ and S2, respectively. To simplify the picture, we combine the $\Sigma_{2}$ and S2 modes and compute the 2-dimensional PES by changing $R^{-}_{5}$ and $\Sigma_{2}$+S2. The resulting contour plot is shown in the left panel of Fig. 3b, and exhibits four minima that correspond to the Pbam structure not (c) – which is a manifestation of collaborative couplings between the $R^{-}_{5}$, $\Sigma_{2}$ and S2 modes. The energy of these minima relative to the cubic structure (center) is 311 meV/f.u., being 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 $\Sigma_{2}$+S2, respectively. Improved estimation of the barrier heights are obtained by the CI-NEB method, yielding 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 very 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. 3b depicts the PES associated with the changes of $R^{-}_{5}$ and $a^{0}a^{0}c^{-}$+$\Gamma^{-}_{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^{-}$+$\Gamma^{-}_{4}$ (or $a^{0}a^{0}c^{+}$+X${}^{-}_{5}$) is also collaborative, and there exist low-barrier 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 $\Delta E=$ 0), the R3c phase with an energy that is only 1 meV/f.u. higher, and $\Delta E=$ 32 meV/f.u. for the Pnma phase. We will discuss the couplings that gives rise to the AFE ground state later in the text, but here one should note that the lower energy of Pbam than R3c depends on the precise structure and chemistry. Slight change of structure or chemical composition (e.g., Ti doping) can affect the energy hierarchy, giving rise to a FE R3c ground state Ayyub et al. (1998); Boldyreva et al. (2007); Mani et al. (2015a); Woodward et al. (2005).

Although the R3c and Pnma phases do not occur in bulk PZO at any temperature, their minima are likely to be visited near $T_{A}$ or above. To see this, we calculate the 0K barrier from Pbam to R3c by the CI-NEB method, which is 57 meV/f.u., and the barrier is 53 mev/f.u. from Pbam to Pnma. Interestingly, the barriers between equivalent Pbam structures is of similar magnitudes, i.e., 51 (or 58) meV/f.u. for reversing $R^{-}_{5}$ (or $\Sigma_{2}$+S2).

This explains the occurrence of the intermediate phase in a very narrow temperature range, and the sensitive dependence on impurities, domain orientation, etc. Tennery (1965, 1966); Goulpeau (1967); Liu (2018), as it only occurs when the system samples a fraction of the energy minima. It also brings to different possible scenario regarding the possible intermediate structure near the AFE transition temperature $T_{A}$ (Fig. 3c). 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 Tennery (1965, 1966); Goulpeau (1967); Scott and Burns (1972); Whatmore and Glazer (1979). 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. 3b), 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.

So far we have demonstrated that deep local minima and the existence of low-barrier energetic paths bypassing the cubic positions are crucial to understand the finite-temperature structures in PZO. Now, 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 sm ); In particular, the three energy terms related to the couplings between the local dipoles and AFD distortions are given by:

[TABLE]

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 sm ), tri-linear, and bi-linear coupling between the local mode vector ${\bf u}$ and the pseudo vector characterizing tilting ${\bm{\omega}}$.

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$\times$4$\times$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$\times$16$\times$2 supercell, which has a long modulation corresponding to $\bm{q}_{\text{LM}}=(2\pi/a)(3/16,3/16,0)$ not (d). 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 Burkovsky et al. (2017). To see the effect of each type of coupling on structural stability, we individually vary the coefficients $D$ (related to $u^{2}\omega^{2}$-type of coupling), $D^{\prime}$ (corresponding to the $u\omega^{2}$ type of couplings), and $D^{\prime\prime}$ (associated with the $u\omega$ 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^{\prime}$ ) favors the R3c (or Pnma) phase, as it is known before in other perovskites Kornev and Bellaiche (2009); Bellaiche and Íñiguez (2013). PZO is especially interesting that the strength of the bi-linear coupling $u\omega$ is found to play an important role to render the Pbam phase ground state, as recently proposed by Patel et al. Patel et al. (2016), in particular for perovskites with A-site cations that tend to move off-center, e.g., Pb2+ or Bi3+ with lone-pair electrons. These findings also reveal that the $u\omega^{2}$ tri-linear coupling that was proposed in Refs. Hlinka et al. (2014); Íñiguez et al. (2014); Prosandeev et al. (2014) 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) not (e).

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 PbZrO3. We thus hope that our findings provide a better understanding of antiferroelectric perovskites, and suggest possible design of new antiferroelectrics having Pbam or related structures, by searching for novel systems with large bi-linear coupling.