# Transformation-induced plasticity in CeO 2 -ZrO 2 ceramics: Atomic-scale insights using a deep neural network potential

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

► To cite this version:

Jin-Yu Zhang, Gaël Huynh, Tristan Albaret, Fu-Zhi Dai, David Rodney. Transformation-induced plasticity in CeO 2 -ZrO 2 ceramics: Atomic-scale insights using a deep neural network potential. Acta Materialia, 2025, 285, pp.120661. 10.1016/j.actamat.2024.120661. hal-04935211

HAL Id: hal-04935211

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

Submitted on 7 Feb 2025

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>

Full length article

# Transformation-induced plasticity in CeO₂-ZrO₂ ceramics: Atomic-scale insights using a deep neural network potential

Jin-Yu Zhang<sup>a,b,∗</sup>, Gaël Huynh<sup>a,†</sup>, Tristan Albaret<sup>a</sup>, Fu-Zhi Dai<sup>c,d</sup>, David Rodney<sup>a,‡,∗∗,1</sup>

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

## ARTICLE INFO

### Keywords:

Zirconia

Machine learning potential

Phase transformation

Molecular dynamics

Plasticity

## ABSTRACT

Zirconia (ZrO₂) ceramics exhibit remarkable mechanical properties, namely transformation-induced plasticity (TRIP), shape memory, and superelasticity. However, understanding the complex atomic-scale processes controlling these phenomena is challenging. Here, we introduce a deep neural network interatomic potential to accurately predict phase transformations in both pure and CeO₂-doped ZrO₂ ceramics. Through molecular dynamics simulations, we examine the mechanical responses of tetragonal CeO₂-ZrO₂ single crystals and polycrystals under uniaxial loading. At odds with the traditional understanding of the TRIP effect in zirconia ceramics, our simulations reveal that the classical stress-induced tetragonal-to-monoclinic phase transformation often involves an intermediate phase, which can either be a tetragonal phase resulting from ferroelastic switching or an orthorhombic Pbc2₁ phase. Consequently, the TRIP effect may be hindered if the applied compression is unfavorable for the intermediate phase, even if the final martensite phase has a high Schmid factor. Our findings on polycrystalline compression underscore the importance of grain boundaries as nucleation sites. Combined with the complex internal stress distribution, this leads to the formation of all three monoclinic lattice correspondences through complex pathways that are analyzed in detail.

## 1. Introduction

Brittleness significantly limits the widespread use of ceramics in structural applications, despite their desirable properties of hardness and corrosion resistance. Zirconia (ZrO₂) ceramics represent a notable exception due to transformation-induced plasticity (TRIP), a phenomenon that enables substantial plastic strain (up to 7% in ceria (CeO₂) doped ZrO₂ [1]), which largely enhances ductility [2–5]. This phase transformation, from a tetragonal (space group P4₂/nmc) to a monoclinic (space group P2₁/c) phase, is characterized by three potential lattice correspondences (A, B, and C), each with two orientation relationships (ORs) as shown in Fig. 1. This figure is based on the definition established by Kelly and Rose [6], with the monoclinic unit cell viewed along the [010]ₘ direction for clearer visualization of the ORs, habit planes, and the macroscopic shear directions. For ORs C-1, B-1, and A-1, where the *a*ₘ and *b*ₘ axes are parallel to the tetragonal axes, the habit plane is close to (001)ₘ, and the macroscopic shear

vector close to [100]ₘ [6]. Meanwhile, for ORs C-2, B-2, and A-2, where the *c*ₘ and *b*ₘ axes are parallel to tetragonal axes, the habit plane is close to (100)ₘ, with a macroscopic shear vector near [001]ₘ. Previous theoretical and density functional theory (DFT) studies suggest that correspondence C is favored due to its minimal Bain strain [6] and a direct transformation path with a single energy barrier (green path in Fig. 2) [7], facilitating the nucleation of the martensitic transformation. In contrast, correspondence B, which has an indirect transformation path through either a ferroelastic switching of the tetragonal phase (succession of the blue and green lines in Fig. 2) or an orthorhombic phase with space group Pbc2₁ (succession of the red and purple lines in Fig. 2) [7,8], exhibits the smallest middle eigenvalue of the transformation stretch tensor, λ₂ [6]. This property is important as it is closely related to interface compatibility and transformation hysteresis [9]. The Bain strain of the tetragonal-to-monoclinic transformation contains both a volume expansion and a shear close to a ⟨100⟩ᵢ{001}ᵢ slip system. However, shear dominates, allowing researchers to analyze the stress

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

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

<footnote><sup>c</sup> School of Materials Science and Engineering, University of Science and Technology Beijing, Beijing, 100083, People's Republic of China</footnote>

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

<footnote>* Corresponding author at: Institut Lumière Matière, Université Claude Bernard Lyon 1, Villeurbanne, 69622, France.</footnote>

<footnote>** Corresponding author.</footnote>

<footnote>E-mail addresses: zhangjinyuchengqi@gmail.com (J.-Y. Zhang), david.rodney@univ-lyon1.fr (D. Rodney).</footnote>

<footnote n="1">¹ David Rodney was an Editor of the journal during the review period of the article. To avoid a conflict of interest, David Rodney was blinded to the record and another editor processed this manuscript.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. The illustration of the tetragonal-to-monoclinic phase transformation (the tetragonal cell contains 12 atoms, i.e., a double cell), showcasing the three distinct lattice correspondences (A, B, and C), together with both corresponding ORs.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. Minimum energy paths predicted in pure ZrO₂ by DFT, reproduced from Ref. [7]. *a*, *b*, and *c* denote the axes in Pbc2₁ phase (orthorhombic phases in this study are named according to their space groups). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

effect using Schmid factors [10]. The latter carries less information than the shape strain and, for instance, does not allow to distinguish between A and C correspondences [11]. However, in the following, for a simple and straightforward analysis, we will use Schmid factors.

The factors that control the selection of these different lattice corre-
spondences are central to understanding the TRIP effect of ZrO₂ ceram-
ics. Prior studies present different experimental observations regarding
the correspondences. Zeng et al. [12] investigated the crystallogra-
phy of both thermally- and stress-induced transformations in micron-
sized zirconia single crystals, identifying correspondence B in thermal
transformations and correspondence C in stress-induced transforma-
tions, respectively. Moreover, micropillar compression tests across var-
ious loading directions generally support correspondence C for stress-
induced transformations, with near [101]ᵢ orientation exhibiting the

lowest critical stress and highest Schmid factor associated with corre-
spondence C [10]. In contrast, dislocation slip instead of phase trans-
formation is favored near [110]ᵢ, which exhibits a high Schmid factor
for correspondence B, although exceptions involving phase transfor-
mations exist [10]. The situation is further complicated in polycrys-
talline ZrO₂, where both correspondences B and C were reported in
thermally-induced transformations [13–18], with B being more promi-
nent [13,14]. Similar findings were reported for stress-induced trans-
formations [13,14]. Correspondence A was typically observed only in
conjunction with other correspondences, suggesting a minor role in the
TRIP effect [16,17].

Molecular dynamics offers a valuable alternative to explore mi-
crostructural evolutions during the martensitic transformation in single-
[19–21], bi- [22] and poly-crystalline zirconia [21,23]. However, the
empirical potentials used in previous studies often struggle to accurately capture the stability and energetics of the ZrO₂ polymorphs observed experimentally, as pointed out in Ref. [24]. Although the latest empirical potentials (GAOP1 and GAOP2 in Ref. [24]) successfully reproduce the stability order of ZrO₂ phases, they exhibit rather large energy deviations from DFT calculations, such as approximately 20 meV/atom for the cubic phase. Moreover, the Young's moduli reported in previous simulations are significantly higher than experimental values [10] that are typically less than 300 GPa. For example, the apparent Young's moduli along the [001] direction are approximately 600 GPa in Fig. 8 of Ref. [20] and several thousand GPa in Fig. 2 of Ref. [19]. Moreover, Yin et al. [25] systematically evaluated various empirical potentials for cubic yttria-stabilized zirconia, including the potentials used in Refs. [19,20]. They found that the Young moduli exhibit relative errors exceeding 50% compared to experiments. Furthermore, previous simulations lack a detailed crystallographic analysis of the TRIP effect, particularly in terms of lattice correspondences during the tetragonal-to-monoclinic transformation and the potential role of ferroelastic switching. These factors are crucial for understanding the activation of the TRIP mechanism under different loading conditions.

In contrast, machine learning interatomic potentials can accurately reproduce the fundamental properties of zirconia, typically with an energetic deviation of less than 2 meV/atom from DFT results [7]. Various machine learning interatomic potentials have been developed for pure ZrO₂, effectively describing thermally-induced phase transformations [26–29], as well as transformations driven by both thermal and pressure changes [7]. However, machine learning interatomic potentials often face challenges in large-scale atomistic simulations due to their high computational cost. To overcome this, deep potentials (DP) utilize a compression technique for neural networks [30], enhancing efficiency by an order of magnitude and enabling simulations of systems with millions of atoms [31–34]. Also, we note that most existing machine learning interatomic potentials primarily focus on pure ZrO₂ and neglect dopants like CeO₂ that are crucial to stabilize the tetragonal phase at ambient temperature. An exception is the potential developed by Peng et al. [35], but the latter is mainly limited to applications in the cubic phase at high Ce composition (>50%, with % indicating mol% hereafter when referring to Ce composition) at ambient pressure. Very recently, Rettenmaier et al. [36] developed a CeO₂-ZrO₂ deep potential, which has a good prediction of the tetragonal-to-monoclinic transformation temperatures while the predicted λ₂ value as a function of CeO₂ composition has an opposite trend to experiments. Moreover, neither potential accounts for orthorhombic phases, e.g., Pbc2₁, which may influence the transformation paths from the tetragonal phase [37–39] and may reduce transformation toughening [6,37,38].

