# A deep-neural network potential to study transformation-induced plasticity in zirconia

Jin-Yu Zhang, Gaël Huynh, Fu-Zhi Dai, Tristan Albaret, Shi-Hao Zhang, Shigenobu Ogata, David Rodney

► To cite this version:

Jin-Yu Zhang, Gaël Huynh, Fu-Zhi Dai, Tristan Albaret, Shi-Hao Zhang, et al. A deep-neural network potential to study transformation-induced plasticity in zirconia. Journal of the European Ceramic Society, 2024, 44 (6), pp.4243-4254. 10.1016/j.jeurceramsoc.2024.01.007. hal-04494023

HAL Id: hal-04494023

https://hal.science/hal-04494023v1

Submitted on 7 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

Original article

# A deep-neural network potential to study transformation-induced plasticity in zirconia

Jin-Yu Zhang<sup>a</sup>, Gaël Huynh<sup>a</sup>, Fu-Zhi Dai<sup>b</sup>, Tristan Albaret<sup>a</sup>, Shi-Hao Zhang<sup>c</sup>, Shigenobu Ogata<sup>c</sup>, David Rodney<sup>a,*</sup>

<figure><img src="image_4.png" /></figure>

## ARTICLE INFO

### Keywords:

Zirconia

Machine learning potential

Phase transformation

Molecular dynamics

Plasticity

## ABSTRACT

Zirconia (ZrO₂) ceramics uniquely exhibit transformation-induced plasticity, allowing plastic deformation prior to failure, setting them apart from most other ceramics. However, our understanding of ZrO₂ plasticity is hindered by the challenge of simulating stress-induced atomic-scale phase transformations, owing to lack of an efficient interatomic potential accurately representing polymorphism and phase changes in ZrO₂. In this work, we introduce a novel deep neural network interatomic potential, DP-ZrO₂, constructed using a concurrent-learning approach. DP-ZrO₂ reproduces properties of various ZrO₂ phases, matching their phase diagrams as well as transformation pathways with accuracy comparable to *ab initio* density functional theory. Leveraging DP-ZrO₂, we conducted molecular dynamics simulations of temperature-induced interphase boundary migration and nanocompression. These simulations demonstrate the potential's efficiency and applicability in studying deformation microstructures involving phase transformations in ZrO₂. Our approach opens the door to large-scale simulations under complex loading conditions, which will shed light on the conditions favouring ZrO₂ transformation-induced plasticity.

## 1. Introduction

Zirconia (ZrO₂) ceramics have attracted significant attention in biomedical and energy applications due to their ductility and fracture toughness remarkably larger than other ceramics [1–3]. This unique combination of mechanical properties is closely linked to phase transformations [3,4], which enhance the resistance to crack propagation (transformation-induced toughening) and increase the ductility (transformation-induced plasticity) [5–7]. Experimental and theoretical studies were conducted to optimize the composition and processing parameters of ZrO₂ ceramics, aiming to establish a relationship between properties and transformation microstructures [8–10]. However, this research remained mostly empirical because of our lack of basic understanding regarding the nanoscale processes that control phase transformations in ceramics in general and in ZrO₂ ceramics in particular. Key questions remain, such as the onset of phase transformations under different loading orientations, the kinetics of interface migration during phase transformations, and the mechanisms of transformation propagation across grain boundaries.

Molecular dynamics (MD) is a powerful tool to study phase transformations at the atomic scale. In ZrO₂ ceramics, previous MD simulations

investigated the competition between dislocations and phase transformations under different loading orientations [11,12] and the nanoindentation of tetragonal ZrO₂ (t-ZrO₂) polycrystals [13] or nanopillars with pre-existing defects [14]. However, these studies used an empirical Born–Mayer–Buckingham potential with a Coulombic term, which does not describe accurately the tetragonal-to-monoclinic (m-ZrO₂) transformation, the martensitic transformation induced under stress in ZrO₂. Furthermore, empirical potentials face challenges in describing other experimentally observed ZrO₂ phases, such as cubic ZrO₂ (c-ZrO₂), Pbca-ZrO₂, Pmna-ZrO₂, and Pbc2₁-ZrO₂ [15,16]. These phases can alter the transformation paths from t-ZrO₂ [17–19] and have been reported to reduce transformation toughening [4,17,18]. Therefore, a numerically-efficient interatomic potential that accurately captures the polymorphism in ZrO₂ is necessary for large-scale simulations of deformation in ZrO₂ ceramics.

In recent years, machine-learning interatomic potentials have gained significant attention and have successfully addressed various challenges in different fields. For instance, machine-learning potentials have been developed for ZrO₂ to study piezo- and pyro-electricity [20], phase transitions at ambient pressure [21,22] and oxygen diffusion [23,24].

<footnote><sup>a</sup> Institut Lumière Matière, Université Claude Bernard Lyon 1, Villeurbanne, 69622, France</footnote>

<footnote><sup>b</sup> Artificial Intelligence for Science Institute, Beijing, 100084, People's Republic of China</footnote>

<footnote><sup>c</sup> Department of Mechanical Science and Bioengineering, Osaka University, Osaka, 560-8531, Japan</footnote>

<footnote>* Corresponding author.
*E-mail address:* david.rodney@univ-lyon1.fr (D. Rodney).</footnote>
Among these machine-learning potentials, deep potentials (DP) stand out as promising candidates for the present study. DP has been applied to various systems, achieving accuracy comparable to density functional theory (DFT) calculations. Examples include ultra-high temperature ceramics [25], (de)lithiation processes in silicon anodes [26], and the phase diagrams of ice [27] and gallium [28]. The efficiency of DP can be significantly enhanced by a compression technique [29], enabling large-scale simulations up to millions of atoms [30,31], which is suitable for investigating transformation-induced plasticity in ZrO₂ ceramics.

In this study, we employ a DP approach [32] to construct an atomistic model of ZrO₂ (Section 2), called here Deep-Potential ZrO₂, or simply DP-ZrO₂, that can accurately describe experimentally observed ZrO₂ phases over a wide range of temperatures (0 to 4500 K) and pressures (–5 to 20 GPa). In Section 3, we systematically examine various properties of ZrO₂ polymorphs and investigate the energy barriers controlling phase transformations, with comparisons to DFT results. The phase diagram of DP-ZrO₂ is determined using thermodynamic integration. Finally, we demonstrate the accuracy and efficiency of the DP potential to simulate relevant mechanical properties of ZrO₂ systems by studying for the first time the migration of an experimentally observed interface between t-ZrO₂ and m-ZrO₂, and by performing a model nanocompression under a controlled orientation. We emphasize that the current DP model focuses on pure ZrO₂ systems, while real applications mostly involve doped ZrO₂. Nonetheless, this work represents a mandatory step forward to model ceramics at the atomic scale. The current DP model and corresponding DFT dataset will soon be available at a publicly accessible website [33].

## 2. Computational methods

The machine-learning DP model, proposed by Zhang et al. [32], was employed to train an interatomic potential for ZrO₂ using a dataset obtained from DFT calculations. In the DP model, the total energy of a system is computed as a sum of atomic contributions, $E = \sum_i E_i$. The atomic energies $E_i$ depend on the local environment, $E_i = E(R^i)$, where $R^i = (r_{1i}, r_{2i}, ..., r_{mi})$ represents the relative coordinates of neighbouring atoms around atom $i$ within a cutoff distance $r_c$. The mapping from $R^i$ to $E_i$ is obtained thanks to an embedding and a fitting neural network. The embedding network transforms $R^i$ into a descriptor matrix $D^i$ that preserves the symmetry operations of the system, while the fitting neural network maps $D^i$ into the atomic energy $E_i$ [34].

We used the smooth edition of the DeepPot descriptor [34] as implemented in the DeePMD-kit package [35,36]. The cutoff was set to 7.5 Å with a smoothing function starting at 2.0 Å.

The embedding network is composed of three hidden layers with sizes (25, 50, 100) and the fitting network comprises three hidden layers with sizes (240, 240, 240). The dimension of the descriptor matrix is 16 × 100. During the training process, the following loss function was minimized [35,36]:

$$L(p_e, p_f, p_\xi) = \frac{p_e}{N} \Delta\epsilon^2 + \frac{p_f}{3N} \sum_i |\Delta F_i|^2 + \frac{p_\xi}{9N} \|\Delta\xi\|^2, \quad (1)$$

where $N$ is the number of atoms. $\Delta\epsilon$, $|\Delta F_i|$, and $\|\Delta\xi\|$ denote the differences between DP model and DFT calculations on energy, atomic force, and virial tensors, respectively. $p_e$, $p_f$, and $p_\xi$ are corresponding prefactors and are adjusted during the training process, which is performed in two stages, each consisting of eight million steps. During the first stage, the prefactors of $p_e$, $p_f$, and $p_\xi$ change gradually from 0.02 to 1, 1000 to 1, and 0.02 to 1, respectively. In the second stage, the pretrained DP model is fine-tuned with prefactors $p_e$, $p_f$, and $p_\xi$ fixed to 100, 1, and 10, respectively. With an initial prefactor of 1000, forces are already well predicted at the end of stage 1. Stage 2 is then a fine-tuning step where the accuracy on energies and stresses is improved without degrading the force accuracy [27,37].

