Elucidating the proton-coupled oxygen reduction pathway in protonic ceramic fuel cells

Abstract

The proton-coupled oxygen reduction reaction (PC-ORR) at protonic ceramic fuel cell (PCFC) cathodes involves multiple charges, which are proton, oxygen ion, and electron/electron hole, and its complexity has long impeded unambiguous identification of the reaction pathway and the rate-determining step (RDS). The difficulty is amplified by the prevailing practice of subjectively positing an a priori “most probable” pathway to infer the RDS—a procedure that heightens the risk of decisive bias and error. Consequently, mutually inconsistent pathways have been proposed for ostensibly the same reaction. Here, we present a generalized microkinetic framework that infers the RDS without prior pathway construction. Applying this approach, we resolve the RDS for two widely studied PCFC cathodes, PrBa_0.5Sr_0.5Co_1.5Fe_0.5O_5+δ (PBSCF) and BaCo_0.4Fe_0.4Zr_0.1Y_0.1O_3−δ (BCFZY). PBSCF is limited by vacancy-assisted O₂ dissociation, whereas BCFZY is limited by a proton-coupled OH adsorbates formation step involving adsorbed atomic oxygen and bulk protons. While both exhibit sufficient proton transport and fast bulk diffusion such that surface reactions dominate porous cathode performance, the origin of the contrasting RDSs is traced to their different proton uptake/release mechanisms.

---

Broader context

Understanding the cathodic reaction mechanism in protonic ceramic fuel cells is pivotal to advancing their performance and commercialization. The proton-coupled oxygen reduction reaction (ORR) involves protons, oxygen ions, and holes as charge carriers, leading to a vast array of possible pathways and as many as 225 potential rate-determining steps (RDSs). This intrinsic complexity has hindered clear identification of the true RDS and the rational design of high-activity cathodes. In this work, we demonstrate that determining the RDS without bias is a prerequisite for mechanistic understanding. To achieve this, we developed a pathway-agnostic microkinetic methodology that deduces the RDS directly from experimentally measured reaction orders in oxygen and water partial pressures, avoiding the conventional assumption of a predefined “most probable” pathway. Applying this generalized framework to representative triple-conducting cathodes enabled us to unambiguously identify their respective RDSs and propose strategies for performance enhancement. The methodology provides a universal route to unbiased mechanistic insight, extendable to diverse electrochemical reactions beyond proton-coupled ORR.

1. Introduction

Solid oxide fuel cells (SOFCs), including protonic ceramic fuel cells (PCFCs), have been actively developed toward lower operating temperatures to suppress long-term degradation and improve economic feasibility using inexpensive component materials.¹ However, as the operating temperature decreases, the overall cell performance becomes increasingly limited not by the ohmic resistance of the electrolyte (related to ionic conduction) but by the polarization resistance at the electrolyte/electrode ionic interface, which governs the electrochemical reaction kinetics.² Compared to the hydrogen oxidation reaction occurring at the anode, the (proton-coupled) oxygen reduction reaction (ORR) at the cathode exhibits a significantly higher activation energy. Consequently, the cathodic ORR constitutes the dominant resistive component of the electrochemical device under low-temperature operation. This reaction proceeds through a series of complex multi-step processes.³ For example, in oxygen-ion–conducting SOFCs, the oxygen reduction reaction (eqn (1)) typically involves oxygen molecule adsorption, dissociation, ionization, and incorporation steps.

However, except for the initial adsorption step, the sequence of subsequent processes may vary. Depending on factors such as whether surface oxygen vacancies act as adsorption sites, how many electronic holes are transferred during the intermediate stages, whether the oxygen intermediates adopt superoxide (O₂⁻ or O⁻) or peroxide (O₂²⁻ or O²⁻) configurations, and whether surface ionic conduction of the cathode is involved, numerous combinations of elementary reaction steps become possible.⁴ In fact, for the reaction described in (eqn (1)), as many as 108 distinct elementary reaction pathways are theoretically possible.⁵ Consequently, identifying the definitive reaction pathway and pinpointing the rate-determining step remains an exceptionally challenging task, and extensive efforts have been devoted to elucidating these mechanisms.

In proton-conducting electrolyte-based SOFCs, also referred to as protonic ceramic fuel cells (PCFCs), the incorporation of electrolytes exhibiting superior ionic conductivity enables more facile operation at intermediate-to-low temperatures.^6,7 Owing to this advantage, PCFCs have recently attracted significant research attention. However, as mentioned above, the relatively low performance in PCFCs is often dominated by the sluggish cathodic reaction kinetics.^8,9 In particular, unlike the simple oxygen reduction reaction (ORR) occurring in conventional SOFCs, the cathode reaction in PCFCs involves an additional process in which protons—supplied from either the cathode bulk or surface—couple with oxygen species to form H₂O (eqn (2)).

Accordingly, the number of possible elementary reaction steps increases exponentially. Because it is practically impossible to account for all these possibilities, previous studies have typically relied on the researchers’ assumptions to define specific a priori reaction pathways, within which the rate-determining step (RDS) was subsequently explored.^10–13 This empirical approach has long hindered the systematic understanding of the overall reaction mechanism.

In this study, we aim to establish a theoretical framework for systematically categorizing the possible reaction pathways of the proton-coupled oxygen reduction reaction (PC-ORR) and to discuss experimental methodologies for their identification. Based on this framework, we propose a new a posteriori approach in which the RDS is first identified and the most plausible reaction pathway is subsequently constructed. As representative examples, we demonstrate this strategy for two state-of-the-art PCFC cathode materials—PrBa_0.5Sr_0.5Co_1.5Fe_0.5O_5+δ (PBSCF)¹⁴ and BaCo_0.4Fe_0.4Zr_0.1Y_0.1O_3−δ (BCFZY)¹⁵—by determining their respective RDSs and deriving the most probable reaction pathways.

In particular, a comprehensive analysis of the physicochemical properties of BCFZY has been reported by L. R. Tarutina et al., where the material was shown to exhibit mixed ionic–electronic transport with possible hydrogenation behavior, along with detailed evaluations of total conductivity and individual ionic contributions.⁹ In that study, surface reaction and bulk diffusion coefficients related to oxygen- and hydrogen-involved processes were systematically investigated, providing valuable insights into the transport and kinetic properties of BCFZY. However, the reported values span a wide range, making it difficult to draw a definitive conclusion regarding the governing reaction-limiting process. This indicates that further in-depth and systematic investigations are required to fully clarify the reaction kinetics of BCFZY.

2. Results and discussion

2.1. Reaction pathway of proton-coupled ORR

The possible reaction pathways at the PCFC cathode can be categorized into five representative cases, as illustrated in Fig. 1. Broadly, they can be divided according to whether the cathode material exhibits sufficient proton conductivity. In the first case, when the cathode does not possess significant proton-conducting properties—i.e., protons cannot migrate through the cathode bulk—the electrochemical reaction proceeds near the triple-phase boundary (TPB), where adsorbed oxygen species or adsorbed protons on the cathode surface diffuse and react. This mechanism is analogous to that of a conventional SOFC employing a cathode such as La_1−xSr_xMnO₃ (LSM), which lacks oxygen-ion conductivity.

Fig. 1 Five representative reaction pathways for proton-coupled oxygen reduction reaction (PC-ORR) in PCFC cathodes are governed by the proton conductivity of the cathode material. In the TPB/Near-TPB limited pathway(insufficient σ_proton), the dominant route depends on the relative diffusivity of protons and oxygen species at the cathode surface. In the extended surface pathway (sufficient σ_proton), proton release occurs via either a dehydrogenation mechanism (not involving oxygen vacancies) or a dehydration mechanism (involving oxygen vacancies), the latter of which is further classified based on the relative rates of bulk oxygen diffusion (D_chem,O) and surface oxygen exchange (k_chem,O).

In contrast, when the cathode material exhibits sufficient proton-conducting properties, protons supplied through the electrolyte can migrate via bulk diffusion and participate in the proton-coupled ORR at the cathode surface. Such materials are often referred to as triple-conducting oxides (TCOs) or triple ionic–electronic conductors (TIECs).¹⁶ In this case, the effective reaction region can extend from the TPB to a broader portion of the cathode surface, which is advantageous for reaction kinetics. (Of course, the actual electrochemically active area depends on the interplay between the cathode geometry, proton conduction path length, and surface reaction kinetics.).

