The role of metastability in enhancing water-oxidation activity

While metastability enhanced water-oxidation activity was experimentally reported, the reason behind this effect is still unclear. We determine here, using density functional theory calculations, oxygen evolution reaction overpotentials for a variety of defective (001) surfaces of three different perovskite materials. For all three, we find a large range of overpotentials for different reaction sites including also overpotentials at the top of the activity volcano. Assuming that these sites dominate the apparent catalytic activity, this implies that a large number of geometrically different reaction sites, as they occur when a catalyst is operated at the border of its stability conditions, can lead to a strong enhancement of the apparent activity. This also implies that a pre-treatment of the catalyst creating a variety of different reactive sites could lead to superior catalytic activities for thermodynamically stable materials.


INTRODUCTION
Hydrogen fuel production by (photo)electrochemical water splitting has been intensively studied as a route to convert electrical or solar energy to chemical energy. The bottleneck in (photo)electrochemical water splitting is the oxygen evolution reaction (OER). Previous studies of metal 1 and oxide 2 catalysts proposed a likely OER reaction mechanism to consist of four consecutive protoncoupled one-electron transfer (PCET) steps with reaction intermediates *OH, *O and *OOH as shown in equations 1-4. * Several experimental studies found an inverse correlation between the activity and the stability of catalyst materials for the OER. Markovic and co-workers investigated activity-stability trends of the OER for crystalline and amorphous RuO 2 3,4 as well as for different surface orientations of the perovskite SrRuO 3 5 . Based on the concentration of dissolved metal ions, they found a higher activity for less stable surfaces of both materials and concluded that a highly active OER catalyst should balance activity and stability. A mechanistic understanding of a catalyst's stability in enhancing the OER activity is however still elusive and actively debated. It was suggested that the OER activity is controlled by the density of surface defects rather than by the binding energy of the reaction intermediates on the perfect surface and that it may be linked to changes in oxidation state of Ru ions close to defects 5 . A very recent study on RuO 2 thin films did, however, not find such a fundamental relationship between a material's stability and its activity 6 . While a correlation between the dissolution and the activity for different surface orientations was initially observed, it was also shown that the active OER sites are decoupled from the fastest-corroding Ru sites 6 .
To fully understand the activity-stability relationship one would ideally study the OER on specific reaction sites of defective dissolving materials. Experimentally it is, however, very challenging to characterize the geometry and activity of specific reaction sites since spectroscopic methods probe the average state of a catalyst. Computationally, one could readily calculate the overpotential and hence the activity of a given reaction site with a specific local geometry, but there exists no well-established method to determine the surface structure of a material with a low stability under OER conditions. For materials with a thermodynamically stable bulk, reasonable surface structures can be determined from the surface energy as a function of the applied potential and the pH by means of so-called surface Pourbaix diagrams (SPD) 7 . For thermodynamically unstable materials, on the other hand, one can only investigate the stability of the bulk with respect to competing phases 8,9 , while a SPD would yield the surface with the minimal number of atoms as most stable, indicating that the thermodynamic equilibrium is the fully dissolved state. Nevertheless, in practice the dissolution may be kinetically slow and still allow these materials to be used as catalysts 8 . We refer here to such materials, that are thermodynamically unstable under OER conditions and for which we therefore cannot determine their surface structure by means of SPDs, as metastable materials.
The aim of the present work is to obtain a density functional theory (DFT) description of the OER on surfaces of metastable materials and to reveal how metastability can increase the performance of a catalyst. Computational studies on different oxides have, so far, mainly compared the OER activity of ideal stoichiometric, clean 10 or adsorbate-covered 11 , low-index surfaces which are unlikely to be good approximations for surfaces of a metastable material that may slowly dissolve under OER conditions 8 . Few studies reported the theoretical activity for reaction sites in the vicinity of point defects, steps or kinks 7,12-14 , however without linking such structures to metastability. Our approach to account for metastability is to study the OER on a large variety of different reaction sites by sampling many arXiv:1807.08546v4 [cond-mat.mtrl-sci] 27 Aug 2019 defective surface structures that should be representative for a slowly dissolving metastable material. We study the effect of metastability on the OER for three chemically very different perovskite structured catalysts, which enables us to distinguish between effects resulting from intrinsic material properties and effects of specific reaction sites. We compare LaTiO 2 N (LTON) and SrRuO 3 (SRO), which differ by their d 0 respectively d 4 d-electron occupation and their insulating respectively metallic character. LTON is the best studied member of the family of perovskite oxynitrides that by virtue of their smaller band gap than oxides are considered promising photocatalysts 15 . Calculations 16 and experimentally observed nitrogen loss 17 indicate that LTON is metastable under OER conditions and SRO was also reported to be metastable 5,8 . As these two materials differ not only in their metallic respectively insulating behaviour but also in their cation composition, we complement the set of materials with SrTiO 3 (STO), which shares properties with both LTON and SRO: like SRO it is a pure oxide but like LTON it is a d 0 insulator. Moreover, it shares the same A and B site with SRO and LTON respectively but unlike the other two materials STO was not reported to be metastable.

