First-Principles Design of Rutile Oxide Heterostructures for Oxygen Evolution Reactions

The oxygen evolution reaction (OER) plays a key role in the determination of overall water-splitting rate. Lowering the high overpotential of the OER of transition metal oxides (TMOs), which are used as conventional OER electrocatalysts, has been the focus of many studies. The OER activity of TMOs can be tuned via the strategic formation of a heterostructure with another TMO substrate. We screened 11 rutile-type TMOs (i.e., MO2; M = V, Cr, Mn, Nb, Ru, Rh, Sn, Ta, Os, Ir, and Pt) on a rutile (110) substrate using density functional theory calculations to determine their OER activities. The conventional volcano approach based on simple binding energies of reaction intermediates was implemented; in addition, the electrochemical-step symmetry index was employed to screen heterostructures for use as electrode materials. The results show that RuO2 and IrO2 are the most promising catalysts among all candidates. The scaling results provide insights into the intrinsic properties of the heterostructure as well as materials that can be used to lower the overpotential of the OER.


INTRODUCTION
Green hydrogen production remains a challenge that must be overcome to achieve a hydrogen economy (Turner, 2004). Water electrolysis is one of the approaches toward hydrogen ecofriendliness. It is based on the hydrogen evolution reaction (HER) and oxygen evolution reaction (OER) as cathodic and anodic reactions, respectively. However, the slow kinetics of the OER limit the commercialization of this approach. Oxide systems based on Ir or Ru (i.e., IrO x or RuO x ) are known to have the best OER performances (Lee et al., 2012;Frydendal et al., 2014;Suen et al., 2017). However, these materials are expensive. Therefore, many studies have been conducted to identify a cost-effective alternative oxide material with high activity. For example, transition-metal substitution (García-Mota et al., 2011) and the introduction of oxygen vacancies (Xiao et al., 2020) have been studied as means to control the compositions of expensive materials. However, transitionmetal substitution and oxygen vacancies only locally affect the active sites.
To activate the metal oxide, support materials are mixed with the oxide to enhance the electrical conductivity and activity of the electrocatalytic reaction (Kumar et al., 2016;Qingxiang et al., 2018;Bu et al., 2019;Lu et al., 2019). In addition, research on heterostructures is conducted with the goal of activating the material on the substrate. Modifications are mainly based on the strain effect and charge transfer due to the formation of an interface between the support material and catalyst. For instance, the heterostructure of La 0.5 Sr 0.5 CoO 3-δ and MoSe 2 induces the phase transition of MoSe 2 from the 2H to the 1T phase (Oh et al., 2019). In addition, the charge transfer from the Co ion to Mo ion improves the electrochemical activity. The results of another study showed that the heterostructure of IrO x and SrIrO 3 outperforms iridium or ruthenium oxide systems (Seitz et al., 2016). Based on density functional theory (DFT) calculations, IrO 3 or anatase IrO 2 motifs are formed during strontium leaching in the outermost surface layers of SrIrO 3 and contribute to the high activity. Another strength of the heterostructure is that it can possibly contribute to reducing the use of precious metals by replacing them into low-cost metals while maintaining similar intrinsic activity of each active site (Esposito et al., 2010;Zhou et al., 2014;Wang et al., 2015;Jin et al., 2016). Rutile (TiO 2 ) is a well-known substrate material for metal oxide systems because of its high structural stability and simple structure (Hanaor et al., 2012). It is suitable for the growth of oxide films, facilitating the fabrication of heterostructures with various metal oxides. In addition, rutile affects the catalytic activity of metal oxides and has a high cost effectiveness (Seitsonen and Over, 2010;Stacchiola et al., 2013;Wei et al., 2015;Sun et al., 2016;Li et al., 2017). However, the effect of the rutile support on the activity remains unclear due to difficulties with respect to the characterization of the heterostructure during experiments (Stacchiola et al., 2013).
In this study, we performed DFT calculations to theoretically investigate the effect of the TiO 2 substrate on the heterostructure and thus on the OER. We screened the OERs of heterostructures consisting of various rutile-type metal oxides (i.e., VO 2 , CrO 2 , MnO 2 , NbO 2 , RuO 2 , RhO 2 , SnO 2 , TaO 2 , OsO 2 , IrO 2 , and PtO 2 ) and a TiO 2 substrate. The results show that these rutile heterostructures follow the universal scaling relations of metals and oxides; however, the binding strengths of the O* intermediates increase due to the TiO 2 substrate. The volcano plot and electrochemical-step symmetry index (ESSI) show that the RuO 2 and IrO 2 are the closest to an ideal catalyst. The results of our computational screening provide insights into the effects of support materials on electrocatalytic reactions.

