A new look at the temperature-dependent properties of the antiferroelectric model PbZrO₃: an effective Hamiltonian study

As indicated above, an effective Hamiltonian, H${}_{\text{eff}}$, is developed in the aim of understanding and modeling finite-temperature properties of PbZrO₃ bulk. This H${}_{\text{eff}}$ has the following degrees of freedoms: (1) the local soft modes $\bm{u}_{i}$ in each 5-atom cells $i$, which are proportional to the electric dipoles of that cell Zhong and that are centered on Pb ions here; (2) the $\bm{\omega}_{i}$ pseudo-vectors that describes oxygen octahedral tiltings in the unit cells $i$ Igor2006 and that are centered on Zr ions; $\bm{\omega}_{i}$ is such as its direction is the axis about which the oxygen octahedron of cell $i$ rotates and its magnitude is the angle (in radians) of such rotation; (3) the $\{\eta_{H}\}$ homogeneous strain tensor for which the zero value of its diagonal elements (in Voigt notation) $\eta_{H,1}=\eta_{H,2}=\eta_{H,3}$ is associated with the calculated first-principles-derived lattice constant of the paraelectric cubic state of PbZrO₃ at 0 K ; and (4) $\bm{v}_{i}$ vectors that quantify the inhomogeneous strain at each 5-atom cell $i$  Zhong , and that are centered on Zr ions here.

Following Refs.[NaNbO3Heff, ; CsPbI3Heff, ; BaCeO3Heff, ], the total internal energy contains two different main energies:

where $E^{\text{FE}}$ describes energetics associated with local modes, elastic variables and their interactions, while $E^{\text{tilt}}$ characterizes energetics involving oxygen octahedral tilts and their couplings with local modes and the total $\{\eta_{l}\}$ strain (that is the sum of inhomogeneous and inhomogeneous strains).

As indicated in Ref.[Zhong, ], $E^{\text{FE}}$ can be decomposed into five terms:

where $E^{\text{self}}$ is the local mode self energy, $E^{\text{dpl}}$ pertains to the long-range dipole-dipole interaction, $E^{\text{short}}$ represents the short-range interactions between local modes that go beyond dipole-dipole interactions, $E^{\text{elas}}$ denotes the elastic energy, and $E^{\text{int}}$ describes the interaction between elastic variables and local modes. As proposed in Ref.[Zhong, ], these five energies can be written as follows:

with $\bm{R}_{ij}=\bm{R}_{i}-\bm{R}_{j}$ for which $\bm{R}_{i}$ and $\bm{R}_{j}$ are the lattice vectors locating sites $i$ and $j$, respectively. $\alpha$ and $\beta$ denote Cartesian components along the [100], [010] and [001] directions, that are chosen to span the x-, y-, and z-axes, respectively. The $i$ index runs over all the Pb ions. For any given Pb-site $i$, the $j$ index in $E^{\text{dpl}}$ runs over all the other Pb ions. In contrast, the $j$ index in $E^{\text{short}}$ “only” runs over the first, second and third nearest Pb neighbors of the Pb ion located at the site $i$.

Furthermore, the $J_{ij\alpha\beta}$ coefficients in $E^{\text{short}}$ can be written in the following way for the different nearest neighbor (NN) interactions:

where $\delta$ is the Kronecker symbol and $\hat{\bm{R}}_{ij,\alpha}$ represents the $\alpha$-component of $\bm{R}_{ij}/R_{ij}$. A detailed explanation of the $j_{3}$ coefficient is mentioned in Ref [Zhong, ]. It represents the strength (and sign) of a specific intersite interaction between second nearest neighbor for the local mode.

Regarding the second main energy of Eq. (1), we generalize that of Refs.[NaNbO3Heff, ; CsPbI3Heff, ; BaCeO3Heff, ; SergeyNano, ] by writing:

where the sums over $i$ run over all the Zr sites. The first sum of $E^{\text{tilt}}$ describes the onsite contributions associated with the oxygen octahedral tilts, as proposed in Refs.[Albrecht, ; Igor2007, ; Lisenkov, ; SergeyNano, ]. The second and third terms represent short-range interactions between oxygen octahedral tiltings SergeyNano ; Igor2007 being first-nearest neighbors in the Zr sublattice. Note that $j$ runs over the six Zr sites that are first-nearest neighbors of the Zr site $i$ in the second term while $\omega_{i+\alpha,\alpha}$ in the third term denotes the $\alpha$-component of the $\omega$ pseudo vector at the site shifted from the Zr site $i$ to its nearest Zr neighbor along the $\alpha$ axis. Note also that the $K_{ij\alpha\beta}$ parameters of this second energy are written as:

Here, $k_{2}$ characterizes the strength and sign of a specific interaction between first nearest neighbors for the tilting modes (it is the equivalent of $j_{2}$ shown in Fig. 1 of Ref. [Zhong, ] but for the tilting degrees of freedom rather than the local modes). Moreover and following Ref.[Igor2007, ], the fourth term of $E^{\text{tilt}}$ characterizes the interaction between strain and tiltings. As first proposed in Ref.[SergeyNano, ], the fifth term (that depends on the ${\rm D}_{ij,\alpha\beta}$ coefficients) represents a trilinear coupling between oxygen octahedral tiltings and local modes. The sixth energy, which involves the $E_{\alpha\beta\gamma\delta}$ parameters, characterizes bi-quadratic couplings between such tiltings and modes Igor2007 . Note that the $j$ index runs over the eight Pb ions that are first nearest neighbors of the Zr-site $i$ in these fifth and sixth terms. Finally, the seventh term is the novelty here with respect to the energies provided in Refs.[Albrecht, ; Igor2007, ; Lisenkov, ; SergeyNano, ]. It represents a recently discovered energy that couples oxygen octahedral tiltings and local modes in a bi-linear fashion Kinnaryuw ; Reviewcouplingenergies . It is allowed by symmetry and was proposed to explain the occurrence of complex antiferroelectric, ferrielectric and even incommensurable phases in some materials Kinnaryuw ; Reviewcouplingenergies ; BinTanPZO . Note that the sum over j is about the eight Zr-sites that are nearest neighbors of the Pb-sites i. It is given, in a less compact form than in Eq. (5), by

where $\mathrm{F_{ii}}$ is a material-dependent constant that characterizes the strength of this coupling. The sum over $i$ runs over all the 5-atom cells of the perovskite structure, and the $x$, $y$, and $z$ subscripts denote the Cartesian components of the ${\bf u}_{i}$ vectors and $\boldsymbol{\omega}_{i}$ pseudo-vectors. $\omega_{ilmn}$, with l, m or n being 0 or 1, characterizes the $\omega$ pseudo-vectors located at $a(l\hat{x}+m\hat{y}+n\hat{z}$ from that of site $i$), with $a$ being the 5-atom lattice parameter and $\hat{x}$, $\hat{y}$ and $\hat{z}$ being the unit vectors along the x-, y- and z-axes, respectively. Note that the Zr site $i$ is located at $-a(\frac{1}{2}\hat{x}+\frac{1}{2}\hat{y}+\frac{1}{2}\hat{z})$ with respect to the Pb site $i$. Figure 1 schematizes the coupling terms inherent to the first line of Eq. (II.1). This new coupling term is important in some systems having Pb atoms (Ref [Kinnaryuw, ]), hence relaxors.

All the parameters of this effective Hamiltonian are provided in Table 1. They are determined by conducting first-principles calculations using the local density approximation (LDA) LDA within density functional theory. Cells containing up to 40 atoms are employed, along with the CUSP code Laurent3 and the ultrasoft-pseudopotential scheme David with a 25 Ry plane-wave cutoff. Practically, we use the same pseudopotentials than those of Ref.[DomenicKS, ], and thus consider the following valence electrons: Pb 5d, Pb 6s, Pb 6p, Zr 4s, Zr 4p, Zr 4d, Zr 5s, O 2s and O 2p electrons. It is important to know that these effective Hamiltonian parameters are determined by considering small perturbations with respect to the cubic state, and thus do not involve the full ionic and cell relaxation of any phase Zhong . Consequently, the fact that we used LDA to obtain these parameters does not automatically imply that our energetic results for relaxed structures predicted by the effective Hamiltonian will be closer to those resulting from the direct use of LDA than to those obtained from other functionals such as the generalized gradient approximation (GGA) in the form of revised Perdew-Burke-Ernzerhof pbesolPerdew2008 (PBEsol). This is especially true if employing these two functionals within first-principles calculations provides similar results in term of structure but rather small differences between the total energy of different relaxed phases – as it is, in fact, known in PbZrO₃ bulk HAramberri2021 .

