Application of Cluster Variation and Path Probability Methods to the Tetragonal-Cubic Phase Transition in ZrO2

TL;DR

This paper extends the cluster variation and path probability methods to study the kinetics and equilibrium states of the displacive tetragonal-cubic phase transition in ZrO2, a novel application of these techniques.

Contribution

It introduces a new approach by applying CVM and PPM to displacive phase transitions, viewing atomic displacements as configurational variables.

Findings

Kinetics of phase transition in ZrO2 analyzed.

Equilibrium states of tetragonal and cubic phases characterized.

Methodology bridges replacive and displacive transition studies.

Abstract

Cluster variation method (CVM) and path probability method (PPM) have generally been employed to study replacive phase transitions in alloy systems. Recently, displacive phase transitions have been explored within the realm of replacive phase transition in the CVM theoretical framework by viewing displaced atoms as different atomic species, i.e., by converting a freedom of atomic displacement to a configurational freedom. The same methodology is applied to the PPM calculations in this work, and the kinetics of displacive phase transition from tetragonal to cubic phases in ZrO2 are investigated as well as their equilibrium states.

Tables

Table 1. Table 1: Pair interaction energies in each sub-lattice e i ​ j α ​ α / α ​ β / β ​ β subscript superscript 𝑒 𝛼 𝛼 𝛼 𝛽 𝛽 𝛽 𝑖 𝑗 e^{\alpha\alpha/\alpha\beta/\beta\beta}_{ij} using e A ​ A subscript 𝑒 𝐴 𝐴 e_{AA} , e A ​ B z subscript superscript 𝑒 𝑧 𝐴 𝐵 e^{z}_{AB} , and e A ​ B x , y subscript superscript 𝑒 𝑥 𝑦 𝐴 𝐵 e^{x,y}_{AB} , where x 𝑥 x , y 𝑦 y , and z 𝑧 z denote the x 𝑥 x , y 𝑦 y , and z 𝑧 z axes, respectively.

	$\alpha$–$\alpha$ (or $\beta$–$\beta$)	$\alpha$–$\beta$ (or $\beta$–$\alpha$)
$e_{AA}$ ($=e_{BB}$)	$e_{AA}$	$e_{AA}$
$e_{AB}$ ($=e_{BA}$)	$e_{AB}^z$	$e_{AB}^{x,y}$

Table 2. Table 2: Fitting parameters of Lennard–Jones potentials for each pair.

	$a_{ij}$	$b_{ij}$	$c_{ij}$
	$\left(\text{Ry}\cdot {\text{au}}^7/\text{atom}\right)$	$\left(\text{Ry}\cdot {\text{au}}^{3.5}/\text{atom}\right)$	$\left(\text{Ry}/\text{atom}\right)$
$e_{AA}$	$3.579\times {10}^6$	$2.461\times {10}^3$	$4.226\times {10}^{-1}$
$e_{AB}^z$	$3.599\times {10}^6$	$2.190\times {10}^3$	$4.277\times {10}^{-1}$
$e_{AB}^{x,y}$	$3.629\times {10}^6$	$2.493\times {10}^3$	$4.227\times {10}^{-1}$

Equations

$$\begin{array}{rl}E& =\frac{1}{2}NZ\sum _{i,j}{e}_{ij}{y}_{ij}\\ & =\frac{1}{2}N\sum _{i,j}(e_{ij}^{\alpha \alpha }{y}_{ij}^{\alpha \alpha }+4e_{ij}^{\alpha \beta }{y}_{ij}^{\alpha \beta }+e_{ij}^{\beta \beta }{y}_{ij}^{\beta \beta }),\end{array}$$

$$e_{ij}(r)=\frac{a_{ij}}{r^7}-\frac{b_{ij}}{r^{3.5}}+c_{ij},$$