Despite this potential advantage, the mechanistic details—particularly the proton release process within the cathode and the relative rate limitation between bulk diffusion and surface reaction—remain insufficiently understood, hindering a precise elucidation of the overall reaction pathway.

Specifically, whether the proton release mechanism in the cathode material follows a dehydration reaction or a dehydrogenation reaction has been largely overlooked, and most prior studies have proceeded under the implicit assumption that the reaction follows a hydration/dehydration process. If the cathode operates via a dehydrogenation mechanism, oxygen vacancies within the electrode are not directly involved in the reaction, and the process corresponds precisely to the proton-coupled oxygen reduction reaction (PC-ORR) described in (eqn (2)).¹⁷ However, if the cathode follows a dehydration mechanism, then—considering the continuous operation of the fuel cell—the cathode must undergo not only the dehydration reaction itself but also a separate oxygen reduction reaction sequentially, as shown below.

These two types of proton-release mechanisms have seldom been explicitly distinguished in explaining PCFC cathode reactions, and in many cases, they have even been used interchangeably. For instance, in the dehydration mechanism, two oxygen vacancies are involved in the reaction, whereas in the dehydrogenation mechanism, oxygen vacancies are not required. Some studies, however, have proposed an alternative interpretation in which only a single surface oxygen vacancy participates, serving as the catalytic site for O₂ molecule adsorption and subsequent O₂ dissociation.¹³ In such a case, the first adsorbed atomic oxygen may be released via the dehydrogenation pathway, while the second adsorbed oxygen follows the dehydration pathway. These mixed mechanistic interpretations highlight the necessity of re-examining the proton-release process based on the intrinsic thermodynamic properties of the cathode material.

In addition, when the proton-release process follows the dehydration mechanism, and the bulk oxygen chemical diffusion is significantly faster than the surface oxygen reduction kinetics, an additional reaction pathway can be considered. In this case, the surface oxygen vacancies generated by dehydration are not replenished by oxygen from the gas phase but rather by oxygen species diffusing from the bulk. This alternative pathway underscores the strong interplay between bulk oxygen transport and surface reaction kinetics in determining the overall cathodic mechanism of PCFCs.

2.2. Generalized microkinetic model

Even when the proton-release mechanism and the relative rates of bulk and/or surface reactions are compared to identify the slower process, thereby narrowing the possible reaction pathways in Fig. 1 to one of the five representative cases, the large number of sequential elementary steps still makes it challenging to isolate the true rate-determining step (RDS). In particular, the direct identification of high-temperature intermediate species is experimentally difficult, prompting continuous efforts to perform mechanistic analysis based on macroscopic experimental data.³

One common strategy involves refining microkinetic models using experimentally measured reaction orders with respect to the oxygen partial pressure (pO₂) and the water partial pressure (pH₂O).¹³ Numerous studies have applied such models to analyze the ORR kinetics of SOFC and PCFC cathodes.^{10–13,18,19} However, when the assumed sequential reaction pathway does not accurately represent the actual mechanism, key intermediate steps corresponding to the true RDS may be omitted. As a result, different steps can appear to yield identical reaction orders, leading to potentially misleading mechanistic interpretations.

To address this issue, Chueh and co-workers proposed a framework that classifies the individual elementary steps into preceding, rate-determining, and following categories, and expresses the reaction rate of the RDS in terms of measurable or definable activities—namely, those of pO₂, holes, oxygen vacancies, and protons—based on the number of participating ionic and electronic point defects whose concentrations can be experimentally quantified.^3,5 They demonstrated that comprehensive information on the full reaction sequence is not necessary. However, in fact, the information concerning the steps preceding or following the RDS is not essential for determining the reaction rate. Only the number of point defects that are incorporated into—or directly participate in—the reactant species during the RDS is required. Building upon this approach, we generalize the formulation to the proton-coupled oxygen reduction reaction (PC-ORR) that governs PCFC cathode behavior.

In this framework, the tunable and/or measurable parameters are pO₂, oxygen vacancies , holes (h˙), and protons . The exchange rate of the RDS can thus be expressed as a function of the activities or concentrations of these species. Although electrons (e′) may also participate in the ORR under certain conditions,⁴ most PCFC cathodes are hole-dominant conductors; therefore, electrons are neglected for simplicity. Similarly, while pH₂O is also a tunable parameter, its dependence is implicitly incorporated in the point-defect concentrations of the equilibrium steps following the RDS. Consequently, pH₂O can be omitted as an independent variable without loss of generality. The only exceptional case is when the final product of the proton-coupled ORR, H₂O(g), becomes transport-limited due to insufficient gas-transport pathways in the porous electrode. In this situation, H₂O(g) must effectively be treated as a reactant species, and thus pH₂O must be explicitly considered. Under these considerations, the global reaction rate of the RDS can be expressed as shown in (eqn (4)).

Based on these considerations, the possible exponents for each parameter can be enumerated arithmetically and are summarized in Tables 1–4. In summary,

Table 1 The dependence of reaction rate of RDS on pO₂ for various situations

No.	Surface-adsorbed oxygen species as reactants in RDS	Exponent in pO2	Example
1	Diatomic	1	O^n−2,ad, H2O2,ad
2	Atomic	1/2	O^n−ad, H2Oad
3	None (lattice incorporated)	0

---

Table 2 The dependence of reaction rate of RDS on

No.	Number of as reactant in RDS	Exponent in	Example
1	0	0
2	1	1
3	2	2

---

Table 3 The dependence of reaction rate of RDS on [h˙]

No.	Number of h˙ incorporated in reactant in RDS	Exponent in [h˙]	Example
1	0	0
2	1	−1	O^−2,ad, O^−ad
3	2	−2	O^2−2,ad, H2Oad
4	3	−3	H2O^−2,ad
5	4	n/a (impossible)	–

---

Table 4 The dependence of reaction rate of RDS on

No.	Number of as reactant in RDS	Exponent in	Example
1	0	0	–
2	1	1
3	2	2
4	3	n/a (impossible)	—
5	4	n/a (impossible)	—

---

(1) pO₂: the exponent α (in eqn (4)) is determined by whether the surface-adsorbed oxygen species serving as the reactant in the RDS exist in a diatomic or atomic form. Specifically, α = 1, ½, or 0 corresponds to diatomic O₂ adsorption, atomic adsorption (atomic O), or the lattice incorporated oxygen, respectively.

(2) : the exponent β reflects the total number of oxygen vacancies incorporated into or directly participating in the reactant species during the RDS. When zero, one, or two oxygen vacancies are involved, β = 0, 1, or 2, respectively. For cathode materials following the dehydrogenation mechanism, no oxygen vacancies participate (β = 0), whereas for dehydration-type reactions, two oxygen vacancies can be involved via the simple ORR process (eqn (3.1)). As summarized in Table S1, a total of seven distinct configurations are possible; however, β is determined solely by the number of oxygen vacancies participating during the RDS.

(3) h˙: the exponent γ depends on the total number of electronic holes incorporated into the reactant species during the RDS, and can take values of 0, −1, −2, or −3. Holes generated as products of the RDS are not counted, as they do not affect the reaction order. From the overall reaction (eqn (2)), up to four holes can, in principle, participate. In practice, however, after the incorporation of two or three holes, H₂O species are released, making intermediates containing all four incorporated holes thermodynamically implausible.

(4) : the exponent ω corresponds to the total number of protons incorporated into or participating in the reactant species during the RDS, and can take values of 0, 1, or 2. Although up to four protons can participate in the overall reaction (eqn (2)), once two protons are incorporated, the H₂O species are released; thus, intermediates containing three or four incorporated protons cannot exist. When all possible proton involvement cases—before, during, and after the RDS—are enumerated, a total of 15 configurations can be identified. Nevertheless, only the number of protons incorporated during the RDS determines the reaction rate, and therefore only the three cases (ω = 0, 1, 2) need to be considered.

2.3. Proton uptake mechanism and kinetics: electrical conductivity relaxation

A simple method to distinguish whether a hole–conducting cathode material follows the hydration reaction or the hydrogenation reaction as its proton uptake mechanism is to examine whether the total electrical conductivity (≈hole conductivity) varies with the water partial pressure (pH₂O).²⁰ In the case of hydration/dehydration, the incorporated protons and oxygen vacancies mutually compensate each other's charge, resulting in no change in the hole concentration. In contrast, in the hydrogenation/dehydrogenation mechanism, charge neutrality is maintained between protons and holes, leading to a change in hole concentration. Consequently, as illustrated in Fig. 2a and b, for materials following the hydration mechanism, the equilibrium total conductivity remains essentially constant with increasing pH₂O. Conversely, for materials governed by the hydrogenation mechanism (Fig. 2c and d), the equilibrium total conductivity decreases relative to its initial value at low pH₂O, owing to the reduction in hole concentration.

