Introduction

Photoelectrochemical (PEC) water splitting is one of the most promising ways to convert solar energy into clean hydrogen energy1,2,3,4. Compared to the two-electron hydrogen evolution reaction (HER) at the cathode, the four-electron oxygen evolution reaction (OER) at the anode imposes a larger limitation on the efficiency of PEC water splitting due to its sluggish kinetics5,6,7. Therefore, there is a significant amount of theoretical and experimental research focusing on exploring the OER mechanisms and designing efficient OER catalysts. In addition to high catalytic activity, an ideal OER catalyst should be non-toxic, cost-effective, chemically and mechanically stable, highly conductive, and easy to integrate into the anode8,9. Hematite (α-Fe2O3, Fe2O3 hereafter) is a highly competitive candidate that satisfies all these requirements10.

However, Fe2O3 as a photoanode material still faces many challenges, such as narrow visible light absorption range, high photo-generated electron-hole recombination rate, low photo-generated electron-hole lifetimes, and particularly, large OER overpotentials11,12. To tackle these challenges, large amount of theoretical and experimental research has been conducted, as summarized in several review publications13,14,15,16,17,18,19, and indeed, significant progress has been made11,20,21. However, there is still a large gap existing between the record efficiency and the theoretical limit22. To improve further the performance of PEC water splitting on Fe2O3, it is important to explore the nature of Fe2O3-water interfaces and to propose a reaction mechanism for water oxidation23. Reaction intermediates are crucial for proposing a reaction mechanism, but they are difficult to identify due to the high interfacial complexity under the reaction conditions. Experiments with various in situ probing techniques have been used to resolve the challenge. Until now, however, few intermediates have been determined and only a rough sketch was drawn24,25,26,27,28,29,30. Theoretical simulations can offer details of the reaction mechanism at the atomistic level, and some intermediates were expected to be essential in the process of water (photo)oxidation, for instance, hydroxyl group (*OH), ferryl group (Fe=O), and peroxo (O-O) species involving lattice Os31. However, at present most of these simulations explore the structure of Fe2O3-water interfaces32,33,34,35,36,37 or adopt the classic four-step PCET mechanism to study water (photo)oxidation, which limits the observation of various intermediates 38,39,40,41,42,43,44.

As an example, the bottleneck intermediate of PEC water oxidation on Fe2O3 photoanodes is still under debate among various experimental observations and theoretical simulations. The operando infrared spectroscopy experiment demonstrated that the bottleneck intermediate in the PEC water oxidation was assigned to the ferryl group24, while many ab-initio simulations supported a viewpoint on the bottleneck intermediate of the hydroxyl group38,39,40,43,44. There are also some simulation studies which offered the viewpoint of the bottleneck intermediate consistent with the experiment but through diverse ways, such as comparing simulated and experimental optical absorbance spectrum or signature45,46, employing elaborate simulation models involving iron or oxygen vacancy defects41,47. Since the simplified four-step PCET mechanism was adopted in these simulations, more details of the reaction mechanism are missing. Experimental studies have already demonstrated that hydrogen bonds play an important role in the process of water photo-oxidation25,48. The complex hydrogen bonding network presents also a challenge in defining the structure of Fe2O3-water interface 35,49.

In this work, as shown in Fig. 1, we propose a comprehensive reaction mechanism for OER on Fe-terminated Fe2O3 (0001) surface in the aid of ab-initio molecular dynamics (AIMD) simulations, a computational method based on quantum mechanics without recourse to empirical parameters to describe the inter-atomic interactions50. The AIMD simulation suggests a potential capability to capture the structure of the hydrogen bonding network and to explore possible intermediates of the water oxidation reaction. Also, some intermediates which were adopted in previous study31 appear to have more stable alternatives and excluded from the reaction mechanism, for example merely the *OOH group on Fe-terminated Fe2O3 (0001) surface. Through the reaction mechanism, we assigned the key OER intermediates on the non-solvated and solvated surfaces, which are consistent with the experimental observations24,25,51. We further proposed that the O2 desorption may be a key step to affect the OER dynamics on Fe2O3 surfaces, due to the improvements in exposure of the surface reaction sites and contact of water solvent, which thereby alters the reaction pathway. Our work may help understand the nature of water oxidation on Fe2O3 and design enhanced Fe2O3 photoanode materials for PEC water splitting.

Fig. 1: Schematic diagram of the OER mechanisms on Fe-terminated Fe2O3 (0001) surface.
figure 1

The figure illustrates the intermediates involved in the classic 4-step PECT mechanism, the 10-step non-solvated surface mechanism, and the solvated surface mechanism.

Full size image

Results

OER on Non-solvated Fe-terminated Fe2O3 (0001) Surface

The initial step of OER on the non-solvated Fe-terminated (0001) surface involves adsorption and dissociation of the 1st water molecules, namely equation A and B in Section 3.2. In literature32,33,34, two water adsorption configurations on this surface were proposed, i.e., molecular and dissociative water adsorption. Fig. 2a, b shows the processes of water dissociation on the surface at 300 K with the initial configurations of molecular and dissociative water, respectively. Spontaneous water dissociation was observed and the dissociative adsorption is more energetically favorable than the molecular adsorption with an energy gain of about 0.18 eV. This result is supportive of the calculated small activation energies for water dissociation on this surface33,34,52, but does not fully agree with the observation of polarization-dependent infrared reflection absorption spectroscopy, where the molecular water adsorption existed and was identified as minor species53. Analyzes of the O-H distances on this surface provide insights into the dynamic behaviors of chemical bonding during water dissociation. The details were shown in Supplementary Fig. 1 and Note 1. The HOMO charge density of the dissociative water adsorption was presented in Fig. 3a. The charge density is spread over the slab instead of the surface hydroxyl group, indicating that the lattice Os, rather than the adsorbed O, might be the active site for the following 1st PCET step.

Fig. 2: Intermediates for OER on non-solvated Fe-terminated Fe2O3 (0001) surface determined with AIMD simulations.
figure 2