In the present work, we developed an interatomic potential for ZrO₂ ceramics doped with CeO₂ (Section 2) using a machine learning-based deep potential (DP) model [40], which has resolved challenges across various material systems [41–45]. The DP model is trained based on a DFT dataset constructed via an active learning scheme and further refined through a transfer learning technique. In Section 3, we evaluate the properties of CeO₂-doped ZrO₂ ceramics, such as elastic constants and phase diagrams, comparing them quantitatively against DFT calculations and/or experimental data. The DFT datasets and corresponding DP models are accessible on a public website [46]. Employing the potential, we conduct large-scale simulations on the TRIP effect in single- and poly-crystalline CeO₂-ZrO₂ ceramics (Section 4) under uniaxial loading. These simulations aim to demonstrate the models' capability to reproduce key features such as the TRIP effect and superelasticity, and to explore the development of various lattice correspondences under deformation in CeO₂-ZrO₂ ceramics.

## 2. Computational methods

### 2.1. Setup of the deep potential model

The DP model, implemented in the DeePMD-kit package [47,48], is utilized to accurately represent atomic interactions within the CeO₂-ZrO₂ system. The DP model represents the potential energy (E) of the

system as a sum of individual atomic energies, expressed as $E = \sum_i E_i$. Each atomic energy, $E_i$, is a function of the local atomic environment, defined by $E_i = E(\mathbf{R}_i)$, where $\mathbf{R}_i = (\mathbf{r}_{i1}, \mathbf{r}_{i2}, \dots, \mathbf{r}_{im})$ represents the matrix of the relative coordinates of the neighboring atoms of atom i, confined within a cutoff radius $r_c$. This mapping from $\mathbf{R}_i$ to $E_i$ is achieved through a two-tier neural network architecture comprising an embedding neural network followed by a fitting neural network. Initially, the embedding network transforms $\mathbf{R}_i$ into a feature matrix (descriptor) $\mathbf{D}_i$, ensuring preservation of the system's translational, rotational, and permutational symmetries [49]. Subsequently, $\mathbf{D}_i$ is processed through the fitting network to obtain $E_i$.

For the present study, we employed the smooth version of the Deep-Pot descriptor [49]. The cutoff radius was set to 7.5 Å, incorporating a smoothing function that extends from 2.0 Å to 7.5 Å. The embedding and fitting networks contain three hidden layers, with dimensions (25, 50, 100) and (240, 240, 240), respectively. The feature matrix has a size of 16 × 100. The parameter optimization within the embedding and fitting networks is guided by minimizing the loss function [47,48]:

$$L(p_c, p_f, p_\xi) = \frac{p_c}{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 denotes the number of atoms. Here, $\Delta\epsilon$, $|\Delta F_i|$, and $\|\Delta\xi\|$ represent respectively the discrepancies between the predictions of the DP model and DFT calculations for energies, atomic forces, and virial tensors. The corresponding prefactors $p_c$, $p_f$, and $p_\xi$ are dynamically adjusted in the training and fine-tuning processes. Initially, the system undergoes training over eight million steps with gradual adjustments of $p_c$, $p_f$, and $p_\xi$ from 0.02 to 1, 1000 to 1, and 0.02 to 1, respectively. Subsequently, the DP model undergoes a fine-tuning of eight million steps with fixed prefactors set to 10, 1, and 1 for $p_c$, $p_f$, and $p_\xi$ respectively. This fine-tuning stage is pivotal in enhancing energy accuracy while maintaining the precision of forces and stresses [50,51].

### 2.2. Active learning scheme for training DP-PBEsol

Fig. 3(a) (left side) depicts the workflow of an active learning scheme, employed to efficiently generate a DFT dataset for training the DP model. This scheme, integrated within the DP-GEN (Deep Generator) software [52], operates through an iterative “DPGEN Run” process comprising three stages: exploration, labeling, and training. During the training step, DP-GEN utilizes the DeePMD-kit package to train on a DFT dataset four distinct DP models obtained starting from different initial network parameters. To accelerate the “DPGEN Run” process, each iteration involves 800,000 training steps, omitting the fine-tuning phase. The initial training dataset incorporates structures of pure ZrO₂ and CeO₂-ZrO₂ solid solutions. Notably, DFT data from a prior study [7] on DP-ZrO₂, which demonstrated accuracy comparable to DFT for phase transformations, were included in the initial CeO₂-ZrO₂ training dataset. Besides that, various ZrO₂ structures comprising 96 atoms (monoclinic, tetragonal, cubic, Pbc2₁, Pbca, Pmna, and liquid phases) were generated by molecular dynamics utilizing the DP-ZrO₂ model. Subsequently, in these structures, 1 to 8, 10, 12, 16, 20, 25, and 32 Zr atoms were randomly substituted with Ce atoms to encompass the entire compositional spectrum. The low-index surfaces of monoclinic and tetragonal CeO₂-ZrO₂ phases were similarly generated and included. Additionally, the transition states of the tetragonal-to-monoclinic phase transformation and of ferroelastic switching were included in the initial training dataset based on generalized solid-state nudged elastic band (G-SSNEB) method [53] calculations in 96-atom cells with various Ce contents. Overall, the initial training dataset used to train the four initial CeO₂-ZrO₂ DP models comprised 1160 CeO₂-ZrO₂ and 5057 ZrO₂ structures.

During the exploration stage, DP-GEN uses LAMMPS [54] to conduct hybrid molecular dynamics/Monte Carlo simulations for CeO₂-ZrO₂ systems across various phases, compositions, pressures, and temperatures, thereby extensively probing the configurational space. The
phases and composition spectrum are consistent with those of CeO₂-ZrO₂ structures in the initial training dataset, i.e., 96-atom structures of all the experimentally observed zirconia phases, with random substitution of 1 to 8, 10, 12, 16, 20, 25, and 32 Zr atoms by Ce. Temperatures were progressively raised from 100 K to 4500 K over the course of the DP-GEN iterations, exploring pressures ranging from –5 to 20 GPa. Because zirconia exhibits a wide band gap (for example, 4.7 eV for the amorphous phase [55]), which corresponds to a thermal excitation temperature of approximately 5×10⁴ K), electron ground-state calculations using DFT are suitable for describing high-temperature environments. The low-index surfaces of CeO₂-ZrO₂ monoclinic and tetragonal phases below 500 K were also generated in this step. Representative initial configurations are shown in Fig. 3(b). To efficiently evaluate the predictive accuracy of the four DP models trained in the preceding step, we measured the prediction deviations among DP models, defined as the maximum standard deviation in predicted atomic forces across different DP models. Large model deviations highlight low predictive accuracy for the generated structures, underscoring the need for training dataset refinement. In the present study, the configurations with model deviations from 0.20 eV/Å to 0.40 eV/Å were randomly added into the training dataset. For liquid structures, the range was adjusted to 0.30–0.50 eV/Å. The upper limit of model deviations was set to avoid picking up potential unphysical configurations.

In the labeling phase, the energy, atomic forces, and stresses of the chosen configurations were computed using the Vienna ab initio simulation package (VASP) [56,57], employing the Perdew–Burke–Ernzerhof exchange–correlation functional revised for solids (PBEsol) [58], which is denoted as DFT-PBEsol below. The calculations were set with a plane-wave cutoff of 900 eV and a 3 × 3 × 3 Monkhorst–Pack k-mesh to ensure convergence of the energies, atomic forces, and virials for the distorted crystals and liquid configurations of CeO₂-ZrO₂ systems. To accelerate single point DFT calculations, a sequential two-step strategy was implemented. In the first step, wavefunctions of electrons were brought to a preliminary level of convergence with reduced precision, facilitated by enabling specific computational efficiency settings in VASP (see Supplementary Materials for details). Following this, the calculation was advanced to a high precision, commencing from the initially preconverged wavefunctions. This approach effectively reduces the overall computational time while maintaining accuracy of the results. A total of 12,303 configurations were incorporated into the DFT-PBEsol dataset. The DP model, as outlined in Section 2.1 and henceforth referred to as the DP-PBEsol model, is trained based on the DFT-PBEsol dataset. The validation of the DP-PBEsol model, exhibiting comparable accuracy to the DP-ZrO₂ model [7] and capturing the qualitative aspects of tetragonal phase stabilization by CeO₂, is detailed in the Supplementary Materials.

## 2.3. Transfer learning scheme for training DP-SCAN+rVV10

While the DP-PBEsol model demonstrates accuracy comparable to DFT-PBEsol and reasonable agreement with experimental results, it is important to note the presence of quantifiable deviations. Specifically, in pure ZrO₂, the phase transformation temperatures of m-ZrO₂/t-ZrO₂, t-ZrO₂/c-ZrO₂, and c-ZrO₂/liquid ZrO₂ (L-ZrO₂) under ambient pressure exhibit discrepancies of 200–400 K when compared to experimental data [59] (see also Table S4 in the Supplementary Materials). This is particularly pronounced for the t-ZrO₂ phase, which is experimentally established as the most stable between 1367–2584 K over a range of about 1200 K [59], whereas DP-PBEsol predicts stability within a narrower temperature range of about 700 K (1640–2360 K). Furthermore, the transformation pressure predicted for the high-pressure Pbca phase is 5.2 GPa (Table S4), notably higher than the experimentally determined 3.0 GPa [60]. There is also a divergence in the equilibrium CeO₂ composition at 0 K between the tetragonal and monoclinic phases, with DP-PBEsol predicting roughly 20% (Fig. S7(a)), compared to the experimental value of 16% [61].

Those discrepancies partly stem from the intrinsic accuracy of the PBEsol functional itself, as listed in Table S1, where DFT-PBEsol predicts a higher energy for the tetragonal phase (79.2 meV/f.u.) relative to the experimental value (63 meV/f.u.) [62,63]. In contrast, the DFT prediction (60.8 meV/f.u.) using SCAN+rVV10 [64], i.e., strongly constrained and appropriately normed (SCAN) meta-generalized gradient approximation [65] combined with the van der Waals interaction from rVV10 (the revised Vydrov–van Voorhis non-local correlation functional [66]), is better aligned with experimental data [67]. We refer to the DFT calculations using SCAN+rVV10 as DFT-SCAN+rVV10 hereinafter. Consequently, a broader temperature range for the tetragonal phase is anticipated when using SCAN+rVV10, and the equilibrium CeO₂ composition is expected to be reduced (Fig. S6), converging to the experimental observations. Furthermore, the lower energy calculated for the Pbca phase by DFT-SCAN+rVV10 (23.4 meV/f.u.) compared to DFT-PBEsol (40.9 meV/f.u.) in Table S1 suggests a stabilization of the Pbca phase at a lower pressure, potentially closer to the experimental value [60].

