Introduction

This contribution to the special issue for José Zagal’s 65th birthday combines two prominent aspects of José’s work: electrocatalysis on molecular complexes, and volcano activity correlations [1]. My aim in this paper is to take the canonical example of hydrogen evolution and introduce a new pathway into the classic models considered decades ago by Parsons and Gerischer [2, 3], in order to derive new volcano correlations between activity and system parameters. The reason for doing so is that classic models for the electrochemical hydrogen evolution reaction (HER) on metal electrodes do not always apply to molecular catalysts. In models for the HER on solid electrodes, the formation of the hydrogen molecule often follows a homolytic pathway, assuming the recombination of two identical hydrogen adatoms (Tafel reaction). In molecular electrocatalysis and enzyme catalysis, hydrogen molecule formation typically follows a heterolytic pathway, combining a (negatively charged) hydride and a proton. Costentin et al. have recently published a diagnostic strategy for distinguishing homolytic from heterolytic pathways in molecular electrocatalysis, showing the importance of reaction rate constants and acid concentration [4]. In this paper, I will map the reactivity for heterolytic hydrogen evolution onto two descriptors: the bond strength of the hydride intermediate bound to the catalyst, and the solution pH (or the pKa of the AH/A− proton donor), as opposed to the Parsons-Gerischer model, which applies to homolytic hydrogen evolution, and which has only the hydrogen-catalyst bond strength as descriptor. The current model is an extension of a recently developed theoretical framework for multiple proton-coupled electron transfer (PCET) highlighting the role of pH in electrocatalytic mechanisms which feature (negatively) charged intermediates [5, 6].

Model and results

Homolytic pathway

The model is based on a well-known model for hydrogen evolution introduced more than 50 years ago by Parsons [2] and Gerischer [3]. The model explains the so-called “volcano activity plot” in electrocatalysis [7–12]. Parsons’ model assumes the usual three reactions for the mechanism of hydrogen evolution:

$$ {\mathrm{H}}^{+} + \mathrm{M} + {\mathrm{e}}^{\hbox{-} }\ \leftrightarrows\ \mathrm{M}-{\mathrm{H}}_{\mathrm{ads}} $$
(1)
$$ {\mathrm{H}}^{+} + \mathrm{M}-{\mathrm{H}}_{\mathrm{ads}} + {\mathrm{e}}^{\hbox{-} }\ \leftrightarrows\ {\mathrm{H}}_2 + \mathrm{M} $$
(2)
$$ 2\ \mathrm{M}-{\mathrm{H}}_{\mathrm{ads}}\ \leftrightarrows\ {\mathrm{H}}_2 + 2\ \mathrm{M} $$
(3)

termed Volmer, Heyrovsky, and Tafel reactions, respectively. M-Hads is the adsorbed hydrogen intermediate. On a metal surface, the adsorbed hydrogen is to be considered as uncharged, as confirmed by DFT calculations [13, 14]. The rates for these reactions are expressed as a function of the electrode potential E SHE (referenced with respect to the standard hydrogen electrode (SHE)), the concentration of protons [H+], the coverage of the adsorbed hydrogen θ H , and the bond strength of M-Hads surface bond, ΔG H (to be understood as a free energy):

$$ {v}_1={k}_1^0\left[{H}^{+}\right]\left(1-{\theta}_H\right) \exp \left(\frac{-\alpha F{E}_{SHE}-\beta \varDelta {G}_H}{RT}\right) $$
(4)
$$ {v}_{-1}={k}_{-1}^0{\theta}_H \exp \left(\frac{\left(1-\alpha \right)F{E}_{SHE}+\left(1-\beta \right)\varDelta {G}_H}{RT}\right) $$
(5)
$$ {v}_2={k}_2^0\left[{H}^{+}\right]{\theta}_H \exp \left(\frac{-\alpha F{E}_{SHE}+\left(1-\beta \right)\varDelta {G}_H}{RT}\right) $$
(6)
$$ {v}_{-2}={k}_{-2}^0{p}_{H_2}\left(1-{\theta}_H\right) \exp \left(\frac{\left(1-\alpha \right){E}_{SHE}-\beta \varDelta {G}_H}{RT}\right) $$
(7)
$$ {v}_3={k}_3^0{\theta}_H^2 \exp \left(\frac{2\left(1-\beta \right)\varDelta {G}_H}{RT}\right) $$
(8)
$$ {v}_{-3}={k}_{-3}^0{\left(1-{\theta}_H\right)}^2{p}_{H_2} \exp \left(\frac{-2\beta \varDelta {G}_H}{RT}\right) $$
(9)