$$\begin{split}S=Nk_{B}\;&\mbox{ln}\Bigg{[}\frac{1}{2}\left(\sum_{i}L(x^{\alpha}_{i})+L(x^{\beta}_{i})\right)\\ -\frac{1}{2}&\left(\sum_{i,j}L(y^{\alpha\alpha}_{ij})+4\sum_{i,j}L(y^{\alpha\beta}_{ij})+\sum_{i,j}L(y^{\beta\beta}_{ij})\right)\\ &+\left(\sum_{i,j,k,l}L(z^{\alpha\beta\alpha\beta}_{ijkl})+2\sum_{i,j,k,l}L(z^{\alpha\alpha\beta\beta}_{ijkl})\right)\\ &\;\;\;\;\;\;\;\;\;\;\;-\sum_{i,j,k,l,m,n,o,p}L(w^{\alpha\beta\alpha\beta\alpha\beta\alpha\beta}_{ijklmnop})\Bigg{]}\;,\end{split}$$

$$\left(\frac{\partial F}{\partial w_{ijklmnop}^{\alpha \beta \alpha \beta \alpha \beta \alpha \beta }}\right)=0\text{mbox}and\left(\frac{\partial F}{\partial r}\right)=0.$$

$$\begin{split}\mbox{ln}P_{1}&=\frac{1}{2}N\Big{[}\mbox{ln}\left(\theta\cdot\Delta t\right)\left(\left(X^{\alpha}_{2}+X^{\alpha}_{3}\right)+\left(X^{\beta}_{2}+X^{\beta}_{3}\right)\right)\\ &+\mbox{ln}\left(1-\theta\cdot\Delta t\right)\left(\left(X^{\alpha}_{1}+X^{\alpha}_{4}\right)+\left(X^{\beta}_{1}+X^{\beta}_{4}\right)\right)\Big{]}\;,\end{split}$$

$$\text{mbox}ln{P}_2=-\frac{\Delta E_{\text{mbox}system}}{2k_{B}T},$$

$$\begin{split}\mbox{ln}P_{3}&=N\Bigg{[}\frac{1}{2}\left(\sum_{i}L(X^{\alpha}_{i})+L(X^{\beta}_{i})\right)\\ &-\frac{1}{2}\left(\sum_{i,j}L(Y^{\alpha\alpha}_{ij})+4\sum_{i,j}L(Y^{\alpha\beta}_{ij})+\sum_{i,j}L(Y^{\beta\beta}_{ij})\right)\\ &\;\;\;\;+\left(\sum_{i,j,k,l}L(Z^{\alpha\beta\alpha\beta}_{ijkl})+2\sum_{i,j,k,l}L(Z^{\alpha\alpha\beta\beta}_{ijkl})\right)\\ &\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;-\sum_{i,j,k,l,m,n,o,p}L(W^{\alpha\beta\alpha\beta\alpha\beta\alpha\beta}_{ijklmnop})\Bigg{]}\;,\end{split}$$

$$\theta ={\theta }_0\cdot \text{mbox}exp\left(-\frac{Q}{k_{B}T}\right)$$

$$\begin{split}&\Delta E_{\mbox{\scriptsize{system}}}=\frac{1}{2}N\Big{[}\sum_{i,j}e^{\alpha\alpha}_{ij}\left(y^{\alpha\alpha}_{ij}(t+\Delta t)-y^{\alpha\alpha}_{ij}(t)\right)\\ &\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+4\sum_{i,j}e^{\alpha\beta}_{ij}\left(y^{\alpha\beta}_{ij}(t+\Delta t)-y^{\alpha\beta}_{ij}(t)\right)\\ &\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+\sum_{i,j}e^{\beta\beta}_{ij}\left(y^{\beta\beta}_{ij}(t+\Delta t)-y^{\beta\beta}_{ij}(t)\right)\Big{]}\\ &=\frac{1}{2}N\left[\sum_{i,j}\Delta e^{\alpha\alpha}_{ij}Y^{\alpha\alpha}_{ij}+4\sum_{i,j}\Delta e^{\alpha\beta}_{ij}Y^{\alpha\beta}_{ij}+\sum_{i,j}\Delta e^{\beta\beta}_{ij}Y^{\beta\beta}_{ij}\right],\end{split}$$

$$\eta =\frac{x_A^{\alpha }-x_B^{\alpha }}{2}.$$

Full text

I Introduction

