A statistical understanding of oxygen vacancies in distorted high-entropy perovskite oxides

Domain: biology
Published: 2026-03-31
Editorial outcome: major_revision

Abstract
High-entropy perovskite oxides have emerged as promising electrode materials for solid oxide electrolyzers. However, their compositional complexity makes the formation of oxygen vacancies, which influence properties such as oxygen ionic conductivity and thermal expansion, challenging to predict. Here, we experimentally measure changes in oxygen vacancy concentration for fourteen perovskite oxides with high and low-entropy A-site compositions, finding a dependence on cation size variance in addition to divalent cation fraction. Atomistic simulations using a machine-learned universal interatomic potential reveal cation size mismatches broaden a distribution of vacancy formation energies, shown through statistical thermodynamics to shift bulk formation thermodynamics. Treating oxygen vacancies statistically enables accurate predictions of oxygen vacancy formation compared to traditional models. Practically, increasing the size variance between A-site cations reduces the temperature sensitivity of oxygen vacancy concentrations, making it key for tuning critical properties. More broadly, this study demonstrates statistical treatment of oxygen vacancies is essential for understanding high-entropy perovskite oxides.

Introduction
Perovskite oxides (ABO3−δ), such as LSCF (La0.6Sr0.4Co0.2Fe0.8O3−δ), are widely utilized as anodes for high-temperature solid oxide electrolyzer cells (SOECs), which offer the potential for efficient hydrogen production using renewable energy and heat. However, the performance of SOECs is often limited by anode properties, including low oxygen ion conductivity and high thermal expansion [1]. The oxygen ion conductivity of perovskites relies on a high concentration of oxygen vacancies (δ) to enable efficient ion conduction via vacancy hopping [2]. The concentration of oxygen vacancies also affects other important properties, such as the thermal expansion coefficient, thermochemical stability, and electrical conductivity [3, 4]. As a result, understanding and careful control of oxygen vacancy behavior is critical for optimizing the performance of perovskite oxides in high-temperature SOECs.

Recently, high entropy perovskite oxides (HEPOs) (ABO3−δ), with five or more cations with nearly equal molar amounts in either A- or B-site [5], have attracted great attention for SOECs as they theoretically could have good high-temperature chemical and structural stability due to the high mixing entropy [2]. In addition, high-entropy perovskites were shown to introduce bond distortion [6] and local ordering [7], which could potentially influence the formation of oxygen vacancies. It has been demonstrated that substituting metal cations can change the oxygen vacancy concentration [3, 8, 9]. However, exactly how the oxygen vacancy concentration in HEPOs varies with metal cation composition and temperature remains elusive as the net energy needed to remove a lattice oxygen is significantly complicated by variations from different metal cations and configurational entropy.

Previous research on simple perovskites has provided valuable methodologies and insights for studying oxygen vacancy. For example, experimental studies vary concentrations of one or two metal cations systematically and measure oxygen vacancies using techniques, such as thermo-gravimetry analysis (TGA) [4, 10–14], Columbic titration [9, 15–17], X-ray photoelectron spectroscopy (XPS) [8], or iodometric titration [17]. Then, the vacancies’ formation thermodynamics can be determined by applying concentration data to established defect models [3, 16]. Those studies consistently reveal the importance of the concentration of 2+ cations in the A-site (X2+) for forming high oxygen vacancy concentrations [12, 16]. However, experimentally sweeping every metal cation concentration is impractical for high-entropy materials with five or more dopants, and no such experimental studies have been reported in the literature.

In comparison, there are a few computational studies that aim to predict high-entropy doping effects on the oxygen vacancy concentration by calculating the vacancy formation energy (Ev). Density functional theory (DFT) studies calculated Ev for simple [11] and high-entropy oxides [7, 8, 13] and found cation substitutions could tune Ev as well as where vacancies preferred to form. An empirical study found such DFT-calculated Ev values could be predicted from the oxide’s formation enthalpy, electronegativity, and band gap [18]. However, while it is possible to predict Ev for many high-entropy oxide compositions computationally, these Ev predictions have yet to consistently align with experimentally measured oxygen vacancy concentrations. A recent study by Park et al. [13] calculated Ev with DFT across many oxygen sites for one HEPO and three simpler oxides and found all four had significant deviations in Ev depending on their cation neighbors. To address this variability, Park proposed Ev values exist in a distribution, noting the potential implications of the mean and variance on equilibrium vacancy concentration. Xu et al. [7] also used a statistical approach in their DFT study of complex perovskites to predict oxygen vacancy formation, finding that vacancies preferentially form near cobalt sites.

Both studies suggest the statistics of Ev may be key to predicting and tuning oxygen vacancies, however, the importance of such statistical effects for predictive models is unclear. The mechanisms of Ev statistics has also remained largely unexplored for a wide range of high-entropy oxides compositions, both experimentally and computationally. This gap needs to be filled to provide principles or guiding models for the selection of cation dopants to tune oxygen vacancies, and the critical properties dependent on them, in complex perovskite oxides.

In this work, we establish a combined experimental, computational, and theoretical workflow to study the effect of high-entropy mixed A-sites in perovskite oxides on the concentration of oxygen vacancies. Experimentally, we measured the oxygen vacancy concentration of 11 high-entropy (HE) and 3 low-entropy (LE) perovskites as a function of temperature using TGA, and found an enthalpy-entropy compensation relation where HE samples exhibit both lower enthalpy and entropy of oxygen vacancy formation. We found that the oxygen vacancy concentrations in LSCF-based perovskite oxides are directly affected by two primary variables: the conventional 2+ cation molar ratio in the A-site (X2+), and the newly identified A-site ionic radius variance (σA). Computationally, our high-throughput simulations with a machine-learned interatomic potential show that a increased σA distorts oxygen–B-site bonds, increasing the variance in oxygen bonding energy (σEv). Our statistical thermodynamic analysis shows that this broadening in σEv reduces the vacancy enthalpy (∆Hf) and entropy (∆Sf) of formation, matching experimental trends. We further propose a statistical model of vacancy formation that accurately predicts the vacancy concentrations of LE and HE materials from simulation. Overall, our results highlight A-site size mismatches in the A-site (σA) as a practical lever for tuning oxygen vacancies and their temperature dependence in LSCF-based perovskite oxides. Theoretically, tuning σA via structural distortion or other means could have similar effects in other oxide systems, high-entropy or not.

Results
We study perovskites with LE (2-3 cations) and HE (5+ cations) A-sites taking the form A(Co0.2Fe0.8)O3−δ, derived from LSCF, the state of the art for SOEC anodes [1, 19]. A-sites consist of Lanthanides and Alkaline earth metals and span a set of common compositional descriptors, such as divalent metal cation concentration X2+ (Eq. 19), A-site cation size variance σA (Eq. 1), Goldschmidt tolerance factor t (Eq. 20), and mixing entropy ∆Smix (Eq. 21). We measure samples with both different and same A-site divalent (2+ cation) doping concentrations (X2+) to differentiate divalent doping from other A-site mixing effects.

All 14 perovskite oxides were synthesized using solution combustion synthesis followed by thermal annealing at 800°C and those oxide powders were further processed with high-speed ball-milling to break apart agglomerations (see Methods). Those oxides show a homogeneous distribution of elements, as illustrated by a representative sample HE-2 using scanning electron microscopy (SEM) paired with energy dispersive X-rays (EDX) (Supplementary Figure 1a). The structure of all 14 oxides was characterized with X-ray diffraction (XRD) and found to match the same perovskite phase (orthorhombic Pnma) (Supplementary Figure 1b). While all synthesized HE oxides in the study were single phase, predicting stable high-entropy compositions remains a challenge for other compositions [20–22].

The change in oxygen vacancy concentration at temperature T, ∆δT, was measured using TGA. The weight was stabilized at 300°C under ambient air, then increased from 500 to 1000°C with 1-2 hour holds every 100°C. Stabilized mass readings suggest the vacancy formation reaction was in equilibrium (Supplementary Figure 2). Changes in oxygen vacancy concentration, ∆δT, are offset by δ0 (δT = δ0 + ∆δT, where δ0 = δ(500°C)) to find vacancy concentrations δT. The oxygen vacancy concentrations are shown in Figure 1a for all samples. Values of δ0 are expected to be small and were fit using nonlinear least squares (Table I) to a reaction model described below (Eq. 7), a typical approach for similar studies [4, 10, 12]. The fitted δ0 values range from 0.0005-0.014 and are indeed small compared to the measured concentration increase ranging from 0.035-0.13 seen in Figure 1a. As such, errors in model-fit δ0 values are likely to have marginal effects on the overall vacancy formation trends. The measured δ(T) curve for LE-1 (equivalent to LSCF) in Figure 1a is in agreement with reported values for the similar composition (La0.6Sr0.4FeO3): δ=0.044 vs 0.035 in [23] at 900°C.

Upon close examination of Figure 1a, it can be seen that the HE samples (solid lines) tend to exhibit more oxygen vacancies than LE samples (dashed lines), particularly at temperatures ≤ 900°C. Furthermore, the curves of the HE samples are quasilinear, but the LE samples show more exponential growth of oxygen vacancy concentration with temperature. Even among HE samples, HE-3 and HE-4 have steeper curves than HE-7 or HE-5. Figure 1b shows the oxygen vacancy concentration at 800°C as a function of the mole fraction of 2+ cations in the A-site (X2+), and clearly, larger values of X2+ lead to more vacancies in general for LE and HE samples. This is because X2+ is directly tied to the formation of oxygen vacancies as explained in Section IIB. Despite the important role of X2+, there are still large unexplained variations in the concentration of oxygen vacancy between samples for the same X2+=0.4, indicating that other factors also play an important role. Figure 1c plots the curvature of Figure 1a (in units °C−1) against the A-site size variance σA. The HE samples have consistently smaller curvatures and thus more linear growth of oxygen vacancies. σA is given by Eq. 1 in percent where rA,i is the ionic radius of the ith A-site element out of N elements and rˆA is the average radius.

σA = 100 · sqrt( (sum(rA,i - rˆA)^2) / (sum(rA,i)^2) ) (1)

The importance of X2+ is well known in the literature and will be discussed in the context of defect reaction thermodynamics in the following section. In comparison, the importance of σA in oxygen vacancy thermodynamics is yet to be explored and will be the subject for the remainder of this paper.

This section aims to determine the ∆Gf, ∆Hf, and ∆Sf values for the oxygen formation reaction (Eq. 2) for both the LE and HE samples. In the reaction given by Eq. 2, OxO is lattice oxygen, h. are holes, and V..O are oxygen vacancies in Kroger-Vink notation. When the A-site of a perovskite is doped with a 2+ cation (X2+), its negative charge relative to other 3+ cations will be compensated by forming either a B-site hole (h.) or half of an oxygen vacancy (V..O) to maintain charge neutrality as expressed by Eq. 3. The holes can be thought of as 4+ oxidation states of B-site cations (Fe. or Co.) with a net positive defect charge relative to the typical 3+ state, as the 2+ oxidation state of B-site elements (Fe′, Co′) has negligible concentrations in oxidizing conditions [10, 16]. So long as Eq. 3 is relevant, the maximum concentration of oxygen vacancies is limited by hole charge availability such that: max([V..O])= 1/2[X′2+]. This relation assumes the following: 1) 2+ A-site doping is the dominant form of hole generation (i.e., the intrinsic charge carrier concentration is negligible); 2) only the B-site cations change oxidation states during the oxygen vacancy formation; and 3) holes are not differentiated by B-site element [16, 23, 24]. Moreover, as a perovskite unit cell has three oxygen sites, the sum of vacancies V..O and occupied sites OxO must equal three (Eq. 4).

OxO + 2h. ↔ V..O + 1/2 O2 (2)

[X′2+] = 2[V..O] + [h.] (3)

[OxO] + [V..O] = 3 (4)

The equilibrium constant, Kp,ox for the oxygen formation reaction (Eq. 2) can be expressed as Eq. 5 considering the concentration relations in Eqs. 3-4 and assuming a reference oxygen partial pressure (pO0) of 1 bar. Kp,ox is now solely a function of X2+, pO2, and δT (Eq. 5) which are known or measured with TGA (Figure 1). Similar or identical expressions for the oxygen vacancy equilibrium constant have been used in previous studies [4, 10, 12, 14]. The explicit relation of oxygen vacancies, δT, on the concentration of X2+ in this expression is consistent with that seen in Figure 1b. X2+ is critical because previous studies have shown oxygen vacancy concentration saturates at δ = 1/2X2+ due to polaron depletion until the system reaches highly reducing conditions (pO2 ≤ 10−10 atm for similar materials) [10, 12, 24].

Kp,ox(T) = ( [V..O] * pO2^1/2 ) / ( [OxO] * [h.]^2 ) = ( δT * pO2^1/2 ) / ( (3-δT) * ([X2+] - 2δT)^2 ) (5)

This equilibrium constant can also be expressed as a function of the Gibbs free energy of formation ∆Gf, or equivalently enthalpy of formation ∆Hf and entropy of formation ∆Sf of oxygen vacancy, the specific gas constant R, and temperature T (Eq. 6). The combination of Eqs. 5, and 6 result in Eq. 7 below.

Kp,ox(T) = e^(-∆Gf / RT) = e^(∆Sf / R - ∆Hf / RT) (6)

ln( (δT * pO2^1/2) / ((3-δT)([X2+] - 2δT)^2) ) = ∆Sf / R - ∆Hf / RT (7)

Now, the δT curves measured by TGA for all 14 samples (Figure 2a) can be fit to Eq. 7 using nonlinear least squares to extract ∆Hf, ∆Sf, and δ0 (R2 >0.99 for all samples) and the values are reported in Table I. The fitted values of ∆Hf and ∆Sf for LE-1 are similar to previously reported values for the similar oxide (La0.6Sr0.4FeO3): 90 kJ mol−1 vs 95 kJ mol−1 [23] and 90 kJ mol−1 [24] for ∆Hf and 56 J mol−1 K−1 vs 54 J mol−1 K−1 [23] and 58 J mol−1 K−1 [24] for ∆Sf, validating the methodology.