a and b Adsorption and dissociation of the 1st water molecule, namely equations A and B in Section 3.2; c and d Formation of a hydroxyl group after the 1st PCET step, namely equation C in Section 3.2; e and f Formation of a ferryl and peroxo-like group after 2nd PCET step, namely equations D and E in Section 3.2; g, h, and i Adsorption and restructuring of the 2nd water molecule, namely equations F and G in Section 3.2; j Formation of a superoxo-like species after the 3rd PCET step, namely equation H in Section 3.2; k Formation of an adsorbed O2 molecule after the 4th PCET step, namely equation I in Section 3.2. Source data are provided as a Source Data file. The optimized structures are provided in Supplementary Data 1.

Full size image
Fig. 3: HOMO charge density distribution on the non-solvated surface.
figure 3

a Dissociative water adsorption; b Fe-hydroxylated surface; c *OOH species stabilized by ΔH; d Superoxo-like species. HOMO charge density distribution on solvated surface. e Dissociative water adsorption with adsorption of another water molecule; f Hydroxylated surface with adsorption of another water molecule; g Bi-hydroxyl group. The charge density was visualized with Vesta software82. Minimum energy pathways for (h) formation of the peroxo-like species, i restructuring of the 2nd water molecule to form the *OOH species, and (j) restructuring of the bi-hydroxyl group to form *OOH species. Source data are provided as a Source Data file. The structures of initial, transition, and final states are provided in Supplementary Data 1.

Full size image

The 1st PCET step involves removal of one H from the stable configuration of dissociative water adsorption, namely equation C in Section 3.2. Fig. 2c, d show the formation of a hydroxyl group (*OH intermediate) adsorbed at the surface Fe. The intact hydroxyl group has an energy gain of about 0.7 eV over the dissociative hydroxyl, demonstrating that the stable product of the 1st PCET step is the Fe-hydroxylated surface. This result is consistent with previous theoretical studies41,44 and the experiment where the active site accumulates a large amount of the hydroxyl group before the onset of OER28. Supplementary Fig. 2 shows the analyzes of the O-H distances on this surface, which supports the stable Fe-hydroxylated surface. The HOMO charge density of this surface was shown in Fig. 3b. The HOMO is mainly localized on the O of the hydroxyl group rather than the surface lattice Os, suggesting the capability of the hydroxyl group to capture holes during water oxidation reaction.

The 2nd PCET step leads to formation of a ferryl group (*O in equation D in Section 3.2) or a peroxo-like species (*OΔ in equation E in Section 3.2), which have been reported in previous simulations with the former less stable than the latter31,54. Taking the ferryl group as the initial configuration, our AIMD simulation gives rise to a tilted *O intermediate adsorbing at the surface Fe, see Fig. 2e. The tilted ferryl group is more stable than the vertical one, with the Fe–O bond length shortened by 0.045 Å and an energy gain of about 0.3 eV. This phenomenon may be explained by the tetrahedral crystal field of surface Fe atom coordinated with four O atoms, in which the slight tilt of Fe–O bond alleviates the repulsion interaction between Fe and O orbitals. Starting from the latter configuration, as shown in Fig. 2f, the peroxo-like species remains stable on the surface during the simulation, which is slightly less stable than the tilted ferryl group by 0.05 eV. The analyzes of the O–O distances in the two configurations were shown in Supplementary Fig. 3, and the corresponding dynamic behavior of chemical bonding was described in Supplementary Note 2.

No transition was observed between the two configurations in our AIMD simulations, suggesting that there might be a notable activation energy for such a transition. Fig. 3h shows the minimum energy pathway for the transition between the two configurations. The activation energy for the transition state was 0.83 eV with only one imaginary frequency of 463 cm−1. Gebauer et al. also reported an activation energy of 0.86 eV for this transition state, and assigned the origin of the high energy barrier to the different electronic and spin structures of the two configurations54.

The following step is adsorption and restructuring of the 2nd water molecule on the surface, namely equations F and G in Section 3.2. Molecular adsorption of the 2nd water molecule can be stabilized during the period of simulation, see Fig. 2g. Restructuring of the 2nd water molecule gives rise to a *OOH species binding to the surface Fe and a H bonding with the surface lattice Os (*OOH + ΔH), as shown in Fig. 2h, or a bi-hydroxyl group at the surface Fe, as shown in Fig. 2i. Among the three configurations, the most stable is the molecular water adsorption, then is the *OOH species with ΔH, and the least stable is the bi-hydroxyl. The last one has been reported previously, but therein the bi-hydroxyl configuration is more stable than the molecular adsorption32,33,34. The reason may be assigned to the absence of the O–Os bond in the configuration of molecular water adsorption therein. The dynamic chemical bonding in adsorption and restructuring of the 2nd water molecule was shown in Supplementary Fig. 4 and analyzed in Supplementary Note 3.

The *OOH species is expected to be important for OER, which can be formed by restructuring the absorbed 2nd water molecule or the bi-hydroxyl group. The HOMO charge density for the *OOH configuration was shown in Fig. 3c, which is mainly localized on the *OOH, suggesting that the *OOH species may be oxidized in the following steps. In the former way, the *OOH species can be formed in such a way that one H of the 2nd water molecule binds with the surface lattice Os followed by breaking of the O–Os bond and formation of the *OOH species. The minimum energy pathway for such a transformation was shown in Fig. 3i and the transition state has an activation energy of 1.04 eV with only one imaginary frequency of 874 cm−1. However, in the latter way, restructuring of the bi-hydroxyl group shows a smaller activation energy of 0.63 eV for the transition state, see Fig. 3j, with only one imaginary frequency of 918 cm−1. Furthermore, the HOMO charge density of the bi-hydroxyl configuration in Fig. 3g was localized on the Os of the two hydroxyl groups, suggesting a possibility that the bi-hydroxyl configuration experiences directly the 3rd PCET step.

