Tandem bulk oxygen diffusion and surface reactions in reducible metal oxides control redox cycle dynamics

Domain: chemistry
Published: 2026-04-11
Editorial outcome: minor_revision

Abstract
The interplay between bulk oxygen diffusion and surface reactions in reducible metal oxides is key in heterogeneous catalysts, but direct measurements of oxygen mobility, transient kinetics, and in situ spectroscopies have been lacking. Here, we reveal complex dynamic behavior of ceria-zirconia by H2 using transient kinetics via mass spectrometry and in situ Raman and near-ambient pressure x-ray photoelectron spectroscopies. Molecular dynamics simulations with a machine learning potential delineate competitive oxygen diffusion mechanisms, with an optimal mobility at intermediate reductions. We expose a compensation between vacancy availability and lattice distortion at intermediate to high reductions and Frenkel defects at low reductions, underscoring a potential deficiency of 16O/18O exchange experiments in deducing oxygen mobility. Vacancies in proximity require electron localization on Ce atoms further away. The continuous replenishment of surface oxygen results in a varying reduction rate, with H2 dissociation being the rate-limiting step. Multiscale transient simulations, consistent with experiments, indicate catalysts of potentially spatially varying oxidation states. The approach is broadly applicable to reducible oxide materials.

Introduction
The versatile chemistry of reducible metal oxides enables their application in many thermochemical processes. Ceria-based materials are ubiquitous in heterogeneous catalysis, owing to their highly reactive surface oxygen (vacancies), tunable composition, and flexible environment-dependent surface chemistry.1 Their catalytic applications include CO and hydrocarbon oxidation, which leverage the reactive surface oxygen; redox reactions, such as water-gas shift and methane reforming, which leverage the oxygen vacancies; and organic reactions, which leverage the tunable surface acidic and basic sites.2,3 The high oxygen storage capacity and thermal stability of doped cerias uniquely qualify them as an oxygen conductor in solid oxide fuel cells and an oxygen (vacancy, VO)-storage medium for oxidation and reduction reactions in thermochemical looping in newer applications.4,5 The rapid temperature response of emerging electrified reactors facilitates transient processes, such as chemical-looping water-splitting and methane reforming, where lanthanide-doped cerias and Ce/Zr mixed oxides (CZO) have been demonstrated as active materials.6,7 Ultimately, the design and optimization of these novel reactors for hydrocarbon processing require an understanding of material property and reactivity across the atomic and continuum scales and in variable reaction regimes.

Redox reactions over ceria-based catalysts often follow (reverse) Mars-van Krevelen (MvK)-type mechanisms, involving reactant oxidation by lattice oxygen and the replinishment of oxygen vacancies by a soft oxidant. Despite the popularity of ceria-based catalysts, the kinetic behavior of MvK-type mechanisms, especially regarding the roles of surface and bulk oxygen, is a topic of debate.8 Although oxidation and reduction elementary steps involve both oxygen species, surface properties have typically been the focus of kinetic and mechanistic studies. The high sensitivity of ceria-based materials to oxidizing and reducing environmental conditions necessitates the use of in situ and operando techniques for surface and bulk characterization. Presently, the role of bulk oxygen is addressed through reactor studies in top-down approaches by measuring the rates of oxygen removal reactions at differing degrees of overall reduction and by isotope exchange experiments.9-12

Detailed modeling of bulk and surface properties of ceria-based materials eludes conventional computational methods. Although the high diffusivity of oxygen in ceria-based materials is well-known, modeling the mobility of bulk oxygen requires multitude of dynamical phenomena. At high temperatures, the distortion of the lattice structure due to reduction, the diffusion of oxygen in conjunction with polaron hopping, and the collective and/or off-lattice motion of atoms all influence bulk oxygen mobility.13-16 Although these effects can be captured by ab-initio molecular dynamics (AIMD), the (relatively) costly density functional theory (DFT)-based electronic structure calculations hinder the exploration of phenomena beyond the picosecond timescale. On the other hand, the limited accuracy of conventional DFT methods for the prediction of the energetics of surface redox reactions presents a further challenge for building kinetic models from first-principles.17,18

Here, we combine experimental kinetic and spectroscopic measurements, molecular simulations, and multiscale modeling to develop an atomistic understanding of ceria-based materials’ reduction. CZO was reduced while monitoring the extent of reduction and the possible origin of the removed oxygen. The reduction experiments were paired with in situ spectroscopic measurements to provide further insights into the reduction process. Molecular dynamics (MD) simulations based on a machine learning potential (MPLPs) trained on DFT calculations uncover nanosecond-scale spatiotemporal dynamics of lattice oxygen and unbiasedly address their mobility at experimentally relevant temperatures and degrees of overall reduction. Finally, joint experimental/computational transient kinetic studies resolve the mechanism of CZO reduction and quantify the roles of surface reaction and surface/bulk oxygen exchange.

Results
Kinetic and Spectroscopic Observations of CZO Reduction
The H2 temperature programmed reduction (TPR) spectrum of CZO has two peaks around 545 and 700 °C that are commonly associated with surface reduction and bulk reduction, respectively (Figure 1a).19 The 545 °C peak corresponds to 1125 µmol/g, so lattice O (Ob) is removed from at least the outer 66% mantle of the CZO particle, indicating that, unlike typical literature assignments, Ob diffusion from the bulk to the surface significantly contributes to the first peak. Higher temperatures enable further reduction; the peak at 700 °C accounts for the remaining 502 µmol/g of removed O. By 800 °C, a total of 1627 µmol/g of O has been removed, and the CZO is reduced to an extent of reduction (ξred, see Methods and Section S1) of 0.96. Clearly, labeling the two peaks as surface and bulk reduction is an oversimplification, as the majority of the CZO reduction occurs in the 545 °C peak. In fact, within the first peak, the quantity equivalent to all removable surface O (“first monolayer” in Figure 1a) has already been removed by 396 °C, long before the maximum.