Notably, the model in Eq. 5 is only valid for oxides with extrinsic polaron holes that interact with oxygen vacancies, as is widely accepted for LSCF-like perovskite oxides [19]. This occurs when the divalent doping exceeds oxygen vacancy concentrations (δ ≤ 1/2X2+) which is typical in the oxidizing conditions used in this study (pO2 = 0.21 atm). Furthermore, this model does not differentiate polarons by B-site, in part because polarons prefer forming on iron over cobalt [19, 25] so long as iron sites are available. For lower iron versus cobalt compositions where [Fe] < X2+ polarons could overflow onto cobalt sites affecting oxygen vacancy dynamics.

This section aims to understand differences in ∆Gf(T) between HE and LE samples that lead to the oxygen vacancy concentration differences at various temperatures. The Van’t Hoff diagram (Figure 2a) shows ln Kp,ox (or −∆Gf/RT) vs 1/T (Eq. 6) exclusively for samples with same X2+ = 0.4. ∆Sf and ∆Hf and are the intercept and negative of the slope, respectively. LE and HE samples show different trends; HE samples exhibit a larger Kp,ox (smaller ∆Gf) at lower temperatures, indicating more oxygen vacancies. With increasing temperature, the Kp,ox of HE samples increases more slowly than that of LE samples due to their smaller ∆Hf. This is equivalent to the oxygen vacancy concentration of HE samples increasing linearly, instead of exponentially, with temperature (Figure 1a and 1c).

We conducted a correlation study between the measured ∆Hf and ∆Sf (Table I) and four compositional descriptors (X2+, t, σA, ∆Smix) across n=14 samples (see details in Methods) shown in Figure 2b. Such a correlation study has not been reported previously. First, the cross-correlation of the four compositional descriptors is limited, suggesting the chosen sample compositions effectively span a four-dimensional A-site design space. In addition, ∆Hf and ∆Sf are not strongly correlated with X2+, indicating the dependence seen in Figure 1b is encapsulated by the inclusion of X2+ in the equilibrium constant (Eq. 5). ∆Hf and ∆Sf show moderate correlations with ∆Smix and tolerance factor t, but strong correlation with the A-site size variance σA.

When our experimentally measured ∆Hf and ∆Sf (isolating X2+=0.4) are plotted against σA (Figures 2c and 2d), their clear correlations can be fit with empirical relations included in the plots. When we added ∆Hf and ∆Sf values reported by [14] and [24] for variants of LE-1, they nicely follow the trends. The ∆Hf relationship with σA in Figure 2c (R2=0.95) and ∆Sf relationship in Figure 2d (R2=0.83) are quasilinear or quadratic. Together, we can observe a positive correlation between ∆Hf and ∆Sf in Figure 2e. This empirical correlation between ∆Hf and ∆Sf within a family of materials, known as enthalpy-entropy compensation, has been documented for reactions in other systems, although there is no agreement on the underlying mechanism [26, 27]. Suggested mechanisms include solvent interactions, hidden reaction intermediates [26], and numerical artifacts from data fitting [28], the last of which should not apply as our ∆Gf values are distinct between LE and HE samples (Figure 2a).

A proposed relationship between oxygen vacancies and the A-site size variations is summarized in Figure 2f. Introducing mismatches in the A-site cation size induces distortion in the lattice [6] and strains oxygen bonds [13], which reduces ∆Hf and ∆Sf. In short, a larger σA (e.g. HE samples) leads to a flatter ∆Gf, and hence a slower increase in oxygen vacancy with temperature (Figure 1a). The next section dives into how tuning σA can flatten ∆Gf by lowering ∆Hf and ∆Sf.

The formation enthalpy ∆Hf of oxygen vacancies is a function of the energy to create an oxygen vacancy, Ev, which can be calculated with atomistic simulations.

∆Hf = f(Ev) (8)

Ev = Edefect - Ebulk + 1/2 EO2 (9)

A unique aspect of disordered high entropy materials is that the vacancy energies Ev vary across lattice sites, so we need to estimate ∆Hf and ∆Sf statistically using a distribution of vacancy energies g(Ev). Treating Ev statistically is unconventional and has only recently been considered [7, 13]. We accomplished this by sampling hundreds of oxygen sites in an atomistic simulation using a machine learning universal interatomic potential (MLUIP) provided by Matlantis [29]. Included in the Supplementary Information is a detailed comparison of the Matlantis potential against DFT for calculating oxygen vacancy energy in smaller supercells over a selection of relevant materials (Supplementary Section 4). To control for contributions from divalent doping, only the eight materials in Table II with X2+ = 0.4 are considered. All materials are modeled with an orthorhombic Pnma perovskite phase, which fits well with XRD data (Supplementary Figure 1). Large randomly populated perovskite supercells are generated in the Atomic Simulation Environment [30] to represent LE and HE samples (Supplementary Figure 3). The resulting simulation cell contains 1280 atoms (256 perovskite unit cells), large enough that even the sparsest A-site element is represented within 2% of its stated mole fraction.

Vacancy energy Ev will be related to the B-site oxygen bonding as it is typically stronger than A-site oxygen bonding and dominates the bonding energy in perovskites [31]. Intuitively, variations in A-site size will cause variations in B-site bonding (Figure 3a) [31]. Such octahedra distortions have been directly observed with scanning transmission electron microscopy (STEM) in high-entropy oxides [6, 32]. To quantify this, we first find the energy-minimized structures of the initial supercells of various HE and LE samples with no oxygen vacancies. After energy minimization, we record the bond lengths and bond angles of B-O bonds for all 768 oxygen atoms (see Methods). Variance in bond length is known to affect properties in high-entropy materials and has been used as a predictor of phase stability [20]. As shown in Figure 3b, a large variance in A-site cation size σA leads to higher root mean squared deviations (RMSD) for both the B-O bond lengths and the B-O-B bond angles. This suggests that a large σA leads to a broader distribution of B-site oxygen bonding geometries, hence a broader distribution of energy to create an oxygen vacancy.

The oxygen vacancy formation energy Ev is calculated using Eq. 9, where Edefect and Ebulk are the total bulk defect energy of a relaxed supercell before and after removing one oxygen atom, and EO2 is the energy from forming gaseous oxygen found to be −9.71eV with the Matlantis potential (compared to the DFT tabulated elemental phase reference −9.52 eV [33]). It should be noted that our calculated vacancy formation energy values range from 0 to 2.5 eV, much larger than the mean average error of 0.03 eV relative to DFT the MLUIP reports for disordered systems [29].

Calculated Ev values for HE samples have higher variance compared to LE samples (Table II). g(Ev) represents the probability distribution of a given oxygen site having vacancy formation energy Ev. Curves for g(Ev) in Figure 2c are estimated by sampling Ev for all oxygen sites in the supercell and smoothing the distribution using a kernel density estimation with bandwidth 0.03 eV. Figure 3c shows that the means Eˆv for the two HE and one LE samples are comparable, but the two HE samples have a higher variance, as we expected from Figure 3b. The longer tail of the distribution, where Ev <0.6eV suggests that more oxygen vacancies will be generated at lower temperatures for HE samples. Park et al. suggested a similar effect of high-variance g(Ev) distributions in their DFT study of disordered perovskites, although experimental verification of this model was inconclusive [13]. Notably, the average sampled vacancy energy Eˆv across all HE samples was comparable the average for LE samples (1.19 vs 1.18 eV) while their average vacancy energy standard deviations σEv diverged significantly (0.19 vs 0.27 eV) suggesting Ev variance is key to differentiating LE and HE vacancy behavior (Table II). Next, we will use statistical thermodynamics to understand how the broadening of g(Ev) in HE samples leads to observed differences in ∆Hf and ∆Sf.

Oxygen vacancies can be modeled as the exchange of lattice oxygen with an external reservoir at a fixed temperature, described by the grand canonical ensemble. Each site either binds to an atom or not, so they can be treated as ideal Fermions following Fermi-Dirac statistics. In this ensemble, the probability of site i being vacant, pi, is given by Eq. 10 where ϵi is the free energy of an oxygen vacancy on site i, µv is the chemical potential of the oxygen vacancies in the system, and kB is the Boltzmann constant. The oxygen vacancy fraction x, equivalent to the vacancy concentration δ over three, is the probability of vacancy averaged over N oxygen sites (Eq. 11).

pi = 1 / (1 + e^((ϵi - µv) / kBT)) (10)

x = δ / 3 = 1/N * sum(pi) (11)

This vacancy probability pi is key to the statistical calculation of ∆Hf and ∆Sf. Finding pi (Eq. 10) depends on the vacancy free energy at each site ϵi, which we will treat as a distribution with average free energy ϵˆ and variance σϵ^2 (Eq. 12) due to small perturbations from the mean.

σϵ^2 = 1/N * sum((ϵi - ϵˆ)^2) (12)

Given this problem setup, the objective of this derivation is to find the first-order effects of energy perturbations due to variance σϵ^2 on predicted values of ∆Hf and ∆Sf. The first-order effects are derived in detail in the Supplementary Section 3. In short, by expanding around first-order terms of σϵ^2, we inverted Eqs. 10-11 to arrive at an expression for oxygen vacancy chemical potential µv as a function of vacancy fraction x and σϵ^2 (Eq. 13). Eq. 13 conforms to that for the regular solution theory [34]: the correction to the ideal mixing given by the rightmost term, of scale σϵ^2/(kBT), is the first-order contribution of energy dispersion. When energy broadening σϵ is sufficiently large, this additional term can be a substantial correction which has not been reported before and can be prominent for high entropy materials.

µv = ϵˆ + kBT ln(x / (1-x)) + (x - 1/2) * (σϵ^2 / kBT) (13)

In equilibrium, the chemical potential of lattice oxygen equals that of oxygen in the environment. Constitutive assumptions outlined in the Supplementary Section 3 connect the free energy of vacancies ϵi and vacancy chemical potential µv to vacancy energy Ev and gaseous oxygen formation entropy SO2 (a tabulated value). We assume that variations in free energy are primarily due to variations in the vacancy energy, that is, σϵ^2 = Var[Ev] = σEv^2. We define the distribution of oxygen vacancy energies with mean Eˆv and variance σEv^2 as g(Eˆv, σEv^2). Equating chemical potentials and collecting terms results in expressions for ∆Hf and ∆Sf as they are used in the law of mass action (Eqs. 5-6). This results in Eqs. 14 and 15 which are valid for a small σEv and small vacancy fractions (x). Full details, including expressions valid for larger vacancy fractions (x), can be found in the Supplementary. Notably, this derivation is applicable to oxygen vacancies in any material, and it has the potential to be extendable to other types of defects with variable formation energies.

The first derived expression, Eq. 14, predicts oxygen vacancy formation enthalpy ∆Hf including the first-order effects from broadening g(Eˆv, σEv^2). This expression is a departure from the assumption of uniform vacancies, where ∆Hf = Eˆv [8, 11, 13, 35], and is relevant for disordered or high-entropy materials, where σEv > 0. This correction can be large: for σEv = 0.2eV, our typical value for HE materials, this correction reduces ∆Hf by 45 kJ mol−1 at 1000 K compared to experimental values of ∆Hf=90 kJ mol−1 for LE materials. At the limit of small vacancy concentrations, broadening g(Eˆv, σEv^2) will always be lower ∆Hf. This is an intuitive result because g(Eˆv, σEv^2) will have a larger tail in the low Ev region; more low Ev oxygen sites will be available for vacancies. This is consistent with Park’s numerical results that shifting σEv can shift the equilibrium oxygen vacancy concentration [13]. The analytical expression in Eq. 14 gives an interpretable starting point for researchers to harness disorder-induced energy broadening to tune defect thermodynamics.

∆Hf ≈ Eˆv - σEv^2 / (kBT) (14)

∆Sf ≈ 1/2 SO2 - σEv^2 / (2kBT^2) - SO,vib (15)

The second derived expression, Eq. 15, is a departure from common assumptions [7, 11, 13], and describes the effects of oxygen vacancy energy broadening on ∆Sf. Broadening g(Eˆv, σEv^2) always reduces the effective change in entropy. This is because introducing oxygen sites that favor vacancies (low Ev) increases the likelihood of configurations where those sites are vacant, lowering the configurational entropy. Skewing pi away from the uniform case, where pi = pˆv always lowers configurational entropy. In short, having preferred sites for oxygen vacancies results in less configurational entropy, an effect propagated to a smaller ∆Sf. This is similar to how high-entropy materials with preferred cation configurations have lower configurational entropy, reducing their phase stability [21, 22]. In Eq. 15, we also considered the phonon (vibrational) entropy contribution by an oxygen atom in the lattice, which was normally assumed to be small. Our estimated value of SO,vib based on simulations described in the Supplementary Section 6 is not small, but about 25 J mol−1 K−1 (or 3R) for both LE and HE materials.

We now have a theoretical basis for explaining the experimentally observed enthalpy-entropy correlation (compensation mechanism) for oxygen vacancies (Figure 2e). Eqs. 14 and 15, found by approximating the first-order effects of σEv from first-principles, correctly predict HE samples with higher σEv will exhibit lower ∆Hf and ∆Sf. In fact, Eqs. 14 and 15 can be combined in Eq. 16 to predict a ∆Hf - ∆Sf scaling relation when Eˆv, SO2, and SO,vib are constant.

∆Hf = 2T∆Sf + (Eˆv - TSO2 + 2TSO,vib) (16)

Eq. 16 suggests the scaling between ∆Hf vs ∆Sf from perturbing σEv is approximately 2T. This prediction can be compared with our experimental results because simulated Eˆv values are comparable between LE and HE samples (Table II). The slope of the line of best fit in Figure 2e, known as the compensation or isoequilibrium temperature [27], is found to be 1120K, a comparable order of magnitude to the predicted slope of 2T (1500-2500K).