The 3rd PECT step involves removal of one H from the *OOH species, namely equation H in Section 3.2. Previous theoretical studies have reported that the *OOH species is the product of the 3rd PECT step31,40. However, our AIMD simulation suggests that only the *OOH species on the non-solvated surface can not be stabilized at room temperature, as shown in Fig. 2j. Starting with the *OOH intermediate without ΔH, there formed simultaneously a superoxo-like species with a magnetic moment of ~1.1 μB adsorbed at the surface Fe and a H binding to the surface lattice Os (ΔH), giving an energy gain of 0.68 eV. Meanwhile, a hydrogen bond was formed between the superoxo-like species and the H, which contributes to the energy gain. Supplementary Fig. 5a shows that the H of *OOH binds to the surface lattice Os within 0.5 ps, giving rise to a O–H bond with the average length around 0.97 Å and a O–O bond with the average length around 1.3 Å. There are two other hypothsized configurations for the product of the 3rd PCET step with each having a O–Os bond. However, they are energetically unfavorable, as shown in Supplementary Fig. 5b, d. Based on this observation, we inferred that the 3rd PCET step took place via oxidation of the H of *OOH rather than the ΔH, as expected from the HOMO charge density distribution shown in Fig. 3c. The HOMO charge density of the superoxo configuration was shown in Fig. 3d, which is mainly localized on the superoxo-like species and the surface Fe.

Then, the 4th PCET step removes the H bonded with the surface lattice Os in the superoxo configuration, namely equation I in Section 3.2, which is followed by desorption of the O2 molecule (equation J in Section 3.2). Figure 2k shows that after removing the H, the superoxo configuration can not be stabilized any longer and the configuration of an O2 molecule adsorbed at the surface Fe was formed simultaneously. Supplementary Fig. 6a shows that the O–O bond length in the adsorbed O2 molecule is around 1.25 Å. The adsorbed O2 configuration was also considered as the initial state for another AIMD simulation, as shown in Supplementary Fig. 6b, which is stable during the period of simulation.

Based on the reaction intermediates identified with the AIMD simulations, Fig. 4 depicts the reaction mechanism for OER on the non-solvated Fe-terminated Fe2O3 (0001) surface. The reaction pathway involves ten reaction steps (equation A to J), which are made up of four electrochemical steps and six non-electrochemical steps. All the electrochemical steps without external electric biases, equations C, D, H, and I in Fig. 4a, are uphill processes in Gibbs free energy. Equation D, dehydrogenation of the hydroxyl group adsorbed at the surface Fe, has the highest thermodynamic energy barrier of 2.10 eV and thus an overpotential of 0.87 eV, which is consistent with previous calculations done by Hellman et al38. and Liao et al39. Among the six non-electrochemical steps, equations A, B, and F are exothermic processes and equations E and J are endothermic processes. Equation G depends on the configuration of its reactant, namely an endothermic process for molecular adsorption of the 2nd water molecule or an exothermic process for the bi-hydroxyl configuration, as illustrated by the dashed lines in the middle part of Fig. 4. The endothermic steps are of particular importance to the overall reaction. Since the electrochemical steps can be modulated by external electric biases, the non-electrochemical steps can only be thermally activated. The non-electrochemical endothermic steps thus constitute the rate-limiting step under PEC operation conditions, as shown in the free energy diagram for OER with an external electric bias of 2.10 V in Fig. 4b. The last two steps in the reaction pathway could also be exchanged, namely desorption of the O2 molecule before the 4th PCET step, as illustrated by the dashed line in the right part of Fig. 4. In this case, however, the desorption of the O2 molecule has a larger thermodynamic energy barrier of 0.59 eV.

Fig. 4: Reaction free energy diagrams for OER on non-solvated Fe-terminated Fe2O3 (0001) surface with the error bars marked, i.e., the standard deviation of the DFT energy for a set of the snapshots sampled from the simulation trajectory of each intermediate.
figure 4

a Without electric bias. b With an electric bias of 2.10 V. Source data are provided as a Source Data file.

Full size image

The non-electrochemical step equation E, formation of the O–Os bond between the adsorbed O and surface lattice Os, has a very small thermodynamic energy barrier of 0.05 eV but a considerable activation energy of 0.83 eV. However, formation of the bi-hydroxyl group from the reaction of the tilted ferryl group with the 2nd water molecule, namely from equation D directly to F, is exothermic without noticeable activation energy. Equation G of restructuring the 2nd water molecule to the *OOH species has a small thermodynamic energy barrier of 0.15 eV, but a larger activation energy of 1.04 eV. In contrast, equation G of restructuring the bi-hydroxyl group to the *OOH species is an exothermic process with a smaller activation energy of 0.63 eV. The adsorbed 2nd water molecule is energetically more stable than the bi-hydroxyl group, but kinetically unfavorable to form the *OOH species. As a result, the reaction pathway of D → F (the bi-hydroxyl group) → G is preferred to that of D → E → F (the absorbed 2nd water molecule) → G in case that the water molecules can readily approach the surface. Otherwise, the tilted ferryl group and the adsorbed 2nd water molecule can be deemed as the key intermediates for the OER if the water solvent is not available timely on the non-solvated Fe-terminated Fe2O3 (0001) surface. The former intermediate has been verified with the experimental assignment of the Fe = O group as a key intermediate in PEC water oxidation reaction24 and also agrees with the simulation analysis of in-situ optical absorbance spectra46. Some experimental clues support the existence of the second intermediate. From the operando IR absorption spectra reported by Omid Zandi et al24., the peak at 743 cm−1 can be interpreted as the bending vibration (rocking) of water molecules in the *OΔ + *OH2 configuration, as detailed in Supplementary Note 4. Besides, Yuchao Zhang et al.25 reported that a peak around 1550 cm−1, attributed to the scissoring vibration of water molecules, becomes more prominent with increasing potential.

OER on Solvated Fe-terminated Fe2O3 (0001) Surface

