Linear scaling relations have led to an understanding of trends in catalytic activity and selectivity of many reactions in heterogeneous and electro-catalysis. However, linear scaling between the chemisorption energies of any two small molecule adsorbates is not guaranteed. A prominent example is the lack of scaling between the chemisorption energies of carbon and oxygen on transition metal surfaces. In this work, we show that this lack of scaling originates from different re-normalized adsorbate valence energies of lower-lying oxygen vs higher-lying carbon. We develop a model for chemisorption of small molecule adsorbates within the *d*-band model by combining a modified form of the Newns–Anderson hybridization energy with an effective orthogonalization term. We develop a general descriptor to *a priori* determine if two adsorbates are likely to scale with each other.

## I. INTRODUCTION

Understanding the formation of a bond between a metallic surface and an adsorbate is critical to the field of heterogeneous catalysis.^{1–10} In principle, the bond (free) energies of all intermediates and transition states of a process fully determines the rate, but, in most cases, the number of such bond energies is so large that it becomes very difficult to single out why a given catalyst surface is the optimal catalyst for the reaction studied.^{11} Scaling relations have provided a method to define a few bond energies as descriptors of the catalytic rate.^{11–14} It has been found quite generally that, for a given surface structure, the bond energies of intermediates and transition states scale linearly with a few bond energies, and once these relations have been found, the full kinetics can be described as a function of the descriptors. These scaling relations result in the “volcano” relationships between rates/selectivities and descriptors that are used extensively in heterogeneous catalysis. It provides a theoretical basis for the heuristic Sabatier principle,^{15} which suggests that the optimal catalyst is one with optimal bond energies for relevant intermediate(s). The scaling relations show which intermediate(s) are relevant and the optimal value(s) for the binding energy of these intermediate(s).

Linear or piece-wise linear scaling relationships are observed when two adsorbates bind to the metal through the same atom, such as C*, CH*, $CH2*$, and $CH3*$,^{16,17} where * denotes the surface. However, it is not always observed when the binding atom is different, such as in the case of C* and O*.^{11} There is currently no *a priori* measure to suggest if the chemisorption energies of two arbitrary adsorbates scale with one another.

In this work, we develop a model that allows for an *a priori* estimation of whether the chemisorption energies of two adsorbates scale with each other. We use the Hammer–Nørskov model, which suggests that the trends in adsorption energies of transition metals are related to the bonding between the adsorbate valence states and metal *d*-bands.^{5} The hybridization energy is described using a modified form of the Newns–Anderson model,^{18} combined with an effective expression for the orthogonalization energy. The model demonstrates that the root cause of lack of scaling between two adsorbed species lies in differences of their re-normalized adsorbate valence energies. For instance, O – *p* states are considerably lower in energy than C – *p* states. Through this analysis, we identify the conditions under which carbon and oxygen *do* scale with each other.

## II. CHEMISORPTION ENERGIES OF CARBON AND OXYGEN DO NOT SCALE

In this section, we discuss the lack of linear scaling relations between the chemisorption energies of C* and O*. We use density functional theory (DFT) calculations to obtain these chemisorption energies (see Sec. VI for computational details). In general, for any adsorbate *A*, the chemisorption energy Δ*E*_{A} is defined as

where A* represents A adsorbed on the surface, $*$ is the free surface, and A is the adsorbate in vacuum. Figure 1(a) shows the chemisorption energies of O* plotted against those of C*. It is evident that there is a significant amount of scatter and no linear scaling relations are observed.

Despite the lack of linear scaling, there is broadly an increasing trend for both C* and O* energies in Fig. 1(a). Metals that bind C* strongly (for example, Fe) also bind O* strongly and those that bind C* weakly (Au, Ag) also bind O* weakly. However, due to the noticeable amount of scatter, this general “rule” does not always hold. Consider the example of platinum: For a constant Δ*E*_{C}, Δ*E*_{O} is about 2 eV weaker than what linear scaling would predict. This weaker than expected binding affects the catalytic activity of Pt, where it is considered weak binding for oxygen based reactions^{19,20} but has a lot of interesting carbon chemistry.^{12,17,21–28} Furthermore, the noble metals (Au, Ag, and Cu) appear to tail-off, having significantly weaker O* and C* chemisorption energies.