The equilibrium potentials for reactions 1 and 2 have been incorporated into the forward and back rate constants k 01 , k 0− 1 , k 02 and k 0− 2 ; F, R, and T have their usual meaning. The model assumes a Butler-Volmer-type kinetics for the electrochemical steps, with α the corresponding transfer coefficient (that we will simply assume to be a constant equal to 0.5). Variations in the binding energy of hydrogen impact on the rate constant through the Brønsted-Evans-Polanyi (BEP) coefficient β [2, 12, 15].

Expressions for the total current density can be calculated for the case of a Volmer-Heyrovksy VH mechanism (i.e., reactions 1 and 2) and Volmer-Tafel VT mechanism (i.e., reactions 1 and 3). In the steady-state appromixation, assuming that the back reaction to reaction 2 can be neglected, the VH expression for the current density is as follows:

$$ j=-\frac{2F{k}_2^0\left[{H}^{+}\right] \exp \left(\frac{-\alpha F{E}_{SHE}+\left(1-\beta \right)\varDelta {G}_H}{RT}\right)}{1+\frac{k_{-1}^0}{k_1^0\left[{H}^{+}\right]} \exp \left(\frac{F{E}_{SHE}+\varDelta {G}_H}{RT}\right)+\frac{k_2^0}{k_1^0} \exp \left(\frac{\varDelta {G}_H}{RT}\right)} $$
(10)

(The term [H +] multiplying k o2 was erroneously missing from the expression in ref. [12]). In the pre-equilibrium approximation, the following expression can be derived for the VT mechanism:

$$ j=-2F{k}_3^0{\left(\frac{ \exp \left(\frac{\left(1-\beta \right)\varDelta {G}_H}{RT}\right)}{1+\frac{k_{-1}^0}{k_1^0\left[{H}^{+}\right]} \exp \left(\frac{F{E}_{SHE}+\varDelta {G}_H}{RT}\right)}\right)}^2 $$
(11)

At this point, it is useful to introduce the potential versus the reversible hydrogen electrode, E RHE , accounting for the fact that the activity for hydrogen evolution, at least on solid electrocatalysts, is usually measured versus this potential, to take into account the fact that in a real electrochemical device, the second electrode is ideally a fully reversible electrode at which proton and electron transfer always take place concertedly. The RHE is in fact the equilibrium potential of the HE at the pH of the working solution:

$$ {E}_{SHE}={E}_{RHE}-0.059pH={E}_{RHE}-\frac{RT \ln 10}{F}pH={E}_{RHE}+\frac{RT}{F} \ln \left[{H}^{+}\right] $$
(12)

Inserting this into Eqs.10 and 11 gives for the VH mechanism:

$$ j=-\frac{2F{k}_2^0{\left[{H}^{+}\right]}^{1-\alpha } \exp \left(\frac{-\alpha F{E}_{RHE}+\left(1-\beta \right)\varDelta {G}_H}{RT}\right)}{1+\frac{k_{-1}^0}{k_1^0} \exp \left(\frac{F{E}_{RHE}+\varDelta {G}_H}{RT}\right)+\frac{k_2^0}{k_1^0} \exp \left(\frac{\varDelta {G}_H}{RT}\right)} $$
(13)

and for the VT mechanism:

$$ j=-2F{k}_3^0{\left(\frac{ \exp \left(\frac{\left(1-\beta \right)\varDelta {G}_H}{RT}\right)}{1+\frac{k_{-1}^0}{k_1^0} \exp \left(\frac{F{E}_{RHE}+\varDelta {G}_H}{RT}\right)}\right)}^2 $$
(14)

Note that in the expression for the VT mechanism (Eq.14), pH (or [H +]) is now no longer a parameter. In Eq. 13, pH still plays a role because the Heyrovsky reaction, which involves proton transfer, is considered out-of-equilibrium.

Both expressions predict a volcano-type plot for the current as a function of the binding energy of hydrogen ΔG H , as shown in Fig. 1. As shown in a previous paper, the slope of a Tafel plot (i.e., ln|j| vs. E) yields the transfer coefficient α; and the slope of a volcano plot (i.e., ln|j| vs. ΔG H ) yields the Brønsted-Evans-Polanyi coefficient β [12]. It is important to mention two key assumptions made in generating Fig. 1: (1) the rate-determining step is independent on ΔG H , whereas in reality, the rate-determining step may actually differ on either leg of the volcano; (2) ΔG H is a constant independent of coverage.

Fig. 1
figure 1

Activity volcanoes predicted by Eqs. 13 and 14 for the Volmer-Heyrovsky (VH) mechanism and Volmer-Tafel (VT) mechanism. Current in arbitrary units. Most constants in Eqs. 13 and 14 were set equal to 1. The plot for the VH mechanism is drawn for the limit k o2 →0

Full size image

Heterolytic pathway