Fig. 2 Schematic illustration of surface reaction and charge transport processes and corresponding electrical conductivity relaxation (ECR) behaviors for hydration-type (a and b) and hydrogenation type (c and d) proton uptake mechanisms. Each panel shows typical relaxation behavior under surface reaction-limited (a and c) and bulk diffusion-limited (b and d) cases.

Furthermore, the electrical conductivity relaxation (ECR) behavior observed upon a sudden increase in pH₂O provides insight into whether the overall reaction is dominated by surface reaction kinetics or by bulk diffusion. For hydration-type materials, if the overall process is surface-reaction-limited, the conductivity remains nearly unchanged over time (Fig. 2a). When bulk diffusion dominates, however, differences in the diffusivities of protons and oxygen ions—that is, between proton–hole and oxygen-vacancy–hole pairs—lead to a twofold relaxation curve (Fig. 2b).²¹

In contrast, for hydrogenation-type materials, both the surface-reaction-limited and diffusion-limited cases exhibit monotonic relaxation behavior (Fig. 2c and d). The shape of the relaxation curve differs between these two regimes and can be used to distinguish them by comparison with surface-reaction-limited, bulk-diffusion-limited, or mixed (surface and bulk co-limited) kinetic models.²²

As shown in Fig. 3, the intrinsic proton uptake mechanisms and corresponding kinetics of the PBSCF and BCFZY cathodes were investigated experimentally using electrical conductivity measurements and electrical conductivity relaxation (ECR) techniques. For PBSCF, no noticeable change in total conductivity was observed when the atmosphere was abruptly switched from dry to 3% humidified condition (under a constant pO₂ of 0.17 atm) at 650 °C and 550 °C. This indicates that PBSCF follows the hydration-type mechanism and that its proton uptake process is surface-reaction-limited. One might attribute the pH₂O-independent conductivity of PBSCF to its intrinsically low proton uptake ability. Indeed, as reported previously, the proton concentrations in PBSCF are approximately 0.5 mol% at 550 °C and 0.28 mol% at 650 °C. However, if PBSCF undergoes a hydrogenation process rather than hydration, even such small variations in proton concentration would be expected to induce a total conductivity change on the order of several tens of S cm⁻¹. This is because a change in proton concentration would directly alter the hole concentration, and the hole mobility in PBSCF is relatively high (≈0.7–0.9 cm² V⁻¹ s⁻¹).²⁰ We performed additional ECR measurements under a higher 10% humidified condition (pH₂O = 0.1 atm) to more clearly confirm that the total conductivity is independent of pH₂O, and likewise observed no detectable change in conductivity (Fig. S1), confirming that PBSCF exhibits a surface-reaction-limited hydration process rather than an intrinsically negligible proton uptake capability.

Fig. 3 Conductivity response of PBSCF and BCFZY samples under gas switching from dry to 3% H₂O-humidified atmosphere (indicated by the red line) at 650 °C and 550 °C. (a) PBSCF (pO₂ = 0.17 atm), showing negligible change in conductivity upon proton uptake (b) BCFZY (pO₂ = 0.156 atm), exhibiting a clear decrease in conductivity upon proton uptake.

In contrast, for BCFZY, the total conductivity decreased when the atmosphere was switched from dry to 3% humidified condition (under a constant pO₂ of 0.156 atm), and the conductivity evolved monotonically toward equilibrium (Fig. 3b and 4a). This behavior clearly indicates that the overall proton uptake mechanism of BCFZY involves a significant contribution from the hydrogenation reaction. To be noted, monitoring changes in the total conductivity upon hydration and hydrogenation is intended to indirectly track variations in the hole concentration. Therefore, to eliminate the influence of potentially different hole mobilities among materials, it is more appropriate to consider a normalized metric (e.g., the conductivity change relative to the total conductivity, Δσ/σ) rather than the absolute change in total conductivity.

Fig. 4 (a) Electrical conductivity of BCFZY under dry and 3% H₂O-humidified atmosphere (pO₂ = 0.156 atm) at various temperatures. The conductivity decreases under humidified conditions, providing the evidence of hydrogenation behavior. (b) Surface exchange coefficient (k_chem) and chemical diffusion coefficient (D_chem) of BCFZY as a function of temperature, extracted by fitting ECR data obtained during dehydrogenation.

To further quantify the fraction of protons incorporated via the hydration pathway—that is, the degree of hydration—we applied the analysis protocol previously proposed in the literature.²⁰ In this method, under the assumption that the mobility of charge carriers is independent of partial pressures at a fixed temperature, the degree of hydration is estimated from the measured total conductivity (dominated by holes) and sample mass as functions of both pO₂ and pH₂O (see SI Note S1 for details). Using this protocol, the degree of hydration of BCFZY at 550 °C was determined to be 0.84 ± 0.01, indicating that 16% of the total proton uptake occurs through the hydrogenation reaction.

Regarding the kinetics of the proton uptake process, analysis of the relaxation curves in Fig. 3b revealed that, at 650, 600, and 550 °C, the overall process is co-limited by surface reaction and bulk diffusion. From this analysis, the chemical diffusion coefficient (D_chem) and the surface reaction rate constant (k_chem) associated with the hydrogenation process were extracted (Fig. 4b and Fig. S2).

It is worth noting two important points. First, the characteristic length (defined as L = D_chem/k_chem) was found to be approximately 60–70 µm in the temperature range of 550–650 °C. Considering that porous cathodes in practical electrochemical cells typically have characteristic thicknesses much smaller than this value, the overall proton uptake process in real devices is expected to be surface-reaction-limited. Second, and most remarkably, the surface reaction in the hydrogenation/dehydrogenation process of BCFZY is exceptionally fast. A comparison of surface reaction rate constants (k_chem) for the proton-independent ORR in representative SOFC cathodes highlights this point. At 600 °C, the k_chem value for BCFZY hydrogenation is 1.09 × 10⁻³ cm s⁻¹, whereas those for the proton-independent ORR are 2.34 × 10⁻⁴ cm s⁻¹ for Ba_0.5Sr_0.5Co_0.8Fe_0.2O_3−δ (BSCF)²³ and 7.96 × 10⁻⁵ cm s⁻¹ for PrBa_0.5Sr_0.5Co_1.5Fe_0.5O_5+δ (PBSCF).²⁴ For reference, the k_chem of BCFZY for its proton-independent ORR was measured to be 3.27 × 10⁻⁵ cm s⁻¹ (Fig. S3), which is less than half that of PBSCF. This striking contrast clearly demonstrates that the hydrogenation reaction in BCFZY proceeds orders of magnitude faster than its proton-independent ORR counterpart.

As a result, both PBSCF and BCFZY can be considered to proceed via the extended surface pathway of the proton-coupled ORR, as illustrated in Fig. 1. For PBSCF, the reaction follows solely the dehydration mechanism, whereas for BCFZY, both dehydrogenation and dehydration mechanisms operate concurrently. Although the individual kinetic parameters for these two processes in BCFZY cannot be fully resolved, the exceptionally fast surface reaction kinetics observed for the dehydrogenation pathway suggest that this mechanism may, in fact, dominate under operating conditions. If this is the case, the proton-independent ORR (eqn (3.1)) would not be required, and the overall cathodic reaction could proceed without direct involvement of surface oxygen vacancies.

To express the global reaction rate of the RDS (eqn (4)), it is essential to determine the pO₂ and pH₂O dependencies of the oxygen vacancy, hole, and proton concentrations. Because direct experimental measurement of the surface concentrations of these three defect species is extremely challenging, one must assess whether bulk defect data can be used instead—that is, whether the pO₂ and pH₂O dependencies of the surface defect concentrations can be assumed equivalent to those in the bulk.