However, DFT-SCAN+rVV10 is considerably more computationally intensive compared to DFT-PBEsol, precluding its application for all structures in the DFT-PBEsol dataset. To address this challenge, a transfer learning approach was employed, as implemented in the “DP-GEN Simplify” module of the DPGEN package [52]. This technique was designed to reduce computational demands while preserving the accuracy of DFT-SCAN+rVV10. Fig. 3(a) (right side) illustrates the iterative process of “DPGEN Simplify”, comprising three primary steps: Simplifying, Relabeling, and Transfer Learning. The four DP models generated during the “DPGEN Run” serve as pre-trained models for transfer learning. In this process, the knowledge of DP models learned from DFT-PBEsol dataset is reused by fixing the embedding net, which encodes local atomic configurations into descriptors. Only the fitting network is fine-tuned using a small DFT-SCAN+rVV10 dataset. Initially, 500 structures were randomly selected from the DFT-PBEsol dataset to form the preliminary DFT-SCAN+rVV10 dataset for the initial fine-tuning. This fine-tuning employs the same prefactors as discussed in Section 2.2, but each iteration comprises only 800,000 steps to accelerate the process.

In the simplifying step, the model deviations of the fine-tuned DP models are assessed for structures in the DFT-PBEsol dataset. Structures from solid and liquid phases with model deviations larger than 0.20 eV/Å and 0.30 eV/Å, respectively, were selected. If more than 100 structures are to be picked up, a random subset of 100 was added to the dataset. In the relabeling step, DFT-SCAN+rVV10 calculations are performed on these new structures, with the same k-mesh as before and a plane-wave cutoff of 1200 eV. These newly calculated data augment the DFT-SCAN+rVV10 dataset, upon which the DP models are further refined. The iterations cease when the simplifying step does not yield any additional data. Ultimately, the DFT-SCAN+rVV10 dataset contained 1968 structures, which is 16% of DFT-PBEsol dataset. The resultant DP-SCAN+rVV10 model, derived from eight million steps of fine-tuning from the DP-PBEsol fitting network, demonstrates strong concordance with the DFT-SCAN+rVV10 calculations. As depicted in Fig. S3, this is evidenced by root-mean-square errors (RMSE) of 2.7 meV/atom for energies, 0.2 eV/Å for forces, and 15.0 meV/atom for virial, respectively, of the structures in the training dataset. The precision is comparable to those reported in previously developed potentials for similar systems [7,26,27,68].

# 3. Validation against DFT and experimental data

## 3.1. Lattice parameters and elastic constants

The lattice parameters and energies of ZrO₂ polymorphs and cubic CeO₂ predicted by DFT-SCAN+rVV10 and DP-SCAN+rVV10 are presented in Table S1. Both sets of results show excellent agreement
<figure><img src="image_6.png" /><figcaption>Fig. 3. (a) Workflow of active and transfer learning processes used in training the DP models. (b) Representative initial configurations of crystals, surfaces, and liquid.</figcaption></figure>

with experimental values. The temperature-dependent lattice parame-
ters of low-pressure ZrO₂ phases, as calculated by the DP models, are
depicted in Fig. S4. Both DP-SCAN+rVV10 and DP-PBEsol models align
closely with the experimental data, although DP-SCAN+rVV10 exhibits
marginally larger deviations, approximately 0.4%, for the tetragonal
c<sub>t</sub> and monoclinic c<sub>m</sub> parameters. Comparison of the predicted room-
temperature lattice parameters of CeO₂-doped ZrO₂ phases with ex-
perimental data is shown in Fig. 4, indicating similar magnitude of
deviations. Such a deviation falls within an acceptable range for ma-
chine learning interatomic potentials. A more significant difference
is found on the energetics of ZrO₂ polymorphs in Table S1, espe-
cially Pbca, Pbc2<sub>1</sub>, and t-ZrO<sub>2</sub>. Notably, the DP-SCAN+rVV10 model
demonstrates superior accuracy in determining the energy of tetragonal
ZrO<sub>2</sub> when compared to the DP-PBEsol model in Table S1. This en-
hanced precision is crucial for accurately capturing the thermodynamic
properties of the CeO<sub>2</sub>-ZrO<sub>2</sub> system (see Sections 3.2 and 3.3).

Table S2 details the elastic constants of ZrO₂ phases and c-CeO₂,
as determined by DFT and DP models. While the results from the
DP models generally align with those from DFT calculations, notable
discrepancies are observed in the calculated C₁₁ values for cubic and
tetragonal phases when compared with experimental data. C₁₅, C₂₅,
and C₃₅ values of m-ZrO₂ also exhibit opposite signs between exper-
iments and theoretical calculations. These inconsistencies in specific
elastic components were also noted in prior DFT studies on the elastic
properties of ZrO₂ with various exchange-correlation functionals [70–
73]. The root cause of these deviations remains unclear, potentially
stemming from inherent limitations in the accuracy of DFT calculations.

3.2. *Pressure–temperature phase diagram of ZrO<sub>2</sub>*

Thermodynamic integrations [74] implemented in the dpti (deep
potential thermodynamic integration) package [75] were performed
to calculate the free energy of various ZrO₂ phases and obtain the
transformation temperatures of a series of phase transformations at
ambient pressure, i.e., m-ZrO₂→t-ZrO₂→c-ZrO₂→L-ZrO₂. The detailed
calculation procedure is presented in the Supplementary Materials of
Ref. [7]. The results are shown in Table 1 along with other thermody-
namic properties. An excellent agreement is achieved between the DP-
SCAN+rVV10 predictions and experiments. Similarly, transformation

pressures at 300 K from m-ZrO₂ to Pbca-ZrO₂, and then to Pmna-ZrO₂,
were calculated, yielding pressures of 3.1 GPa for m-ZrO₂→Pbca-ZrO₂
and 8.8 GPa for Pbca-ZrO₂→Pmna-ZrO₂ transformations, respectively.
Notably, the calculated pressure for the m-ZrO₂→Pbca-ZrO₂ trans-
formation closely matches the experimentally observed value of 3.0
GPa [60].

Using the calculated transformation points as initial references,
the coexistence lines between various ZrO₂ phases were determined
using Gibbs-Duhem integration [79]. Fig. 5 presents the phase diagram
constructed using the DP-SCAN+rVV10 model, overlapped with exper-
imental transformation points and coexistence lines reported in the
literature [59,60,80,81]. A comparison between the DP-SCAN+rVV10
model calculations and experimental observations reveals a remark-
able quantitative agreement (Fig. 5). Similar quantitative accuracy
with experiments was only achieved by MLFF-RPA<sup>d</sup> from Ref. [28]
(see Table 1) using a training dataset based on higher-level quantum-
mechanical method than DFT, i.e., random phase approximation (RPA).
However, this potential only covers the monoclinic, tetragonal, and
cubic phases at ambient pressure. A comparison in coverage of ZrO₂
phases by different machine learning potentials is given in Table S3,
showing that only the present potential and the potential from our
previous study [7] can cover all the experimentally observed ZrO₂
phases. However, the DP-ZrO₂ model from our previous study [7]
exhibited a significantly narrower temperature range for t-ZrO₂ and a
small stability region for the Pbca-ZrO₂ [7], as shown in Table 1 and
Fig. 5. This is probably due to the overestimated energies of t-ZrO₂ and
Pbca-ZrO₂ by DFT-PBEsol.

3.3. Thermodynamics of diffusionless tetragonal/monoclinic transforma-
tions in CeO₂-ZrO₂

Firstly, a straightforward validation of solute energies predicted by DFT and the DP model is presented in Table S5. The results obtained with both methods show good agreement, with relative differences within 10%. Additionally, the solute energies in the tetragonal phase are lower than those in the monoclinic phase, consistent with the stabilizing effect of CeO₂ on the tetragonal phase.

Experimental data to validate the CeO₂-ZrO₂ part of the poten-
tial are more scarce. Here we consider the thermodynamics of diffu-
sionless tetragonal/monoclinic transformations that is central to the
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Variation of lattice parameters of tetragonal and monoclinic phases with CeO₂ composition at 298 K. (a) *a*<sub>m</sub>, *b*<sub>m</sub>, *c*<sub>m</sub>, *a*<sub>t</sub>, and *c*<sub>t</sub>; (b) *β* (angle between *a*<sub>m</sub> and *c*<sub>m</sub> axes). Dashed lines represent predictions from the DP-SCAN+rVV10 model and crosses are experimental data [69].</figcaption></figure>

<caption>**Table 1**</caption>

Thermodynamic properties of phase transformations between low-pressure polymorphs of ZrO₂ at 0 Pa predicted by the present DP-SCAN+rVV10 model, the machine-learning force field (MLFF) of Ref. [28], and the DP-ZrO₂ model of Ref. [7] compared to experimental data.

