Kinetics-Based Computational Catalyst Design Strategy for the

nical details of the calculations are specified in the Supporting Information. ... in order to encourage the system to converge to the desired electro...
0 downloads 0 Views 5MB Size
Subscriber access provided by Kaohsiung Medical University

Article

Kinetics-Based Computational Catalyst Design Strategy for the Oxygen Evolution Reaction on Transition Metal Oxide Surfaces Craig Patrick Plaisance, Simeon Dominikus Beinlich, and Karsten Reuter J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b08549 • Publication Date (Web): 01 Nov 2018 Downloaded from http://pubs.acs.org on November 17, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Kinetics-based Computational Catalyst Design Strategy for the Oxygen Evolution Reaction on Transition Metal Oxide Surfaces Craig P. Plaisance,∗,†,‡ Simeon D. Beinlich,† and Karsten Reuter† †Chair for Theoretical Chemistry and Catalysis Research Center, Technische Universit¨at M¨ unchen, Lichtenbergstr. 4, D-85747 Garching, Germany ‡Present Address: Cain Department of Chemical Engineering, Louisiana State University, Baton Rouge, LA 70803, USA E-mail: [email protected]

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Density functional theory was used to examine the oxygen evolution reaction on a large number of active sites formed by doping three different surfaces of Co3 O4 with various 3d transition metal atoms. By combining the scaling and BrønstedEvans-Polanyi (BEP) relationships that are determined for these sites, it is shown that the activity of a site is controlled by the redox potential for oxidation of the site. Based on this, a kinetics-based design strategy is presented for identifying the optimal active site at a given electrode potential. This design strategy is shown to be valid regardless of whether the rate limiting water addition step occurs electrochemically or non-electrochemically, as long as certain conditions are met. Another finding is that the BEP relations are sensitive to the structure of the active site, with sites reacting through μ3 -oxo species having the most favorable relations. Finally, the kinetics-based design strategy is compared with the commonly used thermodynamics-based design strategy and it is shown that the former is able to identify a site with equal or greater activity than the latter at the same computational cost.

Introduction Water electrolysis and CO2 reduction are two of the most promising means of storing energy from renewable but intermittent energy sources such as solar and wind. 1 Both of these processes involve the oxygen evolution reaction (OER) at the anode, which oxidizes water to molecular oxygen in order to supply the protons and electrons required to produce H2 or reduce CO2 at the cathode. However, the sluggish kinetics of the OER require high overpotentials in order to obtain a sufficient current density (0.30–1.00 V for 1 A/cm2 2 ), resulting in a significant loss of efficiency in both processes. 3–5 Additionally, the catalysts at the lower end of this overpotential range are those that operate in an acidic electrolyzer and are composed of expensive IrO2 . Catalysts composed of the cheaper 3d transition metal oxides operating in an alkaline electrolyzer require overpotentials of 0.60–1.00 V to achieve a 2

ACS Paragon Plus Environment

Page 2 of 49

Page 3 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

current density of 1 A/cm2 2 (although it is predicted that an overpotential of 0.35 V would be sufficient if the thickness of the separator between the anode and cathode chambers could be reduced 6 ). Large-scale deployment of water electrolysis and CO2 reduction would at the minimum require an OER catalyst that can produce sufficient current densities at low overpotentials and which is composed of either low-cost metals (e.g. Mn, Fe, Co, Ni) or very small quantities of expensive metals (e.g. Ru, Ir). 3 It is a complex question as to which of these two approaches is more likely to lead to a commercially viable electrolyzer. 6 Here, we choose to focus on the first of these two approaches. There do appear to be several 3d transition metal oxide OER catalysts giving high turnover frequencies (TOF) at reasonable overpotentials. 3–5 One of the highest identified so far is Fe doped amorphous Ni oxide, which gives a TOF of 0.2 /site/s at an overpotential η of 0.3 V. 7 Nonetheless, it would still be desirable to identify catalysts that could achieve this same TOF at a lower overpotential. For instance, reducing the overpotential by half would increase the efficiency from 80% to 90% (calculated as 1.23/(1.23 + η)). More importantly, operating at lower overpotentials could increase the long-term stability of the catalyst. An additional benefit to identifying new OER catalysts is that it would increase the materials base for such catalysts and therefore avoid the geopolitical complications of relying too heavily on certain materials. In the search for new and improved catalysts, computational catalyst design has played an increasingly important role over the last 15 years by helping to understand how catalysts work at the molecular level in order to guide experimental synthesis. 8–10 The most widely employed computational design strategy for the OER catalyst to date is based on an approach that considers only the thermodynamics of reaction intermediates on the active site, 5,11–13 avoiding the need to calculate more computationally demanding activation barriers. It provides an estimate of the relative activity of a site based on the calculation of a single easily computable descriptor, typically the difference between the binding energies of O and OH to the active

3

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

site. Despite its simplicity, this thermodynamics-based design strategy does appear to predict activities that roughly correlate with experimental measurements. 5,14 However, one of the main assumptions underlying this approach is that the activation barriers of all steps are relatively small when these steps have negative free energies. 11 Previously, we have found that two particular steps in the catalytic cycle on Co3 O4 – the reaction of water with an oxygen on the active site to form a hydroperoxo intermediate, called the water addition reaction, and the elimination of O2 from the active site at the end of the catalytic cycle – have significant activation barriers ranging from 0.3–1.0 eV. 15 This raises concerns as to the applicability of the thermodynamics-based screening approach for this material. We therefore proposed an alternative catalyst design strategy for the OER on 3d oxides that is based on a simplified kinetic model rather than a purely thermodynamic model. 16 The kinetic model assumes that water addition occurring by a non-electrochemical pathway is the rate limiting step and that all electrochemical steps are quasi-equilibrated. By observing that a Brønsted-Evans-Polanyi (BEP) relation exists between the activation barrier for water addition and its reaction free energy and that a linear scaling relation exists between O – H and O – OH bond dissociation energies, we were able to correlate the activity of a given active site with the redox potentials of two pre-equilibrium oxidation steps that occur prior to water addition. This led us to formulate an alternative design principle in that the optimal active site at a given electrode potential has to have both of these redox potentials equal to the electrode potential. Our previous kinetics-based catalyst design strategy was based on a limited set of data, considering only four active sites on different surface terminations of pure Co3 O4 . Validating and extending our approach, we here examine the OER on a large number of active sites formed by doping three different surfaces of Co3 O4 with various 3d transition metal atoms. First, we give an overview of the mechanism and kinetics previously found for the OER on Co3 O4 . Then, we discuss the mechanism of the water addition step and find that the form of the BEP relation we used previously does not adequately represent the larger set

4

ACS Paragon Plus Environment

Page 4 of 49

Page 5 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

of active sites. We identify a more suitable form for the BEP relations which suggests that electrochemical and non-electrochemical water addition pathways are kinetically equivalent on most active sites. After discussing the scaling relations between O – H and O – OH bond dissociation energies obtained for the larger set of active sites, this allows us to formulate a new set of design criteria for the optimal active site. Finally, we contrast our kineticsbased approach with the thermodynamics-based approach to better understand the range of validity of the latter.

Computational Methods All electronic structure calculations mentioned in this study were performed using spinpolarized density functional theory (DFT) within the generalized-gradient approximation (GGA) 17 as implemented in the Vienna Ab-initio Simulation Package (VASP). 18 The technical details of the calculations are specified in the Supporting Information. To correct for the over-delocalization of the 3d electrons caused by the electron self-interaction present in GGA functionals, 19 we use the GGA+U method introduced by Anisimov and coworkers. 20–22 The value of U employed was determined with a self-consistent linear-response approach, described in more detail in our previous work. 15 Most of the transition states were determined using the nudged elastic band (NEB) method 23,24 in combination with the dimer method. 25 Due to the higher computational cost of the NEB method, it was first used to roughly converge the minimum energy pathway between the reactant and product state in order to obtain an initial guess for the transition state. This initial guess was then refined to a higher level of convergence using the dimer method. A few of the transition states could not be found using the dimer method because they corresponded to a crossing between two potential energy surfaces rather than a saddle point on a single potential energy surface. These transition states were refined using an intersystem crossing algorithm 26 that we implemented in VASP that is able to locate the

5

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

minimum energy crossing between two potential energy surfaces. Details of this algorithm and further details of the transition state search procedure are given in the Supporting Information. A particularly important issue encountered during DFT calculations on systems with localized electrons is that there are often several meta-stable electronic states. 27 For the 3d transition metal oxides, these different electronic states correspond to different spin states on the 3d metal centers. It is therefore important to ensure that the correct electronic state is obtained in each calculation. This is done using two tools that we have developed. First, in order to encourage the system to converge to the desired electronic state, we initially run a calculation in which the magnetic moments on the transition metal atoms are constrained to particular values associated with that electronic state. This is done using a feature we implemented in VASP that is described in more detail in the Supporting Information. We then use the resulting geometry and electronic structure as an initial guess for an unconstrained calculation. This approach is similar in spirit to the occupation matrix control method developed by Watson and coworkers, 27 but is simpler to operate and usually converges to the correct electronic state. In cases where this approach fails, we resorted to the more powerful occupation matrix control method. We also checked the converged electronic structure of each calculation using the quasiatomic orbital method of Qian et al., 28 which we have also implemented in VASP. 29 This approach transforms the plane wave basis set to a minimal basis of atomic orbitals called quasiatomic orbitals. We can then obtain the orbital population matrix in this basis for each atom and use a procedure outlined in the Supporting Information to confirm that the correct electronic configuration has been obtained. In order to calculate free energy differences for electrochemical steps, we employ the computational hydrogen electrode approach developed by Rossmeisl, Nørskov, and coworkers. 11 For steps in which a molecule of water is consumed, we use hexagonal ice as the reference state when calculating the reaction free energy. To reduce computational effort, zero-point vibrational energy and thermal vibrational contributions to the free energy are neglected