Adsorption and dissociation of the 1st water molecule on the solvated Fe-terminated Fe2O3 (0001) surface were shown in Fig. 5a, b. Compared to the same process on the non-solvated surface, the molecular water adsorption can be stabilized during the period of the simulation. The hydrogen bonding of the solvent water molecule to the surface lattice Os or the adsorbed water molecule to the solvent water molecule is beneficial to stabilize the adsorbed water molecule. However, the dissociated water with adsorption of another water molecule is energetically preferred, giving equation C’ in Section 3.2, justifying the experimental observation of majority of dissociative water adsorption and minority of molecular water adsorption53. This configuration was not reported previously55. The HOMO charge density of this favorable configuration was shown in Fig. 3e, which spreads over the slab instead of the surface hydroxyl or adsorbed water molecule and facilitates removal of the H bonding to the surface lattice Os in the following 1st PCET step.

Fig. 5: Intermediates for OER on solvated Fe-terminated Fe2O3 (0001) surface determined with AIMD simulations.
figure 5

a Solvated molecular water adsorption; b Dissociative water adsorption in conjunction with adsorption of another water molecule, namely equation C’ in Section 3.2; c Formation of a hydroxyl group (*OH + *OH2), namely equation D’ in Section 3.2; d Formation of a bi-hydroxyl group (*OH + *OH), namely equation E’ in Section 3.2; e Formation of a peroxo-like group (*OΔ + *OH2); f Solvated *OOH + ΔH species, namely equation F’ in Section 3.2; g Solvated *OOH species; h Formation of a superoxo-like species; i Formation of an adsorbed O2 molecule. (Supplementary Note 5). Source data are provided as a Source Data file. The optimized structures are provided in Supplementary Data 1.

Full size image

The product of the 1st PCET step on the solvated surface, equation D’ in Section 3.2, was a hydroxyl group in conjunction with an adsorbed water molecule, namely *OH + *OH2, as shown in Fig. 5c. Figure 3f shows the HOMO charge density of the *OH + *OH2 configuration. The HOMO is localized on the O of the hydroxyl group rather than the adsorbed water molecule, suggesting that the hydroxyl group is likely to capture the photo-generated holes for the following oxidation reactions.

The product of the 2nd PCET step, namely *O + *OH2, spontaneously transforms into a bi-hydroxyl group (*OH + *OH), equation E’ in Section 3.2. Figure 5d illustrated this process with the AIMD simulation starting with the tilted ferryl group on the solvated surface. The peroxo-like species in conjunction with an adsorbed water molecule at the surface Fe, namely *OΔ + *OH2, can also be stabilized during the period of simulation, where a solvent water molecule automatically adsorbs at the surface Fe, as shown in Fig. 5e. The *OOH species at the surface Fe with a hydroxyl group of surface lattice Os, namely *OOH + ΔH, can also be stabilized on the solvated surface in our AIMD simulation, as shown in Fig. 5f.

The 3rd PCET step gives the reaction product of the *OOH species or the superoxo-like species (*OO + ΔH). In contrast to the instability of the *OOH species on the non-solvated surface, Fig. 5g shows that the *OOH species can be stabilized on the solvated surface due to the hydrogen bonding to the water solvent. Figure 5h shows the superoxo-like species on the solvated surface, which is energetically preferred to the *OOH species. Interestingly, the superoxo-like species pushes the solvent water molecules away and thus exhibits hydrophobicity. Yuchao Zhang et al. have experimentally verified the existence of superoxo intermediates stabilized with hydrogen bonding25. Here, we suggest that the hydrogen bonding should be from the H binding to the lattice Os instead of from the water solvent, due to the hydrophobicity of the superoxo-like species.

Once the surface H was oxidized, the 4th PCET step, the superoxo-like species spontaneously transforms into an adsorbed O2 molecule, as shown in Fig. 5i. The adsorbed O2 molecule is also hydrophobic. Since desorption of the O2 molecule is an endothermic non-electrochemical step, the O2 molecules can accumulate on the surface and impede the solvent water molecules approaching the surface, thus diminishing the influence of water solvent. This hydrophobic surface rationalizes the existence of the tilted ferryl group in our simulation on the non-solvated surface and also in the operando spectra experiment24,46.

Figure 6a depicted the reaction mechanism for OER on the solvated Fe-terminated Fe2O3 (0001) surface. The reaction pathway is composed of nine reaction steps, four electrochemical steps (equations D’, E’, H, and I) and five non-electrochemical steps (equations A, B, C’, F’, and J). Compared to the non-solvated surface, the water layers change the mechanism for OER on the solvated surface. Firstly, the water solvent participates in the OER at the very beginning. The dissociative adsorption of the 1st water molecule occurs with one solvent water molecule adsorbing at the surface Fe, equation C’, and the adsorbed solvent water molecule remains adsorbing at the surface Fe in the 1st PCET step, equation D’. For this reason, the 2nd PCET step, equation E’, directly leads to formation of a bi-hydroxyl group which then restructures to the *OOH species with a relatively small activation energy, equation F’. In this mechanism, both the ferryl and peroxo-like groups were avoided. Furthermore, the water solvent helps form some sub-stable species, such as the adsorbed 1st water molecule and the *OOH species. Finally, both the superoxo-like species and the adsorbed O2 molecule are hydrophobic, which occupies the surface Fe site and prevents the solvent water molecules approaching the surface. As a result, the mechanism for the OER on the non-solvated surface is also available on the solvated surface if the adsorbed O2 molecules can not be timely removed from the surface.

Fig. 6: Reaction free energy diagrams for OER on solvated Fe-terminated Fe2O3 (0001) surface with the error bars marked, i.e., the standard deviation of the DFT energy for a set of the snapshots sampled from the simulation trajectory of each intermediate.
figure 6

a Without electric bias. b With an electric bias of 1.97 V. Source data are provided as a Source Data file. The optimized structures are provided in Supplementary Data 1.

Full size image