In molecular electrocatalysis, concerted proton-coupled pathways such as the Volmer and the Heyrovsky reactions (1) and (2) are typically replaced by reactions generating a negatively charged hydride intermediate and its subsequent heterolytic recombination with a proton. In such mechanisms, proton and electron transfer do not take place concertedly in every reaction step. To illustrate this, the following heterolytic pathway for hydrogen evolution should be considered, in which the reaction generating the hydride is written as consisting of two elementary steps, taken into account the molecular character of the catalyst “M”:

$$ {\mathrm{M}}^{\mathrm{N}} + {\mathrm{e}}^{\hbox{--} }\ \rightleftarrows\ {\mathrm{M}}^{\mathrm{N}-1} $$
(15)
$$ {\mathrm{M}}^{\mathrm{N}-1} + {\mathrm{H}}^{+} + {\mathrm{e}}^{\hbox{--} }\ \rightleftarrows\ {\mathrm{M}}^{\mathrm{N}-1}\hbox{--} \mathrm{H}\ \leftrightarrow\ {\mathrm{M}}^{\mathrm{N}}-{\mathrm{H}}^{\hbox{--} } $$
(16)
$$ {\mathrm{H}}^{+} + {\mathrm{M}}^{\mathrm{N}}\hbox{--} {\mathrm{H}}^{\hbox{--}}\rightleftarrows {\mathrm{M}}^{\mathrm{N}} + {\mathrm{H}}_2 $$
(17)

The first step expresses the fact the (metal in the) molecular catalyst changes oxidation state before it will bind the hydrogen in the second reaction. Following Costentin et al. [4], the hydrogen bound to the molecular catalyst can be in (or is a superposition of) two “resonant” states: hydrogen and hydride. In contrast to hydrogen on metal surfaces, the hydride state is more typical for molecular catalysts, and we will consider the hydride state for the remainder of the paper. If an electrocatalyst follows this heterolytic mechanism, this has significant consequences for the pH dependence of the current. To demonstrate this, consider that the catalyst MN is in a film on an inert electrode surface, and there is good and fast electronic contact between the electrode and the catalyst film. The reaction rates are again expressed as a function of the electrode potential E SHE or E RHE , the concentration of protons [H+], the fraction of MN-1 catalyst in the hydride state θ H , and the bond strength of the MN─H─ bond, ΔG H (to be understood as a free energy). We will consider the final step (reaction 17) to be rate-determining and hence neglect its back reaction, and therefore, we must consider that in a solution of high pH, the forward reaction 17 still takes place through a reaction with water, i.e.,:

$$ {\mathrm{H}}^{+} + {\mathrm{M}}^{\mathrm{N}}\hbox{--} {\mathrm{H}}^{\hbox{--}}\to\ {\mathrm{M}}^{\mathrm{N}} + {\mathrm{H}}_2 $$
(18)
$$ {\mathrm{H}}_2\mathrm{O} + {\mathrm{M}}^{\mathrm{N}}\hbox{--} {\mathrm{H}}^{\hbox{--}}\to\ {\mathrm{M}}^{\mathrm{N}} + {\mathrm{H}}_2 + {\mathrm{OH}}^{\hbox{--} } $$
(19)

We obtain an expression for the current by making a steady-state approximation for θ H :

$$ j=-\frac{2F\left({k}_{3a}{10}^{-pH}+{k}_{3b}\right) \exp \left(\frac{\left(1-\beta \right)\varDelta {G}_H}{RT}\right)}{1+\frac{1}{C_{N-1}}\frac{1}{K_2} \exp \left(\frac{F\left({E}_{RHE}-{E}_2^0\right)+\varDelta {G}_H}{RT}\right)+\frac{k_{3a}{10}^{-pH}+{k}_{3b}}{C_{N-1}{k}_2^0{\left({10}^{-pH}\right)}^{1-\alpha }} \exp \left(\frac{\alpha F\left({E}_{RHE}-{E}_2^0\right)+\beta \varDelta {G}_H}{RT}\right)} $$
(20)

where the labels 1, 2, 3a, and 3b correspond to reactions equations 15, 16, 18, and 19, resp., with K 2 = k 02 /k 0− 2 , and we have considered that reaction 15 is in equilibrium such that:

$$ {C}_{N-1}={C}_N^0{\left[1+ \exp \left(\frac{F\left({E}_{SHE}-{E}_1^0\right)}{RT}\right)\right]}^{-1}={C}_N^0{\left[1+{10}^{-pH} \exp \left(\frac{F\left({E}_{RHE}-{E}_1^0\right)}{RT}\right)\right]}^{-1} $$
(21)

where C 0 N is the initial concentration of MN in the film so that C N  = C 0 N  − C N−1 , C N−1 being the concentration of MN−1 in the film.