The scaling relation between ∆Hf and ∆Sf has opposing effects on the Gibb’s of formation ∆Gf with the enthalpy term winning. For a small vacancy fraction x, broadening g(Eˆv, σEv^2) with always lower ∆Gf (Eq. 17). Respect to δ(300°C), with simulated values based on the current statistical approach with the full distribution g(Eˆv, σEv^2) (Figure 4e), and compared to a single average (Ev = Eˆv) approach (Figure 4f). Experimental data in Figure 4d are eight samples with (X2+ = 0.4) from Figure 1a. For both approaches, we used sampled values of Ev from the simulation in coordination with the unsimplified Eqs. 10-11, assuming SO2 = 205 J mol−1 K−1 [37] and SO,vib = 25 J mol−1 K−1. We added additional terms for effects of pO2 and the charge carrier holes in the Supplementary Section 5. Oxygen vacancy concentrations are found by solving Supplementary Eq. 35 numerically. This model expands one the model suggested by Park et al. [13] by coupling to vacancy formation to polaron destruction, considering oxygen vibrational entropy (SO,vib ̸= 0), and increasing site sampling with a machine-learned interatomic potential. Figures 4d-f clearly show that the statistical model better matches experimental values and sample trends. Low-temperature overestimates of LE samples (Figure 4e) could be due to insufficient activation energy to achieve equilibrium. Importantly, the statistical model accurately captures the oxygen formation behavior vs. temperature in terms of the curvatures of LE and HE materials. For example, the model correctly predicts HE-5 and HE-7 to exhibit the largest vacancy increases at low temperatures, with HE-7 beginning to plateau at high temperatures due to polaron hole depletion when δ → 1/2X2+ (see Section IIB). In contrast, the conventional oxygen vacancy prediction model (Figure 4f) fails to correctly differentiate LE and HE sample behavior.

This work provides the first broad experimental validation of the statistical vacancy model proposed by Park [13] and Xu [7], demonstrating that g(Eˆv, σEv^2) more accurately captures oxygen-vacancy behavior than a single mean Eˆv. Our results confirm that modeling variations between oxygen sites is crucial for predicting the unique behavior of oxygen vacancies in our high-entropy perovskite oxides. In addition, our theoretical analysis (Section IIE) provides physical insight into how vacancy-energy distributions influence oxygen vacancy formation dynamics. Finally, comparison across 14 A-site compositions identifies σA as a key composition-property descriptor for guiding future HEPO design.

Figure 4g schematically summarizes this paper’s findings. Across the 14 materials tested, two compositional parameters appear to dominate oxygen vacancy behavior: mole fraction of A-site divalent cations X2+ and A-site size variance σA. The A-site divalent cations mole fraction X2+ sets an upper limit for oxygen vacancy concentration in oxidizing environments as it provides available charge per Eq. 3 [10, 12]. The A-site size variations σA, as a structural descriptor, affect the vacancy energy variations σEv, which have derivable consequences on the vacancy thermodynamics. A larger σA increases σEv, which stretches the vacancy formation curve horizontally (on a δ − T plot) and results in a less temperature-dependent oxygen vacancy change (e.g., smaller slopes). This suggests that properties closely related to oxygen vacancy concentrations, such as thermal expansion, oxygen ionic conductivity, and electrical conductivity, will also have less temperature dependence with a larger σA.

In summary, we investigated the dependence of oxygen vacancy concentration on temperature for eleven HE and three LE, LSCF-based perovskites through a combined experimental, computational, and theoretical workflow. Our experimental results show that the HE perovskites generally have higher concentrations of oxygen vacancy, but their concentration increases relatively more slowly than that of LE samples. When those experimental data were used to extract the oxygen vacancy formation ∆Hf, ∆Sf, and Kp,ox, they show strong correlations with both the conventional divalent doping concentration X2+ and the newly identified A-site size variance σA.

To better understand this result, we conducted a computational investigation using a machine learning interatomic potential and found that A-site size variations distort oxygen-B-site bonds, adding variance to the oxygen bonding energy. We further derived oxygen vacancy formation thermodynamics with a statistical treatment of oxygen vacancies, described by the vacancy energy mean Eˆv and variance σEv^2. The statistical model enabled experimental estimates of σEv which align with simulations. The statistical model also predicted vacancy concentrations that align remarkably well with experiments compared to a conventional model. These results show that the A-site divalent cations mole fraction X2+ sets an upper limit for oxygen vacancy concentration, and the A-site size variation σA affects the vacancy concentration dependence on temperature. A larger σA increases σEv, which stretches the vacancy formation curve horizontally (on a δ − T plot) and results in a less temperature-dependent oxygen vacancy change (e.g., smaller slopes).

Practically, those results suggest that the A-site size variance could be a simple knob to control properties closely related to oxygen vacancy concentrations, such as thermal expansion, oxygen ionic conductivity, and electrical conductivity, especially their temperature-dependent behaviors. Although this study focuses on A-site cation substitution in LSCF-based perovskite oxides in oxidizing conditions, similar strategies for tuning the variance in formation energy may be theoretically extended to control defect-related properties in other high entropy materials.

Methods
The solution combustion synthesis begins with stoichiometric mixing of metal nitrate precursors from Thermo Fisher Scientific: La(NO3)3 (99.9%), Sr(NO3)2 (≥99%), Ca(NO3)2 (≥99%), Ba(NO3)2 (≥99%), Nd(NO3)3·6(H2O) (99.9%), Sm(NO3)3·6(H2O) (99.9%), Gd(NO3)3·6(H2O) (99.9%), Y(NO3)3·6(H2O) (99.9%), Fe(NO3)2· 9(H2O) (≥98%), and Co(NO3)2·6(H2O) (≥98%). Precursors were dissolved in 100mL of deionized water to form a 1M solution of precursor cations in a 500mL Pyrex beaker. Sucrose (Sigma-Aldrich Omnipur) is added as a fuel in a fuel-oxidizer ratio of 1.0 (1M sucrose). The solution is stirred with a magnetic stir rod and dehydrated at 50°C on a hot plate for >12 hours or until solution forms a gel. The rod is removed and temperature raised to 350°C to initiate the combustion reaction. When the gel had fully burned, about 30 minutes, the product is removed and ground with a mortar and pestle. The powders then are placed in an alumina crucible and annealed at 800°C for 90 minutes with a 4°C min−1 ramp rate. Powders collected from the solution combustion synthesis were processed with high-energy ball-milling. Oxide powders (∼10g) were dissolved in ˜30mL of hexane (99.9% Fischer Chemical) and placed in a 50mL stainless steel vial in a Rietsch Cryomill. The powder underwent 3 hours of room-temperature wet-milling at 30Hz. Powders were collected after drying the resulting slurry on a 50°C hot plate.

The crystal structure (see Supplementary Figure 1) of the synthesized samples was investigated by X-ray diffraction (XRD, PANalytical Empyrean) with a Cu source (Kα, 1.54˚A of the wavelength). The morphologies of synthesized powders were examined using scanning electron microscopy (SEM) on an FEI Magellan 400 XHR, operating at an accelerating voltage of 5 kV and a beam current of 25 pA. SEM/EDS was conducted using the same instrument with an accelerating voltage of 20 kV and a beam current of 3.2 nA.

TGA is a common approach to measuring oxygen vacancies [4, 10–12, 14] and was found to be in agreement with iodometric titration and columbic titration methods [17]. The TGA protocol starts with a 1.5 hour 900°C preheating step to ensure mass stability. The 50-70mg samples are cooled and reheated with 1-2 hour holds every 100°C from 500°C to 1000°C in a Setaram LABSYS evo TGA-DSC (see Supplementary Figure 2). Experiments were performed in air (pO2=0.21 atm) to be consistent with high pO2 present for SOEC air electrode. TGA mass changes at each temperature ∆mT (mg) relative to the initial mass m0 (mg) were used to calculate the oxygen non-stoichiometry changes ∆δT for a material with a theoretical molar mass Mw with Eq. 22.

∆δT = (∆mT / m0) * (Mw / 16) (22)