Phase transition is an important and interesting phenomenon in materials science. There are a variety of phase transitions in metallic alloy systems (e.g., magnetic, configurational, and structural phase transitions), and they have been broadly investigated using various computational approaches, such as molecular dynamics simulation, Monte Carlo method, phase field model, and other statistical mechanics methods, such as the cluster variation method (CVM) kikuchi1951theory which is introduced below.

In the statistical mechanics approach, phase stability is analyzed from the free energy, which is generally composed of energy and entropy terms; for example, Helmholtz free energy is defined as $F=E-TS$, where $E$ is the energy, $S$ is the entropy, and $T$ is the temperature. The CVM is one of the most reliable approaches to formulate configurational entropy by considering a broad range of short-range correlations between different atomic species. The CVM has traditionally been combined with an electronic band structure calculation, or density functional theory (DFT), to determine various materials’ properties at finite temperatures without using any experimental data (which is the so-called first-principles CVM mohri2013first ).

Although a rigid lattice is assumed in conventional CVM calculations ignoring local atomic displacements, the continuous-displacement CVM (CDCVM) kikuchi1998space enables incorporating local atomic displacements by introducing “quasi-lattice points” around a Bravais lattice point and viewing the atoms displaced into the quasi-lattice point as different atomic species. The idea of CDCVM (i.e., conversion from a freedom of local atomic displacements to configurational freedom) has been extended to magnetic freedoms yamada2019conversion and collective atomic displacements mohri2013first ; kiyokane2018modelling . In the latter application, for example, the displacive phase transition between tetragonal and cubic phases in ZrO2 is investigated within the scope of replacive transition by regarding upward- and downward-shifted oxygen atoms as different atomic species located on a cubic lattice point.

This idea of conversion of freedom has also been applied to the path probability method (PPM) kikuchi1966path , which is a natural extension of the CVM to the time domain, by the authors of this work yamada2019atomistic . Although the main atomic migration process in alloy systems is by a vacancy mechanism, a spin-flipping mechanism (or a direct-exchange mechanism) is generally assumed in the PPM because of the huge computational burden of the vacancy mechanism. By treating various kinetic paths in the PPM as cluster probabilities in the CVM (i.e., by converting a freedom of kinetic paths to a configurational freedom), the computational cost in the PPM calculation is significantly reduced, and an atomic relaxation process in Ni3Al ordered phase was successfully explored by explicitly considering the vacancy mechanism yamada2019atomistic .

In the present study, the kinetics from tetragonal to cubic phase in ZrO2 as well as their equilibrium states are explored within the CVM and PPM frameworks using the aforementioned methodologies; displacive phase transition is studied within the realm of replacive transition, and path variables in the PPM are treated as cluster probabilities in the CVM. This paper is organized as follows. The methodologies of CVM and PPM are, respectively, described in Secs. II.1 and II.2. The calculated equilibrium states and relaxation paths are shown and discussed in Sec. III. Finally, results of tetragonal–cubic phase transitions in ZrO2 using CVM and PPM calculations with CDCVM are summarized in Sec. IV.

II Theory

II.1 Cluster Variation Method

The tetragonal–cubic transition in ZrO2 has been investigated using the first-principles CVM by one of the authors in reference mohri2013first , where upward- and downward-shifted oxygen atoms along the tetragonal axis are viewed as different atomic species (which are hereafter denoted as A and B atoms, respectively) with fixed zirconium positions. In this methodology, while the tetragonal phase can be described by alternatively displacing oxygen atoms in opposite directions, the cubic phase is realized when the periodicity of the oxygen displacements disappears; thus, this displacive transition can be calculated within the theoretical framework of replacive transitions, more specifically as an order–disorder phase transition mohri2013first ; kiyokane2018modelling . The same methodology and calculation procedure are followed in the present work, and only the essential points are briefly described below.

By considering a cubic structure composed of O atoms in the environment of Zr atoms (see Fig. 1), the total energy, $E$, of the structure is given within the 1st-nearest-neighbor model as

[TABLE]

