Rational design/screening of catalysts has long been an important target in heterogeneous catalysis ^{[1, 2, 3]}, since the traditional trial-and-error method is inadequate to meet the demands of the rapid development of the catalyst industry. However, the lack of understanding of the catalytic mechanism at the atomic level greatly limits the efficient design of heterogeneous catalysts. Fortunately, the popularization of density functional theory (DFT) tools in the past two decades and surface reaction kinetics allows us to theoretically evaluate a catalyst’s activity and thus facilitates the overall searching process ^[1]. Among that, the volcano curve, which was introduced by Balandin ^{[4, 5]}, plays a key role in current descriptor-based catalyst screening framework ^[6]. Empirically, a volcano-shaped curve (first ascending and then descending) can usually be obtained upon examining the variation of the catalytic activity for a certain reaction with the position of the transition metals in the periodic table ^{[6, 7]}.

Figure 1 shows a typical volcano curve, which was firstly drawn by Trasatti for the hydrogen evolution reaction (HER) ^[8]. It is clear that the optimal metal catalyst for the HER is located near the peak of the volcano curve. Scientifically, the volcano curve can serve as a comprehensive basis for understanding the catalytic activities of different catalysts or active sites; technically, it is useful guide in high-throughput screening or design of excellent catalysts ^{[9, 10]}.

This concept, however, did not prevail until the last 10 to 20 years, possibly owing to the lack of understanding about the reaction mechanism at the microscopic level, including the key thermodynamic and kinetic reaction information. Fortunately, the use of DFT calculations and quantum chemical softwares with massively parallel computing techniques facilitates the accumulation of relevant information. Firstly, DFT calculations have become a powerful tool to investigate catalytic processes ^{[12, 13]}. The binding energies of intermediates, reaction barriers, surface structures, and so on can be routinely obtained with reasonable accuracy. However, some challenges still remain, such as the ab initio simulation of liquid-solid interfaces and proper consideration of the inherent errors in DFT ^{[14, 15]}. Secondly, with the necessary energies obtained by DFT calculations, kinetic analysis has been increasingly used to assess the performance of catalysts and determine the key rate-limiting factor ^{[7, 16]}. In short, the kinetics serves as a bridge between the macroscopic behavior of catalysts and the microscopic reaction pathways. Scheme 1 illustrates the hierarchical relation among DFT calculations, kinetics, and catalyst design/screening. Nørskov et al. theoretically calculated the rate of ammunition synthesis using the above framework, and the results were within a factor of 3 to 20 of the experimental rate ^[17]. Their result quantitatively demonstrated the feasibility and reliability of catalyst screening using DFT calculations and kinetics. With these developments, the volcano curve has become the main tool in descriptor-based catalyst screening. Typically, the adsorption energy (or other related parameters) should be located near the peak position of the volcano curve. Such an idea has been successfully applied to many catalytic systems, and has shown to be a simple but powerful approach to assist the theoretical design of catalysts.

The Brønsted-Evans-Polanyi (BEP) relation is a convenient tool to qualitatively analyze the volcano curve for certain catalytic reactions ^{[6, 18, 19, 20, 21, 22, 23, 24]}. There is generally a universal linear relationship between the reaction barrier E_a and the reaction enthalpy ΔH. Based on the BEP relation, Norskov et al. made a significant advance in obtaining the volcano curves of a series of catalytic reactions, and proposed the “chemical window” concept by exploring the volcano curve: the optimum catalyst should be the catalyst with a binding energy approximately in the range −2 to −1 eV ^[19]. Considering the BEP relation, Hu et al. proposed a two-step model to explore the volcano curve taking both adsorption and desorption processes into consideration ^{[25, 26, 27]}. As well as the BEP relation, the linear scaling relation between surface species further facilitates the process of catalysts screening ^{[28, 29]}. Interestingly, combining BEP with the scaling relation, 3D volcano curves can be obtained by scaling all of the energies into two surface species ^{[30, 31]}. For instance, Hu et al. described a 3D volcano curve of CO hydrogenation by plotting its turnover frequency (TOF) with respect to the chemisorption energies of C and O, revealing the importance of multiphase catalysts because of some specific scaling relations between surface species ^[30].