To assess the kinetics of CZO reduction at various temperatures, we reduced the CZO isothermally by flowing H2 and monitored ξred (Figure 1b). At 500 °C, the H2O formation rate starts at a maximum and gradually declines (a sharp initial rise due to reactor startup was omitted here for brevity). This unexpected behavior reflects the combined dynamics of surface reaction with H2 and bulk O (Ob) diffusion, which replenishes surface vacancies (VO). While various surface and mass transfer effects can affect H2O formation, the H2O generated during reduction exceeds the estimated amount of removable surface O (Os), which is calculated from the particle surface area and O surface density (Eq. S2). This indicates that the particle surface would rapidly become O deficient without continuous replenishment from the bulk (Eq. S5). Specifically, at a constant H2 flow, one would expect the H2O rate to decrease continuously with time as Os declines. Upon removing O from the surface layer, the rate may drop further due to diffusion limitations, with a square root dependence on time, as indicated by textbook analysis of the diffusion equation. The rate normalized by the number of removable Os, 𝑁𝑠(𝑡), determined from an overall O balance and the H2O detected via mass spectrometry (MS), should reflect the intrinsic rate constant of O removal and remain static. However, the data (Figure 1b) indicates a continuous decline in the normalized rate of H2O formation (𝑘𝐻2𝑂) with time on stream or ξred. 𝑘𝐻2𝑂 can be calculated by assuming either surface reaction limitations or bulk diffusion limitations, but assuming diffusion limitations produces unphysical values of 𝑘𝐻2𝑂 (Methods and Section S1).

To probe the relation between bulk and surface oxygens, we characterized the reduction of CZO using in situ Raman spectroscopy, leveraging its moderate surface sensitivity. Upon reduction of the CZO, the longitudinal optic (2LO) mode at 1200 cm-1 decreases in intensity while the feature at 2110 cm-1, corresponding to an electronic transition of Ce3+, increases in intensity (Figure S2a).7,20-22 Interestingly, the 2LO overtone and Ce3+ features coexist when CZO is partially reduced, reflecting a balance between oxidized and reduced states. When the gas flow alternates between 10% H2/Ar and pure Ar, the CZO appears to reoxidize despite the lack of inlet oxidant flow (Figure S2b). This suggests that Ob diffusion from the bulk to the surface replenishes surface VO, and that the delay in oxidation observed by Raman may reflect diffusion limitations. However, the lasers in Raman spectroscopy may penetrate hundreds of nanometers deep into the sample, although this is dependent on factors including, but not limited to, sample condition and laser wavelength.23-25 A CZO particle is approximately 32 nm in diameter (see XRD and BET analyses in Section S1), so the Raman laser is capable of penetrating entire CZO particles and potentially lacking surface sensitivity. Additionally, CZO absorbance increases as it reduces, changing penetration depth and Raman scattering intensity at the affected vibrational modes (2110 cm-1) due to the resonance Raman effect (Figure S2c). A control experiment with inert SiO2 instead of CZO also shows that H2 remains in the cell during the Ar purge, and the residual H2 could explain the delay in oxidation observed by Raman (Figure S3). The observed dead volume and possible temperature inhomogeneities can make in situ spectroscopy cells unsuitable for direct comparison to the reactors used in kinetic studies; criteria for operando cells and designs were explored in other works.26-28 A preliminary kinetic analysis can help determine whether diffusion barriers are present and if they are discernible within the spatiotemporal limitations of Raman spectroscopy.

In analyses of particle diffusion processes, two extreme cases are often employed (see discussion in Section S1). In the fast diffusion limit, the mobile species diffuses through the bulk faster than the characteristic time scales of surface reactions to keep up with surface transformation, and the bulk concentration profile is flat. In the slow diffusion limit, surface reactions are much faster and instantly consume species arriving from the bulk, inducing significant bulk concentration gradients. Previous works have effectively demonstrated that Ob diffusion in CZO is faster than the surface reduction reaction, but they rely on indirect kinetic measurements without direct probing of Ob concentrations within CZO particles.9,10 Techniques such as electron energy loss spectroscopy29 or x-ray photoelectron spectroscopy (angle-resolved, depth profiling)30,31 can directly measure Ce3+ and O concentrations, but these characterizations are rarely performed at relevant reaction conditions and lack the spatial and temporal resolution to track Ob diffusion. Monitoring CZO reduction using near ambient pressure x-ray photoelectron spectroscopy (NAP-XPS) resulted in qualitatively similar results to in situ Raman spectroscopy; though, the results are not sufficient for diffusion analyses. The operating pressure of 1 mbar H2 is still far from typical reaction conditions, and although there is good time resolution, this set of data has no depth resolution, the addition of which would prevent effective time resolution (Figure S4). We instead turn to computational methods to fill in the gaps beneath the surface.

Molecular Simulations of Bulk Oxygen Mobility
The multiple possible mechanistic hypotheses regarding the role of Ob, which may form based on kinetic and spectroscopic observations, raise questions regarding its mobility under reaction conditions. To understand how temperature and VO concentration influence the self-diffusivity of O in the CZO bulk, we performed MD simulations at five ξred levels: 0.06, 0.13, 0.25, 0.38, and 0.50, corresponding to the presence of 4, 8, 16, 24, and 32 VO in the simulation cell, respectively (Methods). Simulations spanned temperatures of 1100-1500 °C to enable sufficient sampling of O diffusion rare events. Figure 2a presents the activation energies (Ea) and pre-exponential factors (Do) extracted from Arrhenius fits, 𝐷 = 𝐷𝑜 exp(−𝐸𝑎/𝑅𝑇), at each ξred. When ξred  0.25, Ea depends weakly on ξred, and Do is also relatively constant. The substantial increase of both Ea and Do at higher ξred is reminiscent of the compensation effect; though, overall, the variations indicate a complex relationship between ξred and diffusivity. These sharp nonlinearities also reveal a deficiency in using the rate of 16O/18O exchange as an indicator of Ob mobility: since the measurements must be made under oxidant-rich environments, the inferred diffusivity only corresponds to that in the low ξred regime where diffusion proceeds through Frenkel defects. To more clearly resolve the dependence of diffusivity on ξred at a fixed temperature, we performed additional simulations at intermediate reductions at 1200 °C. As shown in Figure 2b, the diffusivity exhibits a clear power-law dependence on ξred, with a best-fit in the form of 𝐷 = 1.2×10−6𝜉𝑟𝑒𝑑 0.72 𝑐𝑚2/𝑠 . It is important to note that this expression is temperature-dependent due to the variation of Do and Ea for each reduction. At 1200 oC, the behavior is autocatalytic: the more vacancies, the faster the diffusion.