Figures 1(b) and 1(c) show the chemisorption energies of O* and C* plotted against the *d*-band center, a commonly used descriptor in heterogeneous catalysis.^{5} The tailing-off behavior of the chemisorption energies at large negative *ϵ*_{d} values is similar to that seen in Fig. 1(a).

Scatter in Fig. 1(a) is markedly reduced when each row of transition metals is plotted separately as shown in Figs. 1(d)–1(f). Mid-row transition metals for each series display linear behavior (Fe, Co, and Ni for 3*d*; Ru, Rh, and Pd for 4*d*; and Os, Ir, and Pt for 5*d*), with the final elements of each row (Cu, Ag, and Au) falling off a linear fit.

These observations raise two fundamental questions for understanding scaling relations:

Why do the noble metals (Cu, Ag, and Au) deviate from the linear scaling relationships constructed along a given row of transition metals [such as in Figs. 1(d)–1(f)]?

What is the source of the scatter in scaling relations across rows of the Periodic Table [such as in Fig. 1(a), deviation in red, blue, and yellow points]?

In the following section, we present a model of chemisorption, which is able to explain these observations.

## III. MODEL OF CHEMISORPTION

In this section, we develop a model for the chemisorption of mono-atomic adsorbates that includes hybridization and orthogonalization contributions. Following the Hammer–Nørskov *d*-band model,^{5} we split the *total* chemisorption energy into the aforementioned components as

The hybridization component of the energy, *E*_{hyb}, is incorporated through the Newns–Anderson model

where Δ_{d} represents a continuous, *ϵ*-dependent coupling element, *V*_{ak}, between the metal *d*-density of states and the adsorbate state. In line with Ref. 18, it is taken to be a semi-ellipse centered around *ϵ*_{d} with a width *w*_{d}. Δ_{0} is a constant contribution applied to approximate the broadening of the adsorbate states into resonances upon interaction with *sp*-states of the metal, Λ is the Hilbert transform^{18} of Δ_{d} + Δ_{0}, and *ϵ*_{a} is the renormalized energy level of the adsorbate, after interaction with the *sp* states of the metal [see the supplementary material (Sec. S2) for further details] and *ϵ*_{f} is the Fermi level (set to 0). The factor of 2 in both terms represents degeneracy of the spin up and down states considered in the model. The adsorbate projected density of states is given by

where the *ϵ* dependence of Δ and Λ has been omitted for clarity. The occupancy of the adsorbate state is given simply by $na=\u222b\u2212\u221e\u03f5f=0\rho aa(\u03f5)d\u03f5$.

We obtain the orthogonalization penalty to the chemisorption energy by comparing the eigenenergies from the narrow-band limit of the Newns–Anderson model with a two-level problem that includes a constant coupling term *S*,^{5,29}

where *f* is the filling of the metal states, given by

Note that spin-polarization can be considered in *E*_{chem} by performing a sum over individual spins in Eq. (2), as in Ref. 18, along with different *E*_{ortho} values for each spin. In this work, we will instead deal with the spin-paired scenario, as our goal is to determine the factors that cause lack of scaling in Fig. 1.

With the above functional form for *E*_{hyb} and *E*_{ortho}, we now parameterize the sum of the two, i.e., *E*_{chem} [through Eq. (1)]. We choose the O – *p* and C – *p* states to have *ϵ*_{a} = −5 and −1 eV, respectively.^{30} The projected densities of states resulting from Eq. (3) with the chosen *ϵ*_{a} for the O* and C* also qualitatively agree with our DFT calculated projected density of states shown in Fig. 2(a) (see Sec. S12), where the localized O* states (shown in orange) are consistently lower-lying in energy as compared to those of C* (shown in purple). Note that the adsorbate projected density of states of both O* and C* show similar features to experimental x-ray photoemission spectra, such as localized states on either side of the *d*-band from Ref. 31 for Cu and Ni, which further validates our approach of using a single renormalized energy for O* and C*.

