On the origin of multihole oxygen evolution in haematite photoanodes

The oxygen evolution reaction (OER) plays a crucial role in (photo)electrochemical devices that use renewable energy to produce synthetic fuels. Recent measurements on semiconducting oxides have found a power law dependence of the OER rate on surface hole density, suggesting a multihole mechanism. In this study, using transient photocurrent measurements, density functional theory simulations and microkinetic modelling, we have uncovered the origin of this behaviour in haematite. We show here that the OER rate has a third-order dependence on the surface hole density. We propose a mechanism wherein the reaction proceeds by accumulating oxidizing equivalents through a sequence of one-electron oxidations of surface hydroxy groups. The key O–O bond formation step occurs by the dissociative chemisorption of a hydroxide ion involving three oxyl sites. At variance with the case of metallic oxides, the activation energy of this step is weakly dependent on the surface hole coverage, leading to the observed power law. The multihole mechanism of the oxygen evolution reaction on semiconductor electrodes has been hard to elucidate due to a lack of atomic-scale structural characterization of the material interface. Using pulse voltammetry and simulations of α-Fe2O3 photoanodes, this study predicts the chemical origin of the third-order rate dependence on holes.

The oxygen evolution reaction (OER) plays a crucial role in (photo) electrochemical devices that use renewable energy to produce synthetic fuels. Recent measurements on semiconducting oxides have found a power law dependence of the OER rate on surface hole density, suggesting a multihole mechanism. In this study, using transient photocurrent measurements, density functional theory simulations and microkinetic modelling, we have uncovered the origin of this behaviour in haematite. We show here that the OER rate has a third-order dependence on the surface hole density. We propose a mechanism wherein the reaction proceeds by accumulating oxidizing equivalents through a sequence of one-electron oxidations of surface hydroxy groups. The key O-O bond formation step occurs by the dissociative chemisorption of a hydroxide ion involving three oxyl sites. At variance with the case of metallic oxides, the activation energy of this step is weakly dependent on the surface hole coverage, leading to the observed power law.
The oxygen evolution reaction (OER) is the bottleneck in the (photo) electrochemical splitting of water into molecular hydrogen and oxygen, an enabling process for the production of fuels from sustainable energy sources. This four-electron half reaction (2H 2 O → O 2 + 4H + + 4e − , E° = 1.23 V versus the reversible hydrogen electrode (RHE)) is a key step in both natural photosynthesis and its synthetic counterpart. The mechanistic similarity of these two routes is an open question 1,2 . In plants, the OER is catalysed by the oxygen-evolving complex (OEC) in photosystem II, a Mn 4 CaO x cluster with a distorted cubane structure. Remarkable advances in the experimental characterization of the structure of the OEC have permitted an atomic-level understanding of the sequence of steps leading to the formation of oxygen. The OEC undergoes a sequence of four light-induced oxidations whereby the Mn centres store four oxidizing equivalents before the formation of the O-O bond, the rate-limiting step, takes place 3 . The energy cost of water oxidation is mitigated by storing the oxidizing equivalents at different Mn sites, while simultaneously releasing protons to prevent charge accumulation, before performing the four-electron chemistry to circumvent stepwise water oxidation 4 . Although the mechanism by which the O-O bond forms is not yet conclusively known, there are two models that have been extensively discussed in the literature: (1) nucleophilic attack of a water molecule or hydroxide ion on a Mn-oxo entity and (2) direct coupling between two oxo (or oxyl) groups, with ab initio simulations 5 and recent experiments 6 favouring the latter.
At variance with the natural process, our understanding of the synthetic OER route using heterogeneous catalysts is hampered by the difficulty in achieving an atomic-level characterization of the electrode/electrolyte interface under operation. Simulations, mostly based on density functional theory (DFT), have contributed remarkably to our knowledge of the OER mechanism. Most of our current understanding is based on a model consisting of a sequence of one-electron proton-coupled electron transfer (PCET) steps taking place at a metal Article https://doi.org/10.1038/s41929-022-00845-9 As the concentration of surface holes stored at oxyl groups increases approximately linearly with bias, the OER rate is exponential in both the applied bias and the concentration of surface holes, at variance with the power law behaviour for metal oxide photocatalysts.
In spite of the different dependence of the OER rate on surface holes (power law versus exponential), these recent studies stress that, both in metallic oxide electrocatalysts and in semiconducting oxide photocatalysts, the accumulation of holes at the surface of the catalyst is key, in analogy with the OER promoted in plants by the OEC 2 .
In this work we combined pulse voltammetry, DFT calculations of the full solid/liquid interface and microkinetic modelling to obtain an atomic-level understanding of the mechanism responsible for the multihole character of the OER on α-Fe 2 O 3 . We found that holes are accumulated at the solid/liquid interface through one-electron PCET steps leading to the formation of Fe-oxo groups with a character intermediate between Fe(III)-oxyl and Fe(IV)-oxo. In alkaline environments, O-O bond formation proceeds by the facile oxidation of solvent hydroxide ions to form a surface-bound superoxo intermediate. To operate, this elementary step requires the presence of three nearby Fe-oxo/oxyl groups. However, once these groups are present, the activation energy is weakly dependent on their concentration, resulting in a power law dependence of the OER rate on the surface hole concentration.