METHODS
Spin-polarized DFT calculations were performed with the projector-augmented wave (PAW) method (Blochl, 1994) and Vienna Ab initio Simulation Package (VASP) (Kresse and Furthmiiller, 1996). The electron exchange-correlation energy was treated within the generalized gradient approximation (GGA) and the Perdew-Burke-Ernzerhof functional (Perdew et al., 1997). To determine the trend of the OER activity, the DFT + U method within Liechtenstein's approach (Anisimovdag et al., 1997) was used and the following correction parameters were employed: U 4.95 eV for Ti, 2 eV for V, 7.15 eV for Cr, 6.63 eV for Mn, 3.32 eV for Nb, 6.73 eV for Ru, 5.97 eV for Rh, 5.91 eV for Ir, and 6.25 eV for Pt (Xu et al., 2015). The energy cutoff for the plane wave basis set was set at 520 eV. The geometry was optimized using the residual minimization method and the direct inversion in the iterative subspace method (RMM-DIIS) algorithm until the net force on each atom was below 0.02 eV·Å −1 , and the total energy was 10 -6 eV per atom. Dipole slab corrections were also applied to all slab model calculations. The k-point sampling of the Brillouin zone was done with a 4 × 4 × 1 for bulk calculations and 6 × 6 × 8 for slab calculations.