As in Ref. 5, we assume that the coupling elements are proportional to those from the linear muffin tin orbital (LMTO) framework (denoted as *V*_{sd}),^{32} i.e., $Vak2=\beta Vsd2$, where *β* only depends on the adsorbate species. Analogous to Ref. 33, we write *V*_{sd} as

where *η* are structure factors, which do not depend on bond-length between adsorbate and metal, *M*_{a} depends on only the adsorbate states.^{34}^{,}$Md=s2ld+1\Delta l12$, where *l*_{d} and *l*_{a} are the quantum numbers of the metal *d*-states (*l*_{d} = 2) and adsorbate *p* states (*l*_{a} = 2), respectively, *s* is the Wigner–Seitz radii, and Δ_{l} is the Anderson-width (taken from bulk calculations).^{35} In this work, we set *r* = *s* to incorporate the varying bond length between the adsorbate and metal into the coupling elements (see Sec. S3 for parameters).

We also assume that the overlap matrix element is proportional to the coupling element, *S* = −*αV*_{ak}, where *α* is another constant that is identical for all metals.

The only free parameters in the model for a given adsorbate are *α* and *β*, which are independent of the metallic surface species. We obtain these parameters by fitting *E*_{chem} to the DFT energies in Fig. 1. Table I shows the least-squares errors fitted *α* and *β* for C* and O*. Note that *E*_{chem} is multiplied by 3 to account for three bonds due to the preferred adsorption site being the hollow site (see Sec. S4 for ontop site fitting). Figures 2(b) and 2(c) show parity plots for fitted chemisorption energies from the model vs computed DFT energies. Fitting parameters for CO* based on this model, and a comparison with Ref. 5 is shown in Sec. S5.

Adsorbate . | α (eV^{−1})
. | β (eV^{2})
. | Δ_{0} (eV)
. |
---|---|---|---|

C* (ϵ_{a} = −1 eV) | 0.094 | 1.44 | 0.1 |

O* (ϵ_{a} = −5 eV) | 0.062 | 4.55 | 0.1 |

Adsorbate . | α (eV^{−1})
. | β (eV^{2})
. | Δ_{0} (eV)
. |
---|---|---|---|

C* (ϵ_{a} = −1 eV) | 0.094 | 1.44 | 0.1 |

O* (ϵ_{a} = −5 eV) | 0.062 | 4.55 | 0.1 |

Note that the simultaneous addition of hybridization and orthogonalization components to *E*_{chem} has already been incorporated through the Newns–Anderson–Grimley model of chemisorption^{36} and, more recently, has been applied to study chemisorption and catalytic reactions.^{10,37} Within this framework, *E*_{chem} is obtained from the Anderson Hamiltonian using a non-orthonormal basis, i.e., the obtained energies from the model implicitly include the effects of overlap.

In this work, we instead choose a simpler representation of *E*_{chem} through Eq. (1) and subsequent parameterization of *α* and *β* to describe trends in chemisorption. As mentioned in Ref. 36, the same *E*_{chem} may be obtained with and without explicit consideration of overlap, if *V*_{ak} is changed accordingly. By parameterizing Eq. (1), as done in this section, we effectively alter *V*_{ak} for each adsorbate (through *β*) while accounting for the repulsive interactions through *E*_{ortho}. This parameterization allows us to model trends with a simpler model which functions in the same way as the Newns–Anderson–Grimley model. Equivalency between the two models is demonstrated in Sec. S6, where the parameterization procedure described above is performed with the Newns–Anderson–Grimley model instead of Eq. (1), yielding qualitatively similar behavior in terms of trends in variation of *E*_{chem}, *E*_{hyb}, and *E*_{ortho} with *ϵ*_{d}.