This H${}_{\text{eff}}$, along with its parameters, is then used in Monte-Carlo (MC) computations on $12\times 12\times 12$ supercells (which contains 8,640 atoms) for different temperatures. Typically, 40,000 MC sweeps are conducted for each considered temperature, with the first 20,000 sweeps allowing the system to be in thermal equilibrium and the next 20,000 sweeps being employed to obtain statistical averages. However, near phase transitions, a larger number of MC sweeps is needed to have converged results, namely up to 1 million (with the first half used for reaching thermal equilibrium and the next half to extract statistical averages). It is also important to know that a temperature-dependent pressure is added in these Monte-Carlo simulations, with this pressure, P, acting on the strains to produce the following energy:

where $\eta_{H,i}$, with $i$ = 1, 2 and 3, being the diagonal elements of the homogeneous strain tensor in Voigt notations Voigt1910 (note that the shear strains are neglected in Equation (8) because they are either null for the cubic paralectric phase or found to be rather small for the other phases encountered in the simulations). Such dependence is assumed here to be linear, as proposed in Ref.[tintePvsT, ], and is given, in GPa, by: $P=-5.489655172-0.001034483\times T$, where $T$ is the temperature in Kelvin with the two coefficients having been numerically found by trying to reproduce the experimental lattice constants at 10 K and 300 K.

In order to determine which structural phases are predicted from the effective Hamiltonian simulations, the following quantities are extracted from the outputs of the MC simulations at any investigated temperature:

where ${\bf k}$ are vectors belonging to the cubic first Brillouin zone, $\alpha$ denotes Cartesian components and the sums run over all the 5-atom sites. Typically, for the local modes, we look at the ${\bf k}$-vectors at the zone-center (for the polarization) but also at the $\Sigma$-point (for complex antipolar displacements associated with the Pbam state) that is defined as $\Sigma=(\frac{1}{4},\frac{1}{4},0)$ in $2\pi/a$ units. Such latter point is also investigated for the tilting of the oxygen octahedra, in addition to the R-point that is given by $R=(\frac{1}{2},\frac{1}{2},\frac{1}{2})$ still in $2\pi/a$ units. Note that a non-zero $\omega_{R,\alpha}$ characterizes antiphase tilting about the $\alpha$-axis while a finite $\omega_{\Sigma,\beta}$ should also happen for the Pbam ground state with $\beta$ being different from $\alpha$.

We also used direct first-principles calculations to check some unexpected predictions arising from the use of the effective Hamiltonian and that will be discussed in the Results’ section. Practically, we use the generalized gradient approximation (GGA) within the revised Perdew-Burke-Ernzerhof functional (PBEsol) pbesolPerdew2008 , as implemented in the VASP package Kresse1999 , for these first-principles computations. The projector augmented wave (PAW) Blochl1994 ; Kresse1999 is also applied to describe the core electrons, and we consider the Pb ($5d^{10}$$6s^{2}6p^{2}$), Zr ($4s^{2}$$4p^{6}5s^{2}4d^{2}$) and O ($2s^{2}$$2p^{4}$) as valence electrons with a 520 eV plane-wave cutoff.