Experimental determination of the reaction order
To probe the dependence of the OER rate on the hole concentration, we prepared an α-Fe 2 O 3 photoanode by electrodeposition on a F:SnO 2 substrate followed by thermal annealing. In brief, the α-Fe 2 O 3 film had a nominal thickness of ~200 nm with an estimated particle size of 50-100 nm ( Supplementary Fig. 1). Raman and X-ray diffraction analysis ( Supplementary Fig. 2) showed a typical α-Fe 2 O 3 photoanode in excellent agreement with the literature. Further details of the preparation procedure and physicochemical characterization can be found in the Methods and Supplementary Note 1.
Although a number of techniques have been proposed and used to quantify the (surface) hole states [20][21][22] , we used the most direct method, that is, transient photocurrent measurements under chopped illumination (450 nm laser), integrating the cathodic current transients. The current-voltage curves of the photoanode ( Supplementary Fig. 3) in the dark and under static and chopped illumination agree well with the literature 18,23 . Also, the photocurrent increases linearly with light intensity ( Supplementary Fig. 4), as expected. The dependence of the OER rate on the hole concentration is summarized in Fig. 1a centre, whereby water is converted into hydroxy, oxo, hydroperoxo (and possibly superoxo) intermediates before O 2 is released, with the key O-O bond formation step taking place by the nucleophilic attack of water (or a hydroxide ion in an alkaline environment) on an oxo intermediate 1,[7][8][9] . Focusing on the free-energy cost of individual one-electron oxidations, this model has been enormously successful in explaining trends in metal and metal oxide electro-and photocatalysts, identifying the strength of the interactions of water oxidation intermediates with surface metal sites as the key descriptor of catalytic activity.
There is growing evidence, however, of systems where the OER proceeds by different mechanisms. Some oxides, including some perovskites, promote the OER through the oxidation of lattice oxygen [10][11][12] , while other oxides favour the direct coupling of two oxo ligands rather than the formation of a peroxo intermediate through nucleophilic attack 13 . A striking departure from the traditional sequence of PCET steps, however, is the recent discovery of the multihole character of the photoelectrochemical OER on a number of semiconducting oxides. The group of Durrant used transient absorption spectroscopy (TAS) on semiconducting metal oxide photoanodes, such as haematite (α-Fe 2 O 3 ), detecting a correlation between the spectroscopic fingerprints of electron holes accumulated near the solid/liquid interface and the OER rate. In particular, they found that the OER rate obeys a power law in the concentration of surface holes, where J ph is the measured photocurrent, k OER is the OER rate constant, [h + ] is the concentration of surface holes and α is the reaction order with respect to the hole concentration. Their experiments showed a first-order power law at low hole content and a third-order dependence once a sufficient concentration of holes was reached 14 . In another study, TiO 2 , BiVO 4 and WO 3 were found to display the same behaviour 15 . Together with other reports 16,17 , these measurements show that the multihole character of the OER is a general feature of metal oxide photocatalysts.
These recent findings suggest that there is a close analogy between the mechanisms of the photoelectrochemical OER on semiconducting oxides and on the OEC of photosystem II in plants, with hole accumulation preceding the slow O-O bond formation. At variance with the latter system, however, achieving an accurate structural characterization of the active site and detecting the reaction intermediates of the heterogeneous OER have proved daunting challenges, but substantial progress in this direction has been enabled by the development of in operando spectroscopic techniques. Zandi and Hamann 18 , in particular, detected the vibrational fingerprints of Fe(IV)-oxo in α-Fe 2 O 3 . Our recent work on IrO 2 has shown that the OER proceeds in this case by the nucleophilic attack of water on a surface oxyl, with an activation energy that decreases linearly with the concentration of surface holes 19 . and the whole range of the laser power were used to induce varying concentrations of hole charges, giving rise to OER currents that span approximately four orders of magnitude over the four potentials. The data are plotted as log-log plots of the OER current versus the integrated charge (hole number). Because the integrated charge represents the difference in the steady-state charges (light on or light off), the hole number includes surface holes as well as holes generated in the space charge layer (SCL) upon illumination. Thus, we measured surface and (partially) bulk holes. Although in the cathodic transient the steady-state hole concentration created under illumination is reduced to the steady-state level in the dark, the current transients do not count electron-hole recombination events occurring under steady-state illumination. The results suggest a near third-order dependence at all potentials, with some notable and statistically significant differences. At low potential, and thus at the lowest hole concentration, the reaction order is greater than 3. This is at odds with the first-order dependence proposed at low hole densities by Durrant and co-workers 14,15 . At 1.01 and 1.26 V, as well as at 1.51 V in the upper hole concentration range, the order is close to 3 (2.5-2.85). A decrease in the reaction order was observed at 1.51 V under low light intensity.
Not only can current transients be integrated, they can also be fitted to an equation to extract double-layer capacitance and charge (Supplementary Note 2 and Supplementary Fig. 5). The double-layer charge represents the portion of the total charge that shows an exponential dependence in the current transients and excludes diffusion-limited processes. Hence, we argue that the double-layer charge is a close approximation to the number of surface holes. In line with this assessment, the dependence of the OER rate on the surface hole density (see the derivation in Supplementary Note 3) follows closely that of the total charge. Now, however, the bending of the 1.51 V data is not observed, and analysis indicates that the bending is caused by the decreasing proportion of diffusion-limited charge in the low light intensity regime. The diminished influence of diffusion-limited charge in the double-layer charge reveals that the OER rate is nearly third order over the entire current range irrespective of the applied bias (Fig. 1b). It is worth noting that we recently used a similar methodology (pulse voltammetry) to determine the rate-determining step of the electrochemical OER on iridium oxides 19 . There, the OER rate showed an exponential dependence on the charge concentration. Here, the data do not show an exponential dependence but a third-order power law.