Building on these observations, we used the complete set of Arrhenius parameters to construct a diffusivity heat map across the entire temperature range of interest (Figure S12) as well as the estimated diffusion timescales for a spherical CZO particle 30 nm in diameter (Figure 2c). The heat map reveals a maximum for our operating conditions of 600 °C at ξred = 0.25. Note that a volcano-type dependence of O diffusivity on the overall degree of reduction was recently reported based on a regression of the rate of re-oxidation reactions measured in the bulk diffusion-limited kinetic regime.32 The maximum at intermediate ξred most likely arises from competing effects. At low ξred (below ~0.1), diffusion is limited by the lack of vacant sites, while at higher ξred (above 0.35), lattice distortions increase the energy barrier for diffusion. This distortion is visible in the O-O radial distribution functions (gOO) and top-view structures shown in Figure 2d, where an increased density appears between the first and second coordination shells at ξred = 0.5 compared to the pristine case. The broadening of intermediate range order is accompanied by the expected diminished first peak, indicating a disruption of the well-defined local structure. Together, these features are consistent with increased lattice distortion supporting the notion that significant structural disorder occurs at higher vacancy concentrations. The resulting compensation effect between site availability and lattice distortion governs the mechanism of diffusion at intermediate to high ξred, where Frenkel defects play less of a role.

To further understand the mechanism of Ob diffusion, we analyzed the MD trajectories in detail, focusing on the dynamics of VO. We identified the initial positions of VO by locating vacant O sites (tetrahedral interstitials) in the initial frame. To construct a trajectory for VO, we track the O atoms that move into vacant sites, effectively treating VO as if they had moved in the opposite direction. The plot (Figure S13) shows the cross-correlation of VO‘s velocities as a function of pairwise distance, with distances binned into coordination shells that align with the major peaks in gOO. We use the absolute value of velocity components to capture co-directional motion along the same Cartesian axis. The correlated VO mobility is minimal in the first coordination shell (0-3.25 Å) and reaches its maximum in the third shell (5-6.75 Å). The third shell corresponds to intermediate VO separations, with every other lattice site occupied by O in each direction. For the second coordination shell (3.25-5 Å), a moderate correlation is observed at all correlation times, with a maximum at the initial time, indicating brief coordinated motion of closely spaced VO. The decay of the second shell correlation at short times indicates a tendency for VO to move away from one another. Overall, the correlated VO mobility corroborates a uniform distribution of VO in a maximally reduced system (Ce0.5Zr0.5O1.75), as clusters of VO are less energetically favorable. VO in proximity require electron localization on Ce atoms further away, due to the unfavorable nature of Ce2+ sites. These results indicate that enhanced diffusion occurs at intermediate VO concentrations, which favor spacing for VO migration.

Multimodal Kinetic and Mechanistic Investigations of CZO Reduction
With an understanding of the diffusivity of Ob across a broad range of temperatures and ξred, we return to explore the mechanism of CZO reduction. We first assess the appropriateness of calculating 𝑘𝐻2𝑂 in the fast diffusion limit, akin to the traditional turnover frequency definition of the catalytic rate. Using the radius of a CZO particle from BET analysis (16 nm) as the characteristic length scale and a diffusivity of 10-9 cm2/s (Figure S12), the diffusion timescale, 𝜏𝑑𝑖𝑓𝑓 = 𝐿2/𝐷, is ~0.003 s, while the reaction timescale, 𝜏𝑟𝑥𝑛 = 1/𝑘𝐻2𝑂(𝑡=0), is ~40 s, at least four orders of magnitude higher and further exacerbated as diffusivity increases with CZO reduction up to a maximum (Figure 2c and S12). The analysis of 𝑘𝐻2𝑂 in the fast-diffusion limit is well-justified, corroborating with the kinetic analysis evaluated independently from the molecular simulations (Section S1).

Reduction experiments conducted with different H2 partial pressures and at different temperatures reveal insights into reduction kinetics and surface species behavior. Based on the value of 𝑘𝐻2𝑂 at various ξred and their maxima (𝑘𝐻2𝑂,𝑚𝑎𝑥), H2 pressure has a reaction order of 0.88 ± 0.05 (Figure 3a and S5a). A reaction order close to one indicates H2 dissociation is likely the rate-determining step. On the other hand, Arrhenius relations constructed using 𝑘𝐻2𝑂 at various ξred and 𝑘𝐻2𝑂,𝑚𝑎𝑥 shows that the Ea of CZO reduction is 99 ± 8 kJ/mol (Figure 3b and S5b). Unexpectedly, 𝑘𝐻2𝑂 does not remain constant during reduction. There is an initial inflection in 𝑘𝐻2𝑂, which is more gradual at higher temperatures, followed by a maximum and a drop that is more pronounced at 450 °C.