Figure 6b shows the four electrochemical steps become downhill in Gibbs free energy with an applied bias of 1.97 V, giving an overpotential of 0.74 eV which is smaller than that on the non-solvated surface. Among the five non-electrochemical steps, only equations C’ and J are somewhat endothermic with small thermodynamic energy barriers of 0.06 and 0.31 eV and without noticeable activation energies. Therefore, the mechanism on the solvated surface should proceed more readily than that on the non-solvated surface. To enhance the OER activity on Fe2O3, timely removal of the surface O2 molecules is required, which can be achieved by ultrasound or heating techniques. This observation has been demonstrated experimentally in other OER catalytic systems56,57. The improvement in the reaction activity stems from not only exposure of more surface reaction sites but also hydrophilicity of catalyst surfaces which leads to switching in the reaction mechanisms.

In conclusion, using the reaction intermediates identified through AIMD simulations, we propose detailed reaction mechanisms for OER on both the non-solvated and solvated Fe-terminated Fe2O3 (0001) surfaces. The mechanism on the non-solvated surface consists of ten reaction steps, i.e., four electrochemical steps and six non-electrochemical steps, with the highest thermodynamic energy barrier of 2.10 eV. In this mechanism, the tilted ferryl group and the adsorbed 2nd water molecule were deemed as the key intermediates for the OER, in case that the water solvent cannot be timely available during the reaction. Formation of a peroxo-like species (O–Os bond between the adsorbed and surface lattice oxygen) from the tilted ferryl group requires a kinetic activation energy of 0.83 eV, which rationalizes the experimental observation of accumulation of the ferryl groups under PEC water oxidation condition. Moreover, formation of the *OOH species by restructuring the adsorbed 2nd water molecule requires a higher kinetic activation energy of 1.04 eV. The existence of this intermediate can be speculated from several experimental clues. On the other hand, the water solvent changes the OER mechanism on the solvated surface, with the highest thermodynamic energy barrier of 1.97 eV. The water solvent not only participates in the OER at the very beginning but also helps form some sub-stable species, such as the adsorbed 1st water molecule and the *OOH species. In this mechanism, the bi-hydroxyl group can be deemed as the key intermediates. Formation of the *OOH species by restructuring the bi-hydroxyl group requires a lower kinetic activation energy of 0.63 eV. Consequently, the mechanism on the solvated surface is likely to proceed more readily than that on the non-solvated surface. However, our simulation suggests that both the superoxo-like species and the adsorbed O2 molecule are hydrophobic, which prevents the solvent water molecules approaching the surface and thus makes the mechanism on the non-solvated surface possible under real reaction conditions. Therefore, we expect that desorption of the surface O2 molecules may be a key factor controlling the OER dynamics on Fe2O3 surfaces. Timely removal of the surface O2 molecules not only exposes more surface reaction sites but also switches the reaction mechanism by improving the surface hydrophilicity.

This work provides a potential approach for proposing OER mechanisms on the Fe-terminated Fe2O3 (0001) surface by identifying the key reaction intermediates with AIMD simulations at the level of DFT + U. We expect that this approach may be applicable to proposing the mechanisms of similar reactions on the surfaces of other materials, such as transition metal oxides with strong electronic correlation.

Methods

Computational Setups

Structural Models: The Fe-terminated α-Fe2O3 (0001) surface was modeled by a symmetrical slab which was cleaved from the bulk crystal structure. The slab is composed of ten layers of Fe atoms and five layers of O atoms, with a vacuum layer of 15 Å in the z direction, as shown in Fig. 7a. The antiferromagnetic ground state was reserved for the slab and the magnetic moments of Fe atoms were set consistently with our previous work58. To mimic the water layers, referring to previous studies59,60, eight water molecules were introduced in the vacuum layer, see Fig. 7b. Two of these water molecules will participate in the OER on each surface of the slab.

Fig. 7: Structural models for a non-solvated and b solvated Fe-terminated Fe2O3 (0001) surface.
figure 7

Structural models for (a) non-solvated and (b) solvated Fe-terminated Fe2O3 (0001) surface. c Evolution of the potential energy of the fully hydroxylated O-terminated Fe2O3 (0001) surface with several snapshots sampled from a set of local minima in the AIMD trajectory. d Schematic illustrations of the shallow potential energy surface of a hydrogen bond system, where the energy barriers can be readily overcome in the AIMD simulation at room temperature. Source data are provided as a Source Data file. The models are provided in Supplementary Data 1.

Full size image

A larger simulation cell with 2 × 2 periodicity in the surface plane and 32 water molecules in the condensed layer is used to check the effects of the cell size and the number of water molecules on the key intermediates in OER, for instance, the *O and *OΔ + *OH2 intermediates on the non-solvated Fe-terminated Fe2O3 (0001) surface and the *OH2, *OH + *OH2 + ΔH, and *OO + ΔH intermediates on the solvated Fe-terminated Fe2O3 (0001) surface. The simulation results are shown in Supplementary Figs. 7 and 8. On both the non-solvated and solvated surfaces, no new intermediates are found and the properties of the intermediates remain unchanged, for instance, the hydrophobicity of *OO + ΔH on the solvated surface, indicating that the 1 × 1 unit cell with 8 water molecules constituting the condensed layer is sufficient for the simulations in this work.

Computational Details: The electronic structure was calculated using Vienna Ab-initio Simulation Package (VASP) in the framework of spin-polarized plane-wave density functional theory (DFT)61,62. The rotationally invariant PBE + U 63 approach was used to account for the strong correlation among localized Fe 3 d electrons. With U = 5 eV and J = 1 eV, a number of ground and excited state properties of Fe2O3 have been successfully predicted, including lattice constants, bond length, band gap64, formation and hopping of electron small polarons58,65, and even electron-hole relaxation and recombination dynamics66,67,68,69. It has been reported in literature that the Fe2O3-water interface and other Fe-related catalyst-water interface can also be described successfully with DFT + U 59,70,71. Although the hybrid functional such as PBE0 and HSE appears effective to treat the self-interaction error, the computational cost is prohibitive for AIMD simulations. All-electron frozen core projector augmented wave (PAW) potentials72 were employed to treat the electron-core interaction. The valence configurations of Fe and O atoms were 3d64s2 and 2s22p4, respectively. The Fe PAW potential with 3p semi-core electrons was tested, and its effect can be safely neglected, as shown in Supplementary Table 1. The energy cutoff of 450 eV was used to expand the electronic wave function. Gaussian smearing method was used with a smearing width of 0.05 eV. The bulk structure optimization was performed using a Monkhorst−Pack k-point grid mesh of 4 × 4 × 2, in which both the unit cell and atomic positions were allowed to relax. Surface geometry optimization was performed using a 4  × 4 × 1 Monkhorst−Pack k-point grid mesh, while a denser k-point grid mesh of 8 × 8 × 1 was used for electronic structure calculation.