As discussed in SI Note S2, cathode materials typically employed in PCFCs are heavily doped and highly redox-active, meaning that the electrical potential difference between the bulk and surface is relatively insensitive to changes in gas partial pressure. Fleig et al. employed near-ambient pressure X-ray photoelectron spectroscopy (NAP-XPS) to study three p-type conducting oxides, which were La_0.6Sr_0.4CoO_3−δ (LSC), La_0.6Sr_0.4FeO_3−δ (LSF), and SrTi_0.7Fe_0.3O_3−δ (STF), and showed that the electrical potential difference between the bulk and the surface is independent of overpotential and pO₂ under oxidizing conditions.²⁵ Therefore, although the absolute defect concentrations at the surface may differ from those in the bulk, their partial-pressure dependencies can be considered nearly identical.¹³ This assumption allows the use of bulk defect concentration dependencies to represent surface behavior.

Since both PBSCF and BCFZY are hole-dominant conductors, the local pO₂ and pH₂O dependencies of the hole concentration ([h˙]), assuming constant mobility, can be determined experimentally by measuring variations in total electrical conductivity under small perturbations of pO₂ and pH₂O. The underlying assumption here is that the mobility of holes is independent of both partial pressures. For BCFZY, the dependencies of the hole concentration on pO₂ and pH₂O can be obtained from equilibrium electrical conductivity data shown in Fig. 3b (pH₂O variation at 550 °C) and Fig. S4b (pO₂ variation at 550 °C), respectively. Subsequently, using the equilibrium relationships described by eqn (3.1) and (3.2), the corresponding pO₂ and pH₂O dependencies of oxygen vacancies and protons can be derived as follows. (see SI Note S5 for details).

For PBSCF, the corresponding partial-pressure dependencies of the defect concentrations at 550 °C can be derived from the same experimental data previously reported in our earlier study.²⁰

2.4. Electrochemical analysis

A 2 × 2 cm² PCFC single cell was fabricated using a tape-based process, consisting of a Ni–BaCe_0.4Zr_0.4Y_0.1Yb_0.1O_3−δ (BCZYYb) anode, a BCZYYb electrolyte, and either a PBSCF or BCFZY cathode. The overall fabrication procedure, photographs of the PCFCs, and cross-sectional microstructures of both cells are shown in Fig. S12.

The electrochemical characteristics of the PCFCs employing PBSCF and BCFZY cathodes were analyzed by electrochemical impedance spectroscopy (EIS) at 550 °C under various conditions: (i) as a function of pO₂ at the cathode, (ii) as a function of pH₂O at the cathode, and (iii) as a function of pH₂ at the anode. Fig. 5 presents the resulting impedance spectra in Nyquist plots: (a and d) show the effect of varying pO₂ and (b and e) the effect of pH₂O at the cathode (with a constant flow of 3% humidified H₂ at the fuel electrode), while (c and f) display the effect of pH₂ at the anode (with a constant flow of dry air at the cathode).

Fig. 5 Electrochemical impedance spectra (Nyquist plots) of PCFCs with PBSCF (top row) and BCFZY (bottom row) cathodes measured at 550 °C under different gas conditions. (a and d) Effect of varying pO₂ at the cathode, (b and e) effect of pH₂O at the cathode (with a constant flow of 3% humidified H₂ at the anode), and (c and f) effect of pH₂ at the anode (with a constant flow of dry air at the cathode). Nyquist plots correspond to PBSCF (a–c) and BCFZY (d–f) cathodes, respectively.

To further resolve the individual electrochemical processes, a distribution of relaxation time (DRT) analysis was performed. In all cases, five distinct peaks (denoted as R1–R5) were observed in the regularization parameter γ(τ) as a function of frequency, as shown in Fig. 6, indicating that the total polarization resistance consists of at least five contributions, each corresponding to a distinct elementary reaction step.

Fig. 6 Distribution of relaxation time (DRT) spectra of PCFCs with PBSCF (top row) and BCFZY (bottom row) cathodes measured at 550 °C under various gas conditions. (a and d) Effect of varying pO₂ at the cathode, (b and e) effect of pH₂O at the cathode (with a constant flow of 3% humidified H₂ at the anode), and (c and f) effect of pH₂ at the anode (with a constant flow of dry air at the cathode). Five distinct relaxation processes (R1–R5) are identified, indicating that the total polarization resistance consists of multiple elementary reaction steps.

When the pO₂ at the cathode was increased, both the PBSCF- and BCFZY-based cells exhibited a pronounced decrease in polarization resistance in the low-frequency region (frequency (=1/τ) ≈ 10¹–10² Hz) (Fig. 6a and d). In contrast, varying the pH₂ at the anode produced negligible changes in the impedance spectra for either cell (Fig. 6c and f). A notable distinction between the two cathode materials emerged upon varying the pH₂O at the cathode. For PBSCF, the polarization resistance remained virtually unchanged with increasing pH₂O (Fig. 6b), whereas for BCFZY, the polarization resistance in the low-frequency region (1/τ ≈ 10¹–10² Hz) decreased markedly with higher pH₂O, indicating enhanced cathodic kinetics.

To quantitatively evaluate these trends, the area-specific polarization resistances (ASRs) corresponding to R1–R5 were calculated, as shown in Fig. 7. The calculation was based on the principle that the fractional area of each relaxation peak (R1–R5) with respect to the total area represents the ratio of each process's polarization resistance to the total polarization resistance. Consistent with Fig. 6, the R1 component (centered near 10 Hz) accounted for the largest portion of the total polarization for both cells, and its variations with pO₂ and pH₂O were most pronounced.

Fig. 7 Area-specific polarization resistances (ASRs) of individual relaxation processes (R1–R5) derived from DRT analysis of PCFCs employing PBSCF and BCFZY cathodes at 550 °C under different gas conditions. (a) Dependence on pO₂ at the cathode, (b) dependence on pH₂O at the cathode (3% humidified H₂ at the anode), and (c) dependence on pH₂ at the anode (dry air at the cathode). The fractional area of each relaxation component (R1–R5) represents its relative contribution to the total polarization resistance for each cell.

By plotting the derived polarization resistances against pO₂ and pH₂O, the corresponding partial-pressure exponents (A, B) in (eqn (7)) were determined, as presented in Fig. 8.

where R(i) denotes the area-specific polarization resistance of component i. From the exponents extracted in Fig. 8, the reaction rate of the dominant PBSCF cathodic process R1 is proportional to pO₂^0.24±0.03·pH₂O^0.00±0.01. Cathodic processes R1 and R2 in BCFZY are proportional to pO₂^0.29±0.02·pH₂O^0.44±0.10 and pO₂^0.07±0.09·pH₂O^0.06±0.03, respectively. Next, we determine the exponents of the defect activities in (eqn (4)) by using the partial-pressure dependencies of the relevant point defects derived in (eqn (5)) for BCFZY and (eqn (6)) for PBSCF to obtain the identical numerical values of exponents A and B. With those dependencies fixed, and referencing each material's proton-uptake mechanism (hydration vs. hydrogenation), the admissible combinations listed in Tables 1–4 allow us to assign (α, β, γ, ω) in , and, among the feasible scenarios, identify the most plausible RDS.

Fig. 8 Dependence of area-specific polarization resistances (ASRs) for individual relaxation components (R1–R5) on pO₂ and pH₂O at 550 °C for PCFCs with PBSCF (a and b) and BCFZY (c and d) cathodes (oxygen electrode side). (a and c) log (ASR) as a function of log(pO₂), showing the influence of oxygen partial pressure on each relaxation component. (b and d) log (ASR) as a function of log (pH₂O), highlighting the dependence of polarization processes on water vapor content at the cathode.

Within the temperature, pO₂, and pH₂O ranges investigated in this work, both PBSCF and BCFZY exhibited negligible pH₂O dependences of the hole and oxygen vacancy concentrations (eqn (5) and (6)), whereas only the proton concentration showed a strong pH₂O dependence. Consequently, if the RDS is an elementary step that does not include protons as reactants, the overall rate is expected to show negligible pH₂O dependence, because the concentrations of the reactant point defect species entering the RDS are themselves nearly independent of pH₂O. In other words, it is not the mere presence or absence of protons in the RDS that dictates the pH₂O dependence. Rather, the pH₂O dependence of the point defect species participating in the RDS as reactants ultimately determines the reaction order with respect to pH₂O.

Based on this understanding, from the distinct pH₂O dependencies observed for PBSCF and BCFZY, it is possible to first determine whether protons are involved in the RDS. For PBSCF, the absence of any measurable pH₂O dependence indicates that the RDS is proton-independent, corresponding to one of the steps in the proton-independent ORR pathway. In contrast, for BCFZY, the clear dependence on pH₂O confirms that the RDS involves proton participation.