The DFT dataset was prepared using a concurrent-learning scheme implemented in DP-GEN (Deep Generator) [38]. The DP-GEN software

iteratively explores the configurational space through three steps: training, exploration, and labelling. In the training step, DP-GEN invokes the DeePMD-kit package to train four DP models with different activation functions based on an initial DFT dataset. The fine-tuning stage is skipped to accelerate the DP-GEN iteration. The initial dataset for the first training step consists of 638 ZrO₂ configurations with random distortions and atomic displacements. The configurations are supercells of 96 atoms for the six experimentally observed ZrO₂ crystal structures, namely monoclinic, tetragonal, cubic, Pbca, Pmna, and Pbc2₁ (see Fig. 1). Additionally, the initial dataset includes structures generated by the generalized solid-state nudged elastic band (G-SSNEB) method [39] to better capture the transition states involved in the phase transformations.

In the exploration step, DP-GEN uses LAMMPS [40] to perform MD simulations at various temperatures and pressures to explore configurational space. The temperature is incrementally increased from 100 K to 4500 K during the DP-GEN iterations, and the explored pressures range from –5 to 20 GPa. To capture surface effects, surface configurations corresponding to (100)m, (010)m, (001)m, {100}t, (001)t, {101}t, and {110}t at temperatures below 500 K are also included into the dataset. For the explored configurations, all the four DP models will give predictions of atomic forces. If a configuration is not well described by the current DP models, they will provide different predictions and exhibit a large model deviation, which is defined as the maximum standard deviation of the predicted atomic forces. In the present study, the configurations with model deviations ranging from 0.15 eV/Å to 0.35 eV/Å are randomly selected and added to the training dataset. In total, 5057 configurations are included into the dataset. For liquid structures, the range is adjusted to 0.30–0.50 eV/Å. An upper limit of the model deviation is set to avoid adding possible unphysical configurations to the dataset.

In the labelling step, DP-GEN invokes the Vienna ab initio simulation package (VASP) [41,42] to calculate the energy, atomic forces, and stresses of the selected configurations. The Perdew–Burke–Ernzerhof exchange–correlation functional revised for solids (PBEsol) [43] was chosen for the calculations. This functional was shown to be a good trade-off between accuracy and numerical cost [44], representing cohesion in ZrO₂ much better than the generalized gradient approximation (GGA) for a numerical cost much lower than the strongly constrained and appropriately normed parameterization (SCAN) of meta-GGA [45]. A plane-wave cutoff of 900 eV and 3 × 3 × 3 Monkhorst–Pack k-mesh are set to ensure convergence of the energies, atomic forces, and stresses within 1 meV/atom, 1 meV/Å, and 1 meV/atom, respectively, for both distorted crystals and liquid structures.

## 3. Validation

### 3.1. Structural and vibrational properties

The DP model shows a good agreement with DFT results with root-mean-square errors (RMSE) of 3.1 meV/atom and 0.2 eV/Å for energy and force components, respectively, comparable to previous potentials for similar systems [20,21,46] (see Fig. 2).

Table 1 compares the lattice parameters and energies of ZrO₂ polymorphs at 0 K, demonstrating a relative error of about 0.1% for lattice parameters and energy deviations of less than 1 meV/atom. Moreover, the predicted elastic constants of ZrO₂ crystals exhibit an excellent agreement, typically within 10 GPa, as summarized in Table 2. Fig. 3 displays the results of DP and DFT calculations for the equation of states, revealing a remarkable agreement. In particular, with both DFT and DP, the most stable crystals follow the sequence m-ZrO₂ → Pbca-ZrO₂ → Pmna-ZrO₂ as the volume decreases. These calculations are consistent with high-pressure experiments that observed phase transformations from monoclinic to Pbca and from Pbca to Pmna with increasing hydrostatic pressures [47,48].
<figure><img src="image_5.png" /><figcaption>**Fig. 1.** (a) Monoclinic, (b) Pbca, (c) Pbc2₁, (d) tetragonal, (e) cubic and (f) Pmna crystalline phases. Atoms are coloured according to their coordination number. Zr and O atoms are represented by large and small atoms, respectively. *a*-axis and *b*-axis of the unit cells are parallel to the horizontal direction and viewing direction, respectively.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Comparison of the energy (a) and forces (b) computed by DFT and predicted by the DP model. Insets show error distributions.</figcaption></figure>

We employed the Phonopy package [49] to evaluate the vibrational properties of ZrO₂ polymorphs. Fig. 4 displays the phonon spectra of the low-pressure phases, while Supplementary Materials present high-pressure polymorphs. To ensure accurate phonon calculations, the initial configurations consisting of conventional cells made of 12 atoms, with the exception of the Pbca polymorph, which consisted of a 24-atoms cell, were relaxed with both the DP model and DFT until the norm of all forces was below 10⁻⁶ eV/Å. To ensure convergence with respect to system size, 96 atoms supercells are used for the finite displacement calculations. Due to the sensitivity of the second-order derivatives of the energy, achieving an accurate description of phonon dispersion through the training process was challenging, but we can see in Fig. 4 that the DP potential reproduces very well the DFT data. We also note that both DFT and DP correctly predict the instability of the cubic phase at 0 K [50].

The lattice parameters evolution with temperature of the monoclinic, tetragonal and cubic phases is shown in Fig. 5. These results were extracted from 300 ps long simulations using 6 × 6 × 6 supercells (2592 atoms) in the isothermal–isobaric ensemble at zero pressure. The equilibrium phase at each given temperature was taken as starting point. The linear thermal expansion coefficients are in excellent agreement with the experiments.

With the DP model and in the experiments, the lattice parameters are discontinuous at the m-ZrO₂/t-ZrO₂ transition, implying that the transition is of first-order as expected for a martensitic transformation. As shown in Verdi et al. [21], near the t-ZrO₂/c-ZrO₂ transition temperature, the system undergoes frequent fluctuations between both phases. However, by computing separately the lattice parameters of the tetragonal and cubic phases at each temperature (as discussed below), we find that the lattice parameters are discontinuous at the transition
<table><caption>**Table 1**
Lattice parameters and energy of different ZrO₂ phases calculated by DFT and DP (f.u. stands for formula unit)</caption><thead><tr><th></th><th></th><th>a (Å)</th><th>b (Å)</th><th>c (Å)</th><th>β (°)</th><th>ΔE (meV/f.u.)</th></tr></thead><tbody><tr><td rowspan="2">m-ZrO<sub>2</sub></td><td>DFT</td><td>5.133</td><td>5.215</td><td>5.300</td><td>99.6</td><td>0.0</td></tr><tr><td>DP</td><td>5.136</td><td>5.216</td><td>5.295</td><td>99.7</td><td>0.0</td></tr><tr><td rowspan="2">Pbca-ZrO<sub>2</sub></td><td>DFT</td><td>10.055</td><td>5.252</td><td>5.075</td><td>90.0</td><td>40.9</td></tr><tr><td>DP</td><td>10.036</td><td>5.253</td><td>5.079</td><td>90.0</td><td>39.8</td></tr><tr><td rowspan="2">Pbc2<sub>1</sub>-ZrO<sub>2</sub></td><td>DFT</td><td>5.055</td><td>5.260</td><td>5.075</td><td>90.0</td><td>52.2</td></tr><tr><td>DP</td><td>5.051</td><td>5.261</td><td>5.071</td><td>90.0</td><td>52.4</td></tr><tr><td rowspan="2">t-ZrO<sub>2</sub></td><td>DFT</td><td>5.074</td><td>5.074</td><td>5.178</td><td>90.0</td><td>79.2</td></tr><tr><td>DP</td><td>5.073</td><td>5.073</td><td>5.178</td><td>90.0</td><td>77.6</td></tr><tr><td rowspan="2">c-ZrO<sub>2</sub></td><td>DFT</td><td>5.068</td><td>5.068</td><td>5.068</td><td>90.0</td><td>140.4</td></tr><tr><td>DP</td><td>5.069</td><td>5.069</td><td>5.069</td><td>90.0</td><td>139.1</td></tr><tr><td rowspan="2">Pmna-ZrO<sub>2</sub></td><td>DFT</td><td>5.558</td><td>3.321</td><td>6.478</td><td>90.0</td><td>202.2</td></tr><tr><td>DP</td><td>5.555</td><td>3.323</td><td>6.478</td><td>90.0</td><td>199.7</td></tr></tbody></table>

<caption>**Fig. 3.** Equation of states of different ZrO₂ phases. Solid lines and crosses indicate the results by DP-ZrO₂ model and DFT, respectively.</caption>

temperature, suggesting a first-order transition. This is in agreement with experimental data but in contradiction with the results proposed by Verdi et al. [21].

