The Sabatier Principle in Electrocatalysis: Basics, Limitations, and Extensions

The Sabatier principle, which states that the binding energy between the catalyst and the reactant should be neither too strong nor too weak, has been widely used as the key criterion in designing and screening electrocatalytic materials necessary to promote the sustainability of our society. The widespread success of density functional theory (DFT) has made binding energy calculations a routine practice, turning the Sabatier principle from an empirical principle into a quantitative predictive tool. Given its importance in electrocatalysis, we have attempted to introduce the reader to the fundamental concepts of the Sabatier principle with a highlight on the limitations and challenges in its current thermodynamic context. The Sabatier principle is situated at the heart of catalyst development, and moving beyond its current thermodynamic framework is expected to promote the identification of next-generation electrocatalysts.


INTRODUCTION
Electrocatalysis is gaining widespread attention as a critical technology to enhance the sustainability of human society. Core technology for the hydrogen economy, such as fuel cells and water electrolysis, rely on efficient electrocatalysts to perform the interconversion between water and hydrogen (Mehta and Cooper, 2003). The electrochemical synthesis of next-generation fuels, derived from atmospheric CO 2 rather than from fossil fuels, is also under intensive investigation . Even the carbon footprint of the artificial nitrogen cycle can be improved using electrocatalysis, because the Haber-Bosch process used to synthesize ammonia and chemical fertilizers requires vast amounts of hydrogen, which are mainly produced from fossil fuels (Smith et al., 2020). Therefore, the development of electrocatalysts, especially those based on earth-abundant elements, is currently a central topic in chemistry and energy research.
However, the parameter space to search for new catalytic materials is massive. For example, oxides, such as the dimensionally stable anodes in the chlor-alkali industry (Karlsson and Cornell, 2016) or the iridium oxides in electrolyzers (Carmo et al., 2013), are a widely used family of catalysts. Based on the periodic table of elements, a chemist may choose from approximately 40 elements with which to make a single metal oxide. However, if the material candidates are expanded to binary, ternary, or even more complex mixed oxides, the possibility of materials expands massively (40 4 > 2 million), especially if dopants are considered. Furthermore, every elemental composition has multiple possible atomic orientations in the bulk and the surface (Ulissi et al., 2017). Surface reconstruction, which can potentially occur during catalysis, are an additional contribution to the huge parameter space of materials (Fabbri et al., 2017;Li et al., 2019), leading to a combinatorial explosion of the possible candidate materials. In order to identify next-generation catalysts within this massive parameter space, a design principle which connects the chemical structure to catalytic activity is needed.
Nowadays, this role is fulfilled by the Sabatier principle (Sabatier, 1913;Che, 2013). This principle states that "an ideal catalyst must bind to the reactant at an intermediate strength which is neither too weak nor too strong." This concept is based on the notion that if the bond is too weak, the catalyst and the reactant will hardly interact with each other, whereas if the bond is too strong, the reactant will not desorb from the catalyst surface, effectively inhibiting further reactions. This principle was first proposed in 1913 by the Nobel laureate, Paul Sabatier, based on empirical observations (Sabatier, 1913;Che, 2013). At that time, there was no method to directly obtain the binding energy, and therefore, other experimentally accessible parameters, such as hydride (Trasatti, 1972) or chelate formation energies (Rootsaert and Sachtler, 1960), were used to quantitatively assess the interaction between the catalyst and the reactant. A century later, however, it is now the norm to directly evaluate the binding energy of intermediate states using in silico computation Rossmeisl et al., 2005;Rossmeisl et al., 2007b). In particular, the widespread accessibility of density functional theory (DFT) has made it almost every-day practice to predict the activity of a potential candidate prior to performing actual experiments (Greeley et al., 2006). Recent computational studies are further supplemented by machine-learning techniques to reduce computational costs (Back et al., 2019a,b;Li et al., 2020;Ulissi et al., 2016Ulissi et al., , 2017. In the following section, we will provide a brief, conceptual overview of the Sabatier principle in its current form, which has been widely employed to evaluate catalysts in silico. In section "CHALLENGES OF THE SABATIER PRINCIPLE AND APPROACHES FORWARD, " its challenges and limitations will be highlighted, along with ongoing approaches to overcome these bottlenecks. The main question addressed in this review is "what is the (theoretical) requirement for active electrocatalysts?", not "which material satisfies these conditions?" For the reader with an interest in the progress of state-of-the-art materials, we would like to recommend several excellent reviews published recently (Xiao et al., 2015;Tahir et al., 2017;Song et al., 2020;Theerthagiri et al., 2020).

THE SABATIER PRINCIPLE AND ITS THERMODYNAMIC INTERPRETATION
The Ideal Free-Energy Landscape The free-energy difference between the initial (reactant) and final (product) states is defined by thermodynamics and is independent of the actual catalyst. When the applied electrode potential (U) is the same as the half-cell potential U 0 , the free energies of the reactant and product are equal, and the system is in equilibrium. The reaction can proceed once the electrode potential is shifted away from the half-cell potential, such that the free energy of the product state becomes lower than that of the reactants. The difference between the applied potential and the half-cell potential is known as the applied overpotential, η = U -U 0 . A positive overpotential is required for anodic (oxidation) processes, such as the chlorine evolution reaction (CER), whereas a negative overpotential is necessary for cathodic (reduction) reactions, such as the hydrogen evolution reaction (HER). However, the word "overpotential" is used to refer to absolute values in the following sections for simplicity, and hence "a larger overpotential" should be interpreted as the reaction becoming more thermodynamically favorable.
Within this framework, the role of the catalyst is to tune the free-energy landscape between the reactant and the product in a way that is beneficial for the catalytic process (Koper, 2011b). For a two-step reaction, such as the CER (Trasatti, 1987;Exner, 2020b) or HER (Gerischer and Gerischer, 1956;Parsons, 1958;Nørskov et al., 2005;Zeradjanin et al., 2016;Exner, 2020e), there is only one reaction intermediate, whose free energy depends on the catalyst surface to which it is adsorbed to (Figure 1A; Koper, 2011bKoper, , 2013a. The free energy of this reaction intermediate with respect to the reactant at equilibrium potential is denoted as G RI . Currently, the general consensus of an ideal free-energy landscape is one where this reaction intermediate has the same free energy as the reactant and the product at equilibrium. In other words, G RI = 0 is the criterion for an ideal catalyst (Koper, 2011b(Koper, , 2013aLaursen et al., 2012), and G RI can be considered as a "descriptor" for catalytic activity (Greeley and Nørskov, 2007;Peterson and Nørskov, 2012;Dubouis and Grimaud, 2019).
The reason behind why thermoneutral bonding ( G RI = 0) is considered ideal is shown in Figure 1. If the reaction intermediate is bound too weakly ( G RI > 0), the first step is thermodynamically unfavorable (Figure 1A, blue lines). On the other hand, if the intermediate is bound too strongly ( G RI < 0), the second step is thermodynamically unfavorable (Figure 1B, green lines). If we assume that the efficiency of the overall reaction is determined by the most thermodynamically unfavorable step, a thermoneutral landscape where no elementary step is thermodynamically unfavorable corresponds to the ideal situation ( Figure 1A, red lines). This is often discussed quantitatively using the concept of the thermodynamic overpotential (η TD ), introduced by Nørskov, Rossmeisl, and co-workers . η TD is defined as the minimum overpotential necessary to make all elementary reaction steps either thermoneutral or exergonic (Rossmeisl et al., 2005;Rossmeisl et al., 2007b). A smaller value of η TD implies higher activity, and η TD is commonly used to compare and rationalize the activity of various materials. In the case of a two-step reaction, η TD = | G RI |/e with e = elementary charge, and therefore, a thermoneutral landscape ( G RI = 0) yields η TD = 0 V. For this reason, the thermodynamically ideal catalyst is shows the free-energy landscape for a catalyst that binds the reaction intermediate weakly ( G RI > 0 at equilibrium, black). By applying an overpotential of η = η TD , all elementary reaction steps become thermoneutral or downhill in free energy (red). sometimes referred to as the "zero overpotential catalyst" (Koper, 2011b(Koper, , 2013aCraig et al., 2019). This discussion is based on purely thermodynamic considerations and does not explicitly differentiate between outer-sphere and inner-sphere electron transfer mechanisms. However, all electrocatalytic reactions involve the rearrangement of chemical bonds, suggesting they can be considered as inner-sphere processes (Bard, 2010;Schmickler and Santos, 2010).
Nowadays, the free-energy landscape of an electrocatalyst can be calculated using computational techniques, such as density functional theory (DFT). Nørskov and coworkers introduced the computational hydrogen electrode (CHE) model, which substitutes the free energy of an electron-proton pair with that of half a hydrogen molecule at U = 0 V vs. RHE (reversible hydrogen electrode) . Using this method, it is possible to determine the free energies of reaction intermediates as a function of the applied electrode potential and pH by a posteriori analysis. Once the free energies of reaction intermediates are determined by DFT at zero electrode potential and pH zero conditions, they can be shifted to any electrochemical condition based on the CHE model (Calle-Vallejo and Koper, 2012). In principle, this method can also be used to determine the free energies of charged reaction intermediates (Craig et al., 2019), although in practice, most studies focus on non-charged intermediates (Calle-Vallejo and Koper, 2012) due to the large computational errors for charged intermediate states on solid surfaces. Regardless of these limitations, the CHE model and the thermodynamic overpotential (η TD ) are by far the most popular way to evaluate electrocatalytic activity in silico. To date, many experimental trends have been rationalized and new materials have been identified based on this simple framework (Hinnemann et al., 2005;Greeley et al., 2006;Diaz-Morales et al., 2016;Seitz et al., 2016).

Practical Application and Linear Scaling Relationships
From a thermodynamic viewpoint, the free-energy landscape should be completely flat at equilibrium to obtain the ideal catalyst. This notion holds true even when the mechanism involves more than two elementary steps. In practice, however, it is much easier to tune the landscape for two-electron transfer reactions, such as the HER and CER, compared to reactions which transfer more electrons, such as the oxygen evolution and reduction reactions (OER and ORR: 2 . This is reflected in the experimental literature, as some HER/HOR catalysts, such as Pt, are so active that practically no overpotential is visible in the cyclic voltammogram and the current density increases exponentially as soon as an overpotential is applied (Auinger et al., 2011;Durst et al., 2014). Other examples include a Pd-based catalyst which also exhibits almost no overpotential to catalyze the two-electron transfer reaction between CO 2 and formic acid (Kortlever et al., 2015a). In contrast, even the most active OER/ORR catalysts require overpotentials of several hundred mV to obtain current densities in the order of mA/cm 2 (Suntivich et al., 2011a,b;Diaz-Morales et al., 2016;Seitz et al., 2016;Chen et al., 2017). The reason why it is difficult to realize thermodynamically ideal catalysts for reactions with more than two electrons transferred has been rationalized due to the socalled scaling relationships, which impose physical limitations on how the free-energy landscape can be tuned. Let us take the OER as an example.
The OER is a four-electron process, in which water molecules are transformed into gaseous oxygen according to equation (1): Rossmeisl et al. (2005Rossmeisl et al. ( , 2007b calculated the OER free-energy landscape according to a pathway, which assumes the OH, O, and OOH adsorbates as reaction intermediates. Each step was assumed to proceed via concerted proton-electron transfer steps. This was originally due to technical reasons, namely to maintain charge neutrality (cf. Section "Stepwise Proton-Electron Transfer and Charged Intermediates") (Rossmeisl et al., 2005(Rossmeisl et al., , 2007b, but this mechanism is still adopted in the recent literature (Halck et al., 2014;Busch et al., 2016;Dickens and Nørskov, 2017;Hajiyani and Pentcheva, 2018).
In equations (2) -(5), S denotes the active site of an OER catalyst. In order to realize an ideal catalyst, every step must be thermoneutral at the equilibrium potential. In other words, when the free-energy change of the i th step is denoted as G i , the criterion of a thermodynamically ideal catalyst is: G 1 = G 2 = G 3 = G 4 = 0 eV at U = 1.23 V vs. RHE (Figure 2, blue lines). This results in G 1 = G 2 = G 3 = G 4 = 1.23 eV at U = 0 V vs. RHE (Figure 2, black lines), which corresponds to the usual conditions of DFT calculations within the CHE framework.
In reality however, no electrode material that fulfills the requirements of a thermodynamically ideal catalyst has been identified so far. This is because the binding energies of structurally similar intermediates are difficult to tune independently. For example, the binding energies of the OH and OOH adsorbates are known to exhibit a linear correlation. This relationship is often referred to as a linear scaling relationship (Man et al., 2011;Viswanathan et al., 2012), and is given by equation (6): The number 3.2 eV is the value that has been reported for planar metal oxides (Rossmeisl et al., 2005;Man et al., 2011;Viswanathan et al., 2012). The values in more recent studies span a range of about 2.8 eV -3.4 eV (Calle-Vallejo et al., 2013;Viswanathan and Hansen, 2014;Tao et al., 2019;Exner, 2020a), suggesting that different classes of materials may follow different scaling relationships.  Due to the linear scaling relationship between the OH and OOH adsorbates ( G 2 + G 3 = 3.2 eV), it is not possible to adjust both G 2 and G 3 to 1.23 eV simultaneously (blue and green pathways). The minimum thermodynamic (theoretical) overpotential is realized when G 2 = G 3 = 1.6 eV; that is, the free-energy penalty due to the scaling relations is shared equally between these two elementary steps. (B) The thermodynamic overpotential is plotted as a function of G 2 in the form of a volcano plot.
In any case, the scaling relationship directly implies that a thermodynamically ideal catalyst cannot be realized. For example, if a certain electrocatalyst exhibits ideal thermodynamics for the second step ( G 2 = 1.23 eV at U = 0 V), the scaling relationships indicate that the third step cannot be thermodynamically ideal due to: G 3 = 3.2 eV -1.23 eV = 1.97 eV ( Figure 3A, green line). On the other hand, setting G 3 to 1.23 eV leads to a scenario where G 2 is unfavorable ( Figure 3B, blue line). In other words, optimizing either G 2 or G 3 leads to the other deviating from the thermodynamic ideal. If the most thermodynamic unfavorable reaction step dictates catalytic activity, the best compromise is to distribute the free-energy penalty from the scaling relationship into equal parts; that is, G 2 = G 3 = 3.2 eV / 2 = 1.6 eV (Man et al., 2011). In this case, the thermodynamic overpotential amounts to η TD = 1.60 V -1.23 V = 0.37 V ( Figure 3A, red line), which is significantly smaller than η TD = 1.97 V -1.23 V = 0.74 eV which results by making either the second step or third step thermodynamically ideal ( Figure 3A, blue and green lines). In terms of the volcano plot, electrocatalysts with optimum G 2 or G 3 are located at the leg of the volcano because they bind the O adsorbate either too strongly or too weakly, respectively, such that the other step becomes thermodynamically prohibitive. This leads to a minimum thermodynamic overpotential, also denoted as minimum theoretical overpotential, of 0.37 V for the OER (Figure 3B). This value is roughly consistent with the experimental overpotential necessary to reach a current density of 10 mA/cm 2 , a benchmark for the solar-cell community (Suntivich et al., 2011b;Diaz-Morales et al., 2016;Seitz et al., 2016). However, some recent materials have overpotentials below 200 mV to sustain of 10 mA/cm 2 Su et al., 2018;Retuerto et al., 2019;Hao et al., 2020), indicating that η TD is not necessarily a quantitatively accurate descriptor of electrocatalytic activity (Exner, 2021a). This discrepancy is most likely because the thermodynamic overpotential and the experimental overpotential have different definitions, and thus have no direct physical correlation: while η TD is a purely thermodynamic indicator of activity based on the free-energy landscape at equilibrium (zero overpotential), the experimental overpotential η must be evaluated at some non-zero current density, and is thus always influenced by the kinetics. However, the concept of η TD does offer at least a qualitative explanation why the OER is a more difficult process compared to the HER and other two-electron processes, where the thermodynamics is not restrained by any scaling relationships.

Examples for the Application of the Sabatier Principle in Practice
The Sabatier principle and its thermodynamic interpretation have been successfully applied for the development of new catalysts. For example, Greeley et al. used G RI to screen alloy catalysts for the HER (Greeley et al., 2006), and identified a Bi-Pt surface alloy with higher electrocatalytic activity than platinum. Similarly, Hinnemann et al. (2005) identified MoS 2 nanoparticles as a potential HER electrocatalyst by applying the Sabatier principle in its thermodynamic form. The work of Seitz et al. is another example, where the highly active IrO x /SrIrO 3 catalyst for the OER was predicted based on the framework of scaling relations and η TD discussed in the previous section (Seitz et al., 2016). These works, among others, contributed to the fact that in silico screening procedures became extremely popular in electrocatalysis and, nowadays, are broadly used for the discovery of catalytic materials (Yao et al., 2019;Yang et al., 2020).
More recently, the rise of data science has promoted the usage of the Sabatier principle for materials' screening (Basdogan et al., 2020a;Toyao et al., 2020). For example, Ulissi et al. used a neural network-based surrogate model to reduce the number of explicit DFT calculations by an order of magnitude, allowing them to extensively screen active sites for CO 2 reduction (Ulissi et al., 2017).
Machine-learning has also been beneficial to obtain more accurate binding energies across wider configurations (Back et al., 2019a).

Application of the Sabatier Principle to Enzymes, Molecular Catalysts, and Beyond
The Sabatier principle has also been successfully applied to enzymes and homogeneous catalysts. For example, biological OER and ORR enzymes were compared with their artificial counterparts, and their higher activity was rationalized by their lower η TD (Rossmeisl et al., 2007a;Kjaergaard et al., 2010). These enzymes deviate from the linear scaling relationships reported for planar metal oxides (Man et al., 2011). Therefore, this finding has stimulated various studies aiming to find strategies to break scaling relations (Li and Sun, 2016;Pegis et al., 2017;Huang et al., 2019). However, breaking the scaling relationship in itself is not a sufficient condition to realize active materials, because it does not guarantee a lower η TD as the recent literature has shown (Govindarajan et al., 2018;Exner, 2020f;Piqué et al., 2020).
More recently, Kari et al. (2018) demonstrated experimentally that enzymes also follow a Sabatier-type trend. The Michaelis-Menten constant K m is a quantitative measure of the substrateenzyme interaction and was used as a substitute for the binding energy to construct volcano plots of various cellulases. The maximum rate was observed at intermediate values of K m , suggesting that an intermediate interaction strength between the catalyst and the reactant may also be important to enhance enzymatic activity.
Scaling relations and volcano plots have also been introduced into the field of organometallic chemistry. For example, the thermodynamics of Suzuki coupling catalysts was rationalized based on linear scaling relationships and free-energy changes (Busch et al., 2015(Busch et al., , 2016Wodrich et al., 2018). Catalysts known to be active from experiments appear near the top of the volcano plot, suggesting the Sabatier principle may also help in the understanding of organometallic catalysts.
The Sabatier principle was also transferred to the field of lithium-ion batteries in order to study the intercalation of lithium-ions into spinel lithium titanate electrodes (Exner, 2018. While the conventional Sabatier principle only considers activity by analyzing binding energies, stability was introduced as a second performance parameter to rationalize experimental trends by the construction of activity-stability volcano plots. The volcano plots indicate that an intermediate lithium binding strength is the best compromise between enhanced energy density while maintaining stability, thus building a bridge between the catalysis and battery science communities (Exner, 2018. This may be the reason why the idea of an intermediate binding strength in conjunction with the volcano notion has found entrance for the screening of electrode materials in batteries (Pande and Viswanathan, 2019). The advanced framework of activity-stability volcano curves can also be employed to electrocatalytic processes such as the CER (Exner, 2019a).

CHALLENGES OF THE SABATIER PRINCIPLE AND APPROACHES FORWARD
The simplicity of the Sabatier principle and its thermodynamic interpretation has led to its widespread application so far. At the same time, however, this simplicity implies limitations where theory does not match with experiments. This can be due to issues at multiple levels.
For example, the binding energy estimated from the CHE model may not reflect the true binding energy without a more extensive treatment of factors such as solvation or the effective electrochemical potential in the double layer. In this case, the discrepancy between theory and experiments is a matter of accuracy which can be resolved by sufficient technical advancement. The Sabatier principle itself is not the problem if a more accurate binding energy can yield predictions consistent with experiments.
In other cases, the issue may be rooted deeper in the theory. For example, the Sabatier principle has its basis on the thermodynamics because the binding energy is in essence an equilibrium constant between the adsorption and desorption of a reactant. However, thermodynamics alone cannot explain how quickly electron-transfer reactions occur. Conventionally, kinetic factors such as activation barrier heights are assumed to scale universally with the thermodynamics. However, if this assumption is not valid, activity trends can no longer be predicted from thermodynamic criteria such as the thermodynamic overpotential, binding energies, or scaling relations. This is an issue which is beginning to attract more attention, stimulating recent microkinetic studies (Ooka and Nakamura, 2019;Exner, 2020c,e) which have attempted to expand the Sabatier principle (Exner, 2019d(Exner, , 2021b by the incorporation of overpotential and kinetic factors. Finally, there are some chemical processes which the Sabatier principle does not consider at all. The Sabatier principle predicts the intrinsic activity of a material, but in some cases, the observed electrocatalytic behavior is not purely due to the intrinsic activity. For example, mass transport is well known to influence electrode kinetics, such as in the case of the CO 2 reduction reaction (Ringe et al., 2020). The electrochemically active surface area is another consideration which cannot be determined by DFT calculations due to the usage of simplified slab models.
Proper comparison between theoretical and experimental activities requires calculation of the full free-energy landscape including the kinetics, which is computationally expensive. This is the reason why theoretical and experimental activities are often evaluated at the level of material trends by comparing the thermodynamics (theory) with the kinetics (experiment). A more rigorous and critical verification based on a direct quantitative comparison would be beneficial to improve the accuracy of theoretical predictions. In the following sections, we have attempted to sort the existing challenges of the Sabatier principle and the CHE model in order to provide a guideline of possible approaches beyond the current thermodynamic framework.

Influence of the Electrolyte: Beyond the CHE Model
Electrocatalysis is in general a phenomenon which is also dependent on the electrolyte. The most famous example for pH dependence may be the HER on platinum, which is more active in acid compared to alkaline conditions (Sheng et al., 2010). However, the CHE model treats the electrolyte ion concentrations of protons and hydroxides based on the Nernst equation, and thus, the thermodynamic overpotential η TD is independent of the pH. As η TD is the common theoretical measure for electrocatalytic activity, this is a direct contradiction with experiments. Other ions besides protons and hydroxides in the electrolyte may also influence electrocatalytic activity Garcia et al., 2019;Shinagawa et al., 2019), yet such factors are not considered in the CHE model.
Although the influence of the electrolyte can be ascribed to multiple factors, at least part of the reason is due to the binding energy being electrolyte dependent. The electrolyte determines the double-layer structure surrounding the adsorbate (Ledezma-Yanez et al., 2017), which has direct implications on the stability of intermediates on the electrode.
The double layer structure is also influenced by the electrochemical potential, and thus the binding energy may also deviate from the Nernstian behavior assumed in the CHE model when the electrode potential or pH is changed (Hörmann et al., 2019(Hörmann et al., , 2020. There are several recent studies which have attempted to explain the pH dependence of Pt in the HER by considering the binding energies to be pH dependent (Kunimatsu et al., 2007;Sheng et al., 2010Sheng et al., , 2015Zheng et al., 2016;Ledezma-Yanez et al., 2017;Zhu et al., 2020). Yet, such effects are not considered in most theoretical studies, which tacitly assume the binding energy to be a material-specific value which is pH independent (Nørskov et al., , 2009. The above-highlighted challenges could potentially be overcome if the electrolyte solution, consisting of all ions and solvent molecules, is explicitly modeled in the theoretical framework. This, however, requires computationally expensive ab initio molecular dynamics (AIMD) simulations (Kenmoe et al., 2018;Sakong and Groß, 2020). The cell sizes that can be treated by AIMD (Groß, 2020) so far are too small to simulate the dependency of binding energies on electrolyte concentrations (Groß, 2020). Although a steep increase in the computational power is needed to accurately determine the binding energies as a function of electrolyte composition, the effect of the electrolyte solution on the binding energy of an intermediate is significant, nonetheless. For example, explicit consideration of a solvent bilayer was found to change the DFT-obtained binding energy of adsorbed hydroxyl (OH ads ) on Pt(111) by −0.5 eV, and assuming the same solvation configuration for different metals resulted in errors exceeding 0.2 eV (He et al., 2017). This is the same order of magnitude compared to typical DFT errors (Nørskov et al., 2009), indicating that proper modeling of solvation within the available methods is indispensable for the computational design of electrocatalysts.
Currently, continuum solvation models such as Vaspsol (Mathew et al., 2014(Mathew et al., , 2019, COSMO (Klamt and Schüürmann, 1993), or Jaguar (Bochevarov et al., 2013) are commonly used to describe solvation at the electrode/electrolyte interface (Abidi et al., 2020;Basdogan et al., 2020b;Groß, 2020). However at least in the case of the intermediates relevant to the ORR and CO 2 reduction reactions on Cu, Au, and Pt, the widely used continuum solvation models did not improve the accuracy of binding energies compared to the respective calculations in vacuum (Heenen et al., 2020). Due to the complexity and computational costs of AIMD, solvation models by the inclusion of machine-learning techniques have been developed (Basdogan et al., 2020b), which are a promising alternative to continuum solvation models.
Another shortcoming of the CHE model refers to the fact that this approach relies on a canonical ensemble, in which the charge is fixed. Experimental conditions, however, correspond to a grand canonical ensemble where the electrode potential is fixed instead of the charge. Recent progress in this field put forth grand canonical models beyond the CHE method (Hansen and Rossmeisl, 2016;Hörmann et al., 2019Hörmann et al., , 2020Abidi et al., 2020), allowing the electrode potential to be explicitly included into the DFT calculations (Govind Rajan et al., 2020;Groß, 2020). Yet, a limitation of this approach is that the outcome may severely depend on the chosen value for the proton/hydrogen work function, for which values between 3.9 eV and 4.7 eV have been reported in the literature (Busch et al., 2020;Bhattacharyya et al., 2021).

Influence of Surface Charge
Another challenge in the accurate evaluation of the binding energy is the influence of the surface charge density on the electrode. Adsorbates on the electrode surface are usually charged due to the partial charge transfer in the chemisorption process. For instance, halogen ions adsorbed to a metal surface still carry a significant fraction of its negative charge after adsorption (Ávila et al., 2020) leading to the formation of surface dipole moments. The free energy of these surface dipoles depends on the local electric field, which is further modulated by the electrode surface charge. Therefore, the binding energy of reaction intermediates usually depends on the electrode surface charge. Furthermore, the electrode surface charge determines the local reaction conditions, such as the reactant concentration and the driving force at the reaction plane, thus having a pronounced effect on electrocatalytic activity. A recent example is the peroxodisulfate reduction on Pt, which has been shown to be inhibited when the surface charge density is negative (Martnez-Hincapié et al., 2018). This result could be explained by incorporating effects of the electrochemical double layer on electron transfer kinetics described using the Marcus-Hush-Chidsey theory (Zhang and Huang, 2020).
In order to incorporate double layer effects, an explicit model of the interface between the electrode and the electrolyte is necessary. When a reaction intermediate is adsorbed to the electrode surface, this is usually concurrent with the desorption of a solvent molecule. The displacement of the solvent molecule is rarely modeled explicitly due to the fact that conventional DFT studies are commonly performed in vacuum or with implicit solvation models. However, it has been proposed that this process is a possible origin of the pH-dependent HER activity of Pt (Sheng et al., 2015;Zheng et al., 2016;Ledezma-Yanez et al., 2017). If the chemisorbed intermediate possesses a dipole moment such as OH ads , this intermediate contributes to the electrical field at the electrode surface (Ávila et al., 2020). The electrical field is dependent on the coverage of the intermediates, which is in turn dependent on the electrochemical potential, thus making the modeling of the double layer's potential dependency a nontrivial task.
The normal potential of zero charge (PZC) is defined for ideally polarizable electrodes, at which variation of the metal potential only alters the distribution of the electron spillover without resulting in electron transfer from or onto ions and molecules in the solution. This is not the case for electrocatalytic interfaces, at which ions and molecules may adsorb onto the electrode surface, forming chemical bonds. The chemisorption dramatically changes the structure of the electric double layer. One known effect is associated with surface dipole moments that arise from partially charged chemisorbates. The chemisorptioninduced surface dipole moments generate a significant potential drop on the electrode surface. This additional surface potential drop will change the surface charging behavior of the electrode, resulting in a second PZC in the potential region where chemisorption occurs.
In the case of Pt (111) in 0.1 M HClO 4 , the normal PZC is around 0.3 V vs RHE, and shifting the electrode potential to more positive potentials should lead to an increase in the surface charge ( Figure 4A). However, at potentials above 0.6 V vs. RHE, water dissociation leads to the chemisorption of negatively charged OH − , forming an inward surface dipole moment ( Figure 4B, the direction of a dipole moment is from negative to positive charge). This effectively neutralizes the potential-dependent increase of the surface charge, resulting in a nonmonotonic surface charging relation as can be seen from the second PZC in the OH ads potential region (Huang et al., 2016b;Huang et al., 2018).
In general, when there is a net dipole change upon displacing the solvent with an adsorbate, there is a possibility for a second PZC, and therefore, this phenomenon is most likely not unique to the combination of platinum and OH intermediates (Frumkin and Petrii, 1975;Martnez-Hincapié et al., 2018). In addition to the ORR (Huang et al., 2016a;Zhang et al., 2019;Zhou et al., 2020), the double layer model above has been used successfully to model the formic acid oxidation (Zhu and Huang, 2019), and peroxodisulfate reduction (PDSR) reactions (Zhang and Huang, 2020) at Pt(111).

Stepwise Proton-Electron Transfer and Charged Intermediates
Another challenge to obtain accurate binding energies is the involvement of charged intermediate species. The CHE model avoids charged intermediates by assuming that the proton and electron are always transferred together in a single step, also denoted as concerted proton-electron transfer (CPET). This is partially justified because the CPET mechanism is in general thermodynamically more favorable compared to stepwise proton-electron transfer (SPET) (Hammes-Schiffer and Stuchebrukhov, 2010;Warren et al., 2010;Koper, 2013c). However, there is experimental evidence suggesting decoupled electron-proton transfer steps are involved in reactions such as CO 2 reduction on copper (Kortlever et al., 2015b) and gold (Wuttig et al., 2016), or the ORR on gold (Rodriguez and Koper, 2014). If so, this implies that some intermediates are charged, and evaluating their free energies accurately is necessary to construct the free-energy landscape for the electrocatalytic process. In fact, from a purely thermodynamic standpoint, the pH dependence of many electrocatalytic reactions already suggests the participation of SPET within the reaction mechanism (Hammes-Schiffer and Stuchebrukhov, 2010; Koper, 2013b;Sakaushi et al., 2020). Furthermore, the theoretical model by Koper suggests SPET may even be beneficial in some cases because modulating the driving force of electrochemical and non-electrochemical steps separately may lead to higher activity or selectivity (Koper, 2013b). The pH dependent selectivity regulation of nitrite reduction on MoS 2 is a direct support to this theoretical model (He et al., 2018. SPET pathways have also been proposed for other electrocatalytic processes, including the OER and ammonia oxidation reaction (Katsounaros et al., 2016;Saveleva et al., 2018;Exner and Over, 2019).
To date, the energetics of SPET pathways has been assessed theoretically for some molecular catalysts, and it has been shown that charged intermediates are indeed preferred during reactions such as CO 2 reduction (Göttle and Koper, 2017) or the OER (Craig et al., 2019). Thorough treatment of CPET and SPET pathways by quantum chemical considerations can be found in the works of Hammes-Schiffer and coworkers, who have developed computational protocols for the calculation of concerted or decoupled proton-electron transfers for biological enzymes and bio-inspired molecules (Goings and Hammes-Schiffer, 2020;Sayfutyarova and Hammes-Schiffer, 2020;Warburton et al., 2020). Although direct calculation of charged intermediates using DFT tends to be avoided due to their large errors, moving beyond the conventional picture of uncharged intermediate states within the CHE model is of particular importance to the electrocatalysis community. First steps for molecular catalysts have been made, whereas the description of charged species on solidstate surfaces is still in its early stages (Kastlunger et al., 2018;Lindgren et al., 2020).

Coverage Dependence of the Binding Energy
When evaluating the binding energy of intermediate species, it is also important to account for the chemical environment of the adsorbate. For example, the hydrogen binding energy on Pt(111) was shown to shift from strong-binding to weak-binding when the coverage was increased (Skúlason et al., 2010). The hydrogen binding energy also increases with the coverage of CO on metals such as Cu or Mo, which has been used to explain why Cu favor CO 2 reduction of HER when the two are in competition (Zhang et al., 2014). The CER on RuO 2 (110) is another example which is influenced by the coverage (Exner et al., 2016). The active site under typical CER conditions consists of a surface oxygen atom, O ads , on which chloride anions from the electrolyte solution adsorb under formation of an OCl ads precursor (Exner et al., 2014). The binding strength of chlorine to the underlying oxygen atoms is stronger if there is a neighboring OH group instead of an oxygen atom, due to the attractive interaction between the adsorbed chlorine and the neighboring hydrogen atom in the OH group.
In order to calculate the binding energy more accurately, it is desirable to resolve not only the active site but also its surrounding environment under typical reaction conditions. One approach is to use DFT to construct surface Pourbaix diagrams (Hansen et al., 2008;Sumaria et al., 2018;Vinogradova et al., 2018). These diagrams indicate the thermodynamically most favorable surface structure as a function of the applied electrode potential and pH (Gossenberger et al., 2020), and by assuming that this surface corresponds to the active surface, it is possible to predict the coverage of intermediates and adsorbates. However, this approach is purely thermodynamic in nature, and cannot be used if the surface structure in experiments deviates from the most energetically favorable one. The gap between thermodynamics and kinetics is discussed in detail in the following section, and some of the methods discussed therein may allow the construction of more accurate Pourbaix diagrams which include not only thermodynamic but also kinetic information.

Thermodynamics and Kinetics
The BEP Relationship The previous section has focused on challenges to obtain accurate binding energies. However, even if the free energy of reaction intermediates could be obtained without any error, there is still a conceptual gap between the Sabatier principle and the activity in real experiments. This is because the traditional interpretation of the Sabatier principle is based on pure thermodynamic considerations, whereas experiments always contain the influence of kinetics.
A framework based on pure thermodynamic considerations has historically been rationalized by the Brønsted (Bell) -Evans -Polanyi (BEP) relationship, which states that the free-activation energy ( G RI # ) is linearly correlated with the reaction free energy ( G RI ) Cheng et al., 2008). In other words, G RI # = β G RI , where β denotes the so-called BEP coefficient (0 < β < 1) (Van Santen et al., 2009).
An example for a two-electron reaction, such as the HER, is shown in Figure 5. For a weak-binding catalyst (Figure 5A), the first step should have a higher transition-state free energy compared to the second because the first step is uphill in free energy. In contrast, for a strong-binding catalyst (Figure 5B), the BEP relation predicts that the second elementary reaction step would have a higher transition-state free energy than the first step. The optimum situation corresponds to a thermoneutral catalyst ( Figure 5C) because in this case, the transition-state free energies of both elementary reactions are at the same height, and the highest transition state within the entire landscape would be lower than that of the weak-binding and strong-binding catalysts. In general, the positive nature of β ensures that a thermodynamically unfavorable step would have unfavorable kinetics due to a higher transition state. Based on this logic, the thermodynamically ideal (zero-overpotential) catalyst should also be the kinetically optimum catalyst, in agreement with the discussion in Section 2.1. Following the pioneering work of Parsons from the 1950s (Parsons, 1951(Parsons, , 1958, the reaction intermediates preceding the highest transition state can be considered to be in quasi-equilibrium, and the rate of the entire reaction (Figure 5) is governed by the free energy difference between the intermediate with smallest free energy and the transition state with highest free energy .
On the other hand, there are examples in the literature where the binding energy and the activity do not seem to correlate. For example, thermodynamic DFT calculations analyzed in terms of the BEP relationship for the ORR predict an activity trend of Pt > Ir > Rh, although the experimental trend is Pt > Rh > Ir (Figure 6A; Zhou et al., 2020). This discrepancy persists despite the experiments being performed using single crystals and even after diffusion effects are compensated, suggesting FIGURE 5 | The free-energy landscape of a two-step electrocatalytic process in equilibrium. The heights of transition states indicated by # were drawn assuming a constant BEP coefficient: (A) weak-binding catalyst; (B) strong-binding catalyst; (C) thermoneutral catalyst. The highest transition state, shown in red, occurs at the reaction step with the least favorable thermodynamics. Therefore, the thermodynamically ideal catalyst (C) should also be the kinetically ideal catalyst. that the binding energy alone cannot predict the activity of electrocatalysts. A discrepancy between theory and experiment is observed even within different single crystal facets of Pt ( Figure 6B; Gómez-Marín et al., 2014). There are also examples where the activity trend is reversed when the overpotential is increased; that is, a material which is less active than another electrocatalyst at low overpotentials can become more active at large overpotentials (Exner, 2019d). These discrepancies are mainly because the experimental activity is evaluated by the current density (kinetics) whereas theoretical trends are evaluated based on the thermodynamics (cf. Section "The Ideal Free-Energy Landscape") (Exner, 2021b).

Validity of the BEP Relationship
First, it is important to note that the BEP relationship is in direct contradiction to the Marcus theory of electron transfer, and therefore, in some cases, it does not depict the free-energy landscape accurately. Although the BEP relationship states a ). This is a direct consequence of how these two theories modeled the free-energy landscape: while Brønsted, Evans, and Polanyi assumed the free energy to be a linear function of the reaction coordinate, Marcus assumed it to be quadratic. As these assumptions are fundamental to both theories, the applicability of these theories depends on the validity of these presumptions, which should ideally be tested on a system-by-system basis. As a general guideline, the discrepancy between the two is less apparent near equilibrium as the nonlinearity, predicted by the Marcus theory, is less important. However, the predictions from the two theories diverge when the overpotential is increased. The most iconic example is the Marcus inverted region (Hammes-Schiffer, 2009;Parada et al., 2019), where once the overpotential is increased beyond a certain point, the free energy of activation starts to increase. This corresponds to a negative BEP coefficient, which is a direct contradiction to the assumption of 0 < β < 1.
As for the value of the BEP coefficient, β = 0.5 is a frequent assumption, although there is no physical basis for why β should be 0.5 or why it should be independent of the material (Koper, 2011a;Akhade et al., 2018). There is also no basis for why β should be constant throughout the various elementary steps of a multistep electrocatalytic reaction. β = 0.5 is only a convenient assumption for theoretical analysis and should be treated as such.

Overpotential Effects on the BEP Relationship
So far, the entire discussion connecting thermodynamic and kinetics via the BEP relationship is based on the free-energy landscape at equilibrium (zero overpotential). Experimentally, however, electrocatalytic turnover can only be observed if a non-zero overpotential is applied to make the overall reaction exergonic (cf. Section "The Ideal Free-Energy Landscape"). Therefore, there is a discrepancy in the condition between theory (η = 0) and experiments (η > 0) (Exner, 2019a(Exner, , 2020c(Exner, , 2019b. How can the applied overpotential be introduced in the analysis of the free-energy changes (thermodynamics) within the reaction mechanism? For a two-electron process, Exner suggested analyzing the potential dependency of the reaction intermediates in the free-energy landscape, | G RI (η)| (Exner, 2019d(Exner, , 2021b. | G RI (η)| is an overpotential-dependent activity descriptor, reaching its maximum at the overpotential η target that fulfills | G RI (η target )| = 0 eV. As a consequence, thermoneutral bonding of the reaction intermediate at the targeted overpotential, η target , has been identified as a design criteria for active electrocatalysts (Exner, 2020d). This is in contrast to the conventional approach which aims for thermoneutral bonding at zero overpotential (cf. Section "The Ideal Free-Energy Landscape"), and has been coined as the extended Sabatier principle (Exner, 2021b).
Accounting for the overpotential dependency in predicting activity trends appears to be especially important for the evaluation of active catalysts located near the top of the volcano plot (Exner, 2020a). The conventional volcano analysis with the thermodynamic overpotential as the activity descriptor predicts IrO 2 > RuO 2 in both CER and OER (Briquet et al., 2017;Krishnamurthy et al., 2018), although once the overpotential dependence of the free-energy landscape is considered, the experimental activity trend of RuO 2 > IrO 2 could be reproduced (Exner, 2020d). This suggests that the link between the thermodynamics and kinetics via the BEP relationship should not be set at η = 0 V, but rather at the target overpotential (η > 0 V) based on typical reaction conditions in experiments. While the targeted overpotential may be about 50 -100 mV for a twoelectron process such as the HER and CER, a larger value of 300 -400 mV is more realistic for sluggish four-electron processes such as the OER (Exner, 2019c).
A major consequence of the extended Sabatier principle is the shift in the optimum binding energy with increasing driving force (Exner, 2019d(Exner, , 2020c. In Section "The BEP Relationship, " it was demonstrated that the thermodynamically ideal freeenergy landscape for a two-electron process is one which is completely flat at equilibrium ( G RI = 0 eV at η = 0 V). However, this may no longer be ideal if an overpotential is applied. The criterion for the optimum landscape according to the overpotential-dependent description refers to thermoneutral bonding of the reaction intermediate at the target overpotential (Exner, 2020c). Applying the CHE scheme, the free energy of the reaction intermediate at zero overpotential is described by the relation G RI = e|η target | > 0 eV, and hence the optimum binding energy of the reaction intermediate shifts to weak bonding with increasing overpotential (Figure 7), in which G RI of about 100 -200 meV has been recognized as the ideal situation for a two-electron process. It is noteworthy that only the extended Sabatier principle can explain the high activities of Pt ( G RI = 0.21 eV), MoS 2 ( G RI = 0.20 eV), and Mo 2 C ( G RI = 0.13 eV) in the HER (Tang and Jiang, 2016;Handoko et al., 2019;Lindgren et al., 2020) as well as the high activity of RuO 2 ( G RI = 0.13 eV) in the CER (Exner, 2019b(Exner, , 2021b. In contrast, the conventional Sabatier approach would not consider these materials as highly active electrocatalysts, because they are not located at the top ( G RI = 0 eV) of the volcano plot at η = 0 V. So far, the overpotential dependence for the free energy of the reaction intermediate relating to a simple two-electron process has been discussed. In case of a four-electron process, the situation is more complex since three reaction intermediates need to be considered in the analysis (cf. Section "Practical Application and Linear Scaling Relationships"). The common approach to assess electrocatalytic activity for multiple-electron processes relies on the notion of the thermodynamic overpotential, η TD , assuming that the highest free-energy change at the equilibrium potential governs the catalytic performance. It was already recognized by Calle-Vallejo and coworkers that the assessment of a single free-energy change to approximate electrocatalytic activity is too simplistic (Govindarajan et al., 2018;Piqué et al., 2020). These authors put forth an alternate activity descriptor, the electrochemical-step symmetry index (ESSI), which in contrast to η TD analyzes the entire free-energy landscape at zero overpotential. As such, the ESSI is a more balanced activity descriptor than η TD and provides better results in the sorting of highly active electrocatalysts (Exner, 2020f;Piqué et al., 2020). Yet, the ESSI does not contain any overpotential effects in its framework. Recently, Exner derived an overpotentialdependent activity descriptor for multiple-electron processes, denoted G max (η) (Exner, 2020g). In a similar fashion to the descriptor | G RI (η)| for a two-electron process, G max (η) analyzes the thermodynamics of the reaction intermediates at the target overpotential, but at the same time makes use of the freeenergy span model (Kozuch and Shaik, 2011;Kozuch, 2012) by assessing the free-energy difference between the intermediate with smallest free energy and the intermediate with highest free energy in the free-energy landscape. This is achieved by a dedicated procedure of renumbering the electron-proton transfer steps with increasing overpotential, as illustrated in Figure 8.
The concept of G max (η) as activity descriptor goes far beyond the assessment of a single free-energy change, as encountered with η TD . It is noteworthy that the application of G max (η) to OER electrocatalysts predicts that the volcano apex is situated at about G 2 = 1.40 eV (Exner, 2020g). Following the discussion in Section "Practical Application and Linear Scaling Relationships, " the thermodynamic analysis at zero overpotential in terms of η TD purports G 2 = 1.60 eV as the optimum situation in the OER, indicating a shift of the volcano top by about 200 meV. Thus, the displacement of the volcano apex with increasing driving force is in the same order of magnitude as observed for a twoelectron process (Figure 7), whereas the direction of the shift is toward stronger bonding of the oxygen adsorbate (smaller value of G 2 ). The move of the volcano top to stronger oxygen FIGURE 8 | (A) The free-energy span model is applied to determine the distance between the intermediate with smallest free energy and the intermediate with highest free energy in the OER landscape, described by the descriptor G max . By accounting for the overpotential dependency of the elementary reaction steps, indicated by brown arrows in a), the free-energy diagram is translated to an applied overpotential of 300 mV in (B) to evaluate G max (η). Herein, the electron-proton transfer steps are renumbered in that the reaction always commences from the intermediate with smallest free energy. The peculiarity of the overpotential-dependent descriptor G max (η), albeit of thermodynamic nature, is that this activity measure directly scales with the kinetics of the reaction, which is governed by the transition state with highest free energy (G rds # ) in the free-energy diagram. Reproduced with permission from Exner (2020g). Copyright 2020 American Chemical Society.
bonding is further corroborated by a unifying screening approach which considers not only binding energies but also the applied overpotential, rate-determining steps, and catalytic symmetry, and reports G 2 = 1.30 eV as the optimum situation at an OER overpotential of 400 meV (Exner, 2019c). Consideration of the applied overpotential for the approximation of electrocatalytic activity by the thermodynamics is of crucial importance when analyzing free-energy changes to comprehend trends in a class of materials. The concept of the extended Sabatier principle by connecting the thermodynamics to the kinetics at the target overpotential of the reaction is recognized as a valuable tool to distinguish between active and inactive electrode materials, particularly since it has been shown that the approximation of electrocatalytic activity by a single free-energy change, as encountered with the notion of η TD , is insufficient (Exner, 2021c). Yet, this framework cannot entirely close the gap between theory and model experiments. This is also related to the fact that the kinetics is only partly covered in this theory, motivating the following discussion on kinetic effects.

Kinetic Effects Beyond the BEP Relationship
Under the assumption that the reaction intermediates preceding the transition state with highest free energy are in quasiequilibrium (Parsons, 1951(Parsons, , 1958, the rate of an electrocatalytic process is governed by the transition state with the highest free energy in the free-energy landscape (cf. Figure 8). It should be noted that the transition-state free energy consists of a coverage term [(i) → (iii) in Figure 8A] as well as a freeenergy barrier [(iii) → # 3 in Figure 8A], in which the transition state with highest free energy can change upon increasing the overpotential (cf. Figure 8B). However, such a switch in the transition state with highest free energy cannot be predicted solely by thermodynamic considerations, leading to potential discrepancies between thermodynamics and kinetics.
The transition state with highest free energy in the free-energy landscape is connected with the experimentally measurable current density, governing electrocatalytic activity (Parsons, 1951(Parsons, , 1958Exner et al., 2018). Experiments commonly assess the relation between the current and electrode potential by measuring the Tafel plot, in which the applied overpotential is depicted as a function of the logarithm of the current density. This graph facilitates quantifying the so-called Tafel slope (Hu et al., 2004;Guidelli et al., 2014;Shinagawa et al., 2015), which is a measure for the electrocatalyst's sensitivity to the applied electrode potential. The Tafel slope, b, indicates by how many mV the applied overpotential needs to be enhanced to yield an increase of the current density by one order of magnitude: b = ∂η ∂ log 10 j In equation (7), η and j denote the applied overpotential and current density, respectively. Some catalysts are extremely sensitive to the applied overpotential. For instance, a Tafel slope of 40 mV/dec. is obtained for HER on Pt in acid at low overpotentials, whereas the Tafel slope increases to about 120 mV/dec. at higher overpotentials. In alkaline solutions, a single Tafel slope with b = 120 mV/dec. is observed (Sheng et al., 2010). In general, a smaller Tafel slope is desirable because a slight increase in the overpotential can significantly increase the current and thus, the observed electrocatalytic activity. The different Tafel slopes of 40 mV/dec. and 120 mV/dec. can be ascribed to differences in the location of the highest transition state in the free-energy landscape (Parsons, 1951(Parsons, , 1958. The Tafel slope can be expressed in terms of the apparent transfer coefficient, β, which is given by Exner et al. (2018) In case of β = 0.5, the Tafel slope amounts to 120 mV/dec., whereas for β = 1.5 a Tafel slope of 40 mV/dec. is obtained. This finding directly indicates that different transition states govern the rate of the HER over platinum in acidic or alkaline environment at small overpotentials. The apparent transfer coefficient is related to the transfer coefficient for the highest transition state (α) and the number of electrochemical reaction steps consisting of a charge transfer preceding it (n) as follows: Commonly, the assumption α = 0.5 is used, even if deviations from this notion have been reported in the literature (Fang and Liu, 2014;Rojas-Carbonell et al., 2018;Exner and Over, 2019). Now, it becomes clear that n = 0 results in β = 1 /2, and thus b = 120 mV/dec., whereas n = 1 yields ß = 1.5 and b = 40 mV/dec. As a result, when the first or second elementary reaction step governs the rate, the Tafel slope amounts to 120 mV/dec. or 40 mV/dec., respectively. It shall be noted that these different Tafel slopes can be obtained even if the thermodynamic landscape is the same, which is illustrated in Figure 9 for a two-electron process.
Recently, there have been several attempts to bridge the gap between the thermodynamics and kinetics, at least for twoelectron processes such as the HER. For example, a microkinetic study by Koper has explicitly shed light on the relationship between thermodynamic and kinetic free-energy landscapes, and has shown how the thermodynamic landscape may overestimate the catalytic rate due to unfavorable kinetics (Koper, 2013a). A recent study by Ooka and Nakamura demonstrates that the optimum binding energy in terms of maximizing reaction rates may not be thermoneutral as expected from the traditional thermodynamic framework (cf. Section "The Ideal Free-Energy Landscape"). Instead, the optimum binding energy is dependent on the kinetic rate constants, and higher rates are achieved if the binding energy of the reaction intermediate shifts to weak or strong bonding with increasing driving force (Ooka and Nakamura, 2019). Exner has taken the analysis one step further to propose the direction of the binding energy's shift by combining microkinetics with the characteristics of the free-energy diagram in electrocatalysis, pointing out that a shift to strong bonding can be fairly excluded (Exner, 2020e). Rather, displacement of the optimum binding energy to weak bonding is observed for a two-electron process, coinciding with the theory of the extended Sabatier principle (cf. Section "Validity of the BEP Relationship") (Exner, 2019d(Exner, , 2020c. In addition to these and many other studies based on microkinetic considerations in the steady-state approximation (Wang et al., 2009;Holewinski and Linic, 2012;Marshall and Vaisson-Béthune, 2015;Shinagawa et al., 2015;Tao et al., 2019;Mefford et al., 2020;Tichter et al., 2020;, Dauenhauer and coworkers directly calculated the timedependent kinetics of electrochemical processes, proposing that an oscillation in the driving force may lead to an enhancement in the reaction rate compared to the steady state (Ardagh et al., 2019;Shetty et al., 2020). This is because temporal oscillations allow a single catalytic material to exhibit multiple free-energy landscapes, which promote different parts of the reaction.

Chemical Processes Outside the Scope of the Sabatier Principle
Finally, in this section, we will consider factors which influence the observed electrocatalytic performance in ways which are beyond the scope of the Sabatier principle. The Sabatier principle is useful to predict the intrinsic activity of a material. However, in some cases, the intrinsic activity is not the limiting factor of the catalysis.

Mass Transport and Reactant Availability
A prominent example where mass transport and reactant availability significantly influence the observed electrocatalytic reaction is CO 2 reduction Ooka et al., 2017;Lum and Ager, 2018). It has been demonstrated both experimentally and theoretically that the availability of CO 2 , protons, and other electrolyte ions at the electrode heavily influences the selectivity of this reaction (Hori et al., 1989;Murata and Hori, 1991;Varela et al., 2016;Pérez-Gallent et al., 2017;Ringe et al., 2020). While it is difficult to understand the effect of each ion directly from experiments, it is possible to resolve trends of different ions qualitatively using theoretical approaches (Ringe et al., 2019). This requires multi-scale modeling bridging different time and length scales, which accounts for factors such as surface charging, the potential of zero charge, and the double layer capacitance (cf. Section "Challenges to Accurately Predict the Binding Energy").
Another difficulty of CO 2 reduction is the change of the pH during the reaction. For example, the bicarbonate ion inside the electrolyte, either added intentionally or formed as a result from the equilibrium with CO 2 , acts as a pH buffer, thus influencing the pH Gupta et al., 2006). The buffering effect changes as the reaction progresses because there is a gradual accumulation of products inside the reactor. Although gaseous products, such as CO, CH 4 , and C 2 H 4 , have no buffering capacity, oxygenated products such as HCOOH and acetate also have protonation equilibria of their own, and hence are expected to influence the pH. As CO 2 reduction is known to be highly sensitive to the pH (Bumroongsakulsawat and Kelsall, 2014;Schouten et al., 2014;Gao et al., 2015;Birdja et al., 2019), this effect would directly impact electrocatalytic activity and selectivity. It is possible to steer the selectivity of electrochemical reactions by manipulating the local concentrations of substrates, as has been shown experimentally in the case of CO 2 reduction vs. HER (Ooka et al., 2017) or the CER vs. OER (Vos et al., 2018;Wintrich et al., 2019).

Surface Area
The surface area and the availability of active sites on the electrocatalyst surface are other factors which influences the observed activity in ways independent from the intrinsic activity. Although the number of active sites can be estimated theoretically from simplified slab models, there is no direct access to the electrochemically active surface area. From an experimental point of view, the surface area is a critical factor to design new electrodes because the activity is usually evaluated based on the geometric area. For this reason, many reports of active materials are based on porous substrates, such as carbon felt or Pt/Ti mesh, compared with glassy carbon or conductive glass substrates (Park et al., 2012;Xie et al., 2017;Kang et al., 2019).
The surface area may also indirectly influence diffusion at the vicinity of the electrode as a rough surface would disrupt the diffusion layer. Although most studies so far assume the concentration of chemical species within the diffusion layer to be dependent only on the distance from the electrode, a more accurate prediction may be obtained by taking the roughness of the electrode into account (Pajkossy and Nyikos, 1989;Louch and Pritzker, 1993;Kant, 1994). This is especially important for cascade reaction systems, where the intermediate must reach the next catalyst before it diffuses into the bulk solution (Lum and Ager, 2018;Gurudayal et al., 2019). So far, both the electrochemically active surface area and the roughness of the electrode are difficult to describe using DFT, indicating another limitation for the comparison between experiment and theory, which is particularly pronounced when surface blocking or roughness effects govern electrocatalytic activity.

CONCLUSION
The present review has summarized the key concepts of the Sabatier principle, which is widely used today to guide catalyst development. From traditional materials science to the most advanced machine learning studies, there is a general consensus that the binding energy of reaction intermediates is a critical parameter to predict electrocatalytic activity. This descriptor-based approach has also been recently transferred to the surrounding fields of homogeneous catalysis and enzyme catalysis.
Since the concept of the binding energy and the Sabatier principle are central in the current electrocatalysis theory, it is important to highlight not only its usefulness, but also its challenges and shortcomings in a critical light. This may help not only to prevent misleading analyses, but also to clarify what steps need to be taken next. To this end, we have attempted to accentuate at which level issues can arise. Some matters are technical in nature and can be overcome if the accuracy of binding energy calculations is enhanced. On the other hand, some aspects are not an issue of computational accuracy and are more fundamental in nature, such as the bridge between the thermodynamics and kinetics by the Brønsted -Evans -Polanyi relationship or the influence of electrolyte, surface charging and mass transport effects, which require the application of multiscale frameworks. While the simplicity of the CHE model has led to its widespread success and application, more comprehensive approaches which include a grand canonical description of the electrochemical double layer as well as an explicit treatment of the overpotential and its influence on kinetic barriers may be the next step forward. Increasing the model complexity can be promoted considering the emergence of new techniques such as machine learning, which can be used to reduce computational costs.
The theoretical advancement based on the CHE model and the Sabatier principle in its thermodynamic form has been a major driving force behind the advancement of electrocatalysis in the last fifteen years. Building upon and expanding our current framework is critical to promote the identification of next-generation electrocatalysts, and to realize a truly sustainable society.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

FUNDING
HO is grateful to the support from the RIKEN Special Postdoctoral Researcher Program and to the grant-in-aid for scientific research (No. 20K15387) from the Japan Society for the Promotion of Science. JH appreciates the support from the National Natural Science Foundation of China (No. 21802170). KE gratefully acknowledges funding from the Alexander von Humboldt Foundation. KE is associated with the CRC/TRR247: "Heterogeneous Oxidation Catalysis in the Liquid Phase" (Project number 388390466-TRR 247). KE is associated with the RESOLV Cluster of Excellence, funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC 2033 -390677874 -RESOLV. The contribution of KE is based upon the work from COST Action 18234, supported by COST (European Cooperation in Science and Technology).