COMPUTATIONAL APPROACH
For all three materials we investigate both the AO and BO 2 terminated (001) surface, assuming that for a metastable or unstable material both terminations will be exposed as dissolution proceeds. For the oxynitride we also have to consider the anion order, which is different in the surface and the bulk: LaTiO 2 N assumes a cis anion order in the bulk to maximize the overlap between the N 2p and the Ti 3d orbitals 18 but the (001) surface prefers a non-polar trans anion order 19 . We therefore perform calculations on the trans surface but to account for the metastability-induced dissolution of the thin trans layer 20 also investigate the OER on a cis ordered surface. As we find no significant difference between the two anion orders regarding OER reaction free energies, we do however not distinguish between them in our results.
Our stoichiometric and clean surface slabs contain four surface cation sites and we create defective surfaces with between one and four nitrogen, titanium, ruthenium, lanthanum or strontium vacancies in the surface layer (see Fig. 1a for an example). We do not explicitly consider oxygen vacancies, since these are unlikely to exist on the surface under oxidizing OER conditions 11,21 . Occasionally, during surface relaxation -in particular for high cation-vacancy concentrations -either O 2 or N 2 recombine and desorb. In these cases, we stop the calculation and remove the desorbed molecules before continuing the calculation, which can in effect lead to oxygen vacancies. Since potentials relevant for (photo)electrochemical water splitting often lead to formation of oxidising surface species 11 , we also cover surface metal sites with either zero or one monolayer of OH or O in top position in the case of Ti and Ru as well as in the bridge position for La and Sr (see Fig. 1a, as well as ESI † Figs. S1 and S2 for examples of different surface structures). For these surfaces we calculate the OER free energy profile on symmetrically inequivalent titanium, ruthenium, lanthanum and strontium reaction sites in the topmost formula-unit layer -meaning that for example on a defective BO 2 terminated surface also an A atom could become the reaction site.
We consider here only the conventional OER mechanism with four PCET steps (see equation 1-4) that is depicted in Fig. 1b) as we want to study the influence of different reaction sites irrespective of the reaction mechanism. This conventional mechanism does normally not include oxygen ions that are part of the slab. However, in the vicinity of defects and adsorbates it is not always clear if an oxygen originated as part of the adsorbate or the slab and whether the conventional mechanism depicted in Fig. 1b can be applied or not. Therefore, we also include the OER mechanism shown in Fig. 1c which is mathematically equivalent to the conventional mechanism shown in Fig. 1b but operates on an oxygen deficient surface. This so-called lattice oxygen evolution 22 was experimentally observed [23][24][25] and also explained by basic thermodynamic concepts 26 . Alternative lattice oxygen evolution mechanisms have been proposed 27,28 , which are however mathematically different from the conventional mechanism (equations 1-4) and were therefore not included in this study.
In total we calculate the OER on 770 symmetrically inequivalent reaction sites (327 trans LTON, 135 cis LTON, 176 STO, 132 SRO) on 106 LTON, 35 STO and 39 SRO defective (001) surface models. The larger number of reaction sites on the oxynitride stems from anioninduced symmetry breaking. For roughly two thirds of the reaction sites either the *OOH or *OH intermediate was not stable and we exclude these cases from our analysis.
We calculate the change in free energy (∆G) of the four reaction steps (1)-(4) using the computational hydrogen electrode (CHE) 29 , where the energy of a proton and an electron equals half the energy of a hydrogen molecule. As the theoretical overpotential is not dependent on the pH or the potential 2 we perform our calculations at standard conditions (pH=0, T = 298.15 K) and U = 0 V. Zero-point energies (ZPE) and entropies (S) of the reaction intermediates were included as detailed elsewhere 30 . In contrast to other studies, where ZPE was calculated for reaction intermediates in different environments (bridge vs. top site) 30 , we always use the ZPE at the top site. Given that for defective surfaces a full ZPE evaluation would be computationally prohibitively expensive and that previously reported changes in ZPE were minor 30 , this approach will yield correct trends for defective surfaces.
We estimate the activity of a specific reaction site by the largest step in its OER free energy profile: ∆G 0 1−4 being the change in free energy of the four OER steps (equations 1-4). The calculated overpotential is given by: where 1.23 V is the potential needed to make all ∆Gs equal to zero for an ideal catalyst. The adsorption free energies of the intermediate species ∆G O , ∆G OH and ∆G OOH were calculated as follows: where free energies include changes in ZPE and S while n and m are stoichiometric coefficients that preserve the number of atoms on both sides of the respective reaction. We further compare the calculated overpotential η OER with the overpotential η OER D predicted by the single descriptor ∆G 0 2 : It was established that there exists a universal scaling relation between the adsorption energies of the reaction intermediates *OH and *OOH, their difference being approximatively 3.2 eV for metal and oxide surfaces irrespective of the reaction site 31 . Since the overpotential for oxides and metals is generally determined by either step 2 (∆G 0 2 ) or 3 (∆G 0 3 = 3.2eV − ∆G 0 2 ), the former is often a suitable descriptor of the OER activity that determines the overpotential 2 : where ∆G 0 2 is obtained via the difference in adsorption energies.