To further confirm the nature of the t-ZrO₂/c-ZrO₂ transformation, we conducted longer simulations of 10 ns to precisely measure the volume and enthalpy of both c-ZrO₂ and t-ZrO₂ phases in the vicinity of the transformation point (see Figs. 6(a) and (b)). It is worth noting that the system fluctuates between c-ZrO₂ and t-ZrO₂ phases near the transformation temperature, consistent with previous findings [21]. This fluctuation is directly visualized by the probability distribution of the instantaneous lattice parameters ($L_x, L_y, L_z$) (Fig. 6(c)) at different timesteps of MD simulations, which indicates that the system explores the c-ZrO₂ phase and the three t-ZrO₂ variants. In contrast to the previous study in Ref. [21], where a statistical average was computed across all structures, we adopted a different approach. Here, we classify structures at different timesteps as either c-ZrO₂ ($c/a \approx 1.00$) or t-ZrO₂ ($c/a \approx 1.02$) based on their instantaneous $c/a$ values and compute separate statistical averages for each phase. We defined $c$ as the longest lattice parameter at a given timestep and $a$ as the average of the other two lattice parameters. Specifically, a structure is identified as t-ZrO₂ if $c/a > 1.01$ and as c-ZrO₂ if $c/a < 1.01$. Our results clearly reveal small differences between both phases in terms of volume and enthalpy around the transformation temperature, supporting the conclusion that

this is indeed a first-order phase transformation. The same methodology was applied in Fig. 5 to compute the temperature-dependence of the average lattice parameters.

We further examined the free energy surface associated with the c-ZrO₂/t-ZrO₂ transformation based on the distribution of lattice parameters presented in Fig. 6(c). In this analysis, we considered the lattice parameters ($L_x, L_y, L_z$) as the collective variables, and we calculated the free energy surface in terms of lattice parameters as [54]:

$$F(L_x, L_y, L_z) = -kT \log(p(L_x, L_y, L_z)). \quad (2)$$

Here, $p(L_x, L_y, L_z)$ is the normalized probability distribution of the system. We can see that the lattice parameters of the system approximately lie in the (111) plane ($L_x + L_y + L_z = \text{constant}$), since the system vibrates around its average volume and the volume difference between both phases is very small (about 0.06% at 2390 K). Therefore, the dimension of the collective variables is reduced by projecting the data points onto the (111) plane and the probability distribution of the system on this plane is computed. Using Eq. (2), the calculated free energy surface at 2390 K is reported in Fig. 6(d), which shows a free energy barrier of about 0.1 meV/atom. The free energy barrier of the transformation always exists around the transformation temperature (Fig. 6(e)) and is much smaller compared with other typical first-order transitions in ZrO₂ (see Section 3.2), which indicates that the t-ZrO₂/c-ZrO₂ transformation is a weak first-order phase transformation. The energy barrier and time during which the system resides in a specific phase are anticipated to rise with system size, potentially surpassing the timescale of MD simulations. This aligns with the previously observed overheating of the tetragonal phase, as documented in Table S1 in Ref. [55].

Other physical properties at high temperatures, including the diffusion coefficient and heat capacity (which shows a Schottky anomaly in c-ZrO₂), were computed with the present potential. They are reported in Supplementary Materials.

## 3.2. Minimum energy paths of phase transformations

Zirconia can transform between different phases through diffusion-less transformations. The two most studied are the ferroelastic switching between t-ZrO₂ variants [56,57] and the martensitic transformation between the t-ZrO₂ and m-ZrO₂ phases [3,4]. At zero applied stress and temperature, the various phases are separated by energy barriers that can be determined using saddle-point search methods. Here, we used G-SSNEB method [39] to evaluate the energy barrier for transformations between different phases and/or transformation variants using both DFT and the DP model. Calculations were performed to compute the minimum-energy path between an initial and a final unit cell made of 12 atoms with periodic boundary conditions. We thus consider here homogeneous transformations. The case of heterogeneous transformations will be addressed in Section 3.5. The result of the G-SSNEB calculations is shown in Fig. 7.

The first barrier on the left corresponds to the ferroelastic switching between two tetragonal variants of ZrO₂, i.e. a transformation of the type ($a_t, b_t, c_t) \rightarrow (c_t, a_t, b_t)$. As seen in Fig. 7, the transformation occurs between two variants that have the same energy (taken as reference for all G-SSNEB calculations) and the DP model reproduces the energy barrier obtained with DFT with an error of 1 meV/f.u.. The second barrier (in green) corresponds to a martensitic transformation from a tetragonal unit cell ($a_t, b_t, c_t$) into a monoclinic unit cell ($a_m, b_m, c_m$). This transformation has an orientation relationship of the so-called C-type [58], where the c-axis of the tetragonal phase transforms into the c-axis of the monoclinic phase. We see in Fig. 7, that both phases are separated by a single energy barrier slightly higher than that of ferroelastic switching (16.49 meV/f.u. compared to 16.30 meV/f.u.). We note that with an LDA functional, the DFT calculations predict a very low ferroelastic switching barrier (10 meV/f.u.) and a significantly lower
<table><caption>**Table 2**
Elastic constants (GPa) calculated by DFT and DP.</caption><thead><tr><th></th><th></th><th>C<sub>11</sub></th><th>C<sub>22</sub></th><th>C<sub>33</sub></th><th>C<sub>12</sub></th><th>C<sub>13</sub></th><th>C<sub>23</sub></th><th>C<sub>44</sub></th><th>C<sub>55</sub></th><th>C<sub>66</sub></th><th>C<sub>15</sub></th><th>C<sub>25</sub></th><th>C<sub>35</sub></th><th>C<sub>46</sub></th></tr></thead><tbody><tr><td rowspan="2">c-ZrO<sub>2</sub></td><td>DFT</td><td>548</td><td></td><td></td><td>105</td><td></td><td></td><td>69</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>DP</td><td>542</td><td></td><td></td><td>115</td><td></td><td></td><td>71</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">t-ZrO<sub>2</sub></td><td>DFT</td><td>473</td><td></td><td>308</td><td>133</td><td>69</td><td></td><td>36</td><td></td><td>77</td><td></td><td></td><td></td><td></td></tr><tr><td>DP</td><td>470</td><td></td><td>303</td><td>125</td><td>63</td><td></td><td>37</td><td></td><td>76</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">Pbc2<sub>1</sub>-ZrO<sub>2</sub></td><td>DFT</td><td>364</td><td>380</td><td>373</td><td>129</td><td>120</td><td>153</td><td>122</td><td>87</td><td>84</td><td></td><td></td><td></td><td></td></tr><tr><td>DP</td><td>362</td><td>372</td><td>385</td><td>123</td><td>113</td><td>149</td><td>129</td><td>105</td><td>91</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">Pbca-ZrO<sub>2</sub></td><td>DFT</td><td>357</td><td>420</td><td>356</td><td>171</td><td>130</td><td>133</td><td>98</td><td>90</td><td>126</td><td></td><td></td><td></td><td></td></tr><tr><td>DP</td><td>368</td><td>416</td><td>374</td><td>171</td><td>140</td><td>140</td><td>108</td><td>94</td><td>130</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">Pmna-ZrO<sub>2</sub></td><td>DFT</td><td>458</td><td>356</td><td>409</td><td>159</td><td>193</td><td>164</td><td>65</td><td>98</td><td>127</td><td></td><td></td><td></td><td></td></tr><tr><td>DP</td><td>456</td><td>358</td><td>400</td><td>168</td><td>193</td><td>152</td><td>75</td><td>98</td><td>137</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">m-ZrO<sub>2</sub></td><td>DFT</td><td>310</td><td>397</td><td>252</td><td>169</td><td>91</td><td>152</td><td>86</td><td>88</td><td>130</td><td>48</td><td>-5</td><td>14</td><td>-9</td></tr><tr><td>DP</td><td>311</td><td>404</td><td>270</td><td>172</td><td>117</td><td>153</td><td>96</td><td>95</td><td>124</td><td>49</td><td>-1</td><td>24</td><td>-6</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Phonon dispersion curves of (a) m-ZrO₂, (b) t-ZrO₂ and (c) c-ZrO₂. Red and blue lines refer to DFT and DP, respectively.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Pressure-free temperature-dependent lattice parameters for the low-pressure polymorphs of ZrO₂. Filled regions represent the standard deviation, while crosses and dots refer to experimental data and the DP model, respectively. Experimental data from Refs. [51–53]. The vertical dashed lines highlight the disparity between experimental transition temperatures and the values obtained from MD simulations.</figcaption></figure>

energy difference between t-ZrO₂ and m-ZrO₂ phases [59], whereas the GGA functional PBE predicts a high ferroelastic switching barrier (29.2 meV/f.u.) and an extremely low t-ZrO₂→m-ZrO₂ (2.6 meV/f.u.), as well as a large energy difference between both phases [60].