Specifically, for PBSCF, as summarized in Table 5, two possible cases were identified for the R1 component—considered the most likely RDS—whose reaction rate exponents with respect to pO₂ and pH₂O closely match the experimentally derived values. Based on the corresponding defect concentration exponents, possible reaction expressions were derived for each case. Comparison of these two cases (Case I and Case II) reveals that the key distinction lies in whether holes participate in the RDS.

Table 5 Exploration of possible cases (Case I and II) yielding reaction exponents (α, β, γ, ω) consistent with the experimentally derived rate expression of the RDS (R1) in PBSCF, along with the corresponding possible reaction expressions derived from these exponents

R1 in PBSCF							Possible reaction expressions
	α	β	γ	ω	A	B
Case I	1	2	−1	0	0.224	0
Case II	1	2	0	0	0.300	0
Experimental result	n/a				0.24 ± 0.03	0.00 ± 0.01	—

---

Because the hole concentration in PBSCF exhibits only a weak dependence on pO₂ (∝pO₂^0.076), it is difficult to conclusively identify which case is more plausible solely by comparing exponent values. Indeed, the literature presents differing interpretations. D. Poetzsch et al. suggested that, given the much lower bonding energy of peroxide species (∼1.5 eV) compared to the O–O double bond (∼5 eV), the peroxide form is more likely to exist.¹³ In contrast, D. Chen et al. experimentally demonstrated that in Pr-doped ceria, the RDS of the ORR corresponds to the dissociation of neutral molecular oxygen adsorbates in which electrons do not participate.⁵ In any cases, both share a common feature: the RDS involves the dissociation of diatomic oxygen molecules, with two oxygen vacancies participating in the reaction. This diatomic oxygen dissociation process thus represents the most plausible rate-determining step for the PBSCF cathode. From our atomistic simulation results (see the following Section S2.5), the pathway involving charged oxygen intermediates is energetically more favorable than with hole-free dissociation of neutral molecular oxygen.

As summarized in Table 6, analysis of the R1 component of BCFZY—where the RDS involves proton participation—yielded three possible cases whose reaction rate exponents with respect to pO₂ and pH₂O closely match the experimentally observed values. Among these, two chemically feasible reaction expressions (Cases I and III) were derived. The key distinction between Case I and Case III lies in whether oxygen vacancies participate in the ORR: in Case III, an oxygen vacancy serves as the adsorption site for diatomic oxygen. If the reaction indeed followed Case III, then the advantage of the dehydrogenation mechanism—namely, that the overall reaction could proceed without requiring a separate proton-independent ORR—would no longer hold, since oxygen vacancies would again be directly involved.

Table 6 Exploration of possible cases yielding reaction exponents (α, β, γ, ω) consistent with the experimentally derived rate expression of the RDSs (R1 and R2) in BCFZY, along with the corresponding possible reaction expressions derived from these exponents

R1 in BCFZY					pO2^A·pH2O^B		Possible reaction expressions
	α	β	γ	ω	A	B
Case I	1/2	0	−1	1	0.250	0.500
Case II	1/2	1	−1	1	0.222	0.498	No valid reaction step
Case III	1	1	−3	1	0.250	0.500
Experimental result	n/a				0.29 ± 0.02	0.44 ± 0.10	—

---

R2 in BCFZY					pO2^A·pH2O^B		Possible reaction expressions
	α	β	γ	ω	A	B
Case I	1/2	0	−2	0	0.028	0.002	O^2−ad → O^2−ad,TPB
Case II	1/2	1	−2	0	0	0	No valid reaction step
Experimental result	n/a				0.07 ± 0.09	−0.06 ± 0.03	—

---

Further analysis of the R2 component, which exhibits a polarization magnitude comparable to R1, revealed that it corresponds to the surface diffusion of atomic oxygen adsorbates. Considering that the R1 and R2 processes are sequential reaction steps, Case I—where oxygen vacancies are not involved—is more plausible than Case III. This conclusion is consistent with the findings of M. Shang et al., who first analyzed the RDS of BCFZY cathodes and similarly identified a proton-involved surface reaction mechanism as the dominant pathway.²⁶

Furthermore, by applying the methodology established in this study to 600 °C, we confirmed that the RDS for both materials remains identical to that observed at 550 °C (see SI Note S4).

In summary, for PBSCF, the rate-determining step (RDS) corresponds to the dissociation of diatomic oxygen molecules involving surface oxygen vacancies, whereas for BCFZY, the RDS involves the reaction between adsorbed atomic oxygen and bulk protons to form OH adsorbates. Both materials possess sufficient proton transport capability and exhibit rapid bulk diffusion kinetics; thus, in practical porous cathodes, the surface reaction governs the overall cathodic process in both cases. Nevertheless, the fundamentally different RDSs originate from their distinct proton uptake/release mechanisms.

Under fuel-cell operating conditions, PBSCF follows a dehydration-type mechanism, which necessitates an additional sequential proton-independent ORR. Because this ORR step is kinetically slower than the dehydration reaction itself, it becomes the rate-limiting process. This interpretation is consistent with the approximately 14-fold smaller surface reaction rate constant (k_chem) for the proton-independent ORR in PBSCF compared with the hydrogenation reaction in BCFZY.

These insights suggest distinct strategies for enhancing cathodic activity. For hydration-type materials such as PBSCF, improving performance requires increasing the surface oxygen vacancy concentration and introducing catalysts or dopants that facilitate diatomic oxygen dissociation. In contrast, for hydrogenation-type materials such as BCFZY, optimizing performance involves enhancing the degree of hydrogenation and bulk proton conductivity, thereby suppressing the surface-diffusion limitation of adsorbed atomic oxygen.

While the potential for an RDS shift under overpotential remains, the theoretical framework by Z. Guan et al.²⁷ implies that the current methodology retains its interpretability under bias. However, practical evaluations face inherent constraints, particularly electronic leakage in the electrolyte. As noted by C. Duan et al.²⁸ and further confirmed by H. Sumi et al.²⁹ such leakage can complicate analysis under overpotential. In this regard, integrating complementary strategies, such as the isotope effect analysis proposed by Y. Okuyama et al.¹² would enable a more definitive identification of the RDS under such operating conditions (A detailed discussion on this matter is provided in SI Note S6).

2.5. Atomistic simulations

To gain further insights into the potential RDS steps in Tables 5 and 6, we carried out quantum-mechanical level atomistic simulations using a universal machine-learning interatomic potentials (MLIPs).^30–32 Among universal MLIPs, we used SevenNet-Omni, which demonstrating state-of-the-art accuracy in reproducing adsorption energetics and reaction barriers for catalyst systems.³³ We further validated its applicability to PBSCF and BCFZY by benchmarking representative reaction energetics against DFT reference data (see Fig. S8).

We first investigate the reaction energies of potential RDS reaction on PBSCF (Table 5: Case I: hole-participating dissociation vs. Case II: hole-free dissociation). Fig. 9a shows the adsorption energies of O₂ adsorbed at one oxygen vacancy in the presence of a neighboring vacancy, together with the corresponding Bader charge of the adsorbed O₂ (q_ad). Various vacancy configurations were examined, and for clarity we present the configurations exhibiting the highest q_ad values (closest to neutral O₂) as well as the lowest-energy configurations, for both Co–Fe-terminated and Co–Co-terminated surfaces (see Fig. S9 for all configurations). The adsorption energy becomes more negative with increasing charge transfer, indicating that charged oxygen species are energetically more stable than neutral molecular oxygen. In addition, we find that the activation barrier for the transition from neutral-like O₂ to a superoxide-like species is very small (0.03 eV; not shown), suggesting that neutral molecular oxygen is unlikely to persist as a kinetically relevant intermediate. Finally, we evaluated the activation barrier for the subsequent oxygen dissociation reaction involving charged oxygen intermediates, as shown in Fig. 9b, obtaining a value of 0.655 eV. (The top-view configurations are shown in Fig. S10b and c) Taken together, these results suggest that oxygen dissociation on PBSCF is more consistent with a pathway involving charged oxygen intermediates than with hole-free dissociation of neutral molecular oxygen, thereby providing computational support for Case I over Case II.