To understand the mechanism of CZO reduction at the atomistic level, we studied the dissociation of H2 and the removal of surface O using DFT calculations. Specifically, we considered the (111) facet of CZO as it was previously determined to be the most stable one.33 Surface reduction initiates through the heterolytic dissociation of H2 across a surface Zr-O bond, facilitated by the stronger Lewis acidity of the surface Zr4+ cation compared to Ce4+ [See Figure S14 for site nomenclature and Figure S15 for intermediates (IS) and transition states (TS)]. Following dissociation, the H atop Zr can migrate to a nearby O, reducing two surface Ce4+ to Ce3+. Two migration pathways were explored; the first (direct) pathway was the migration to subsurface OCe,ss, which elevates OCe,ss close to the surface and forms a H-OCe,ss + H-OZr,s pair. This pair is strongly stabilized by hydrogen bonding, and its decomposition to H2O and OCe,ss occurs barrierlessly. Finally, the chemisorbed H2O desorbs to create an O vacancy, VO,Zr,s, completing the OZr,s +H2 → VO,Zr,s +H2O redox half-cycle. The alternative pathway is the migration of H to OZr,s, with a similar barrier as that of the direct pathway. Although the produced H-OZr,s + H-OZr,s pair does not involve H bonding, it is more stable than H-OCe,ss + H-OZr,s, which can be subsequently formed by crossing a small barrier.

To account for the influence of experimental conditions on the surface reduction energetics, we constructed a pressure-corrected free energy diagram (Figure 4a). Under a representative environment (T = 500 C, PH2 = 10−1 atm, and PH2O = 10−2 atm), the initial surface reduction of CZO clearly follows two macroscopic timescales, corresponding to the dissociation of H2 and the desorption of H2O. These timescales are controlled by the stability of the H migration TS, of the H-OZr,s pair product, of the TS of H2O desorption (see Section S4 and ref 34), and of the VO,Zr,s product. Further refinement of the energetics of key states through single point calculations at the HSE06+D3BJ level (red and blue markers in Figure 4 representing the use of 25% and 15% exact exchange, respectively) reveals that the initial H2 activation sequence is the most kinetically relevant, as the net free energy barrier relative to gaseous H2 is much higher than that from the H-OZr,s + H-OZr,s pair to the TS of H2O desorption.

Although the kinetics of the initial surface reduction can be deduced from the free energy diagram, the role of Ob is not yet clear. Since the self-diffusion of Ob occurs rapidly (Figure S12), the Ob concentration profile in the particle interior should be uniform, and the catalyst’s response should be controlled by surface reduction reactions.9,10 Indeed, considering a diffusivity of ~10-9 cm2/s and an O-O site spacing of ~2.5 Å, one estimates the O/V hopping frequency to be ~107 s-1. This effective single-event rate constant far exceeds that of H2 dissociation. To capture the multiscale (simultaneous surface and bulk) reduction of a CZO nanoparticle, we developed a model treating the surface and a uniform bulk as distinct phases in equilibrium (see Section S4, Figure S16 and Eq. S6-S22).

Isothermal kinetic simulations of a single CZO particle reveal that the reduction of CZO occurs at a timescale of minutes and depends positively on PH2 (Figure 5b and Figure S17). At T = 500 C, PH2 = 10−1 atm, and PH2O = 10−2 atm, the single-particle 𝑘𝐻2𝑂,𝑝 (𝑘𝐻2𝑂,𝑝, see Eq. S21 and S22) reaches a maximum of 0.05 s-1 after the rapid quasi-equilibration of H2O desorption (Figure S18). In line with the free energy pathway, H2 dissociation is clearly the rate-determining process: the elementary step is always far from quasi-equilibrium, and the maximum 𝑘𝐻2𝑂,𝑝 (𝑘𝐻2𝑂,𝑚𝑎𝑥,𝑝) depends linearly on PH2 (Figure S19), in agreement with experimental measurements. Appreciable formation of VO begins at ~10 s, where the removal of OZr,s and Ob occur simultaneously, in agreement with spectroscopic observations. As more VO form over the course of reduction, the coverage of HOZr,s gradually grows as it is controlled by PH2O and the re-adsorption of H2O on VO; nevertheless, it remains low (Figure 5b). A greater PH2 results in a decline in 𝑘𝐻2𝑂,𝑚𝑎𝑥,𝑝 (Figure S20). For the bimolecular H2 activation process (heterolytic dissociation to a homolytic product), as O is gradually removed, 𝑘𝐻2𝑂,𝑝 slowly declines as fewer OZr,s sites remain to react with H2. Despite the amount of O removed already exceeding ¼ of OZr,s sites (removable surface O; see Section S4) by 21 s, 𝑘𝐻2𝑂,𝑝 only slightly diminishes due Ob continuously replenishing OZr,s. Like at 500 C, 𝑘𝐻2𝑂,𝑝 at 450 and 600 C also quickly reaches 𝑘𝐻2𝑂,𝑚𝑎𝑥,𝑝 early on and only slowly decays (Figure 5c). At constant PH2 and PH2O, the coverage of HOZr,s decreases with temperature, as H2O re-adsorption becomes progressively endergonic (Figure 5c). Based on the calculated 𝑘𝐻2𝑂,𝑚𝑎𝑥,𝑝, we estimate an Ea of 101 kJ/mol (Figure S19b), in agreement with the experimentally measured 99 kJ/mol (Figure S5).

With a working model of isothermal reduction kinetics in hand, we revisit the temperature programmed reduction (TPR) of a single CZO particle. When the free energy change of the bulk/surface O exchange (GO) was set to 0 eV (same as that in isothermal reduction simulations), a single H2O formation peak was observed at 542 C, in agreement with experimentally observed 545 C (Figure S21). At the onset of H2O desorption, O removal occurs simultaneously on the surface and in the bulk, where the amount equivalent to all removable surface O is removed by 404 C, in agreement with the experimentally determined 396 C. The plateau observed in the experiment after the main peak is possibly due to a multimodal distribution of VO formation energies in realistic CZO particles;35 computational studies have shown that VO formation energies depend on the local Ce/Zr distribution.33,36-38 TPR simulations performed with GO = 0.25 eV results in a flatter peak at 773 C, and a combination of spectra computed with GO = 0 and 0.25 eV strongly resembles the experimentally observed one (Figure S22).