Interestingly, if we look at another orientation relationships, for example the B-type where $(a_t, b_t, c_t) \rightarrow (c_m, a_m, b_m)$, i.e. the c-axis of the t-ZrO₂ phase becomes the b-axis of the m-ZrO₂ phase, or the A-type where $(a_t, b_t, c_t) \rightarrow (b_m, c_m, a_m)$, i.e. the c-axis of the t-ZrO₂ phase becomes the a-axis of the m-ZrO₂ phase, we find in both DFT and DP calculations that there is no direct path between the t-ZrO₂ and m-ZrO₂. Instead, the t-ZrO₂ phase undergoes first a ferroelastic switch to align its c-axis with the final c-axis of the m-ZrO₂ phase and then the unit cell undergoes a C-type transformation. The energy path then corresponds to the succession of the blue and green paths in Fig. 7.

Another path for the B orientation relationship was proposed by Trolliard et al. [61]. This path is also indirect with the initial t-ZrO₂ phase first transforming into an intermediate orthorhombic phase Pbc2₁ of orientation $(a_0, c_0, b_0)$ before transformation into m-ZrO₂ with orientation $(a_m, c_m, b_m)$. The complete transformation $(a_t, b_t, c_t) \rightarrow (a_0, c_0, b_0) \rightarrow (a_m, c_m, b_m)$ is represented as the succession of the red and purple lines in Fig. 7. While the energy barrier for the Pbc2₁ to monoclinic transformation is very high compared to others (77.6 meV/f.u.), the energy barrier from t-ZrO₂ to Pbc2₁ is smaller and comparable to the barrier of the t-ZrO₂ → m-ZrO₂ path. Given the ~4.5% volume increase along the t-ZrO₂ → m-ZrO₂ path, the t-ZrO₂ → Pbc2₁ transformation may be favoured considering the local stress-strain state.

Finally, there is another possible path between the t-ZrO₂ and m-ZrO₂ phases, shown as the succession of the red and yellow paths, and going through a Pbc2₁ phase. As seen in Fig. 7, the initial t-ZrO₂ phase transforms into a Pbc2₁ unit cell with orientation $(c_0, a_0, b_0)$ then back to a tetragonal phase with a c-axis parallel to the final c-axis of the monoclinic phase (i.e. an indirect ferroelastic switching through a Pbc2₁ phase) and then the tetragonal phase undergoes a C-type transformation towards the final monoclinic phase. Interestingly, the ferroelastic switching of the intermediate tetragonal phase only appears when a strict convergence criterion of $10^{-3}$ eV/Å is used, which may explain why this intermediate state has not been evidenced in a previous study [60].

Although there already exist multiple paths for the t-ZrO₂ → m-ZrO₂ transformation at zero pressure, we note that the energy barriers are expected to evolve as a function of the stress-strain state in the system. This point has already been quoted by Guan et al. [60] and should be investigated in further studies. We also note that other transition paths exist but were not discussed as they involve very high energy barriers. Overall, in all the cases presented here, we see that the DP model reproduces very accurately the transition paths and corresponding barriers compared to the DFT calculations.

### 3.3. Phase diagram of ZrO₂

Fig. 8 presents the phase diagrams of ZrO₂ obtained from previous experiments [47] and predicted by the DP model. The computation
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Analysis of the phase transition between c-ZrO₂ and t-ZrO₂. The volume (a) and enthalpy (b) of c-ZrO₂ and t-ZrO₂ near the transformation temperature. (c) Distribution of lattice parameters at 2390 K. (d) Free energy surface of the t-ZrO₂/c-ZrO₂ transformation at 2390 K. (e) Free energy curves between t-ZrO₂ and c-ZrO₂ at different temperatures along the dashed line in (d).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Minimum energy paths predicted by DFT (empty circles) and DP-ZrO₂ (solid lines).</figcaption></figure>

covered a temperature range of 0–3700 K and a pressure range of 0–20 GPa. As a first step, we employed thermodynamics integration [62] (details in Supplementary Materials) to determine transformation temperatures at zero pressure. We obtained transformation temperatures for m-ZrO₂ → t-ZrO₂, t-ZrO₂ → c-ZrO₂, and c-ZrO₂ → liquid ZrO₂

equal to 1610 K, 2400 K and 2630 K, respectively. These values agree reasonably well with experimental data (1367 K [63], 2584 K [64], and 2983 K [55], respectively). Using the calculated transformation points as starting points, we then conducted Gibbs-Duhem integrations [65] to obtain coexistence lines up to triple-points. We note that the Pbc2₁ phase does not appear in the calculated phase diagram, which we verified by performing thermodynamics integrations at different pressures. We confirmed that the Pbc2₁ phase is never the most stable phase in the phase diagram at prescribed pressures and temperatures.

While a quantitative difference between the experimental and DP model phase diagrams is observed in Fig. 8, it is worth noting that, to the best of our knowledge, this is the first time that a complete phase diagram of ZrO₂ has been obtained using an interatomic potential. One reason is that no interatomic potential developed so far was able to predict the stability of the various phases of ZrO₂. The predicted transformation sequences, m-ZrO₂ → t-ZrO₂ → c-ZrO₂ → liquid with increasing temperature at 0 Pa and m-ZrO₂ → Pbca-ZrO₂ → Pmna-ZrO₂ with increasing pressure at low temperatures, perfectly agree with the sequences observed experimentally. The deviation of the calculated phase diagram from the experimental one is likely due to inaccuracies of the PBEsol exchange-correlation functional used in the DFT dataset. In comparison to more accurate results [44] obtained using the van
<figure><img src="image_11.png" /><figcaption>Fig. 8. Phase diagrams of ZrO₂ obtained experimentally (a) (reproduced from Ref. [47]) and calculated with the present DP model (b).</figcaption></figure>

der Waals corrected SCAN (SCAN-rVV10) [66], the PBEsol predicts a higher energy for t-ZrO₂, potentially accounting for the narrower region of t-ZrO₂ at 0 Pa in the calculated phase diagram. Previous DFT calculations using PBE functional also predicted a smaller Pbca region at low temperatures compared to experiments [67].

Another obvious difference between the experimental and DP phase diagrams is that the t-ZrO₂/c-ZrO₂ coexistence line depends significantly on pressure experimentally, while it is almost flat numerically. The reason is that experimentally, there is a volume increase from t-ZrO₂ to c-ZrO₂ of about 0.55% [55] whereas the DP model yields a value of only 0.06% at 2400 K (Fig. 6(a)), resulting in a very weak pressure dependence of the transition temperature. Given that DFT calculations using SCAN-rVV10 are time-consuming and resource-intensive, additional efforts are needed to enhance the DP model and achieve better experimental accuracy.

In addition to deviations arising from simulation results, it is important to note that debates exist regarding transformation pressures of high-pressure phases observed in experimental studies. Experimental measurements of transformation pressures conducted at room temperature may exhibit discrepancies due to kinetic effects. The studies supporting the reversibility of phase transitions [48,68] reported lower transformation pressures, i.e., 3.1 GPa for m-ZrO₂ → Pbca-ZrO₂ transformations [68] and 12.5 GPa for Pbca-ZrO₂ → Pmna-ZrO₂ transformations [48], compared with the values shown in Fig. 8(a). Moreover, in experiments by Ohtaka et al. [48], Pmna-ZrO₂ was found stable in the high pressure and high temperature region (up to 1800 °C at 24 GPa), which should indicate that the slopes of phase boundaries between the Pmna and t-ZrO₂ phases and between Pmna and c-ZrO₂ phases are high. This result is inconsistent with the experimental phase diagram in Ref. [47] used in Fig. 8(a) but agrees well with the present DP-phase diagram.

### 3.4. Temperature-induced migration of a t-ZrO₂/m-ZrO₂ interface

With our accurate DP model at hand, we can explore the structure and migration of the t-ZrO₂/m-ZrO₂ interfaces, pushing the boundaries of our model out of the training dataset and further validating its applicability in complex interface simulations. The orientation relationship between the two phases and the interface normal are calculated based on the theory of martensite crystallography [70,71]. A coherent interface can be obtained when [010]ₘ and [001]ₜ are corresponding vectors

<figure><img src="image_12.png" /><figcaption>Fig. 9. Comparison between a simulated t-ZrO₂/m-ZrO₂ interface structure and an experimental scanning transmission electron microscopy image from Ref. [69]. An inverse Fourier transform was performed on the experimental image to remove the noise. For clarity, only Zr atoms are shown.</figcaption></figure>