Modelling the spectroscopic fingerprints of surface holes
To understand the role of surface hole concentration on the OER mechanism over haematite, we turned to computational modelling, performed with DFT calculations. We focused on the (0001) facet, as this is the most frequently exposed surface in crystalline haematite minerals 24 . In contact with water, the surface is hydroxylated, and the two most favourable surface structures are the hydroxylated forms of the O and Fe terminations 25 . Previous work by Nguyen et al. 26 indicated that the O termination is catalytically more active, hence we focused on this structure. We further discuss the choice of this structural model for the haematite/water interface in Supplementary Note 4.
Holes photogenerated in the SCL are driven by band bending to the solid/liquid interface, where they can promote the oxidation of the surface hydroxy groups. To model this process we considered PCET steps leading to the dehydrogenation of the surface: *OH − + h + + OH − → *O 1− + H 2 O, where the * symbol indicates surface-bound species and h + indicates a valence band hole. The free-energy change associated with the PCET steps was evaluated using the computational hydrogen electrode method 27 . Solvation effects were accounted for by first-principles molecular dynamics (FPMD) of the explicit solid/liquid interface using the approach introduced by Wang et al. 28 . The position of the valence band maximum (VBM) on the electrochemical scale, E VBM = 1.91 V RHE , was estimated by the three-step approach introduced by Guo et al. 29 , and we assumed that the quasi-Fermi level of the holes at the solid/liquid interface coincides with the valence band edge of haematite (see Supplementary Note 6 for further details). The resulting surface Pourbaix diagram, shown in Fig. 2a, indicates that the dehydrogenation process covers a potential window of ~0.4 V, between 1.53 and 1.95 V. Under illumination, assuming the thermodynamic driving force for the PCET steps of dehydrogenation is given by the position of the VBM edge 30 , ignoring kinetic Article https://doi.org/10.1038/s41929-022-00845-9 effects and assuming these steps to be quasi-equilibrated, the surface should be almost completely dehydrogenated. We found that in the 1 to 2/3 monolayer (ML) hydrogen coverage range, the cost of dehydrogenation is approximately coverage-independent, whereas at lower hydrogen coverage, the cost increases, indicating that oxidizing the surface through PCET becomes increasingly difficult as the coverage of hydrogen is reduced (Supplementary Table 3). The removal of a proton and electron from the surface can formally lead to the creation of either Fe(IV)-O 2− or Fe(III)-O 1− surface species, depending on whether the hole localizes predominantly on the Fe or O atom. Our previous work 31 indicated that the localization of holes in bulk haematite depends strongly on the exchange and correlation functional adopted in the DFT calculations; hybrid functionals such as PBE0(α X ) (ref. 32 ), in particular, tend to favour localization at the oxygen site for values of α X larger than 0.30, while smaller values of α X lead to localization at the Fe site. To model the electronic structure of haematite, here we adopted the PBE0 functional with α X = 0.12. We obtained a fundamental gap of 2.47 eV and an optical gap of 2.07 eV, in good agreement with experimental estimates (2.6 ± 0.4 and 2.0-2.2 eV; refs. [33][34][35]. With this functional we found that the surface hole is shared equally between the Fe and O ions in an antibonding orbital, as shown in Fig. 2b. To simplify the notation, we will refer to this species as *O 1− bound in a bridging position between two Fe atoms. In the case of the removal of an electron without proton abstraction, the hole is localized at the Fe ion only, leading to a formal Fe(IV) subsurface species (Fig. 2c).
In operando TAS experiments showed a correlation between the presence of long-lived surface holes and the OER rate 36 . The spectroscopic fingerprint of these surface holes consists of a broad absorption band peaking at around 650 nm and extending into the infrared region to over 1,000 nm. The surface Pourbaix diagram shown in Fig. 2a indicates that holes reaching the solid/liquid interface promote the dehydrogenation PCET steps, leading to the formation of *O 1− . To explore whether these surface species are responsible for the TAS fingerprint, we performed time-dependent density functional theory (TDDFT) simulations with hybrid functionals (PBE0 with α X = 0.12) to compute the optical spectrum of the α-Fe 2 O 3 (0001)/water interface before and after a PCET step. The difference between the computed absorption spectra yields the spectroscopic fingerprint of a surface hole, trapped in an *O 1− intermediate. The spectrum shown in Fig. 2e displays a broad feature in the 600-1,000 nm region, with a maximum at 643 nm, in excellent agreement with TAS experiments 36 .
A similar simulated spectrum was obtained when we removed only an electron, leading to an Fe(IV) subsurface species (Fig. 2f). In this case, the maximum of the spectrum is at 626 nm, blueshifted by 17 nm with respect to the spectrum of *O 1− . Bulk holes also introduce a weak peak at around 650 nm in an otherwise featureless and very broad spectrum that extends from the absorption edge to about 2,500 nm (Fig. 2g). Our calculations therefore suggest that the formation of *O 1− surface species is thermodynamically favourable, and their optical spectrum matches closely the experimental feature at around 650 nm that correlates with the OER rate. Subsurface and SCL Fe(IV) holes, however, could also make a contribution to the absorption spectrum in the same energy range, if they survive recombination long enough to be probed by TAS measurements.