A number of AIMD simulations were carried out using VASP with the simulation time of 10~20 ps, which depends on the system under investigation, and the time step of 1.0 fs. The temperature was maintained at 300 K using Nosé–Hoover thermostat73 within the canonical ensembles (NVT). As shown in Fig. 7c, the AIMD simulation on the fully hydroxylated O-terminated Fe2O3 (0001) surface appears to rapidly converge to stable surface structures, which, however, may be challenging to achieve using a uniform sampling method39. The capability of the AIMD simulation at room temperature to search for the global minima in the potential energy surface of a hydrogen bonding system is expected to stem from easy formation and breakage of hydrogen bonds at the oxide-water interface, since the hydrogen bond energies are usually low (20 kJ·mol−1), see Fig. 7d. In each of the AIMD runs, several snapshots were sampled from a set of local minima in the trajectory and then were structurally optimized.

The minimal energy pathways were calculated using the climbing image nudged elastic band (CI–NEB) method74 and the Dimer method as implemented in VASP of the VTST version75,76. Four image points were used to construct the reaction pathway. Frequency calculations were performed to verify the reliability of the transition states searched.

Mechanisms of Water oxidation

Different from the simplified four-steps PCET mechanism, we proposed a more detailed reaction mechanism based on the intermediates from the AIMD simulations. The reaction steps constituting the most probable pathway on the non-solvated surface were listed as follows:

$${{{\rm{H}}}}_{2}{{\rm{O}}}+\ast \to*{{{\rm{OH}}}}_{2}$$
(A)
$$\ast {{{\rm{OH}}}}_{2}{+}^{\Delta }\to*{{\rm{OH}}}+{}^{\Delta }{{\rm{H}}}$$
(B)
$$\ast {{\rm{OH}}}+{}^{\Delta }{{\rm{H}}}\to*{{\rm{OH}}}{+}^{\Delta }+{{{\rm{H}}}}^{+}+{{{\rm{e}}}}^{-}$$
(C)
$$\ast {{\rm{OH}}}\to*{{\rm{O}}}+{{{\rm{H}}}}^{+}+{{{\rm{e}}}}^{-}$$
(D)
$$\ast {{\rm{O}}}{+}^{\Delta }\to*{{{\rm{O}}}}^{\Delta }$$
(E)
$${{{\rm{H}}}}_{2}{{\rm{O}}}+\ast {{{\rm{O}}}}^{\Delta }\to*{{{\rm{OH}}}}_{2}+\ast {{{\rm{O}}}}^{\Delta }$$
(F)
$$\ast {{{\rm{OH}}}}_{2}+\ast {{{\rm{O}}}}^{\Delta }\to*{{\rm{OOH}}}+{}^{\Delta }{{\rm{H}}}$$
(G)
$$\ast {{\rm{OOH}}}+{}^{\Delta }{{\rm{H}}}\to*{{\rm{OO}}}+{}^{\Delta }{{\rm{H}}}+{{{\rm{H}}}}^{+}+{{{\rm{e}}}}^{-}$$
(H)
$$\ast {{\rm{OO}}}+{}^{\Delta }{{\rm{H}}}\to*{{\rm{OO}}}{+}^{\Delta }+{{{\rm{H}}}}^{+}+{{{\rm{e}}}}^{-}$$
(I)
$$\ast {{\rm{OO}}}\to {{{\rm{O}}}}_{2}+\ast$$
(J)

Reaction I and J can be exchanged.

$$\ast {{\rm{OO}}}+{}^{\Delta }{{\rm{H}}}\to \ast+{}^{\Delta }{{\rm{H}}}+{{{\rm{O}}}}_{2}$$
(I’)
$${}^{\Delta }{{\rm{H}}}{\to }^{\Delta }+{{{\rm{H}}}}^{+}+{{{\rm{e}}}}^{-}$$
(J’)

* stands for the surface Fe site and Δ for the surface lattice O site. *OH2, *OH, *O, *OOH, and *OO stand for the reaction intermediates adsorbed at the surface Fe, ΔH for the H of water adsorbed at the surface lattice O, and *OΔ for the O of water bonded with both surface Fe and lattice O. The Gibbs free energy for each reaction step was calculated according to the following equations.