Although the kinetic trends observed in single-particle simulations agree with experimental observations, finer differences remain in the sub-unity order of H2, the slow 𝑘𝐻2𝑂 ramp-up at 600 C, and the slight deviation from Arrhenius-type behavior (Figure 3 and S5b). These observations may be due to the consumption of H2 and accumulation of H2O in the reactor, as a lower PH2 and a higher PH2O both negatively affect 𝑘𝐻2𝑂,𝑝. To examine effects beyond the kinetic limitations, we simulated the reduction of CZO in a packed bed reactor (see Section S4, Eq. S23 to S27, and Table S2). Isothermal simulations reveal that, although the reactor 𝑘𝐻2𝑂,𝑟 (𝑘𝑟𝐻2𝑂) reaches its peak at the beginning of reduction at 450 and 500 C, 𝑘𝑟𝐻2𝑂 grows slowly at 600 C (Figure 5d). Indeed, while the time required to reach the maximum 𝑘𝑟𝐻2𝑂 (𝑘𝑟𝐻2𝑂,𝑚𝑎𝑥) initially decreases up to ~500 C due to faster 𝑘𝐻2𝑂,𝑝 (Figure 5c,e), it sharply increases at 550 C. Axial concentration profiles at 600 C show that at short times (t  10 s), CZO at the reactor entrance quickly consumes H2, hindering the formation of H2O downstream (Figure S23). When 𝑘𝑟𝐻2𝑂 approaches 𝑘𝑟𝐻2𝑂,𝑚𝑎𝑥 (t ~ 300 s), the axial profile of 𝑘𝐻2𝑂,𝑝 becomes uniform, resulting in a near-linear profile of H2 and H2O concentrations. This behavior contrasts from that at 500 C: without overconsumption, 𝑘𝐻2𝑂,𝑝 quickly becomes uniform; though, the gradual conversion of H2 along the reactor still slows 𝑘𝐻2𝑂,𝑝. Ultimately, reactor effects decrease Ea with increasing temperature and lower the order of H2 (Figure S24).

Methods
Experimental Methods

Materials
Cerium zirconium oxide (CZO, 50 mol% Ce) was procured from Sigma Aldrich. H2 (5.0), O2 (4.7), He (5.0), and Ar (5.0) were procured from Keen Gas.

Characterization and Experiments
X-ray diffraction (XRD) patterns were measured using a Bruker D8 XRD with a Cu Kα source. Average crystallite sizes were estimated using the Scherrer equation. BET surface areas were measured using N2 at 77 K with a Micromeritics ASAP 2020 instrument. Fresh CZO samples were used as is. Pre-reduced CZO samples were reduced at 800 °C for 3 h in 100 mL/min of 10% H2/He, then passivated in 1% O2/N2 at room temperature for 1 h prior to XRD or BET analyses.

Temperature programmed reduction (TPR) was done by monitoring H2O formation with a MS. 500 mg of CZO were pelletized and sieved into 60-100 mesh particles and placed in a quartz tube using quartz wool. The TPR was conducted in a laboratory-scale reactor (Lindberg Blue M, Thermo Scientific) starting with pristine CZO and 100 mL/min of 10% H2/He flow. All outlet lines were heated to 110 °C to avoid H2O condensation.

Isothermal reduction was performed similarly to the TPR. 1.8 g of pristine CZO were loaded in the same tube and reactor then heated to 450-600 °C in 100 mL/min of 10% O2/He at 10 °C/min, purged with pure He for 30 min at the target temperature, then exposed to 5-20% H2/He. All flow rates were maintained at 100 mL/min total flow.

In situ Raman spectra were acquired in a Linkam CCR1000 stage using a HORIBA Raman Spectrometer with a green-line solid-state laser (532 nm). 40 mg of CZO or SiO2 were pretreated in-situ with 10% O2/Ar at 500 °C for 30 min. The samples were then purged with pure Ar then reduced in 10% H2/Ar at 500 °C for 40 min or alternated between flows of 5.5 min of 10% H2/Ar and 30 min of 100% Ar. All flowrates were maintained at 50 mL/min of total flow.

In situ UV-visible spectra were acquired in an in-house-built low-void-volume environmental cell using an Avantes ULS2048CL-EVO UV-vis spectrometer equipped with an AvaLight-XE light source and an adjustable UV/VIS/NIR high-temperature probe. BaSO4 was used as a reference material for the background spectrum of the sample absorbance spectra. The CZO sample was calcined in 10% O2/Ar for 30 min at 500 °C, purged with Ar, then reduced in 10% H2/Ar for 30 min at 500 °C.

Near-ambient pressure X-ray photoelectron spectroscopy (NAP-XPS) spectra were collected at Brookhaven National Laboratory using a differentially pumped hemispherical analyzer (SPECS Phoibos) equipped with a monochromatic Al source (1486.6 eV). The powder sample was pressed into a thin layer onto a copper sheet. The CZO samples were pre-treated in situ with 2 mbar of O2 at 300 °C for 20 min, then cooled to 50 °C before evacuating the chamber and repressurizing with 1 mbar of H2. After H2 pressure was stabilized, the CZO was heated to 500 °C at 20 °C/min, and spectra were collected every 5 min at 500 °C under 1 mbar of H2. The XPS spectra were analyzed using CasaXPS.46

Extent of Reduction
The extent of reduction (ξred) is defined as the fraction of removable oxygen atoms that have been removed. A maximum of 12.5% of the lattice oxygen (OL) in CZO can be removed by reduction, defined based on the conversion of all Ce4+ to Ce3+. As such, ξred is equal to zero in a pristine CZO and one when 12.5% of its oxygen is removed. This fraction applies to the CZO bulk and surface.