The lattice parameters have been determined by X-ray diffraction using a high resolution diffractometer with a Copper radiation ($\lambda$=Cu${}_{\textrm{K}_{\alpha,1}}$=1.5405 Å) issued from an 18 kW-rotating anode from room temperature to 750K using a furnace with a resolution better than 0.1 K. The X-ray diagrams were analyzed fitted with a pseudo-tetragonal unit cell due to the overlap of the peaks related to the $a$ and $b$ lattice parameters of the orthorhombic unit cell. The lattice parameters at 10 K and 300 K have been obtained by Rietveld analysis (Jana software) Rietveld on results obtained from neutron diffraction performed at the Laboratoire Léon Brillouin (beam line 3T2 , $\lambda$ = 1.2252 Å). The lattice parameters obtained at room temperature from neutron diffraction and from X-ray diffraction are in very good agreement. The orthorhombic unit cell parameters (a${}_{\textrm{o}}$, b${}_{\textrm{o}}$, c${}_{\textrm{o}}$) are presented in Fig.(3b) in a pseudo-tetragonal (a${}_{\textrm{pt}}$, c${}_{\textrm{pt}}$) setting where a${}_{\textrm{o}}=\sqrt{2}$a${}_{\textrm{pt}}$, b${}_{\textrm{o}}\approx 2\sqrt{2}$a${}_{\textrm{pt}}$ and c${}_{o}=2$c${}_{\textrm{pt}}$. At room temperature the orthorhombic cell parameters are a${}_{\textrm{o}}$=5.878Å, b${}_{\textrm{o}}$=11.783Å, c${}_{\textrm{o}}$=8.228Å  while, at 10 K, they are a${}_{\textrm{o}}$=5.878Å, b${}_{\textrm{o}}$=11.784Å, and c${}_{\textrm{o}}$=8.197Å.

Let us now report the predictions of the effective Hamiltonian with the parameters indicated in Table 1 and the total energy described in Equation (1), with the additional feature that the second line of Equation (7) is dropped out. This dropping is made because we numerically found that incorporating all the terms of Eq.(7) provides unphysical solutions for tilting arrangements (such as tiltings associated with the $X$-point of the first Brillouin zone, which is given by $(0,0,\frac{1}{2})$ in $2\pi/a$ units) while keeping the first eight terms of Eq.(7) gives rise to a selected Pbam ground state (Note that dropping all terms of Eq.(7) will not yield such latter ground state).

Figure 2(a) reports the temperature evolution of the non-zero tilting-related quantities, namely $\omega_{R,x}=\omega_{R,y}$ and $\omega_{\Sigma,z}$, while Figure 2(b) shows that the only significant local mode quantity is $u_{\Sigma,x}=u_{\Sigma,y}$, as obtained by heating up the system from the ground state (note that identical results are obtained by cooling down the system, and that the H${}_{\text{eff}}$ calculations also provide a $S_{4}$ mode, but for which the weight is rather small, namely less than 0.32%). The fact that all these quantities are finite at small temperature does characterize a Pbam ground state.

As the temperature increases up to $\simeq$ 650 K, $u_{\Sigma,x}=u_{\Sigma,y}$ decreases while $\omega_{\Sigma,z}$ is barely affected and $\omega_{R,x}=\omega_{R,y}$ (anti-phase oxygen octahedra tilts) gets reduced significantly until vanishing through a jump. A first-order transition to what we call here the $\Sigma$-phase, and that is only characterized by non-zero $u_{\Sigma,x}=u_{\Sigma,y}$ and $\omega_{\Sigma,z}$ (complex tilts along z direction), occurs. Such latest phase still has the Pbam symmetry with space group no. 55 (as indicated by the FINDSYM software FINDSYM ; HTStokes2005 ), therefore indicating that our predicted transition around 650 K is isostructural in nature VRajan2005isostrucutral ; Tian2018 ; Kennedy1999Isostructuralphasetransition . Note that isostructural transitions are usually of first-order, which is consistent with the jump of $\omega_{R,x}=\omega_{R,y}$ seen here at 650 K SergeyPZO2014 . Interestingly, it is experimentally known that the Pbam ground state phase disappears around 511 K, which is not so far away from our predicted 650 K – especially taking into account that effective Hamiltonians have the tendency to not accurately predict transition temperatures (while qualitatively reproducing phase transition sequences, as demonstrated, e.g., by Ref.[Zhong, ]). However, we are not aware of any previous work mentioning that the Pbam ground-state phase transforms into another Pbam state (that is, the $\Sigma$ phase here) under heating, therefore making our predictions provocative and novel. Within this $\Sigma$ phase, Figs 1a and 1b indicate that both $\omega_{\Sigma,z}$ and $u_{\Sigma,x}=u_{\Sigma,y}$ decrease with temperature and then vanish through a jump at about 1300 K. A first-order transition to cubic paraelectric $Pm\bar{3}m$ is thus predicted to happen, according to our effective Hamiltonian simulations.