through the phase transformation and the small mismatch (about 0.7%) between them is accommodated elastically [4]. Based on the calculated orientation relationship and interface normal, a coherent interface is constructed and simulated. In Fig. 9, we considered as an example a so-called B-1 orientation relationship, where [100]ₘ || [010]ₜ, [010]ₘ || [001]ₜ, (001)ₘ || (100)ₜ. A slight prestrain (less than 1%) is imposed on the crystals to satisfy the periodic boundary condition. The interface normal is parallel to (0.110 0 0.994)ₘ || (0.994 0.110 0)ₜ and the simulated interface is composed of terraces parallel to (001)ₘ || (100)ₜ and atomic steps with a step vector [001]ₘ or [100]ₜ. The simulated interface structure at 300 K is superimposed on an atomic image obtained experimentally using scanning transmission electron microscopy in yttria-stabilized zirconia [69] in Fig. 9. The comparison can only be qualitative but the very good agreement between atomic positions in the atomistic model and the experiment is an indication of the quality of the DP potential. In both cases, the interface has a highly coherent
<figure><img src="image_13.png" /><figcaption>Fig. 10. Shear-coupled interface migration at 300 K of a t-ZrO₂/m-ZrO₂ interface with a B-type orientation relationship.</figcaption></figure>

structure, with identical misorientations on both sides of the interface. Deviations of atom positions can be seen locally, which may result from local strains in the experimental sample and the difference in lattice parameters between the simulation (pure ZrO₂) and the experiment (3 mol% Y₂O₃ doped ZrO₂).

The migration behaviour of the t-ZrO₂/m-ZrO₂ interface is further simulated in Fig. 10 by performing an MD simulation at 300 K without applied stress. The interface migration is thus driven by the free energy difference between both phases. The migration rate is around 500 m/s at room temperature, which is a typical value during an athermal martensitic transformation. We see that the interface simply translates in the simulation cell, without altering its structure or leaving behind crystalline defects. The development of a shear and a dilatation of the simulation cell are clearly visible in Fig. 10 accompanying the interface migration. They are consequences of the shear and dilatation components of the Bain transformation between both phases, which form the physical origin of transformation-induced plasticity in ZrO₂.

Interestingly, we see that t-ZrO₂ can undergo a B-type transformation in presence of a t-ZrO₂/m-ZrO₂ interface, while there is no direct path for B-type transformation at the level of unit cells, as discussed in Section 3.2. Observation of the B-type transformation in the present simulation as well as in experiments thus probably results from the heterogeneous structure of the bi-crystal near the interface. The coherent interface was shown to stabilize a gradual path between the two phases across the interface, as demonstrated in a prior study on graphite-to-diamond transformations [72]. This is reinforced by the high coherency of the t-ZrO₂/m-ZrO₂ interface in the B-type transformation. We have modelled other interfaces with different lattice orientations, resulting also in highly mobile and highly coherent interfaces with shear-coupled migration.

### 3.5. Transformation-induced plasticity

The relevant applications of zirconia-based systems as structural materials rely on the martensitic transformation from the tetragonal to the monoclinic phase, which depends on the local stress and strain state. In this section, we test the DP-ZrO₂ model on the nanocompression of a thick tetragonal zirconia pillar along the [101], direction by a planar indenter. The three-dimensional system with size 15×15×15 nm³ is represented in Fig. 11 and contains approximately 320,000 atoms.

Periodic boundary conditions (PBC) are applied perpendicular to the compression axis (along x and y directions). The positions of the atoms in the bottom layers are fixed in the [101], direction of compression (z direction), but are allowed to relax parallel to the surface plane. The compression is performed at 0 K by monitoring the position of the planar indenter through a series of steps of 0.1 Å (0.07% engineering strain compression) until a compressive strain of 10%. Here, we define the strain as the vertical displacement of indenter over the height of the pillar. The indenter is represented by a purely repulsive potential acting on the atoms that have a vertical z-coordinate larger than the fixed indenter position. Static calculations are performed to avoid dynamical effects related to the large strain rates applied in MD simulations. At each compression step, the atomic positions and the cell dimensions in the surface plane (x, y) are relaxed to ensure $|σ_{xx}| ≈ |σ_{yy}| ≤ 5$ MPa. Technically, the simulation proceeds through a series of FIRE [73] and conjugate gradient minimizations, consisting of 1000 steps each, until a maximum atomic force of $2.5 × 10⁻⁴$ eV/Å is reached and the stress convergence criterion is simultaneously reached. All calculations were performed with the LAMMPS code [40] including the DeePMD-kit plugin.

The t-ZrO₂ phase is deformed elastically below 4.6% compression. At 4.6% compressive strain, the system undergoes a tetragonal to monoclinic phase transition, keeping its translational symmetry along the y direction (initial tetragonal b-axis [010],). During the transformation, the pillar shrinks in the vertical direction and detaches from the indenter. From 4.6% to 6.6% compression, the indenter moves down but remains separated from the pillar. The stress is thus zero. At 6.6% compression, the pillar reforms a contact with the indenter and after 6.6% compression, the monoclinic pillar is deformed elastically. We show in Fig. 12 five 2D snapshots seen along b-axis that illustrate various stages of the t-ZrO₂ → m-ZrO₂ transformation at 4.6% compression.

Fig. 12(a) shows the initial configuration of the tetragonal structure at $ε_{zz} = ε_{[101]} = 0\%$. The green colour indicates atoms in a tetragonal environment, while the atoms in grey at the bottom and close to the indenter have a local environment that deviates from the crystalline bulk due to the free surfaces. Fig. 12(e) displays the atomic configuration after 4.6% compression, where the martensitic transformation has taken place. The red colour signifies a complete t-ZrO₂ → m-ZrO₂ transformation towards a C-2 variant characterized by the orientation
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Simulation of nanocompression along the [101] direction. (a) Simulation setup. (b) Stress–strain curve. The stress is zero in the strain range between 4.6% and 6.6%, indicating that the indenter is no longer in contact with the pillar as a result of the phase transformations.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Atomic snapshots of the stress-induced t-ZrO₂ → m-ZrO₂ transformation during [101]ₜ nanocompression. (a) the relaxed t-ZrO₂ at zero strain; (b)–(e) are snapshots taken in a single relaxation at 4.6% compression. For visualization purposes, only Zr atoms are shown. Red, green, and blue atoms represent the m-ZrO₂, t-ZrO₂, and Pbc2₁ phases, respectively, while yellow atoms belong to a t-ZrO₂ phase that has undergone ferroelastic switching. Grey atoms have an unidentified environment.</figcaption></figure>

relationship [001]ₘ || [001]ₜ, [010]ₘ || [010]ₜ, (100)ₘ || (100)ₜ. The occurrence of this variant can be rationalized by considering the plastic strain produced by each of the 12 possible variants during compression along the [101]ₜ axis.

Analysis of the Bain transformation matrix reveals that only four variants, i.e. the C-1, C-2 and two A-type variants, exhibit a compressive strain along the [101]ₜ direction. They all produce a strain of approximately ~7%. If formed, they lead to stress-relaxation and are thus favoured by the applied strain. All the other variants show positive strains in the [101]ₜ direction, which correspond to a length increase

in the direction of compression. They are therefore not likely to appear here.

Further information on the nature of the possible variants can be inferred from the G-SSNEB results presented in Section 3.2, which show that, without an interface as in previous section, there is no direct path from the tetragonal structure to A or B monoclinic variants. Moreover the A variants have the largest coherency strain among all variants. C-1 and C-2 are therefore expected to be the favourite candidate orientations for this t-ZrO₂ → m-ZrO₂ transformation, in good agreement with the final C-2 structure obtained by the DP-ZrO₂ model.
Apart from the analysis of the initial and final structures, the details of the t-ZrO₂ → m-ZrO₂ transformation in Figs. 12(b–d) show other interesting features concerning the intermediate structures obtained in the course of the minimization process. The yellow colour in Fig. 12(c) corresponds to a ferroelastic switch ($a_t, b_t, c_t$) → ($b_t, c_t, a_t$) with the c tetragonal axis rotating perpendicular to the image plane. This ferroelastic switch initiates in the upper-right part of Fig. 12(c) in the vicinity of both the surface and the t-ZrO₂/m-ZrO₂ interface. This transformation produces a small negative strain and its localization suggests that it can be triggered by local defects.

Fig. 12(d) shows that the initial ferroelastic switch is followed by a t-ZrO₂ → Pbc₂₁ transformation. The orientation relationship is ($b_t, c_t, a_t$) → ($a_o, b_o, c_o$) with the long direction $b_o$ perpendicular to the image. From this Pbc₂₁ phase, a transition to the surrounding C-2 monoclinic variant is then available by following the same path as in Fig. 7 (purple curve). This last Pbc₂₁ → m-ZrO₂ transformation is characterized by the lattice correspondence ($a_o, b_o, c_o$) → ($a_m, b_m, c_m$). It is worth noting that the Pbc₂₁ phase has been experimentally observed under specific conditions in pure zirconia, such as growth on substrates [74], in nanoparticles [75], or in doped zirconia [17]; it has also been proposed from DFT calculations [60] as an intermediate structure in the formation of the monoclinic B variant, which is precisely what is described here through the successive transformations: ($b_t, c_t, a_t$) → ($a_o, b_o, c_o$) → ($a_m, b_m, c_m$).