Computational Methods
Density Functional Theory (DFT) Calculations of Bulk Ce/Zr Mixed Oxide (CZO)
All DFT calculations of CZO were performed using the Vienna ab-initio Simulation Package (VASP), version 5.4.1.47-50 The exchange-correlation energy was calculated using the Perdew-Burke-Ernzerhof (PBE) functional.51 The DFT-D3 method with the Becke-Johnson damping function (D3BJ) was used to calculate the van der Waals dispersion energy.52,53 Hubbard-type on-site self-interaction corrections (DFT+U) were included using Dudarev’s method,54 where an effective value of U (Ueff) of 5.0 eV and 4.0 eV were applied to the Ce 4f and Zr 4d states, respectively.33,55,56 The core electrons were described using the projector augmented-wave method.57,58 In DFT calculations, a (2×2×2) supercell of the bulk CZO unit cell was used, where the Brillouin zone was sampled at the  point. The plane-wave kinetic energy cutoff was set to 520 eV. The electronic structure in each self-consistent field (SCF) cycle was considered converged when the difference in energies between consecutive SCF steps fell below 10-6 eV. Atomic positions in structural optimization calculations were considered converged when the Hellmann-Feynman forces on relaxed atoms fell below 0.05 eV/Å.

Construction of Machine-learning Potentials (MLPs) for Molecular Dynamics (MD) Simulations of Bulk CZO
We use the moment tensor potentials (MTPs) as part of version 2 of the MLIP software package.59,60 The potentials were configured using a cutoff radius of 4.5 Å and an MTP level of 12 (i.e., including all invariant moment-tensor basis functions with level ≤ 12) with 8 radial basis functions in 3 channels. Three different levels of oxide reduction were explored for training our MLP: pristine CZO, 1 oxygen vacancy (Ce0.5Zr0.5O1.97), and 4 oxygen vacancies (Ce0.5Zr0.5O1.88). For each level of reduction, 250 configurations were generated by randomly perturbing the structure from its relaxed geometry up to a maximum of 0.25 Å in the x, y, and z directions. This provides an initial training set for active learning with MD. Forces and energies for each structure were obtained from stationary single-point DFT calculations without geometric optimization.

The initially trained MTPs were refined through an active learning loop using MD simulations performed with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS).61 The temperature was controlled using a chain of Nosé-Hoover thermostats, with a set temperature of 1800 ºC and a step size of 1 fs.62-64 Training using the active-learning framework of the method implemented in MLIP-2 does not use a fixed training/test split; instead, model reliability is controlled by the D-optimality extrapolation grade. For more details, see Ref 65. For active learning, we collected structures from MD above an extrapolation grade of 2 and exited MD if the grade was above 10. DFT energies and forces were computed for new structures encountered during active learning and appended to the training set. The process continues in a loop until no structures exceed an extrapolation grade of 2 during a nanosecond-long trajectory. Finally, results of all three different levels of reduction, totaling 3635 structures, were used to train a single final MTP (see Figure S7 for accuracy and Figure S8 for sample extrapolation grades).

For the computation of oxygen self-diffusivities, we utilized the Einstein relation based on the mean squared displacement of the lattice oxygen atoms.66,67
𝐷 ≈ 1/6𝑡 〈|𝐫(𝑡)−𝐫(0)|2〉 (1)

MD calculations were performed in the NVT ensemble at the desired temperature for 10 ns of equilibration, followed by 40 ns of production in the NVE ensemble (Figure S10 and S11). A larger (4×4×4) supercell was used for improved sampling.

We compute a pairwise cross-correlation function for vacancy mobility, denoted as Cr(𝜏) , as shown in Eq. (2) A normalized form of this correlation is given in Eq. (3).
Cr(𝜏) = 1/𝑁𝑟 ∑∑〈|𝐯𝑖(0)∙𝐯𝑗(𝜏)|〉 (2)
𝑖 𝑗≠𝑖
𝐶̃𝑟(𝜏) = Cr(𝜏) / max (Cr(𝜏)) (3)

Here, 𝑁𝑟 is the number of vacancies in the coordination shell r, 𝐯𝑖 is the velocity of vacancy i, and the assignment of distances is based on the vacancy pair separation at correlation time 𝜏 = 0.

DFT Calculations of CZO Surface Reduction and Simulations of Reduction Kinetics
DFT calculations of CZO surface reduction were also performed using VASP version 5.4.1. The electronic energies were calculated using the PBE functional with the same corrections (PBE+U+D3BJ) as stated above. One electron wavefunctions were expanded using a basis set of plane waves with kinetic energy up to 500 eV. The (111) facet of CZO was simulated using a nine-layer-thick p(2×2) slab, where periodic images in the z direction were separated by a 20 Å thick vacuum layer (more details are provided in Section S3 and previous work33). The Brillouin zone was sampled at the  point. Unless otherwise stated, spin polarization was used in slab calculations. Dipole corrections were applied in the z direction for all slab calculations.68 The same convergence criteria were used for the SCF cycles and structural optimizations. During structural relaxation, atoms in the bottommost tri-layer were fixed to bulk positions. Structural optimizations of transition states (TS) were performed using the climbing image nudged elastic band (CI-NEB) and quasi-Newton algorithms.69,70 Vibrational frequency calculations were performed using the finite difference method with a step size of 0.015 Å.