6

ACS Paragon Plus Environment

Page 6 of 49

Page 7 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

for surface calculations. Further details and justification of this approach are given in the Supporting Information.

Results and Discussion Mechanism and Kinetics of the OER on Co3 O4 In our previous work 15 we determined that the OER on Co3 O4 proceeds on an active site consisting of two redox centers (M – O)a and (M – O)b that are responsible for oxidizing water. Each center consists of an octahedrally coordinated 3d metal cation and a protonatable oxygen atom that is accessible to the electrolyte. This two-center site was found to carry out the OER on several different surface terminations of Co3 O4 by the generalized catalytic cycle depicted in Scheme 1. It begins by the oxidation of the (M – O)a center to form a reactive oxo species in an electrochemical step, followed by a second electrochemical step in which the (M – O)b center is also oxidized, n+1 − −− (Mn −OH)a + OH− ) −* − (M −O)a + e

,

{1}

m+1 (Mm −OH)b + OH− − −O)b + e− )− −* − (M

.

{2}

In the next step, an O – O bond is formed between Oa (the oxygen of (M – O)a ) and a molecule of water from the electrolyte to form a hydroperoxo intermediate. This step, called water nucleophilic addition or the acid/base mechanism in the literature 30 (or simply water addition), can occur either by a non-electrochemical pathway 3 or by an electrochemical pathway 30 (the latter pathway bypasses step 2),

(Mn+1 −O)a + (Mm+1 −O)b + H2 O −−→ (Mn −OOH)a + (Mm −OH)b (Mn+1 −O)a + OH− −−→ (Mn −OOH)a + e− 7

ACS Paragon Plus Environment

.

,

{3} {30 }

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 49

We will see later that these two pathways should be kinetically equivalent on most active sites. Following water addition, the hydroperoxo intermediate is electrochemically oxidized to a superoxo and the (M – O)b center is again electrochemically oxidized, n − (Mn −OOH)a + OH− − )− −* − (M −O2 )a + e

{4}

,

m+1 −− −O)b + e− (Mm −OH)b + OH− ) −* − (M

{5}

.

Finally, the superoxo transfers an electron to Mb (the metal cation in (M – O)b ) to eliminate O2 in a non-electrochemical step, followed by dissociation of a second molecule of water to complete the catalytic cycle,

(Mn −O2 )a + (Mm+1 −O)b + H2 O −−→ (Mn −OH)a + (Mm −OH)b + O2

.

{6}

Alternatively, it may also be possible for the O2 elimination step to occur electrochemically and bypass step 5,

(Mn −O2 )a + OH− −−→ (Mn −OH)a + O2 + e−

.

{60 }

Although the exact details of each step were found to vary between the different surface terminations of Co3 O4 that were examined, this general catalytic cycle remained unchanged. In our previous work 15,16 as well as the current work, we assume that only the water addition (3 and 30 ) and O2 elimination (6) steps are kinetically relevant, while the electrochemical oxidation steps are fast and quasi-equilibrated. This is justified by the fact that these latter steps only involve simple proton transfer, with no other bonds being formed or broken. Calculations have shown that the free energy barriers for such reactions are 0.1– 0.2 eV in the thermodynamically favorable direction, 31–33 significantly less than the barriers

8

ACS Paragon Plus Environment

Page 9 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

H2O + e−

OH−

O=O

OH−

H2O

6

H2O + e−

OH−

6'

e−

OH−

5 H2O + e−

1

e−

O=O

3'

2 3

4

H2O

OH− OH−

H2O + e−

Scheme 1: Generic catalytic cycle for the the oxygen evolution reaction on a transition metal oxide surface. Numbers in the figure correspond to the reaction step equations in the main text. Steps assumed to be quasi-equilibrated are depicted with a circle through the reaction arrows. of 0.3–1.0 eV we have found for the water addition step. 15,16 We also found that in most circumstances the water addition step is rate limiting so that the O2 elimination step could be neglected. 15 In this case, the kinetics are governed by the two quasi-equilibrated electrochemical steps 1 and 2, followed by the rate limiting water addition step 3 to form OOH (we will consider electrochemical water addition by step 30 later). With this simplification, the turnover frequency (TOF) can be written as a function of the electrode potential η, 15 k3    ,   e(η − η1 )  e(η − η2 ) 1 + exp − 1 + exp −  kB T kB T

TOF = 



(1)

where k3 is the intrinsic rate constant for the water addition step, while η1 and η2 are the equilibrium potentials for the two redox steps. In deriving this rate expression (shown in the Supporting Information), it is assumed that the site can exist in only four states – both (M – O)a and (M – O)b centers reduced, one reduced and the other oxidized, and both 9

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

oxidized. In all four of these states, the metal cations and oxygen anions composing the remainder of the surface in which the active site is embedded are held in fixed oxidation and protonation states. In reality, the oxidation and protonation states of these surrounding atoms will change with the electrode potential, leading to a more complex dependence of the TOF on the electrode potential as seen in our previous work. 15 As we will discuss in the final section of the manuscript, however, these modifications to the TOF should not change the design principles that we eventually formulate. The ideal OER catalyst should be able to achieve a high TOF while operating at low electrode overpotentials, thus maximizing both current density and efficiency. Our previous work, 16 however, indicated that it may not be possible to optimize both of these characteristics simultaneously. It was observed that a linear relation exists between the logarithm of the intrinsic rate constant k3 for water addition and the redox potentials η1 and η2 of the pre-equilibrium electrochemical oxidation steps. A detailed quantum chemical explanation for this relationship was given, essentially that the O – H bond broken during the oxidation step and the O – OH bond formed during water addition are both controlled by the energy to localize a hole on the reactive oxo of the (M – O)a center, leading to a linear scaling relation between the energies of the two processes. This scaling relation, combined with a linear BEP relation between the activation barrier and reaction energy for water addition, means that we cannot have a catalyst with a low barrier that operates at a low overpotential. As mentioned in the introduction, the relationship just described was based on a limited set of data, consisting of active sites on only four different surface terminations of undoped Co3 O4 . In this work, we now examine the water addition step on 48 different active sites constructed by doping three of these four surface terminations of Co3 O4 with different combinations of 3d transition metal cations. As illustrated in Figure 1, the resulting active site geometries differ in the coordination number of the oxygen Oa that forms the O – O bond. A site where Oa is coordinated to three metal cations (μ3 -Oa ) exists on the (311) surface, a site where Oa is coordinated to two metal cations (μ2 -Oa ) exists on the (110)-A surface, and

10

ACS Paragon Plus Environment

Page 10 of 49

Page 11 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

a)

b) Ob

Ob

Oa Ma

Ma Oa

Ma

Oa

Mb

Mb

Ob

Mb

c)

(311) µ³-oxo

(110)-A µ²-oxo

(110)-B η-oxo

Figure 1: Active site geometries on the (311), (110)-A, and (110)-B surfaces of Co3 O4 showing (a) the top-down view of the surface (with the active site cluster indicated with darker colored atoms) and the active site cluster in the (b) reactant state and (c) transition state for water addition. Oxidized metal centers are indicated in light green for the +4 oxidation state and dark green for the +5 oxidation state.

11

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

a site where Oa is a terminal oxo (η-Oa ) exists on the (110)-B surface. The (001) surface that we examined previously was not considered in this study because it also has a μ3 -Oa and is thus similar to the (311) surface. The active site on each of the three surfaces features the two redox-active metal cations (Ma and Mb ) that are central to the above described OER mechanism. The doped active sites were constructed by substituting different 3d transition metal cations (V, Cr, Mn, Fe, Co, and Ni) at these two redox-active centers. In the next section, we discuss the mechanism for the water addition step on these doped active sites and show that a qualitatively different BEP relation is more suitable than the one we obtained previously.

Mechanism and Activation Barriers for the Water Addition Step As the water addition step is likely the rate limiting step in the OER catalytic cycle, it should lie at the core of any design strategy for the active site. Specifically, it is necessary to understand how the activation barrier of this step is related to the thermodynamic properties of the active site. In our previous work, 16 we assumed that this is described by a standard BEP relation in which the activation barrier is linearly correlated to the reaction free energy of water addition. However, with the more extensive set of active sites examined in the present work, we find that a qualitatively different type of BEP relation is better suited for expressing the activation barrier. First, we discuss the mechanism of water addition on the doped active sites examined in this study. It was found to be essentially the same as the mechanism we found previously on undoped Co3 O4 , 16 shown in Scheme 2. The reaction begins by transfer of a proton from a water molecule to the Ob atom of the active site, accompanied by transfer of an electron to an empty 3d orbital on Mb . This leads to an intermediate state consisting of a hydroxyl radical, which may or may not be a local minimum on the potential energy surface. The hydroxyl radical then approaches the Oa atom, and at some point a second electron is transferred from the antibonding orbital of the nascent O – O bond into an empty 3d orbital on Ma . 12

ACS Paragon Plus Environment

Page 12 of 49

Page 13 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Following this, the O – O bond continues to form until the reaction step is complete. R

O−O bond formation H+/e− transfer TS

OH−

H2O + e−

I

P



Scheme 2: Mechanism for the water addition reaction by non-electrochemical (top) and electrochemical (bottom) pathways, consisting of quasi-equilibrated H+ /e – transfer steps followed by a common O – O bond formation step. The H+ /e – transfer steps occur between a reactant state (R or R’) and a hydroxyl radical intermediate (I) that is common to both pathways. The kinetically relevant transition state common to both pathways occurs during the O – O bond formation step in which the intermediate goes to the product state (P). Since the water addition reaction involves two steps, H+ /e – transfer and O – O bond formation, there can in principle be two transition states. We found this to be the case for most of the sites on the (110)-A surface, although the O – O bond formation transition state was found to be the higher in energy of the two for all sites where both barriers were calculated. On all sites examined on the (311) surface, only a single transition state was found, corresponding to the O – O bond formation step. A single transition state was also found for all sites examined on the (110)-B surface, however in most cases it corresponded to H+ /e – transfer and O – O bond formation occurring concertedly. In contrast to what we found in our previous work on undoped Co3 O4 , a single BEP relation of the form ∆G‡3 = α∆G3 + β