In this paper, the Matlantis PFP estimator v7.0.0, trained on data without the Hubbard U correction, is used to perform structural relaxation using the LBFGS algorithm. We conducted atomistic simulations using an MLUIP provided by Matlantis [41]. Its stated mean average error of 0.03eV for disordered systems is small compared to calculated vacancy formation energy values in the range of 0.15-1.4 eV [29]. The structural analysis determined B-site bonding to oxygen with 1.1˚A radius collision spheres and used a relaxation force tolerance of f <0.05eV˚A−1. For vacancy formation energy, the program iterates over removing 768 oxygen sites in a 1280 atom supercell (Supplementary Figure 3) randomly populated with cations in stoichiometric ratios with no initial oxygen vacancies.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript investigated the statistical nature of oxygen vacancy formation in high-entropy perovskite oxides (HEPO). Systematic HEPO samples were synthesized and measured to provide experimental evidence for this important feature in high-entropy oxides. The computational work introduced the distribution of oxygen vacancy formation energy to the formation entropy and enthalpy. Despite the claim that a combined experimental, computational, and theoretical workflow has been established for oxygen vacancies in A-site mixed HEPO, the manuscript reads like two separate papers. The authors never put experimental data and modeling data next to each other. This is a missed opportunity to directly calibrate the modeling results with experiments and learn from the comparison. Careful examination of the experimental data and modeling results showed many disagreements, requiring detailed discussions. Table 1 and Table 2 contain modeled oxygen vacancy formation energy, enthalpy, and entropy, but they are not compared, nor do they agree with each other. Fig 1a shows the measured oxygen vacancy concentration increasing with temperature, Fig 4c shows the modeled oxygen vacancy concentration vs temperature. The shapes of these curves do not agree with each other. Fig 2a shows the ln(K) with 1/T in experiments, Fig 4b shows the modeled data. But the shape of HE disagrees with each other. Fig 2e shows the oxygen vacancy formation enthalpy and entropy from experiments, Fig 3d shows the modeled data. However, the experiment showed a curved trend, and the modeling results showed a linear trend. With these many discrepancies, it is hard to believe the model had captured all the essential physics in the experiments. The authors should comment on the reasons for the agreement and disagreement, so further improvements can be made.", "reviewer_2": "The paper reports on work using experimental synthesis, atomic-scale modeling, and statistic thermodynamics to determine the role of variation in the cation sublattice on the formation of vacancies on the oxygen sublattice in a high-entropy oxide (HEO) system. The resulting model could be used to design HEO systems with a desired concentration of oxygen vacancies, making them more useful for applications where oxygen transport is important. - What are the noteworthy results? The most noteworthy result is the statistical thermodynamics component that predicts the experimentally observed enthalpy-entropy correlation for oxygen vacancies in experimentally relevant systems. The finding using atomistic simulations with machine learned potentials that having a fraction of oxygen sites that are strained lowers the enthalpy of formation for oxygen vacancies and decreases their configurational entropy is not surprising or especially new, although it does inform the model described above and provides important quantification of these effects using an inexpensive and accurate computational method. The authors do a nice job emphasizing that while averages between LE and HE compositions are similar (Figure 3c and 3e), the distribution of these vacancy formation energies can vary and impact resulting predictions of future HEOs. - Will the work be of significance to the field and related fields? How does it compare to the established literature? If the work is not original, please provide relevant references. Some aspects of the thermodynamic work around bond length and entropy is reminiscent of other published work (Sarker, P. et al. Nat Commun 9, 4980 (2018). https://doi.org/10.1038/s41467-018-07160-7, Pitike, K.C. et al. Chemistry of Materials 32(17), 7507-7515 (2020). https://doi.org/10.1021/acs.chemmater.0c02702, Sivak, J.T. et al. Phys. Rev. Lett. 134, 216101 (2025). https://doi.org/10.1103/PhysRevLett.134.216101) that are not cited here. These other works focus more on the single-phase formation of HEOs, and the quantitative prediction and understanding of oxygen vacancy formation in high-entropy oxide materials discussed in this paper will influence the field and the many group investigating HEOs. - Does the work support the conclusions and claims, or is additional evidence needed? The work supports the conclusions and claims made by the authors sufficiently. - Are there any flaws in the data analysis, interpretation and conclusions? - Do these prohibit publication or require revision? No significant flaws have been identified. - Is the methodology sound? Does the work meet the expected standards in your field? There are many machine-learned potentials available, and the authors did not develop the one used here. Most new ones such as the one used in this work provide reasonable results for work such as done here, but it is difficult to judge exactly. While the authors comment on the accuracy of the MLUIP from reference 25 with respect to the vacancy formation enthalpy values observed here, a strong argument would be made if a more explicit benchmarking comparison was performed for the compositions explored here. Obviously one could only do this with DFT with much smaller supercells, however it would strength the argument for these compositions explored in this paper, as the MAE from reference 25 explores a much wider chemistry (~45 elements). The thermodynamic and experimental work appear to be exemplary and in-line with expected standards in the field. - Is there enough detail provided in the methods for the work to be reproduced? While the authors show the importance of considering the B-O bond length in this study, they do not comment on the relaxation tolerances when using the MLUIP. This detail should be added to the text to ensure other scientists can reproduce their results. Looser/tighter relaxation tolerances could change the findings, therefore we find this information important to include explicitly. The raw data generated from the Matlantis MLUIP denoted as “computational supplemental information” supplied with this manuscript provided sufficient detail to reproduce the analysis. Other notes: Figure 2 is distracting just showing LE and HE. It would be nicer if a more scientifically descriptive labeling system was used. Additionally, the data markers in the legends are quite small (e.g., the blue star is difficult to see). In the Methods, the authors state that the Matlantis MLUIP is trained on data without Hubbard U corrections, yet Co- and Fe-oxide systems commonly require such corrections within DFT calculations. Could the authors comment on how this might affect their predictions from the MLUIP here?", "reviewer_3": "This study discusses the effect of the ionic radius of different A-site element species, in addition to the valence states of the acceptor elements, on the oxygen deficiency energy of perovskite oxides using experimental and computational approaches. It is understood that the presence of many elements at the A-site leads to a larger probability density distribution and smaller oxygen deficiency energy. These phenomena are expected to be observed as the number of acceptor elements at the A site increases. In other words, although the effect of multi-element substitution is described, it is not clear that this is a high entropy effect that shows specific characteristics in materials with five or more substituents. In the paper, the authors addressed that “ A computational investigation using a machine learning interatomic potential finds that these A-site size variations lead to distortion of the oxygen-B-site bonds, which modulate oxygen bonding energy.” This clearly explains the nonstoichiometry behavior of the perovskite-type oxides. The oxygen defects in perovskite oxides vary greatly depending on the B-site elements. Therefore, the phenomenon discussed in this study is a result of the oxygen defect energy of perovskite oxides with a composition of Co0.2Fe0.8 at the B site. In other words, the content of this study is only applicable to LSCF-based materials, whereas the title of this paper gives the impression that it is applicable to all high-entropy materials, which is misleading. It is assumed to some extent that the inclusion of multiple A-site elements will cause local structural changes in the perovskite structure due to the large ionic radius of the A-site elements, resulting in changes in oxygen deficiency energies. In the claims of this study, no new findings can be found regarding the high entropy effect in high entropy materials created by more than five substitutions. I believe Nature Communications top journal requires a certain new challenge and trial, such as a new concept, new experimental technique, drastic enhancement of properties, theoretical understandings of difficult problems, etc. In light of the criterion, the salient result of the manuscript is maybe a coating method proposal, and this work would be appreciated as such a study on the fabrication method of the cathode materials. Therefore, I do not recommend this manuscript publication in Nature Communications, and I would offer the authors to submit it to other related journals.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "RE: Nature Communications 1st Revision\nOriginal Manuscript ID: NCOMMS-25-34369\nTitle: \"A statistical understanding of oxygen vacancies in distorted high-entropy oxides\"\nAuthors: Adam Potter; Yifan Wang; Kiran Hamkins; Dongjae Kong; Yuzhe Li; Jian Qin;\nXiaolin Zheng\nResponse to Reviewers’ Comments: We thank all reviewers for their thorough and\ncritical reviews. We have addressed the comments below.\nOverview\nWe noted three primary suggestions from the reviewers and responded with the following\ncontributions:\nSuggestion 1: Further analyzing and comparing experimental and simulation\nresults\n• Response: Detailed analyses are included in the revised Figure 4 comparing\nexperiments and simulations along several fronts, with the new Section 2F\nproviding discussion. The new Supplemental Section S5 gives methods for\npredicting vacancy concentrations (results in Figure 4d-f) better than the\nconventional model. We estimate the average contribution of phonon entropy in\nthe new Supplemental Section S6.\nSuggestion 2: Benchmarking the universal interatomic potential against DFT\n• Response: We performed additional DFT calculations of vacancy energy for\nthree materials (LE and HE) and compared them with the universal interatomic\npotential in the new Supplemental Section 4.\nSuggestion 3: Clarifying the universality of results to other materials and the role,\nif any, of mixing entropy\n• Response: We clarified the scope of conclusions and any causal relationship\nbetween mixing entropy, A-site size variance, vacancy energy variance, and\nvacancy concentration through revisions to the title, abstract, introduction,\nSection 2E, Section 2F, and the conclusion.\nReviewer comments are in black, our response to reviewer comments will be in blue and\nrevised portions of the manuscript in red here and in the main text and supplemental.\nReviewer 1\nThis manuscript investigated the statistical nature of oxygen vacancy formation in high-\nentropy perovskite oxides (HEPO). Systematic HEPO samples were synthesized and\nmeasured to provide experimental evidence for this important feature in high-entropy\noxides. The computational work introduced the distribution of oxygen vacancy formation\nenergy to the formation entropy and enthalpy. Despite the claim that a combined\nexperimental, computational, and theoretical workflow has been established for oxygen\nvacancies in A-site mixed HEPO, the manuscript reads like two separate papers.\nComment 1: The authors never put experimental data and modeling data next to each\nother. This is a missed opportunity to directly calibrate the modeling results with\nexperiments and learn from the comparison. Careful examination of the experimental data\nand modeling results showed many disagreements, requiring detailed discussions.\nThis comment will be addressed together with all comments below.\nComment 1a: Table 1 and Table 2 contain modeled oxygen vacancy formation energy,\nenthalpy, and entropy, but they are not compared, nor do they agree with each other.\nComment 1b: Fig 1a shows the measured oxygen vacancy concentration increasing with\ntemperature, Fig 4c shows the modeled oxygen vacancy concentration vs temperature.\nThe shapes of these curves do not agree with each other.\nComment 1c: Fig 2a shows the ln(K) with 1/T in experiments, Fig 4b shows the modeled\ndata. But the shape of HE disagrees with each other.\nComment 1d: Fig 2e shows the oxygen vacancy formation enthalpy and entropy from\nexperiments, Fig 3d shows the modeled data. However, the experiment showed a curved\ntrend, and the modeling results showed a linear trend.\nWith these many discrepancies, it is hard to believe the model had captured all the\nessential physics in the experiments. The authors should comment on the reasons for the\nagreement and disagreement, so further improvements can be made.\nWe agree with the reviewer that direct comparison of experimental and\ncomputational/theoretical results will make a strong case for the proposed oxygen\nvacancy model. Specifically, the reviewer suggested direct comparison of the following\ntrends: formation enthalpy (ΔH) and entropy (ΔS) values, oxygen vacancy concentration\nf f\nas a function of temperature, log(K ) vs 1/T (Van’t Hoff trend), and ΔH versus ΔS\np,ox f f\ncompensation trends, and asked to explain the discrepancies in some of the trend curves'\nshapes. We agree with these comments and have contributed the following work to better\nintegrate the experimental and computational work and clarify discrepancies.\nFirst, we have added Section 2F, dedicated to directly comparing experimental and\ncomputational results and discussing discrepancies. This section discusses the newly\nrevised Figure 4. Previously, Figure 4 illustrated a numerical example based on our\nderived model, providing curves independent of the simulations described in Figure 3.\nHowever, this distinction was unclear and unexpected in the flow of the paper. We have\nreplaced much of this numerical example with plots directly comparing enthalpy (ΔH) and\nf\nentropy (ΔS) values, changes in oxygen vacancy concentration, and vacancy energy\nf\ndispersion (σ ). This directly addresses concerns from Reviewer 1.\nEv\nIt should be noted that the simulations presented in Figure 4 and Figure 3 have been\nrerun in response to comments from Reviewer 2. The only difference is that we used the\nupdated version of the machine-learning model and included the Hubbard U correction.\nSimulation predictions of ΔH and ΔS in Figure 4a-b show excellent agreement with\nf f\nexperiments. In Figure 4c, we introduce the idea of estimating σ with experimental fitting\nEv\nand demonstrate agreement with simulated values. Finally, using our new simulation\nresults and a vibrational entropy offset, we estimate in the Supplemental (Section S6) that\nwe achieve compelling predictions of oxygen vacancy concentrations (Figure 4d-f). Using\nthe full set of samples vacancy energy values (Figure 4e) matches experiments (Figure\n4d) far better than using an average vacancy energy (Figure 4e) as is done traditionally.\nWe also provide possible explanations for why predictions in Figure 4 may deviate from\nexperiments in Section 2F.\nSection 2F\nWe now compare the theoretically derived ∆H from Eq. 13 and ∆S from Eq. 14 with\nf f\nvalues extracted from experiments for all LE and HE samples (Figures 4a and b). Since\nthe experimental ˆE and σ values were extracted from the temperature range of 500-\nv Ev\n1000◦C, Table I), we used an intermediate temperature of 800◦C when using Eqs. 13 and\n14. We applied a Gaussian fit to g( ˆE , σ2 ) sampled data in simulation (Table II) to find\nv Ev\nvalues for ˆE and σ , and used S = 25 J/mol-K (or 3k ) estimated with simulations\nv Ev O,vib B\ndetailed in Supplemental Section 6. This value, averaged over multiple materials, is\nsimilar to a previous DFT- calculated S ≈ 2 for LaFeO at 800◦C [34]. Simulations\nO,vib kB 3\npredict ∆H values well (Figure 4a) with a strong correlation of R2 > 0.8. Calculated ∆S\nf f\nvalues (Figure 4b) have good quantitative agreement with experimental values, especially\nfor LE samples. Overall, Eqs. 13 and 14 exhibit good predictive power, while factors such\nas sensible enthalpy, vacancy-vacancy coupling effects, experimental uncertainty,\nvibrational entropy, and accuracy of the Matlantis potential could account for remaining\ndiscrepancies. To further validate the thermodynamically derived Eq. 17, we used this\nequation to extract the value of σ from experimental data and compared it with the\nEv\nsimulated σ in Table II (Figure 4c). Specifically, we use experimentally determined ∆δT\nEv\nvalues, a tabulated value for S [35], and S =25 J/mol-K (Supplemental Section 6) to\nO2 O,vib\nfit values of σ , ˆE , and δ . Figure 4c shows σ values fit to 17 closely match the\nEv v 0 Ev\nsimulation sampled values (Table II) and both values of σ increase with increasing σ ,\nEv A\nindicating that the A-site size variation causes a vacancy formation energy distribution.\nFinally, we compare experimentally measured oxygen vacancy concentration change,\n∆δ (Figure 4d), with respect to δ(300◦C), with simulated values based on the current\nT\nstatistical approach with the full distribution of vacancy energies from simulations\n(E =g( ˆE , σ2 )),(Figure 4e), and a classical average value (E = ˆE ) approach (Figure\nv v Ev v v\n4f). Experimental data in Figure 4d are eight samples of (X =0.4) from Figure 1a. For\n2+\nboth approaches, we used sampled values of Ev from the simulation in coordination with\nthe unsimplified Eqs. 9-10, assuming S =205 J/mol-K [35] and S =25 J/mol-K.\nO2 O,vib\nWe\nadded additional terms for effects of pO and electron hole charge coupling (Eq. 2) with\n2\nfurther details in the Supplementary Section 5. Oxygen vacancy concentrations are found\nby solving Eq. S35 numerically. Figures 4d-f clearly show that the statistical model not\nonly matches the experimental values well, and the low-temperature overestimates of LE\nsamples (Figure 4e could be due to insufficient activation energy to achieve equilibrium.\nImportantly, the statistical model accurately captures the oxygen formation behavior vs.\ntemperature in terms of the curvatures of LE and HE materials. For example, the model\ncorrectly predicts HE-5 and HE-7 to exhibit the largest vacancy increases at low\ntemperatures, with HE-7 beginning to plateau at high temperatures. In contrast, the\nconventional oxygen vacancy prediction model (Figure 4f) fails to correctly differentiate\nLE and HE sample behavior. Since Park [13] and Xu [7] introduced the importance of\nusing g( ˆE , σ2 ) instead ˆE to predict oxygen vacancies, this work is the first report to\nv Ev v\nuse experimental data to illustrate the benefits of the statistical model over the other. Our\nresults confirm that modeling variations between oxygen sites is crucial for predicting the\nunique behavior of oxygen vacancies in our high-entropy perovskite oxides. Figure 4g\nschematically summarizes this paper’s findings. Across the 14 materials tested, two\ncompositional parameters appear to dominate oxygen vacancy behavior: mole fraction of\nA-site divalent cations X and A-site size variance σ . The A-site divalent cations mole\n2+ A\nfraction X sets an upper limit for oxygen vacancy concentration in oxidizing\n2+\nenvironments as it provides available charge (Eq. 2). The A-site size variations σ , as a\nA\nstructural descriptor, affect the vacancy energy variations σ , which have derivable\nEv\nconsequences on the vacancy thermodynamics. A larger σ increases σ , which\nA Ev\nstretches the vacancy formation curve horizontally (on a δ − T plot) and results in a less\ntemperature- dependent oxygen vacancy change (e.g., smaller slopes). This suggests\nthat properties closely related to oxygen vacancy concentrations, such as thermal\nexpansion, oxygen ionic conductivity, and electrical conductivity, will also have less\ntemperature dependence with a larger σ .\nA\nReviewer 2\nThe paper reports on work using experimental synthesis, atomic-scale modeling, and\nstatistic thermodynamics to determine the role of variation in the cation sublattice on the\nformation of vacancies on the oxygen sublattice in a high-entropy oxide (HEO) system. The\nresulting model could be used to design HEO systems with a desired concentration of\noxygen vacancies, making them more useful for applications where oxygen transport is\nimportant.\n- What are the noteworthy results?\nThe most noteworthy result is the statistical thermodynamics component that predicts the\nexperimentally observed enthalpy-entropy correlation for oxygen vacancies in\nexperimentally relevant systems.\nThe finding using atomistic simulations with machine learned potentials that having a\nfraction of oxygen sites that are strained lowers the enthalpy of formation for oxygen\nvacancies and decreases their configurational entropy is not surprising or especially new,\nalthough it does inform the model described above and provides important quantification\nof these effects using an inexpensive and accurate computational method. The authors\ndo a nice job emphasizing that while averages between LE and HE compositions are\nsimilar (Figure 3c and 3e), the distribution of these vacancy formation energies can\nvary and impact resulting predictions of future HEOs.\n- Will the work be of significance to the field and related fields? How does it compare to\nthe established literature? If the work is not original, please provide relevant references.\nComment 1: Some aspects of the thermodynamic work around bond length and entropy\nis reminiscent of other published work (Sarker, P. et al. Nat Commun 9, 4980 (2018).\nhttps://doi.org/10.1038/s41467-018-07160-7, Pitike, K.C. et al. Chemistry of Materials\n32(17), 7507-7515 (2020). https://doi.org/10.1021/acs.chemmater.0c02702, Sivak, J.T.\net al. Phys. Rev. Lett. 134, 216101 (2025).\nhttps://doi.org/10.1103/PhysRevLett.134.216101) that are not cited here. These other\nworks focus more on the single-phase formation of HEOs, and the quantitative prediction\nand understanding of oxygen vacancy formation in high-entropy oxide materials\ndiscussed in this paper will influence the field and the many group investigating HEOs.\nWe thank the reviewer for pointing out relevant literature that could strengthen the context\nof our work. In response, we have referenced these studies at various points throughout\nthe paper.\nIn the context of synthesizing high-entropy materials (Section 2A)\n“While all synthesized HE oxides in the study were single phase, predicting stable high-\nentropy compositions remains a challenge for other compositions (Sivak, 2025) (Pitike,\n2020).”\nIn the context of using bond length (and angle) variance as a descriptor of structural\ndistortion (Section 2D)\n“Variance in bond length is known to affect properties in high-entropy materials and has\nbeen used as a predictor of phase stability (Sivak, 2025).”\nIn the context of how preferred configurations reduces the configurational entropy of\noxygen vacancies (Section 2E)\n“This is similar to how high-entropy materials with preferred cation configurations have\nlower configurational entropy, reducing their phase stability (Pitike, 2020) (Sarker, 2018).”\n- Does the work support the conclusions and claims, or is additional evidence needed?\nThe work supports the conclusions and claims made by the authors sufficiently.\n- Are there any flaws in the data analysis, interpretation and conclusions? - Do these\nprohibit publication or require revision?\nNo significant flaws have been identified.\n- Is the methodology sound? Does the work meet the expected standards in your field?\nComment 2: There are many machine-learned potentials available, and the authors did\nnot develop the one used here. Most new ones such as the one used in this work provide\nreasonable results for work such as done here, but it is difficult to judge exactly. While the\nauthors comment on the accuracy of the MLUIP from reference 25 with respect to the\nvacancy formation enthalpy values observed here, a strong argument would be made if\na more explicit benchmarking comparison was performed for the compositions explored\nhere. Obviously, one could only do this with DFT with much smaller supercells, however it\nwould strength the argument for these compositions explored in this paper, as the MAE\nfrom reference 25 explores a much wider chemistry (~45 elements). The thermodynamic\nand experimental work appear to be exemplary and in-line with expected standards in the\nfield.\nWe agree that a proper benchmarking of the MLUIP against DFT would strengthen the\nsimulation results. In response, we have performed a thorough benchmarking of the\nMLUIP’s ability to calculate oxygen vacancy energies with respect to DFT. As the\nreviewer mentioned, we had to use smaller 100 atom cells for DFT calculations in one LE\nand two HE materials chosen to span all relevant elements. We calculated the vacancy\nenergy of three oxygen atoms in each and benchmarked against the MLUIP. Overall, we\nfound the Matlantis potential suitable for studying oxygen vacancies in high-entropy\nmaterials in this study.\n“The results shown in Figure S4b show strong correlations between the DFT and\nMatlantis predicted oxygen vacancy energy over multiple sites for all materials. The R2\nscore, a measure of explained variance, was R2=0.89 for the two HE materials and R2\n=0.47 when LE-1 is included. The mean average error, MEA, is nontrivial across all\nmaterials (MEA=0.41 eV), much of this is due to a consistent underestimation of E for\nv\nLE-1. For just HE-5 and HE-9, the MAE=0.24 eV, much less than the 2.3 eV range of\nDFT-calculated vacancy energies. The Matlantis potential underestimates E for LE-1\nv\n(and possibly other LE materials) by about 0.5 eV compared to DFT. In effect, accounting\nfor this offset would shift the predicted oxygen vacancy formation curves for LE materials\nto higher temperatures (Figure 4e in main text), improving oxygen vacancy predictions\nrelative to experiments (Figure 4d). “ *\n*More of Section S4 in Supplemental Information\n- Is there enough detail provided in the methods for the work to be reproduced?\nComment 3: While the authors show the importance of considering the B-O bond length\nin this study, they do not comment on the relaxation tolerances when using the MLUIP.\nThis detail should be added to the text to ensure other scientists can reproduce their\nresults. Looser/tighter relaxation tolerances could change the findings, therefore we find\nthis information important to include explicitly.\nWe agree that this is an important parameter and have included details in the methods\n(Section 4D).\n“The structural analysis determined B-site bonding to oxygen with 1.1A radius collision\nspheres and used a relaxation force tolerance of 0.05 eV/A.”\nThe raw data generated from the Matlantis MLUIP denoted as “computational\nsupplemental information” supplied with this manuscript provided sufficient detail to\nreproduce the analysis.\nComment 4: Figure 2 is distracting just showing LE and HE. It would be nicer if a more\nscientifically descriptive labeling system was used.\nWe agree and have removed these labels from Figures 2c-e and Figures 3d, as we can\nsee how these labels can be distracting on scatter plots. For the schematic Figure 2f and\nFigure 3e we have replaced “LE” and “HE” labels with “Low σ ” and “High σ ” as they are\nA A\nmore specific and descriptive. This also helps clarify that the described phenomenon is\nnot inherently tied to materials with high cation mixing entropy (high entropy materials),\njust more likely in such materials.\nSee Figures 2c-f and 3d-e\nComment 5: Additionally, the data markers in the legends are quite small (e.g., the blue\nstar is difficult to see).\nWe have increased the marker size in the legend in Figures 2c-d.\nSee Figure 2c-d\nComment 6: In the Methods, the authors state that the Matlantis MLUIP is trained on\ndata without Hubbard U corrections, yet Co- and Fe-oxide systems commonly require\nsuch corrections within DFT calculations. Could the authors comment on how this might\naffect their predictions from the MLUIP here?\nWe agree that the Hubbard U correction is commonly used with these transition metals.\nIn response, we have recalculated all vacancy energy distributions using a new version\nof the machine-learning potential trained on DFT calculations, including the Hubbard\ncorrection. Figure 3-4 and Table 2 have been updated accordingly. DFT calculations for\npotential benchmarking (Comment 2) also use the Hubbard U correction.\nSee Figures 3-4 and Table 2 as well as numerical references in the text.\nReviewer 3\nThis study discusses the effect of the ionic radius of different A-site element species, in\naddition to the valence states of the acceptor elements, on the oxygen deficiency energy\nof perovskite oxides using experimental and computational approaches.\nComment 1: It is understood that the presence of many elements at the A-site leads to\na larger probability density distribution and smaller oxygen deficiency energy. These\nphenomena are expected to be observed as the number of acceptor elements at the A\nsite increases. In other words, although the effect of multi-element substitution is\ndescribed, it is not clear that this is a high entropy effect that shows specific characteristics\nin materials with five or more substituents.\nThis comment will be addressed in conjunction with Comment 3, as both regard concerns\nthat the mechanism described for oxygen vacancies is not specific to high-entropy\nmaterials.\nComment 2: In the paper, the authors addressed that “A computational investigation\nusing a machine learning interatomic potential finds that these A-site size variations lead\nto distortion of the oxygen-B-site bonds, which modulate oxygen bonding energy.”\nThis clearly explains the nonstoichiometry behavior of the perovskite-type oxides. The\noxygen defects in perovskite oxides vary greatly depending on the B-site elements.\nTherefore, the phenomenon discussed in this study is a result of the oxygen defect energy\nof perovskite oxides with a composition of Co0.2Fe0.8 at the B site. In other words, the\ncontent of this study is only applicable to LSCF-based materials, whereas the title of this\npaper gives the impression that it is applicable to all high-entropy materials, which is\nmisleading.\nThe reviewer notes the possibility that the described mechanism for oxygen vacancies\nmay not only be specific to perovskites, but also specific to the B-site composition used\nfor all materials in the study. This would mean the observed mechanism is limited to one\nfamily of materials, constraining the scope and impact of the paper.\nWe believe that the implications in this paper extend beyond LSCF-based materials to\nimpact the broader field of high-entropy oxides. We would like to clarify the scope of our\nconclusions both here and through revisions to the text. In our paper, we argue for two\npoints: (1) A-site cation substitution in LSCF-based perovskites distorts the lattice to\ncreate a distribution of oxygen bonding energies and (2) broadening this distribution has\nmeasurable effects on the statistics of the oxygen vacancy thermodynamics. Point #1 (A-\nsite size variance → oxygen bonding variance) is a material design strategy specific to\nperovskites. It is likely that A-site size variations affect oxygen bonding for perovskites\noutside those with B-site (Co Fe ), but testing this hypothesis is outside the scope of\n0.2 0.8\nthe paper. Point #2 (oxygen bonding variance → lower enthalpy/entropy of formation) is\na result that comes directly from statistical models that are not structure specific (Section\n2E and S3 in Supplemental). These models suggest that introducing bonding energy\nvariance in any materials has effects on their vacancy formation thermodynamics. When\nwe apply this structure-ambivalent model to simulations, we find close agreement with\nexperiments, shown in the new Figure 4. Point (2) could suggest a strategy for tuning\noxygen vacancies in a far broader group of materials.\nIn short, we observe trends in a specific family of materials (LSCF-based perovskites)\nthat suggest physics with the potential to influence defect-engineering in highly distorted\n(including high-entropy) materials.\nWe recognize that the current title implies a specific generalization to all high entropy\nmaterials in the paper, so we have changed the title to specify the perovskite system.\nTitle: “A statistical understanding of oxygen vacancies in distorted high-entropy perovskite\noxides”\nTo better define the scope and generalizability of our findings, we have revised the\nintroduction, discussion, and conclusion. All red excerpts are new revisions. In the\nintroduction, we refine the scope of our contributions to 1) noticing a trend in specific\nfamily of materials 2) developing a theoretical framework to explain and better predict\nsuch differences and 3) noting that the theoretical framework could be extendable to other\noxides.\nIntro: “Overall, our results highlight A-site size mismatches in the A-site (σ ) as a practical\nA\nlever for tuning oxygen vacancies and their temperature dependence in LSCF-based\nperovskite oxides. Theoretically, tuning σ via structural distortion or other means could\nEv\nhave similar effects in other oxide systems, high-entropy or not.”\nIn Section 2E, we note that this theoretical framework outlining the effects of vacancy\nenergy variance could be applicable to other systems, although it is not tested here.\n“Notably, this derivation is applicable to oxygen vacancies in any material, and it has the\npotential to be extendable to other types of defects with variable formation energies.”\nIn our conclusion, we further clarify in our Conclusion that varying A-site cation size is\nclearly to perovskites, but the larger theme of modulating σ has the potential to be widely\nEv\napplicable.\n“Practically, those results suggest that the A-site size variance is a simple knob to control\nproperties closely related to oxygen vacancy concentrations, such as thermal expansion,\noxygen ionic conductivity, and electrical conductivity, especially their temperature-\ndependent behaviors. Although this study focuses on A-site cation substitution in\nperovskites, similar strategies for tuning the variance in formation energy may be\ntheoretically extended to control defect-related properties in other high entropy materials.”\nComment 3: It is assumed to some extent that the inclusion of multiple A-site elements\nwill cause local structural changes in the perovskite structure due to the large ionic radius\nof the A-site elements, resulting in changes in oxygen deficiency energies. In the claims\nof this study, no new findings can be found regarding the high entropy effect in high\nentropy materials created by more than five substitutions.\nWe agree with this statement: we draw no specific connections between A-site mixing\nentropy and oxygen vacancies. In Figure 2b, we found a limited to no correlation between\nA-site mixing entropy and oxygen vacancy thermodynamics. Instead, we describe the\nconsequences of introducing a continuous distribution of oxygen vacancy energies, a\nscenario most relevant to disordered high-entropy materials.\nAs an example, consider a metal oxide with a small concentration of a cation dopant.\nThere are approximately two possible oxygen binding energies; the oxygen is either\nadjacent to the dopant or only atoms in the bulk. Thus, the oxygen vacancy energy\ndistribution primarily consists of two low-variance peaks, and because they are weighted\nexponentially, we typically just consider the lowest energy site (likely next to the dopant).\nIn contrast, a high-entropy system has so many combinations of cation neighbors that the\ndistribution is essentially continuous, so we can no longer just consider the lowest value.\nWhile this effect is not specific to high-entropy materials, it is most relevant to systems\nwith many dopants at higher concentrations as it results in more oxygen-neighbor\ncombinations, and thus more possible oxygen vacancy energies. The vacancy energy\nvariance, σ affects both our low and high entropy samples, but results in Table 2\nEv ,\nsuggest the effect is stronger in high-entropy samples (our model suggests effects on ΔH\nf\nand ΔS scale with σ 2). This explains the distinction between LE and HE samples\nf Ev\nobserved experimentally (Figures 1a, 2a). Our thermodynamic model does not\ndiscriminate by composition or structure, so it is possible that similar statistics apply to\nany material with randomized distortion, whether it be technically high-entropy or not.\nWe recognize the discussion in our paper repeatedly differentiates samples as “low-\nentropy\" versus “high-entropy\", thus inviting an incorrect connection to A-site mixing\nentropy. In addition, we failed to clarify that the relevance to high entropy oxides stems\nfrom high-entropy oxides tending to have more structural (and energetic) distortion. These\nare key distinctions we must better communicate. We have made the following revisions\nto the abstract, conclusion, and discussion section 3C. In addition, we have also revised\nour “LE” vs “HE” labeling system.\nWe found the “LE” and “HE” labeling conventions misleading, so we have removed them\nin Figures 2c-e and 3d. The schematics in Figures 2f and 3e have also been relabeled to\nuse “Low σ ” and “High σ ” to avoid association with the formal definition of “high-entropy”\nA A\nsystems.\nSee Figures 2-3\nIn the Introduction, we clarify that the phenomenon is not specific to high-entropy\nmaterials.\nTheoretically, tuning σ via structural distortion or other means could have similar effects\nEv\nin other oxide systems, high-entropy or not.”\nIn Section 2F, and throughout the manuscript, we attribute the effect to structural\ndistortions, not a high-entropy effect.\n“The A-site size variations σ , as a structural descriptor, affect the vacancy energy\nA\nvariations σ , which have derivable consequences on the vacancy thermodynamics. A\nEv\nlarger σ increases σ , which stretches the vacancy formation curve horizontally (on a δ\nA Ev\n− T plot) and results in a less temperature-dependent oxygen vacancy change (e.g.,\nsmaller slopes). This suggests that properties closely related to oxygen vacancy\nconcentrations, such as thermal expansion, oxygen ionic conductivity, and electrical\nconductivity, will also have less temperature dependence with a larger σ .”\nA\nIn Section 2F, we note that the described phenomenon is necessary to explain differences\nbetween LE and HE samples, while not specific to HE samples.\n“Our results confirm that modeling variations between oxygen sites is crucial for predicting\nthe unique behavior of oxygen vacancies in our high-entropy perovskite oxides.”\nWe also revise conclusions to be applicable to “complex perovskite oxides” instead of\n“high-entropy materials”. The end of our revised abstract reads:\n“Altogether, this study not only shows that the A-site cation size variance is an important\nparameter in tuning properties of HEPOs, but also demonstrates that statistical treatment\nof the vacancy thermodynamics is critical for understanding high entropy materials.”\nI believe Nature Communications top journal requires a certain new challenge and trial,\nsuch as a new concept, new experimental technique, drastic enhancement of properties,\ntheoretical understandings of difficult problems, etc. In light of the criterion, the salient\nresult of the manuscript is maybe a coating method proposal, and this work would be\nappreciated as such a study on the fabrication method of the cathode materials. Therefore,\nI do not recommend this manuscript publication in Nature Communications, and I would\noffer the authors to submit it to other related journals.\nIn this major revision of the manuscript, we have rerun all simulations, directly compared\nexperimental and simulation results, benchmarked our potential against DFT, and found\nan appropriate model to accurately predict oxygen vacancy concentrations of LE and HE\nmaterials directly from simulations. We believe Section 2F, Figure 4, Supplemental\nSections S4-6, and scattered edits in other sections are substantial contributions to this\nmanuscript making it suitable for publication in Nature Communications.\nReviewer #4\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports.\nThis is part of the Nature Communications initiative to facilitate training in peer review and\nto provide appropriate recognition for Early Career Researchers who co-review\nmanuscripts.\nThis reviewer does not provide additional feedback to address."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I am glad that the authors have made great efforts to compare experiments and modeling in the revised manuscript. These comparisons revealed large discrepancies between theory and experiments, despite some agreements when the vacancy concentration is low and the temperature is below 1000C. 1. The oxygen vacancy concentration variation with Temperature. In experiments, Fig. 1a and Fig 4d show either a linear increase or an accelerated (slope increases with T) increase. These two trends are replotted schematically in Fig 2f based on experimental data. However, the modeling results showed different trends. The increase with T slows down at higher temperatures. The schematic drawing in Fig. 4g has a totally different trend from experiments and shows plateaued oxygen vacancy concentrations. 2. The range in which this model is validated must be discussed. 3. There are some overclaims about the modeling. The statistics model (Eq. 9, 10, 11) is generally similar to the model in Park et al. (Ref. 13, npj computational materials, 2023). The authors did not discuss what is new in the model and what has already been developed before. 4. Park et al. considered interacting vacancies and found that it is extremely important with increasing temperature and increasing oxygen vacancies. The authors never commented on whether they have considered interacting vacancies, why or why not consider the interactions. 5. This modeling difference may contribute to the validity of the model, especially the vacancy concentration at higher temperatures. In fact, the oxygen vacancy concentration increases with temperature, and the slope also increases with temperature in Park’s paper, consistent with the experimental trends and schematics discussed in comment 1. Overall, I don’t think a model has to capture experiments completely. However, not stating its limitations and only focusing on the agreements with experiments may lead to biased predictions.", "reviewer_2": "The authors have addressed the reviewer comments well. The revised paper is consequently stronger than the original paper. I now recommend that the paper be accepted for publication.", "reviewer_3": "Concerns regarding the revised manuscript, for which I request further correction or clarification, are as follows: 1. The authors estimate the oxygen deficiency using TGA. However, a holding time of only 1–2 hours during the TGA measurements is insufficient for the system to reach equilibrium. Therefore, comparison with calculations that assume equilibrium is not valid. The oxygen nonstoichiometry should be determined under equilibrium conditions in order to allow meaningful discussion. 2. In Eq. 2, charge compensation for A-site 2+ doping is restricted to B-site holes and oxygen vacancies, without distinguishing between different B-site cations. This treatment does not reflect the electronic states of Co and Fe in the Co0.2Fe0.8 composition. The effect of the Co/Fe ratio on B-site hole formation should be examined. For example, if the composition were Co0.5Fe0.5, the results could differ significantly. It is therefore important to clarify the valid range and limitations of the present study. 3. Although the manuscript describes a transition from the LE/HE notation to Low σA/ High σA, the figures, legends, and text are not yet fully consistent. A comprehensive review and unification of terminology throughout the manuscript is required. 4. The computational results discuss distortions of the B–O bonds and distributions of site-specific bonding energies. Experimental validation of local structural changes around the B-site is necessary. Techniques such as EXAFS, TEM, PDF, or Raman spectroscopy are available, and the manuscript would be strengthened by incorporating or at least discussing such experimental observations.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "RE: Nature Communications 2nd Revision\nOriginal Manuscript ID: NCOMMS-25-34369\nTitle: \"A statistical understanding of oxygen vacancies in distorted high-entropy\nperovskite oxides\"\nAuthors: Adam Potter; Yifan Wang; Kiran Hamkins; Dongjae Kong; Yuzhe Li; Jian Qin;\nXiaolin Zheng\nResponse to Reviewers’ Comments: We thank all reviewers for their thorough and\ncritical reviews. We have addressed the comments below.\nWe noted three primary suggestions from the reviewers and responded with the\nfollowing contributions:\nSuggestion 1: Model limitations and assumptions\nThe authors should clearly state the model limitations, validation regime, and\ndifferences from previous works, particularly Park et al. (Comments 1.1-5, 3.2)\n• Response: We agree with these comments and have revised the text to clearly\nidentify 1) limiting model assumptions, 2) the compositional and environmental\nrange used to validate the model, 3) and novel contributions to modeling and\nanalysis. We believe this better qualifies which of our results are material-specific\nand which could have wider implications.\nSuggestion 2: Experimental evidence of bond distortion\nDiscussion around lattice distortion and B-O bond variations would be enriched by\nincorporating direct experimental evidence. (Comment 3.4)\n• Response: Multiple studies have observed bond distortions due to high-entropy\ncation mixing using STEM. We have added these references to our discussion.\nSuggestion 3: Ensuring TGA experiments are in equilibrium\nThe authors should prove the oxygen vacancy formation reaction is in equilibrium to\nsupport their thermodynamic analysis. (Comment 3.1)\n• Response: Experimental TGA data suggests rapid weight stabilization in\nresponse to temperature steps suggesting oxygen exchange has reached\nequilibrium. Estimates for the oxygen diffusion rates in LSCF also suggest a\nrelaxation time on the scale of seconds.\nReviewer comments are in black, our response to reviewer comments will be in blue and\nrevised portions of the manuscript in red here and in the main text and supplemental.\nI am glad that the authors have made great efforts to compare experiments and\nmodeling in the revised manuscript. These comparisons revealed large discrepancies\nbetween theory and experiments, despite some agreements when the vacancy\nconcentration is low and the temperature is below 1000C.\nComment 1: The oxygen vacancy concentration variation with Temperature. In\nexperiments, Fig. 1a and Fig 4d show either a linear increase or an accelerated (slope\nincreases with T) increase. These two trends are replotted schematically in Fig 2f based\non experimental data. However, the modeling results showed different trends. The\nincrease with T slows down at higher temperatures. The schematic drawing in Fig. 4g\nhas a totally different trend from experiments and shows plateaued oxygen vacancy\nconcentrations.\nThe reviewer comments that while experimental data suggests linear or exponential\ngrowth in δ, the model predictions in Figure 4e and the schematic in Figure 4g suggest\nsaturation at higher temperatures. This saturation is a physical predication of our model\nand is consistent with other works. Although our experiments do not reach the full\nsaturation plateau, the onset of this behavior is observable. We have revised the text to\nclarify those points.\nThe model predicts saturation because the primary mechanism for oxygen vacancy\nformation (O + 2h → V + ½ O ) relies on the reduction of polaron holes (i.e.,\nx O 2\nFe4+→Fe3+, less favorable Fe3+→Fe2+). The concentration of the polaron holes is finite\nand determined by 2+ dopant concentrations. Once the holes are consumed, the\nvacancy concentration plateaus. This saturation effect is well-documented in literature.\nIn similar 2+ doped perovskite oxides, oxygen vacancy formation has been observed to\nstagnate at δ = ½ [X ] unless the pO2 is very low (< 10-10 atm) or possibly if\n2+\ntemperature is high (>1200C) (Bae et al., 2019; Søgaard et al., 2007). This was also\nobserved for (La Sr )(Cr Ti )O for y=0.1/0.2/0.3, (La Sr )(Cr Mn )O , and\n0.7 0.3 1-y y 3 0.75 0.25 0.5 0.5 3\n(La Sr )(Cr Fe )O (Oishi et al., 2008). Our model explicitly incorporates this\n0.75 0.25 0.5 0.5 3\nphysical limit. The term “(X -2δ)2” in the denominator (Eq 4) represents the diminishing\n2+\nconcentration of holes, resulting in saturation as of δ in Figure 4e and Figure 4g. This\ntreatment has also been used in models elsewhere (Mizusaki et al., 1985; Oishi et al.,\n2008) and has successfully predicted the effects of Sr concentration (2+ doping) on\noxygen vacancy (Bae et al., 2019).\nAlthough our experiments do not reach the full saturation plateau, the onset of this\nbehavior is observable. In Figure 4d, the oxygen vacancy formation curves of HE-5 and\nHE-3 show a slight but distinct deceleration (downward curvature) in vacancy formation\nabove 800oC. This marks the beginning of saturation. The model prediction in Figure 4e\nsuccessfully captures this deceleration for HE-5. Meanwhile, our model correctly\n2\npredicts linear growth in HE-2 and HE-4 through 1000oC. Figure 4g extends the model\npredication to temperatures beyond our experimental range (gray region), resulting in\nbehavior consistent with the literature, but unvalidated within this work.\nWe have added related discussion on vacancy saturation and associated assumptions\nin Section IIB and Section IIF.\n“ 𝑋 is critical because previous studies have shown oxygen vacancy concentration\n2+\n1\nsaturates at δ = 𝑋 due to polaron depletion until the system reaches highly reducing\n2+\n2\nconditions (pO ≤ 10−10atm for similar materials)”\n2\n1\n“…due to polaron hole depletion when δ → 𝑋 (see Section IIB) …”\n2+\nWe agree that we should clearly state the validation range of our model. The model is\nvalidated for (1) high oxygen partial pressure (pO =0.21 atm), (2) constant B-site\n2\n(Co Fe ) and (3) low to moderate oxygen vacancy (δ<0.15). Our model assumes that\n0.2 0.8\ninteractions between oxygen vacancies and extrinsic charge (polaron holes) are due to\ncation doping; it is only relevant for similar doped systems (X > 0). For example, the\n2+\nmodel could fail under reducing environments (including deep vacuum) or higher\ntemperatures when the system enters a region on the Brouwer diagram where oxygen\nvacancies no longer interact with polaron holes (Mizusaki et al., 1985). Previous work\nsuggests this happens to La Sr FeO , a similar material to LE-1, when pO <10-10 atm\n1-x x 3 2\n(Bae et al., 2019; Søgaard et al., 2007).\nWe have added those validation ranges and limitations of our model in Section IIF and\nthe Conclusion.\n“Notably, the model in Eq. 4 is only valid for oxides with extrinsic polaron holes that\ninteract with oxygen vacancies, as is widely accepted for LSCF-like perovskite oxides\n\\cite{jiang_development_2019}. This occurs when the divalent doping exceeds oxygen\n1\nvacancy concentrations (δ ≤ 𝑋 ) which is typical in the oxidizing conditions ( pO =\n2+ 2\n2\n0.21 atm) used in this study. Furthermore, this model does not differentiate polarons by\nB-site, in part because polarons prefer forming on iron over cobalt (Jiang,\n2019)(Hartmann, 2023) so long as iron sites are available. For lower iron versus cobalt\ncompositions where [Fe] < 𝑋 polarons could overflow onto cobalt sites affecting\n2+\noxygen vacancy dynamics.”\n3\n“The above model has only been validated for perovskite oxides in air with Fe Co B-\n0.8 0.2\nsites and divalent doping of 0.17 ≤ 𝑋 ≤ 0.5, however, the theoretical framework\n2+\nbehind energetically heterogeneous oxygen sites could be relevant to other complex\noxides.”\nComment 3: There are some overclaims about the modeling. The statistics model (Eq.\n9, 10, 11) is generally similar to the model in Park et al. (Ref. 13, npj computational\nmaterials, 2023). The authors did not discuss what is new in the model and what has\nalready been developed before.\nThis work is an extension of the excellent research of Park et al., as mentioned in\nparagraph 4 of the Introduction, who brought attention to oxygen vacancy energy\nvariance in distorted perovskites. As an inspiration to this work, it is important to\ndifferentiate contributions of our work beyond what has been developed by Park et al.\nBelow we list our work’s general contributions, which will be clarified in the Introduction\nand Conclusion, and specific differences to the model, which will be added to Section\nIIF. We believe much of our contributions lie mostly in the first category.\nBroad Contributions\n• Park et al. (Park et al., 2023a) suggest a statistical model (labeled as “HE”) in\ncontrast to a traditional model (“LE”) because high variance in vacancy energies\nfrom distortion could lead to measurable differences in behavior. This is a notable\ninnovation, however, experimental comparisons in Figure 5 do not clearly\nmotivate one model over the other for differentiating materials. As such, Park et\nal. introduced vacancy energy variance as a potentially important parameter but\ndoes not prove its necessity for predicting differences between materials. Our\nstudy identifies experimental trends in the thermodynamics of 14 materials which\ncan be explained theoretically by vacancy energy variance (Figures 2c-e, Figure\n4a-b). We also find our statistical model better matches experimental Δδ values\ncompared to the traditional model for predicting differences among 8 materials\n(Figurers 4d-f). This provides strong evidence that including vacancy energy\nvariance is necessary for accurate predictions.\n• We built on the numerical analysis of Park et al. (Figure 6 in (Park et al., 2023b))\nby analytically deriving the effects of vacancy energy variance on formation\nthermodynamics from statistical mechanics. This provides physical intuition into\nhow vacancy energy statistics affect vacancy formation enthalpy and entropy\n(Eqs. 13-16). These derived expressions correctly predict a decrease in\nformation enthalpy and entropy with increasing vacancy energy variance (Figures\n4a-b) and were even used to estimate vacancy energy variance experimentally\n(Figure 4c).\n• Finally, by measuring 14 samples over a range of A-site compositions we can\ndefine composition-property relations. Park et al. suggested vacancy energy\nvariance increases with the number of A-site cations among their 4 materials, but\nthat the “E distribution does not show any strong correlation with the local A-site\nv\ncharge distribution or the structural distortions.” In contrast, we found A-site size\nvariance to be a more relevant predictor suggesting it matters which cations are\nmixed, not just how many. This results in clear guiding principles for tuning\noxygen vacancies using X and A-site size variance, easily calculated from\n2+\ncomposition.\nOxygen Vacancy Model Differences\n• One difference between our model and that proposed by Park et al. is the use of\ndefect models to describe 2+ cation concentration and polaron holes. This\ndistinction relates to the term (X -2δ)2 in our model and is discussed further in\n2+\nresponses to Comments 1, 2, and 4.\n• Additionally, we include contributions of oxygen vibrational entropy which are\ntypically dismissed as negligible. Our calculations, detailed in the Supplemental,\nfound the loss of lattice vibrational entropy from an oxygen removal is ~25 J/mol-\nK which is significant (scaling the equilibrium constant by e-25/8.314). Differences in\noxygen vibrational entropy between materials (see Figure S5) are cited as a\npossible source of model error in Figure 4b.\n• The last difference is we use a universal interatomic potential to sample values of\nE instead of DFT. This is a methodology choice resulting in superior oxygen site\nv\nsampling (N=768 vs N=96) at the cost of reduced E accuracy relative to DFT.\nv\nWhile the cost to accuracy is not clear, it is possible a predictive framework using\nuniversal interatomic potentials could also be more accessible to researchers\nthan DFT.\nWe have added explicit comparisons with Park et al. in Section IIF and better\nhighlighted specific contributions listed above in the Introduction and Conclusion.\n“However, the importance of such statistical effects for predictive models is unclear. The\nmechanisms of 𝐸 statistics has also remained largely unexplored for a wide range of\nv\nhigh-entropy oxides compositions, both experimentally and computationally.”\n“This model expands one the model suggested by Park et al. by coupling to vacancy\nformation to polaron destruction, considering oxygen vibrational entropy ( 𝑆 ≠ 0 ),\nO,vib\nand increasing site sampling with a machine-learned interatomic potential.”\n“Since Park (Park, 2023) and Xu (Xu, 2024) introduced using 𝑔(𝐸̂,σ2 ) instead 𝐸̂ to\nv Ev v\npredict oxygen vacancies, this work is the first to use broad experimental data to\nsuggest the superiority of the statistical model. Our results confirm that modeling\n5\nvariations between oxygen sites is crucial for predicting the unique behavior of oxygen\nvacancies in our high-entropy perovskite oxides. Other contributions of this work include\nthe theoretical analysis in Section IIE providing physical intuition and guiding principles\nfor how vacancy energy distributions affect bulk formation thermodynamics. Finally,\ncomparisons of 14 A-site compositions identified σ parameter as a key composition-\n𝐴\nproperty relation for guiding future material design. ”\nComment 4: Park et al. considered interacting vacancies and found that it is extremely\nimportant with increasing temperature and increasing oxygen vacancies. The authors\nnever commented on whether they have considered interacting vacancies, why or why\nnot consider the interactions. This modeling difference may contribute to the validity of\nthe model, especially the vacancy concentration at higher temperatures. In fact, the\noxygen vacancy concentration increases with temperature, and the slope also increases\nwith temperature in Park’s paper, consistent with the experimental trends and\nschematics discussed in comment 1.\nWe agree with the reviewer that vacancy-vacancy interactions can play a role, as noted\nin Section IIF, paragraph 1. However, the referenced paper applies Boltzmann statistics\nto each site independently, whereas we adopted the Fermi-Dirac statistics. The latter\naccounts for the equilibrium among different sites by enforcing a common chemical\npotential, which is more appropriate for finite concentrations but complicates the direct\ncomparison. It is also unclear to us how a linear shift in vacancy energy captures\nvacancy-coupling effects, since a formal first-order treatment would require inclusion of\nthe quadratic term in the excess free energy - the same order as the contribution our\nanalysis attributes to vacancy-energy dispersion (𝜎 ). This implies that the strength of\nEv\nvacancy-vacancy interactions cannot be cleanly separated from the effects of vacancy-\nenergy dispersion. In our model, the fitted dispersion parameter matches the calculated\ndispersion (Figure 4C), suggesting that vacancy-vacancy interactions are minor or, at\nmost, comparable in magnitude to the dispersion contribution. While Park et al. (2023)\nalso emphasized the importance of the distribution of vacancy formation energies for\nthermodynamics, their discussion remained qualitative and assumed independent\nvacancies following a Boltzmann distribution, an assumption valid only in the classical\nlimit of 𝑇 → ∞ or 𝑛 → 0. In contrast, our work starts from a rigorous derivation based on\nchemical potentials and Fermi-Dirac statistics. A central result is that the variance of the\nenergy distribution (𝜎2 ) introduces well-defined quadratic terms in both the vacancy-\n𝐸𝑣\nformation enthalpy and entropy (Eqs. 13–14), terms not accounted for in Park et al.\nThese quadratic corrections provide key physical insight into non-ideal mixing in high-\nentropy perovskites, analogous to the quadratic activity corrections in regular-solution\ntheory and reveal how mixing enthalpy and entropy contribute to oxygen-vacancy\nformation in high-entropy oxides.\nmay account for some of what is referred to as “vacancy interactions”. Park et al. argue\nthat when vacancies exceed the dilute limit (δ>0.05) their large polarons repel each\nother, slowing further vacancy formation. This assumes each oxygen vacancy creates\nlarge (less localized) negative polarons due to charge exchange. In contrast, we\nassume our system has divalent dopants (this is true for SLFO, BLFO, NSBLFO, but not\nNLFO in Park et al.) which create positive polaron holes localized on B-sites. Then,\nforming an oxygen vacancy don’t create negative polarons, but destroys positive\npolarons; this is a common treatment for the defect modeling in similar perovskites (Bae\net al., 2019; Choi et al., 2014; Oishi et al., 2008). When polaron holes are treated as a\ndepleting reactant in the law of mass action (via the term (X -2δ)-2 in Eq. 4), the model\n2+\nslows oxygen vacancy formation with increasing concentration. This has a similar effect\non limiting vacancy formation to the interaction term used in Park et al. Furthermore, our\ncurrent model’s predictions closely match experiments in magnitude unlike the “HE-DIL”\nmodel in Park et al (Supplementary Figure S5 in their paper). Overall, the argument is\nsimilar: in Park et al. oxygen vacancies are accompanied by a relative “excess” charge\nrepelling additional vacancies, in this defect model oxygen vacancies are accompanied\nby a “depletion” in charge making it harder to form more vacancies. Our polaron\nmodeling does not account for all interaction mechanisms, but it may capture the\noxygen vacancy coupling with polaron coulombic repulsion described in the Introduction\nof Park et al.\nComment 5: Overall, I don’t think a model has to capture experiments completely.\nHowever, not stating its limitations and only focusing on the agreements with\nexperiments may lead to biased predictions.\nWe do not want to overclaim our model capabilities and have made various edits to the\nConclusion and final paragraph of the Introduction to better define the studies scope.\nAs per the comments above, we have revisited the text to minimize model\noverstatements.\n“…in LSCF-based perovskite oxides… “\n“…LSCF-based perovskite oxides in oxidizing conditions…”\nThe authors have addressed the reviewer comments well. The revised paper is\nconsequently stronger than the original paper. I now recommend that the paper be\naccepted for publication.\nWe thank this reviewer for their positive feedback!\nConcerns regarding the revised manuscript, for which I request further correction or\nclarification, are as follows:\nComment 1: The authors estimate the oxygen deficiency using TGA. However, a holding\ntime of only 1–2 hours during the TGA measurements is insufficient for the system to reach\nequilibrium. Therefore, comparison with calculations that assume equilibrium is not valid.\nThe oxygen non-stoichiometry should be determined under equilibrium conditions in order\nto allow meaningful discussion.\nIn this work, the equilibrium assumption is critical as it justifies the study of\nthermodynamics instead of kinetics. The reviewer is requesting further justification that\nthe oxygen exchange with the environment is in equilibrium. Our first evidence comes from\nthe rapid changes, then steady holds, in sample weight in TGA. In the supplemental\ninformation, Figure S2 shows the LE-1 (LSCF) mass response to temperature steps after\nsubtracting the empty crucible response. Shown below is the mass response for LE-1 (left)\nand (HE-2) right as temperature is increased stepwise. Even below 600C, temperature\nsteps cause the sample weight to rapidly drop, indicating vacancy formation, followed by a\nsteady response suggesting oxygen exchange reached equilibrium.\nThis rapid response is expected for oxygen conductors, especially in powders where the\ncharacteristic length for oxygen diffusion is on the micron scale. For example, the oxygen\ndiffusion coefficient in LSCF at 600C is measured to be approximately D=10-11 m2/s\n(Endler-Schuck et al., 2015). Our SEM analysis suggests our powders are as large as L=10\num. This suggests the characteristic relaxation time for the diffusion process is on the\nscale of τ=L2/D, or 10 seconds, in agreement with TGA. Altogether, we feel that the\nequilibrium assumption is valid.\nWe have added the figure above for HE-2 to Figure S2 in the supplemental.\nWe cannot make explicit claims about materials with varying B-site compositions as we\ndid not test this parameter experimentally and it is outside the scope of the study.\nHowever, we can provide justification for why our model does not distinguish Fe from Co in\nCo Fe and provide an argument for the range of Co-Fe ratios we might expect our model\n0.2 0.8\nto be relevant.\nIn our Co Fe system, oxygen vacancy formation is coupled to the Co/Fe electronic state\n0.2 0.8\nthrough polaron destruction: O x + 2h → V O + ½ O 2. This is equivalent to reducing\nCo4+/Fe4+ as oxygen is removed. If the energy to destroy an Fe-localized polaron is different\nthan a Co-localized polaron it could affect the oxygen vacancy formation thermodynamics.\nHowever, we believe this effect is minimal because the system’s permanent holes prefer\nto form on Fe. A review of experimental LSCF studies suggests that polaron holes tend to\nlocalize on Fe, while Co contributes additional charge through disproportionation (Jiang,\n2019). A DFT study also found polarons prefer localizing on Fe (Hartmann et al., 2023).\nWhile this cannot guarantee no Co localized holes, it suggests holes are likely to localize\non Fe so long as Fe sites are available. Furthermore, high electrical conductivity (>100\nS/cm for LE-1 (Jiang, 2019)) suggests if polarons did form on higher energy Co sites,\nthey could readily migrate to lower energy Fe sites.\nHowever, this understanding suggests a possible compositional threshold beyond which\nour model may be less accurate. When perovskites with Co/Fe B-sites satisfy [A2+]>[Fe]\nthere could be more polaron holes than Fe sites and polaron holes may be forced to\n“overflow” onto higher energy Co sites. This could create two regions of oxygen vacancy\nbehavior: the first where vacancy formation reduces Co, and another where they reduce\nFe. None of our materials satisfy [A2+]>[Fe], but it could be relevant to note in the text that\nour model may fail when extended to this compositional region. We will clarify this\ncondition in Section IIB.\nThis potential limitation on generalizability to B-site composition has been added to\nSection 2B.\n“Notably, the model in Eq. 4 is only valid for oxides with extrinsic polaron holes that\ninteract with oxygen vacancies, as is widely accepted for LSCF-like perovskite oxides\n\\cite{jiang_development_2019}. This occurs when the divalent doping exceeds oxygen\n10\n1\nvacancy concentrations (𝛿 ≤ 𝑋 ) which is typical in the oxidizing conditions ( pO =\n2+ 2\n2\n0.21 atm) used in this study. Furthermore, this model does not differentiate polarons by\nB-site, in part because polarons prefer forming on iron over cobalt (Jiang, 2019)\n(Hartmann, 2023) so long as iron sites are available. For lower iron versus cobalt\ncompositions where [Fe] < 𝑋 polarons could overflow onto cobalt sites affecting\n2+\noxygen vacancy dynamics.”\nComment 3: Although the manuscript describes a transition from the LE/HE notation to\nLow σ / High σ , the figures, legends, and text are not yet fully consistent. A comprehensive\nA A\nreview and unification of terminology throughout the manuscript is required.\nIn the previous round of revisions, we worked toward distinguishing the “LE/HE” and “Low\nσ / High σ “notations in the figures and text. The difference can be difficult to\nA A\ncommunicate; while samples were synthesized to span LE and HE compositional\ncomplexities, the study later concludes oxygen vacancy behavior is better distinguished by\nA-site size variance (Low σ / High σ ). Our ultimate takeaway is based on the A-site\nA A\nvariance, but LE and HE roughly correlate with A site variance making it easier for readers\nidentify samples.\nAs such, we feel it is best to label samples with LE and HE, while also using LE vs HE to\ndescribe initially observed trends in Figure 1. As soon as we introduce σ and identify it as a\nA\nkey descriptor for oxygen vacancy trends in Figure 2, all further figures and discussion will\nbe in terms of σ .\nA\nThe figures and text terminology have been revised as per the reasoning above.\nComment 4: The computational results discuss distortions of the B–O bonds and\ndistributions of site-specific bonding energies. Experimental validation of local structural\nchanges around the B-site is necessary. Techniques such as EXAFS, TEM, PDF, or Raman\nspectroscopy are available, and the manuscript would be strengthened by incorporating or\nat least discussing such experimental observations.\nOur model attributes the variation in oxygen vacancy formation energy to the local\nstructural and chemical heterogeneity around the B-site. This assumption is a\nfundamental concept in the study of high-entropy materials and is strongly supported by\ndirect experimental observations in the literature. For example, Su, et al used STEM to\ndirectly observe structural fluctuations in high-entropy perovskite oxides, including\noctahedral tilt, as a consequence of cation mixing. They note, “the existence of\ncompositional fluctuations along with significant composition-dependent distortion of\n11\nmetal-oxygen bonds is observed” (Su et al., 2022). Recently, Zou, et al also used STEM to\ncharacterize variations in cobalt octahedral tilt and Co-O bond lengths in perovskite oxides\nwith high-entropy A-sites (Zou et al., 2025). This discussion will be added to Section IID to\nstrengthen claims of B-O distortion.\nThese experimental studies and a brief discussion have been added to discussion Section\nIID1.\n“Such octahedra distortions have been directly observed with scanning transmission\nelectron microscopy (STEM) in high-entropy oxides (Su, 2022) (Zou, 2025).”\nBae, H., Hong, J., Singh, B., Srivastava, A. K., Joo, J. H., & Song, S.-J. (2019). Investigations\non Defect Equilibrium, Thermodynamic Quantities, and Transport Properties of La\n0.5\nSr FeO . Journal of The Electrochemical Society, 166(4), F180–F189.\n0.5 3-δ\nhttps://doi.org/10.1149/2.0311904jes\nChoi, M.-B., Jeon, S.-Y., Singh, B., Yoo, Y.-S., Hwang, J.-H., & Song, S.-J. (2014). Correlation\nbetween defect structure and electrochemical properties of mixed conducting\nLa0.1Sr0.9Co0.8Fe0.2O3−. Acta Materialia, 65, 373–382.\nhttps://doi.org/10.1016/j.actamat.2013.11.006\nEndler-Schuck, C., Joos, J., Niedrig, C., Weber, A., & Ivers-Tiffée, E. (2015). The chemical\noxygen surface exchange and bulk diffusion coefficient determined by impedance\nspectroscopy of porous La0.58Sr0.4Co0.2Fe0.8O3 − δ (LSCF) cathodes. Solid State\nIonics, 269, 67–79. https://doi.org/10.1016/j.ssi.2014.11.018\nHartmann, C., Geneste, G., Saravanabavan, K., & Laurencin, J. (2023). Investigation of the\nReaction Mechanisms in La1-XSrxCo1-YFeyO3 Electrode by DFT Calculations. ECS\nTransactions, 111(6), 1239. https://doi.org/10.1149/11106.1239ecst\nJiang, S. P. (2019). Development of lanthanum strontium cobalt ferrite perovskite\nelectrodes of solid oxide fuel cells – A review. International Journal of Hydrogen\nEnergy, 44(14), 7448–7493. https://doi.org/10.1016/j.ijhydene.2019.01.212\nMizusaki, J., Yoshihiro, M., Yamauchi, S., & Fueki, K. (1985). Nonstoichiometry and defect\nstructure of the perovskite-type oxides La1−xSrxFeO3−°. Journal of Solid State\nChemistry, 58(2), 257–266. https://doi.org/10.1016/0022-4596(85)90243-9\nOishi, M., Yashiro, K., Sato, K., Mizusaki, J., & Kawada, T. (2008). Oxygen nonstoichiometry\nand defect structure analysis of B-site mixed perovskite-type oxide (La, Sr)(Cr,\nM)O3−δ (M=Ti, Mn and Fe). Journal of Solid State Chemistry, 181(11), 3177–3184.\nhttps://doi.org/10.1016/j.jssc.2008.08.015\nPark, J., Xu, B., Pan, J., Zhang, D., Lany, S., Liu, X., Luo, J., & Qi, Y. (2023a). Accurate\nprediction of oxygen vacancy concentration with disordered A-site cations in high-\nentropy perovskite oxides. Npj Computational Materials, 9(1), 1–13.\nhttps://doi.org/10.1038/s41524-023-00981-1\nPark, J., Xu, B., Pan, J., Zhang, D., Lany, S., Liu, X., Luo, J., & Qi, Y. (2023b). Accurate\nprediction of oxygen vacancy concentration with disordered A-site cations in high-\nentropy perovskite oxides. Npj Computational Materials, 9(1), 29.\nhttps://doi.org/10.1038/s41524-023-00981-1\nSøgaard, M., Vang Hendriksen, P., & Mogensen, M. (2007). Oxygen nonstoichiometry and\ntransport properties of strontium substituted lanthanum ferrite. Journal of Solid\nState Chemistry, 180(4), 1489–1503. https://doi.org/10.1016/j.jssc.2007.02.012\nSu, L., Huyan, H., Sarkar, A., Gao, W., Yan, X., Addiego, C., Kruk, R., Hahn, H., & Pan, X.\n(2022). Direct observation of elemental fluctuation and oxygen octahedral\ndistortion-dependent charge distribution in high entropy oxides. Nature\nCommunications, 13(1), 2358. https://doi.org/10.1038/s41467-022-30018-y\nZou, G., Liu, H., Liu, T., Han, S., Tan, Z., Zhang, S., Guo, Y., Yu, J., Zhang, X., Lu, F., Song, Y.,\nWang, G., & Bao, X. (2025). High-entropy-induced CoO 6 octahedral distortion for\n13\nboosted oxygen evolution reaction at high temperature. Energy & Environmental\nScience, 18(21), 9478–9489. https://doi.org/10.1039/D5EE01370D\n14"}, {"round": 3, "reviewer_comments": {"reviewer_1": "The authors have addressed the reviewer's concerns.", "reviewer_3": "The results in Figure S2 indicate that the oxygen‐vacancy concentration appears to relax on a very short timescale. However, oxygen diffusion in perovskites is generally not sufficiently fast to account for such an immediate response. This raises the possibility that the observed rapid mass changes are influenced by buoyancy effects associated with temperature variations, rather than solely by oxygen vacancy formation. When the furnace temperature changes, the gas density and therefore the buoyant force acting on the sample change instantaneously, producing an apparent mass change in the TGA signal. The subsequent gradual drift would correspond to the true relaxation of oxygen nonstoichiometry toward equilibrium at each temperature. To verify whether buoyancy contributes to the rapid mass responses observed in LE-1 and HE-2, the authors should perform a control TGA experiment using an oxygen‐stoichiometry‐invariant material, such as Al₂O₃, under the same temperature‐step protocol. This would allow the buoyancy contribution to be quantified and properly separated from the intrinsic oxygen exchange behavior."}, "author_response": "RE: Nature Communications 3rd Revision\nOriginal Manuscript ID: NCOMMS-25-34369\nTitle: \"A statistical understanding of oxygen vacancies in distorted high-entropy\nperovskite oxides\"\nAuthors: Adam Potter; Yifan Wang; Kiran Hamkins; Dongjae Kong; Yuzhe Li; Jian Qin;\nXiaolin Zheng\nResponse to Reviewers’ Comments: We thank all reviewers for their thorough and\ncritical reviews. We have addressed the comments below.\nOverview\nWe have performed the experiment suggested by Reviewer 3 and found the effects of\nbuoyancy are minimal and should not affect out TGA results.\nThere are no new revisions to the manuscript main text or supplemental information.\nReviewer comments are in black, our response to reviewer comments will be in blue and\nrevised portions of the manuscript in red here and in the main text and supplemental.\nThe authors have addressed the reviewer's concerns.\nReviewer 3\nThe results in Figure S2 indicate that the oxygen‐vacancy concentration appears to\nrelax on a very short timescale. However, oxygen diffusion in perovskites is generally\nnot sufficiently fast to account for such an immediate response. This raises the\npossibility that the observed rapid mass changes are influenced by buoyancy effects\nassociated with temperature variations, rather than solely by oxygen vacancy formation.\nWhen the furnace temperature changes, the gas density and therefore the buoyant\nforce acting on the sample change instantaneously, producing an apparent mass\nchange in the TGA signal. The subsequent gradual drift would correspond to the true\nrelaxation of oxygen nonstoichiometry toward equilibrium at each temperature.\nTo verify whether buoyancy contributes to the rapid mass responses observed in LE-1\nand HE-2, the authors should perform a control TGA experiment using an oxygen‐\nstoichiometry‐invariant material, such as Al₂O₃, under the same temperature‐step\nprotocol. This would allow the buoyancy contribution to be quantified and properly\nseparated from the intrinsic oxygen exchange behavior.\nWe recognize that there are several factors that could affect crucible weight upon\nheating such as drag forces or buoyancy. We have taken steps to properly control these\nfactors to isolate oxygen loss. The effects of drag are controlled when we subtract the\nidentical TGA run using an empty crucible. To test the effects of buoyancy, we have\nperformed the experiment described by the reviewer and confirmed that the mass\nchange we observed for LE and HE is not due to buoyancy.\nTo determine if weight changes are due to oxygen loss or buoyancy effects, we\nperform TGA on 50-60mg each of LE-1, HE-2, and alumina. For each experiment, we\ncorrected the baseline by subtracting the data of an empty crucible measured using the\nsame protocol, which was shown to effectively minimize the buoyancy effect (Buoyancy\nEffect of TGA Experiment – Cement Science, 2017). For the alumina test, we\nmaintained the same temperature ramp rate to compare buoyancy effects, if any, but\nshortened the temperature hold. If buoyancy was a dominant factor, we should observe\nlarge weight changes for the alumina, which is the same sample mass and similar\ndensity (4.0 g/cm3 vs 5.2 g/cm3 for LE-1). We observe no significant weight changes,\nonly small fluctuations, from 500-1000 for alumina. These small fluctuations (0.05%) are\nnormal, but they are clearly significantly less than the apparent mass drop (0.5-1%) for\nLE-1 and HE-2. Hence, we conclude that buoyancy does not play a large role in our\nTGA data. We agree that the oxygen exchange is surprisingly fast (< 5 minutes), this is\n2\ndue to the small particle size (~300-500 nm in SEM) since the relaxation time scales\nwith L2.\nWorks Cited\nBuoyancy Effect of TGA Experiment – Cement Science. (2017, June 21).\nhttps://www.cementscience.com/2013/04/buoyancy-effect-of-tga-experiment.html\n3"}, {"round": 4, "reviewer_comments": {"reviewer_3": "In the buoyancy measurements, the measurement time at each temperature is approximately 0.2 hours, which is relatively short. Because of the short duration, the buoyancy of the gas at each temperature may not be accurately reflected. In the actual experiments, measurements are typically conducted for about two hours at each temperature. Therefore, I propose that the buoyancy measurements should likewise be performed under steady-state conditions."}, "author_response": "RE: Nature Communications 4th Revision\nOriginal Manuscript ID: NCOMMS-25-34369\nTitle: \"A statistical understanding of oxygen vacancies in distorted high-entropy\nperovskite oxides\"\nAuthors: Adam Potter; Yifan Wang; Kiran Hamkins; Dongjae Kong; Yuzhe Li; Jian Qin;\nXiaolin Zheng\nResponse to Reviewers’ Comments: We thank all reviewers for their thorough and\ncritical reviews. We have addressed the comments below.\nOverview\nWe have performed the experiment suggested by Reviewer 3 with longer holds at each\ntemperature and found the effects of buoyancy are minimal and should not affect out TGA\nresults.\nThere are no new revisions to the manuscript main text or supplemental information.\nReviewer comments are in black, our response to reviewer comments will be in blue and\nrevised portions of the manuscript in red here and in the main text and supplemental.\nIn the buoyancy measurements, the measurement time at each temperature is\napproximately 0.2 hours, which is relatively short. Because of the short duration, the\nbuoyancy of the gas at each temperature may not be accurately reflected. In the actual\nexperiments, measurements are typically conducted for about two hours at each\ntemperature. Therefore, I propose that the buoyancy measurements should likewise be\nperformed under steady-state conditions.\nWe have repeated the experiment on a longer timescale.\nMethods are identical to the previous test apart from longer holds at each\ntemperature for the alumina sample. To determine if weight changes are due to oxygen\nloss or buoyancy effects, we perform TGA on 50-60mg each of LE-1, HE-2, and\nalumina. For each experiment, we corrected the baseline by subtracting the data of an\nempty crucible measured using the same protocol. For the alumina test, we maintained\nthe same temperature ramp rate to compare any buoyancy effects. If buoyancy was a\ndominant factor, we should observe large weight changes for the alumina, which is the\nsame sample mass and similar density (4.0 g/cm3 vs 5.2 g/cm3 for LE-1).\nOver a 12-hour test with approximately 2 hour holds at each temperature, the\nalumina exhibits no significant weight changes, only small fluctuations, from 500-1000C\nfor alumina. The alumina weight fluctuated as much as 0.017% from its mean over the\ntest after an initial small drop upon heading to 500C. This is significantly less than the\napparent mass drop (0.5-1%) for LE-1 and HE-2. Hence, we conclude that buoyancy\ndoes not play a large role in our TGA data and weight changes can be attributed to loss\nof oxygen.\n2\n3"}]

Decision letter