RESULTS AND DISCUSSION
We begin our analysis by comparing the calculated overpotentials η OER with the overpotentials determined via the single descriptor η OER D and find in general a good agreement for the majority of the calculated overpotentials (see Fig. 2a). Most remarkable is the fact that we find for all three materials a continuous distribution of overpotentials within a large range of the descriptor ∆G 0 2 , including also the top of the volcano. The complex dissolution kinetics, make it difficult to determine with certainty which sites will be predominantly present on a dissolving surface, yet we can attempt to distinguish between reactive sites occurring more likely (filled circles in Fig. 2a) and less likely (empty circles in Fig. 2a) based on thermodynamic and kinetic considerations: From surface Pourbaix diagrams (see ESI † section S2 for details) we see, in agreement with experimental studies 32,33 , that A sites dissolve more readily and therefore B sites represent likelier OER reactive sites. Moreover, the Bell-Evans-Polanyi principle 34,35 implies that the formation of B-site vacancies is also kinetically less likely than the formation of A-site vacancies. Since under OER conditions surfaces are normally oxidized 11 , we consider B sites oxidized either by an OH or O species to be the thermodynamically and kinetically most likely OER reaction sites. For LTON we consider the trans anion order to be more likely due to electrostatic considerations 19 and for STO we only consider the O covered defect free TiO 2 termination as likely since the material is stable. When considering only the overpotentials of these likely reactive sites we still find a continuous distribution of overpotentials, many of them at the top of the activity volcano. Since many of the less likely (less stable) sites also have overpotentials close to the top of the volcano, this implies, in agreement with recent experiments 6 , that there exists no correlation between the stability of a single site and its catalytic activity and that the activity-stability relations must have a different origin.
Given the large number of likely sites close to the top of the volcano, we can assume that some of these will at least be transiently present. Since we omitted combinations of different defects types and considered only zero or full adsorbate coverage, the actual variety of reaction sites should be even larger. Together with alternative reaction mechanisms that could be active, we thus expect even more sites with small overpotentials on a metastable catalyst. Since the activity of a catalyst depends exponentially on the overpotential, even a small number of highly-active sites can dominate the apparent activity 14,36,37 . We therefore postulate that in operando, surface dissolution can lead to the appearance of highly active sites, resulting in an enhanced apparent activity of the catalyst. In other words, the less a material is stable, the higher should be the diversity of reaction sites and the higher the chance of finding highly active reaction sites, which we believe to be the origin of the activity-stability relations. We stress that at this point the proposed activity improvement requires the material to be metastable under the chosen conditions. However, metastable mate-rials such as LTON and SRO can of course only benefit from this proposed activity enhancement if the dissolution is slow enough for the material to still sustain extended periods of operation.
While previous studies [3][4][5] suggest that active dissolution of the material results in the best activities, our results suggest that an activity enhancement can also be achieved by an initial dissolution (preconditioning) of the catalyst that creates the variety of reaction sites but that the dissolution does not have to continue for high OER activities. This agrees with recent experimental findings for RuO 2 thin films, which show that while the surface orientation with the initially highest dissolution rate results in the best activities, an ongoing dissolution is not necessary for highly active surfaces 6 .
The similarity of the overpotential distribution of LTON, SRO and STO implies further that the high activity does not originate purely from the bulk electronic structure of the catalyst as is often assumed, but that reaction sites with geometries that enable a lower overpotential can form on any of the investigated materials. This is supported by the experimental finding that perovskite oxide electrocatalysts with structural flexibility can lead to superior activity 38 . We were not able to correlate catalytic activity with geometrical descriptors or simple structural features such as the reaction site species, angles and bond lengths of reaction intermediates or the adsorbate coverages to mention only a few. We also find no correlation with global properties such as the dipole moment of the slab or the nom-inal charge of the atoms. However, in agreement with other studies 13,39,40 , a preliminary analysis that will be reported elsewhere 41 hints at a correlation of ∆G 0 2 with the electronic structure and in particular the centre of the O2p band of the oxygen adsorbate.
We continue by predicting the minimum overpotentials of the three materials. From the relation between the adsorption energies of the *OH and *OOH intermediates, we predict minimum overpotentials of 0.23 V for STO, 0.26 V for SRO and 0.40 V for LTON, the two former being smaller than the 0.37 V predicted by the scaling relations, while the latter is larger. Despite these differences, we do for STO and SRO not find actual reaction sites with overpotentials smaller than 0.37 V. This is in line with the observation that breaking the universal scaling relation between *OOH and *OH is a necessary but not sufficient condition to optimize the OER 42 . Further, while there are only small differences between the two oxides that show a rather good agreement between the calculated η OER and the single-descriptor overpotentials η OER D , for LTON we find a not insignificant number of reaction sites that show large deviations from η OER D .
A more detailed analysis of the reaction free energies of the four charge transfer steps ∆G 0 i for LTON (see Fig.  2b), reveals that while for most reactions ∆G 0 2 or ∆G 0 3 (blue and red data points) have the largest free-energy change, there are some reactions where ∆G 0 4 (green data points) has the largest free energy change and is the limiting step. This implies that for the oxynitride the step from *OOH to O 2 is energetically less favourable while at the same time the formation of *OH is more favourable. This difference with respect to oxides 2 suggests a stronger adsorption of *OOH on the oxynitride and leads to larger calculated overpotentials compared to the single-descriptor overpotentials η OER D for some of the OERs.
From the volcano plot in Fig. 2a we can see a further difference between the two oxides and the oxynitride: While for the oxynitride the descriptor value ∆G 0 2 ranges from -2 eV to 6 eV, the two oxides seem to have a much smaller range of ∆G 0 2 and their overpotentials are therefore more concentrated around the top of the volcano. The spread in ∆G 0 2 is important for the apparent activity. A catalyst with an activity of the perfect surface far from the top of the volcano would benefit from a large spread to increase the number of highly active sites. An already good catalyst located not far from the top of the volcano, on the other hand, would benefit from a small spread to have a large proportion of highly active sites. We can explain the difference in the ∆G 0 2 range between oxides and oxynitrides by analysing the charge transfer towards different reaction intermediates. To do so, we look at the change of average Löwdin charges 43 of the anions for different reaction steps. We find the best correlations of the Löwdin charges with ∆G 0 2 if we include all anions and not only surface atoms. We compute differences between element-averaged Löwdin charges for the different reaction steps: where a is an anion element (here O or N), i and j are two consecutive reaction intermediates, N i/j a is the number of a atoms in the surface structure with the i/j intermediate and l i/j x is the Löwdin charge for atom x in the structure with intermediate i/j, while the sum runs over all atoms of element a.
As we show in Fig. 3a, there is a direct correlation between the descriptor ∆G  Fig. 3b we find an inverse correlation between these charge differences, suggesting a charge transfer from N to O. From the colour coding of the data points, we see that ∆G 0 2 is generally smaller when nitrogen loses more electrons (large positive ∆L * OH− * O N ), while the charge variation on the oxygen atoms is small. From these observations we propose that nitrogen can act as an electron reservoir: The reaction from *OH to *O is favoured if nitrogen ions provide electrons, resulting in a lower descriptor value ∆G 0 2 , while reaction sites, where this charge transfer is not possible, result in larger ∆G 0 2 s. More generally, this implies that a more flexible valence-band electronic structure, for example in mixed-anion materials, can yield a larger variety of ∆G 0 2 s.