Mechanism of the OER
Having established that photogenerated holes driven to the α-Fe 2 O 3 (0001)/water interface can accumulate at the surface through the formation of *O 1− species, we next investigated the elementary OER steps, with the goal of understanding the effect of surface hole coverage. In the following, we will refer to the removal of a H atom as the creation of a surface hole.
We started by considering the conventional mechanism of the OER, consisting of a series of PCET steps, with O-O bond formation taking place via the formation of an *OOH − intermediate through the nucleophilic attack of a reactant molecule (a hydroxide ion in the case of alkaline conditions, a water molecule in acidic conditions) on a surface *O 1− intermediate 1,8,12,37 These four steps are PCET reactions and their free-energy cost, ΔG, was evaluated using the computational hydrogen electrode method 27 (9); bottom). The gold, red and white balls represent iron, oxygen and hydrogen atoms, respectively. The blue spheres represent the oxygen and hydrogen atoms that react with the surface oxygen atoms, and the green spheres represent the dehydrogenated surface oxygen atoms.
Article https://doi.org/10.1038/s41929-022-00845-9 ΔG of each step depends only on the adsorption energy of the various intermediates, not on whether the steps are considered in an acidic or alkaline environment 37 . The free-energy cost of each step is reported in Supplementary Fig. 11 as a function of the degree of surface oxidation, that is, the coverage of *O 1− . In agreement with previous studies 38 , we found that the potential limiting step is the first dehydrogenation (reaction (1)), that is, the formation of the *O 1− intermediate, which therefore determines the electrochemical potential that holes must have for all elementary steps to become exergonic. Importantly, when evaluated in a 4 × 4 unit cell, step (1) has a free-energy cost that is nearly constant up to 1/3 ML coverage of *O 1− . Moreover, pairs of *O 1− species interact attractively, with a nearest-neighbour interaction (−0.06 eV) that is more than twice the thermal energy available at room temperature, so even at low coverage pairing is expected (Supplementary Note 10). Note that we did not include the empty site * (that is, an oxygen vacancy) explicitly as a reaction intermediate as our FPMD simulations indicated that a nearby water molecule or hydroxide ion would be adsorbed spontaneously in just 2-3 ps. Moreover, in the presence of both an oxygen vacancy and *O 1− , our FPMD simulations indicated that water is adsorbed and dissociates spontaneously, suggesting that O 2 release can also take place according to the following elementary step: As it is thermodynamically favourable for *O 1− to be present in pairs, hydroperoxo can also be formed by a chemical reaction. This consists of the dissociation of a water molecule into *OOH − and *OH − adsorbates (reaction (6)). When three *O 1− are in close proximity, a new path becomes available, the chemical dissociation of a water molecule, leading to the formation of a superoxo *OO with formal charge 1− and two *OH − species (reaction (7)): Reactions (6) and (7) are both exergonic and, when compared at equal coverage of *O 1− , the formation of the superoxo is always thermodynamically favoured ( Supplementary Fig. 11).
We have therefore identified three mechanisms leading to the formation of the O-O bond: the standard PCET step, leading to a hydroperoxo through nucleophilic attack (reaction (2)), and the dissociation of a water molecule, leading to either a hydroperoxo (reaction (6)) or a superoxo (reaction (7)) intermediate. As the thermodynamics is favourable in all cases, we needed to examine their kinetics to determine which path dominates. To locate the transition state and compute the activation energy for the O-O bond formation step, we performed nudged elastic band (NEB) calculations in explicit solvent. The simulations were performed in a 2 × 2 charge neutral cell with the atomic coordinates taken from a snapshot extracted from an equilibrated FPMD simulation. The results of the NEB calculations are shown in Fig. 3.
To evaluate the kinetics of the PCET step leading to the hydroperoxo intermediate in an alkaline environment (reaction (2)), we included both a solvated OH − and Na + counterion in the explicit solvent. The hole was introduced by doping the haematite slab with a divalent Mg atom, substituting a trivalent Fe atom. The activation energy was calculated to be 0.43, 0.21 and 0.13 eV in the presence of one, two and three *O 1− , respectively.
The dissociation of water at two *O 1− sites leading to *OOH − and *OH − (reaction (6)) has a high activation energy of 1.40, 1.00 and 0.90 eV in the presence of two, three and four *O 1− in the 2 × 2 unit cell. When the dissociation of water leads to the formation of *OO − and two *OH − (reaction (7)), a path enabled by the presence of at least three *O 1− species, we found an almost identical transition state to the hydroperoxo case, involving the cleavage of an O-H bond, which also results in high activation energies (0.98 and 0.91 eV in the presence of three and four *O 1− , respectively). At sufficiently high surface hole coverage, the hydroperoxo intermediate is not stable and spontaneously donates the proton and an electron to the surface, leading to the more stable superoxo intermediate, in agreement with the thermodynamics reported in Supplementary Fig. 11. Moreover, we can see in Fig. 3 that the accumulation of surface holes has a large impact on the activation energy of O-O bond formation by water dissociation, but even on a fully oxidized surface the barrier remains substantial at 0.40 eV.
As O-O bond formation is generally acknowledged to be the rate-limiting step of the OER on haematite 16,39 , the barriers to the dissociation of a water molecule (~1 eV at coverages of around 3 h + nm −2 ) show poor agreement with the experimentally determined apparent activation energies on α-Fe 2 O 3 : 0.30 and 0.06 eV at low (<1 h + nm −2 ) and high (>1 h + nm −2 ) surface hole coverage, respectively, measured by Mesa et al. 15 , and 0.1-0.2 eV determined by Zhang et al. 40 . As discussed earlier, the kinetics of reaction (2), involving the hydroxide ion, is much more favourable. For this reason we also examined paths where the dissociative adsorption of OH − , rather than water, is responsible for O-O bond formation. Considering partial atomic charges, we found that the formation of the hydroperoxo (where both oxygen atoms are in the 1− oxidation state) is accompanied by the reduction of a surface *O 1− to *O 2− . Similarly, the formation of the superoxo is accompanied by the formation of both an *OH − and an *O 2− . We can therefore write these two steps as follows: The activation energies of these reactions are substantially lower than those involving water, and in better agreement with the experimental apparent activation energies 15,40 . We found that the barrier for the formation of the hydroperoxo is 0.23 and 0.16 eV in the presence of two and three surface holes, respectively. The barrier for the formation of the superoxo is 0.16 eV in the presence of three surface holes and decreases with increasing coverage, down to 0.02 eV on the fully oxidized surface. At a sufficiently high coverage of *O 1− , the hydroperoxo intermediate spontaneously evolves into a superoxo and a surface hydroxy group. As shown in Fig. 3, the dependence of the activation energy on the surface hole coverage is much weaker in reactions involving hydroxide ions than in those involving water.
These simulations show that O-O bond formation can take place at an isolated *O 1− site by the PCET step shown in reaction (2), leading to an *OOH − intermediate. As soon as pairs or trios of *O 1− sites are available, chemical reactions involving the (dissociative) adsorption of a hydroxide ion can also take place, with activation energies that are substantially lower than the single-site PCET mechanism. We also examined reaction (2) in further detail in the presence of more than one *O 1− site and the results are presented in Supplementary Fig. 13.
As we were performing static calculations without finitetemperature sampling of the degrees of freedom of the solvent, these estimates could be strongly affected by the choice of solvent structure (based on a single snapshot extracted from a FPMD simulation). Although this is common practice [41][42][43] , to validate our approach we performed a simulation based on thermodynamic integration; we used the O-O bond length as the reaction coordinate and sampled all the degrees of freedom of the system at 350 K by FPMD ( Supplementary  Fig. 16). We focused on the dissociation of a hydroxide ion in the presence of three *O 1− , leading to a superoxo intermediate. We obtained an activation energy of 0.19 eV and a reaction free energy of −0.73 eV, to be compared with 0.16 eV and −0.89 eV obtained by our static approach. This shows that our approach based on NEB calculations is reasonably accurate for the calculation of both activation and reaction energies. The difference in the results obtained by the two approaches (0.03 eV Article https://doi.org/10.1038/s41929-022-00845-9 for barriers and 0.16 eV for reaction free energies) can be taken as a rough estimate of the numerical accuracy of our approach.

Microkinetic model
To model the full OER we considered a simplified steady-state microkinetic model within the framework of transition-state theory. Our goal was to explore whether the elementary steps described in the previous section were sufficient to capture the qualitative features of the experimental measurements, namely a power law dependence of the rate on the coverage of surface holes and apparent activation energies in the order of 0.1-0.3 eV.
Our starting point was the microkinetic model of Wang et al. 28 , who investigated the photoelectrochemical OER on TiO 2 . A key feature of this model is the consideration among the reagents of a concentration of photogenerated valence band holes, C(h + ), in the SCL (SCL-holes) that reach the surface. In their model, Wang et al. investigated the OER rate as a function of C(h + ) using some experimental estimates of concentration. We adopted the same strategy, but used experimental data on the coverage of surface holes stored in *O 1− intermediates, Θ(*O 1− ), rather than the concentration of SCL-holes. Note that our simplified model does not describe the generation, recombination and diffusion of holes in the SCL, nor the spatial extent of the SCL, nor the gradient of hole concentration. As a consequence, the values of the SCL-hole concentration C(h + ) used in our model cannot be compared with an experimental estimate of this quantity, at variance with surface holes stored in *O 1− .
The reactions included in the model are listed in Table 1, and the reaction and activation energies are shown in Supplementary  Table 7. The first nine reactions (1)-(9) correspond to the elementary steps discussed in the previous section. In reaction (10), a surface *O 2− , generated in the chemical steps (8) and (9), is oxidized to *O 1− : The kinetic equations were solved under the conditions of mole fractions χ H 2 O = 1 and χ O 2 = 10 −7 at pH 8 and 300 K. The choice of pH was motivated by the fact that the experimental point of zero charge (PZC) is estimated to be at ~7.9 ± 1.7 (ref. 44 ). As our DFT simulations were performed on a neutral surface, our first-principles estimates of activation and reaction energies were obtained under PZC conditions.
Note, however, that we did not include in our microkinetic model the acid-base chemistry of the surface, which could describe charging at values of pH different from the PZC. The pH is considered only through the concentration of OH − in the expression of the rates of the elementary steps.
To describe the dependence of the activation energies of the O-O bond formation steps on the surface coverage of *O 1− , we linearly fitted the data obtained from NEB calculations. For the corresponding reaction free energies of steps (1)-(7), we linearly fitted the data displayed in Supplementary Fig. 11. The reaction free energies of steps (8) and (9) were obtained from the initial and final states of the NEB calculations, adding zero-point and finite-temperature corrections. We also explored the effect of using the difference between the initial and final states of the NEB calculations for steps (2), (6) and (7) as reaction energies, finding negligible differences (Supplementary Note 24). The reaction free energy of reaction (10) was obtained from a total energy calculation, as described in detail in Supplementary Note 19.
We define the turnover frequency as the rate at which O 2 is released (r OER ), that is, the sum of the rates of steps (4) and (5). Plotting r OER versus the coverage of *O 1− (Supplementary Fig. 18), we found that the data can be well fitted with a power law dependence: r OER = k OER [Θ(*O 1− )] α , with α = 3.20. This reaction order is consistent with the fact that, as shown in Supplementary Table 8, the O-O bond forms predominantly through step (9), which involves three neighbouring *O 1− sites. Small deviations from the power law behaviour can be seen at *O 1− coverages of around 0.1 ML. The activation energies of step (9) were evaluated in a 2 × 2 cell with *O 1− coverages exceeding 3/12 = 0.25 ML, which is considerably higher than both the maximum hole coverage obtained in our measurements (0.10 ML; Supplementary Note 3) and the maximum *O 1− coverage obtained using the microkinetic model (0.13 ML). To explore the behaviour of the activation energy of step (9) at low coverages, we exploited the reaction energies obtained in the 4 × 4 cell and the linear relationship between activation energy and reaction energy, as discussed in Supplementary Note 13. We found that the dependence of the activation energy of step (9) on Θ(*O 1− ) is lower at low *O 1− coverage than at high coverage (0.10 versus 0.22 eV per monolayer). The resulting microkinetic model yields a reaction order of 3.07, as shown in Fig. 4b. In the following we will discuss the results obtained with this refined model.   The reversibility is expressed as k b /k f , where k f is the forward rate and k b is the backward rate. The activation energy for hole diffusion, E h diff , was set as 0.29 eV, the experimental bulk value. We report results obtained at two SCL-hole concentrations, C(h + ) = 1 and 0.01 ML. Θ(*O 1− ) is the surface hole coverage, E app is the apparent activation energy, (C(h + )) is the reaction order with respect to SCL-holes and (Θ(*O 1− )) is the reaction order with respect to surface holes.
Article https://doi.org/10.1038/s41929-022-00845-9 To obtain the apparent activation energy, following Mesa et al. 15 , we fitted the dependence of k OER versus 1/T, obtaining E app = 0.18 eV at C(h + ) = 1 ML (Fig. 4c). To better interpret the outcome of the microkinetic model, we report in Table 1 the rates of the elementary steps, their reversibility and the degree of rate control (DRC), assuming the activation energy for hole diffusion, E h diff , to be the experimental value of 0.29 eV (ref. 45 ). We considered two SCL-hole concentrations, C(h + ) = 1 and 0.01 ML (we report the results of simulations with different values for C(h + ) and E h diff in Supplementary Notes 20-23). We can see in both cases that O-O bond formation takes place by step (9), involving the hydroxide ion, while steps (6) and (7), involving the water molecule, play a negligible role. Also step (2), the PCET reaction leading to the formation of the *OOH − intermediate, has a very low rate. The release of O 2 takes place by step (5), involving the dissociative adsorption of a water molecule.
At a high concentration of SCL-holes, the reaction with the highest DRC is step (1), suggesting the limiting factor is the slow diffusion of holes. At a low concentration of SCL-holes, the reaction with the highest DRC is step (9), suggesting that the limiting factor is the low coverage of *O 1− necessary for O-O bond formation. The involvement of the hydroxide ion in the steps with the highest DRC leads to both an OER rate and an OER rate constant with a power law dependence on the concentration of OH − , that is, they are exponential in pH ( Supplementary Fig. 20). We show in Supplementary Fig. 21 and Supplementary Note 22 that when step (10), which does not involve OH − , is the step with the highest DRC, the OER rate is essentially independent of the pH, giving a sublinear dependence. This can be achieved by increasing the barrier of step (10), for example, using the experimental estimate of Gardner et al. 46 of 0.46 eV for the hole diffusion barrier rather than the experimental estimate of 0.29 eV taken from Warnes et al. 45 . This choice of parameter represents a scenario in which hole diffusion through the SCL is slower than the surface reactions; this is a situation that could also be compatible with very efficient SCL-hole recombination, resulting in a high effective barrier for hole diffusion.
When we computed the reaction order with respect to the concentration of SCL-holes, C(h + ), rather than considering surface holes, we obtained a reaction order of approximately 1 when step (1) has the highest DRC and is therefore not quasi-equilibrated, which happens at high SCL-hole concentrations (Table 1). However, we obtained a reaction order of approximately 3 when step (1) is quasi-equilibrated, as Θ(*O 1− ) and C(h + ) are linearly related; this happens at low SCL-hole concentrations. The first order in SCL-holes is compatible with experiment as the OER rate is first order with respect to the light intensity ( Supplementary Fig. 4).
Overall we can conclude that this elementary microkinetic model is able to reproduce the main features of the OER on haematite seen experimentally, namely a power law dependence of r OER on Θ(*O 1− ) with reaction order 3, a low apparent activation energy, an OER rate that increases with pH and an approximately linear increase in rate with light intensity.
In Supplementary Note 21 we also discuss the results of a simplified model where we neglected the dependence of the activation and reaction energies on the coverage of *O 1− and arbitrarily enabled or blocked selected reactions to investigate the effects on the reaction order. Including only the four PCET steps (reactions (1)-(4)) in the model, that is, considering the traditional single-site model normally adopted in the simulation of the OER on metal and oxide surfaces 7,8 , we obtained an OER rate that is second order in Θ(*O 1− ). This is because step (1) is quasi-equilibrated, hence Θ(*O 1− ) and C(h + ) are linearly related.
Step (2) has a DRC of 0.99 and its rate depends on the product of Θ(*O 1− ) and C(h + ), hence the second-order behaviour. Including all ten steps but blocking step (8), we obtained a reaction order of 2.95, while blocking step (9) gave a reaction order of 2.00, showing that when these two reactions are enabled, the reaction order depends on their relative rates.
The possible mechanisms emerging from our microkinetic model are summarized in Fig. 5. Figure 5a illustrates the single-site mechanism consisting of a sequence of four PCET steps, with the O-O bond-forming step occurring by nucleophilic attack. Figure 5b,c illustrates multisite mechanisms, where the PCET step in reaction (1) leads to the accumulation of surface holes before O-O bond formation, which takes place by a chemical reaction. Using our computed values of activation and reaction energies, the microkinetic model supports a third-order mechanism in which step (9), the dissociative adsorption of a hydroxide ion on three neighbouring *O 1− , leads to O-O bond formation.

Discussion
High-valent Fe-oxo species are well-known intermediates in a number of important reactions: haem and non-haem Fe(IV)-oxo complexes activate dioxygen to perform a wide range of oxidation and C-H activation reactions 47 , Fe(IV)-oxo is the key intermediate in the Fenton reaction at around neutral pH (ref. 48 ) and is also formed as an intermediate in the OERs performed in the presence of a number of Fe-based molecular catalysts [49][50][51] . Ligands play an important role in determining the oxidizing power of mononuclear Fe(IV)-oxo complexes 47 . The implications of the precise location of the hole in non-haem Fe(IV)-oxo complexes, namely whether the intermediate should be described as Fe(IV)-oxo or Fe(III)-oxyl, have been recently discussed by Ye and Neese 52 in the context of C-H bond activation, and have also previously been investigated by Decker et al. 53 . While these studies focused on mononuclear Fe-oxo intermediates, binuclear complexes have also been investigated. In these systems, terminal oxygen atoms, Fe(IV)=O, were found to be more reactive than bridging oxygen atoms, Fe(IV)-O-Fe(IV), in both the oxygen transfer reaction and the hydrogen abstraction reaction 54 .
In the context of heterogeneous OER, Zandi and Hamann detected Fe(IV)-oxo in haematite using operando vibrational spectroscopy 18 , and this intermediate was assumed to be formed by the oxidation of Fe(III)-OH in a number of mechanistic studies of the OER on haematite 15,38 . Our model, where the oxidation of Fe(III)-OH to the (formal) Fe(III)-oxyl intermediate is the key step in the accumulation of holes at the surface, is therefore in line with existing knowledge of this system.
The key prediction of our model is the power law dependence of the OER rate on the concentration of surface holes, which is fully consistent with our experimental measurements as well as with earlier experimental findings. This behaviour contrasts with results reported for metallic oxides such as IrO x , where the OER rate depends exponentially on the concentration of surface holes 19 . Our calculations have revealed that the origin of this difference is the dependence of the activation energy for the formation of the O-O bond on the concentration of surface holes: in both materials, the activation energy decreases linearly with increasing hole coverage, but the slope in the case of haematite is extremely small (Fig. 3a), to the extent that the deviation from linearity in the log-log plot in Fig. 4b is not detectable. As activation and reaction energies are approximately linearly related (Brønsted-Evans-Polanyi (BEP) principle, see Supplementary Fig. 13b), the problem can equivalently be recast into a different behaviour of the activation energy for O-O bond formation versus the reaction energy. As shown in Supplementary  Fig. 13b, for the dominant step (9) involving the dissociative adsorption of OH − , we found an activation energy that is weakly dependent on the reaction energy (BEP slope <0.1), a signature of an early transition state. The O-O bond formation step in the case of IrO x involves the dissociative adsorption of a water molecule, which displays a large BEP slope of around 0.9 (ref. 19 ), a signature of a late transition state.
The experimental literature 15,17,40,55,56 consistently reports photocurrents, at fixed potential on the RHE scale, that increase with pH. While the results of Iandolo and Hellman 55 and of Liu et al. 56 suggest an approximately linear increase in photocurrent with pH, Zhang et al. 40 reported a sharp increase above pH 10, followed by a slower increase above pH 13. The trends described in the other reports are less obvious, but nonetheless consistent with increasing OER rates under more alkaline conditions. While Zhang et al. attributed this behaviour to the direct involvement of OH − in alkaline conditions, the other authors reported a mitigation of Fermi level pinning under increasingly alkaline conditions. Our results support the interpretation of Zhang et al. 40 . We stress, however, that the possible involvement of surface states in Fermi level pinning and the mitigation of such an effect with increasing pH cannot be described by our model. Note that the surface states detected by Liu et al. 56 through impedance measurements inversely correlate with the OER rate, at variance with the holes detected by both optical spectroscopy and transient photocurrent measurements, suggesting that the different experimental techniques probe different species.
When the step with the highest DRC involves a hydroxide, our model predicts an OER rate that has a power law dependence on the concentration of OH − , that is, an exponential dependence on pH. As shown in Supplementary Fig. 21, this behaviour is critically affected by the relative rates of the elementary steps considered. It is possible that in experiments displaying a linear dependence of the photocurrent on pH, the step with the highest DRC does not involve OH − . This could be due to a number of reasons. In conditions more alkaline than the PZC, the surface is partially deprotonated, hence holes will reach the surface by step (10) rather than by step (1). Another possibility is that, regardless of the pH, hole migration through the SCL might be the step  with the highest DRC. This would be consistent with a reaction such as step (10), namely a step independent of the concentration of OH − , being the step with the highest DRC. To accurately model the effects of pH one should consider both protonation and deprotonation of surface oxygen sites, study how the activation and free energy of the elementary steps depend on pH, and include the diffusion of holes through the SCL, features that are not present in our model. Our calculations have revealed that the origin of the approximate third-order behaviour is the favourable kinetics of step (9), the direct formation of the superoxo intermediate upon adsorption of a hydroxide ion. This step requires the presence of a surface oxyl active site as well as two nearby surface holes to accommodate the release of two electrons from the hydroxide ion, hence the third-order behaviour. The sequence of PCET steps at neighbouring Fe(III)-OH sites (reaction (1)), or electron transfer steps in highly alkaline conditions (reaction (10)), are the steps that enable this process, leading to the accumulation of the three surface holes necessary for the adsorption, dissociation and oxidation of the incoming hydroxide ion.
Our findings suggest an analogy with the OEC promoting the OER in biological systems, where holes are stored at different metal-oxo sites in the PCET steps that precede O-O bond formation. The key notion of redox potential levelling 57 in nature's OEC, enabled by the PCET steps, is mirrored in haematite by the near independence of the free-energy cost of step (1) on the surface hole concentration.
We monitored the dependence of the OER current on the charge stored at the interface between haematite photoanodes and the electrolyte through transient photocurrent measurements. We found that the OER rate exhibits an approximate third-order power law dependence on the surface hole density over four orders of magnitude of photocurrent. This is in agreement with previous measurements based on time-resolved optical spectroscopy, which detected third-order behaviour at high surface hole coverage; the same measurements detected first-order behaviour at low coverage, which was not observed in our measurements.
Our calculations have revealed that holes, upon reaching the surface, can promote a PCET step, oxidizing an *OH − group. Holes form *O 1− surface species that have a strong absorption peak at around 650 nm, in line with experimental evidence. Subsurface Fe(IV) holes, as well as holes in the SCL, however, can also contribute to the absorption spectrum in the same energy range.
To understand the origin of the power law dependence of the OER rate on the surface hole concentration, we considered several possible mechanisms involving one, two and three *O 1− intermediates. We monitored the dependence of the activation and reaction energy on the surface hole concentration and built a self-consistent microkinetic model to account for such dependence. We found that the reaction proceeds through a sequence of one-electron PCET steps that stores three holes at the surface, enabling a chemical reaction in which a hydroxide ion is adsorbed dissociatively, leading to the formation of a superoxo intermediate. This model is consistent with the experimental data, displaying an approximate third-order OER rate in the concentration of surface holes, a low apparent activation energy (0.18 eV), a rate that increases approximately linearly with light intensity and that increases at more alkaline conditions.
Our work shows that the multihole character of the OER on haematite stems from a multisite mechanism that marks a clear departure from previous models. We have therefore established a close analogy between the accumulation of oxidizing equivalents in the OEC in photosystem II that precedes the rate-determining step of the OER 3 and the multihole character of the OER on metal oxide photocatalysts 15 .

Sample preparation
The α-Fe 2 O 3 photoanode was prepared by a two-step wet thermochemical route 58 . Anodic electrodeposition was carried out in a plating solution of 16 mM (NH 4 ) 2 Fe(SO 4 ) 2 , 2.4 M (NH 4 ) 2 SO 4 and 2 mM H 2 SO 4 in ultrapure water. The solution pH was adjusted to 7.5 by addition of KOH. A F:SnO 2 coated glass substrate purchased from Solaronix was cleaned in 2-propanol and ultrapure water and connected in the plating solution as the working electrode in a three-electrode configuration with a Ag/AgCl electrode and a Pt wire serving as the reference and counter electrode, respectively. The substrate was biased at a potential of -0.1 V versus Ag/AgCl and a total charge of 10 μAh cm −2 was deposited. Afterwards, the as-deposited films were cleaned in ultrapure water, dried in a stream of N 2 and annealed for 20 min at 800 °C in air.

Photoelectrochemical measurements
Photoelectrochemical measurements were carried out in a home-built polytetrafluoroethylene cell 30 . Analytical grade 1 M KOH solution purchased from Roth was used as electrolyte. The pH of the electrolyte was 13.8. A Ag/AgCl electrode and a Pt wire served as reference and counter electrode, respectively. Sample illumination was provided by a home-built set-up composed of Thorlabs components: an LP450-SF15 450 nm monochromatic optical fibre laser diode in combination with a TTC001 temperature controller, an LDM9LP diode mount and an SH05 optical shutter. The timing and control of the experiments were achieved with a BioLogic SP200 bipotentiostat in combination with EC-Lab software. The electrode in contact with the electrolyte and the aperture of the light source both had an area of 0.196 cm 2 . Transient photocurrent measurements under chopped illumination were conducted at 0.76, 1.01, 1.26 and 1.51 V RHE . The whole range of laser power, from the lowest settings at which the diode emits light to the maximum power of the laser, was used to vary the incident light intensity.

Characterization methods
UV-visible optical absorption spectra were recorded with an Agilent Cary 5000 UV-Vis-NIR spectrometer. The recordings were executed in transmission mode, with a bare F:SnO 2 substrate serving as baseline. No notable reflective losses were found on either side of the sample. Raman spectra were recorded using a Raman microscope assembled by S&I Spectroscopy & Imaging. A TriVista TR557 spectrometer provided spectral dispersion and a Princeton Instruments PyLoN:2kBUV UV-enhanced CCD (charge-coupled device) camera was used for signal detection. The fundamental wavelength was 633 nm at a focused continuous-wave incident power of 1.24 mW. X-ray diffraction measurements were performed in Bragg-Brentano geometry on a Bruker AXS D8 Advance II Theta/Theta diffractometer using Ni-filtered Cu Kα radiation and a position-sensitive energy-dispersive LynxEye silicon strip detector. Samples were positioned on a larger, supporting glass slide inside a height-adjustable sample holder. A focused ion beam (FIB) lamella was prepared in a Thermo Scientific FIB/SEM DualBeam Helios NanoLab G3 UC device. Transmission electron microscopy images were measured with a Thermo Scientific Talos F200X microscope equipped with an X-FEG field emission gun.

DFT calculations
DFT calculations were performed with the Quantum ESPRESSO code 59-61 using Optimized Norm-Conserving Vanderbilt (ONCV) pseudopotentials 62 with a planewave cut-off of 80 Ry. We used the revised van Voorhis (rVV10) exchange and correlation functional 63 with the b parameter set to 9.3 (ref. 64 ), using a Hubbard correction for the Fe ion of 4.3 eV. The α-Fe 2 O 3 (0001) surface was modelled using 2 × 2 and 4 × 4 surface unit cells with slab geometries including 6 Fe atomic layers separated by at least 10 Å of vacuum. The Brillouin zone was sampled at the gamma point. Minimum energy paths and transition states were identified through climbing-image NEB calculations 65 using 10-12 images. Vibrational frequencies were computed from the atomic displacements (0.01 Å) of the reaction intermediates. Goedecker-Teter-Hutter pseudopotentials 67 , a plane-wave cut-off of 500 Ry, the rVV10 exchange and correlation functional with the b parameter set to 9.3, and a Hubbard correction for the Fe ion of 3.3 eV (ref. 68 ). We used deuterium masses for H atoms, a time step of 1 fs and a temperature of 350 K, imposed by the canonical sampling velocity rescaling thermostat 69 . The simulations were performed on a 2 × 2 unit cell with hydroxylated termination, modelling the aqueous solvent with 93 water molecules. The starting geometry was taken from the work of Blumberger and co-workers 70 . The FPMD simulations to obtain the geometries for the optical spectra were performed using the CP2K code with the PBE0 functional with α = 0.12 and a truncation cut-off of 4 Å for the exchange operator 71 and the cpFIT3 Auxiliary Density Matrix Methods (ADMM) 72 and MOLOPT-DZVP-SR basis sets. To reduce the computational cost of these hybrid functional FPMD calculations, the solvent was modelled with 47 water molecules. The other parameters were identical to the DFT+U simulations described above.

Simulations of the optical spectra
Optical spectra were simulated by time-dependent density functional perturbation theory (TDDFPT) using the hybrid functional implementation of TDDFPT in the CP2K code and computing 300 excited states for bulk haematite, 600 excited states for the pristine surface and 800 excited states in the presence of holes. The spectra of surfaces were obtained by time-averaging spectra of individual snapshots extracted every 500 fs from ~15 ps FPMD, applying a Gaussian broadening of 0.05 eV. The spectra of bulk haematite were simulated using a 3 × 3 × 3 unit cell and snapshots extracted every 100 fs from an ~3 ps FPMD. The optical gap was estimated from a Tauc plot of time-averaged absorption spectra. All the parameters used in these calculations were identical to those used in the FPMD simulations with the PBE0 functional.

Microkinetic model
The Matlab code for solving the microkinetic model was written by J.-F. Chen and is available at the following public repository: https://github. com/JianFuChen/Kinetic-solver. The code was used by Wang et al. 28 to study the OER on TiO 2 . We modified the code to include our ten elementary steps and the dependency of the activation and reaction energies on the hole coverage, which required the rate equations to be solved self-consistently. Further details are given in Supplementary Note 18.

Data availability
The data used to generate Figs. 1-4 are available from the Open Research Data Repository of the Max Planck Society, https://doi.org/10.17617/3. XYOGI0 (ref. 73 ). The data presented in the Supplementary Information are available from the corresponding authors upon reasonable request.