In order to determine if these rather surprising predictions of the H${}_{\text{eff}}$ can be realistic, Fig. 3 reports the computed $\chi_{11}$, $\chi_{22}$ and $\chi_{33}$ diagonal elements of the dielectric tensor as well as their average $\chi_{av}=\frac{1}{3}(\chi_{11}+\chi_{22}+\chi_{33})$ as a function of temperature. One can see that $\chi_{33}$ and thus $\chi_{av}$ adopt very large values at the temperatures around which the known Pbam ground state transforms into the novel $\Sigma$ phase, which is consistent with the experimental observation of a large dielectric response around the temperatures at which the known Pbam ground state disappears RomanFayeRef11 ; RomanFayeRef14 ; RomanFayeRef46 ; RomanFayeRef48 ; RomanFayeRef55 ; GShirane1951DEResponse ; RomanFayeRef69 ; RomanFayeRef70 ; RomanFayeRef80 ; RomanFayeRef81 . Strikingly too, the coefficient $C$ in the fitting of $\chi_{33}$ and $\chi_{av}$ by $\frac{C}{T-T_{0}}$ is predicted to be 1.81$\times$10⁵ K and 1.49$\times$10⁵ K, respectively, in the $\Sigma$ phase, which is of the same order of magnitude than the experimental values ranging between 1.36$\times$10⁵ K and 2.07 $\times$10⁵ K (see Refs.[RomanFayeRef46, ; RomanFayethesis, ] and references therein) for temperatures at which the known Pbam ground state has vanished.

Note that we found that the increase of $\chi_{av}$ and $\chi_{33}$ when the temperature gets reduced in the $\Sigma$ phase towards the known Pbam state can be thought to originate from the fact that there is a $P4mm$ ferroelectric state (with a polarization along the $z$-axis and no oxygen octahedral tilting) that is very close in free energy for these temperatures. In fact, reducing the $F_{ii}$ parameter in Table I to a certain critical value of 0.019 makes the intermediate phase in-between the known Pbam state and the $Pm\bar{3}m$ cubic state becoming $P4mm$ rather than $\Sigma$. Note also that the large dielectric response shown in Fig. 3 implies that the zone-center mode should soften when approaching the known Pbam phase from above, which has been indeed observed OstapchukPolarandCentralmodePZO2001 ; TagantsevPZO . It is also worth mentioning that the proximity in energy of the ferroelectric P4mm state (or other ferroelectric states) could also explain the experimental electric-field- or defect-driven single ferroelectric P-E loop observed above the known Pbam phase HLiu2011 , or stabilized under epitaxial strain LPintilie2007 .

Let us now take a look at the behavior of the pseudocubic lattice parameters as a function of temperature. Figure. 4a reports our predictions from the use of the presently developed effective Hamiltonian (obtained from the strain outputs) while Fig. 4b displays our own experimental results along with those of Ref.[Fujishita2003, ]. One can see that, within the known Pbam state, several features of the measurements are well reproduced, namely the $a$ and $b$ pseudocubic lattice constants are basically independent of the temperature while the $c$ pseudocubic lattice parameter significantly increases when heating PZO. The jump in the $c$ lattice constant seen in Fig. 4a at around 650 K confirms the first-order character of the isostructural transition between the known Pbam state and the $\Sigma$ phase. What is remarkable is that, within the error bars (the range of error bars for the $a$ and $b$ pseudocubic lattice parameter is from 0.0013 to 0.007 Å for the whole temperature range, while it is less than 0.0038 Å for the $c$ lattice parameter), the $a$, $b$ and $c$ parameters are identical within our predicted $\Sigma$ phase, which may thus give the (wrong) impression that PZO is cubic if solely focusing on pseudocubic lattice parameters or strains. One can thus think that this $\Sigma$ phase is a cubic state in disguise, when focusing on pseudocubic lattice parameters. In fact, the predicted temperature evolution of $\frac{1}{3}(a+b+c)$ adopts a thermal expansion that is very close in the $\Sigma$ and paraelectric $Pm\bar{3}m$ phases, and that is characterized by an average coefficient of 7.97 $\times$ 10^-6 K^-1, –which is in good agreement with the present measurements giving 7.4$\times$ 10^-6 K^-1 (note that such coefficient is numerically found to be 7.77 $\times$ 10^-6 K^-1 in our $\Sigma$ phase versus 8.36 $\times$ 10^-6 K^-1 in the simulated $Pm\bar{3}m$ state). This facts suggests, once again, that the present simulations can be realistic. On the other hand, comparing Fig. 4a and 4b tells us that the jump in the $c$ pseudocubic lattice parameter when the known Pbam state disappears is more pronounced in experiments than in our present calculations. Such quantitative discrepancy can be due to the facts that this disappearance happens at higher temperature in the simulations and/or that the computations use relatively small $12\times 12\times 12$ supercells (implying that the first-order character of that transition can be reduced by the increase in critical temperature and/or finite-size effects).