Despite the progress in catalyst design by constructing the volcano curve, some basic issues about the volcano curve are yet to be completely understood. One may ask: what are the important factors that determine the volcano curves in heterogeneous catalysis? What is the essential origin of the volcano curve? Historically, the issue could be understood by the Sabatier principle ^{[32, 33]}, which states that an excellent catalyst should have a moderate binding ability. If adsorption is too weak, reactant will not adsorb on the surface, while if adsorption is too strong, intermediate desorption is difficult. The idea of moderate binding revealed by the Sabatier principle is no doubt illumining. However, the physical picture resulting in the volcano curve mainly relys on experience and commonsense, lacking clear explanation and solid evidence. It has also been suggested that the catalyst surface could be blocked by the intermediates if the binding is too strong, which may slow the overall reaction ^{[34, 35]}. Therefore, to obtain a clear picture of the origin of the volcano curve in heterogeneous catalysis, it is necessary to quantitatively demonstrate and reevaluate the volcano curve. Here, we attempt to carry out a self-consistent kinetic analysis to unveil the essential causes of the volcano curve, as well as the related implications in understanding the activity variation mechanism and assisting catalyst screening.

Including both adsorption and desorption processes, the two-step model is believed to capture the essential of catalytic reactions, and can generally serve as a simplified model to describe the activity trend kinetically ^{[25, 26, 27, 36]}. This model is usually written as

where ^* denotes the free site on the catalyst surface, and R, I, and P represent the reactant, intermediate, and product, respectively. Fig. 2 shows a schematic illustration of the corresponding energy profile for this model reaction. E_ad,R (E_ad,P) and E_R^dis (E_P^dis) indicate the reaction enthalpy and energy barrier for the adsorption (desorption) process, which are defined as E_ad,R = E_I - E_R and E_ad,P = E_I- E_P( = E_ad,R- ΔH). Generally, they are correlated according to the linear BEP relation:

In other words, dE^dis/dE_ad= α. The slope α indicates the reaction type with a limited value between 0 and 1, while the intercept term b varies with the local configuration of the catalyst surface. For the activity trend of a specific reaction occurring on a similar surface structure, b is a less interesting term.

With these prerequisites, we performed an analytical analysis following the steady state approximation within the microkinetic framework ^[16], with the aim of revealing the origin of the volcano curve. Based on the De Donder relation ^{[37, 38]}, the reaction rate of each step can be expressed as

Here, θ_*(θ_I) stands for the coverage of surface free sites (intermediates); z_i is the reversibility (z₁ = θ_I/(K_eq1R’θ_*) and z₂= P’ θ_*/(K_eq2θ_I)), which satisfies z_tot = z₁z₂= P’/(K_eq^totR’). R’ and P’ are the relative pressure at the standard state (P_R/P^o, P_P/P^o) of the reactant and product, respectively, while k₁(k₂) and K_eq1(K_eq2) are the rate constant and equilibrium constant of the adsorption (desorption) process, respectively. Specifically, the rate constant k can usually be obtained by transition state theory:

where ∆H^≠,o and ∆S^≠,o are the differences of enthalpy and entropy between the transition state and initial state at standard conditions, respectively, and k_B, h, and T are the Boltzmann constant, Planck constant, and Kelvin temperature, respectively. It should be noted that the calculated activation and chemisorption energies using DFT are actually the total energy change (∆E) at 0 K. To determine of rate constant k at a certain temperature, we apply ∆H^≠,o ≈ ∆E and the small entropy change term (∆S^≠,o) can be ignored for the surface reaction. However, for the adsorption (desorption) process of the gaseous molecule, ∆S^≠,o must be involved. In the two-step model, this gives

where S_R (S_P) represents the standard molar entropy of the reactant (product).

Along with the surface coverage conservation condition (θ_*+ θ_I = 1) and Eq.(3), the expression of coverage (θ_*) and turnover frequency (TOF = r_sur= r₁ = r₂) under the steady state (dθ_I/dt = r₁ - r₂ = 0) can be obtained:

Applying the BEP relation, the rate would inherently become a function solely depending on E_ad,R (i.e. r_sur= f(E_ad,R)) for a specific reaction. Then, we can obtain the analytical expression of r_sur with respect to E_ad,R on the surface:

in which z_tot = z₁z₂, and the reversibility z₁ and z₂ can be obtained for further analysis. Usually, z_i → 1 indicates that step i arrives at equilibrium, while z_i → 0 indicates the step is rate limiting:

To investigate the role of the catalyst surface in the formation of the volcano curve, we also examined a hypothetical gas-phase reaction R(g) → I(g) → P(g), which corresponds to the surface reaction and has the same energy profile as Fig. 2. With similar kinetic derivation, the rate expression in the gas phase (r_gas) can be obtained:

It is worth mentioning that such a gas phase reaction may not exist because the energies of gas-phase intermediates are usually higher than those of the reactants, and there may be no evident BEP correlation between E_a and E^dis in the gas phase. However, such a strategy could provide a reference for uncovering the inherent mechanism that controls the activity trend by comparing r_sur and r_gas.

To understand the activity trend of heterogeneous catalysis using a fully self-consistent method, we partially differentiated lnr_sur with respect to E_ad,R to obtain ∂lnr_sur/∂E_ad,R, which will be zero at the maximum activity, resulting in the formation of the volcano curve. Based on Eq.(8), differentiation of the reaction rate can be solved analytically as Eq.(12).

From Eq.(12), A + B + C + D determines the whole variation behavior of lnr_sur, because the rest of the term is always less than zero. Alternatively, if a volcano curve for lnr_sur exists (mathematically corresponding to an extreme point), A + B + C + D must equal zero at some point (i.e., ∂lnr_sur/∂E_ad,R = 0). We subsequently prove this in a self-consistent manner. Firstly, at the point we have (A + D)/(B + C) = -1, rearranging each term gives the simplified relationship

The right-hand side of Eq.(13) is always less than zero. Therefore, the left-hand side is also less than zero. Solving this inequality, an interesting variation range of E_ad,R emerges:

These boundaries of the upper inequality can be rearranged as μ_P + RTln[(1 - α_R)/α_R] and μ_R + RTln[(1 - α_P)/α_P], which is actually the same as the chemical potential range previously suggested for a good catalyst ^[25]. With these two boundaries, we will prove that a zero point of A + B + C + D must exist using intermediate value theorem. Here, we choose the upper inequality in Eq.(14) for demonstration. From this inequality:

Then, at the left boundary, substituting E_ad,R= ΔH - S_PT - RTln[α_R/P’/(1- α_R)] into A + B + C + D (denoted as f(E_ad,R)) gives:

Based on Eq.(16), it is not difficult to determine that the last part (in the square brackets) is less than zero, i.e., f(ΔH - S_PT - RTln[α_R/P’/(1 - α_R)])< 0. Furthermore, together with Eq.(12), this indicates that:

Similarly, on the right boundary, it corresponds to:

It is a common sense that r_sur should be a continuous function because it describes the real reaction. According to intermediate value theorem, there thus must be a zero point for ∂lnr_sur/∂E_ad,R between ΔH - S_PT - RTln[α_R/P’/(1- α_R)] and - S_RT -RTln[α_P/R’/(1- α_P)]. Moreover, because ∂lnr_sur/∂E_ad,R is greater than zero on the left of the zero point, and less than zero on the right, such an extreme point essentially corresponds to a local maximum of r_sur. In other words, plotting r_sur versus E_ad,R will always lead to a volcano-like curve. It is worth noting that the second situation of Eq.(14) would result in the same conclusion. In short, we have proven the existence of the volcano curve for the heterogeneous catalytic reaction in a completely analytical way.

After proving the existence of the volcano curve, one question naturally arises: how should we determine its origin and what is the key factor determining this phenomenon? To address these issues, we then differentiated r_gas for comparison:

It is easy to notice that ∂lnr_gas/∂E_ad,R is always less than zero, indicating a continuously descending reaction rate rather than volcano-like behavior. It is quite surprising that there is no volcano curve for the gas phase reaction while the only difference between r_sur and r_gas is the incorporation of the catalyst in r_sur. To illustrate and verify this result, we subsequently performed numerical simulations of log(r) versus E_ad,R. From Fig. 3, the right-hand sides of the curves are almost identical but the left-hand sides are completely different. For the gas-phase curve, the reaction rate continuously decreases and no peak exists (blue line). However, in the presence of catalyst, r shows a typical volcano-like shape (red line) with a maximum at about −1.4 eV, which is in agreement with chemical window proposed by Nørskov et al.^[6]