Fig. 9 (a) Representative atomic structures of oxygen molecule adsorption. The q_ad indicated the atomic charges of O–O moiety from obtaining Bader charge analysis, and the E_ads is adsorption energy. (b–d) Representative reaction pathways and their atomic configurations from computational results: (b) Oxygen molecule dissociation process involving oxygen vacancy on the defective PBSCF surface. (c) Proton transfer mechanism for BCFZY. Proton in subsurface of BCFZY slab exposes on the surface, and moves to oxygen adsorbate. (d) HO₂ dissociation mechanism on the defective BCFZY surface. Dashed line circle indicates the oxygen vacancy, which involves the ORR pathway. O_react(blue) indicate oxygen atoms which directly take part in the reactions.

For BCFZY, we calculated the reaction energies and activation barriers of (Case I) and (Case III), respectively in Fig. 9c and d (see Table 6 for details). (The top-view configurations are shown in Fig. S10c and d) For the first pathway, the reaction proceeds via two sequential steps. First, a subsurface oxygen migrates to the surface with an activation barrier of 0.041 eV. Subsequently, this oxygen species transfers to the adsorbed oxygen (O_ad) with an activation barrier of 0.126 eV. In contrast, dissociation of the HO₂ species at the same active site exhibits a substantially higher activation barrier of 0.533 eV. These results indicate that the first pathway is kinetically more favorable than the second. This finding is consistent with our interpretation in the previous subsection that the R1 and R2 processes are sequential reaction steps and that oxygen vacancies are unlikely to be directly involved in the rate-determining step for BCFZY.

3. Conclusion

In this study, we established a systematic theoretical and experimental framework for elucidating the reaction pathway of the proton-coupled oxygen reduction reaction (PC-ORR) in PCFC cathodes. Using PBSCF and BCFZY as representative examples, we identified their RDSs by correlating intrinsic thermodynamic and kinetic properties with a generalized microkinetic model. This newly proposed approach provides a fundamental basis for a unified understanding of PCFC cathode reactions and for designing activity-enhancement strategies. Beyond PCFCs, the methodology presented here can be broadly applied to air electrodes in SOFCs and SOECs, offering valuable guidance for developing high-efficiency electrochemical energy-conversion devices.

4. Experimental section

4.1. Fabrication of protonic ceramic fuel cells and electrochemical analysis

For the fuel electrode (anode) support layer (FSL), fuel electrode functional layer (FFL), and electrolyte layer (EL), tapes were prepared using the tape casting method. The FSL tape consisted of NiO (Mechema, Taiwan), BaCe_0.4Zr_0.4Y_0.1Yb_0.1O_3−δ (BCZYYb) (Terrafuelcell, Korea), and poly(methyl methacrylate) (PMMA, Sunjin Beauty Science, Korea) pore-forming agent in a weight ratio of 61 : 34 : 5, dissolved in solvents with dispersant, binder, and plasticizer. The FFL tape shared an identical composition with the FSL tape, excluding the PMMA pore former, and the EL tape followed the same procedure as the FSL tape. The tape casting process for FSL, FFL, and EL layers was conducted by FCT (Fine Ceratech, Korea), resulting in approximate thicknesses of 85 µm for FSL, 9–25 µm for FFL, and 14 µm for EL.

The fuel electrode/electrolyte assembly was fabricated by sequential stacking and roll calendering. Approximately seven sheets of FSL and one sheet of FFL tape were first stacked, with Mylar films attached to the top and bottom surfaces, and then pressed using a roll calender at 75 °C between the upper and lower rolls, with the roll gap adjusted to approximately 2% of the total thickness of the stacked tapes and films. The resulting laminated FSL/FFL body served as the fuel electrode substrate. Subsequently, the EL tape was placed on top of this laminated body and pressed with the same method at 75 °C, with the roll gap adjusted to approximately 3%, to achieve uniform adhesion between the FFL and EL layers. The final laminated FSL/FFL/EL body was trimmed along the edges to a size of 7.5 × 7.5 cm², sintered at 1450 °C for 4 h, and then cut into 2 × 2 cm² specimens for electrochemical evaluation.

The PrBa_0.5Sr_0.5Co_1.5Fe_0.5O_5+δ (PBSCF) oxygen electrode (cathode) ink was prepared by blending the PBSCF powder (Kceracell, Korea) with α-terpineol as the solvent, KD-6 (3 wt% relative to the powder) as the dispersant, BH-3 (3.5 wt%) as the binder, and dibutyl phthalate (DBP, 3.5 wt%) as the plasticizer using a planetary milling machine. The resulting ink was then screen printed onto the sintered electrolyte (effective area: 1 × 1 cm²) and sintering at 950 °C for 5 hours.

The BCFZY oxygen electrode ink was prepared using BCFZY powder (Terra Fuel Cell, Korea) following the same fabrication procedure as that described above for the PBSCF ink.

The cross-sectional images of the sintered single cell with the oxygen electrode were examined using a scanning electron microscope (SEM, Inspect F50, FEI) to evaluate the microstructure. Electrochemical impedance measurements were conducted at 550 °C in a lab-made testing station. The cell was placed between Inconel 600 alloy metal interconnectors and compression-sealed using a glass-ceramic sealant. Au mesh and Ni foam served as the current collectors for the oxygen and fuel electrodes, respectively. Impedance spectra were analyzed under open-circuit voltage (OCV) conditions using a frequency response analyzer and a potentiostat (Solartron 1260/1287) with an amplitude of 10 mV over a frequency range of 10⁶ to 10⁻¹ Hz.

The electrode reactions were investigated by analyzing impedance spectra as a function of gas compositions at both the oxygen and fuel electrodes. At the oxygen electrode, (i) the oxygen partial pressure (pO₂) was controlled within −1.20 ≤ log(pO₂/atm) ≤ −0.68 for the PBSCF cell and −1.13 ≤ log(pO₂/atm) ≤ −0.68 for the BCFZY cell by adjusting the mixing ratio of air and N₂, and (ii) the water vapor partial pressure (pH₂O) was controlled within −1.54 ≤ log(pH₂O/atm) ≤ −0.64 for the PBSCF cell and −1.60 ≤ log(pH₂O/atm) ≤ −0.70 for the BCFZY cell by altering the mixing ratio of water vapor and non-humidified (dry) air. Throughout these variations at the oxygen electrode, 3% H₂O-humidified H₂ was continuously supplied to the fuel electrode side.

For the gas variation at the fuel electrode, the hydrogen partial pressure (pH₂) was controlled within −0.54 ≤ log(pH₂/atm) ≤ −0.01 for the PBSCF cell and −0.47 ≤ log(pH₂/atm) ≤ −0.01 for the BCFZY cell by modifying the mixing ratio of H₂ and N₂. Throughout these variations at the fuel electrode, non-humidified air was consistently supplied to the oxygen electrode side.

4.2. Preparation of BCFZY and PBSCF samples for ECR measurements

BCFZY powder was synthesized via the glycine-nitrate process (GNP) with a glycine-to-nitrate molar ratio of 0.95 : 1. Stoichiometric amounts of Ba(NO₃)₂ (≥99.0%, Sigma-Aldrich), Co(NO₃)₂·6H₂O (≥98.0%, Sigma-Aldrich), Fe(NO₃)₃·9H₂O (≥98.0%, Sigma-Aldrich), zirconyl nitrate solution (35 wt% in diluted nitric acid, Sigma-Aldrich), Y(NO₃)₃·6H₂O (≥99.8%, Sigma-Aldrich) were dissolved in deionized water. After complete dissolution, glycine was added to the solution. The mixed precursor solution was heated at 80 °C to evaporate the solvent and form a gel, which was then combusted at 500 °C to yield a fine BCFZY powder.

The resulting powder was ball-milled in ethanol, dried, and subsequently calcined at 800 °C in a box furnace. To obtain a fine and uniform powder, planetary wet milling was conducted in ethanol, followed by drying prior to pellet formation. The powder was uniaxially pressed at 15 MPa into a rectangular mold (20 mm × 15 mm) for 1 min, followed by cold isostatic pressing (CIP) at 200 MPa for 3 min to enhance densification. The green pellets were sintered at 1275 °C to obtain dense rectangular bodies with approximately 95% of the theoretical density, as measured using the Archimedes method.

After cutting the sintered pellets to the required dimensions using a low-speed diamond saw, its surfaces were ground and polished to produce flat surfaces with minimal roughness. This preparation step ensured well-defined sample geometry for 2D analysis, enabling precise evaluation of surface kinetics. The final specimens had dimensions of 0.80 × 0.80 × 14.35 mm³.

