# Fileset

[Lu_PhysrevMaterials_2026.pdf](https://mdr.nims.go.jp/filesets/391f107f-1c74-4289-8c65-f639495f272f/download)

## Creator

[Anh Khoa Augustin Lu](https://orcid.org/0000-0003-4702-0933), Naoki Maekawa, Akane Ikeda, Koji Shimizu, [Hiroshi Masuda](https://orcid.org/0000-0003-1032-8790), [Hidehiro Yoshida](https://orcid.org/0000-0002-0759-8258), [Satoshi Watanabe](https://orcid.org/0000-0002-8069-6938)

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Study of the ion mobility in defect-laden                    <math>                      <msub>                        <mi>ZrO</mi>                        <mn>2</mn>                      </msub>                    </math>                    under an electric field using neural network with predictions for Born effective charges](https://mdr.nims.go.jp/datasets/563b1fb9-a76b-4e14-865e-962152439fdd)

## Fulltext

Study of the ion mobility in defect-laden ${\rm ZrO}_2$ under an electric field using neural network with predictions for Born effective chargesPHYSICAL REVIEW MATERIALS 10, 066001 (2026)Study of the ion mobility in defect-laden ZrO2 under an electric field using neuralnetwork with predictions for Born effective chargesAnh Khoa Augustin Lu ,1,2,* Naoki Maekawa ,1 Akane Ikeda ,1 Koji Shimizu ,3Hiroshi Masuda ,1 Hidehiro Yoshida ,1 and Satoshi Watanabe 1,†1Department of Materials Engineering, The University of Tokyo, Tokyo 113-8656, Japan2Research Center for Materials Nanoarchitectonics (MANA), National Institute for Materials Science (NIMS), Tsukuba 305-0044, Japan3Materials DX Research Center, National Institute of Advanced Industrial Science and Technology, Tsukuba 305-8568, Japan(Received 18 August 2025; revised 12 March 2026; accepted 4 May 2026; published 2 June 2026)Unusual mass transport behavior in tetragonal ZrO2 ceramics has attracted attention under flash events inducedby strong electric fields. However, this observation cannot be attributed solely to Joule heating, suggesting theimportance of understanding the ion behaviors associated with defective states under a strong electric field.Previous studies have studied the impact of an external electric field but were typically limited to fixed formalcharges for the ions. In this work, to incorporate the response of ions to an electric field, we calculate Borneffective charges, and use them in addition to the energy and forces to train neural network potentials. Ourmolecular dynamics simulations using trained models show that under an applied electric field, the diffusivityof oxygen ions is enhanced in defect-laden ZrO2 with a preexisting oxygen vacancy, which could be associatedwith the observed unusual mass transport behavior. This is a milestone towards the accurate description ofdefect-laden materials under an applied electric field.DOI: 10.1103/jcsd-dbl2I. INTRODUCTIONHigh-strength structural ceramics are known for their hard-ness, wear resistance, and chemical stability. However, theyare typically brittle and have limited plastic deformability.Among them, tetragonal zirconia (ZrO2) is one of the mostwidely used examples, with applications in mechanical tools,medical implants, oxygen sensors, and thermal barrier coat-ings [1–4].Metastable phases of zirconia can be maintained at roomtemperature by adding stabilizers such as calcium, magne-sium, cerium, or yttrium [5]. When stabilized with a moderateamount of yttria (Y2O3), yttria-stabilized tetragonal zirconiapolycrystal (Y-TZP) is known for its combination of mechan-ical strength and fracture toughness [6]. Recently, unusuallyrapid sintering and plastic deformation behaviors of Y-TZPunder strong electric fields have attracted attention from re-searchers [7–11].Cologna et al. reported that when Y-TZP is placed underan electric field, the sintering process can be completed in afew seconds at a low temperature of 850 ◦C instead of severalhours at 1450 ◦C required for conventional sintering [12]. Therapid sintering is associated with the flash event, where the*Contact author: LU.Augustin@nims.go.jp†Contact author: watanabe@cello.t.u-tokyo.ac.jpPublished by the American Physical Society under the terms of theCreative Commons Attribution 4.0 International license. Furtherdistribution of this work must maintain attribution to the author(s)and the published article’s title, journal citation, and DOI.electrical conductivity nonlinearly increases, and has beenobserved in oxide ceramics [8], suggesting the enhancementof oxygen mobility in tetragonal ZrO2. It is known that theflash phenomenon is triggered by thermal runaway becauseof the combination of Joule heating and an increase in elec-trical conductivity. However, Raj calculated the temperatureincrease because of Joule heating in a Y-TZP sample duringthe flash phenomenon and showed that it was several hun-dred degrees lower than the temperature required to sinterthe sample in a few seconds [7]. Therefore, the rapid masstransport cannot be explained by the effect of Joule heatingalone. Raj also implied that the action of the electric field andhigh temperature can significantly increase the mass transportrate by generating a large number of Frenkel pair defects. Mo-tomura and coworkers observed that instead of brittle fracture,large plastic deformation was obtained in Y-TZP under a flashevent [13]. They experimentally demonstrated the nonthermaleffect of the flash event in accelerating plastic deformation.Itoh et al. detected an excessive concentration of chargedoxygen vacancies in 3Y-TZP after flash sintering using photo-luminescence spectroscopy [14]. These studies highlight theimportance of understanding the ion behaviors in defectivestates, which can be related to the mass transport behaviorspromoted nonthermally under flash events.Several experimental studies suggested the importance ofthe response of the system to an external electric field in flashsintering phenomena [12,13]. Xu and coworkers theoreticallystudied oxygen diffusion in zirconia by performing moleculardynamics simulations with a Coulomb-Buckingham potential[15]. They incorporated the effect of the electric field byadding an electrostatic force F = qe, where q denotes theformal charge of each ion [16]. These simulations, based on2475-9953/2026/10(6)/066001(11) 066001-1 Published by the American Physical Societyhttps://orcid.org/0000-0003-4702-0933https://orcid.org/0009-0002-4893-1974https://orcid.org/0009-0009-5640-8201https://orcid.org/0000-0001-5622-9582https://orcid.org/0000-0003-1032-8790https://orcid.org/0000-0002-0759-8258https://orcid.org/0000-0002-8069-6938https://ror.org/057zh3y96https://ror.org/026v1ze26https://ror.org/01703db54https://crossmark.crossref.org/dialog/?doi=10.1103/jcsd-dbl2&domain=pdf&date_stamp=2026-06-02https://doi.org/10.1103/jcsd-dbl2https://creativecommons.org/licenses/by/4.0/ANH KHOA AUGUSTIN LU et al. PHYSICAL REVIEW MATERIALS 10, 066001 (2026)TABLE I. Details of the training data.Structure type Materials Project ID Symmetry group Defect Charge DatapointsCubic 1565 Fm3̄m Pristine 0 1978Cubic 1565 Fm3̄m O vacancy +2 1942Tetragonal 2574 P42/nmc Pristine 0 1999Tetragonal 2574 P42/nmc O vacancy +2 2000Tetragonal 754 403 I41/amd O vacancy +2 597Monoclinic 2858 P121/c1 O vacancy +2 1587Total 10 103empirical potentials, suggested that diffusion is mediated bydefect formation and migration. However, the accuracy ofsuch empirical potentials remains uncertain, and the use of afixed formal-charge model neglects the polarization response.In recent years, the emergence of high-accuracy neural net-work potentials (NNPs), pioneered by Behler and Parrinellohigh-dimensional neural network potentials [17], has deeplymodified the landscape of molecular dynamics simulations,enabling near-DFT accuracy calculations for a fraction of thecost. Following this breakthrough, several models have beenemerging, such as DeepMD [18], SchNet [19], or NequIP[20], which is the basis for Allegro [21] and SevenNet [22].While NNPs have enabled an accurate description ofinteratomic interactions, their prediction accuracy may beinsufficient for the cases where the long-range electrostaticinteractions play a crucial role. To address this, Zhang et al.proposed an extension of DeepMD, DPLR, that adds a deepneural network to predict distributions of spherical Gaus-sian charges located at ionic and valence-electron Wanniercentroids [23]. On the other hand, Shimizu and coworkersdeveloped a neural network architecture, which directly pre-dicts the Born effective charge tensor relevant to a fixed fielddirection [24] at each time step. Their method does notconsider the long-range electrostatic interaction explicitly, incontrast to DPLR. We can expect, however, that their methodcan describe the ion movements under an applied electricfield with sufficient accuracy since the Born effective chargeis directly related to the force change owing to an appliedelectric field.In this work, we performed first-principles moleculardynamics simulations (FPMD) and density functional per-turbation theory (DFPT) of zirconia. As zirconia in its pureform exists with monoclinic symmetry (P21/c) at room tem-perature, with a transition to tetragonal symmetry (P42/nmc)at 1170 ◦C, then cubic symmetry (Fm3̄m) at 2370 ◦C and amelting point of 2716 ◦C [1], we included these phases tobuild a dataset to train neural network potentials with a goodaccuracy on energy, forces and Born effective charges. Ourmolecular dynamics simulations confirm previous results thatthe presence of oxygen vacancies enhances the diffusivity, butmore importantly, reveal that the applied electric field furtherincreases it in the direction of the electric field. While thepresent study focuses on pure ZrO2, the insight into oxygenion mobility under an electric field and in the presence ofoxygen vacancies provides qualitative guidance for the morecomplex yttria-stabilized zirconia.Our results are an important milestone towards the un-derstanding of the properties of ceramics under an electricfield. More generally, their application to ionic systems shouldprovide valuable insights into their response to an appliedelectric field.II. METHODOLOGYA. First-principles calculationsFirst-principles molecular dynamics (FPMD) calculationswere performed to generate atomic configurations to build thetraining set for training machine learning potentials (MLPs)using the Vienna Ab initio Simulation Package (VASP) [25,26].We used the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [27] and the projector-augmented-wave(PAW) method [28] for describing the interaction betweeninner-shell electrons and nuclei. The energy cutoff for planewaves was set to 520 eV, spin polarization was considered,and we used Monkhorst-Pack k-point sampling [29]. We useda 2 × 2 × 2 grid for cubic ZrO2 (Fm3̄m), 3 × 3 × 2 and3 × 3 × 1 for the two tetragonal phases of ZrO2 (P42/nmc andI41/amd), and 2 × 2 × 2 for monoclinic ZrO2 (P121/c1). TheMaterials Project ID [30] for each initial structure is listed inTable I. For molecular dynamics simulations, the time stepwas set to 1 fs. For each structure, 100-ps simulations atconstant volume and temperature (NVT ensemble) were per-formed. A snapshot was taken every 1 ps. For each snapshot,density functional perturbation theory (DFPT) was used tocalculate the Born effective charges [31]. The energy conver-gence criterion was set to 0.0003 eV for FPMD calculationsand 10−6 eV for DFPT calculations.For each structure, temperature values were 1300, 1500,1700, and 1900 K, and structures with strain values of −2%,−1%, 0%, 1%, and 2% were built. For structures with an oxy-gen vacancy, we considered +2 charged vacancies by addinga homogeneous background charge in VASP. In total, the train-ing dataset for the proposed MLPs contained 10 103 ZrO2snapshots. The details can be found in Table I.B. Neural networksIn this work, we develop neural network (NN) potentialsin order to describe not only the interatomic interactions,but also the impact of an applied electric field, followingthe methodology introduced by Shimizu and coworkers [24].Specifically, a conventional NN of the Behler-Parrinello type[17], hereby denoted by NNP, is dedicated to the calculationof energy and forces without an electric field, and a NN isdedicated to the prediction of Born effective charges for eachion and is denoted by BEC-NN.066001-2STUDY OF THE ION MOBILITY IN DEFECT-LADEN … PHYSICAL REVIEW MATERIALS 10, 066001 (2026)In the formalism of NNP, atomic forces can be calculatedby applying the chain ruleFj = − ∂U∂u j= − ∂U∂Gv∂Gv∂u j, (1)where U is the total energy and Gv are atom-centered sym-metry functions (ACSFs) introduced by Behler and Parrinello[17].The Born effective charge is a tensor that quantifies how themacroscopic polarization in a material changes in response toatomic displacements. It reflects the coupling between atomicmotion and an external electric field and is defined by [32].Z∗i j = �∂Pi∂u j= ∂Fi∂Ej, (2)where � is the cell volume, Pi is the macroscopic polarization,u j are the atomic coordinates, Fi are the atomic forces, andEi is the electric field. It is expressed in units of elementarycharge e. The subscripts i and j represent the x, y, and zdirections.The formalism of NNP was adapted to predict the Borneffective charge tensor by first replacing U by −�e Pi and Fjby Z∗i j and by preserving directional information in the inputsvia vector atomic fingerprints (VAF) [33,34],V 1,αi =∑jRαi jRi je−η(Ri j−Rs ) fc(Ri j ),V 2,αi = 21−ζ∑j∑k(Ri j + Rik )α{1 + cos(θi jk − θs)}ζ× e−η(Ri j +Rik2 −Rs )2fc(Ri j ) fc(Rik ), (3)where η and ζ are width parameters. Importantly, Ri j andRik are the atomic vectors of atom i with atom j and atomk, respectively. θi jk is the angle between atoms i, j, and kat vertex i. α represents the direction, either x, y, or z. Rsand θs define the peak positions and fc is the cutoff function,defined infc(Ri j ) ={12[cos(πRi jRc) + 1]if Ri j � Rc0 if Ri j > Rc, (4)where RC is the cutoff distance. In this work, the cutoff dis-tance is set to 7 Å. Here, we consider the Born effectivecharge along a specific direction, assuming an electric fieldalong the z direction, in other words, �E = Ez. In this case, theexpression of the relevant components is given asZ∗z j = �∂Pz∂V zv∂V zv∂u j. (5)Therefore, the present BEC-NN model requires only VAFwith α = z as its input. Assuming that the forces are a linearfunction of the electric field, the external forces acting on theions can be expressed as a function of the BEC, F BEC−NNj = Z∗z jEz. (6)Finally, the total forces are calculated as the sum of theexternal forces and the values obtained by the BEC-NN,F Totalj = F NNPj + F BEC−NNj . (7)FIG. 1. (a) Top view and (b) side view of the initial structure ofpristine ZrO2 (750 atoms/cell).This methodology enables accurate simulations of ion dynam-ics under an electric field.It is worth noting that in recent years, several groups havedeveloped methods to predict tensor quantities such as theBorn effective charge (BEC) tensor using equivariant graphneural networks [35]. Such a method has been developedby Kutana et al. to predict the BEC [36]. They adopted theDFPT data obtained in the present work as one of the testcases to demonstrate the validity of their method. Their modelhas a slightly better accuracy than ours, but at the same timedemonstrates the validity of our BEC-NN.C. Molecular dynamics simulations under an electric fieldMolecular dynamics simulations using the developed neu-ral networks were performed with LAMMPS [37] using theinterface developed by Shimizu and coworkers [24]. Calcula-tions were performed in the NPT ensemble using a time stepof 1 fs, where the temperature and pressure were controlledusing Nosé-Hoover style non-Hamiltonian equations of mo-tion proposed by Shinoda et al. [38,39]. For each moleculardynamics simulation, the mean-squared displacements (MSD)were calculated over 1 ns after an initial equilibration periodof 100 ps.We built a 5 × 5 × 5 supercell of pristine tetragonal ZrO2containing 750 atoms, as illustrated in Fig. 1. One oxygenatom was removed to build the supercell with a vacancy,which thus contains 749 atoms.FIG. 2. Comparison between DFT and NNP on (a) the energyand (b) the forces.066001-3ANH KHOA AUGUSTIN LU et al. PHYSICAL REVIEW MATERIALS 10, 066001 (2026)FIG. 3. Comparison between DFPT and BEC-NN on the Born effective charges along the z direction.The diffusion coefficient is calculated using the expres-sion proposed by Kemp et al. [40] for field-dependent ionmobilities, which result from a time-dependent contribu-tion from drift and a contribution from diffusion, as writtenin below Eq. (8). An example is presented in Fig. S10 withinthe Supplemental Material [41]. In the present work, we useEq. (8) in the direction with a nonzero electric field component(z), and below Eq. (9) in other directions,〈x2i〉 = (vd,it )2 + 2D∗i t . (8)FIG. 4. Hopping events, counted along three directions.066001-4STUDY OF THE ION MOBILITY IN DEFECT-LADEN … PHYSICAL REVIEW MATERIALS 10, 066001 (2026)TABLE II. Fitted diffusion coefficient (in cm2/s) along (top) xdirection (D∗x ) and (bottom) z direction (D∗z ) for oxygen ions. Forcases marked by *, no hopping was observed and the MSD curvehad a zero slope.D*x Pristine With defectT (K) No field 0.05 V/Å No field 0.05 V/Å1000 0.0∗ 0.0∗ 3.93 × 10−8 2.55 × 10−81300 0.0∗ 0.0∗ 9.37 × 10−8 7.47 × 10−81500 2.43 × 10−7 2.85 × 10−7 1.93 × 10−7 5.36 × 10−61700 2.49 × 10−6 4.68 × 10−6 3.93 × 10−6 4.51 × 10−61900 5.99 × 10−6 7.72 × 10−6 5.15 × 10−6 1.12 × 10−5D*z Pristine With defectT (K) No field 0.05 V/Å No field 0.05 V/Å1000 0.0∗ 0.0∗ 1.06 × 10−8 4.15 × 10−81300 0.0∗ 0.0∗ 3.52 × 10−8 1.61 × 10−71500 2.12 × 10−7 1.32 × 10−7 1.15 × 10−7 1.26 × 10−51700 2.87 × 10−6 1.52 × 10−5 3.86 × 10−6 1.27 × 10−51900 4.69 × 10−6 3.36 × 10−5 4.73 × 10−6 6.01 × 10−5In absence of an electric field in a given direction, the firstterm vanishes, and the equation reduces to Einstein’s relationin one dimension, 〈x2i〉 = 2D∗i t . (9)For comparison purpose, we also performed simulationsunder identical conditions using a Buckingham-Coulomb po-tential developed by Xu et al. [16].III. RESULTSA. Development of neural networks1. Conventional neural network potential (NNP)The conventional neural network potential was constructedusing a network architecture of 66 input nodes, two hiddenlayers with 10 nodes and one output node for each atomicspecies: [66-10-10-1]. The dataset is split with a 90%–10%repartition into a training set and a validation set. The par-ity plots are presented in Fig. 2. The root-mean-squarederror (RMSE) of the energy and the forces were evaluatedat 5.4 (5.1) meV/atom and 0.105 (0.099) eV/Å, respectively,for the training (validation) set.2. Neural network for Born effective charges (BEC-NN)The neural network for predicting Born effective charges(BEC-NN) was constructed with a similar architecture, butunlike the NNP, the descriptor is VAF. The neural networkconsists of two hidden layers with 20 nodes and one outputnode for each atomic species: [94-20-20-1]. The RMSE ofthe training and validation sets are, respectively, 0.0938 e and0.0936 e, leading to good accuracy, as shown in Fig. 3. Adiscussion of the performance of the BEC-NN with respect tothe number of input nodes for VAF can be found in Fig. S1 andTable S-I within the Supplemental Material [41]. The diagonalelement Z∗zz shows that the effective charges can strongly differfrom the formal charges, i.e., +4 for Zr and −2 for O, high-lighting the importance of predicting accurate charges in orderto appropriately describe the ion displacements in moleculardynamics simulations. (More details can be found in Fig. S2within the Supplemental Material [41].) Similar significantfluctuation of the Born effective charges were also seen in thecase of Li3PO4 [24].B. Molecular dynamics simulationsFor both pristine ZrO2 and defect-laden ZrO2, we per-formed molecular dynamics simulations at 1000 K, 1300 K,1500 K, 1700 K, and 1900 K, with an electric field of 0.0V/Å, 0.05 V/Å, and 0.1 V/Å along the z direction. Additionalcalculations with an electric field of 0.01 V/Å and 0.02 V/Åwere also performed for structures with defects. The intensityof the electric field used in our simulations is higher than inexperimental setups in order to accelerate the dynamics ofthe system. Experimental observations of flash sintering spantimescales on the order of seconds, which are inaccessible tomolecular dynamics simulations that typically operate in theFIG. 5. Comparison between MSD of O atoms in (a) pristine ZrO2 and (b) defect-laden ZrO2.066001-5ANH KHOA AUGUSTIN LU et al. PHYSICAL REVIEW MATERIALS 10, 066001 (2026)Å Å ÅÅ Å ÅFIG. 6. MSD of defect-laden ZrO2 under an electric field of 0.05 V/Å. (a) MSD for O atoms and (b) MSD for Zr atoms.nanosecond-to-microsecond regime. For each simulation, wefirst allowed the system to relax for the first 100 ps, then themean-squared displacement was monitored for 1 ns.First, we observed that over 1 ns, all structures remainedstable except for defect-laden ZrO2 at 1900 K and 0.1 V/Å,where void formation was observed. As such, a high electricfield is likely destructive, the discussion hereafter will focuson the results with an electric field up to 0.05 V/Å.1. ZrO2 without electric fieldFor pristine ZrO2 without an electric field, a higher temper-ature leads to an increase in the mean-squared displacement(MSD) of O atoms, especially at the highest temperatures(1700 K and 1900 K) (see Fig. S3 within the SupplementalMaterial [41]). In contrast, the MSD profiles of Zr atoms,while slightly increased, reach a plateau. Zr atoms thereforeform a rigid framework. The MSD profiles do not depend onthe direction (see Fig. S4 within the Supplemental Material[41]). At 1700 K, it was observed that after 136 ps, an O atommoved to an interstitial position (Oi) and created a Frenkeldefect. The position of the created vacancy (VO) was sub-sequently taken by a neighboring O atom. Snapshots of thisprocess are shown in Fig. S5 within the Supplemental Ma-terial [41]. After the occurrence of this initial defect, hoppingfrom O atoms to the defect position is facilitated, and the MSDfurther grows quickly. At 1900 K, this occurred immediatelyafter the simulation started. At 1500 K, the same phenomenonwas observed after 364 ps, and is illustrated in Fig. S6 withinthe Supplemental Material [41]. For this Frenkel defect, wecalculated a formation energy of 1.80 eV. Using this valueas a minimum theoretical estimate, considering that 1500 Kand 1700 K correspond to thermal energies of 0.129 eVand 0.146 eV, respectively, that there are 216 O atomsper atomic cell, and a standard atomic vibration frequencyof 1013/s, we calculated the expected shortest waiting timeτ to be around 500 ps and 100 ps, respectively, in line withour observed times of 364 ps and 136 ps. Unsurprisingly,this phenomenon was not observed at lower temperatures. Toinvestigate the hopping events more in detail, we monitoredatomic displacements every 5 ps and considered that if thedisplacement along one axis was more than 1.5 Å, then onehopping event has occurred and is accounted for; hoppingevents towards negative x (resp. y, z) and positive x (resp. y, z)are counted separately. The net count, obtained by calculatingthe difference between the number of positive and negativehopping events, represents the global trend. The results, dis-played in Fig. 4, show that the number of hopping eventsslowly increases at first, then accelerates after the first defect isgenerated. Note that only hopping of O atoms was found. Analternative visual representation of the hopping event countis also shown in Fig. S7a within the Supplemental Material066001-6STUDY OF THE ION MOBILITY IN DEFECT-LADEN … PHYSICAL REVIEW MATERIALS 10, 066001 (2026)FIG. 7. Hopping events in defect-laden ZrO2 at 1500 K under an electric field of 0.05 V/Å along the z direction.[41]. In all directions, the number of hopping events towardspositive and towards negative values tend to cancel each other,leading to low net hopping counts that mostly hover aroundthe neutral value.One oxygen atom is removed to create the initial structure,leading to a 749-atom structure. The presence of an oxygenvacancy creates a site for neighboring O atoms to hop into,which enhances diffusion, as shown in Table II. In compar-ison with pristine ZrO2, hopping is also observed at lowertemperatures, as shown in Fig. 5. As expected, the diffusioncoefficient is higher for the defect-laden structure compared topristine ZrO2 at low temperature, confirming that preexistingoxygen vacancies promote the diffusion of O atoms.2. Impact of an applied electric fieldNext, we investigate the impact of an applied electric fieldon defect-laden ZrO2. This is possible through the predictionof the Born effective charges with the BEC-NN. We consid-ered an electric field of 0.05 V/Å in the z direction. As canbe seen in the MSD plots shown in Fig. 6, the MSD alongthe z direction is greatly increased by the applied electricfield, therefore enhancing the diffusion of O ions. It is note-worthy that using a fixed formal charge model would largelyunderestimate this phenomenon, as discussed in Figs. S8 andS9 within the Supplemental Material [41].This is reflected in the hopping event counts as the nethopping along z drifts towards negative values over time,as shown in Fig. 7. An alternative visualization is providedin Fig. S7b within the Supplemental Material [41]. The dif-fusion coefficient along the z direction shows a significantincrease under an applied electric field compared to the caseswithout an electric field, as reported in Table II. This is high-lighted by the trajectory lines, as shown in Fig. 8 for the first400 ps at 1500 K. In pristine ZrO2, only small displacementsare visible [Fig. 8(a)]. When an oxygen vacancy is present,displacements increase in all directions compared to the pris-tine case, as the defect allows hopping to occur [Fig. 8(b)].Finally, when an electric field is present, the displacements arefurther increased, with hopping towards negative z (downwarddirection) favored, as negatively charged O ions tend to movein the direction opposite the applied electric field [Fig. 8(c)].At 1500 K and 1700 K, an increase in the diffusion coefficientin defect-laden ZrO2 is also observed in directions perpendic-ular to the applied electric field (Table II). This observation066001-7ANH KHOA AUGUSTIN LU et al. PHYSICAL REVIEW MATERIALS 10, 066001 (2026)FIG. 8. Trajectory lines (gray lines) over 400 ps for trajectories at 1500 K. Zr (resp. O) atoms are colored in green (resp. red). (a) Pristine,no field, (b) With O vacancy, no field, and (c) With O vacancy, 0.05 V/Å.may be associated with non-zero nondiagonal elements of theBorn effective charge tensor, as can be seen in Fig. 3 for Z∗zxand Z∗zy.The diffusion coefficient D∗, calculated from the MSD, isshown in Table II. Under no electric field, D∗ is calculatedusing a linear fit and Einstein’s relation. With an electric fieldalong the z direction, the diffusion coefficient along the xdirection is calculated by linear fit while the diffusion coeffi-cient along the z direction is calculated using Eq. (8). Moredetails including vd can be found in Table S-II within theSupplemental Material [41]. At low temperatures (1000 K),no significant changes are observed with an electric field. Athigher temperatures, the diffusion coefficient along the direc-tion of the applied field D∗z increases sharply. For instance,in Fig. 9, in defect-laden zirconia at 1500 K and 1700 K, anincrease in the diffusion coefficient with respect to the electricfield is observed, and this increase is higher in the direction ofthe electric field D∗z .C. Nudged elastic band (NEB) calculationsThe energy barrier for migration of O atoms from one siteto another along the z direction is evaluated via nudged elasticband calculations, considering a structure with one O vacancyand 14 intermediate structures between the initial and finalstructures. The energy barrier evaluated by DFT (0.67 eV)was consistent with the value for a +2 charged oxygen defectreported by Eichler (0.61 eV) [42]. The profile given by ourneural networks shows a barrier of height 0.58 eV, which isconsistent (within 0.1 eV) with DFT results (more details inFig. S11 within the Supplemental Material [41]). The BECpredicted by the BEC-NN is consistent with DFPT calcu-lations along the NEB images, showing good transferability(Fig. S12 within the Supplemental Material [41]). The calcu-lated BEC is typically around +5.0 for Zr and –2.7 for O, butcan reach +4.5 to 5.4 for Zr and –2.7 to –1.0 for O, which isa significant deviation from the nominal values of +4.0 for Zrand –2.0 for O, and the values vary with the atomic positionsof the ions. Therefore, the assumption of constant nominalcharges is insufficient for accurately determining the forcesarising from an applied electric field. We then calculated theenergy barrier for a larger 383-atom cell, which was foundto be higher at 0.65 eV. The energy barrier for the smallsystem without electric field is underestimated, likely owingto size effects from the small cell size. As can be seen inTable III, the energy barrier is lowered by the applied electricTABLE III. Calculated energy barriers for oxygen hopping alongthe z direction under varying electric field strengths, using thetrained neural network models. Forward and backward barriers cor-respond to transitions from initial to final and final to initial states,respectively.Electric field Forward Backward Barrier(V/Å) barrier (eV) barrier (eV) difference (eV)0.00 0.65 0.65 0.000.02 0.59 0.64 −0.050.05 0.49 0.62 −0.13066001-8STUDY OF THE ION MOBILITY IN DEFECT-LADEN … PHYSICAL REVIEW MATERIALS 10, 066001 (2026)(a)(b)FIG. 9. Evolution of the O ion diffusion coefficient of defect-laden zirconia with respect to the applied electric field at (a) 1500 Kand (b) 1700 K.field, from 0.65 eV without field to 0.59 eV at 0.02 V/Å and0.49 eV at 0.05 V/Å. We also calculated the energy barrierfor two alternative modes oriented at 45° in the x-y plane,but with a small displacement along z of ±0.54Å. The barrierwithout field is around 0.32–0.33 eV, but depending on thez component of the displacement, can be either reduced to0.25 eV at 0.05 V/Å (field-assisted) or increased to 0.43 eV at0.05 V/Å (field-opposed). More details on the energy profilescan be found in Figs. S13 and S14 within the SupplementalMaterial [41].IV. CONCLUDING REMARKSIn this study, we proposed a neural network potential in-cluding prediction of the Born effective charges to study theion dynamics in ZrO2 under an electric field, based on densityfunctional theory and density functional perturbation theorycalculations.Our simulations show that in pristine ZrO2, at 1500 K andabove, defects appear and drive the dynamics of oxygen ions.The prior presence of an oxygen vacancy allows bypassingthe initial time needed to observe defects and increases thediffusion coefficient of oxygen ions.Under an applied electric field along the z direction, a netmovement of oxygen ions along this direction is observed,which is reflected by the MSD profiles as well as the atomichopping counts. The diffusion coefficient along this directionis greatly enhanced by the presence of the electric field andcould be associated with the unusual mass transport behaviorsunder flash events. The enhanced oxygen diffusion along zunder an applied field suggests measurable ionic conductivity.While not computed in the present work, our model outputs(energies, forces, and Born effective charges) could be used toobtain Green-Kubo [43,44] or Nernst-Einstein estimates forthe conductivity in future work.Our findings demonstrate the value of a tailored neuralnetwork architecture and provide atomistic insight into defect-mediated ion transport in zirconia under an electric field. Suchprocesses are relevant to high-temperature and field-drivenphenomena that are known to influence the macroscopic be-havior of high-strength ceramics under direct current. Furtherinvestigations should explore the role of the stabilizing com-pounds (such as yttria) as well as other phases of zirconiaand their impact on the defect chemistry. Improvements tothe neural network architecture, for instance to include allcomponents of the Born effective charge tensor [36], shouldfurther improve the accuracy. In the future, more complexsystems requiring more configurations in datasets, for instanceincluding explicitly Y, will be studied.ACKNOWLEDGMENTSThis work was supported by the Japan Society for thePromotion of Science (JSPS) Kakenhi Grant No. 24K01284and the Japan Science and Technology Agency (JST) CRESTGrant No. JPMJCR1996. Some of the calculations used in thisstudy were performed using the computer facilities at ISSPSupercomputer Center and Information Technology Center,The University of Tokyo.DATA AVAILABILITYThe data that support the findings of this article are openlyavailable [45], embargo periods may apply.[1] J. R. Kelly and I. Denry, State of the art of zirconia for dentalapplications, Dent. Mater. 24, 289 (2008).[2] K. Matsui, H. Yoshida, and Y. Ikuhara, Review: Microstructure-development mechanism during sintering in polycrystallinezirconia, Int. Mater. Rev. 63, 375 (2018).[3] H. Masuda, K. Morita, M. Watanabe, T. Hara, H.Yoshida, and T. Ohmura, Ferroelastic and plasticbehaviors in pseudo-single crystal micropillars ofnontransformable tetragonal zirconia, Acta Mater. 203, 116471(2021).066001-9https://doi.org/10.1016/j.dental.2007.05.005https://doi.org/10.1080/09506608.2017.1402424https://doi.org/10.1016/j.actamat.2020.11.013ANH KHOA AUGUSTIN LU et al. PHYSICAL REVIEW MATERIALS 10, 066001 (2026)[4] D. Tejero-Martin, M. Bai, J. Mata, and T. Hussain, Evolutionof porosity in suspension thermal sprayed YSZ thermal barriercoatings through neutron scattering and image analysis tech-niques, J. Eur. Ceram. Soc. 41, 6035 (2021).[5] J. Chevalier, L. Gremillard, A. V. Virkar, and D. R. Clarke,The tetragonal–monoclinic transformation in zirconia: Lessonslearned and future trends, J. Am. Ceram. Soc. 92, 1901 (2009).[6] F. Kern and B. Osswald, Mechanical properties of an extremelytough 1.5 mol% yttria-stabilized zirconia material, Ceramics 7,1066 (2024).[7] R. Raj, Joule heating during flash-sintering, J. Eur. Ceram. Soc.32, 2293 (2012).[8] H. Yoshida and T. Yamamoto, Fundamentals and futureprospects of flash sintering of advanced ceramics, J. Jpn. Soc.Powder Powder Metall. 64, 523 (2017).[9] H. Yoshida and Y. Sasaki, Low temperature and high strainrate superplastic flow in structural ceramics induced by strongelectric-field, Scr. Mater. 146, 173 (2018).[10] Y. Sasaki, K. Morita, T. Yamamoto, K. Soga, H. Masuda, and H.Yoshida, Electric current dependence of plastic flow behaviorwith large tensile elongation in tetragonal zirconia polycrystalunder a DC field, Scr. Mater. 194, 113659 (2021).[11] K. Wang, G. Chen, Q. Wang, X. Fu, and W. Zhou, Unusualelectrode-dependent deformation of 3Y-TZP induced by weakelectric current in oxygen-lean atmosphere, Scr. Mater. 205,114220 (2021).[12] M. Cologna, B. Rashkova, and R. Raj, Flash sintering ofnanograin zirconia in <5 s at 850◦C, J. Am. Ceram. Soc. 93,3556 (2010).[13] H. Motomura, D. Tamao, K. Nambu, H. Masuda, and H.Yoshida, Athermal effect of flash event on high-temperatureplastic deformation in Y2O3-stabilized tetragonal ZrO2 poly-crystal, J. Eur. Ceram. Soc. 42, 5045 (2022).[14] A. Itoh, T. Tokunaga, A. Kodaira, H. Yoshida, and T.Yamamoto, Variation of photoluminescence intensity depend-ing on the timing of electric field application during isothermalflash sintering for 3 mol%Y2O3-ZrO2 polycrystal, Ceram. Int.48, 28712 (2022).[15] R. A. Buckingham, The classical equation of state of gaseoushelium, neon and argon, Proc. R. Soc. Lond. Ser. A 168, 264(1938).[16] W. Xu, A. Maksymenko, S. Hasan, J. J. Meléndez, and E.Olevsky, Effect of external electric field on diffusivity and flashsintering of 8YSZ: A molecular dynamics study, Acta Mater.206, 116596 (2021).[17] J. Behler and M. Parrinello, Generalized neural-network repre-sentation of high-dimensional potential-energy surfaces, Phys.Rev. Lett. 98, 146401 (2007).[18] L. Zhang, J. Han, H. Wang, R. Car, and E. Weinan, Deep poten-tial molecular dynamics: A scalable model with the accuracy ofquantum mechanics, Phys. Rev. Lett. 120, 143001 (2018).[19] K. T. Schütt, H. E. Sauceda, P.-J. Kindermans, A. Tkatchenko,K.-R. Müller, SchNet–A deep learning architecture formolecules and materials, J. Chem. Phys. 148, 241722 (2018).[20] S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa,M. Kornbluth, N. Molinari, T. E. Smidt, and B. Kozinsky,E(3)-equivariant graph neural networks for data-efficient andaccurate interatomic potentials, Nat. Commun. 13, 2453 (2022).[21] A. Musaelian, S. Batzner, A. Johansson, L. Sun, C. J Owen,M. Kornbluth, and B. Kozinsky, Learning local equivariantrepresentations for large-scale atomistic dynamics, Nat.Commun. 14, 579 (2023).[22] Y. Park, J. Kim, S. Hwang, and S. Han, Scalable parallelalgorithm for graph neural network interatomic potentials inmolecular dynamics simulations, J. Chem. Theory Comput. 20,4857 (2024).[23] L. Zhang, H. Wang, M. C. Muniz, A. Z. Panagiotopoulos,and R. Car, Weinan E., A deep potential model with long-range electrostatic interactions, J. Chem. Phys. 156, 124107(2022).[24] K. Shimizu, R. Otsuka, M. Hara, E. Minamitani, and S.Watanabe, Prediction of Born effective charges using neuralnetwork to study ion migration under electric fields: Applica-tions to crystalline and amorphous Li3Po4, Sci. Technol. Adv.Mater.: Methods 3, 2253135 (2023).[25] G. Kresse and J. Hafner, Ab initio molecular dynamics for liquidmetals, Phys. Rev. B 47, 558 (1993).[26] G. Kresse and J. Furthmüller, Efficient iterative schemes forab initio total-energy calculations using a plane-wave basis set,Phys. Rev. B 54, 11169 (1996).[27] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gra-dient approximation made simple, Phys. Rev. Lett. 77, 3865(1996).[28] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B50, 17953 (1994).[29] H. J. Monkhorst and J. D. Pack, Special points for Brillouin-zone integrations, Phys. Rev. B 13, 5188 (1976).[30] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S.Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, K. A.Persson, The Materials Project: A materials genome approachto accelerating materials innovation, APL Mater. 1, 011002(2013).[31] S. Baroni, P. Giannozzi, and A. Testa, Green’s-function ap-proach to linear response in solids, Phys. Rev. Lett. 58, 1861(1987).[32] X. Gonze and C. Lee, Dynamical matrices, Born effectivecharges, dielectric permittivity tensors, and interatomic forceconstants from density-functional perturbation theory, Phys.Rev. B 55, 10355 (1997).[33] V. Botu and R. Ramprasad, Adaptive machine learning frame-work to accelerate ab initio molecular dynamics, Int. J.Quantum Chem. 115, 1074 (2015).[34] W. Li and Y. Ando, Comparison of different machine learningmodels for the prediction of forces in copper and silicon diox-ide, Phys. Chem. Chem. Phys. 20, 30006 (2018).[35] S. Falletta, A. Cepellotti, A. Johansson, C. W. Tan, M. L.Descoteaux, A. Musaelian, C. J. Owen, and B. Kozinsky, Uni-fied differentiable learning of electric response, Nat. Commun.16, 4031 (2025).[36] A. Kutana, K. Shimizu, S. Watanabe, and R. Asahi, Represent-ing Born effective charges with equivariant graph convolutionalneural networks, Sci. Rep. 15, 16719 (2025).[37] S. Plimpton, Fast parallel algorithms for short-range moleculardynamics, J. Comput. Phys. 117, 1 (1995).[38] G. J. Martyna, D. J. Tobias, and M. L. Klein, Constant pres-sure molecular dynamics algorithms, J. Chem. Phys. 101, 4177(1994).[39] W. Shinoda, M. Shiga, and M. Mikami, Rapid estimation ofelastic constants by molecular dynamics simulation under con-stant stress, Phys. Rev. B 69, 134103 (2004).066001-10https://doi.org/10.1016/j.jeurceramsoc.2021.04.020https://doi.org/10.1111/j.1551-2916.2009.03278.xhttps://doi.org/10.3390/ceramics7030070https://doi.org/10.1016/j.jeurceramsoc.2012.02.030https://doi.org/10.2497/jjspm.64.523https://doi.org/10.1016/j.scriptamat.2017.11.042https://doi.org/10.1016/j.scriptamat.2020.113659https://doi.org/10.1016/j.scriptamat.2021.114220https://doi.org/10.1111/j.1551-2916.2010.04089.xhttps://doi.org/10.1016/j.jeurceramsoc.2022.04.055https://doi.org/10.1016/j.ceramint.2022.06.185https://doi.org/10.1098/rspa.1938.0173https://doi.org/10.1016/j.actamat.2020.116596https://doi.org/10.1103/PhysRevLett.98.146401https://doi.org/10.1103/PhysRevLett.120.143001https://doi.org/10.1063/1.5019779https://doi.org/10.1038/s41467-022-29939-5https://doi.org/10.1038/s41467-023-36329-yhttps://doi.org/10.1021/acs.jctc.4c00190https://doi.org/10.1063/5.0083669https://doi.org/10.1080/27660400.2023.2253135https://doi.org/10.1103/PhysRevB.47.558https://doi.org/10.1103/PhysRevB.54.11169https://doi.org/10.1103/PhysRevLett.77.3865https://doi.org/10.1103/PhysRevB.50.17953https://doi.org/10.1103/PhysRevB.13.5188https://doi.org/10.1063/1.4812323https://doi.org/10.1103/PhysRevLett.58.1861https://doi.org/10.1103/PhysRevB.55.10355https://doi.org/10.1002/qua.24836https://doi.org/10.1039/C8CP04508Ahttps://doi.org/10.1038/s41467-025-59304-1https://doi.org/10.1038/s41598-025-01250-5https://doi.org/10.1006/jcph.1995.1039https://doi.org/10.1063/1.467468https://doi.org/10.1103/PhysRevB.69.134103STUDY OF THE ION MOBILITY IN DEFECT-LADEN … PHYSICAL REVIEW MATERIALS 10, 066001 (2026)[40] D. Kemp, A. Tarancón, and R. A. De Souza, Recipes for su-perior ionic conductivities in thin-film ceria-based electrolytes,Phys. Chem. Chem. Phys. 24, 12926 (2022).[41] See Supplemental Material at http://link.aps.org/supplemental/10.1103/jcsd-dbl2 for details on the parameters used to fit theBEC-NN, the MSD without electric field, a visual represen-tation of the formation of a Frenkel defect, trajectory lineswithout electric field, MSD profiles under strong electric field,comparison with a fixed formal charge model, MSD along dif-ferent directions without electric field, and a table with all fitteddiffusion coefficient values.[42] A. Eichler, Tetragonal Y-doped zirconia: Structureand ion conductivity, Phys. Rev. B 64, 174103(2001).[43] M. S. Green, Markoff random processes and the sta-tistical mechanics of time-dependent phenomena. II. Ir-reversible processes in fluids, J. Chem. Phys. 22, 398(1954).[44] R. Kubo, Statistical-mechanical theory of irreversible pro-cesses. I, J. Phys. Soc. Jpn. 12, 570 (1957).[45] A. K. A. Lu, BEC-NN, URL, https://github.com/AugustinLu/BEC-NN.066001-11https://doi.org/10.1039/D2CP01335Ehttp://link.aps.org/supplemental/10.1103/jcsd-dbl2https://doi.org/10.1103/PhysRevB.64.174103https://doi.org/10.1063/1.1740082https://doi.org/10.1143/JPSJ.12.570https://github.com/AugustinLu/BEC-NN