From the above discussion, both the analytical and numerical results indicate that there is a volcano curve for the catalytic reaction but not in the gas phase. If the volcano curve results from a large energy barrier for the intermediate to adsorb or desorb, there is no reason for the volcano curve to disappear in the gas phase. Essentially, it is the decreasing number of free sites on the catalyst surface that results in the left-hand side of the curve being “dragged” down to form a volcano curve (red line). Chemically speaking, the “reaction site” in the gas phase can be deemed to be infinite, thus the volcano curve does not form, while a typical heterogeneous catalytic system usually has a limited number of active sites. This is the reason why on the right-hand side, where the adsorption is weak and there is enough room to accommodate intermediates both on the catalyst surface and in the gas phase, the curves of r_gas and r_sur coincide. For r_gas, the accumulation of intermediates outweighs the augmentation of the desorption barrier, and its value will continuously increase. For r_sur, the magnitude of the free active site rapidly decreases when E_ad,R decreases, which is much more pronounced than the increase of the rate constant, making the curve of r_sur bend down on the left-hand side. Therefore, we propose that the self-poisoning effect accompanying heterogeneous catalysis is the fundamental factor that causes the volcano-shaped activity trend.

Interestingly, we can obtain a clearer physical picture of the generation of the volcano curve by further mathematical analysis. From Eqs.(6),(7) , and (11), there is a simple relationship between r_gas and r_sur:

This indicates that r_sur can be simply expressed by r_gas modified by the free sites, further showing the importance of the free sites in determining the activity trend of heterogeneous catalysis. Taking the natural logarithm of both sides, the differential with respect to E_ad,R can be expressed as

We can determine many interesting implications from this expression. Combining Eqs.(9),(10) and introducing the reversibility z₁ and z₂, the gas-phase term ∂lnr_gas/∂E_ad,R in Eq.(21) can essentially be reformulated as

Similarly, ∂lnθ_*/∂E_ad,R can be expressed as

With these expressions, we can obtain a deeper insight into the slopes of the curves in Fig. 3. On the far left-hand side, adsorption is strong and desorption would be the rate-limiting step, i.e., z₁≈ 1, z₂ ≈ z_tot, and θ_*≈ 0. Then, the slopes of lnr_sur and lnr_gas can be obtained from Eqs.(21)-(23):

Similarly, on the far right-hand side, adsorption is weak and can be regarded as the rate-limiting step, i.e., z₁ ≈ z_tot, z₂ ≈ 1 and θ_*≈ 1:

One can clearly see that the slope of lnr_right is the same for the surface and gas-phase reactions, but that of lnr_left is different. For surface reactions, the slope is negative for lnr_sur,right and positive for lnr_sur,left, resulting in a volcano curve. However, r_gas continuously decreases with a different negative slope on both sides (-α_P/RT and -α_R/RT). These characteristics also agree with the behavior of the curves in Fig. 3. Moreover, Eqs.(22) and (23) can generally be simplified in another way:

in which we assume that α_R = α_P = α for simplicity because α_R and α_P often have a similar magnitude ^[26]. ∂lnr_gas/∂E_ad,R is a constant, which again accounts for the absence of the volcano curve for the gas-phase reaction. However, ∂lnθ_*/∂E_ad,R varies with θ_*, and remains correlated with E_ad,R. Combining Eqs.(21) and (26) gives

Quantitatively, α is a constant between 0 and 1 for a given reaction, whereas (1 - θ_*) ranges from 0 (when E_ad,R is very positive) to 1 (when E_ad,R is very negative). As a result, ∂lnr_sur/∂E_ad,R will be positive on the left-hand side and negative on the right-hand side, leading to a volcano curve. It seems that the competition between ∂lnr_gas/∂E_ad,R and ∂lnθ_*/∂E_ad,R leads to a volcano curve, which further emphasizes the role of the incorporation of the surface term (∂lnθ_*/∂E_ad,R). Interestingly, integrating the equation ∂lnθ_*/∂E_ad,R = (1 - θ_*)/ RT gives the expression of θ_*with respect to E_ad,R: θ_*= 1/(1+ exp((-E_ad,R + C)/RT)), where C is a constant. Thus: ∂lnθ_*/∂E_ad,R = 1/(RT(1 + exp((-E_ad,R + C)/RT))), whose value changes from 0 to 1/RT as E_ad,R changes from being very weak to very strong, indicating that lnθ_*is nonlinear and rapidly decreases as the adsorption strength increases.