PBSCF samples for ECR was fabricated via tape casting method using a slurry prepared from commercial PBSCF powder (Kceracell, Korea), as described in previous study.²⁰ The resulting specimens had a rectangular geometry with dimensions of 10.4 × 5.4 × 0.36 mm³.

4.3. Electrical conductivity relaxation (ECR) measurements setup & method

ECR measurements were conducted using a conventional four-probe DC method. Two current leads were positioned at the outer ends of the bar-shaped sample, and two voltage leads were placed near the center. Each wire was tightly bound to the sample surface to ensure reliable electrical contact.

A custom-built jig was designed to keep a stable atmosphere environment and allow quick changes between gas atmospheres. The jig consisted of a small-volume alumina end cap, a six-hole alumina rod, and a metallic holder with gas ports. The alumina endcap formed an enclosed space to stably maintain the desired gas composition around the sample, and its compact volume enabled rapid changes in gas atmosphere. The six-hole alumina rod contained four Pt wires for electrical connections and an S-type thermocouple positioned near the sample. Its multi-channel structure physically separated all leads, preventing contact between the thermocouple and Pt wires. The alumina assembly was mounted on a metallic holder equipped with gas inlet and outlet ports to ensure stable and continuous gas flow during ECR measurements. The sample was heated in a furnace to the target temperature, and the temperature was monitored by a thermocouple. Electrical conductivity was continuously recorded in real time using DC current source (6220, Keithley, USA) and a digital multimeter (2000, Keithley, USA).

ECR measurements were carried out to investigate proton uptake and ORR kinetics at 550, 600 and 650 °C. For the proton uptake experiment, a dry gas with a pO₂ of 0.156 atm was prepared by blending high-purity argon and oxygen. To achieve the humidified condition (3% H₂O), a gas mixture initially containing 16.07% O₂ was passed through a water saturator immersed in a chiller bath held at 25 °C, yielding a humidified gas with pO₂ matched to the dry condition (0.156 atm). For ORR kinetics measurements, gas mixtures with pO₂ values of 0.156 atm and 0.123 atm were used. The total gas flow rate was maintained at 300 sccm for all conditions.

For PBSCF, ECR measurements were conducted under slightly different conditions, as described in a previous study.²⁰ Proton uptake was investigated by switching the gas atmosphere between dry and 3% humidified conditions at a constant pO₂ of 0.17 atm. To evaluate ORR kinetics, measurements were conducted under pO₂ conditions of 0.17 and 0.08 atm. All oxygen partial pressure conditions were precisely controlled using a pO₂ measurement sensor (Nanoionics, Korea).

At the beginning of each ECR test, the sample was stabilized under the designated temperature and initial gas atmosphere until the conductivity stabilized at a steady value. After switching to the new gas condition, conductivity was continuously monitored until it remained stable over time. To enable rapid gas switching, a 4-way valve was installed upstream of the jig. Two gas mixtures were continuously supplied in parallel - one directed into the jig and the other vented just before the inlet. The vented gas was continuously monitored using a pO₂ gas sensor. Switching the valve allowed immediate redirection of the gas flow without delay, enabling fast atmosphere exchange within the jig.

4.4. ECR measurements curve fitting

The electrical conductivity relaxation curves were analyzed based on a two-dimensional (2D) diffusion model that incorporates both bulk diffusion and surface exchange processes. This analytic expression, derived from Fick's second law and incorporating appropriate initial and boundary conditions, and is given as follow³⁴where the β_n are the positive roots of (eqn (7.2)), and dimensionless length L is defined as the ratio k_chem/D_chem.

In the above expression, σ(t) is the measured electrical conductivity at time t, while σ(0) and σ(∞) correspond to the initial and final equilibrium conductivities, respectively. The term a denotes the half-thickness of the sample. The parameters D_chem and k_chem represent the chemical diffusion coefficient and the surface exchange coefficient, respectively. The parameter t₀ is the starting time of conductivity relaxation process.

Initial parameter fitting was conducted using MATLAB scripts to estimate key transport parameters for the conductivity relaxation model. The fitting parameters included the σ(∞), D_chem and five β_n values (where n = 1 to 5), while the σ(0) and a were treated as fixed parameter. To ensure robust optimization, the fitting was performed 3000 times using randomly generated initial value of k_chem and D_chem, from which corresponding β_n values were determined by solving (eqn (7.2)). For each trial, the simulated conductivity curve was compared to experimental data, and the parameter set yielding the lowest residual error was selected as the best fit.

4.5. Computational methods

All DFT calculations are performed useing the Vienna ab initio simulation package (VASP).^35,36 The interactions between the electrons and ions were described using the projector-augmented-wave (PAW) method.^37,38 The exchange–correlation energy was modeled using the Perdew–Burke–Ernzerhof (PBE) type generalized gradient approximation (GGA) functional.³⁹ We considered the onsite Coulomb interaction in Co 3d and Fe 3d orbitals within the GGA+U scheme.^40,41 The effective U values of Co and Fe were set 3.32 and 5.3 eV, respectively.^42,43 We used Grimme's DFT-D3 method for van der Waals corrections.⁴⁴ The kinetic energy cutoff was set at 500 eV, and a Γ-centered 3 × 3 × 3 k-point grid was used for the geometrical optimization in PBSCF and BCFZY bulk structures. The convergence tolerance was set at an energy of 10⁻⁶ eV for electronic structures. In the slab structure, the kinetic energy cutoff was set at as same as the bulk structure calculation and a Γ-centered 3 × 3 × 1 k-point grid was used. Bader charge analysis⁴⁵ was employed to calculate the atomic charge for adsorbed structures.

We used the pre-trained machine learning interatomic potential (MLIP), SevenNet,³³ with atomic simulation environment (ASE) package⁴⁶ for investigating the reaction energies activation energy barriers. The nudged elastic band (NEB) calculations were performed using the MLIP framework with the climbing-image NEB (CI-NEB) method^47,48 to calculate the activation energy barrier. The reaction path was discretized using 7–9 intermediate images (excluding the initial and final states) and a spring constant of 0.1 eV Å⁻². The NEB optimization was converged when the maximum force on any image was below 0.01 eV Å⁻¹.

The optimized atomic structure of PBSCF (BCFZY) bulk state was obtained with lattice parameters of 7.72 × 7.69 × 7.69 ų (8.19 × 8.16 × 8.24 ų) from DFT calculations (see Fig. S4c). The surface index, which is (100) for both, and the surface terminations are followed by ref. 42 and 43. A vacuum layer of 15 Å is inserted for each simulation cell. The bottom two layers are fixed during the optimizations and activation-energy calculations. The slab structures of PBSCF and BCFZY are provided in Fig. S11.

As the surface oxygen vacancies play a key role in the reaction pathway for PBSCF, defective structures containing surface oxygen vacancies were considered. The Gibbs free energy of defective structures for vacancy formation was calculated as:

ΔG_vac = E_defective − E_defect-free + µ_O2(T, p) (8)

where E_defective and E_defect-free are the total energies of defective PBSCF slab with two oxygen vacancies and defect-free PBSCF slab, respectively. Here, µ_O2(T, p) denotes the chemical potential of molecular oxygen. The oxygen chemical potential was calculated by combining DFT total energies with temperature-dependent thermodynamic corrections obtained from the NIST Chemistry WebBook.⁴⁹ The oxygen chemical potential was constructed as:where Δµ_O2(T, p) was taken from NIST. The temperature (T) and oxygen partial pressure (p) were set to 823 K and 0.17 atm, respectively, in consideration of the experimental condition. The three representative atomic configurations of the defective PBSCF are presented in Fig. S5a, along with the corresponding Gibbs free energy. Moreover, the adsorption energy was calculated as follows:

E_ads = E_adsorbed − (E_defective + E_O2) (10)

where E_adsorbed is the total energy of adsorbed structure. The E_O2 is obtained from the most stable gas-phase oxygen molecule.

Conflicts of interest

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

Data availability

Additional data that are not available online can be obtained from the corresponding author upon reasonable request.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5ee06170a.

Acknowledgements

This work was supported by National R&D Program through the National Research Foundation of Korea (NRF) funded by Ministry of Science and ICT (RS-2025-02310288), and the National Research Council of Science & Technology (NST) grant by the Korea government (MSIT) (No. GTL24051-500).

References