The DP-ZrO₂ model thus provides valuable information on the mechanical properties of zirconia, showing that an intermediate phase can participate in the tetragonal to monoclinic transformation in pure zirconia. The results given by this machine learning potential are also consistent with both DFT results and elementary structural analysis. However, a distinct advantage of the DP-ZrO₂ model is its computational cost, which makes it suitable to study large systems of 320,000 atoms and more, far beyond the limitations of standard DFT calculations. At this system size, the numerical efficiency can be compared to a classical Born–Mayer–Buckingham potential [76–78] commonly used for zirconia systems. In this case, the average time per iteration is only twice larger for the DP-ZrO₂ model when running the same simulation on two graphic processing units (Nvidia Tesla V100 32 GB).

## 4. Conclusion

In this study, we developed a neural network interatomic potential to simulate temperature-induced and stress-induced phase transformations in ZrO₂. The accuracy of the potential was validated on both zero-temperature properties and the phase diagram of various ZrO₂ phases with temperatures up to 3700 K and pressures up to 20 GPa. The potential exhibits a high level of agreement with DFT calculations and experimental data. Using the novel neural network potential, the experimentally observed t-ZrO₂/m-ZrO₂ interface was successfully modelled, revealing a shear-coupled migration mode, which aligns well with the principles of martensite crystallography. The nanocompression of a t-ZrO₂ system showing transformation-induced plasticity involving multiple transformation paths was also simulated. Our study presents a precise and efficient approach for modelling microstructure evolution involving phase transformations in ZrO₂. This potential will be central to perform large-scale simulations of plasticity in ZrO₂ in order to understand the physical conditions favouring transformation-induced plasticity in single as well as polycrystalline ZrO₂ ceramics.

### CRediT authorship contribution statement

**Jin-Yu Zhang:** Interatomic Potential Development, Investigation, Visualization, Writing – original draft. **Gaël Huynh:** Investigation, Visualization, Writing – original draft. **Fu-Zhi Dai:** Methodology, Resources. **Tristan Albaret:** Investigation. **Shi-Hao Zhang:** Methodology. **Shigenobu Ogata:** Methodology. **David Rodney:** Conceptualization, Supervision, Result analysis, Funding acquisition.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

We would like to thank Dr. Xiaoyang Wang at Institute of Applied Physics and Computational Mathematics (China) for many useful discussions on phase diagram calculations. D. Rodney, T. Albaret & G. Huynh acknowledge funding from the ANR NANOTRIP (ANR-21-CE08-0019) and computational resources from the Jean-Zay machine at the Grand Équipement National de Calcul Intensif (GENCI) (Grant No. A0130810637). The support from National Key Research and Development Program of China (No. 2022YFB3709000) is also acknowledged. The support from the CNRS/IN2P3 Computing Center (Lyon - France) is gratefully acknowledged for providing computing and data-processing resources needed for this work. S.H.Z. and S.O. were supported by the JSPS Postdoctoral Fellowships for Research in Japan (Standard), the Grant-in-Aid for JSPS Research Fellow (Grant No. 22F22056) and the JSPS KAKENHI, Japan (Grant No. JP22KF-241), and used computational resources of supercomputer Fugaku provided by the RIKEN Center for Computational Science (Project IDs: hp230205 and hp230212), the large-scale computer systems at the Cybermedia Center, Osaka University, and the large-scale parallel computing server at the Center for Computational Materials Science, Institute for Materials Research, Tohoku University. S.O. acknowledges the support by the Ministry of Education, Culture, Sport, Science and Technology of Japan (Grant Nos. JPMXP1122684766, JPMXP1020230325, and JPMXP1020230327), and the support by JSPS KAKENHI (Grant Nos. JP23H00161 and JP23K20037).

## Appendix A. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.jeurceramsoc.2024.01.007.

## References

[1] J. Chevalier, L. Gremillard, A.V. Virkar, D.R. Clarke, The tetragonal-monoclinic transformation in zirconia: Lessons learned and future trends, J. Am. Ceram. Soc. 92 (9) (2009) 1901–1920, http://dx.doi.org/10.1111/j.1551-2916.2009.03278.x.

[2] Y.-W. Chen, J. Moussi, J.L. Drury, J.C. Wataha, Zirconia in biomedical applications, Expert Rev. Med. Devices 13 (10) (2016) 945–963, http://dx.doi.org/10.1080/17434440.2016.1230017.

[3] J. Chevalier, A. Liens, H. Reveron, F. Zhang, P. Reynaud, T. Douillard, L. Preiss, V. Sergio, V. Lughi, M. Swain, N. Courtois, Forty years after the promise of <<ceramic steel?>>: Zirconia-based composites with a metal-like mechanical behavior, J. Am. Ceram. Soc. 103 (3) (2020) 1482–1513, http://dx.doi.org/10.1111/jace.16903.

[4] P.M. Kelly, L.R. Francis Rose, The martensitic transformation in ceramics — Its role in transformation toughening, Prog. Mater. Sci. 47 (5) (2002) 463–557, http://dx.doi.org/10.1016/S0079-6425(00)00005-0.

[5] T.K. Gupta, F.F. Lange, J.H. Bechtold, Effect of stress-induced phase transformation on the properties of polycrystalline zirconia containing metastable tetragonal phase, J. Mater. Sci. 13 (7) (1978) 1464–1470, http://dx.doi.org/10.1007/BF00553200.

[6] M. Imariouane, M. Saâdaoui, G. Denis, H. Reveron, J. Chevalier, Low-yttria doped zirconia: Bridging the gap between strong and tough ceramics, J. Eur. Ceram. Soc. 43 (11) (2023) 4906–4915, http://dx.doi.org/10.1016/j.jeurceramsoc.2023.04.021.

[7] M. Li, B. Tunca, B. Van Meerbeek, J. Vleugels, F. Zhang, Tough and damage-tolerant monolithic zirconia ceramics with transformation-induced plasticity by grain-boundary segregation, J. Eur. Ceram. Soc. 43 (5) (2023) 2078–2092, http://dx.doi.org/10.1016/j.jeurceramsoc.2022.11.069.

[8] E.L. Pang, C.A. McCandler, C.A. Schuh, Reduced cracking in polycrystalline ZrO₂-CeO₂ shape-memory ceramics by meeting the cofactor conditions, Acta Mater. 177 (2019) 230–239, http://dx.doi.org/10.1016/j.actamat.2019.07.028.

[9] E.L. Pang, G.B. Olson, C.A. Schuh, The mechanism of thermal transformation hysteresis in ZrO₂-CeO₂ shape-memory ceramics, Acta Mater. 213 (2021) 116972, http://dx.doi.org/10.1016/j.actamat.2021.116972.
[10] E.L. Pang, G.B. Olson, C.A. Schuh, Low-hysteresis shape-memory ceramics designed by multimode modelling, Nature 610 (7932) (2022) 491–495, http://dx.doi.org/10.1038/s41586-022-05210-1.

[11] N. Zhang, M. Asle Zaeem, Competing mechanisms between dislocation and phase transformation in plastic deformation of single crystalline yttria-stabilized tetragonal zirconia nanopillars, Acta Mater. 120 (2016) 337–347, http://dx.doi.org/10.1016/j.actamat.2016.08.075.

[12] J. Zhang, M. Zhang, L. Deng, J. Jin, P. Gong, X. Wang, Mechanical behavior of tetragonal zirconia nanopillars subjected to uniaxial loading: A molecular dynamics study, Mech. Mater. 151 (2020) 103666, http://dx.doi.org/10.1016/j.mechmat.2020.103666.

[13] J. Zhou, Z. Jiao, J. Zhang, Z. Zhong, Nanoindentation of single-crystal and polycrystalline yttria-stabilized zirconia: A comparative study by experiments and molecular dynamics simulations, J. Alloys Compd. 878 (2021) 160336, http://dx.doi.org/10.1016/j.jallcom.2021.160336.

[14] N. Zhang, M. Asle Zaeem, Effects of twin boundaries and pre-existing defects on mechanical properties and deformation mechanisms of yttria-stabilized tetragonal zirconia, J. Eur. Ceram. Soc. 40 (1) (2020) 108–114, http://dx.doi.org/10.1016/j.jeurceramsoc.2019.09.017.

[15] P.K. Schelling, S.R. Phillpot, D. Wolf, Mechanism of the cubic-to-tetragonal phase transition in zirconia and yttria-stabilized zirconia by molecular-dynamics simulation, J. Am. Ceram. Soc. 84 (7) (2001) 1609–1619.

[16] J. Yu, R. Devanathan, W.J. Weber, Unified interatomic potential for zircon, zirconia and silica systems, J. Mater. Chem. 19 (23) (2009) 3923–3930.

[17] E.H. Kisi, C.J. Howard, R.J. Hill, Crystal structure of orthorhombic zirconia in partially stabilized zirconia, J. Am. Ceram. Soc. 72 (9) (1989) 1757–1760, http://dx.doi.org/10.1111/j.1151-2916.1989.tb06322.x.