A study of Wang et al ^[36]. gave a good example to quantitatively understand the conclusion we made above. Choosing NO oxidation as an example, they performed a careful kinetic analysis of the complete reaction. As shown in Fig. 4, the results numerically indicate that log(r) is determined by the sum of the linear log(k₁₊) and curved 2log(θ_*) terms with some simplification, where k₁₊ is the forward rate constant of O₂ dissociative adsorption. Additionally, they proposed a method to identify the peak position of log(r) assuming the slopes of log(k₁₊) and 2log(θ_*) to be opposite, suggesting the role of the nonlinear character of log(θ_*) in forming the volcano curve, because two straight lines should also give rise to a straight line. Their numerical analysis validates what we analytically determined above. The slope of log(k₁₊) is approximately −α/RT in this circumstance, and the shape of the lnθ_* curve in Fig. 4 (red line) agrees well with Eq.(26). Here, it is again worth pointing out that, because ∂lnθ_*/∂E_ad,R ranges from 0 to 1/RT (when adsorption is very weak or strong), somewhere in the curve the slope must equal α/RT, considering that the value α in the BEP relation is well-known to be in the range 0-1. In short, the nonlinear character of the lnθ_*curve can be considered to be an inherent property of heterogeneous catalysis that contributes to the presence of the volcano curve.

More interestingly, from Eq.(27), when r_sur reaches a maximum, we get θ_*^opt = 1 - α. With this equation, we can estimate the coverage of active sites for the best catalyst. Typically, if α ≈ 0.5, then θ_*^opt = 0.5, which agrees with common sense. Similar concepts have been mentioned elsewhere ^{[26, 39, 40]}, here we would like to emphasize the implications of this equation (θ_*^opt = 1 - α) in understanding the activity trend for rationally screening catalyst with several examples.

Firstly, this equation might be able to explain why the optimal catalyst is different for ammonia synthesis and its reverse decomposition process. Experimentally, the N adsorption strength of the optimal catalyst for ammonia synthesis is always higher than that in its decomposition, which has been previously investigated by Boisen et al. and Kozuch et al. ^{[41, 42]} In our model, this phenomenon can be interpreted as follows. For NH₃ synthesis, the dissociative adsorption of N₂ is the dominant process. As the value of the parameter α ≈ 0.9 for N₂ dissociation ^[27], the optimal catalyst should be with θ_*^opt = 0.1. On the other hand, for NH₃ decomposition, NH₃ dissociation dominates as its concentration increases and α ≈ 0.3 ^[27], corresponding to θ_*^opt = 0.7. Therefore, it is easy to understand that ammonia synthesis requires an early transition metal (e.g., Fe and Ru) with a strong interaction with N₂ to guarantee a small number of free sites ^[43], while more noble metals (e.g., Ru, Co, and Ni) with less strong adsorption ability, which facilitates larger θ_*^opt, are more applicable for NH₃ dissociation ^[44].

Secondly, this equation could provide qualitative insight into the origin of different catalytic reactions with a different optimal catalyst. As a specific case, the optimal catalysts for NH₃ synthesis and CO methanation are usually early transition metals, while those for CO or NO oxidation tend to be noble metals. Differing from the first example, the activation of the key reacting species, i.e., diatomic molecular dissociation of N₂, CO, and O₂, have been proven to obey the same BEP relation on metal surfaces ^[19], and thus would correspond to almost the same θ_*^opt value from the equation θ_*^opt = 1 - α. Consequently, a similar dissociation energy E_ad^opt is essential, as demonstrated in Fig. 5. In this situation, we plotted the general trend of the dissociative adsorption energies of N₂, CO, and O₂ with the metal position in the periodic table (Fig. 5), which clearly shows the location of the best catalyst by the intersection with the horizontal line at a given θ_*sup>opt. The red line corresponding to N₂ and CO dissociation is higher than that of O₂ dissociation (blue line), reflecting a known fact that N₂ and CO are harder to dissociate than O₂. Therefore, from the intersection point (Fig. 5), it is evident that the optimal catalyst for CO or NO oxidation should be more noble than that for NH₃ synthesis or the CO methanation reaction. It is worth noting that with a reasonable approximation in our equation, there should be a range for E_ad^opt (gray shaded region in Fig. 5).