1. Z. Gao, L. V. Mogni, E. C. Miller, J. G. Railsback and S. A. Barnett, Energy Environ. Sci., 2016, 9, 1602–1644.
2. S. B. Adler, Chem. Rev., 2004, 104, 4791–4844.
3. W. C. Chueh and S. M. Haile, Annu. Rev. Chem. Biomol. Eng., 2012, 3, 313–341.
4. R. Merkle and J. Maier, Phys. Chem. Chem. Phys., 2002, 4, 4140–4148.
5. D. Chen, Z. Guan, D. Zhang, L. Trotochaud, E. Crumlin, S. Nemsak, H. Bluhm, H. L. Tuller and W. C. Chueh, Nat. Catal., 2020, 3, 116–124.
6. H. Iwahara, Solid State Ionics, 1996, 86, 9–15.
7. K.-D. Kreuer, Annu. Rev. Mater. Res., 2003, 33, 333–359.
8. H. An, H.-W. Lee, B.-K. Kim, J.-W. Son, K. J. Yoon, H. Kim, D. Shin, H.-I. Ji and J.-H. Lee, Nat. Energy, 2018, 3, 870–875.
9. L. R. Tarutina, M. A. Gordeeva, D. E. Matkin, M. T. Akopian, G. N. Starostin, A. V. Kasyanova, A. P. Tarutin, N. A. Danilov, I. A. Starostina and D. A. Medvedev, Chem. Eng. J., 2024, 490, 151615.
10. F. He, T. Wu, R. Peng and C. Xia, J. Power Sources, 2009, 194, 263–268.
11. S. M. Choi, H. An, K. J. Yoon, B.-K. Kim, H.-W. Lee, J.-W. Son, H. Kim, D. Shin, H.-I. Ji and J.-H. Lee, Appl. Energy, 2019, 233, 29–36.
12. Y. Okuyama, K. Kasuga, M. Shimomura, Y. Mikami, K. Yamauchi, T. Kuroha and H. Sumi, J. Power Sources, 2023, 586, 233647.
13. D. Poetzsch, R. Merkle and J. Maier, J. Electrochem. Soc., 2015, 162, F939.
14. S. Choi, C. J. Kucharczyk, Y. Liang, X. Zhang, I. Takeuchi, H.-I. Ji and S. M. Haile, Nat. Energy, 2018, 3, 202–210.
15. C. Duan, J. Tong, M. Shang, S. Nikodemski, M. Sanders, S. Ricote, A. Almansoori and R. O’Hayre, Science, 2015, 349, 1321–1326.
16. M. Papac, V. Stevanović, A. Zakutayev and R. O’Hayre, Nat. Mater., 2021, 20, 301–313.
17. D. Poetzsch, R. Merkle and J. Maier, Faraday Discuss., 2015, 182, 129–143.
18. R. Merkle and J. Maier, Angew. Chem., Int. Ed., 2008, 47, 3874–3894.
19. R. De Souza, Phys. Chem. Chem. Phys., 2006, 8, 890–897.
20. S. Im, M. A. Berk, S. Yang, B.-K. Kim, K. J. Yoon, J.-W. Son, J.-H. Lee and H.-I. Ji, J. Mater. Chem. A, 2022, 10, 16127–16136.
21. H.-I. Yoo, J.-Y. Yoon, J.-S. Ha and C.-E. Lee, Phys. Chem. Chem. Phys., 2008, 10, 974–982.
22. C.-R. Song and H.-I. Yoo, Solid State Ionics, 1999, 120, 141–153.
23. D. Chen and Z. Shao, Int. J. Hydrogen Energy, 2011, 36, 6948–6956.
24. J. Kim, S. Im, S. H. Oh, J. Y. Lee, K. J. Yoon, J.-W. Son, S. Yang, B.-K. Kim, J.-H. Lee and H.-W. Lee, Sci. Adv., 2021, 7, eabj8590.
25. A. Nenning, A. K. Opitz, C. Rameshan, R. Rameshan, R. Blume, M. Hävecker, A. Knop-Gericke, G. N. Rupprechter, B. Klötzer and J. R. Fleig, J. Phys. Chem. C, 2016, 120, 1461–1471.
26. M. Shang, J. Tong and R. O'Hayre, RSC Adv., 2013, 3, 15769–15775.
27. Z. Guan, D. Chen and W. C. Chueh, Phys. Chem. Chem. Phys., 2017, 19, 23414–23424.
28. C. Duan, R. Kee, H. Zhu, N. Sullivan, L. Zhu, L. Bian, D. Jennings and R. O’Hayre, Nat. Energy, 2019, 4, 230–240.
29. H. Sumi, H. Shimada, K. Watanabe, Y. Yamaguchi, K. Nomura, Y. Mizutani and Y. Okuyama, ACS Appl. Energy Mater., 2023, 6, 1853–1861.
30. J. Behler and M. Parrinello, Phys. Rev. Lett., 2007, 98, 146401.
31. C. Chen and S. P. Ong, Nat. Comput. Sci., 2022, 2, 718–728.
32. I. Batatia, P. Benner, Y. Chiang, A. M. Elena, D. P. Kovács, J. Riebesell, X. R. Advincula, M. Asta, M. Avaylon, W. J. Baldwin, F. Berger, N. Bernstein, A. Bhowmik, F. Bigi, S. M. Blau, V. Cărare, M. Ceriotti, S. Chong, J. P. Darby, S. De, F. Della Pia, V. L. Deringer, R. Elijošius, Z. El-Machachi, E. Fako, F. Falcioni, A. C. Ferrari, J. L. A. Gardner, M. J. Gawkowski, A. Genreith-Schriever, J. George, R. E. A. Goodall, J. Grandel, C. P. Grey, P. Grigorev, S. Han, W. Handley, H. H. Heenen, K. Hermansson, C. H. Ho, S. Hofmann, C. Holm, J. Jaafar, K. S. Jakob, H. Jung, V. Kapil, A. D. Kaplan, N. Karimitari, J. R. Kermode, P. Kourtis, N. Kroupa, J. Kullgren, M. C. Kuner, D. Kuryla, G. Liepuoniute, C. Lin, J. T. Margraf, I.-B. Magdău, A. Michaelides, J. H. Moore, A. A. Naik, S. P. Niblett, S. W. Norwood, N. O’Neill, C. Ortner, K. A. Persson, K. Reuter, A. S. Rosen, L. A. M. Rosset, L. L. Schaaf, C. Schran, B. X. Shi, E. Sivonxay, T. K. Stenczel, C. Sutton, V. Svahn, T. D. Swinburne, J. Tilly, C. van der Oord, S. Vargas, E. Varga-Umbrich, T. Vegge, M. Vondrák, Y. Wang, W. C. Witt, T. Wolf, F. Zills and G. Csányi, J. Chem. Phys., 2025, 163, 184110.
33. J. Kim, J. You, Y. Park, Y. Lim, Y. Kang, J. Kim, H. Jeon, S. Ju, D. Hong and S. Y. Lee, arXiv, 2025, preprint, arXiv:2510.11241 DOI:10.48550/arXiv.2510.11241.
34. C. B. Gopal and S. M. Haile, J. Mater. Chem. A, 2014, 2, 2405–2417.
35. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558.
36. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
37. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
38. D. Hobbs, G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 62, 11556.
39. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
40. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. Humphreys and A. P. Sutton, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 1505.
41. V. I. Anisimov, F. Aryasetiawan and A. Lichtenstein, J. Phys.: Condens. Matter, 1997, 9, 767.
42. Y. Zhu, Z. He, Y. Choi, H. Chen, X. Li, B. Zhao, Y. Yu, H. Zhang, K. A. Stoerzinger and Z. Feng, Nat. Commun., 2020, 11, 4299.
43. Y. Wang, Z. Wang, K. Yang, J. Liu, Y. Song, J. Li, Z. Hu, M. J. Robson, Z. Zhang and Y. Tian, Adv. Funct. Mater., 2024, 34, 2404846.
44. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
45. G. Henkelman, A. Arnaldsson and H. Jónsson, Comput. Mater. Sci., 2006, 36, 354–360.
46. A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer and C. Hargus, J. Phys.: Condens. Matter, 2017, 29, 273002.
47. G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9978–9985.
48. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
49. P. Linstorm, J. Phys. Chem. Ref. Data, Monogr., 1998, 9, 1–1951.

---

Footnote

† Equally contributed author.

---