[18] D.B. Marshall, M.R. Jarnes, J.R. Porter, Structural and mechanical property changes in toughened magnesia-partially-stabilized zirconia at low temperatures, J. Am. Ceram. Soc. 72 (2) (1989) 218–227, http://dx.doi.org/10.1111/j.1151-2916.1989.tb06104.x.

[19] S. Liu, W. Hu, Y. Zhang, J. Xiang, F. Wen, B. Xu, J. He, D. Yu, Y. Tian, Z. Liu, Metastable adaptive orthorhombic martensite in zirconia nanoparticles, J. Appl. Crystallogr. 47 (2) (2014) 684–691, http://dx.doi.org/10.1107/S1600576714003331.

[20] R. Ganser, S. Bongarz, A. von Mach, L. Azevedo Antunes, A. Kersch, Piezo- and pyroelectricity in zirconia: A study with machine-learned force fields, Phys. Rev. A 18 (5) (2022) 054066, http://dx.doi.org/10.1103/PhysRevApplied.18.054066.

[21] C. Verdi, F. Karsai, P. Liu, R. Jinnouchi, G. Kresse, Thermal transport and phase transitions of zirconia by on-the-fly machine-learned interatomic potentials, npj Comput. Mater. 7 (1) (2021) 1–9, http://dx.doi.org/10.1038/s41524-021-00630-5.

[22] P. Liu, C. Verdi, F. Karsai, G. Kresse, Phase transitions of zirconia: Machine-learned force fields beyond density functional theory, Phys. Rev. B 105 (6) (2022) L060102.

[23] C. Wang, A. Tharval, J.R. Kitchin, A density functional theory parameterised neural network model of zirconia, Mol. Simul. 44 (8) (2018) 623–630.

[24] S.-H. Guan, C. Shang, Z.-P. Liu, Resolving the temperature and composition dependence of ion conductivity for yttria-stabilized zirconia from machine learning simulation, J. Phys. Chem. C 124 (28) (2020) 15085–15093.

[25] F.-Z. Dai, B. Wen, H. Xiang, Y. Zhou, Grain boundary strengthening in ZrB₂ by segregation of W: Atomistic simulations with deep learning potential, J. Eur. Ceram. Soc. 40 (15) (2020) 5029–5036, http://dx.doi.org/10.1016/j.jeurceramsoc.2020.06.007.

[26] F. Fu, X. Wang, L. Zhang, Y. Yang, J. Chen, B. Xu, C. Ouyang, S. Xu, F.-Z. Dai, W. E., Unraveling the atomic-scale mechanism of phase transformations and structural evolutions during (de)lithiation in Si anodes, Adv. Funct. Mater. 33 (37) (2023) 2303936, http://dx.doi.org/10.1002/adfm.202303936.

[27] L. Zhang, H. Wang, R. Car, W. E., Phase diagram of a deep potential water model, Phys. Rev. Lett. 126 (23) (2021) 236001, http://dx.doi.org/10.1103/PhysRevLett.126.236001.

[28] H. Niu, L. Bonati, P.M. Piaggi, M. Parrinello, Ab initio phase diagram and nucleation of gallium, Nature Commun. 11 (1) (2020) 2654, http://dx.doi.org/10.1038/s41467-020-16372-9.

[29] D. Lu, W. Jiang, Y. Chen, L. Zhang, W. Jia, H. Wang, M. Chen, DP compress: A model compression scheme for generating efficient deep potential models, J. Chem. Theory Comput. 18 (9) (2022) 5559–5567, http://dx.doi.org/10.1021/acs.jctc.2c00102.

[30] D. Lu, H. Wang, M. Chen, L. Lin, R. Car, W. E, W. Jia, L. Zhang, 86 PFLOPS deep potential molecular dynamics simulation of 100 million atoms with ab initio accuracy, Comput. Phys. Comm. 259 (2021) 107624, http://dx.doi.org/10.1016/j.cpc.2020.107624.

[31] J. Li, Q. An, Quasiplastic deformation in shocked nanocrystalline boron carbide: Grain boundary sliding and local amorphization, J. Eur. Ceram. Soc. 43 (2) (2023) 208–216, http://dx.doi.org/10.1016/j.jeurceramsoc.2022.10.014.

[32] L. Zhang, J. Han, H. Wang, R. Car, W. E., Deep potential molecular dynamics: A scalable model with the accuracy of quantum mechanics, Phys. Rev. Lett. 120 (14) (2018) 143001, http://dx.doi.org/10.1103/PhysRevLett.120.143001.

[33] AIS square, https://www.aissquare.com/.

[34] L. Zhang, J. Han, H. Wang, W. Saidi, R. Car, W. E., End-to-end symmetry preserving inter-atomic potential energy model for finite and extended systems, in: Advances in Neural Information Processing Systems, vol. 31, Curran Associates, Inc., 2018, pp. 4441–4451.

[35] H. Wang, L. Zhang, J. Han, W. E., DeePMD-Kit: A deep learning package for many-body potential energy representation and molecular dynamics, Comput. Phys. Comm. 228 (2018) 178–184, http://dx.doi.org/10.1016/j.cpc.2018.03.016.

[36] J. Zeng, D. Zhang, D. Lu, P. Mo, Z. Li, Y. Chen, M. Rynik, L. Huang, Z. Li, S. Shi, Y. Wang, H. Ye, P. Tuo, J. Yang, Y. Ding, Y. Li, D. Tisi, Q. Zeng, H. Bao, Y. Xia, J. Huang, K. Muraoka, Y. Wang, J. Chang, F. Yuan, S.L. Bore, C. Cai, Y. Lin, B. Wang, J. Xu, J.-X. Zhu, C. Luo, Y. Zhang, R.E.A. Goodall, W. Liang, A.K. Singh, S. Yao, J. Zhang, R. Wentzcovitch, J. Han, J. Liu, W. Jia, D.M. York, W. E, R. Car, L. Zhang, H. Wang, DeePMD-kit v2: A software package for deep potential models, J. Chem. Phys. 159 (5) (2023) 054801, http://dx.doi.org/10.1063/5.0155600.

[37] T. Wen, L. Zhang, H. Wang, W. E, D.J. Srolovitz, Deep potentials for materials science, Mater. Futures 1 (2) (2022) 022601, http://dx.doi.org/10.1088/2752-5724/ac681d.

[38] L. Zhang, D.-Y. Lin, H. Wang, R. Car, W. E., Active learning of uniformly accurate interatomic potentials for materials simulation, Phys. Rev. Mater. 3 (2) (2019) 023804, http://dx.doi.org/10.1103/PhysRevMaterials.3.023804.

[39] D. Sheppard, P. Xiao, W. Chemelewski, D.D. Johnson, G. Henkelman, A generalized solid-state nudged elastic band method, J. Chem. Phys. 136 (7) (2012) 074103, http://dx.doi.org/10.1063/1.3684549.

[40] A.P. Thompson, H.M. Aktulga, R. Berger, D.S. Bolintineanu, W.M. Brown, P.S. Crozier, P.J. in 't Veld, A. Kohlmeyer, S.G. Moore, T.D. Nguyen, R. Shan, M.J. Stevens, J. Tranchida, C. Trott, S.J. Plimpton, LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Comm. 271 (2022) 108171, http://dx.doi.org/10.1016/j.cpc.2021.108171.

[41] G. Kresse, J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci. 6 (1) (1996) 15–50, http://dx.doi.org/10.1016/0927-0256(96)00008-0.

[42] G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54 (16) (1996) 11169–11186, http://dx.doi.org/10.1103/PhysRevB.54.11169.

[43] J.P. Perdew, A. Ruzsinszky, G.I. Csonka, O.A. Vydrov, G.E. Scuseria, L.A. Constantin, X. Zhou, K. Burke, Restoring the density-gradient expansion for exchange in solids and surfaces, Phys. Rev. Lett. 100 (13) (2008) 136406, http://dx.doi.org/10.1103/PhysRevLett.100.136406.

[44] W. Mayr-Schmölzer, J. Planer, J. Redinger, A. Grüneis, F. Mittendorfer, Many-electron calculations of the phase stability of ZrO₂ polymorphs, Phys. Rev. Res. 2 (4) (2020) 043361, http://dx.doi.org/10.1103/PhysRevResearch.2.043361.

[45] J. Sun, A. Ruzsinszky, J.P. Perdew, Strongly constrained and appropriately normed semilocal density functional, Phys. Rev. Lett. 115 (3) (2015) 036402, http://dx.doi.org/10.1103/PhysRevLett.115.036402.

[46] J. Wu, Y. Zhang, L. Zhang, S. Liu, Deep learning of accurate force field of ferroelectric HfO₂, Phys. Rev. B 103 (2) (2021) 024108, http://dx.doi.org/10.1103/PhysRevB.103.024108.

[47] P. Bouvier, E. Djurado, G. Lucazeau, T. Le Bihan, High-pressure structural evolution of undoped tetragonal nanocrystalline zirconia, Phys. Rev. B 62 (13) (2000) 8731–8737, http://dx.doi.org/10.1103/PhysRevB.62.8731.