Furthermore, while the orthogonalization model presented in this section is simple and parameterizable, in Secs. S6 and S7, we show that other orthogonalization models yield near identical behavior of *E*_{chem}, with qualitatively similar behavior of variation along *ϵ*_{d} for *E*_{hyb} and *E*_{ortho}.

Having fit *α* and *β* to chemisorption energies from DFT, we compute the variation of *E*_{chem} for O* and C* with *ϵ*_{d} for 3*d*, 4*d*, and 5*d* series of transition metals shown in Figs. 2(d) and 2(e), respectively. Note that, for a fixed row, the chemisorption energy is a function of only *ϵ*_{d}, as $Vak2$ and *w*_{d} for each row of the Periodic Table are functions of *ϵ*_{d}.

The model reproduces two important features present in the DFT calculated energies. First, the variation of *E*_{chem} for O* and C* with *ϵ*_{d} in Figs. 2(d) and 2(e) shows nearly linear behavior up to *ϵ*_{d} ≈ −3 eV. At more negative *ϵ*_{d} values, *E*_{chem} decreases, consistent with the DFT computed trends in Figs. 1(b) and 1(c). Second, the lack of scaling between different rows of transition metals is also seen in the model chemisorption energies. Figure 2(f) shows *E*_{chem} for O* and C*, plotted against each other. Similar to the trends from DFT calculations in Figs. 1(b) and 1(c), taken together the points *do not* scale, but they *do* scale when each transition metal row is plotted separately, up to the weak binding noble metals.

Having developed and parameterized the model, we now use it to explain the variation of the DFT energies with *ϵ*_{d} in Sec. IV.

## IV. ROOT CAUSE OF LACK OF LINEAR SCALING LIES IN DIFFERENCES OF RE-NORMALIZED ENERGIES BETWEEN ADSORBATES

The lack of scaling between *E*_{chem} of O* and C* arises from differences in how *E*_{hyb} and *E*_{ortho} of the two adsorbates vary with respect to *ϵ*_{d}. Previous studies have demonstrated that either hybridization^{5} or orthogonalization^{38} dominate for certain adsorbates and metal surfaces. In this section, we generalize these observations through the model in Sec. III and describe (1) why noble metals Cu, Ag, and Au fall off otherwise linear scaling relationships, (2) why scaling of transition metals is better performed across a fixed row. These two observations come directly from computed energies in Sec. II. We trace the origin of both these effects to differences in *ϵ*_{a} values of the two adsorbates.

### A. Saturation of hybridization energies leads to fall-off of Cu, Ag, and Au from linear scaling

Figures 3(a) and 3(b) show *E*_{hyb} of O* and C* vs *ϵ*_{d}, which has three distinct regions. At positive *ϵ*_{d} values, *E*_{hyb} increases approximately linearly with *ϵ*_{d}. At intermediate *ϵ*_{d} values (−2.5 < *ϵ*_{d} < 0 eV for 4*d* and 5*d* metals and −1.6 < *ϵ*_{d} < 0.9 eV for 3*d* metals), *E*_{hyb} decreases with *ϵ*_{d}. Finally, at large, negative *ϵ*_{d} values (−3.5 eV for 4*d* and 5*d* metals and −1.9 eV for 3*d* metals), *E*_{hyb} → 0, i.e., it saturates completely for O*, while just starting to saturate for C*.

The difference in the variation of *E*_{hyb} with *ϵ*_{d} at large, negative *ϵ*_{d} (−3.5 eV for 4*d* and 5*d* metals and −1.9 eV for 3*d* metals) gives rise to scatter in linear scaling relations between O* and C* along a *fixed* row of transition metals. For O*, the variation of *E*_{chem} in this large, negative *ϵ*_{d} region is determined solely by *E*_{ortho} as *E*_{hyb} → 0. For C*, all regions of *ϵ*_{d} are jointly determined by the contributions of *E*_{hyb} and *E*_{ortho}. This mismatch of energy contributions, *E*_{ortho} alone for O* and *E*_{hyb} + *E*_{ortho} for C*, explains why metals such as Cu, Ag, and Au in the large, negative *ϵ*_{d} region fall off linear scaling relationships.