RESULTS AND DISCUSSION
To theoretically investigate the OER activities of rutile-type heterostructures, we considered the heterostructures of 11 rutile-type oxides (i.e., MO 2 , where M V, Cr, Mn, Nb, Ru, Rh, Sn, Ta, Os, Ir, and Pt) with a TiO 2 substrate. To accurately illustrate the OER activity, the following magnetic structures were used for all calculations according to Xu et al.'s work: nonmagnetic (NM) for TiO 2 , NbO 2 , RuO 2 , RhO 2 , IrO 2 , and PtO 2 , and ferromagnetic (FM) for CrO 2 and MnO 2 (Xu et al., 2015). Furthermore, the ground state magnetic configurations of the other candidate models were identified as FM for VO 2 and NM for SnO 2 , TaO 2 , and OsO 2 , respectively (Supplementary Figure S1). Surface models, that is, six-layer stoichiometric slabs, were built by using the 2 × 1 supercell of the optimized unit cell. The vacuum of ∼15 Å was applied in the (110) direction, which is the most stable facet of rutile-type oxides ( Figure 1A; Kung, 1989). Note that a tri-layer consisting of oxygen-metal-oxygen (O-M-O) atomic layers was considered to be a single layer in our slab models. The top four layers of the slab models were allowed to fully relax. The two layers at the bottom were fixed to represent the bulk state. For heterostructures, the top three layers of the TiO 2 slab model were replaced with MO 2 layers while maintaining the cell dimensions of the TiO 2 substrate. The coordinatively unsaturated sites (CUS) of the metal atoms at the top surface were considered to be the adsorption sites for each intermediate (i.e., OH*, O*, and OOH*) for the OER. To examine the oxygen coverage effect on the OER activity, the pristine and fully O* covered surfaces were representatively compared (i.e., denoted as 2O b and 2O b 2O c , where the subscripts "b" and "c" represent the bridge sites and CUS, respectively).
The overall OER process consists of four elementary steps involving proton-coupled electron transfer (PCET; Figure 1B) (Hammes-Schiffer, 2015; Costentin and Savéant, 2017). In this study, we followed the conventional four-electron pathways with respect to the OER of rutile heterostructures, which can be described as follows: The adsorption free energy was calculated using the following equation: where ΔE is the binding energy of each reaction intermediate; ΔZPE is the zero-point vibrational enthalpy; −TΔS is the entropic correction at room temperature, −eU is the energy shift by the electrode potential, where U is the electrode applied potential relative to the standard hydrogen electrode (SHE) and e is the elementary charge transferred; and k B T ln a H + is used as correction for the free energy of H + ions, where a H + is the activity of the proton. In this study, we considered the standard conditions for Gibbs free energy calculations where E(OH * ), E(O * ), and E(OOH * ) represent the total energies of the slab models for each adsorbate, E(p) is the total energy of the bare slab, and E(H 2 O) and E(H 2 ) represent the total energies of an isolated water molecule and hydrogen gas, respectively. The differences in the Gibbs free energy (ΔG) of each step were calculated as follows: Finally, the theoretical overpotential of the OER (η OER ) can be calculated as follows: where U eq indicates the equilibrium potential of the OER (i.e., 1.23 V vs. reversible hydrogen electrode). Eqs 1-13 assume that thermodynamics of the reaction steps is a valid descriptor for reaction kinetics based on the Brønsted-Evans-Polanyi (BEP) relations , which refer to that a free-energy change in the transition states (kinetics) follows a change in the reaction heat (thermodynamics). Note that the connection between thermodynamics and kinetics is not always established (Kuo et al., 2017;Kuo et al., 2018), which requires further kinetic experiments or microkinetic modeling based on the activation barrier calculations for all plausible transition states.
We examined the scaling relations between the adsorption free energies of the reaction intermediates (i.e., OH*, O*, and OOH*; Supplementary Table S2) for all rutile-type catalysts of interest ( Figure 2A). The binding energies of OOH* and OH* are linearly correlated, with an offset of 3.20 eV. Note that our scaling trend is   (Rossmeisl et al., 2007).
To analyze the substrate effect of TiO 2 on the heterostructures, we divided the scaling relations into two groups, that is, (110) surfaces of MO 2 and MO 2 /TiO 2 heterostructures ( Figure 2B,C). Notably, the binding energies of OOH* and O* on the MO 2 /TiO 2 (110) surface are more concentrated than those of MO 2 . The scaling relation between OOH* and OH* of MO 2 /TiO 2 (red solid line) shows an increased intercept by 0.18 compared with the corresponding scaling relation of MO 2 , implying slightly weakened interactions between OOH* and the surface. Note that the scaling relations between the binding strengths of O* and OH* species on MO 2 and MO 2 /TiO 2 (blue solid lines) apparently differ. The slopes of MO 2 /TiO 2 (110) heterostructures are closer to 2 (i.e., double bonding nature of O*) than those of MO 2 (110) surfaces (Rossmeisl et al., 2007). This is due to the intensified binding strengths of O* intermediates in heterostructures compared with MO 2 surfaces (Divanis et al., 2020). This indicates that the TiO 2 substrate generally stabilizes the O* intermediates, which might lead to a decrease in η OER of rutile oxides involving weakly bound O* intermediates.
Based on the scaling relations between the binding energies of the reaction intermediates, a volcano plot was constructed, as  shown in Figure 3A. We chose ΔG O * − ΔG OH * (denoted as ΔG 2 ) as a descriptor, which is commonly used to predict the OER activity of the 4 e − reaction Krishnamurthy et al., 2018). The plot shows that the O* bonds of PtO 2 , OsO 2 , TaO 2 , and RhO 2 on the TiO 2 substrate are stronger than those without the TiO 2 substrate, and the activities change along the volcano curve. In the cases of both oxygen coverages of PtO 2 / TiO 2 and RhO 2 /TiO 2 , ΔG 2 , which is the potential-determining step for each system, was reduced, exhibiting the improved OER activities (i.e., smaller η OER ) compared to those of PtO 2 and RhO 2 , respectively. In the cases of OsO 2 , TaO 2 , and NbO 2 , the bonds of the O* intermediates are too strong in the presence of the TiO 2 substrate, resulting in a decreased OER activity. Similar to the previous studies (Rossmeisl et al., 2007;Man et al., 2011), RuO 2 and IrO 2 were both identified as the active OER catalysts among the MO 2 candidates (i.e., η OER < ∼ 0.5 V, Figure 3A). In the case of 2O b 2O c of RuO 2 , the OER activity maintains similarity in the presence of TiO 2 substrate, which is located near the top of the volcano. In general, the TiO 2 substrate stabilizes the adsorbates (i.e., OOH*, O*, and OH*) on 2O b 2O c of RuO 2 , which appears as stronger adsorption free energies by 0.54-0.70 eV ( Supplementary Table S1). Nevertheless, η OER is almost unchanged (i.e., Δη OER −0.01 V) because the free energy change in the potential determining step (i.e., ΔG OOH* − ΔG O* ) remains similar (Supplementary Table S2). In addition, it is also noteworthy that η OER for 2O b 2O c of IrO 2 is relatively similar, still exhibiting higher activity than other candidates except RuO 2 .
The changes in η OER induced by the TiO 2 substrate are summarized in Figure 3B. As the value becomes more negative, the activity of MO 2 /TiO 2 improves compared with that of MO 2 . Notably, 2O b 2O c of CrO 2 , MnO 2 , RuO 2 , RhO 2 , and PtO 2 shows an enhanced activity. In addition, 2O b of RhO 2 and PtO 2 also shows an improved activity on the TiO 2 substrate. However, regardless of the coverage, the TiO 2 substrate decreases the OER activity of VO 2 , OsO 2 , and IrO 2 . In particular, 2O b 2O c of VO 2 and 2O b of SnO 2 exhibit significantly increased overpotentials (Δη OER of 0.53 and 0.57 V, respectively).
To determine the promising free-energy regime for the OER, we further analyzed the free energies of the steps in terms of the ESSI, as suggested by Calle-Vallejo and coworkers (Govindarajan et al., 2018). The ESSI is an energetic descriptor that indicates the degree of similarity with an ideal catalyst, where all OER steps are perfectly symmetric at 1.23 eV. The ESSI is defined by the following equation and is only applied to steps with ΔG i (i 1, 2, 3, 4) ≥ 1.23 eV (denoted as ΔG p i ): ESSI 1 n n i ΔG p i − 1.23).
(14) Figure 4A shows η OER of each model as a function of the ESSI, representing a good linear correlation. 2O b of RuO 2 is the closest model to the ideal catalyst, followed by 2O b 2O c of RuO 2 and RuO 2 /TiO 2 . The 2O b of PtO 2 appears on the η OER ESSI (red dashed line) and has a null bar with respect to the ESSI, which means that only a single step is greater than 1.23 V (i.e., ΔG 2 3.13 eV). In the presence of the TiO 2 substrate, ΔG 2 for 2O b of PtO 2 is reduced, leading to a decrease of the ESSI and η OER . Meanwhile, 2O b 2O c of PtO 2 lies relatively far from the line of η OER ESSI, which has a wide bar with respect to the ESSI, and thus corresponds to a good candidate for optimization (Govindarajan et al., 2018). In practice, the TiO 2 substrate on 2O b 2O c of PtO 2 works for enhancing the OER activity with a decrease in both η OER and ESSI.
Based on the scaling relation between the ESSIs and η OER 's, we determined a promising group of OER catalysts using the ESSI-ΔG 2 activity map introduced by Exner (Exner, 2019; Figure 4B). The activity map is used to identify OER candidates by adjusting the ESSI threshold. The ESSI threshold (< 0.53 V) was determined from Supplementary Figure S2 by applying η OER of IrO 2 , which is the conventional OER catalyst. The free-energy regime was set to 1.20 eV < ΔG 2 < 2.06 eV by applying the standard deviation of ± 0.43 eV of the scaling relationship between ΔG 2 and ΔG 3 (Supplementary Figure S3) on 1.63 eV. The median value of 1.63 eV in the free-energy regime is determined by assuming the threshold electrode potential as the point where the experimental Tafel slope exceeds 59 mV/dec (i.e., η OER > 0.4 V), accounting for kinetics (Exner and Over, 2019). The green-colored area includes 2O b of RuO 2 , 2O b 2O c of RuO 2 , 2O b 2O c of RuO 2 /TiO 2 , 2O b of IrO 2 , 2O b 2O c of IrO 2 , 2O b 2O c of IrO 2 /TiO 2 , and 2O b 2O c of RhO 2 , which are the most promising candidates. Particularly, a portion of promising candidates for the OER (i.e., RuO 2 , RuO 2 /TiO 2 , IrO 2 , IrO 2 / TiO 2 , and RhO 2 /TiO 2 ) are more evidently classified on the ESSI-ΔG 2 map ( Figure 4B), while they are somewhat deviated from the apex in the volcano plot ( Figure 3A). This implies that conventional volcano analysis does not guarantee to predict the most active OER catalyst, whereas the ESSI-ΔG 2 activity map, which is based on the kinetic scaling relations, is a more robust descriptor for the OER (Exner, 2019;Exner and Over, 2019). The 2O b 2O c surface of IrO 2 exhibits high symmetry of the reaction steps comparable to RuO 2 . On the 2O b and 2O b 2O c surfaces of RhO 2 /TiO 2 , the values of ΔG 2 decrease compared to those of RhO 2 . Accordingly, the free-energy differences of the other steps are adjusted to compensate for the decrease in ΔG 2 , while maintaining ΔG 2 as a potential-determining step with reduced η OER (cf. Supplementary Table S2). Next, 2O b of IrO 2 /TiO 2 belongs to the yellow-colored region, which needs to be reconsidered for activity optimization (Govindarajan et al., 2019). The candidates in the region of ESSI > 0.53 V (i.e., redhighlighted region in Figure 4B) are classified as an inferior group, showing poor OER activities over all ranges of ΔG 2 due to the highly asymmetric free-energy changes of the OER intermediates (Exner, 2019).

CONCLUSION
We screened a variety of rutile oxide heterostructures based on a TiO 2 substrate using scaling relations and relevant descriptors to identify a promising OER catalyst. The scaling relations between the reaction intermediates demonstrate that the rutile-type MO 2 heterostructures follow the universal scaling relationship of metal oxides. In addition, the TiO 2 substrate stabilizes the O* bond on the (110) metal oxide surface. Based on the conventional volcano plot, RuO 2 and IrO 2 are found to be highly active OER catalysts as previously reported. Based on the ESSI descriptor, the superior activity of the RuO 2 can be attributed to the high symmetry of the reaction steps. Furthermore, based on the ESSI-ΔG 2 activity map, the candidates can be classified into an optimum group, a second promising group of OER catalysts with potential for optimization, and an inferior group that does not require particular attention.
The results of our computational screening using the scaling relations of rutile-type heterostructures provide valuable insights into the effect of the support material on the overpotential and thus guidelines for the design of a promising OER catalyst.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.