<table><thead><tr><th></th><th></th><th>T<sub>0</sub> (K)</th><th>ΔH (kJ/mol)</th><th>ΔV (%)</th></tr></thead><tbody><tr><td rowspan="4">m-ZrO<sub>2</sub>/t-ZrO<sub>2</sub></td><td>Experiment</td><td>1367 [59]</td><td>5.43 ± 0.31 [76]</td><td>—</td></tr><tr><td>DP-SCAN+rVV10</td><td>1369</td><td>5.67</td><td>—</td></tr><tr><td>MLFF-RPA<sup>d</sup> in Ref. [28]</td><td>1415</td><td>6.69</td><td>—</td></tr><tr><td>DP-ZrO<sub>2</sub> in Ref. [7]</td><td>1610</td><td>—</td><td>—</td></tr><tr><td rowspan="4">t-ZrO<sub>2</sub>/c-ZrO<sub>2</sub></td><td>Experiment</td><td>2584 [59]</td><td>3.4 ± 3.1 [77]</td><td>0.55 ± 0.09 [78]</td></tr><tr><td>DP-SCAN+rVV10</td><td>2616</td><td>6.9</td><td>0.39</td></tr><tr><td>MLFF-RPA<sup>d</sup> in Ref. [28]</td><td>2585</td><td>5.14</td><td>—</td></tr><tr><td>DP-ZrO<sub>2</sub> in Ref. [7]</td><td>2400</td><td>—</td><td>—</td></tr><tr><td rowspan="3">c-ZrO<sub>2</sub>/L-ZrO<sub>2</sub></td><td>Experiment</td><td>2983 [59]</td><td>55 ± 7 [78]</td><td>—</td></tr><tr><td>DP-SCAN+rVV10</td><td>2964</td><td>59.9</td><td>—</td></tr><tr><td>DP-ZrO<sub>2</sub> in Ref. [7]</td><td>2630</td><td>—</td><td>—</td></tr></tbody></table>

shape memory and superelastic properties of CeO₂-ZrO₂ ceramics. This validation involved conducting thermodynamic integrations [74] at different fixed CeO₂ compositions to compute the transformation temperatures of diffusionless tetragonal/monoclinic phase transformation. Note that we did not include configurational entropy since we consider here a diffusionless transformation, in accordance with a previous thermodynamic study [82]. The results are shown in Fig. 6, showcasing that the model predictions of equilibrium temperatures and latent heats align closely with experimental data [9,61] over a large range of CeO₂ concentrations.

Additionally, Fig. 7 presents the lattice parameters and λ₂ values at equilibrium temperatures under different lattice correspondences, as a function of the CeO₂ composition. The calculated lattice parameters are influenced by both composition and the prediction accuracy of transformation temperatures. Fig. 7 indicates that the DP-SCAN+rVV10 model captures a trend consistent with experimental findings [9]. We also calculated the λ₂ values for different correspondences, equal to bₘ/cₜ for correspondence B and bₘ/aₜ for correspondences A and C. Notably, the λ₂ value of correspondence B approaches one as the CeO₂ composition increases, since bₘ and cₜ approach each other. This trend agrees with the prior study on thermal transformations hysteresis in CeO₂-ZrO₂ ceramics [9]. Overall, these comparisons effectively highlight the DP-SCAN+rVV10 model capability to faithfully replicate the thermodynamic behavior of CeO₂-ZrO₂ ceramics. This consistency is not only comparable to DFT calculations but extends to the level of precision observed experimentally.

# 4. Applications

## 4.1. Transformation-induced plasticity in a quasi-2D CeO₂-ZrO₂ polycrystal

In this section, we use the DP model to simulate the complex microstructural evolutions that take place during the uniaxial deformation of a polycrystalline sample. Employing the DP-SCAN+rVV10 model by using LAMMPS [54] combined with the DeePMD-kit [47,48], we investigate the TRIP effect in tetragonal polycrystals of CeO₂-ZrO₂ ceramics. Initially, a quasi-2D tetragonal polycrystal oriented normal to the [010]₁ direction (parallel to the y-axis) is subjected to uniaxial compression along z-axis, to showcase the DP model's ability to capture key features of the TRIP effect. The quasi-2D simulation facilitates the visualization of the transformation process and allows for the use of larger grain sizes, as in a recent simulation of deformation-induced martensitic transformation in polycrystalline iron [83]. Using Atomsk's polycrystal mode [84], a periodic polycrystalline configuration with six grains was created with dimensions of 50 nm × 6aₜ × 50 nm (x × y × z), comprising approximately 680,000 atoms. To simulate the superelasticity phenomenon, a 14% CeO₂ composition was adopted to ensure that the monoclinic and tetragonal phases have almost the same free energy at 300 K (Fig. 6). Relaxation of the polycrystal was accomplished through an iteration of FIRE [85] and conjugate gradient minimization steps, each consisting of 1000 iterations, until reaching a
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Pressure-temperature phase diagram of ZrO₂ polymorphs predicted by DP-SCAN+rVV10 (red solid lines) and DP-ZrO₂ [7] (cyan solid lines) compared to experimental data. Experimental transformation points (crosses) are from Ref. [59] (m/t, t/c, and c/L) and Ref. [80] (Pbca/Pmna). Experimental boundary lines of forward (dashed lines) and inverse transformations (dotted lines) are from Ref. [81] (m/t), Ref. [60] (m/Pbca and t/Pbca at 3–8 GPa), and Ref. [80] (t/Pbca at 6–12 GPa and t/Pmna). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** (a) Transformation temperature and (b) latent heat of the diffusionless tetragonal/monoclinic transformation in the CeO₂-ZrO₂ system predicted by DP-SCAN+rVV10 (dashed line) while crosses are experimental data from Ref. [9] (blue crosses) and Ref. [61] (red crosses). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

maximum atomic force of $2.5 \times 10^{-4}$ eV/Å. Following relaxation, the system was equilibrated at 300 K for 30 ps, then uniaxially compressed along the z-axis at a strain rate of $2 \times 10^{-2}$/ns (the same strain rate is used in following simulations). A template matching method was applied to identify the different tetragonal, monoclinic, and Pbc2₁ variants. The Zr-O bonding vectors of different variants, transformed from the initial tetragonal grain via known lattice correspondences, are defined as reference templates to capture the approximate transformed crystalline orientations. Correspondence between bonds in the simulation and a reference template is established by matching the closest bonding vectors. The correct variant and its precise orientation are then determined by minimizing the deviation between the bonds observed in the simulation and those in the reference templates. For details, readers

are referred to the Supplementary Materials. Starting from the initial tetragonal phase (dark blue in Fig. 8), ferroelastic switching leads to two tetragonal variants (lighter blue). For each lattice correspondence of transformations from tetragonal to monoclinic and Pbc2₁ phases, the four variants are respectively shown with a different color lightness, as indicated in the colormap of the figures. The identified variants can be confirmed by visual inspection of their atomic structures from low-index directions, compared to the reference atomic structures provided in Supplementary Materials.

The compression-induced microstructural evolution is illustrated in Fig. 8 and Supplementary Movie 1, while the corresponding stress-strain curve is plotted in Fig. 9(a). In Fig. 8(a), all the tetragonal grains are numbered and the orientations of the *a*- and *c*-axes are shown,
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Lattice parameters of tetragonal and monoclinic phases and $\lambda_2$ values of the tetragonal/monoclinic transformation for different lattice correspondences at transformation temperatures as a function of CeO$_2$ composition. Lines and crosses are predictions by DP-SCAN+rVV10 and experimental data [17], respectively.</figcaption></figure>

revealing that the loading directions of Grains 1 to 3 are close to [001]$_t$, whereas those of Grains 4 and 5 are near [101]$_t$, and Grain 6 is close to [100]$_t$. Initially, ferroelastic switching of tetragonal phase (($a_t, b_t, c_t$)→($c_t, b_t, a_t$)) is observed in Grains 1 and 2 under compression, as shown in Fig. 8(b), and subsequently in Grain 3 (Fig. 8(c)). Following this, monoclinic phase nucleation with correspondence C occurs within and along the grain boundaries of Grains 4 and 5, respectively, as shown in Fig. 8(b). These nuclei evolve into martensite plates with OR C-2 in Grain 4 and OR C-1 in Grain 5, respectively. Subsequently, a martensite plate forms in Grain 1 along the boundary with Grain 5, while a similar martensitic transformation occurs in Grain 2 along the boundary with Grain 4 (Fig. 8(c)), exemplifying an autocatalytic transformation that propagates across multiple grains [6]. In the case of Grain 6, the monoclinic phase nucleates at the grain boundaries and grows with OR C-1, as depicted in Supplementary Movie 1. It should be noted that for the martensitic transformations after ferroelastic switching in Grains 1 to 3, the martensitic transformation corresponds to ORs C-1 and C-2 with respect to the grain after ferroelastic switching, but corresponds to ORs A-1 and A-2 with respect to the original tetragonal orientation, respectively. While martensite transformations for those grain orientations are not observed in the experimental micropillar compression [10], the transformation-induced stresses from neighboring grains may help martensite nucleation in the polycrystal in contrast with a single crystal.

At a later stage of the transformation (Fig. 8(d)), the martensite plate thickens, predominantly via a propagation perpendicular to the habit plane once the martensite tip has crossed the entire grain. Additionally, retained tetragonal phases are observed in each grain, likely due to grain boundary constraints and the limited extent of TRIP at this stage. Notably, some monoclinic phases exhibit two habit planes, as shown in the boxed region in Fig. 8(d). A closer examination reveals that the retained tetragonal phase is sheared in order to form an OR C-1 variant along one habit plane and an OR C-2 with the other. Both ORs may thus be present in a single grain. The transition from a single facet to two facets in the late stage of a phase transformation was also observed in a recent molecular dynamics simulation of NiTi shape memory materials [86].

Subsequently, the compressed polycrystal was gradually unloaded by reducing the applied strain with the same strain rate until the stress was back to zero, a process detailed in Supplementary Movie 2. During unloading, inverse transformations from monoclinic to tetragonal were

observed in each grain, ultimately restoring a fully tetragonal polycrystalline structure. The inverse transformations occurred mainly through the migration of the existing interfaces between monoclinic and the retained tetragonal phases. However, no reverse ferroelastic switching was observed, likely due to the minor shape changes being accommodated elastically. The grains thus recovered the orientation they had either before deformation or after ferroelastic switching. The observed residual strain following the initial cycle of loading and unloading may thus be attributed to this absence of inverse ferroelastic switching. Further investigation through two additional cycles of loading and unloading reveals an overlap of the stress-strain curves, indicating a pronounced superelastic behavior, a characteristic prominently observed in CeO$_2$-ZrO$_2$ ceramics at high CeO$_2$ compositions [1].