$$\Delta {G}_{{{\rm{A}}}}=\frac{1}{2}{E}_{*{\mbox{O}}{{\mbox{H}}}_{2}}-\frac{1}{2}\left({E}_{*}+2{E}_{{{{\rm{H}}}}_{2}{{\rm{O}}}}\right)+\Delta {\mbox{ZP}}{{\mbox{E}}}_{{{\rm{A}}}}-{\rm{T}}\Delta {\rm{S}}_{{{\rm{A}}}}$$
(1)
$$\Delta {G}_{{{\rm{B}}}}=\frac{1}{2}{E}_{(*{\mbox{OH}}+{}^{\Delta }{{\rm{H}}})}-\frac{1}{2}{E}_{(*{\mbox{O}}{{\mbox{H}}}_{2}+\Delta )}+\Delta {\mbox{ZP}}{{\mbox{E}}}_{{{\rm{B}}}}-{\rm{T}}\Delta {\rm{S}}_{{{\rm{B}}}}$$
(2)
$$\Delta {G}_{{{\rm{C}}}}=\frac{1}{2}\left({E}_{\left(*{\mbox{OH}}{+}^{\Delta }\right)}+{E}_{{{{\rm{H}}}}_{2}}\right)-\frac{1}{2}{E}_{\left(*{\mbox{OH}}+{}^{\Delta }{{\rm{H}}}\right)}+\Delta {\mbox{ZP}}{{\mbox{E}}}_{{{\rm{C}}}}-{\rm{T}}\Delta {\rm{S}}_{{{\rm{C}}}}-{\rm{eU}}$$
(3)
$$\Delta {G}_{D}=\frac{1}{2}\left({E}_{*{\mbox{O}}}+{E}_{{{{\rm{H}}}}_{2}}\right)-\frac{1}{2}{E}_{*{\mbox{OH}}}+\Delta {\mbox{ZP}}{{\mbox{E}}}_{{{\rm{D}}}}-{\rm{T}}\Delta {\rm{S}}_{{{\rm{D}}}}-{\rm{eU}}$$
(4)
$$\Delta {G}_{{{\rm{E}}}}=\frac{1}{2}{E}_{*{{{\rm{O}}}}^{\Delta }}-\frac{1}{2}{E}_{\left(*{\mbox{O}}{+}^{\Delta }\right)}+\Delta {\mbox{ZP}}{{\mbox{E}}}_{{{\rm{E}}}}-{\rm{T}}\Delta {\rm{S}}_{{{\rm{E}}}}$$
(5)
$$\Delta {G}_{{{\rm{F}}}}=\frac{1}{2}{E}_{\left(*{\mbox{O}}{{\mbox{H}}}_{2}+\ast {{{\rm{O}}}}^{\Delta }\right)}-\frac{1}{2}\left({E}_{*{{{\rm{O}}}}^{\Delta }}+2{E}_{{{{\rm{H}}}}_{2}{{\rm{O}}}}\right)+\Delta {\mbox{ZP}}{{\mbox{E}}}_{{{\rm{F}}}}-{\rm{T}}\Delta {\rm{S}}_{{{\rm{F}}}}$$
(6)
$$\Delta {G}_{{{\rm{G}}}}=\frac{1}{2}{E}_{\left(*{\mbox{OOH}}+{}^{\Delta }{{\rm{H}}}\right)}-\frac{1}{2}{E}_{\left(*{\mbox{O}}{{\mbox{H}}}_{2}+\ast {{{\rm{O}}}}^{\Delta }\right)}+\Delta {\mbox{ZP}}{{\mbox{E}}}_{{{\rm{G}}}}-{\rm{T}}\Delta {\rm{S}}_{{{\rm{G}}}}$$
(7)
$$\Delta {G}_{{{\rm{H}}}}=\frac{1}{2}\left({E}_{\left(*{\mbox{OO}}+{}^{\Delta }{{\rm{H}}}\right)}+{E}_{{{{\rm{H}}}}_{2}}\right)-\frac{1}{2}{E}_{\left(*{\mbox{OOH}}+{}^{\Delta }{{\rm{H}}}\right)}+\Delta {\mbox{ZP}}{{\mbox{E}}}_{{{\rm{H}}}}-{\rm{T}}\Delta {\rm{S}}_{{{\rm{H}}}}-{\rm{eU}}$$
(8)
$$\Delta {G}_{{{\rm{I}}}}=\frac{1}{2}\left({E}_{\left(*{\mbox{OO}}{+}^{\Delta }\right)}+{E}_{{{{\rm{H}}}}_{2}}\right)-\frac{1}{2}{E}_{\left(*{\mbox{OO}}+{}^{\Delta }{{\rm{H}}}\right)}+\Delta {\mbox{ZP}}{{\mbox{E}}}_{{{\rm{I}}}}-{\rm{T}}\Delta {\rm{S}}_{{{\rm{I}}}}-{\rm{eU}}$$
(9)
$$\Delta {G}_{{{\rm{J}}}}=\frac{1}{2}\left({E}_{*}+2{E}_{{{{\rm{O}}}}_{2}}\right)-\frac{1}{2}{E}_{*{\mbox{OO}}}+\Delta {\mbox{ZP}}{{\mbox{E}}}_{{{\rm{J}}}}-{\rm{T}}\Delta {\rm{S}}_{{{\rm{J}}}}$$
(10)
$$\Delta {G}_{{{\rm{I}}}^{\prime}}=\frac{1}{2}\left({E}_{\left( \ast+{}^{\Delta }{{\rm{H}}}\right)}+2{E}_{{{{\rm{O}}}}_{2}}\right)-\frac{1}{2}{E}_{\left(*{{\mbox{OO}}}+{}^{\Delta }{{{\rm{H}}}}\right)}+\Delta {{\mbox{ZP}}}{{\mbox{E}}}_{{{\mbox{I}}^{\prime}} }-{\rm{T}}\Delta {\rm{S}}_{{{\rm{I}}}^{\prime}}$$
(11)
$$\Delta {G}_{{{\rm{J}}^{\prime}}}=\frac{1}{2}\left({E}_{\Delta }+{E}_{{{{\rm{H}}}}_{2}}\right)-\frac{1}{2}{E}_{{}^{\Delta }{{\rm{H}}}}+\Delta {\mbox{ZP}}{{\mbox{E}}}_{{\mbox{J}}{\prime} }-{\rm{T}}\Delta^{\prime} -{\rm{eU}}$$
(12)

E is the single-point energy. U denotes the applied potential. ZPE is the zero-point energy and S is the vibrational entropy. The ZPE corrections and entropy contributions were obtained from vibrational frequencies77. The entropy contributions for adsorbed species are defined by considering only their vibrational motion, as described in Eq. (13)38:

$$S={k}_{B}{\mathrm{ln}}{\prod }_{i=1}^{N}\frac{1}{1-\exp \left(\frac{{{\hslash }}{\omega }_{X,i}}{{k}_{B}T}\right)}$$
(13)