To obtain more accurate surface redox energetics for the analysis of the free energy pathways and parameterization of kinetic models, single-point calculations were performed for selected structures using the HSE06 hybrid density functional with dispersion energy calculated by the D3BJ method (more details are provided in Section S3).71,72 These structured were deemed important from kinetic analysis (see main text). Free energies of intermediates (IS) and TS were calculated using standard statistical mechanical equations.73 Kinetic parameters for surface reactions were obtained using transition state theory (TST), while those for adsorption steps were calculated using collision theory.74 The multiscale kinetics of CZO reduction were evaluated at both the single-particle and the reactor scale (more details are provided in Section S4).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Comment by co-reviewer:\nThis work combines transient kinetics, in situ Raman spectroscopy, NAP-XPS, MLP-MD, and DFT to investigate bulk\noxygen diffusion in CZO and the associated surface reactions. The research topic is interesting; however, the overall novelty\nof the manuscript is insufficient, and the manuscript contains significant technical issues. I therefore cannot recommend\npublication in its current form.\n1. The description of Figure 1b as “reaches its peak and gradually declines” seems misleading; the rate appears to be\nmaximal essentially at the earliest measured time. More importantly, this behavior alone does not support the conclusion that\nthe reaction proceeds via bulk O replenishing the surface O, since it could also arise from different types of surfaces, and\nmass-transfer effects.\n2. As the authors note, in situ Raman spectroscopy has limitations, its penetration depth may cover the entire ~32 nm\nparticle, reduction-induced changes in absorbance and resonance effects can alter peak intensities, and residual H2 may\nremain during the Ar purge. Therefore, the observed delayed reoxidation cannot uniquely demonstrate that it is surface\noxidation. Given that the authors calculated the H2O formation rate at 500 °C, a systematic comparison of H2O formation-\nrate trends across temperatures would help strengthen the mechanistic interpretation. In addition, the equation S1 appears\ninconsistent with the fitted relationship in Figure 2d, and this needs to be clarified.\n3. The manuscript assumes a removable oxygen fraction of 12.5% at 800 °C. Please provide the original reference and/or\nsupporting experimental evidence.\n4. The oxygen diffusivity is obtained from the MTP. However, a force RMSE of 0.30 eV/Å raises concerns regarding the\nreliability of quantitative diffusion predictions. Since the subsequent fast-diffusion-limit assumption and the multiscale kinetic\nmodel depend on this diffusivity, independent validation and/or an uncertainty (sensitivity) analysis is required.\n5. The training set includes only three reduction levels. Although active learning is employed, the manuscript does not report\nstatistics of the extrapolation grade for each reduction level and temperature, yet the authors extrapolate the model to\nsystems beyond these reduction states, which does not seem well justified.\n6. The MD simulations were conducted at 1100–1500 °C, yet the manuscript reports diffusivity at lower temperatures. Are\nthese values obtained via extrapolation? If so, how is the validity of this extrapolation assessed?\nComment by reviewer:\nI agree with the previous assessment that the manuscript in its current form is too preliminary for publication in Nature\nCommunications. To further strengthen this conclusion, I would like to offer several additional points for the authors'\nconsideration:\n(1) The proposed oxygen diffusion pathways and their associated activation barriers lack validation through Density\nFunctional Theory (DFT) calculations, which is essential for mechanistic credibility.\n(2) The computational models and analysis focus predominantly on (111) surfaces. Given that the experimental system\ninvolves nanoparticles, which inherently expose multiple facets, the discussion should be extended to include other relevant\ncrystal surfaces.\n(3) A microkinetic analysis based on results is necessary to quantitatively validate the proposed reaction mechanism.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.", "reviewer_3": "The paper \"Tandem Bulk Oxygen Diffusion and Surface Reactions in Reducible Metal Oxides Control Redox\nCycle Dynamics\" describes a study of oxygen mobility and surface reactions in ceria–zirconia catalysts. Using transient kinetics, in\nsitu spectroscopies, and machine-learning–based molecular dynamics simulations, the authors uncover complex, reduction-\ndependent oxygen diffusion mechanisms. They find that oxygen mobility is highest at intermediate reduction levels due to a\nbalance between oxygen vacancies and lattice distortions, while different defect types dominate at low and high reductions,\nchallenging conclusions drawn from isotope exchange experiments. The results show that surface oxygen is continuously\nreplenished, making hydrogen dissociation the rate-limiting step, and suggest that catalysts can exhibit spatially varying\noxidation states. The paper addresses a relevant problem and reports a methodology that is broadly applicable to other\nreducible oxide materials. In this respect, the work presented is definitely at the level required for publication in Nature\nCommunications.\nGiven my background in electronic structure theory, I will restrict my comments to this aspect of the work. The authors have\nperformed molecular dynamics simulations using a machine-learning potential to investigate oxygen diffusion processes,\ninteractions with H₂, and water desorption leading to reduced samples. Particularly commendable is the use of state-of-the-\nart electronic structure methods, such as DFT+U and the hybrid HSE06 functional. This level of theory is essential, as the\nreduction of ceria involves the formation of localized Ce³⁺ polarons that cannot be properly described using semilocal\nfunctionals. Overall, the computational work is of high quality and provides a valuable and necessary complement to the\nexperimental observations.\nThere are, however, two aspects that warrant further consideration. First, the analysis of H₂ reaction and H₂O desorption is\nbased exclusively on the CZO (111) facet. This choice is reasonable given that the average particle size is approximately 30\nnm. Nonetheless, even particles of this size can exhibit surface defects and low-coordinated sites, whose influence\nbecomes increasingly important as particle size decreases. It has been shown that small zirconia nanoparticles behave\nquite differently from extended zirconia surfaces with respect to H₂ dissociation and water desorption. In particular, while\nheterolytic H₂ splitting (Zr–H⁻/O–H⁺) is energetically favored on extended surfaces, homolytic splitting, leading to the\nformation of two OH groups and the transfer of two electrons to Zr³⁺ sites, becomes preferred on zirconia\nnanoparticles. Similarly, H₂O desorption is favored over O₂ desorption on ZrO₂ nanoparticles, in contrast to extended surfaces (see\nhttps://doi.org/10.1021/acs.jpcc.6b05984 and https://doi.org/10.1021/acsomega.7b00799). In light of these findings, it would\nbe appropriate to include a discussion of the potential role of surface defects and nanoparticle size in shaping the observed\nkinetic behavior.\nThe second point is more exploratory in nature. The authors consider a wide range of reduction levels for the CZO\nnanoparticles, up to 50%. The discussion, however, focuses primarily on kinetics, particularly on the relative diffusion\nbarriers. It would be interesting to address the accompanying changes in the electronic structure. To what extent are the\nexcess electrons localized, and how does their degree of localization evolve with increasing reduction? What is the\nelectronic nature of highly reduced CZO particles? While one may expect strongly reduced ceria–zirconia to exhibit\nincreased electronic conductivity, it remains unclear whether the system retains an insulating or semiconductor with\npolaronic character, or whether it approaches metallic behavior. Given the level of theory employed, the existing\ncomputational results should enable the authors to at least qualitatively address this question.", "reviewer_4": "The manuscript presents a comprehensive and compelling study that integrates experiments and MLP-based MD\nsimulations to unravel the complex interplay between oxygen diffusion and surface reactions during the reduction of CZO.\nThe application of MLP-based MD overcomes the significant limitations of experimental techniques and conventional DFT\ncalculations. The key finding is that the optimal oxygen mobility occurs at intermediate reduction levels due to a trade-off\nbetween vacancy availability and lattice distortion, providing atomic insights to redox kinetics. In summary, this work is of\nhigh quality and well supported by data. However, this paper is probably published, but revision is needed.\n1. The force threshold for structural optimization is only 0.05 eV/Å in the present DFT study, while this value is usually\nsuggested as 0.01~0.02 eV/Å. Why did the authors use a lower convergence here? Will it bring any unconvergent energetic\nresults?\n2. More details of the MLP should be provided. For example, how to construct the training set and testing set? What is the\nnumber of training set? What kind of machine learning model was used? What is the parameter setting in the ML model?\n3. Figure S7 only demonstrates a general picture for energy and force accuracy. However, whether the MLP can accurately\ndescribe the energy profile of Vo diffusion identical to NEB result using DFT should be demonstrated. Other details of the\nMLP performance (e.g., R2 values of training set and testing set data, and so on) should also be provided.\n4. The active learning loop is terminated when the extrapolation grades of all structures are less than 2 during 1-ns\nsimulation, however, the MLP-based MD simulations are performed in NVT for 10 ns. How to guarantee the accuracy of MD\nresults beyond 1ns, is it possible to explore unknown structures during 10-ns NVT? In addition, the energy profiles of 40-ns\nNVE simulations should be extra provided.\n5. In the TPR of a single CZO particle, is the equilibrium between bulk and surface oxygen ( GO° = 0 eV) relative to the\nextent of reduction (ξred)?\n6. The formula D=Aξ_red^B is fitted only based on data at 1200 °C. It is not stated whether it is applicable to other\ntemperatures. It is recommended to supplement with a discussion of its temperature dependence or validation under broader\nconditions.\n7. The diffusion coefficient is fitted from mean squared displacement through the Einstein equation. Thus, the MSD plots\nshould be provided in SI to confirm the fitting quality.\n8. The values of activation energy in Figure 2a are fitted from Arrhenius equation. Are the fitted values consistent with the\ndiffusion barriers from NEB calculations?"}, "author_response": "We thank co-Reviewer #1 and Reviewer #2 for their collective feedback. Below, we have addressed their\ncomments point-by-point and edited the manuscript (MS) and supporting information (SI) accordingly.\n1. The description of Figure 1b as “reaches its peak and gradually declines” seems misleading; the rate\nappears to be maximal essentially at the earliest measured time. More importantly, this behavior alone\ndoes not support the conclusion that the reaction proceeds via bulk O replenishing the surface O, since it\ncould also arise from different types of surfaces, and mass-transfer effects.\nAuthors’ reply:\nWe thank the reviewer for this comment. The quoted text incorrectly describes Figure 1b, so the text has\nbeen changed. Regarding the point about bulk O replenishing surface O, it is true that surface and mass\ntransfer effects can (and likely do) affect H O formation. However, H O formation must occur at the\n2 2\nsurface, and the amount of O at the particle surface can be estimated from the particle surface area and O\nsurface density as described in Equation S2. If bulk O is not replenishing the surface O, the particle\nsurface would rapidly become O deficient (Equation S5). The amount of H O generated during reduction\n2\nexceeds the amount of removable O at the surface.\nOn page 4 of the MS, text was changed as follows:\n“At 500 °C, the H O formation rate starts at a maximum and gradually declines (a sharp initial rise due to\n2\nreactor startup was omitted here for brevity).”\n2. As the authors note, in situ Raman spectroscopy has limitations, its penetration depth may cover the\nentire ~32 nm particle, reduction-induced changes in absorbance and resonance effects can alter peak\nintensities, and residual H may remain during the Ar purge. Therefore, the observed delayed reoxidation\n2\ncannot uniquely demonstrate that it is surface oxidation. Given that the authors calculated the H O\n2\nformation rate at 500 °C, a systematic comparison of H O formation-rate trends across temperatures would\n2\nhelp strengthen the mechanistic interpretation. In addition, the equation S1 appears inconsistent with the\nfitted relationship in Figure 2d, and this needs to be clarified.\nAuthors’ reply:\nWe thank the reviewer for this comment. H O formation is temperature dependent and can be\n2\ncharacterized by an Arrhenius relation in Figure S5. Equation S1 has been corrected to match Figure S2d.\nOn page 4 of SI #1, Equation S1 was corrected:\n𝜉 −𝜉 = −0.00188𝑒'.))*∗(-./01!2-./01\")\n!"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have made substantial changes to improve the quality of the manuscript. However, placing critical information in\nthe Supporting Information, particularly in Section S1, hinders reader comprehension. I would like to recommend the\npublication of this paper in Nature Communications after they make relevant changes.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.", "reviewer_3": "I went through the answers given by the authors to the various comments raised by the referees. As far as I am concerned I\nconsider the answers well formulated and the paper suitable for publication in Nat. Commun.", "reviewer_4": "The authors have answered all my questions and modified the manuscript according to my comments."}, "author_response": "We thank co-Reviewer #1 and Reviewer #2 their collective feedback and their recommendation for the\npublication of our manuscript. In our revision, we have added a passage alongside the description of the\nexperimentally measured rate of H O formation to state its implications on the apparent oxygen removal\n2\nkinetics.\nOn page 4 of the manuscript:\n“…vacancies (V ). While various surface and mass transfer effects can affect H O formation, the H O\nO 2 2\ngenerated during reduction exceeds the estimated amount of removable surface O (O ), which is calculated\ns\nfrom the particle surface area and O surface density (Eq. S2). This indicates that the particle surface would\nrapidly become O deficient without continuous replenishment from the bulk (Eq. S5). Specifically, at a…”"}]

Decision letter