CONCLUSIONS
In summary, we have shown for three electronically different perovskite materials that the structural variety of reaction sites on dissolving heterogeneous catalysts will, for some sites, lead to overpotentials close to the top of the activity volcano. Given that reaction sites with low overpotentials will dominate the apparent activity, this implies that, independent of the bulk electronic properties, operating a heterogeneous catalyst at the border of its stability region can lead to an increased apparent activity due to the increased variety of geometrically different reaction sites. Also preconditioning catalysts under metastable conditions before operation under stable conditions should result in the same activity enhancement. While we find the single-descriptor approach to generally work well for these highly defective surfaces, we observe that the strong binding of the *OOH and *OH intermediates on LaTiO 2 N can lead to different limiting steps and larger deviations from the descriptor-based overpotential than for oxides. We find only minor differences in activity between metallic SrRuO 3 and insulating SrTiO 3 but show that a flexible valence band structure as it occurs in LaTiO 2 N yields a larger variety of ∆G 0 2 values.

METHODS
We determine energies of the various adsorbate covered, defective surfaces by density functional theory (DFT) calculations using the Quantum ESPRESSO 44 package at the GGA+U level of theory with the PBE 45 exchange-correlation functional and a Hubbard U 46 of 3 eV applied to the titanium 3d states.
We do not apply a Hubbard U on ruthenium 4d orbitals as this setup results in the best agreements with experimentally measured magnetic moments and the density of states at the Fermi energy. We use ultrasoft pseudopotentials 47 with La(5s,5p,5d,6s), Sr(4s,4p,5s), Ti(3s,3p,3d,4s), Ru(4d,5s,5p), O(2s,2p) and N(2s,2p) electrons as valence states to describe the interaction between electrons and nuclei and perform spin-polarized calculations in the case of SrRuO 3 . The cutoff for the plane-wave basis set is 40 Ry and 320 Ry for the kinetic energy and the augmented density respectively for LaTiO 2 N and SrTiO 3 , while for SrRuO 3 we use slightly higher cutoffs of 50 Ry and 500 Ry respectively. We start our calculations from a 40-atom pseudo-cubic perovskite cell and create asymmetric 1 × 1 × 2 (001) surface slabs (see Fig. 1a) with at least 10Å vacuum, two fixed atomic layers at the bottom of the slab and a dipole correction in the vacuum layer 48 . The Brillouin zone is sampled with a 4 × 4 × 1 Monkhorst-Pack 49 k-point grid for LaTiO 2 N and SrTiO 3 and a 6 × 6 × 1 grid for SrRuO 3 . We relax ionic positions until forces converge below 0.05 eVÅ −1 and total energies change by less than 1.4 · 10 −5 eV.