Interestingly, near the temperatures at which the known $Pbam$ state disappears, it has been observed in Ref [71] (see Fig. 2.9 there) that the $(\frac{1}{2},\frac{3}{2},\frac{1}{2})$ diffractions peak typically associated with antiphase oxygen octahedral tiltings disappears at a temperature lower than the $(\frac{1}{4},\frac{3}{4},0)$ diffraction peak that typically characterizes antiparallel displacements of Pb associated with the $\Sigma$ k-point. Such features appear to be qualitatively consistent with our predictions shown in Figs 1a and 1b. Note that our predicted $\Sigma$-phase can be of short-range or broken into small domains in grown samples, which may explain why it has not been directly reported and reported yet. It is also interesting to realize that, in addition to a zone-center soft-mode, a central mode, with a frequency being very low, has been experimentally reported for temperatures above the known Pbam state OstapchukPolarandCentralmodePZO2001 ; TagantsevPZO ; JHKo2013Modesoftintermediatephase2013 . Such feature may be representative of dynamical jumps between different $\Sigma$-phases, that are with positive or negative $\omega_{\Sigma,z}$, $\omega_{\Sigma,y}$ and $\omega_{\Sigma,x}$. Such dynamical jumps will render cubic the overall structure of PZO, while not affecting $\frac{1}{3}(a+b+c)$.

Let us now check if such $\Sigma$ phase is also seen as a low-energy state within direct first-principles calculations. Table 2 reports the energy of such $\Sigma$ phase, but also those of the known Pbam and $R3c$ states, choosing the zero of energy to correspond to the cubic paraelectric $Pm\bar{3}m$ state. One can first see that the known Pbam is the state with the lowest energy but by only a small amount of about 1 meV/f.u. with respect to $R3c$. Very interestingly, Table 2 also tells us that $\Sigma$ is indeed a low-energy state in PZO, with its decrease in energy being about 235.6 meV/f.u. with respect to the cubic state, to be compared with 278.8 meV/f.u. for the known Pbam phase. In other words, having finite $\omega_{\Sigma,z}$ and $u_{\Sigma,x}=u_{\Sigma,y}$ (like in the $\Sigma$ and known Pbam phases) brings about 85% of the energy of the known Pbam state with respect to $Pm\bar{3}m$ state, with the other 15% being due to $\omega_{R,x}=\omega_{R,y}$. Such percentages highlight the importance of the bilinear couplings of Eq. (7), which are responsible for the finite values of $\omega_{\Sigma,z}$ and $u_{\Sigma,x}=u_{\Sigma,y}$ in both the $\Sigma$ and known Pbam states, as well as their low energies. Table 3 provides the crystal structure and atomic positions of the known Pbam and the presently discovered $\Sigma$ phase, respectively, as predicted by direct first-principles calculations at 0 K, in the hope they can be used in the future to investigate and/or revisit structural phases of PZO. One can see from Table 3 that the pseudo lattice parameters of the $\Sigma$ phase are not close to each other at 0K, while they are at high temperatures (see Fig 4a).