[48] O. Ohtaka, H. Fukui, T. Kunisada, T. Fujisawa, K. Funakoshi, W. Utsumi, T. Irifune, K. Kuroda, T. Kikegawa, Phase relations and equations of state of ZrO₂ under high temperature and high pressure, Phys. Rev. B 63 (17) (2001) 174108, http://dx.doi.org/10.1103/PhysRevB.63.174108.

[49] A. Togo, I. Tanaka, First principles phonon calculations in materials science, Scr. Mater. 108 (2015) 1–5, http://dx.doi.org/10.1016/j scriptamat.2015.07.021.

[50] A. Kuwabara, T. Tohei, T. Yamamoto, I. Tanaka, Ab initio lattice dynamics and phase transformations of ZrO₂, Phys. Rev. B 71 (6) (2005) 064301.

[51] P. Aldebert, J.-P. Traverse, Structure and ionic mobility of zirconia at high temperature, J. Am. Ceram. Soc. 68 (1) (1985) 34–40, http://dx.doi.org/10.1111/j.1151-2916.1985.tb15247.x.

[52] R.N. Patil, E.C. Subbarao, Axial thermal expansion of ZrO₂ and HfO₂ in the range room temperature to 1400°C, J. Appl. Crystallogr. 2 (6) (1969) 281–288, http://dx.doi.org/10.1107/S0021889869007217.

[53] S.P. Terblanche, Thermal-expansion coefficients of yttria-stabilized cubic zirconias, J. Appl. Crystallogr. 22 (3) (1989) 283–284, http://dx.doi.org/10.1107/S0021889888013937.

[54] O. Valsson, M. Parrinello, Variational approach to enhanced sampling and free energy calculations, Phys. Rev. Lett. 113 (9) (2014) 090601, http://dx.doi.org/10.1103/PhysRevLett.113.090601.

[55] Q.-J. Hong, S.V. Ushakov, D. Kapush, C.J. Benmore, R.J.K. Weber, A. van de Walle, A. Navrotsky, Combined computational and experimental investigation of high temperature thermodynamics and structure of cubic ZrO₂ and HfO₂, Sci. Rep. 8 (1) (2018) 14962, http://dx.doi.org/10.1038/s41598-018-32847-8.

[56] A.V. Virkar, R.L.K. Matsumoto, Ferroelastic domain switching as a toughening mechanism in tetragonal zirconia, J. Am. Ceram. Soc. 69 (10) (1986) C-224-C-226, http://dx.doi.org/10.1111/j.1151-2916.1986.tb07341.x.
[57] C.-J. Chan, F.F. Lange, M. Rühle, J.-F. Jue, A.V. Virkar, Ferroelastic domain switching in tetragonal zirconia single crystals—Microstructural aspects, J. Am. Ceram. Soc. 74 (4) (1991) 807–813, http://dx.doi.org/10.1111/j.1151-2916.1991.tb06929.x.

[58] W. Kriven, W. Fraser, S. Kennedy, Martensite crystallography of tetragonal zirconia, in: Advances in Ceramics, Vol. 3, American Ceramic Society, 1981, pp. 82–97.

[59] C. Carbogno, C.G. Levi, C.G. Van De Walle, M. Scheffler, Ferroelastic switching of doped zirconia: Modeling and understanding from first principles, Phys. Rev. B 90 (14) (2014) 144109, http://dx.doi.org/10.1103/PhysRevB.90.144109.

[60] S.-H. Guan, X.-J. Zhang, Z.-P. Liu, Energy landscape of zirconia phase transitions, J. Am. Chem. Soc. 137 (25) (2015) 8010–8013.

[61] G. Trolliard, D. Mercurio, J.M. Perez-Mato, Martensitic phase transition in pure zirconia: a crystal chemistry viewpoint, Z. Kristallogr. 226 (3) (2011) 264–290, http://dx.doi.org/10.1524/zkri.2011.1340.

[62] C. Vega, E. Sanz, J.L.F. Abascal, E.G. Noya, Determination of phase diagrams via computer simulation: Methodology and applications to water, electrolytes and proteins, J. Phys.: Condens. Matter 20 (15) (2008) 153101, http://dx.doi.org/10.1088/0953-8984/20/15/153101.

[63] C. Wang, M. Zinkevich, F. Aldinger, The zirconia–hafnia system: DTA measurements and thermodynamic calculations, J. Am. Ceram. Soc. 89 (12) (2006) 3751–3758, http://dx.doi.org/10.1111/j.1551-2916.2006.01286.x.

[64] A. Navrotsky, L. Benoist, H. Lefebvre, Direct calorimetric measurement of enthalpies of phase transitions at 2000 °–2400 °C in yttria and zirconia, J. Am. Ceram. Soc. 88 (10) (2005) 2942–2944, http://dx.doi.org/10.1111/j.1551-2916.2005.00506.x.

[65] D.A. Kofke, Gibbs-Duhem integration: A new method for direct evaluation of phase coexistence by molecular simulation, Mol. Phys. 78 (6) (1993) 1331–1336, http://dx.doi.org/10.1080/00268979300100881.

[66] H. Peng, Z.-H. Yang, J.P. Perdew, J. Sun, Versatile van der waals density functional based on a meta-generalized gradient approximation, Phys. Rev. X 6 (4) (2016) 041005, http://dx.doi.org/10.1103/PhysRevX.6.041005.

[67] W. Wang, Z. Liang, X. Han, J. Chen, C. Xue, H. Zhao, Mechanical and thermodynamic properties of ZrO₂ under high-pressure phase transition: A first-principles study, J. Alloys Compd. 622 (2015) 504–512, http://dx.doi.org/10.1016/j.jallcom.2014.08.114.

[68] S. Block, J. a. H. Da Jornada, G.J. Piermarini, Pressure-temperature phase diagram of zirconia, J. Am. Ceram. Soc. 68 (9) (1985) 497–499, http://dx.doi.org/10.1111/j.1151-2916.1985.tb15817.x.

[69] K. Matsui, K. Nakamura, M. Saito, A. Kuwabara, H. Yoshida, Y. Ikuhara, Low-temperature degradation in yttria-stabilized tetragonal zirconia polycrystal: Effect of Y³⁺ distribution in grain interiors, Acta Mater. 227 (2022) 117659, http://dx.doi.org/10.1016/j.actamat.2022.117659.

[70] M.S. Wechsler, D.S. Lieberman, T.A. Read, On the theory of the formation of martensite, Trans. AIME 197 (11) (1953) 1503–1515.

[71] J.K. Mackenzie, J.S. Bowles, The crystallography of martensite transformations II, Acta Metall. 2 (1) (1954) 138–147, http://dx.doi.org/10.1016/0001-6160(54)90103-0.

[72] K. Luo, B. Liu, W. Hu, X. Dong, Y. Wang, Q. Huang, Y. Gao, L. Sun, Z. Zhao, Y. Wu, Y. Zhang, M. Ma, X.-F. Zhou, J. He, D. Yu, Z. Liu, B. Xu, Y. Tian, Coherent interfaces govern direct transformation from graphite to diamond, Nature 607 (7919) (2022) 486–491, http://dx.doi.org/10.1038/s41586-022-04863-2.

[73] J. Guénolé, W.G. Nöhring, A. Vaid, F. Houllé, Z. Xie, A. Prakash, E. Bitzek, Assessment and optimization of the fast inertial relaxation engine (fire) for energy minimization in atomistic simulations and its implementation in lammps, Comput. Mater. Sci. 175 (2020) 109584, http://dx.doi.org/10.1016/j.commatsci.2020.109584.

[74] G. Trolliard, R. Benmechta, D. Mercurio, Pure orthorhombic zirconia islands grown on single-crystal sapphire substrates, Acta Mater. 55 (17) (2007) 6011–6018, http://dx.doi.org/10.1016/j.actamat.2007.07.009.

[75] S. Liu, W. Hu, Y. Zhang, J. Xiang, F. Wen, B. Xu, J. He, D. Yu, Y. Tian, Z. Liu, Metastable adaptive orthorhombic martensite in zirconia nanoparticles, J. Appl. Crystallogr. 47 (2) (2014) 684–691, http://dx.doi.org/10.1107/S1600576714003331.

[76] X. Li, B. Hafskjold, Molecular dynamics simulations of yttrium-stabilized zirconia, J. Phys.: Condens. Matter 7 (7) (1995) 1255–1271, http://dx.doi.org/10.1088/0953-8984/7/7/007.

[77] P. Koziatek, J. Barrat, D. Rodney, Short-and medium-range orders in as-quenched and deformed SiO₂ glasses: An atomistic study, J. Non-Crystall. Solids 414 (2015) 7–15.

[78] M. Asle Zaeem, N. Zhang, M. Mamivand, A review of computational modeling techniques in study and design of shape memory ceramics, Comput. Mater. Sci. 160 (2019) 120–136, http://dx.doi.org/10.1016/j.commatsci.2018.12.062.