where kB is the Boltzmann constant, ℏ is the Planck constant, ω is wavenumber, i represents a specific vibrational mode, and N is the total number of vibrational modes for the adsorbed specie X. The entropic contributions for gaseous molecules are taken from standard thermodynamics tables78. The entropic contribution of H2O was calculated in the gas phase under a pressure of 0.035 bar, because at this pressure gas phase is in equilibrium with liquid water at 300 K. The ZPE corrections and entropies of other intermediates were calculated by the vibrational analysis at T = 300 K using VASPKIT79.

Reaction steps C-G were replaced by C’-F’ on the solvated surface:

$$\ast {{\rm{OH}}}+{}^{\Delta }{{\rm{H}}}+{{{\rm{H}}}}_{2}{{\rm{O}}}\to*{{\rm{OH}}}+\ast {{{\rm{OH}}}}_{2}+{}^{\Delta }{{\rm{H}}}$$
(C’)
$$\ast {{\rm{OH}}}+\ast {{{\rm{OH}}}}_{2}+{}^{\Delta }{{\rm{H}}}\to*{{\rm{OH}}}+\ast {{{\rm{OH}}}}_{2}{+}^{\Delta }+{{{\rm{H}}}}^{+}+{{{\rm{e}}}}^{-}$$
(D’)
$$\ast {{\rm{OH}}}+\ast {{{\rm{OH}}}}_{2}\to*{{\rm{OH}}}+\ast {{\rm{OH}}}+{{{\rm{H}}}}^{+}+{{{\rm{e}}}}^{-}$$
(E’)
$$\ast {{\rm{OH}}}+\ast {{\rm{OH}}}\to*{{\rm{OOH}}}+{}^{\Delta }{{\rm{H}}}$$
(F’)

The Gibbs free energy for each reaction step was calculated using Eqs. (14)–(17):

$$\Delta {G}_{{{\rm{C}}^{\prime}}}=\frac{1}{2}{E}_{\left(*{\mbox{OH}}*{\mbox{O}}{{\mbox{H}}}_{2}+{}^{\Delta }{{\rm{H}}}\right)}-\frac{1}{2}\left({E}_{\left(*{\mbox{OH}}+{}^{\Delta }{{\rm{H}}}\right)}+2{E}_{{{{\rm{H}}}}_{2}{{\rm{O}}}}\right)+\Delta {\mbox{ZP}}{{\mbox{E}}}_{{\mbox{C}}^{\prime} }-{\rm{T}}\Delta {\rm{S}}_{{{\rm{C}}^{\prime}}}$$
(14)
$$\Delta {G}_{{{\rm{D}}^{\prime}}}= \frac{1}{2}\left({E}_{\left(*{\mbox{OH}}*{\mbox{O}}{{\mbox{H}}}_{2}{+}^{\Delta }\right)}+{E}_{{{{\rm{H}}}}_{2}}\right)\\ -\frac{1}{2}{E}_{\left(*{\mbox{OH}}*{\mbox{O}}{{\mbox{H}}}_{2}+{}^{\Delta }{{\rm{H}}}\right)}+\Delta {\mbox{ZP}}{{\mbox{E}}}_{{\mbox{D}}^{\prime} }-{\rm{T}}\Delta {\rm{S}}_{{{\rm{D}}^{\prime}}} -{\rm{eU}}$$
(15)
$$\Delta {G}_{{{\rm{E}}^{\prime}}}=\frac{1}{2}\left({E}_{*{\mbox{OH}}*{\mbox{OH}}}+{E}_{{{{\rm{H}}}}_{2}}\right)-\frac{1}{2}{E}_{*{\mbox{OH}}*{\mbox{O}}{{\mbox{H}}}_{2}}+\Delta {\mbox{ZP}}{{\mbox{E}}}_{{\mbox{E}}^{\prime} }-{\rm{T}}\Delta {\rm{S}}_{{{\rm{E}}^{\prime}}} -{\rm{eU}}$$
(16)
$$\Delta {G}_{{{{\rm{F}}}^{\prime}}}=\frac{1}{2}{E}_{\left(*{{\mbox{OOH}}}+{}^{\Delta }{{\rm{H}}}\right)}-\frac{1}{2}{E}_{*{{\mbox{OH}}}*{{\mbox{OH}}}}+\Delta {{\mbox{ZP}}}{{\mbox{E}}}_{{{\mbox{F}}}^{\prime} }-{\rm{T}}\Delta {\rm{S}}_{{{\rm{F}}^{\prime}}}$$
(17)

The calculated ZPE and entropy contributions for each intermediates are shown in Supplementary Tables 2 and 3. Since DFT at the level of PBE presents an over-binding nature for the O2 molecule80,81, the Gibbs free energy of the reaction step involving O2, namely reaction J and I’, was corrected according to Eqs. (18) and (19) in terms of the Gibbs free energy of 4.92 eV for the overall water splitting.

$$\Delta {G}_{J}=4.92-\left(\Delta {G}_{{{\rm{A}}}}+\Delta {G}_{{{\rm{B}}}}+\Delta {G}_{{{\rm{C}}}}+\Delta {G}_{{{\rm{D}}}}+\Delta {G}_{{{\rm{E}}}}+\Delta {G}_{{{\rm{F}}}}+\Delta {G}_{{{\rm{G}}}}+\Delta {G}_{{{\rm{H}}}}+\Delta {G}_{{{\rm{I}}}}\right)$$
(18)
$$\Delta {G}_{{{\rm{I}}}^{\prime}}=4.92-\left(\Delta {G}_{{{\rm{A}}}}+\Delta {G}_{{{\rm{B}}}}+\Delta {G}_{{{\rm{C}}}}+\Delta {G}_{{{\rm{D}}}}+\Delta {G}_{{{\rm{E}}}}+\Delta {G}_{{{\rm{F}}}}+\Delta {G}_{{{\rm{G}}}}+\Delta {G}_{{{\rm{H}}}}+\Delta {G}_{{{\rm{J}}}^{\prime}} \right)$$
(19)