,

(2)

is inadequate to capture the changes in activation barrier ∆G‡3 between the different sites (in 13

ACS Paragon Plus Environment

The Journal of Physical Chemistry

terms of the reaction free energy ∆G3 ), as can be seen in the left panel of Figure 2. Using a different BEP relation (different values of α and β) for each of the active site geometries ((311), (110)-A, and (110)-B) yields a better fit, however there are still many sites where deviations greater than 0.1 eV exist between the calculated and BEP barriers. The reason that the activation barrier is not well described by a BEP relation, at least on the (311) and (110)-A sites, is related to the existence of the intermediate state along the reaction trajectory. This intermediate arises because the H+ /e – transfer and O – O bond formation occur sequentially rather than concertedly on these sites. As the transition state for O – O bond formation occurs between the intermediate and product states, there is no reason to expect that the barrier should be well correlated with the free energy difference between the reactant and product states. Rather, a BEP relation should exist between the activation barrier and reaction free energy with respect to the intermediate state. 2.4

2.0 (311) (110)-A top (110)-A bottom (110)-B

2.2 2.0 1.8

1.8

1.6 1.4

α ������������� β ������������� eV

α �������������

∆G‡3' (eV)

∆G‡3 (eV)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 49

α �������������

β ������������� eV

β ������������� eV

1.2

1.6

1.4

1.0 0.8

1.2

0.4

α �������������

α �������������

0.6

β ������������� eV

β ������������� eV

-1.0

-0.5

0.0

0.5

1.0

1.0

1.5

∆G3 (eV)

0.5

1.0

1.5

∆G3' (eV)

Figure 2: (left) Standard BEP relations for the non-electrochemical water addition pathway and (right) modified BEP relations corresponding to the electrochemical water addition pathway. The sites on the (110)-A surface are divided into two sets (top and bottom) depending on the location of the metal center that gets reduced first (explained in the Supporting Information). Note that the y-axis range in the left panel is twice the range in the right panel, i.e. the deviations are larger than they seem compared to those in the right panel.

14

ACS Paragon Plus Environment

Page 15 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Unfortunately, the intermediate state is not a defined local potential energy minimum on the (311) sites so that it is not possible to calculate its free energy on the basis of the static NEB calculations. However, as the intermediate consists of a hydroxyl radical having only non-bonded interactions with the active site, it should be reasonable to use in its place a surrogate structure in which the hydroxyl radical is in vacuum. This is equivalent to splitting up the water addition reaction into the following two components:

(Mm+1 −O)b + H2 O −−→ (Mm −OH)b + OH , (Mn+1 −O)a + OH −−→ (Mn −OOH)a

.

{3a} {3b}

The first step consists of the electrochemical oxidation of liquid water to gas phase hydroxyl (described by the redox potential ηH2 O , calculated to be 2.23 V with respect to the OER equilibrium potential) coupled with the electrochemical reduction of the active site by the reverse of reaction 2 (with redox potential η2 ). It therefore has a reaction free energy ∆G3a of ∆G3a = e(ηH2 O − η2 ) .

(3)

The second step involves the formation of an O – O bond between the hydroxyl radical and the active site, having a reaction free energy ∆G3b equal to the negative of the dissociation free energy of the O – OH bond,

∆G3b = −∆GO−OH

.

(4)

It is this second step that corresponds to the intermediate-to-product part of the reaction pathway on which the BEP relation should be based. This BEP relation correlates the activation barrier ∆G‡3b for step 3b with the free energy of this step – although for reasons that will become apparent in the next section, a more physically meaningful expression is

15

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 49

obtained if we shift the origin of the BEP relation to lower energies on both axes by the quantity eηH2 O , ∆G‡3b + eηH2 O = α(∆G3b + eηH2 O ) + β

.

(5)

Defining ∆G‡30 = ∆G‡3b + eηH2 O and ∆G30 = ∆G3b + eηH2 O , this becomes ∆G‡30 = α∆G30 + β

.

(6)

The activation barrier with respect to the reactant state is finally obtained by adding the free energy of step 3a in eq (3) to ∆G‡3b , leading to a modified BEP relation for the water addition reaction, ∆G‡3 = α∆G30 + β − eη2

.

(7)

The BEP relations obtained from fitting eq (6) to the energetics calculated for the doped sites on the (311) and (110)-A surfaces are shown in the right panel of Figure 2. We can see that the activation barrier for sites on the (311) surface is well described by this modified BEP relation, with a maximum deviation of 0.06 eV, significantly smaller that for the original BEP relation (0.13 eV). The fit is less satisfactory for sites on the (110)-A surface, although the maximum deviation (0.11 eV) is less than when using the original BEP relation in eq (2) (0.27 eV). Using the modified BEP relation for sites on the (110)-B surface resulted in a worse fit than the original BEP relation. Therefore, we use the original BEP relation for this surface, which has a maximum deviation of 0.27 eV. The reason that the modified BEP relation is not suitable for sites on the (110)-B surface is that H+ /e – transfer and O – O bond formation occur concertedly on these sites, rather than sequentially as they do for sites on the other two surfaces. The large deviations in the BEP-predicted activation barriers from the calculated barriers on the (110)-A and (110)-B surfaces (and the smaller deviations on the (311) surface as well) are likely to originate from two effects 34,35 – (1) the mixture of other electronic configurations into the electronic structure of the transition state other than those corresponding to 16

ACS Paragon Plus Environment

Page 17 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

the reactant, intermediate, and product configurations, and (2) variations in the resonance integral between the reactant, intermediate, and product configurations. Without performing additional electronic structure analysis, it is difficult to specify the exact origin of these deviations. It is worth mentioning, though, that if the origin of these deviations could be identified, it would provide an additional lever for tuning the active site performance. The modified BEP plot in the right panel of Figure 2 also shows that the use of different BEP parameters for the (311) and (110)-A surfaces is warranted, as these parameters are very different for the two surfaces. The BEP slope α is 0.58 ± 0.03 for the (311) sites but 0.34 ± 0.04 for the (110)-A sites. Also of consequence, the BEP intercept β is lower for the (311) site (0.97 ± 0.04 eV) than for the (110)-A site (1.41 ± 0.04 eV), which amounts to a difference in the reaction rate constant at 298 K of more than seven orders of magnitude. Thus, we can say that the BEP relations exhibit a high degree of structural sensitivity, with a more favorable relation existing for the (311) site compared to the (110)-A site. As with the deviations discussed in the previous paragraph, however, it is difficult to specify the origin of the structural sensitivity in the absence of extensive electronic structure analysis; however, it is likely related at least in part due to the coordination environment of the reactive oxo Oa .

Equivalence of Electrochemical and non-Electrochemical Water Addition Pathways The fact that the transition state for water addition on the (311) and (110)-A surfaces occurs after the H+ /e – transfer is complete (c.f. Scheme 2) has a very important consequence. Since all states preceding the transition state should be in quasi-equilibrium, this means that the intermediate state with the hydroxyl radical in Scheme 2 is quasi-equilibrated with the resting state of the active site. This entails that only the free energy of the intermediate state is important – not the path by which the site arrives at this intermediate state – as long as the highest energy transition state is associated purely with O – O bond formation. 17

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 49

Therefore, it would be equally valid to consider that the intermediate state is obtained by the electrochemical process shown in the bottom pathway of Scheme 2 in which a water molecule dissociates into a hydroxyl radical at the active site, with a proton going to the electrolyte and an electron going to the electrode. This would be consistent with the electrochemical water addition step in reaction 30 ,

(Mn+1 −O)a + OH− −−→ (Mn −OOH)a + e−

.

{30 }

Essentially, this electrochemical water addition pathway has the same transition state (O – O bond formation) as the non-electrochemical pathway, the only difference being in the trajectory prior to the transition state. As a result, the turnover frequency in this case is independent of which water addition pathway we choose for the catalytic cycle. For the purposes of catalyst design, it turns out to be more intuitive to consider the electrochemical pathway for water addition rather than the non-electrochemical pathway. Another advantage of using the electrochemical pathway is that it is the pathway used most commonly in thermodynamics-based OER catalyst screening studies (c.f. refs 5,12–14), which facilitates comparison between the two approaches. Since the electrochemical water addition pathway is thermodynamically and kinetically equivalent to electrochemical oxidation of the (M – O)b center (step 2) followed by nonelectrochemical water addition (step 3), the activation barrier ∆G‡30 and reaction free energy ∆G30 are given by ∆G‡30 = ∆G‡3 + eη2

,

(8)

∆G30 = ∆G3 + eη2

.

(9)

It turns out that these are the same activation barrier and reaction free energy used in the BEP relationship in eq (6) from the previous section. Thus, the modified BEP relation we 18

ACS Paragon Plus Environment

Page 19 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

used in eq (6) is actually the BEP relation for water addition by the electrochemical pathway, explaining why we chose that particular definition as well as the label 30 . Also, since this is an electrochemical step, we assign it a redox potential η3 , where eη3 = ∆G30 . In addition to the primary electrochemical water addition pathway we have discussed so far that proceeds from an active site where the (M – O)b center is in the reduced state, a secondary electrochemical pathway also exists that instead proceeds from an active site where the (M – O)b center is in the oxidized state. On the undoped (311) surface, this secondary pathway has an activation barrier at the OER equilibrium potential of 1.82 eV, significantly higher than the activation barrier of 1.30 eV for the primary pathway. The reason for the higher activation barrier of the secondary pathway could be due to two factors. First, when the (M – O)b center is reduced, the proton bound to the oxygen can form a strong hydrogen bond with the hydroxyl oxygen in the transition state, as can be seen in Figure 1. The second reason is that the increased oxidation state on the metal atom results in a stronger stabilization of the oxo in the (M – O)a center, making it less reactive for O – O bond formation. Since the activation barrier of the secondary electrochemical pathway is so much higher than for the primary pathway, it should not be kinetically relevant except when the electrode potential is on the order of 0.5 V greater than the redox potential η2 . As will be seen later, a site for which η2 is significantly less than the electrode potential will be a very inefficient catalyst for the OER. As such, we can neglect this secondary electrochemical pathway for any reasonably efficient OER catalyst.