Tensile loading of the quasi-2D polycrystal at 14% CeO$_2$ was also conducted, as illustrated in Fig. S9 and Supplementary Movie 3, with the stress-strain curve shown in Fig. 9(b). The TRIP phenomena observed during compression are mirrored in the tensile mode, showing phase transformations to the monoclinic phase under ORs C-1 and C-2. Supplementary Movie 4 presents additional results of tensile loading along the x-axis (horizontal), demonstrating the occurrence of ferroelastic switching under tension in Grains 1 and 2 where the loading direction is near [100]$_t$. Compared with compression, the phase transformation under tensile loading exhibits a lower critical stress and is more complete, likely due to the positive volume expansion associated with the tetragonal-to-monoclinic transformation. Once the martensitic transformation finishes, a noticeable increase in stress precedes crack nucleation and intergranular propagation, along with an incomplete inverse transformation back to the tetragonal phase due to stress relief at the surfaces.

Overall, the martensitic transformation within CeO$_2$-ZrO$_2$ ceramics unfolds in a progressive manner, forming a uniform C variant (potentially after ferroelastic switching) with possible retained tetragonal regions near grain boundaries. This process effectively accommodates the external strain and thereby enhances the ceramic plasticity. In contrast, when subjecting pure ZrO$_2$ polycrystals to compression, the transformation is primarily driven by the large chemical driving force associated with the difference in cohesive energy between the tetragonal and monoclinic phases, rather than by the applied external stress. This leads to a rapid and complete martensitic transformation immediately beyond a lower critical stress, with a self-accommodated microstructure of martensite plates as shown in Fig. S10 and Supplementary Movie 5. In summary, the quasi-2D simulations
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Transformation-induced plasticity in a quasi-2D polycrystal at 300 K during compressive loading at (a) $\epsilon = 0\%$, (b) $\epsilon = 2.4\%$, (c) $\epsilon = 3.3\%$, and (d) $\epsilon = 4.5\%$. The orientation of all quasi-2D grains is normal to [010]$_t$, defining the viewing perspective. Zr atoms are colored according to the variant they belong to. Different variants of the same lattice correspondence are shown with different lightness. Unidentified Zr atoms are colored white. Small gray atoms denote oxygen. Orientation markers: pink and cyan arrows indicate directions parallel to [100], and [001], respectively, while red and blue arrows align with [100]$_m$ and [001]$_m$. These orientation markers are used consistently throughout this study. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Stress-strain curve of a quasi-2D polycrystal at 300 K under (a) cyclic compressive loading and unloading, and (b) tensile loading.</figcaption></figure>

discussed in this section demonstrate the DP-SCAN+rVV10 model's capability to accurately capture the key aspects of the mechanical response observed in CeO₂-ZrO₂ ceramics. This includes phenomena such as ferroelastic switching, autocatalytic/self-accommodating martensitic transformation, superelastic behavior, and fracture.

## 4.2. Tetragonal-to-monoclinic transformations in nanopillar under $\langle 110 \rangle_t$ loading

While quasi-2D polycrystal simulations offer insights into the TRIP effect in CeO₂-ZrO₂ ceramics, they fall short of capturing the full