Supplementary information for
The role of metastability in enhancing water-oxidation activity

S2. OCCURRENCE OF SURFACE STRUCTURES -PARTIAL SURFACE POURBAIX DIAGRAMS
To qualitatively assess the likelihood for different defective surfaces to occur, we calculate partial surface Pourbaix diagrams meaning that we analyse the potential and pH dependent occurrence of different vacancy types independently of each other. A conventional surface Pourbaix diagram would only show that our metastable materials are not thermodynamically stable for relevant potentials and pHs. However, by comparing the stability of surfaces with single defects to the perfect stoichiometric surface we obtain qualitative information as to which surface structures are more likely to be present than others.
For these calculations we consider defect formation accompanied by formation of solvated ionic species. Since the energy of these solvated ionic species is difficult to calculate within DFT, we use the approach of Persson et al. 1 , where the chemical potentials are calculated with respect to experimental reference energies, allowing to combine DFT energies and Gibbs free energies of solvated ions. We also apply a correction to the energy of the oxygen atom to correct for the well-known overbinding of O 2 within DFT. We obtain the correction of 0.68 eV via fitting the experimental formation enthalpies of metal oxides (no transition-metal oxides) to our DFT energies 2 . Further, DFT+U calculations impose challenges in calculating formation enthalpies of reactions involving compounds with both localized and delocalized electronic states. Therefore, we also add a correction to every element to which a Hubbard U was applied 3 such as to correctly reproduce the formation enthalpy of TiO 2 (for Ru no correction is applied as we use no Hubbard U in our calculations). As we do not include surface adsorbates in this analysis, we do not include zero-point energies of adsorbates. Finally, to be able to determine the partial surface Pourbaix diagram with a numerical solver, we prevent the dissolution of the surface by imposing the presence of atoms that are part of the so-called reference set. Figure S3 shows the obtained partial surface Pourbaix diagrams where we consider only the stoichiometric perfect surface (stoi) and a single vacancy. We do not consider oxygen vacancies as these are likely to be spontaneously healed by O* adsorbates under OER conditions 4 . From Figures S3 a)-c), we see that at a given pH La vacancies (depending on the pH via formation of various solvated species) on the LaTiO 2 N surface start appearing already at potentials below the equilibrium potential for water oxidation. Nitrogen vacancies start forming slightly above the water oxidation potential, whereas Ti vacancies require the highest potential to form. A similar trend is observed for SrRuO 3 in Figures S3 d)-e), where Sr vacancies are formed at significantly lower potentials than Ru vacancies. From this analysis it thus seems that for our set of materials A-site vacancy formation occurs at smaller potentials than B-site vacancy formation, while the potential for N vacancy formation is similar to the B-site. This agrees with experimental observations that the B-site in perovskite oxides is generally thermodynamically more stable than the A-site in contact with an aqueous solvent 5,6 .
The thermodynamic stability is generally correlated with the kinetic stability via the Bell-Evans-Polanyi principle 7,8 , which states that kinetic activation barriers (E a ) obey the relation E a = E 0 + α∆H, where E 0 is a reference energy for a given class of reactions, ∆H is the thermodynamic reaction enthalpy and α characterises the position of the transition state along the reaction coordinate. This principle implies that reactions with a larger thermodynamic energy difference (∆H) will also have a larger activation barrier. Applied to our dissolving surfaces, this implies that the formation of B-site vacancies is both thermodynamically and kinetically less likely than the formation of A-site vacancies. We thus consider that surfaces with B reaction sites are more likely to exist on both thermodynamic and kinetic grounds than surfaces with A reaction sites.