Simplified Kinetic Model of the OER Even though changing the definition of the water addition step from a non-electrochemical pathway to an electrochemical pathway effectively does not change the turnover frequency of the catalytic cycle in eq (1), it is more intuitive to write the rate expression in terms of

19

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 49

the rate constant k30 of the electrochemical water addition pathway, 

 eη   k30 exp   kB T TOF = 







 .

(10)

  e(η − η1 )  e(η − η2 ) 1 + exp   1 + exp − kB T kB T

Here, k30 is defined as the electrochemical water addition rate constant at the equilibrium OER electrode potential and is related to the rate constant of the non-electrochemical water addition pathway by 



 eη2  k30 = k3 exp −  , kB T

(11)

Both rate constants are related to the activation barriers via the Eyring equation,   kB T ∆G‡ k= exp − . h kB T

(12)

The exponential factor of η in the numerator of the above rate expression accounts for the dependence of the rate constant on the electrode potential. Interestingly, this implies that the symmetry factor for electrochemical water addition should be close to unity, rather than 0.5, when the transition state is associated purely with O – O bond formation. This is because the electrochemical part of the pathway (H+ /e – transfer) is already complete prior to the transition state for O – O bond formation. Thus, it should not be possible to tell from an experimentally measured Tafel slope whether water addition is occurring electrochemically or non-electrochemically – nor does it matter, since the turnover frequency is independent of how the hydroxyl radical in the intermediate is formed because this occurs prior to the transition state common to both pathways. The TOF given by eq (10) is plotted as a function of the electrode potential in Figure 3 for three sites on the (311) surface – a site where both metal centers are doped with Ni 20

ACS Paragon Plus Environment

Page 21 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(Ni/Ni), a site where Mb is doped with Cr (Co/Cr), and a site where Ma is doped with Fe (Fe/Co). The TOF on each site exhibits three regimes, each having a different dependence on the electrode potential. At low electrode potentials, the logarithm of the TOF increases linearly with respect to the potential with a slope corresponding to a transfer coefficient of two (or a Tafel slope of 30 mV/dec). When η > η1 for the particular site, the slope of the TOF decreases by half (corresponding to a transfer coefficient of one or a Tafel slope of 60 mV/dec) and then becomes zero for η > η2 . These three regimes exist due to changes in the resting state of the active site according to the redox processes 1 and 2. The transition from the first regime to the second at η = η1 is due to the oxidation of the (M – O)a center according to 1. Prior to this transition, the (M – O)a center exists in an inactive reduced state with a thermally excited population in the active oxidized state that increases exponentially with the electrode potential. This is responsible for the increased slope of log(TOF) with respect to η in this regime. In the second regime, the (M – O)a center is in the oxidized state so that the TOF is equal to the intrinsic rate of the electrochemical water addition pathway. In the third regime, the (M – O)b center becomes oxidized to a state which is inactive for the primary electrochemical water addition pathway. Water addition must proceed from thermally excited active sites having a reduced (M – O)b center, whose population decreases exponentially with the electrode potential. This cancels out the increase in the rate of electrochemical water addition with increasing electrode potential, so the TOF is constant in this regime. The differences between the three active sites in the plot in Figure 3 are due to differences in the two redox potentials of these sites. The Ni/Ni site has the highest η1 , followed by Co/Cr and then Co/Fe. Therefore, the transition to the second regime occurs at an increasingly higher electrode potential going from the last site to the first. The Ni/Ni site also has higher η2 than the other two sites. This causes the Ni/Ni site to transition to the third regime at higher electrode potentials than the other two sites. Interestingly, there appears to be a linear correlation between the points on each curve

21

ACS Paragon Plus Environment

The Journal of Physical Chemistry

3 0

-1

)

-3