where $N$ is the number of oxygen atoms, $Z$ is the coordination number ($Z=6$ for the cubic structure), $e^{\alpha\alpha/\alpha\beta/\beta\beta}_{ij}$ and $y^{\alpha\alpha/\alpha\beta/\beta\beta}_{ij}$ are, respectively, the 1st-nearest-neighbor pair interaction energies and pair cluster probabilities between $i$ and $j$ atoms in the $\alpha$–$\alpha$/$\alpha$–$\beta$/$\beta$–$\beta$ sub-lattice (see Fig. 1). In the conventional CVM framework, the pair interaction energies, $e^{\alpha\alpha/\alpha\beta/\beta\beta}_{ij}$, are independent of a sub-lattice. In the displacive calculations, however, they may depend on the sub-lattice. In the case considered here, for example, the pair interaction energy between A (up-shifted) and B (down-shifted) atoms is different between $\alpha$–$\alpha$ (or $\beta$–$\beta$) and $\alpha$–$\beta$ (or $\beta$–$\alpha$) sub-lattices, because the atomic displacements lead to different types of displacements (see Fig. 2). To distinguish the pair interaction energies, $e_{AB}$ depending on the sub-lattices, they are, respectively, represented by $e^{z}_{AB}$ and $e^{x,y}_{AB}$ (see Table 1). This difference in the pair energy was not taken into account in the preliminary work mohri2013first , and the calculated tetragonal–cubic transition temperature of approximately 1500 K, was significantly underestimated compared with that of the experimental data, which is approximately 2570 K aldebert1985structure . The calculated transition temperature is improved using the pair interaction energies shown in Table 1, as shown in Sec. III.

The pair interaction energies, $e_{AA}$, $e^{z}_{AB}$, and $e^{x,y}_{AB}$, are extracted from the total energies of the three different oxygen arrangements calculated from the full potential linear augmented plane-wave (FLAPW) jansen1984total , which are shown in Fig. 3 in reference mohri2013first . The extracted pair interaction energies are fitted to a Lennard–Jones type potential:

[TABLE]

where $a_{ij}$, $b_{ij}$, and $c_{ij}$ are the fitting parameters, and $r$ is the interatomic distance. The fitting parameters are shown in Table 2.

The entropy formula for a cubic approximation in the CVM is given by kiyokane2010order

[TABLE]

where $k_{B}$ is the Boltzmann constant, $x^{\alpha/\beta}_{i}$, $y^{\alpha\alpha/\alpha\beta/\beta\beta}_{ij}$, $z^{\alpha\beta\alpha\beta/\alpha\alpha\beta\beta}_{ijkl}$, and $w^{\alpha\beta\alpha\beta\alpha\beta\alpha\beta}_{ijklmnop}$ are the respective cluster probabilities of the point, pair, square, and cubic clusters at the sub-lattices, and $L(x)\equiv x\;\mbox{ln}x-x$.

The free energy of a system, $F=E-TS$, is formulated using Eqs. (1) and (3). An equilibrium state of a system is determined by minimizing the free energy with respect to the cluster probabilities and interatomic distances:

[TABLE]

The minimization is achieved by the natural iteration method (NIM) kikuchi1974superposition .

II.2 Path Probability Method

The huge computational burden of the PPM can be significantly reduced by treating the path variables as cluster probabilities in the CVM. Whereas the vacancy mechanism is required for an atomic diffusion process, as in reference yamada2019atomistic , an atomic displacement can be modeled within the mold of the spin-flipping mechanism.

The PPM framework is based on the assumption that a system takes the most probable kinetic path during a relaxation process. The most probable path is determined by maximizing a path probability function (PPF). The PPF, $P$, is a product of three terms; $P=P_{1}P_{2}P_{3}$, where $P_{1}$ is the probability of spin-flipping, $P_{2}$ is the probability of thermal activation for a state change, and $P_{3}$ is the number of equivalent paths. The logarithmic expressions of these three terms are written as

[TABLE]

[TABLE]

and

[TABLE]

where $\theta$ is the spin-flipping probability, $\Delta t$ is the time step, $X^{\alpha/\beta}_{i}$, $Y^{\alpha\alpha/\alpha\beta/\beta\beta}_{ij}$, $Z^{\alpha\beta\alpha\beta/\alpha\alpha\beta\beta}_{ijkl}$, and $W^{\alpha\beta\alpha\beta\alpha\beta\alpha\beta}_{ijklmnop}$ are the respective path variables of the point, pair, square, and cubic clusters in the sub-lattices, which are defined in Table 3, and $\Delta E_{\mbox{\scriptsize{system}}}$ is a change in internal energy before and after spin-flippings in $\Delta t$. The spin-flipping frequency, $\theta$, and internal energy change, $\Delta E_{\mbox{\scriptsize{system}}}$, are written as follows:

[TABLE]

and

[TABLE]

where $\theta_{0}$ is the attempt frequency, $Q$ is the activation energy, and $\Delta e^{\alpha\alpha/\alpha\beta/\beta\beta}_{ij}$ are the energy changes of pair interaction energies in $\Delta t$. The maximization of the PPF is performed by the NIM ohno2005iteration as well as minimization of the free energy in the CVM in Sec. II.1.

Note that while the subscripts of cluster probabilities and pair interaction energies indicate atomic species (e.g., $x_{A}$ and $e_{AA}$), those of path variables and the change of pair interaction energies represent a type of kinetic path (e.g., $X_{1}$ and $\Delta e_{11}$). This indicates that the path variables are treated in the same way as the cluster variables in CVM. yamada2019atomistic

III Results and Discussion

The degree of tetragonality compared to a cubic phase is evaluated using the long-range order (LRO) parameter, $\eta$, which is defined as mohri2013first

[TABLE]

The LRO parameters at the equilibrium states are shown in Fig. 3. As can be seen from the results, the transition temperature is estimated to be around 2500 K. This is quite close to the experimental value of 2570 K, unlike the transition temperature reported in reference mohri2013first , where appropriate pair interaction energies are not used for the energy term (see Sec. II.1).

It is controversial whether the tetragonal–cubic phase transition in ZrO2 is first-order or second-order fabris2001free . Although the transition is considered to be weakly second-order in the preliminary study mohri2013first , the calculated temperature dependence of the LRO parameter shown in Fig. 3 represents a first-order characteristic, i.e., it discontinuously becomes zero. In the preliminary and present calculations, however, only the displacements of oxygen atoms are considered with the fixed tetragonality, $c/a$, of the lattice composed of zirconium atoms (see Fig. 1). It has been noted that the degree of tetragonality also depends on the temperature and plays a role in the mechanism of the transition fabris2001free . Therefore, a more detailed analysis is required for a determination of the order of the tetragonal–cubic phase transition in ZrO2 within the CVM framework.

The relaxation process for a sample quenched from $T_{0}=2000$ K to $T_{R}=2400$ K is shown in Fig. 4. It can be seen that the change of the LRO parameter takes place at the initial stage of the relaxation process and the final equilibrium state corresponds to the one independently calculated from the CVM, $\eta\approx 0.7$ (see Fig. 3).

Note that the time scales here are normalized by the spin-flipping frequency, $\theta$. Their real time dependence can be determined once the attempt frequency, $\theta_{0}$, and the activation energy, $Q$, are calculated from the electronic band structure calculation, but this is beyond the scope of the present study and will be considered in the future work.

IV Conclusions

The displacive transitions from tetragonal to cubic phases in ZrO2 are investigated using the CVM and PPM. In both methodologies, the idea of CDCVM is employed; i.e., displacive and kinetic freedoms are converted into configurational freedoms. The transition temperature from tetragonal to cubic phase is estimated to be around 2500 K, which is quite close to the experimental data, 2570 K, and the relaxation process is successfully calculated using the PPM.

One of the biggest advantages of the PPM over other kinetic models (e.g., kinetic Monte Carlo and phase field model) is that the calculated relaxation process can be easily scaled using the data determined from an electronic band structure calculation. It has been shown that there is a certain relation between relaxation processes calculated from PPM and phase field model ohno2006critical . Thus, it would be possible to incorporate a real time dependence into the phase field simulations by effectively combining with the PPM calculation.

ACKNOWLEDGEMENT

This study is partly supported by a project (P16010) commissioned by the New Energy and Industrial Technology Development Organization (NEDO), and by the Structural Materials for Innovation of the Cross ministerial Strategic Innovation Promotion Program (SIP) of Japan Science and Technology (JST).