complexity revealed by 3D modeling. A notable limitation of quasi-2D simulations is the constrained periodicity along the *y* direction, which limits the deformation out of the *xz* plane. Moreover, the *y* direction is parallel to [010]*t* for all the tetragonal grains, which restricts the variety of grain orientations and loading directions. Loading along other directions demonstrates non-trivial mechanical properties in previous experiments. An unresolved question is the absence of TRIP effect in micropillars under [110]*t*, compression despite correspondence B exhibiting the highest Schmid factor for the tetragonal-to-monoclinic transformation [10]. This contrasts with the readily observed TRIP effect under [101]*t*, loading directions, as demonstrated in both the quasi-2D simulations and previous micropillar compression tests [10].
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Minimum energy paths of tetragonal-to-monoclinic transformations of correspondence B predicted by the DP-SCAN+rVV10 potential under different stresses. (a) The 24-atom supercell of phases at 12.5% CeO₂ (1 Ce atom and 7 Zr atoms) along the transformation path, including the tetragonal phases before (T) and after (T') ferroelastic switching and the final monoclinic phase (M). The compressive and tensile stresses along (110)ₜ, that favors this transformation path are marked by arrows. Minimum energy paths under different uniaxial [110], compressive stresses and [110], tensile stresses are shown in (b) and (c), respectively.</figcaption></figure>

To explore this issue, we performed compression simulations along (110)ₜ of a cylindrical nanopillar at 300 K, as shown in Supplementary Movie 6. A periodic boundary condition was applied along the loading direction, while free boundary conditions were set in the other directions. The results show that the tetragonal-to-monoclinic phase transformation does not occur, even at 10% compressive strain, aligning with experimental observations. To investigate the origin of this behavior, we conducted finite deformation nudged elastic band calculations [87] to determine the energy barrier for correspondence B's phase transformation under varying compressive stresses along (110)ₜ, as shown in Fig. 10(b). Our results indicate that the energy barrier for ferroelastic switching along the transformation path (Fig. 10(a)) of correspondence B becomes the critical bottleneck, while the energy barrier for monoclinic transformations decreases significantly with increasing compressive stress. This persistent energy barrier of ferroelastic switching can be attributed to its tensile strain component along (110)ₜ, which is unfavorable under external compressive stress.

To further verify the importance of ferroelastic switching along correspondence B's path, we calculated the energy profile under (110)ₜ, tensile stress, as shown in Fig. 10(c). In this case, the tensile stress promotes ferroelastic switching, and both energy barriers decrease with increasing tensile stress. We also performed tensile loading simulations of the nanopillar, as shown in Fig. 11 and Supplementary Movie 7. The results reveal that martensitic plates nucleate at the surface, almost simultaneously with ferroelastic switching. Once nucleation occurs, the phase transformation proceeds via tetragonal/monoclinic interfaces under both correspondences B and C. These findings suggest that the resolved strain of intermediate states along the transformation path plays a critical role, in addition to the final state, in determining whether a TRIP mode is active. This factor should be considered in the future when assessing the operational conditions for TRIP.

The simulation also reveals that the TRIP effect for correspondence B in CeO₂-ZrO₂ ceramics is hindered under compressive loading but

may become active under tensile loading, a phenomenon that warrants experimental validation at the microscale, using methods such as push-to-pull testing. In previous tensile tests on single-crystal micro-beams of 16% CeO₂-doped ZrO₂ [88], the TRIP effect was not observed prior to fracture. However, tensile loading directions near (110)ₜ and (101)ₜ were not tested. Other tensile experiments on polycrystals [89,90] have demonstrated the TRIP effect under tensile stress, highlighting its competition with fracture, which varies with microstructural factors such as grain size. It is also important to note that the current atomistic simulation does not account for potential defects, such as micro-voids, which could initiate fracture. As a result, ceramics under experimental conditions are expected to be more brittle than predicted by our simulations.

### 4.3. Transformation-induced plasticity in a 3D CeO₂-ZrO₂ polycrystal

Although correspondence B is rarely observed in single-crystal micropillar compression experiments, studies on polycrystalline systems reported stress-induced correspondence B through indentation [13,14]. To explore the differences in the activation of correspondence B between single- and poly-crystalline systems under compression, this section focuses on 3D polycrystal simulations. Fig. 12(a) showcases a 3D polycrystal with about 940,000 atoms under full periodic boundary conditions, equilibrated at 300 K following the method outlined in Section 4.1. This configuration consists of six grains, each about 10 nm in diameter. Their loading directions are reported in the standard tetragonal stereographic projection triangle in Fig. 12(f). In contrast with the quasi-2D model, this 3D model includes a spectrum of representative loading directions for the tetragonal phase.

Upon applying compressive loading at a strain rate of 2×10⁻²/ns and 300 K, we observed microstructural evolutions depicted in Fig. 12(a)-(d) and Supplementary Movie 8. Fig. 12(e) highlights these phase
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Transformation-induced plasticity in a nanopillar at 300 K during tensile loading at (a) $\epsilon = 3.0\%$, (b) $\epsilon = 3.4\%$, (c) $\epsilon = 3.7\%$, and (d) $\epsilon = 7.2\%$. (e) Corresponding stress–strain curve. The pillar is cut in the middle normal to the viewing direction [001]. The vertical direction in the figure corresponds to ⟨110⟩. Dark and light blue of metal atoms show tetragonal variants of initial orientation and after ferroelastic switching, respectively. Metal atoms in monoclinic phase with correspondence B to initial tetragonal are colored brown and those at the surface are colored white. Gray atoms denote oxygen. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

transformations on the stress–strain curve. Unlike micropillar compression experiments, which frequently show fracture or slip along loading directions far from [101]ᵢ, especially near [001]ᵢ, [100]ᵢ and [110]ᵢ [10], the present 3D simulations reveal transformations into monoclinic phases within each grain at a high stress level. Fracture is not observed probably due to a lack of micro-cracks and the small grain size of the present polycrystal. Fig. 12(f) shows the ORs or lattice correspondences of monoclinic phases relative to their original tetragonal orientations. Differing from experimental observations of single crystals which report only correspondence C during micropillar compression [10,12], the present simulations predict correspondences A and B, in addition to correspondence C. The transformation process

is analyzed below.

The tetragonal-to-monoclinic transformation first appears in Grain 1, as seen in Fig. 12(b), with an OR C-2 (see also Fig. S11(a) in Supplementary Materials). This OR agrees with observations from quasi-2D simulations, as expected since the loading direction is near [101]ᵢ. Subsequent observations identify monoclinic plates with OR C-2 in Grains 3 and 6. The structures of the tetragonal/monoclinic interfaces inside both grains are available in the Supplementary Materials (Figs. S11(b) and (c)). An earlier monoclinic transformation in Grain 1, relative to Grains 3 and 6, is expected because of a larger Schmid factor for the C correspondence near [101]ᵢ, consistent with the experimental observation of a lower critical stress near [101]ᵢ [10].
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Transformation-induced plasticity in a 3D polycrystal at 300 K during compressive loading at (a) $\epsilon = 0.0\%$, (b) $\epsilon = 5.1\%$, (c) $\epsilon = 6.8\%$, and (d) $\epsilon = 9.7\%$. (e) Corresponding stress–strain curve, labeled with events of phase transformations or ferroelastic switching. (f) ORs of the monoclinic phases formed during compression in the six grains of the polycrystal. The grains are represented in a standard triangle by the stereographic projection of the compression axis with respect to the initial tetragonal grain orientation (double cell notation). Other observed phase transformations and ferroelastic switching (FS) are also indicated.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Formation of ORs A-2 and B-2, and Pbc2₁ proceeded by ferroelastic switching during the compression of Grain 5. The grain is sliced and viewed along [001], at (a) $\epsilon = 0\%$, (b) $\epsilon = 2.2\%$, and (c) $\epsilon = 8.8\%$. Orientation markers: pink and cyan arrows indicate directions parallel to [100], and [001], respectively, while red, green, and blue arrows align with [100], [001], and [001] in monoclinic and Pbc2₁ phases. These orientation markers are used consistently throughout this study. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Formation of OR B-2 proceeded by Pbc2₁ during the compression of Grain 2. The grain is sliced and viewed along [001]ᵢ at (a) $ϵ = 0\%$, (b) $ϵ = 6.1\%$, (c) $ϵ = 8.8\%$.</figcaption></figure>

Subsequent compression triggers monoclinic phases in Grains 2 and 5, characterized by lower Schmid factors associated with correspondence C. It is found that the monoclinic phases exhibit correspondences A and B, which form through indirect transformation paths. Specifically, in Grain 5, the emergence of monoclinic phases with ORs A-2 and B-2 is initiated by a ferroelastic switching of the entire grain, followed by a direct transformation with OR C-2, as depicted in Fig. 13 and Figs. S12(a) and (b). These transformation paths are articulated as $(a_t, b_t, c_t)→(c_t, b_t, a_t)→(c_m, b_m, a_m)$ and $(a_t, b_t, c_t)→(c_t, b_t, a_t)→(c_m, a_m, b_m)$, respectively. The transformations have thus a C correspondence with respect to the ferroelastic-switched grain, but have correspondences A and B with respect to the initial tetragonal grains. While the former path is consistent with observations from quasi-2D simulations, the latter could not take place due to the limited dimension along [010]ᵢ.

Moreover, another indirect transformation path is observed, involving a localized transformation to the Pbc2₁ phase such that $(a_t, b_t, c_t)→$

$(c_t, b_t, a_t)→(b_o, c_o, a_o)$. This polar phase is known for ferroelectric applications of zirconia [91], and has been suggested as an intermediate step in the $t$-ZrO$_2$→$m$-ZrO$_2$ transformation [92]. The $t$-ZrO$_2$→Pbc2$_1$-ZrO$_2$ transformation is notably favored when the loading direction is near $a_t$ ($c_t$ before ferroelastic switching) primarily due to $a_o$ being smaller than $a_t$. A similar transformation to Pbc2$_1$ phase, with an identical lattice correspondence, occurs in Grain 2. Upon experiencing higher strain levels, the Pbc2$_1$ phase in Grain 2 undergoes a further transformation into the monoclinic phase, resulting in the full transformation path $(a_t, b_t, c_t)→(c_o, a_o, b_o)→(c_m, a_m, b_m)$ and the formation of a monoclinic structure with correspondence B, as illustrated in Fig. 14 and Fig. S12(c). The observed transformation paths of correspondence B in Grains 2 and 5 corroborate the geometric analysis proposed by Trolliard et al. [93] and previous DFT calculations (Fig. 2) [7,8].

Fig. 15 shows the nucleation and subsequent growth of monoclinic phases with OR B-2 in Grains 3 and 4. Different from the scenario
<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Formation of correspondence B from grain boundaries during the compression of Grains 3 and 4. Grain 3 is sliced and viewed along [001], at $ϵ = 8.8\%$, where a monoclinic region with OR C-2 is also observed. Grain 4 is sliced and viewed along $\langle 010 \rangle$, at $ϵ = 8.8\%$.</figcaption></figure>

in Grain 5, the transformations here do not start with a ferroelastic switching of the entire grains because this process is either only slightly favorable (in Grain 3) or unfavorable (in Grain 4). Instead, in Grain 3, the monoclinic phase grows from the grain boundaries with a ferroelastic switching front slightly ahead of the monoclinic phase. In Grain 4, which has a significant Schmid factor for correspondence B, the monoclinic phase appears only at a very late stage of the compression process. In contrast to the simulation of single-crystal nanopillar compression along $\langle 110 \rangle_t$, the observations in Grains 3 and 4 highlight the critical role of grain boundaries in lattice correspondences devoid of a direct transformation path. Smaller grain sizes will stabilize the tetragonal phase further [82], building a larger stress inside materials. Nucleation sites from grain boundaries and the internal stress may help the phase transformations to cross the energy barriers. Our simulation results are consistent with the occurrence of fracture during the compression of a micropillar (single crystal) along near $[110]_t$ [10], whereas indentation of a polycrystal induces monoclinic phases with predominantly correspondence B [13,14].

## 5. Discussion and conclusion

In this work, we developed a deep neural network potential using a combination of active and transfer learning techniques. The potential reaches DFT-level accuracy and reproduces experimental thermodynamic properties of CeO₂-ZrO₂ ceramics. The main effect of CeO₂ addition, which is faithfully reproduced by the potential, is to stabilize the tetragonal phase at low temperature. As seen in Fig. 6, the potential predicts a tetragonal-monoclinic transition temperature of about 300 K for zirconia containing 14% CeO₂, in close agreement with experimental data. For this composition, the tetragonal and monoclinic phases have thus nearly the same free energy at 300 K. This property is crucial to simulate transformation-induced plasticity in realistic conditions since, upon deformation, the tetragonal-monoclinic transformation should be triggered only by the bias induced by the applied stress. The situation is different in pure ZrO₂, where the driving force is instead the large free energy difference between both phases. With the addition of CeO₂, the potential also reproduces superelasticity [1] due to reversible tetragonal-monoclinic transformations, a

phenomenon only possible when monoclinic and tetragonal phases have close free energies biased by the applied stress. Moreover, the computational efficiency of the deep potential achieves an impressive computational speed of approximately $1.1 \times 10^{-6}$ s/(atom·step) using four GPU NVIDIA V100 (32 GB) with nearly 1 million atoms, which enables 3D simulations of polycrystalline CeO₂-ZrO₂ ceramics.

Using the deep potential, we studied the activation of tetragonal-to-monoclinic transformations under $\langle 110 \rangle_t$ uniaxial loading. This loading direction is of interest because it exhibits a large Schmid factor for correspondence B while no phase transformation is observed experimentally under compression. The simulation results indicate a strong asymmetry between $\langle 110 \rangle_t$ compression and tension. The compressive strain fails to trigger tetragonal-to-monoclinic transformation because of the increasing energy barrier of ferroelastic switching along the transformation path. Meanwhile, the tensile strain lowers both energy barriers along the transformation path and successfully induces the TRIP effect by correspondence B accompanying ferroelastic switching. The results indicate the important role of resolved strain of the intermediate structures to the TRIP effect, in addition to that of the final structure along the transformation path.

The simulations of 3D polycrystalline CeO₂-ZrO₂ ceramics revealed complex microstructural evolutions during uniaxial compression. Contrary to prior experimental findings in single-crystal micropillars [10] but consistent with the scenario in polycrystals [13,14,16,17], the simulations predict the formation of all three lattice correspondences. The grains well oriented for correspondence C, i.e., near $[101]_t$ with a high Schmid factor, transform easily with this correspondence, consistent with previous micropillar compression [10]. Conversely, transformations with correspondences A and B occur with more difficulty and involve a two-step process with either an intermediate tetragonal variant after ferroelastic switching or the Pbc2₁ phase. The reason, mostly overlooked so far, is the absence of a direct path for the correspondences A and B. As shown in Fig. 2, only correspondence C has a direct path between the tetragonal and the monoclinic phases. Thus, ferroelastic switching, which was studied theoretically in the past independently of the martensitic transformation [94], appears to play a larger role than usually considered in mediating correspondences A or B.

Generally, the transformations observed here can be understood using simple arguments based on the Schmid factor or the shape change
between phases. For instance, ferroelastic switching occurs when the grain is compressed along [001] or stretched along (010)ₜ since $c_t > a_t = b_t$. Transformations are expected with correspondences C near [101]ₜ and B near [110]ₜ, where the Schmid factor for both correspondences is maximum. In the later stages of deformation, the effects of coherency strains from neighboring grains facilitate transformations in ill-oriented grains, which aligns with a recent machine learning-based analysis of twinning in Mg [95].

Grain boundaries emerge as a crucial factor in our simulations, allowing transformations to nucleate at lower stresses compared to single crystals. Here, grain boundaries can accommodate substantial deformations without fracture, likely due to the limited grain size and absence of microcracks. Besides acting as a source of transformation, grain boundaries also act as obstacles to propagation but they can be crossed by transformed regions without breaking (see Fig. 8(c) for instance). Overall, grain boundaries facilitate a more uniform transformation-induced plasticity across various grains, echoing findings from a recent experiment on ultrafine-grained austenitic steels [96] and simulations [97]. In experimental conditions, grain boundaries in ceramics may contain impurities, microcracks, and surfaces. They are thus more brittle and may not accommodate the plasticity produced here, explaining why the ductility of zirconia polycrystals is typically limited to about 2% [98], while we were able to deform the 3D polycrystals up to 10% without fracture.

In conclusion, we believe that the present work, based on an accurate and efficient deep neural network potential, opens the door for more systematic studies of atomic-scale transformations in zirconia ceramics, a process which needs a better understanding to further develop ceramics for structural applications.

### CRediT authorship contribution statement

**Jin-Yu Zhang:** Writing – original draft, Visualization, Investigation, Conceptualization. **Gaël Huynh:** Writing – review & editing, Visualization, Investigation. **Tristan Albaret:** Writing – review & editing, Investigation. **Fu-Zhi Dai:** Methodology. **David Rodney:** Writing – review & editing, Supervision, Investigation, Funding acquisition, Conceptualization.

### 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.

### Acknowledgments

We would like to thank Dr. Yinan Wang at Artificial Intelligence for Science Institute (China) for many useful discussions on transfer learning of DP models. D. Rodney, T. Albaret & G. Huynh acknowledge funding from the ANR NANOTRIP, France (ANR-21-CE08-0019) and computational resources from the Jean-Zay machine at the Grand Équipement National de Calcul Intensif (GENCI), France (Grant No. A0130810637). J.-Y. Zhang acknowledges the funding by projects of JSPS KAKENHI (Grant No. JP23K20037). 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.

### Appendix A. Supplementary data

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

### References

[1] A. Lai, Z. Du, C.L. Gan, C.A. Schuh, Shape memory and superelastic ceramics at small scales, Science 341 (2013) 1505-1508, http://dx.doi.org/10.1126/science.1239745.

[2] 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 (2009) 1901-1920, http://dx.doi.org/10.1111/j.1551-2916.2009.03278.x.

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

[4] 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 (2020) 1482-1513, http://dx.doi.org/10.1111/jace.16903.

[5] A. Liens, H. Reveron, T. Douillard, N. Blanchard, V. Lughi, V. Sergio, R. Laquai, B.R. Müller, G. Bruno, S. Schomer, T. Fürderer, E. Adolfsson, N. Courtois, M. Swain, J. Chevalier, Phase transformation induces plasticity with negligible damage in ceria-stabilized zirconia-based ceramics, Acta Mater. 183 (2020) 261-273, http://dx.doi.org/10.1016/j.actamat.2019.10.046.

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

[7] J.-Y. Zhang, G. Huynh, F.-Z. Dai, T. Albaret, S.-H. Zhang, S. Ogata, D. Rodney, A deep-neural network potential to study transformation-induced plasticity in zirconia, J. Eur. Ceram. Soc. 44 (2024) 4243-4254, http://dx.doi.org/10.1016/j.jeurceramsoc.2024.01.007.

[8] S.-H. Guan, X.-J. Zhang, Z.-P. Liu, Energy landscape of zirconia phase transitions, J. Am. Chem. Soc. 137 (2015) 8010-8013, http://dx.doi.org/10.1021/jacs.5b04528.

[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] X.M. Zeng, A. Lai, C.L. Gan, C.A. Schuh, Crystal orientation dependence of the stress-induced martensitic transformation in zirconia-based shape memory ceramics, Acta Mater. 116 (2016) 124-135, http://dx.doi.org/10.1016/j.actamat.2016.06.030.

[11] M.D. Magalhaes, T. Douillard, H. Reveron, S. Comby-Dassonneville, T.W. Corbelius, S. Meille, M. Texier, J.S. Micha, D. Rodney, O. Thomas, J. Chevalier, *In-situ* laue micro-diffraction during compression tests on Ce-TZP single crystal micropillars, J. Eur. Ceram. Soc. 45 (1) (2025) 116794, http://dx.doi.org/10.1016/j.jeurceramsoc.2024.116794.

[12] X.M. Zeng, Z. Du, N. Tamura, Q. Liu, C.A. Schuh, C.L. Gan, *In-situ* studies on martensitic transformation and high-temperature shape memory in small volume zirconia, Acta Mater. 134 (2017) 257-266, http://dx.doi.org/10.1016/j.actamat.2017.06.006.

[13] G.R. Hugo, B.C. Muddle, R.H.J. Hannink, Crystallography of the tetragonal to monoclinic transformation in ceria-zirconia, Mater. Sci. Forum 34-36 (1988) 165-169, http://dx.doi.org/10.4028/www.scientific.net/MSF.34-36.165.

[14] G.R. Hugo, B.C. Muddle, The tetragonal to monoclinic transformation in zirconia, Mater. Sci. Forum 56-58 (1990) 357-362, http://dx.doi.org/10.4028/www.scientific.net/MSF.56-58.357.

[15] S. Deville, G. Guénin, J. Chevalier, Martensitic transformation in zirconia: Part I. Nanometer scale prediction and measurement of transformation induced relief, Acta Mater. 52 (2004) 5697-5707, http://dx.doi.org/10.1016/j.actamat.2004.08.034.

[16] Y. Wang, F. Xu, R. Gauvin, M. Kong, M. Khan, Z. Liu, Y. Zeng, Growth modes for monoclinic yttria-stabilized zirconia during the martensitic transformation, J. Am. Ceram. Soc. 100 (2017) 4874-4883, http://dx.doi.org/10.1111/jace.15004.

[17] 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.

[18] Y. Wang, M. Winhold, M. Kong, M. Khan, P. Frank, C.H. Schwalb, Y. Zeng, Investigation on the habit plane of martensitic transformation in zirconia coatings, J. Aust. Ceram. Soc. 56 (2020) 257-264, http://dx.doi.org/10.1007/s41779-019-00416-5.

[19] 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.

[20] 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.

[21] 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.
[22] 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.

[23] N. Zhang, M. Asle Zaeem, Effects of grain orientations and pre-existing defects on mechanical properties and deformation mechanisms of polycrystalline yttria-stabilized tetragonal zirconia, Materialia 9 (2020) 100553, http://dx.doi.org/10.1016/j.mtla.2019.100553.

[24] S. Fujii, K. Shimazaki, A. Kuwabara, Empirical interatomic potentials for ZrO₂ and YSZ polymorphs: Application to a tetragonal ZrO₂ grain boundary, Acta Mater. 262 (2024) 119460, http://dx.doi.org/10.1016/j.actamat.2023.119460.

[25] C. Yin, F. Ye, C.Y. Yin, D.R. Ou, T. Mori, An assessment of interatomic potentials for yttria-stabilized zirconia, Appl. Mech. Mater. 492 (2014) 239-247, http://dx.doi.org/10.4028/www.scientific.net/AMM.492.239.

[26] 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 (2021) 1-9, http://dx.doi.org/10.1038/s41524-021-00630-5.

[27] 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. Appl. 18 (2022) 054066, http://dx.doi.org/10.1103/PhysRevApplied.18.054066.

[28] P. Liu, C. Verdi, F. Karsai, G. Kresse, Phase transitions of zirconia: Machine-learned force fields beyond density functional theory, Phys. Rev. B 105 (2022) L060102, http://dx.doi.org/10.1103/PhysRevB.105.L060102.

[29] Y. Deng, C. Wang, X. Xu, H. Li, Machine learning potential for Ab Initio phase transitions of zirconia, Theor. Appl. Mech. Lett. 13 (2023) 100481, http://dx.doi.org/10.1016/j.taml.2023.100481.

[30] D. Lu, W. Jiang, Y. Chen, L. Zhang, W. Jia, H. Wang, M. Chen, DP compression: 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.

[31] 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.

[32] 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.

[33] J. Li, K. Luo, Q. An, Mobile dislocation mediated Hall-Petch and inverse Hall-Petch behaviors in nanocrystalline Al-doped boron carbide, J. Eur. Ceram. Soc. 44 (2) (2024) 659-667, http://dx.doi.org/10.1016/j.jeurceramsoc.2023.09.079.

[34] J.-Y. Zhang, Y.-S. Zhang, F. Mompiou, W.-Z. Zhang, Structures and migrations of interfaces between β precipitates and α' matrix in a Ti-2.6 wt% Mo alloy, Acta Mater. 281 (2024) 120429, http://dx.doi.org/10.1016/j.actamat.2024.120429.

[35] Y. Peng, X.-L. Si, C. Shang, Z.-P. Liu, Abundance of low-energy oxygen vacancy pairs dictates the catalytic performance of cerium-stabilized zirconia, J. Am. Chem. Soc. 146 (2024) 10822-10832, http://dx.doi.org/10.1021/jacs.4c01285.

[36] O. Rettenmaier, J.J. Gabriel, S. Patala, Modeling martensitic transformation temperatures in zirconia-ceria solid solutions using machine learning interatomic potentials, Modelling Simul. Mater. Sci. Eng. 32 (2024) 085023, http://dx.doi.org/10.1088/1361-651X/ad801e.

[37] E.H. Kisi, C.J. Howard, Elastic constants of tetragonal zirconia measured by a new powder diffraction technique, J. Am. Ceram. Soc. 81 (1998) 1682-1684, http://dx.doi.org/10.1111/j.1151-2916.1998.tb02533.x.

[38] 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 (1989) 218-227, http://dx.doi.org/10.1111/j.1151-2916.1989.tb06104.x.

[39] 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 nanoparti-cles, J. Appl. Crystallogr. 47 (2) (2014) 684-691, http://dx.doi.org/10.1107/S1600576714003331.

[40] 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 (2018) 143001, http://dx.doi.org/10.1103/PhysRevLett.120.143001.

[41] L. Tang, Z.J. Yang, T.Q. Wen, K.M. Ho, M.J. Kramer, C.Z. Wang, Short- and medium-range orders in Al90Tb10 glass and their relation to the structures of competing crystalline phases, Acta Mater. 204 (2021) 116513, http://dx.doi.org/10.1016/j.actamat.2020.116513.

[42] H. Tang, Y. Zhang, Q.-J. Li, H. Xu, Y. Wang, Y. Wang, J. Li, High accuracy neural network interatomic potential for NiTi shape memory alloy, Acta Mater. 238 (2022) 118217, http://dx.doi.org/10.1016/j.actamat.2022.118217.

[43] J.-Y. Zhang, Z.-P. Sun, D. Qiu, F.-Z. Dai, Y.-S. Zhang, D. Xu, W.-Z. Zhang, Dislocation-mediated migration of the α/β interfaces in titanium, Acta Mater. 261 (2023) 119364, http://dx.doi.org/10.1016/j.actamat.2023.119364.

[44] R. He, B. Zhang, H. Wang, L. Li, P. Tang, G. Bauer, Z. Zhong, Ultrafast switching dynamics of the ferroelectric order in stacking-engineered ferroelectrics, Acta Mater. 262 (2024) 119416, http://dx.doi.org/10.1016/j.actamat.2023.119416.

[45] Y. Wang, B. Wen, X. Jiao, Y. Li, L. Chen, Y. Wang, F.-Z. Dai, The highest melting point material: Searchied by Bayesian global optimization with deep potential molecular dynamics, J. Adv. Ceram. 12 (2023) 803-814, http://dx.doi.org/10.26599/JAC.2023.9220721.

[46] AIS Square, https://www.aissquare.com/.

[47] H. Wang, L. Zhang, J. Han, W. E, DeePMDKit: 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.

[48] 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 (2023) 054801, http://dx.doi.org/10.1063/5.0155600.

[49] L. Zhang, J. Han, H. Wang, W. Saidi, R. Car, W. E, End-to-end symmetry pre-serving 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.

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

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

[52] 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 (2019) 023804, http://dx.doi.org/10.1103/PhysRevMaterials.3.023804.

[53] D. Sheppard, P. Xiao, W. Cemelewski, 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.

[54] 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.

[55] V. Gritsenko, D. Gritsenko, S. Shaimeev, V. Aliev, K. Nasyrov, S. Erenburg, V. Tapilin, H. Wong, M.C. Poon, J.H. Lee, J.W. Lee, C.W. Kim, Atomic and electronic structures of amorphous ZrO₂ and HfO₂ films, Microelectron. Eng. 81 (2) (2005) 524-529, http://dx.doi.org/10.1016/j.mee.2005.03.056.

[56] 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 (1996) 15-50, http://dx.doi.org/10.1016/0927-0256(96)00008-0.

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

[58] 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 (2008) 136406, http://dx.doi.org/10.1103/PhysRevLett.100.136406.

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

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

[61] M. Yoshima, T. Tsuchuhashi, H. Takashina, M. Kakihana, T. Ikegami, M. Yoshimura, Tetragonal—monoclinic phase transition enthalpy and temperature of ZrO₂-CeO₂ solid solutions, J. Am. Ceram. Soc. 78 (1995) 2225-2228, http://dx.doi.org/10.1111/j.1151-2916.1995.tb08642.x.

[62] E.V. Stefanovich, A.L. Shluger, C.R.A. Catlow, Theoretical study of the stabilization of cubic-phase ZrO₂ by impurities, Phys. Rev. B 49 (1994) 11560-11571, http://dx.doi.org/10.1103/PhysRevB.49.11560.

[63] R. Ackermann, E. Rauh, C. Alexander, Thermodynamic properties of ZrO₂(g), High Temp. Sci. 7 (4) (1975) 304-316.

[64] 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 (2016) 041005, http://dx.doi.org/10.1103/PhysRevX.6.041005.

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

[66] R. Sabatini, T. Gorni, S. de Gironcoli, Nonlocal van der Waals density functional made simple and efficient, Phys. Rev. B 87 (2013) 041108, http://dx.doi.org/10.1103/PhysRevB.87.041108.

[67] 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 (2020) 043361, http://dx.doi.org/10.1103/PhysRevResearch.2.043361.

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

[69] M. Yoshima, T. Hirose, S. Katano, Y. Suzuki, M. Kakihana, M. Yoshimura, Structural changes of ZrO₂-CeO₂ solid solutions around the monoclinic-tetragonal phase boundary, Phys. Rev. B 51 (1995) 8018-8025, http://dx.doi.org/10.1103/PhysRevB.51.8018.
[70] S. Gennard, F. Corà, C.R.A. Catlow, Comparison of the bulk and surface properties of ceria and zirconia by ab initio investigations, J. Phys. Chem. B 103 (1999) 10158–10170, http://dx.doi.org/10.1021/jp9913923.

[71] G. Fadda, L. Colombo, G. Zanzotto, First-principles study of the structural and elastic properties of zirconia, Phys. Rev. B 79 (2009) 214102, http://dx.doi.org/10.1103/PhysRevB.79.214102.

[72] Z. Liang, W. Wang, M. Zhang, F. Wu, J.-F. Chen, C. Xue, H. Zhao, Structural, mechanical and thermodynamic properties of ZrO₂ polymorphs by first-principles calculation, Phys. B: Condens. Matter. 511 (2017) 10–19, http://dx.doi.org/10.1016/j.physb.2017.01.025.

[73] M. Delarmelina, M.G. Quesne, C.R.A. Catlow, Modelling the bulk properties of ambient pressure polymorphs of zirconia, Phys. Chem. Chem. Phys. 22 (2020) 6660–6676, http://dx.doi.org/10.1039/D0CP00032A.

[74] 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 (2008) 153101, http://dx.doi.org/10.1088/0953-8984/20/15/153101.

[75] deepmodeling/dpti, https://github.com/deepmodeling/dpti.

[76] Y. Moriya, A. Navrotsky, High-temperature calorimetry of zirconia: Heat capacity and thermodynamics of the monoclinic–tetragonal phase transition, J. Chem. Thermodyn. 38 (2006) 211–223, http://dx.doi.org/10.1016/j.jct.2005.05.002.

[77] 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 (2005) 2942–2944, http://dx.doi.org/10.1111/j.1551-2916.2005.00506.x.

[78] 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 (2018) 14962, http://dx.doi.org/10.1038/s41598-018-32848-7.

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

[80] 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 (2001) 174108, http://dx.doi.org/10.1103/PhysRevB.63.174108.

[81] E.D. Whitney, Electrical resistivity and diffusionless phase transformations of zirconia at high temperatures and ultrahigh pressures, J. Electrochem. Soc. 112 (1965) 91, http://dx.doi.org/10.1149/1.2423476.

[82] E.L. Pang, G.B. Olson, C.A. Schuh, Role of grain constraint on the martensitic transformation in ceria-doped zirconia, J. Am. Ceram. Soc. 104 (2) (2021) 1156–1168, http://dx.doi.org/10.1111/jace.17526.

[83] T. Shimokawa, K. Fujii, T. Niiyama, Atomic simulation study of the factors affecting nucleation in deformation-induced martensitic transformation in grains and at grain boundaries in pure iron, Acta Mater. 265 (2024) 119629, http://dx.doi.org/10.1016/j.actamat.2023.119629.

[84] P. Hirel, Atomsk: A tool for manipulating and converting atomic data files, Comput. Phys. Comm. 197 (2015) 212–219, http://dx.doi.org/10.1016/j.cpc.2015.07.012.

[85] 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.

[86] G. Plummer, M.I. Mendelev, J.W. Lawson, Molecular dynamics simulations of austenite-martensite interface migration in NiTi alloy, Phys. Rev. Mater. 6 (12) (2022) 123601, http://dx.doi.org/10.1103/PhysRevMaterials.6.123601.

[87] A. Ghasemi, P. Xiao, W. Gao, Nudged elastic band method for solid-solid transition under finite deformation, J. Chem. Phys. 151 (5) (2019) 054110, http://dx.doi.org/10.1063/1.5113716.

[88] X.M. Zeng, P. Ye, H.T. Tan, Z. Du, C.L. Gan, Tensile behavior of tetragonal zirconia micro/nano-fibers and beams in situ tested by push-to-pull devices, J. Am. Ceram. Soc. 105 (9) (2022) 5911–5920, http://dx.doi.org/10.1111/jace.18555.

[89] R.H.J. Hannink, M.V. Swain, Metastability of the martensitic transformation in a 12 mol% ceria-zirconia alloy: I, deformation and fracture observations, J. Am. Ceram. Soc. 72 (1) (1989) 90–98, http://dx.doi.org/10.1111/j.1151-2916.1989.tb05959.x.

[90] K.M. Prettyman, Ferroelastic domain formation and switching as a toughening mechanism in ceria-doped zirconia (Ph.D. thesis), The University of Utah, Salty Lake City, 1991.

[91] X. Li, Z. Liu, A. Gao, Q. Zhang, H. Zhong, F. Meng, T. Lin, S. Wang, D. Su, K. Jin, C. Ge, L. Gu, Ferroelastically protected reversible orthorhombic to monoclinic-like phase transition in ZrO₂ nanocrystals, Nature Mater. (2024) 1–8., http://dx.doi.org/10.1038/s41563-024-01853-9.

[92] L.K. Lenz, A.H. Heuer, Stress-induced transformation during subcritical crack growth in partially stabilized zirconia, J. Am. Ceram. Soc. 65 (11) (1982) C-192–C-194, http://dx.doi.org/10.1111/j.1151-2916.1982.tb10796.x.

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

[94] 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 (2014) 144109, http://dx.doi.org/10.1103/PhysRevB.90.144109.

[95] B. Yang, V. Vassilev-Galindo, J. Llorca, Application of machine learning to assess the influence of microstructure on twin nucleation in Mg alloys, NPJ Comput. Mater. 10 (2024) 1–14, http://dx.doi.org/10.1038/s41524-024-01212-x.

[96] W. Mao, W. Gong, S. Harjo, S. Morooka, S. Gao, T. Kawasaki, N. Tsuji, In situ neutron diffraction revealing the achievement of excellent combination of strength and ductility in metastable austenitic steel by grain refinement, J. Mater. Sci. Technol. 176 (2024) 69–82, http://dx.doi.org/10.1016/j.jmst.2023.07.036.

[97] W.-S. Ko, J. Hur, J.-Y. Hwang, Y.-B. Chun, Atomistic investigation of phase transformations in NiTiCu shape memory alloys, Int. J. Mech. Sci. 274 (2024) 109256, http://dx.doi.org/10.1016/j.ijmecsci.2024.109256.

[98] P.E. Reyes-Morel, J.-S. Cherng, I.-W. Chen, Transformation plasticity of CeO₂-stabilized tetragonal zirconia polycrystals: II, pseudoelasticity and shape memory effect, J. Am. Ceram. Soc. 71 (1988) 648–657, http://dx.doi.org/10.1111/j.1151-2916.1988.tb06383.x.