lo g ( T O F /s

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 49

-6 N i/N i C o /C r F e /C o

-9 -1 2 η = η1

-1 5 0 .0

0 .2

η = η2

,C o /C r

0 .4

0 .6

,C o /C r

0 .8

1 .0

η( V ) Figure 3: TOF versus electrode potential for three active sites on the (311) surface. The two vertical lines indicate the regime transitions for the Co/Cr site at η = η1 and η = η2 . The points indicate the TOF at η = η1 for each of the three sites.

22

ACS Paragon Plus Environment

Page 23 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

corresponding to the transition between the first and second regimes at η = η1 . In fact, this correlation exists for all of the sites on the (311) surface. Its origin is a linear relationship between the activation barrier for the electrochemical water addition pathway and the redox potential η1 . Since the activation barrier is linearly related to the reaction free energy by the BEP relation, this correlation must be due to a second linear relationship between the reaction free energies of electrochemical water addition and oxidation of the (M – O)a center. We examine the nature and origin of this relationship in the next section.

Energetic Relationships between the Water Addition Reaction and Pre-equilibrium Oxidation of the Active Site The correlation observed in the last section between the free energies of electrochemical water addition (step 30 ) and oxidation of the active site (step 1) is due to similarities in the chemical bonds that are broken and formed during these two steps. To see this, we first note that the oxidation in step 1 is equivalent to dissociating an O – H bond followed by electrochemical oxidation of the resulting H atom (H + OH – −−→ H2 O + e – ). Consequently, the free energy of this step, eη1 , can be written as

eη1 = eηH + ∆GO−H

,

(13)

where ∆GO−H is the O – H bond dissociation free energy and ηH is the redox potential for hydrogen atom oxidation (calculated to be −3.65 V with respect to the equilibrium OER potential). We then note that electrochemical water addition is equivalent to the electrochemical oxidation of water (in the form of hydroxide under alkaline conditions) to form a hydroxyl radical (OH – −−→ OH + e – , with redox potential ηH2 O ), followed by binding of the hydroxyl to Oa to form the hydroperoxo product. As such, the free energy of this step can be written as eη3 = eηH2 O − ∆GO−OH 23

,

ACS Paragon Plus Environment

(14)

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 49

where ∆GO−OH is the free energy for dissociation of the O – OH bond in the hydroperoxo. Thus, aside from the site-independent quantities ηH and ηH2 O , the two reaction free energies are determined by the O – H and O – OH bond dissociation free energies. In our previous study of different active sites on undoped Co3 O4 , 16 we found that ∆GO−H and GO−OH are not independent but are related by a linear scaling relationship, due to the similarity between O – OH and O – H bond formation at the quantum chemical level. This scaling relation has the form

∆GO−OH = a∆GO−H + b ,

(15)

and is also found to apply to the more extensive set of doped sites shown in Figure 4. We find that the sites on all surfaces can be described by the same scaling relation, with a slope of 0.99 ± 0.03. This is nearly identical to the slope of 1.00–1.05 obtained in other work 12,13 and the near-unity value is consistent with the similarity between the O – H and O – OH bonding processes at the quantum chemical level. Since the slope is equivalent to unity within the uncertainty of the regression, we use a simpler model in which a is set to unity for all three surfaces, which is drawn in Figure 4. We would like to point out that the mathematical form of the scaling relation we use to correlate O – H and O – OH bond dissociation energies is equivalent to the mathematical form used by Man et al. 13 to correlate the formation free energies of OH and OOH intermediates on the surface, ∆GHOO∗ = ∆GHO∗ + b0

.

(16)

As derived in the Supporting Information, b and b0 are related by

b0 = −b + eηH + eηH2 O + 2.46 eV .

(17)

Using this, the value of b0 corresponding to the scaling relation we obtain in Figure 4 is

24

ACS Paragon Plus Environment

Page 25 of 49

(3 1 1 ) (1 1 0 )-A (1 1 0 )-B

2 .5

∆G ( O - O H ) ( e V )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

2 .0

1 .5

1 .0 a =1 b = −2.42 ± 0.02 e V

3 .0

3 .5

4 .0

4 .5

5 .0

∆G ( O - H ) ( e V ) Figure 4: (left) Scaling relation between the O – H and O – OH dissociation energies.

25

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 49

b0 = 3.46 eV, which is higher than the value of b0 = 3.0–3.2 eV calculated by Rossmeisl et al. 13 and Koper. 36 This is possibly due to the fact that they include zero-point vibrational energy in their calculated free energies and use a different reference for water (in the gas phase at 0.035 bar; we use hexagonal ice for reasons discussed in the Supporting Information). The scaling relation between O – OH and O – H bond dissociation energies leads to an interesting way to think about the mechanics of the water addition step. To see this, we first use the scaling relation between the O – OH and O – H bond dissociation energies in eq (15) to write a relation between the free energies of the water addition and pre-equilibrium oxidation steps, making use of eqs (13), (14),

eη1 + eη3 = eηH + eηH2 O − b = eη1+3

.

(18)

Since ηH , ηH2 O , and b are the same for all sites, the sum of η1 and η3 is a site-independent constant equal to 1.00 V, which we will call η1+3 . Therefore, any increase in the free energy of the pre-equilibrium oxidation step results in an equal decrease in the free energy of the water addition step. The physical interpretation of this is that the active site stores the energy supplied by the electrode potential during the oxidation step and releases this energy to drive the water addition step. The origin of the scaling relation that we identified in our previous work explains how the active site stores this energy. We found that the O – H and O – OH bond dissociation energies are related by the energy to transfer a hole from the redox-active metal center Ma to the oxygen atom Oa that forms the O – O bond with the hydroxyl. 16 Thus, the active site stores the energy supplied by the electrode potential in the form of a high energy hole on this metal center, as depicted in Scheme 3. This hole then removes the electron from the antibonding orbital of the nascent O – O bond during the water addition step. The reason this storage step is necessary is because the electron is unlikely to transfer directly from the O – O antibonding orbital to the electrode, as this would involve electron transfer over a long

26

ACS Paragon Plus Environment

Page 27 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

oxidation

water addition

H2O + e−

e−



1 OH−

𝜂

OH−

𝜂1



1

Scheme 3: During the pre-equilibrium oxidation step (1), the energy supplied by the external potential η is stored in a high energy hole (light blue) on the Ma center. During the water addition step (30 ), this energy is released to remove an electron from the antibonding orbital (light green) of the nascent O – O bond. distance. The electron must first transfer to a state having sufficient overlap with the O – O antibonding orbital, most likely the hole on the metal center. We will see in the next section how this energy storage process controls the overall activity of the site and how this leads to a kinetics-based design principle quite different from the widely used thermodynamics-based design principle. A similar process occurs during non-electrochemical water addition on the (110)-B site, except that the active site stores energy in both oxidation steps 1 and 2, rather than only in step 1. This can be seen by using the scaling relation (eq (15)) to write the free energy of non-electrochemical water addition in terms of both redox potentials η1 and η2 (making use of eqs (9), (13), and (14)), ∆G3 = e(η1+3 − η1 − η2 ) ,

(19)

(η1+3 is the same as in eq (18)). The energy stored during step 2 in the form of a high energy hole on Mb is used to drive the H+ /e – transfer. As with the electrochemical pathway, the energy stored during step 1 in the form of a high energy hole on Ma is used to drive O – O formation by removing an electron from the antibonding orbital of the nascent O – O bond.

27

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 49

A Kinetics-based Catalyst Design Strategy for the OER We now arrive at our main objective, which is to develop a computationally cheap method of identifying an active site for the OER having the maximum TOF at a given electrode potential. Of course, criteria other than the TOF also must be considered, such as the stability of the active site under OER conditions or its abundance on the surface of catalyst nanoparticles, but these are beyond the scope of the present study. In a rigorous approach, identification of the optimal active site would require costly transition state optimizations in order to calculate the rate constants for the water addition step, which is assumed to be rate limiting (an even more rigorous approach would not make this assumption). This can be avoided by using the scaling and BEP relations developed in the previous sections to estimate the activation barrier for water addition in terms of the redox potential(s) of the active site. For sites on the (311) and (110)-A surfaces, the electrochemical and non-electrochemical water addition pathways are kinetically equivalent, so we choose the former since the BEP relations are based on it. The activation barrier of this pathway can be estimated from the redox potential η1 by combining the BEP relation in eq (6) with the relationship between η1 and η3 in eq (18), ∆G‡30 = αe(η1+3 − η1 ) + β

.

(20)

Substituting the numerical values of η1+3 and the BEP parameters gives ∆G‡30 = 0.58(1.00 eV − eη1 ) + 0.97 eV

(21)

∆G‡30 = 0.34(1.00 eV − eη1 ) + 1.41 eV

(22)

for the (311) site and

for the (110)-A site. Using the resulting rate constant (obtained via the Eyring equation) in eq (10) gives the TOF for an active site with redox potentials η1 and η2 operating at a

28

ACS Paragon Plus Environment

Page 29 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

given electrode potential η. For sites on the (110)-B surface, water addition occurs by the non-electrochemical pathway. The activation barrier for this pathway can be obtained in terms of the redox potentials by combining the BEP relation for this pathway in eq (2) with the relation between ∆G3 , η1 , and η2 in eq (19), leading to ∆G‡3 = αe(η1+3 − η1 − η2 ) + β

.

(23)

Substituting for η1+3 and the BEP parameters gives ∆G‡3 = 0.62(1.00 eV − eη1 − eη2 ) + 1.20 eV

(24)

The TOF is then given by eq (1), using the rate constant k3 obtained from ∆G‡3 by the Eyring equation. Note that different values of α and β are used for each of the three surfaces. Using either eq (20) or (23), a plot can be made showing the TOFs at a fixed electrode potential of sites having different values of the redox potentials η1 and η2 . One can then identify the redox potentials leading to the highest TOF, defining the optimal site at that electrode potential. While this activity plot will change depending on the electrode potential, these changes will only involve shifts along the η1 and η2 axes and a scaling of the TOF axis. This can be seen by first normalizing the TOF by a quantity TOF∗ , whose physical relevance will become apparent shortly. On the (311) and (110)-A sites, this quantity is defined as 



 eη  TOF∗ = k3∗0 exp   , kB T

(25)

where k3∗0 is the rate constant for electrochemical water addition for a site with η1 = η, obtained from eq (20) via the Eyring equation. Normalizing the TOF in eq (10) by TOF∗

29

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 49

gives, k30 /k3∗0    .  = TOF∗  e(η2 − η)   e(η1 − η) 1 + exp 1 + exp  −    kB T kB T TOF

(26)

Using eq (20) along with the Eyring equation for k30 and k3∗0 , we can see that the numerator depends only on η1 −η. Since the denominator only depends on η1 −η and η2 −η, a generalized activity plot valid for any value of η can be constructed by plotting this normalized TOF with respect to η1 − η and η2 − η. The activity plot for the (110)-B surface can be generalized to any value of the electrode potential η in a similar way. For this surface, TOF∗ is defined as

TOF∗ = k3∗

,

(27)

where k3∗ is the rate constant for non-electrochemical water addition for a site with η1 = η2 = η, obtained from eq (23) via the Eyring equation. Normalizing the TOF in eq (1) by this quantity also leads to a result that is dependent only on η1 − η and η2 − η, TOF TOF∗

k3 /k3∗    .  e(η1 − η)   e(η2 − η) 1 + exp 1 + exp      kB T kB T

=



(28)

The generalized activity plots defined by eqs (28) and (26) are shown in Figure 5 for the (311) and (110)-B surfaces. The plot for the (110)-A surface is qualitatively very similar to the (311) surface, so it is shown only in the Supporting Information. For all three surfaces, one can see that the TOF approaches a maximum equal to the quantity TOF∗ defined in eq (25) ((311) and (110)-A surfaces) or (27) ((110)-B surface). Thus, TOF∗ is a measure of the TOF of the optimal site on a given surface, essentially the maximum attainable activity, at a fixed electrode potential. 30

ACS Paragon Plus Environment

Page 31 of 49

0.3

0.3

0.2

0.2

0 -1 -2

0.1 η2 − η (V)

0.1 η2 − η (V)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

0.0

-3 -4

0.0

-5

-0.1

-0.1

-0.2

-0.2

-0.3 -0.3

-0.3 -0.3

-6 -7 -8

-0.2

-0.1

0.0

0.1

0.2

0.3

-0.2

-0.1

η1 − η (V)

0.0

0.1

η1 − η (V)

(311) µ³-oxo

0.2

0.3 -9 log(TOF/TOF*)

(110)-B η-oxo

Figure 5: Activity plots showing the logarithm of the TOF for sites with different values of the redox potentials η1 and η2 . For applicability to any electrode potential η, the TOF is scaled according to eqs (26) and (28). The plot for the (110)-A surface is similar to the one for the (311) surface, so it is included in the Supporting Information. The plots for both sites have the usual volcano shape along the η1 axis, with a maximum occurring at approximately η1 = η. The behavior is the same as in our previous work and is related to the transition of the active site from an inactive reduced state to an active oxidized state. When η1 < η and η2 > η, the active site exists primarily in the state from which electrochemical water addition proceeds directly, where the (M – O)a center is oxidized and the (M – O)b center is reduced. Thus, the denominator in eq (10) is unity and the turnover frequency is equal to the intrinsic rate constant for electrochemical water addition. Since this rate constant increases exponentially as η1 increases due to the BEP and scaling relations, a site with a higher value of η1 will have a higher turnover frequency. In contrast, when η1 is above the electrode potential, the active site primarily exists in a state where the (M – O)a center is reduced, with only a small thermally excited population in the oxidized state that increases exponentially as η1 decreases. This increase in the active state population is stronger than the decrease in the rate constant so that the turnover frequency increases as η1 decreases. 31

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Along the η2 axis on the activity plots for the (311) and (110)-A sites (the latter is in the Supporting Information), however, the behavior is very different from a volcano plot. For η2 < η, the TOF increases exponentially with η2 but then becomes constant once η2 > η. As such, the two-dimensional activity plot resembles more the end of a ridge than a volcano. Since the BEP relation for the water addition rate constant is independent of η2 , this behavior is caused by the thermodynamics of the second oxidation step 2. When η2 > η, the (M – O)b center is reduced. This is the state from which the electrochemical water addition pathway proceeds, so that the turnover frequency is independent of η2 . When η2 < η, the (M – O)b center is oxidized so that the turnover frequency decreases with η2 in proportion to the thermal population of sites in which the (M – O)b center is reduced. In contrast, the activity plot for sites on the (110)-B surface has the typical “volcano” shape rather than the “ridge” shape seen for sites on the other two surfaces. This is due to the fact that the water addition reaction is governed by a BEP relationship based on the reaction free energy of the non-electrochemical water addition pathway, which depends on both redox potentials. As discussed already, this arises from concerted proton/electron transfer and O – O bond formation in the transition state on these sites. As a result of the different BEP relation, we get the same dependence in the activity plot on η2 − η as on η1 − η (the plot is symmetric across the line η1 = η2 ). From the above discussion, we can see that the maximum TOF occurs at approximately η1 = η and η2 > η for the (311) and (110)-A surfaces, corresponding to the “ridge” on the activity plot. For the (110)-B surface, it occurs approximately at η1 = η2 = η, corresponding to the peak of the volcano. These conditions define the design criteria for the optimal active site on a given surface at a fixed electrode potential. In practice, the optimal site will have redox potentials within kB T of those specified by the above criteria due to the fact that the activity plot is not symmetric across the line η1 = η (and also η2 = η on the (110)-B plot) unless the BEP slope α = 0.5. Also, the maximum TOF is always slightly less than TOF∗ (within a factor of four) due to the fact that the denominators in the rate expressions are

32

ACS Paragon Plus Environment

Page 32 of 49

Page 33 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

greater than unity for the optimal site. However, these deviations are much smaller than the uncertainty of the BEP and scaling relations so that it is reasonable to neglect them in order to keep the design criteria simpler. The design criterion for the (311) and (110)-A sites has a very important difference from the design criterion in our previous work. There, we found that the activity had the usual volcano shape along both the η1 and η2 axes. 16 This occurs for the same reason that we observe this shape for the (110)-B activity plot in the present work – that we used a BEP relation based on the reaction free energy of the non-electrochemical water addition pathway. As such, the design principle we formulated previously stated that the optimal active site should have both redox potentials equal to the electrode potential. However, the more extensive set of sites examined in the current work reveals that it is more appropriate to base the BEP relation on the electrochemical water addition pathway, which has no dependence on η2 and results in an activity plot that reaches a plateau along the η2 axis. This result has an extremely fundamental impact on the design criteria for the OER active site. Since η2 can have any value above η, we have significantly more freedom in selecting an active site having the maximum activity. The physical origin of the design criterion η1 = η is due to the mechanism discussed in the last section by which the active site stores the energy supplied by the external potential during the pre-equilibrium oxidation step 1 and releases it to drive the water addition step. Because of this, it is helpful to think of the active site as analogous to a spring with η1 being equal to the amount of energy required to fully compress the spring. This is depicted graphically in Scheme 3. When η1 > η, the energy supplied by the external potential is not sufficient to oxidize the active site (compress the spring) so that only a low population of oxidized sites exists. When η1 < η, the active site cannot store all of the energy supplied by the external potential (i.e. the spring is not large enough) so that the rest is lost during the oxidation step. For the (110)-B surface, a similar argument applies to the additional design criterion η2 = η since step 2 also stores energy in the active site to drive the non-

33

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 49

electrochemical water addition reaction. Finally, we would like to compare the differences in the maximum attainable activity (TOF∗ ) for the three surfaces that arise from the structure sensitivity of the BEP and scaling relations. These are plotted for the three surfaces in Figure 6 with respect to the electrode potential. We can see that for any electrode potential, the (311) surface has a significantly higher maximum attainable activity than the other two surfaces. This is due to the more favorable BEP relation for this surface (Figure 2), for which the slope is greater and the intercept is lower than on the other two surfaces. Not only is the maximum attainable activity lowest for the (110)-B surface, but also the design criterion of η2 = η on this surface is more restrictive than η2 > η on the other two surfaces. For both of these reasons, it appears that sites forming the O – O bond to an η-oxo are not good candidates for the OER. Interestingly, most studies of the OER focus on such a site (c.f.

refs 12–14). Instead,

surfaces with sites that form the O – O bond to a μ3 -oxo, such as the (311) surface, appear to be the best candidates to focus the search for an improved OER catalyst. Figure 6 also shows how the calculated TOFs of the actual active sites compare to the TOFs estimated from the BEP and scaling relations. For each site, the TOF calculated at η = η1 is plotted – in other words, at the electrode potential where that particular site has the optimal value of η1 . For sites on the (311) and (110)-A surfaces, deviations of the computed TOFs from the TOFs predicted by the correlations are due to deviations in the BEP and scaling relations. Of particular interest, the undoped Co/Co site on the (311) surface has a TOF about two orders of magnitude higher than the TOF predicted by the correlations. It would be interesting to understand the quantum chemical origin of this deviation to identify if there is another property of the active site in addition to η1 that can be tuned to maximize the activity. For sites on the (110)-B surface, the deviations of the computed TOFs from the TOFs predicted by the correlation are due not only to deviations in the BEP and scaling relations, but also to the more stringent design criteria for sites on this surface (η1 = η2 = η). None of

34

ACS Paragon Plus Environment

Page 35 of 49

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 C o /C o /t 21 C o /C o /b 22 23 C o /C o 24 M n /C o 25 N i/N i 26 C r/C r F e /C o /b 27 F e /C r/b 28 C o /C o 29 F e /F e /b 30 C o /C r C r/F e 31 32 33 C r/C r M n /M n /t C r/C o 34 C r/C o M n / F e / t 35 C r/C r F e /C o F e /C r 36 F e /F e 37 C r/C r M n /M n /b 38 39 40 41 42 V /C r/b 43 44 45 V /C r/t 46 V /M n /t 47 V /M n V /F e /t V V / V / 6: 48V / V / t / M n Maximum attainable TOF for the optimal site on each surface as a function of Figure V /V 49 electrode potential. Points corresponding to the calculated sites at η = η1 are shown (labeled 50 Ma /Mb ), as well as for a site with η1 = 0.50 V at all values of η (black line). 51 52 53 54 55 56 57 58 35 59 ACS Paragon Plus Environment 60

(3 1 1 ) (1 1 0 )-A (1 1 0 )-B

1 0

lo g ( T O F /s

-1

)

5 0

-5

-1 0 -1 5

0 .0

0 .2

0 .4

0 .6

η( V )

0 .8

1 .0

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the actual sites examined on this surface have equal redox potentials η1 and η2 so that even when the first condition η1 = η holds, the second condition η2 = η does not hold and the site is sub-optimal.

Comparison with the Thermodynamics-based Screening Approach The kinetics-based design criteria we presented above differ from the thermodynamics-based design criterion 11–13 used in most computational screening studies of the OER to date (c.f. refs 5,14,36). This latter design criterion is based on the catalytic cycle occurring via the electrochemical water addition and O2 elimination pathways depicted in Scheme 1 (although a different numbering of the steps is often used). The criterion states that on the optimal catalyst, the redox potentials of all four electrochemical steps (1, 30 , 4, and 60 ) should be equal to the equilibrium OER potential. Due to the scaling relations between the O – H and O – OH bond dissociation energies, however, this limit cannot be obtained for a real material. Rather, the optimal metal oxide catalyst is predicted to have the redox potentials of steps 1 and 30 , the electrochemical oxidation and water addition steps, equal to one another (η1 = η3 ). 13 It was found that this occurs when both of these redox potentials are equal to 0.37 V. A somewhat higher value of η1+3 /2 = 0.50 V results from eq (18) based on the scaling relation we obtained in Figure 4. In Figure 6, we have also plotted the TOF of such a catalyst (η1 = 0.50 V and η2  η) as a function of the electrode potential using the rate expression in eq (10) and the BEP and scaling relations for the (311) surface. As one might expect, this site has the same TOF as the optimal (311) site at an electrode potential of η = 0.50 V, however, it has a lower TOF at all other electrode potentials. This is related to the fact that the kinetics-based approach predicts a different optimal site for each electrode potential – in contrast, the thermodynamics-based approach predicts a single optimal site independent of the electrode potential. Perhaps the reason that the thermodynamic-based approach often works in practice is that the optimal redox potential for the second step in this approach is close to the electrode potentials at which current state-of-the-art OER catalysts typically 36

ACS Paragon Plus Environment

Page 36 of 49

Page 37 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

operate (c.f. ref 2). The thermodynamics-based design criteria do not identify the optimal catalyst when η 6= 0.50 V because the underlying approach does not consider that (1) water addition is the rate limiting step under most conditions and that (2) the water addition activation barrier decreases as the redox potential η1 increases. This is shown graphically in Scheme 4, which depicts a series of potential energy surface (PES) diagrams for both the thermodynamicsbased model and the kinetics-based model. The PES diagram for the thermodynamics-based model assumes a small activation barrier of 0.2 eV for each of the two steps (similar to the barrier assumed in ref 37), while the kinetics-based model assumes this barrier only for the oxidation step 1, using the BEP barrier in eq (6) for the water addition step. η1 = 0.60

1

0.60 0.50 0.40

0.50



0.40

(M3+–OH)a

0.60

(M3+–OOH)a

(M4+–O)a

η1 =

0.50 0.40

0.30

η = 0.40 V

(M3+–OH)a

0.70 0.60

0.50 0.40

0.40



0.60

1

0.50

(M3+–OOH)a

0.50

(M4+–O)a

η = 0.50 V

η = 0.60 V

Scheme 4: Potential energy surface diagrams at different electrode potentials η and redox potentials η1 for the kinetically relevant part of the OER according to the (top) thermodynamics-based and (bottom) kinetics-based models. The two steps shown are electrochemical oxidation of (M – O)a (step 1) and electrochemical water addition (step 30 ). We first examine the PES diagrams for the thermodynamics-based approach at three different electrode potentials to see how changes in the redox potential η1 affect the apparent activation barrier. The solid black lines show the PES for a the optimal active site predicted by the thermodynamics-based approach, with η1 = η3 = 0.50 V. The dashed blue lines show how changing the redox potential η1 of the active site (along with a corresponding 37

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

equal and opposite change in η3 due to the scaling relation) affects the PES. This has the effect of changing the free energy of the oxo intermediate while leaving the free energies of the hydroxo and hydroperoxo states unchanged. The most important thing to note is that the lowest apparent activation barrier occurs when the free energy of the oxo intermediate falls between the free energies of the hydroxo and hydroperoxo states. Within this window, defined by η1 = 0.50 V ± |η − 0.50 V|, the apparent activation barrier is insensitive to η1 . Since the site with η1 = η3 = 0.50 V falls within this window for all electrode potentials, this site will always be among the sites giving the maximum TOF regardless of the electrode potential. The PES diagrams in which the water addition barrier is described by eq (6) are significantly different. First, water addition is always the rate limiting step regardless of the electrode potential or the redox potential of the active site (within reasonable bounds). In contrast, the thermodynamics-based approach predicts that either of the two steps could be rate limiting, depending on the values of η and η1 . Second, the intrinsic activation barrier for water addition is now a function of η1 . Together, these have the effect that the site with eta1 = 0.50 V no longer gives the minimum apparent barrier at all electrode potentials. Instead, the site with η1 = η is the optimal site (shown with solid black lines in Scheme 4). On sites with higher redox potentials, the site must first climb along the PES up to the oxo intermediate before water addition can proceed. On sites with lower redox potentials, the site must climb out of the low energy oxo intermediate during water addition. Both cases (shown with dashed blue lines in Scheme 4) lead to a higher apparent barrier. To summarize, the kinetics-based approach predicts that the optimal active site will change depending on the electrode potential while the thermodynamics-based approach predicts a single optimal active site independent of the electrode potential. The former approach predicts that there is a trade-off between activity and efficiency – by going to an active site with a higher redox potential, one can obtain higher activity but at the price of operating at a higher electrode potential. The implication is that the catalyst should be matched to the

38

ACS Paragon Plus Environment

Page 38 of 49

Page 39 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

desired current density of the electrolyzer. In the case of a photoelectrocatalytic cell, where the electrode potential is related to the band gap of the semiconductor, it would be expected that one property of an optimal catalyst would be a redox potential close to this band gap. It is worth mentioning that the conclusions obtained with the kinetics-based approach are dependent on the water addition step being the rate limiting step. The same conclusions would not be obtained from a kinetic model in which all steps have activation barriers that follow identical BEP relations, as is done by Hansen et al. 37 Performing a similar analysis as in Scheme 4 for their model predicts that the optimal site has η1 = 0.5 V when η ≥ 0.5 V, and η1 = η when η ≤ 0.5 V. This is due to the fact that with their model the first step becomes rate limiting, rather than water addition, when η > 0.5 V and η1 > 0.5 V. Finally, we would like to note that in order to estimate the activity of a given active site, both the kinetics-based and thermodynamics-based approaches require the computation of only a single descriptor, the redox potential η1 (technically, that latter approach uses the difference in M – O and M – OH binding energies, 12,13 which differs from η1 by a constant). Therefore, the two approaches have the same level of computational difficulty. Because of this and the fact that the kinetics-based approach is always able to identify a site that has equal or greater activity than the one identified by the thermodynamics-based approach, the former approach is preferable.

Limits of the Kinetics-based Screening Approach While the kinetics-based active site design strategy performs better than the thermodynamicsbased design strategy, it still has limitations due to the assumptions it is based on. The first of these assumptions is that the water addition step is rate limiting. In our previous work, however, we found that the O2 elimination step can also be rate limiting on the undoped (311) surface at lower electrode potentials. 15 A future study should thus examine under which conditions this step becomes rate limiting and what properties of the active site control the rate of this step. 39

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

A related assumption is that we consider only four possible states of the active site when deriving the expressions for the TOF in eqs (1) and (10). This ignores the fact that the active site is embedded in an extended surface that can undergo additional redox transitions. As nearby sites become oxidized at increasing electrode potentials, the energetics of the OER pathway on the active site can be modified. This was seen in our previous work, 15 where it led to a more complex dependence of the TOF on the electrode potential than that shown in Figure 3. Regardless, we do not expect these effects to alter the proposed design criteria. This is because the design criteria are based only on the form of the rate expression in the vicinity of η = η1 so that modifications to the rate expression at electrode potentials far from this value are not relevant to the design criteria. Additionally, we do not expect the BEP and scaling relations to be altered by changes in the state of the surface in which the active site is embedded. Nonetheless, the redox potentials of the active site (η1 and η2 ) on which the design criteria are based can be modified by redox transitions of nearby sites on the surface. Therefore, it is important when calculating the redox potentials of the active site that these neighboring sites exist in the correct state for the electrode potential of interest. Under certain conditions, it could also be possible that one of the electrochemical proton/electron transfer steps (e.g. , step 1) becomes rate limiting. However, we think this is unlikely for the optimal sites where η1 = η. This is because step 1 has a free energy of zero for such a site, so that the activation barrier for this step should be on the order of 0.1–0.2 eV in the absence of a significant electric field across the double layer. This barrier would only be higher if the catalyst surface is negatively charged at the operating potential – i.e. for a catalyst with a very high potential of zero charge at the electrolyte pH. As currently only limited methodology exists to compute barriers of electrochemical reactions at a charged electrochemical interface of a semiconducting transition metal oxide, it is unfortunately not possible to examine this assumption at present. The final assumption of our kinetic model is that the O – O bond formation step occurs non-electrochemically – i.e. the electron from the O – O antibonding orbital transfers locally

40

ACS Paragon Plus Environment

Page 40 of 49

Page 41 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

to a hole on the Ma center of the active site. While we have mentioned already that it is unlikely for this electron to transfer directly to the electrode, as this would involve longrange electron transfer, it may be possible for this electron to transfer to a hole in the valence band of the catalyst surface. This would certainly be the case on conducting oxides such as IrO2 and RuO2 , or on an oxide that has conducting surface states. It could also occur on a polaronic p-type semiconducting oxide that develops an accumulation layer under electrochemical conditions (Co3 O4 is likely such a material based on the position of the valence band maximum 38 ). In these cases, however, the free energy and activation barrier of the O – O bond formation process would still be independent of the electrode potential. This is because the O – O antibonding orbital is inside the potential gradient of the electrochemical double-layer, so that it has a constant energy relative to the Fermi level of the catalyst that is independent of the electrode potential. Thus one could still model it as a non-electrochemical process with a barrier that is independent of the electrode potential, since all charge transfer occurs inside the electrochemical double-layer.

Conclusions By examining a large number of 3d transition metal-doped active sites on three surface terminations of Co3 O4 ((311), (110)-A, and (110)-B), we were able to formulate kinetics-based design principles for identifying the optimal active site to carry out the oxygen evolution reaction at a given electrode potential. The design principles are applicable for any active site that carries out the OER by a catalytic cycle in which a rate limiting water addition step is preceded by a quasi-equilibrated oxidation step. We formulated these design principles based on three observations – the existence of BEP relations relating the activation barrier of the water addition step to the free energy of this step, the existence of a linear scaling relation between O – H and O – OH bond dissociation energies at the active site, and a kinetic equivalence between electrochemical and non-electrochemical water addition.

41

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

For sites on both the (311) and (110)-A surfaces, it was found that the transition state for water addition corresponds to formation of the O – O bond in the resulting hydroperoxo intermediate, which occurs after a proton and electron have been completely transfered from the water to the active site to form a hydroxyl radical intermediate. Since the hydroxyl radical intermediate is formed prior to the transition state, it is in thermal equilibrium with the resting state of the active site. Therefore, we postulated that the pathway by which this proton and electron are transferred is not kinetically relevant – thus non-electrochemical water addition where both are transfered to the active site and electrochemical water addition where the proton and electron are transferred to the electrolyte and electrode, respectively, are kinetically equivalent. We chose to base our kinetic model on the latter as it is more intuitive. In contrast, sites on the (110)-B surface carry out water addition through a transition state that involves proton/electron transfer and O – O bond formation occurring concertedly. As such, electrochemical and non-electrochemical water addition pathways are not kinetically equivalent and we are only able to consider the latter. We obtained a BEP relation for the water addition step on each of the three surfaces. On the (311) and (110)-A surfaces, the BEP relation was written in terms of the electrochemical water addition pathway, while on the (110)-B surface it was written in terms of the nonelectrochemical pathway. Importantly, we found that the BEP relations exhibit a high degree of structure sensitivity between sites on the three different surface terminations that we examined. Sites on the (311) surface, which carry out water addition over a μ3 -oxo, have the most favorable BEP relation and give the largest water addition rate constant for a given value of the redox potential η1 . Sites on the (110)-A and (110)-B surfaces which carry out water addition over μ2 - and η-oxos, respectively, have less favorable BEP relations, with the latter having the least favorable. This is interesting, because the majority of computational catalyst screening studies for the OER on transition metal oxides only consider sites that carry out water addition to an η-oxo. The linear scaling relation that we found between O – H and O – OH bond dissociation

42

ACS Paragon Plus Environment

Page 42 of 49

Page 43 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

energies allowed us to express the water addition free energy in terms of the redox potential η1 of the pre-equilibrium oxidation step, as the former reaction involves O – OH formation while the latter involves O – H bond dissociation to the same oxygen. Combining this with the BEP relation for water addition and the rate expression for the catalytic cycle, we obtained an expression for the turnover frequency as a function of the two redox potentials η1 and η2 of the active site and the electrode potential η. The values of the redox potentials that maximize the TOF for a given electrode potential define the design criteria for the active site. The optimal site on the (311) and (110)-A surfaces has η1 = η and η2 > η, while the optimal site on the (110)-B surface has η1 = η2 = η. The origin of the η1 = η condition on all three surfaces is that the active site stores the energy supplied by the external potential during the oxidation step in the form of a high energy hole on the redox-active transition metal center. This energy is then released to extract an electron from the antibonding orbital of the O – O bond formed during the water addition step. When the redox potential of the oxidation step, η1 , is equal to the electrode potential, the site is able to store all of the supplied energy while still having a large population in the active oxidized state. Finally, we contrasted our kinetics-based catalyst design approach to the widely employed thermodynamics-based approach pioneered by Rossmeisl, Nørskov, and coworkers. The main difference is that the thermodynamics-based approach identifies the same optimal site independent of the operating electrode potential, while the kinetics-based approach identifies a different optimal site for every electrode potential. While both approaches identify the same optimal site for an electrode potential of 0.50 V, the kinetics-based approach is able to identify a more active catalyst for all electrode potentials above and below this. The difference is due to the fact that the kinetics-based approach correctly accounts for the change in the activation barrier for water addition as the free energy of that step changes. Since both approaches require the computation of only a single descriptor for each site – effectively the redox potential η1 – the kinetics-based approach is preferable as it is able to identify a more active catalyst with the same computational difficulty and additionally give an estimate of

43

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the absolute turnover frequency.

Supporting Information Available See Supporting Information for technical details of the DFT calculations, method for locating intersystem crossings, method for constraining atomic magnetic moments, procedure for checking the converged electronic state, calculation of free energy differences, derivation of the rate expression in eq (1), derivation of eq (17), activity plot for the (110)-A surface, structures of surfaces used in the calculations, tables of selected bond distances, and tables of calculated energies. This information is available free of charge via the Internet at http://pubs.acs.org.

Acknowledgement We gratefully acknowledge support from the Solar Technologies Go Hybrid initiative of the State of Bavaria and from the German Research Council (DFG) under grant number RE1509/33-1. We also thank the Leibniz Supercomputing Center for access to supercomputing facilities.

References (1) Ibrahim, H.; Ilinca, A.; Perron, J. Energy storage systems - Characteristics and comparisons. Renew. Sustain. Energy Rev. 2008, 12, 1221–1250. (2) Buttler, A.; Spliethoff, H. Current status of water electrolysis for energy storage, grid balancing and sector coupling via power-to-gas and power-to-liquids: A review. Renew. Sustain. Energy Rev. 2018, 82, 2440–2454. (3) Fabbri, E.; Habereder, A.; Waltar, K.; K¨otz, R.; Schmidt, T. J. Developments and per-

44

ACS Paragon Plus Environment

Page 44 of 49

Page 45 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

spectives of oxide-based catalysts for the oxygen evolution reaction. Catal. Sci. Technol. 2014, 4, 3800–3821. (4) Gal´an-Mascar´os, J. R. Water oxidation at electrodes modified with earth-abundant transition-metal catalysts. ChemElectroChem 2015, 2, 37–50. (5) Hong, W. T.; Risch, M.; Stoerzinger, K. A.; Grimaud, A.; Suntivich, J.; Shao-horn, Y. Toward the rational design of non-precious transition metal oxides for oxygen electrocatalysis. Energy Environ. Sci. 2015, 8, 1404–1427. (6) Schalenbach, M.; Tjarks, G.; Carmo, M.; Lueke, W.; Mueller, M.; Stolten, D. Acidic or alkaline? Towards a new perspective on the efficiency of water electrolysis. J. Electrochem. Soc. 2016, 163, F3197–F3208. (7) Burke, M. S.; Enman, L. J.; Batchellor, A. S.; Zou, S.; Boettcher, S. W. Oxygen evolution reaction electrocatalysis on transition metal oxides and (oxy)hydroxides: Activity trends and design principles. Chem. Mater. 2015, 27, 7549–7558. (8) Nørskov, J. K.; Bligaard, T.; Rossmeisl, J.; Christensen, C. H. Towards the computational design of solid catalysts. Nat. Chem. 2009, 1, 37–46. (9) Greeley, J. Theoretical heterogeneous catalysis: scaling relationships and computational catalyst design. Annu. Rev. Chem. Biomol. Eng. 2016, 7, 605–635. (10) Reuter, K.; Plaisance, C. P.; Oberhofer, H.; Andersen, M. Perspective: On the active site model in computational catalyst screening. J. Chem. Phys. 2017, 146, 040901. (11) Nørskov, J. K.; Rossmeisl, J.; Logadottir, A.; Lindqvist, L.; Kitchin, J. R.; Bligaard, T.; J´onsson, H. Origin of the overpotential for oxygen reduction at a fuel-cell cathode. J. Phys. Chem. B 2004, 108, 17886–17892. (12) Rossmeisl, J.; Qu, Z.-W.; Zhu, H.; Kroes, G.-J.; Nørskov, J. Electrolysis of water on oxide surfaces. J. Electroanal. Chem. 2007, 607, 83–89. 45

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(13) Man, I. C.; Su, H.-Y.; Calle-Vallejo, F.; Hansen, H. A.; Mart´ınez, J. I.; Inoglu, N. G.; Kitchin, J.; Jaramillo, T. F.; Nørskov, J. K.; Rossmeisl, J. Universality in oxygen evolution electrocatalysis on oxide surfaces. ChemCatChem 2011, 3, 1159–1165. (14) Mavros, M. G.; Tsuchimochi, T.; Kowalczyk, T.; McIsaac, A.; Wang, L.-P.; Voorhis, T. V.; Mc Isaac, A.; Wang, L.-P.; Van Voorhis, T. What can density functional theory tell us about artificial catalytic water splitting? Inorg. Chem. 2014, 53, 6386–6397. (15) Plaisance, C. P.; van Santen, R. A. Structure sensitivity of the oxygen evolution reaction catalyzed by cobalt(II,III) oxide. J. Am. Chem. Soc. 2015, 137, 14660–14672. (16) Plaisance, C. P.; Reuter, K.; van Santen, R. A. Quantum chemistry of the oxygen evolution reaction on cobalt(II,III) oxide implications for designing the optimal catalyst. Faraday Discuss. 2016, 188, 199–226. (17) Perdew, J. P.; Burke, K.; Wang, Y. Generalized gradient approximation for the exchange-correlation hole of a many-electron system. Phys. Rev. B 1996, 54, 16533– 16539. (18) Kresse, G.; Furthm¨ uller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B 1996, 54, 11169–11186. (19) Cohen, A. J.; Mori-Sanchez, P.; Yang, W. Insights into current limitations of density functional theory. Science (80-. ). 2008, 321, 792–794. (20) Anisimov, V. I.; Gunnarsson, O. Density-functional calculation of effective Coulomb interactions in metals. Phys. Rev. B 1991, 43, 7570–7574. (21) Liechtenstein, A. I. I.; Anisimov, V. I.; Zaanen, J. Density-functional theory and strong interactions: Orbital ordering in Mott-Hubbard insulators. Phys. Rev. B 1995, 52, R5467–R5470. 46

ACS Paragon Plus Environment

Page 46 of 49

Page 47 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(22) Anisimov, V. I.; Aryasetiawan, F.; Lichtenstein, A. I. First-principles calculations of the electronic structure and spectra of strongly correlated systems: the LDA + U method. J. Phys. Condens. Matter 1997, 9, 767–808. (23) Henkelman, G.; Uberuaga, B. P.; Jonsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J. Chem. Phys. 2000, 113, 9901. (24) Henkelman, G.; Jonsson, H. Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points. J. Chem. Phys. 2000, 113, 9978. (25) Henkelman, G.; Jonsson, H. A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives. J. Chem. Phys. 1999, 111, 7010. (26) Bearpark, M. J.; Robb, M. A.; Bernhard Schlegel, H. A direct method for the location of the lowest energy point on a potential surface crossing. Chem. Phys. Lett. 1994, 223, 269–274. (27) Allen, J. P.; Watson, G. W. Occupation matrix control of d- and f-electron localisations using DFT + U. Phys. Chem. Chem. Phys. 2014, 16, 21016–21031. (28) Qian, X.; Li, J.; Qi, L.; Wang, C.-Z.; Chan, T.-L.; Yao, Y.-X.; Ho, K.-M.; Yip, S. Quasiatomic orbitals for ab initio tight-binding analysis. Phys. Rev. B 2008, 78, 245112. (29) Plaisance, C. P.; van Santen, R. A.; Reuter, K. Constrained-orbital density functional theory. Computational method and applications to surface chemical processes. J. Chem. Theory Comput. 2017, 3561–3574. (30) Betley, T. A.; Wu, Q.; Van Voorhis, T.; Nocera, D. G. Electronic design criteria for OO bond formation via metaloxo complexes. Inorg. Chem. 2008, 47, 1849–1861.

47

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(31) Desai, S. K.; Neurock, M. First-principles study of the role of solvent in the dissociation of water over a Pt-Ru alloy. Phys. Rev. B - Condens. Matter Mater. Phys. 2003, 68, 1–7. (32) Tripkovi´c, V.; Sk´ ulason, E.; Siahrostami, S.; Nørskov, J. K.; Rossmeisl, J. The oxygen reduction reaction mechanism on Pt(111) from density functional theory calculations. Electrochim. Acta 2010, 55, 7975–7981. (33) Stecher, T.; Reuter, K.; Oberhofer, H. First-principles free-energy barriers for photoelectrochemical surface reactions: proton abstraction at TiO2 (110). Phys. Rev. Lett. 2016, 117, 276001. (34) Shaik, S. S. What happens to molecules as they react? A valence bond approach to reactivity. J. Am. Chem. Soc. 1981, 103, 3692–3701. (35) Pross, A.; Shaik, S. S. A qualitative valence-bond approach to organic reactivity. Acc. Chem. Res. 1983, 16, 363–370. (36) Koper, M. T. M. Theory of multiple proton-electron transfer reactions and its implications for electrocatalysis. Chem. Sci. 2013, 4, 2710. (37) Hansen, H. A.; Viswanathan, V.; Nørskov, J. K. Unifying kinetic and thermodynamic analysis of 2 e- and 4 e- reduction of oxygen on metal surfaces. J. Phys. Chem. C 2014, 118, 67066718. (38) Chen, S.; Wang, L.-W. Thermodynamic oxidation and reduction potentials of photocatalytic semiconductors in aqueous solution. Chem. Mater. 2012, 24, 3659–3666.

48

ACS Paragon Plus Environment

Page 48 of 49

Page 49 of 49 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Graphical TOC Entry Activity (311)

Optimal Site

𝜂1 redox Water Addition: Release Energy

e− 𝜂

𝜼𝟐 > 𝜼𝟏 = 𝜼 e− Oxidation: Store Energy

49

ACS Paragon Plus Environment

potential