To quantify the *ϵ*_{d} value at which *E*_{hyb} saturates, we compute the derivative of *E*_{hyb} with respect to *ϵ*_{d}, $Ehyb\u2032$, and determine the *largest* *ϵ*_{d} value at which it tends to zero, which we denote as *ϵ*_{s},

where $Ehyb\u2032$ is given by

where $\Delta d\u2032(\u03f5,\u03f5d)=\u2202\Delta d/\u2202\u03f5d$ and Λ′(*ϵ*, *ϵ*_{d}) = *∂*Λ/*∂ϵ*_{d} (see Sec. S8).

Figures 3(a) and 3(b) (dashed lines) show $Ehyb\u2032$ for both O* and C*, respectively. Consistent with the aforementioned trends in *E*_{hyb}, there are three regions of $Ehyb\u2032$. At positive *ϵ*_{d} values, $Ehyb\u2032$ is negative, in line with linear increasing variation of *E*_{hyb}. At intermediate *ϵ*_{d} values, $Ehyb\u2032$ starts to increase, eventually vanishing at large, negative *ϵ*_{d}. *ϵ*_{s} is indicated by a starred point for each row (in red, blue, and yellow for 3*d*, 4*d*, and 5*d,* respectively). A larger, more positive value of *ϵ*_{s} indicates that saturation of *E*_{hyb} is expected to occur at more positive values of *ϵ*_{d}, as seen in O* as opposed to more negative values for C*.

*E*_{hyb} alone cannot explain the row dependence (3*d*, 4*d*, and 5*d*) of the scaling relations seen in Fig. 1. Figure 3(c) shows the scaling between *E*_{hyb} for O* and C*. For almost the entire range, there is an overlap between the three rows of transition metals (with the exception of the early transition metals), implying that this quantity does not cause the row dependence of scaling seen in the computed energies of Fig. 1(a). To explain this row dependence, we study the other component of *E*_{chem}, i.e., *E*_{ortho}.

### B. Row dependence of scaling stems from orthogonalization energies

Figures 3(d) and 3(e) (solid lines) show *E*_{ortho} as a function of *ϵ*_{d} and dashed lines show the sum of *n*_{a} and *f*, which are multiplicative factors to *E*_{ortho} [Eq. (4)]. At positive *ϵ*_{d} values, *E*_{ortho} increases with decreasing *ϵ*_{d} due to the increase in $Vsd2$, *n*_{a}, and *f*. For *ϵ*_{d} < 0 eV, *E*_{ortho} decreases with decreasing *ϵ*_{d} due to a decrease in $Vsd2$, though *n*_{a} + *f* still tend upward.

Differences in the behavior of *n*_{a} + *f* with *ϵ*_{d} between −3 < *ϵ*_{d} < 0 for O* and C* cause row-dependence of *E*_{ortho} and, hence, *E*_{chem}. Figure 3(f) shows the scaling between *E*_{ortho} of O* and C*; noticeable deviations are seen in linear scaling of *E*_{ortho} but also in the overlap between the rows of transition metals.

Furthermore, differing behavior of *n*_{a} + *f* with *ϵ*_{d} is a consequence of different *ρ*_{aa}(*ϵ*) [determined from Eq. (3)] for different adsorbates, i.e., different *ϵ*_{a}. As an example, Fig. 3(g) shows *ρ*_{aa}(*ϵ*) for C* (purple) and O* (orange) on Rh, a mid-row transition metal. Since *ϵ*_{a} for C* is more positive than O*, the adsorbate projected density of states broadens and splits differently. In the case of C*, *ρ*_{aa} moves toward the upper edge of the *d*-band, resulting in a larger fraction of its peak lying above the Fermi level as compared to O*. This disparate *ρ*_{aa} behavior results in different *n*_{a}, causing differences in *E*_{ortho} [from Eq. (4)] and, hence, row-dependent scaling.

The variations in *E*_{ortho} at *ϵ*_{d} ≈ −2 eV are also caused by changes in *ρ*_{aa}, as the peaks start to cross the Fermi level. We note that this energy variation is less prominent when using the Newns–Anderson–Grimley model, as seen in Sec. S6. However, these variations in energy are of the order of 0.1–0.2 eV and do not change the findings of Fig. 3(f), where scaling of *E*_{ortho} is different by ≈1 eV.

### C. Generalizing findings to determine scaling relations

We now generalize our findings beyond C* vs O* scaling by proposing a descriptor to determine the *ϵ*_{d} value up to which scaling is expected for a *fixed* row, i.e., regions where the chemisorption energies are dictated by both hybridization and repulsive interactions.

In general, if two adsorbates have similar *ϵ*_{s} values, they are likely to scale. Figure 4(a) shows *ϵ*_{s} for a range of re-normalized energy levels *ϵ*_{a}. Since each *ϵ*_{s} value corresponds to a specific *ϵ*_{a} value for each row of transition metals, it suffices to say that two adsorbates with similar *ϵ*_{a} are likely to scale for a given row of transition metal elements. This condition is satisfied for scaling relations between adsorbates binding to the metal with the same atom, such as C* vs CH_{x}* (*x* = 1, 2, 3),^{16} COOH* vs CO*,^{39} or OOH* vs OH*.^{20} It does not hold for adsorbates with very different *ϵ*_{s} and corresponding *ϵ*_{a}, such as C* and O*.

Finally, Fig. 4(b) shows the effect that different *ϵ*_{a} has on the row-dependent scaling behavior. By plotting *E*_{ortho} for O* (for *ϵ*_{a} = −5) against *E*_{ortho} for S* (*ϵ*_{a} = −4 eV), N* (*ϵ*_{a} = −3 eV), and C* (*ϵ*_{a} = −1 eV), we determine the range of *ϵ*_{a} for which row dependent scaling is to be expected. For example, for N* (*ϵ*_{a} = −3 eV), *E*_{ortho} for 3*d*, 4*d*, and 5*d* elements show less mismatch than for the case of *ϵ*_{a} = −1 eV (see Sec. S10 for *ϵ*_{a} for adsorbates), indicating that adsorbates with similar *ϵ*_{a} values are likely to not show row-dependence in their scaling behavior. The range of allowable *ϵ*_{a} values depends on the variation of *E*_{ortho} across and along rows of the Periodic Table. Note that for adsorbates with near-identical *ϵ*_{a} (such as hydrogenated versions of the above monoatomic adsorbates), Fig. 4(b) essentially becomes a parity plot, leading to there being no row-dependence to scaling, as has been observed in a previous work.^{16}

We note in passing that the value of *ϵ*_{s} is dependent on how scaling is constructed. We expect different *ϵ*_{s} and scaling behavior when a scaling line is constructed using just a single metal with varied facets [$Vak2$ is a constant for a fixed metal, rather than =*f*(*ϵ*_{d}), see Sec. S9 for further discussion]. Since scaling lines are, in general, constructed across metals, we emphasize this type of scaling in this work.

### D. Fitting scaling lines in practice

Our analysis of C* and O* scaling, in this section, produces two main insights for the use of scaling lines in constructing activity volcanoes.

When constructing scaling lines for elements of the same row, metals with large, negative

*ϵ*_{d}values, such as noble metals may be left out. The reason is that they are likely to fall off linear scaling lines, as their energies are largely governed by orthogonalization and not hybridization. A more rigorous analysis may include using*ϵ*_{s}to gauge which metals linear scaling should be fit to.Linear scaling is likely to be much improved when all metals belong to the same row of transition metals. The suitability of using multiple rows can be judged based on the

*E*_{ortho}scaling between the two adsorbates.

Note that, in practice, there could be several other reasons for lack of scaling between two adsorbates. Relaxation of the surface and site dependence of the adsorbate might cause scatter on simple transition metal surfaces. For that reason, one should always consider scaling for a fixed surface structure. Adsorbates other than those that are monoatomic could have alternative effects that predominate due to greater degrees of freedom. However, the adsorbates and surfaces we describe in this work are transition metals with the simplest of mono-atomic adsorbates, with chemisorption energies in the dilute coverage limit. The fact that these simple systems do not scale in an obvious manner shows the complexity associated with surface chemical bond formation. Care must be taken while fitting linear scaling relations for transition metals for more complex reactions and intermediates.

## V. CONCLUSION

In this work, we develop a model to understand why chemisorption energies of O* and C* do not follow a linear scaling relation. Following the *d*-band model of chemisorption, our model combines a modified version of the Newns–Anderson hybridization energy with effective terms for the orthogonalization repulsion. We show that lack of linear scaling along a row of the Periodic Table is a consequence of the saturation of the hybridization energy of O* for metals with lower lying *d*-band centers, while the hybridization energy for C* shows little to no saturation at relevant ranges of *ϵ*_{d}. We suggest a measure, *ϵ*_{s}, for the highest *ϵ*_{d} value at which saturation is expected. We quantify the row-dependence of scaling by determining the variation of the orthogonalization energy with the *d*-band center and show that adsorbates with similar *ϵ*_{a} are likely across rows. Finally, we show that linear scaling relations can be recovered in the case of C* and O* by considering metals of the same row of the transition metal series, bar noble metals.

## VI. COMPUTATIONAL METHODS

Density functional theory calculations were carried out using Quantum ESPRESSO.^{40} Core electrons were treated using pseudopotentials from the standard solid-state pseudopotentials (SSSP) Precision pseudopotential database.^{41} Valence electrons were described using plane-waves with kinetic energy and density cutoffs of at least 80 and 600 Ry, respectively. In cases where the recommended cutoffs are greater that these values, the larger of the two is used. Gaussian smearing of width 0.0075 Ry (≈0.1 eV) was used for all calculations. A *k*-point mesh of (4, 4, 1) was used for all surface calculations. The PBE^{42} functional was used to describe exchange and correlation effects. Calculations were considered converged if the maximum force on each atom was less than 10^{−4} Ry bohr^{−1} and the energy difference between successive SCF cycles is less than 2 · 10^{−5} Ry. All calculations were performed *without* spin polarization. Structures were prepared using the atomic simulation environment.^{43} The adsorbate was initially placed at symmetrically inequivalent sites at the start of the relaxation. It is found that adsorbates at the hollow site for all metals have the lowest energy. Provenance of each calculation in this work was stored using AiiDA^{44} and performed using the AiiDA quantumespresso plugin.

The Newns–Anderson expressions were implemented using trigonometric functions from the Python *numpy*^{45} library. Multi-precision C-library *arb*^{46} was used to integrate the projected density of states using the *acb_calc_integrate* function. Decimal precision of 50 was used throughout this work. Least-squares error fitting were performed using the Python *scipy*^{47} library.

## SUPPLEMENTARY MATERIAL

See the supplementary material for the list of symbols, model development, model parameterization details, and extended discussion on parameterization of the Newns–Anderson–Grimley model.

## ACKNOWLEDGMENTS

This research received funding from the European Union’s Horizon 2020 Research and Innovation Program under Grant Agreement No. 851441, SELECTCO2, and from The VILLUM Center for the Science of Sustainable Fuels and Chemicals (9455) from VILLUM FONDEN. Sudarshan Vijay and Karen Chan acknowledge computational resources from PRACE (Project No. 2020235596) and the Juelich Supercomputing Center.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

Python codes used to generate all figures in this manuscript are available at https://github.com/sudarshanv01/newns-anderson. The Newns–Anderson equations and those used in the text for the chemisorption energies and their derivatives are implemented in Python and can be accessed through the following package: https://github.com/sudarshanv01/CatChemi. AiiDA archive files are available at 10.24435/materialscloud:m2-fp.

## REFERENCES

_{2}-CO mixtures catalyzed by platinum: Alkali effects on rates and selectivity

_{2}on Cu and Ni surfaces