It can be verified numerically that for certain values of the parameters, Eq. 20 yields a volcano-shaped curve as function of pH, confirming our previous result that PCET reactions in which decoupled proton-electron transfer takes place may exhibit an optimal pH when activity is considered on the “relevant” RHE potential scale [5],[6]. Figure 2 shows a typical volcano curve for the activity of the reaction as a function of pH, for a given fixed value of ΔG H . In the limit that the third term in the denominator is large compared to the second, one derives the following expression for the optimal pH:

$$ {\left[{H}^{+}\right]}_{top}\approx \frac{\alpha }{1-\alpha}\frac{1}{B};\ p{H}_{top}\approx \log B- \log \frac{\alpha }{1-\alpha } $$
(22)
Fig. 2
figure 2

Activity volcano for the heterolytic pathway as predicted by Eq. 20 for the current (in a.u.) vs. pH. k 3a  = 1, k 3a  = 0.001, \( \exp \left(\frac{F\left({E}_{RHE}-{E}_1^0\right)}{RT}\right) \) =105, ΔG H  = 0, K 2  = 1, C N 0  = 1, \( \exp \left(\frac{F\left({E}_{RHE}-{E}_2^0\right)}{RT}\right) \) = 1, \( \frac{1}{k_2^0} \exp \left(\frac{\alpha F\left({E}_{RHE}-{E}_2^0\right)}{RT}\right) \) =103, α = 0.5

Full size image

where B = \( \exp \left(\frac{F\left({E}_{RHE}-{E}_1^0\right)}{RT}\right) \). Since the model now depends on two parameters, i.e., the catalyst property ΔG H and the electrolyte property pH, we plot in Fig. 3 a two-dimensional activity plot to illustrate that the ideal catalyst now optimizes for both properties. It can be seen from this plot that for very negative ΔG H , the activity vs pH dependence changes shape though it in fact still displays a maximum, albeit for unrealistic values of pH. I do not want to dwell on the detailed shape of the surface shown in Fig. 3 as it may be the result of specific and perhaps not very realistic assumptions (under certain conditions) made in the model and in the steady-state approximation to obtain Eq. 20.

Fig. 3
figure 3

Two-dimensional activity plot for the heterolytic pathway as predicted by Eq. 21. Parameters as in Fig. 2, β = 0.5

Full size image

In volcano-type activity correlations in molecular catalysts, as considered in numerous papers by Zagal [1, 11, 16], it is customary to consider the redox potential of the MN/MN−1 redox couple as a descriptor, E 01 , rather than the bond strength ΔG H . It is expected that there is a relation between the redox potential and its ability to stabilize the hydride, i.e., between E 01 and ΔG H . Since it is assumed that the reduced state of the redox couple binds the hydride, a more positive redox potential (i.e., more stable reduced state) should yield a stronger binding of the hydride. One may postulate assume a simple linear equation such as:

$$ {E}_1^0(H)={E}_1^{00}-\gamma F\varDelta {G}_H\ \mathrm{or}\ \varDelta {G}_H=\frac{E_1^0(H)-{E}_1^{00}}{\gamma F} $$
(23)

to describe this relationship.

In the above model, it was tacitly assumed that the HER takes place in an aqueous electrolyte in which the pH is a measurable and adjustable parameter. In non-aqueous electrolytes, it is rather the concentration of the acid that is varied, together with its pKa, i.e., its ability to donate protons. This yields a corresponding equation for [H +] to be used in the above equations. Also, a second proton donor, such as the water in reaction 19, may not exist in non-aqueous solvents.

As a final comment, let me emphasize that the activity plots shown in Fig. 2 and 3 were obtained by evaluating the current as a function of pH at a fixed potential on the RHE scale. In my view, this is a more sensible reference potential, as it automatically corrects for the pH dependence of the equilibrium potential of the overall PCET reaction, and hence, rates are compared at the same thermodynamic driving force. Activity plots using the SHE as the electrode potential reference do not predict volcano relationships, as explained in detail in ref. [6].

Conclusions

In this paper, I have considered a simple model for the electrocatalytic hydrogen evolution reaction which includes a heterolytic H-H bond making step and derived an expression for the activity of the reaction as a function of the ability to stabilize the hydride intermediate and the proton-donating capability of the electrolyte (i.e., pH or pKa of the acid). The model generates a two-dimensional activity volcano, with an optimal hydride binding energy (similar to the existing models for hydrogen evolution) and with an optimal pH so that the acidity of the electrolyte solution (or equivalently the proton-donating capability of the acid) matches the proton-accepting capability of the catalyst active intermediate. The model highlights that the acid/base character (or generally, the charged nature) of the catalytic intermediates featuring in so many pathways relevant to molecular catalysts is what makes these pathways so sensitive to electrolyte pH, and what often distinguishes them from pathways on heterogeneous catalysts [5].