Finally, it is worth noting that we also considered the catalytic model involving two surface species incorporated in one elementary step (i.e., R(g) + 2^* → 2I^* → P(g) + 2^*). In this case, the model can still be solved analytically, but it is too complicated to show. The numerical simulation result gives a very similar conclusion, further verifying our results.

Focusing on the use of the volcano curve in catalyst screening, we briefly discuss some of the latest examples of electrocatalytic conversion (the hydrogen evolution reaction (HER) and dye-sensitized solar cells (DSCs)).

Traditionally, H₂ evolution from water splitting by solar energy over Pt/TiO₂ is deemed to take place on metallic Pt nanoparticles. However, the effect of other Pt species dispersed on the surface of TiO₂ is usually ignored. Toward the high efficient utilization of platinum in the HER of the water splitting process, Wang and Yang et al. elucidated the underlying mechanism of the HER and clarified the effect of the valence state and particle size of the Pt co-catalyst using the DFT-based volcano curve together with experimental characterization ^{[46, 47, 48]}. From DFT modeling and microkinetic analysis, they suggested that the atomic H adsorption energy can serve as an activity descriptor, and indicated that the activity of oxidized Pt species, highly dispersed Pt subnanoclusters ^[43], and even single atoms embedded in the TiO₂ surface ^[44] are located much closer to the volcano peak than metallic Pt nanoparticles (Fig. 6(a)), revealing the key catalytic role of oxidized Pt species in the Pt/TiO₂ hydrogen evolution photocatalyst ^[46]. In particular, motivated from the above understanding, they also proposed a new oxidized PtO co-catalyst to efficiently promote the HER and inhibit the reverse H₂ oxidation reaction ^[48].

DSCs are a way of efficiently harnessing solar energy. The discovery of efficient non-platinum counter electrode (CE) material to catalyze triiodide electroreduction was important for its large-scale application. To avoid trial-and-error tests, a general framework for screening CE material is urgently required. With a two-step kinetic model, Hou et al. ^[14] successfully used the adsorption energy of iodine as a descriptor to evaluate the activity variation. In Fig. 6(b), the catalysts in the blue area are “good catalysts” with relatively high activity. Yang et al. gave this area a physical picture by investigating the peak position of the volcano curve at adsorption and desorption determining circumstance ^[25]. Guided by this search criterion, they predicted and verified the activity of the different crystal planes of Pt ((111) > (411) > (100)) ^[49], facilitating the shape-controlled synthesis of Pt nanocrystals and reducing the amount of Pt used. More significantly, by high-throughput simple adsorption energy calculations of a vast number of candidate materials, they predicted the highly efficient rust (α-Fe₂O₃) and RuO₂ electrodes ^{[14, 50]}. Furthermore, this criterion also guided the rational modification of inert indium oxide (In₂O₃) by interstitially doping heteroatom N to tune its adsorption strength to the iodine atom ^[51].

In summary, a self-consistent mathematical analysis of the origin of the volcano curve was performed using a two-step kinetic model. We analytically proved the existence of the volcano curve in heterogeneous catalysis from a mathematical perspective, revealing the crucial role of the number of free sites. Conceptually, the rapid occupation of the active sites with increasing adsorption strength, which is defined as the self-poisoning effect, results in a slow adsorption rate and is the essential cause of the volcano curve. Some interesting implications for catalyst screening and practical applications of volcano-curve-based models are also discussed and reviewed. The concept of the volcano curve was proposed half of a century ago. Modern kinetic analysis and the development of computational chemistry have allowed us to obtain a deeper understanding of the volcano curve and facilitated its use in real catalytic processes. However, some aspects of the volcano curve are still not completely understood, despite the discussions given above. Some challenging questions still remain to be answered: Is there a simple routine to determine the optimal adsorption energy? How can we describe the kinetics of a multiphase catalyst? These issues deserve further attention in the future.

Prof. Hu thanks the Chinese Government for the “Thousands Talents” program.