Modelling adsorbate coverage on complex alloy surfaces

Abstract

In order to extend catalysis theory to complex alloys and multiple adsorbates, we have to face the fact that the number of possible surface site–adsorbate combinations gets too large to calculate. We, instead, define rules for adsorbate–adsorbate interactions; specifically, blocking rules in terms of disallowed local adsorbate–adsorbate configurations. We then conduct simple simulations to investigate how different rules entail certain outcomes. For the PdAg intermetallic and PdAg solid solutions, we find that the presence of Ag atoms hinders O* species from covering the whole (111) surface, which is the case for unary Pd(111), and instead allows for adsorbed *OH species. We predict that the adsorbed *OH species improves the oxygen reduction reaction activity because they have adsorption energies at the top of the activity volcano. Experiments can use our results to distinguish between the different possible PdAg(111) alloy surface manifestations, and to better understand adsorbate coverage on complex alloys. Lastly, we use our approach on Ag₁₄Ir₁₆Pd₃₀Pt₁₄Ru₂₆ high-entropy alloys, but find that the choice of adsorbate–adsorbate interaction rules affects the oxygen reduction in less distinguishable ways compared to the binary PdAg alloys.

---

1 Introduction

Adsorption on surfaces is relevant in many contexts and central to the theory of heterogeneous catalysis. Langmuir started the surface science approach to adsorption by using a set of basic “rules”, i.e., assumptions.¹ The rules are: equilibrium exists between the adsorbed and desorbed state of the molecule; all the adsorption sites on the surface are the same; adsorption is limited to one monolayer (ML); and there is no adsorbate–adsorbate interaction. These rules directly give rise to the Langmuir isotherm for adsorption – an analytical expression that describes surface coverage (Fig. 1a).

Fig. 1 (a) Coverage, θ, as function of potential, U, (or chemical potential). The Langmuir isotherm model in blue, , where ΔG_ads is the adsorption free energy. The Frumkin isotherm model in red, , where α is the adsorbate–adsorbate interaction term. The plot uses α = 5k_BT. (b) Examples of local *OH–O*, *OH–*OH, and O*–O* adsorption configurations that could be disallowed by the blocking rules.

The last rule regarding the non-interacting adsorbates is known to be violated on real catalyst surfaces. Actually, the Langmuir rules imply that there is an infinite repulsion between molecules on the same site and no interaction at all between molecules on neighbouring sites. The simplest way of including adsorbate–adsorbate interactions is the Frumkin adsorption model for surface coverage wherein the interaction between adsorbed molecules is approximated by a mean-field linear term, which only depends on the average coverage and not the local adsorbate configuration on the surface (Fig. 1a). A related version of the Frumkin isotherm, the Temkin isotherm, could be deployed for surfaces with heterogeneous adsorption energies. However, the result is similar to the Frumkin isotherm, as the binding energy is assumed to increase linearly, thus implicitly mimicking the effect of the mean-field interaction term. Fig. 1a shows how the inclusion of the interaction term results in a delayed coverage response to the applied potential, which weakens the adsorption energy and correspondingly impacts the catalytic activity and selectivity. Additionally, the inclusion of an interaction term means that the coverage of adsorbed molecules must be solved iteratively; thus, the model loses the simplicity of the Langmuir isotherm, even in the simplest approximation of adsorbate–adsorbate interactions.

In electrocatalysis on uniform surfaces, isotherms resembling those of the Langmuir and Frumkin models are observed, e.g., for underpotential deposition of hydrogen on Pt(100) and Pt(111), respectively.² The Langmuir isotherm also describes *OH adsorption on Pt(111), and the Frumkin isotherm describes intermediates related to oxygen evolution on Ir oxide.^3,4

All the local adsorbate–adsorbate interactions, i.e., nearest neighbour, second nearest neighbour and so on, can be calculated and implemented in Monte Carlo simulations. It is complicated and far from trivial to conduct these simulations, but this will certainly see increased use in the future, especially for studying dynamics. Such simulations can provide accurate answers to specific questions, but it is difficult to get chemical intuition from the results. Another commonly employed technique to study the effect of coverage is the recording of cyclic voltammograms (CVs). As CVs are both achievable experimentally and theoretically, this methodology presents a helpful way of comparing results across different approaches.⁵ When the number of different adsorbates increases, it becomes even more complicated, and if the surface is heterogeneous, as for high-entropy alloys (HEAs), it becomes practically impossible. Thus, understanding the adsorbate–adsorbate interaction on HEAs represents a frontier in the atomic-scale modelling of catalysis.

Our approach is pragmatic; we cannot exhaust the different possible interactions on a HEA surface with density functional theory (DFT) simulations. This limitation arises from the inherent complexity of HEAs, where a single on-top site exhibits 5¹⁰ distinct permutations when accounting for the nearest neighbouring atoms. Therefore, we define a set of rules for adsorption, just as Langmuir did. We then conduct simple simulations to investigate the consequences and predictions of these rules and compare with experiments (if such are available). In most cases, there are no experiments to compare with, but we can still investigate how different rules entail certain outcomes.

For the simple case with only one type of site and one type of intermediate, we can use the Frumkin or Langmuir isotherms directly. However, it is known that spectators can block active sites and reduce the number of sites that participate in catalysis.^6–8 So, if spectators and other intermediates are adsorbing, we have to model how adsorbates block other adsorbates from accessing nearby surface sites. As it is not feasible to use DFT directly on a large surface or calculate all possibilities, we start by estimating the adsorption energies on the different types of sites for the different possible (non-interacting) intermediates. We define this as the gross distribution for non-interacting adsorption of each intermediate. This distribution suggests that each surface atom could bind to more than one intermediate, which is deemed unphysical. At this point, we must postulate rules for how to fill the surface with adsorbates.

We assume that there can be only one intermediate per site, thus the intermediates with the strongest adsorption will block sites for weaker adsorbing intermediates. More specifically, we define blocking rules in terms of disallowed local adsorbate configurations. Examples of such rules could be to disallow adsorbed *OH and O* that share a surface atom, three adsorbed *OH on adjacent surface atoms, or two adsorbed O* that share a surface atom (Fig. 1b). Furthermore, imagine a surface with different types of sites; we can fill the surface from the strongest binding sites, which might then block associated sites for weaker adsorbing intermediates. Even if this seems reasonable, these are rules that we are imposing, and other rules could be imagined.⁹ The resulting distribution after filling the surface according to the imposed rules we name the net distribution. The net distribution determines what sites are accessible for catalysis.

The catalytic activity of a surface is the activity of each site summed over all the sites. Catalysis is a complex phenomenon, and catalytic performance cannot be reduced to a simple easily calculated quantity. Consequently, we make use of descriptors known to correlate with catalytic performance. Adsorption energies have proven to be the most predictive descriptor in catalysis.¹⁰ Adsorption energies fit well with the Sabatier Principle, which states that the optimal catalyst strikes a balance between binding reactants strongly enough to facilitate their transformation and releasing products efficiently. In recent decades, the advent of DFT simulations has enabled the development of a quantitative theory for catalytic activity based on these principles.¹¹ The resulting catalytic activity is often represented as a volcano curve, which establishes a relationship between a descriptor (the adsorption energy, ε, of a reaction intermediate) and catalytic activity. Given a specific value for the descriptor, the volcano curve, A(ε), predicts the corresponding activity.

The activity in the descriptor approach is:

Activity ∝ ∫D(ε)A(ε)dε.

D is the net distribution of adsorption energies, ε is the adsorption energy descriptor, and A is the activity as a function of the adsorption energies or the volcano curve. On a uniform surface all sites have the same adsorption energy, ε₀, and D(ε₀) is then a delta function and the activity becomes A(ε₀) times the number of sites. On surfaces with sites of different adsorption energies, D(ε) can take any shape; however, the area under the distribution corresponds to the total number of occupied sites. Thus, the integral over adsorption energy captures the catalytic contributions from all occupied sites.

In this paper, we investigate a special case where adsorbate–adsorbate interactions and alloying can unlock highly active sites for catalyzing the oxygen reduction reaction that are otherwise inaccessible on the pure element surface. We look at Pd and Ag as unary metals, an intermetallic, and a solid solution of the two. We have limited the analysis to alloys composed of Pd and Ag, as these are reported to outperform both of their associated pure metal alloys in terms of activity, of which Pd performs decently, while Ag performs poorly.¹² Additionally, the chosen alloy polymorphs display unique adlayer patterns due to their structural order (unary and intermetallics) or structural disorder (solid solution), making them ideal surfaces for our analysis. While neither the intermetallic nor the solid solution qualifies as a HEA, from a simulation standpoint, they exhibit a phenomenon also observed on HEA surfaces, where we initially identified it.¹³

In the last part, we extend our analysis to the (111) surface of an Ag₁₄Ir₁₆Pd₃₀Pt₁₄Ru₂₆ HEA. Opposite to the PdAg intermetallic and PdAg solid solution, the simulated potential-dependent ORR currents from the HEA are rather featureless, and applying different blocking rules only scales the resulting curves. This actually makes PdAg alloys a better experimental choice for studying how adsorbate coverage on complex alloy surfaces influences the (electro)catalytic activity.

2 Computational details

Throughout this work we have utilized a graph convolutional neural network as a regression model. The regression model is developed in line with the methodology fully elucidated in earlier works.¹⁴ Briefly, the regression model is based on DFT-calculated *OH and O* adsorption energies at on-top sites and face-centered cubic (FCC) hollow sites, respectively, on FCC-structured HEA (111) surfaces comprising Ag, Ir, Pd, Pt, and Ru. It should be noted that all DFT calculations are performed on adsorbates in isolation, which results in the regression model predicting exclusively non-interacting adsorption energies. The regression model is employed to predict *OH and O* adsorption energies on simulated 96 × 96 × 3 atom-sized slabs, either as unary alloys, intermetallics, or solid solutions. These extended slabs form the foundation for the presented analyses.

The DFT calculations, used for training, testing, and validation,^15,16 were performed with the revised Perdew–Burke–Ernzerhof (RPBE) exchange–correlation functional as implemented in the GPAW code version 1.5.1.^17–19 The wavefunctions were expanded in plane-wave basis sets with an energy cutoff of 400 eV and the Brillouin zones were sampled with 4 × 4 × 1 k-point Monkhorst–Pack grids.^20,21 The structures were set up and manipulated using the atomic simulations environment version 3.17.0 (ASE).^22,23 The Ag–Ir–Pd–Pt–Ru alloys were simulated as 3 × 3 × 5 atom-sized slabs with periodic boundary conditions and 20 Å of vacuum in the z-direction. The applied FCC lattice parameter of the slab was taken as the weighted mean of the top-layer constituent-elements-calculated FCC lattice constants in line with Vegard’s Law.^24,25 The two bottommost layers of the slabs were held fixed to emulate bulk properties and only the remaining layers were allowed to relax until the maximum force on any atom was within a threshold of 0.1 eV Å⁻¹.

The *OH and O* adsorption free energies on the multimetallic alloys were obtained by reference to Pt(111) using the following assumptions:

Assumption (1): the *OH adsorption free energy on Pt(111) is ΔG^Pt_*OH = 0.76 eV. This is an observation from blank Pt(111) cyclic voltammograms (CV) in water.²⁶ This adsorption energy also fits the observation that Pt(111) adsorbs *OH 0.1 eV stronger than the optimum of the ORR activity volcano.²⁷ The optimum of the ORR volcano is fixed by the *OH/*OOH scaling relation, ΔG^opt_*OOH = ΔG^opt_*OH + 3.2 eV, and the formation free energy of two water molecules of −4.92 eV to ΔG^opt_*OH = (4.92 − 3.2) eV/2 = 0.86 eV described in detail elsewhere.^28,29

Assumption (2): on Pt(111), the O* and *OH adsorption free energies display a perfect scaling relation, i.e., ΔG^Pt_O* = 2ΔG^Pt_*OH. This can be inferred from the same CVs in water and from simulations of the free energy contributions of *OH and O* on Pt(111).³⁰

Assumption (3): entropic contributions to the DFT energies to obtain free energies are the same as on Pt(111), i.e., the differences between free energies and DFT energies are assumed to be identical:

ΔG^alloy_*OH − ΔG^Pt_*OH = ΔE^alloy_*OH − ΔE^Pt_*OH and ΔG^alloy_O* − ΔG^Pt_O* = ΔE^alloy_O* − ΔE^Pt_O*.

The above assumptions allow the computation of *OH and O* adsorption free energies as

ΔG^alloy_*OH = ΔE^alloy_*OH − ΔE^Pt_*OH + ΔG^Pt_*OH = (E^alloy_*OH − E^alloy) − (E^Pt_*OH − E^Pt) + 0.76 eV,

G^alloy_O* = ΔE^alloy_O* − ΔE^Pt_O* + ΔG^Pt_O* = (E^alloy_O* − E^alloy) − (E^Pt_O* − E^Pt) + 0.76 eV.

As the potential, U, is increased, the surface is filled with *OH and O*. The free energies for oxidative adsorption of OH and O are ΔG^alloy_*OH(U) = ΔG^alloy_*OH − eU and ΔG^alloy_O*(U) = ΔG^alloy_O* − 2eU. Hence, the competition between *OH and O* for sites, at any potential, can be determined by comparing ΔG^alloy_*OH alloy with ½ΔG^alloy_O*. In an anodic potential sweep, adsorbates fill the surface sites from lowest to highest adsorption energies. For convenience in the following, we define the following *OH and O* adsorption free energies relative to *OH on Pt(111):

ΔΔG^alloy_*OH ≡ ΔG^alloy_*OH − ΔG^Pt_*OH,

ΔΔG^alloy_O* ≡ ½ΔG^alloy_O* − ΔG^Pt_*OH = ½ΔG^alloy_O* − ½ΔG^Pt_O*.

The *OH and O* binding energies were obtained on all unique sites of ∼300 unique Ag–Ir–Pd–Pt–Ru slabs totalling ∼5000 unique adsorption energies, which constituted the training, validation, and test data for the regression model in an 80-10-10 split, respectively. The data was converted to graph-type features using the geometry tools incorporated in ASE to define the adjacency of atoms within a cutoff-radius determined by the constituent elements’ covalent radii. The regression model itself was constructed within the pyTorch framework with an architecture of 3 convolutional layers of 18 dimensions and ReLU as the activation function.³¹ Each training epoch used the AdamW optimizer (betas = (0.9, 0.999), eps = 10⁻⁸, weight_decay = 0.01), the mean squared error (MSE) as loss function, a learning rate of 1 × 10⁻³, and a batch size of 64.³² The rolling mean (±10 epochs) of the validation error was compared with that of the earlier epochs excluding the last 100, and if the rolling mean validation error had not decreased by 1% compared to the prior errors, early stopping was invoked. The model weights that gave rise to the lowest validation error were then chosen for the trained model.

3 Results and discussion

Our goal is to predict the ORR activity of different PdAg(111) alloy surface manifestations. We estimate ORR activity by using the stability of adsorbed *OH/O* species relative to *OH adsorption free energy on Pt(111) and assert that the highest ORR activity is achieved when the *OH/O* adsorption energy fulfills ΔΔG^alloy_*OH/O* = 0.1 eV. This follows the approach described in earlier works.¹³ The main obstacle is that the net distribution of adsorption energies, rather than the less complex gross distribution, determines the ORR activity.

3.1 Gross *OH/O* adsorption energy distributions on different PdAg(111) alloy surfaces

The non-interacting gross *OH/O* adsorption energy distributions on different manifestations of PdAg alloy surfaces are depicted in Fig. 2.³³ The first considered possibility is that the surface is unary Pd(111) or Ag(111), for instance, brought about by surface segregation or element-specific dissolution.^34–38 The gross *OH/O* adsorption energy distributions are distinguished by sharp peaks, where the *OH adsorption energy on Pd(111) is very close to the optimum value (ΔΔG^Pd_*OH = 0.09 eV). However, as will be discussed later, these sites are unavailable due to the strong O* adsorption (ΔΔG^Pd_O* = −0.02 eV) which fully covers the Pd(111) surface. On Ag(111), the coinciding *OH and O* adsorption energies are weak compared to adsorption on Pd(111) and do not contribute to the overall ORR activity in an appreciable amount until a large overpotential is reached.

Fig. 2 Gross *OH/O* adsorption energy distributions and illustration of the atomic configurations at the surface for (a) unary Pd(111) or Ag(111), (b) intermetallic PdAg(111), and (c) equimolar solid solution PdAg(111). In (a) and (b) it is indicated which O binding energy corresponds to which metal.

The second possibility is that the Ag and Pd atoms sit in an ordered pattern at the surface, referred to as an intermetallic. Here, we study the (111) surface cut of a layered bulk PdAg alloy that exposes alternating rows of Pd and Ag atoms (Fig. 2b). The bulk layered PdAg alloy is predicted to be stable by Materials Project; however, experiments suggest that PdAg does not form ordered phases at 25 °C and above.^39,40 Nevertheless, the intermetallic is included, as it showcases the ideal case of adlayer formation and the ligand effect coming together. The intermetallic surface has two distinct on-top adsorption sites on Pd and Ag atoms, respectively, and two distinct FCC hollow adsorption sites consisting of two Pd atoms and one Ag atom, or one Pd atom and two Ag atoms. The adsorption energy distribution still consists of sharp peaks; however, the peaks can be shifted from the unary metal adsorption energies. For instance, the most stable O* adsorption energy is destabilized slightly compared to O* on unary Pd(111), because the O* atom interacts with two reactive Pd atoms and one inert Ag atom instead of three Pd atoms. In contrast, *OH adsorbed on-top the Pd atom has an unaltered adsorption energy compared to unary Pd(111) and is now slightly more stable than O*.

The third possibility is that the Pd and Ag elements are randomly distributed on the surface, resulting in what is referred to as a solid solution (Fig. 2c). This is the situation predicted by most phase diagrams for bulk Pd_xAg_(1−x) alloys.^41–43 The adsorption energy distributions are broadened around the unary metal adsorption energies due to ligand effects, while the distributions can be shifted in the cases where the element composition of the surface differs from that of the bulk due to strain effects (which is not included in our model).^44–48 We note that the PdAg surfaces depicted in Fig. 2 are idealized cases and that the surface could exist in a state in between those illustrated in Fig. 2. For instance, the surface could be enriched in one element without the other being fully absent, or there could be short-range local order in an otherwise disordered surface.^49,50

3.2 Blocking rules for O* and *OH

After having obtained adsorption energies for every possible site on the different PdAg alloy surface manifestations, we need to consider to what extent adsorption of one species blocks adsorption of other species at the neighbouring surface sites. Ideally, we would recalculate the adsorption energies at the neighbouring surface sites after every adsorption event, but the regression model does not account for repulsion between species, while DFT calculations quickly become untenable for alloys with many different surface-site combinations.

Instead, we define two sets of adsorbate-specific blocking rules for O* and *OH, and subsequently examine how the different sets of rules affect the final ORR activity estimate. Specifically, we use an iterative procedure: we begin by placing the most stable species on the surface, then sequentially add the next most stable species that are not excluded by the blocking rules imposed by the previously placed species.

The first set of rules, originally defined in ref. 13, is referred to as mutual blocking. Mutual blocking stipulates that O* can only adsorb in FCC hollow sites, and when it does, it blocks on-top *OH adsorption on the three surface atoms that make up the hollow site. Similarly, *OH adsorbing at on-top positions blocks the three nearest FCC hollow sites. Finally, an adsorbate blocks further adsorption on the same binding site. Fig. 3 explains how the rules are used to place adsorbates on a PdAg solid solution surface.

Fig. 3 Illustration of how to fill the solid solution surface shown in (a) using blocking rules. (b) The most stable species (an O*) is placed first and blocks three on-top *OH sites in mutual blocking (marked by crosses) and an additional six FCC hollow sites in extended blocking (marked by diamonds). (c) The next most stable species (an *OH that has not been blocked) in turn blocks three FCC hollow sites (marked by crosses) in both set of rules. (d) The next most stable species (another *OH) blocks additional FCC hollow sites in both set of rules, while the combination of two adjacent *OH species blocks the two shared on-top neighbour sites in extended blocking. This process continues until there are no unblocked sites left, resulting in the coverage (e) for mutual blocking and (f) for extended blocking. The colour scheme overlaid in (b)–(d) depicts the GNN-predicted adsorption energy of each site, with the dark red regions indicating strong binding sites (more negative adsorption energy) and the blue regions indicating weak binding sites (more positive adsorption energy).

The surface atomic configuration is depicted in Fig. 3a, and in the following figures we overlay the adsorption energies from the regression model as a heat map. The most stable species is an O*, which is placed in Fig. 3b, and blocks three on-top *OH sites (marked by crosses). Next, the most stable species, that is not blocked, is an *OH (Fig. 3c), which, in turn, blocks three FCC hollow sites (marked by crosses). This iterative process can be continued, where the next adsorbate is another *OH (Fig. 3d), until all sites are blocked (Fig. 3e). Importantly, this set of rules does not leave uncoordinated surface atoms and allows for 1 ML coverage of O* or *OH.

The second set of rules, denoted extended blocking, retains the rules from mutual blocking, but adds that an O* species blocks the nearest six FCC sites (diamonds in Fig. 3b). This limits the O* coverage to 1/3 ML, which is often predicted to be the highest possible O* coverage at conditions relevant to ORR.^51,52 It also avoids unphysical coverages, where bulk-metal oxidation becomes more energetically favourable than surface O* adsorption. Extended blocking further adds that two adjacent *OH species block the two shared on-top neighbour sites (diamonds in Fig. 3d), which we do to promote hexagonal patterns of *OH species with up to 2/3 ML coverage. Experimental cyclic voltammetry suggests that Pt(111) forms a surface adlayer with between 0.35 and 0.45 ML coverage of *OH.^53–56 This adlayer is often envisioned as a hexagonal pattern of alternating *OH species (1/3 ML) and H₂O species (1/3 ML) connected by hydrogen bonds, causing a very stable adsorbate configuration.^57–59

Together, the rules also ensure that the maximum coverage of O* (1/3 ML) and of *OH (2/3 ML) yields an equally oxidized surface and therefore that their relative stability is potential independent. This improves on mutual blocking, where we might predict one monolayer of *OH at high potential, even though we know that the more oxidized one monolayer of O* becomes more stable at high potentials. Finally, the O*–O* repulsion and the *OH–*OH interaction are small or unchanged at or below these coverages.^60–62 The weak interactions also mean that it is probably not necessary to include a Frumkin term. The full surface coverage after applying the extended blocking rules (Fig. 3f) has some important differences to that from the mutual blocking rules (Fig. 3e), namely that the surface contains empty sites that do not support adsorbates and therefore do not provide ORR activity. There is also significantly more intermixing of *OH and O* and fewer domains dominated by one species. Finally, hexagonal *OH patterns are visible in the bottom-left part of the surface, whereas the bottom-right part of the surface shows that *OH also sits in line patterns.

Before discussing the net *OH/O* adsorption energy distributions for the different PdAg(111) alloy surface manifestations, we want to emphasize, in detail, how the blocking rules affect catalytic activity. Applying the extended blocking rules to unary Pd(111) results in an ordered O* adsorption pattern with 1/3 ML coverage (Fig. 4, top row). The O* species have low ORR activity because ΔΔG^Pd_O* is −0.02 eV rather than the optimal +0.1 eV. In contrast, if the Pd(111) surface contains a single Ag atom, the ordered O* pattern is interrupted and the extended blocking rules add two *OH species instead of an O* species (Fig. 4, bottom row). The ΔΔG^Pd_*OH is 0.09 eV, so the two *OH species have one order of magnitude higher ORR activity than the combined O* species.

Fig. 4 Coverage of adsorbed species using extended blocking rules for Pd(111) (top row) and Pd(111) with a single Ag atom (bottom row). The presence of *OH species in Pd(111) with a single Ag atom improves the ORR activity by one order of magnitude.

3.3 Net *OH/O* adsorption energy distributions on different PdAg(111) alloy surfaces

We now apply the blocking rules in order to figure out which species are absorbed on the surface with which adsorption energies, i.e., the net *OH/O* adsorption energy distributions. We start with the mutual blocking rules for unary Pd(111) and Ag(111) (Fig. 5a). The O* species are most stable in both cases, and, due to the absence of O*–O* blocking, they can cover the whole surface. The net distribution therefore solely consists of O* species and the coverage, as a function of reactant chemical potential, is obtained by applying the Langmuir isotherm or the Frumkin adsorption model. The difference between net and gross distribution is therefore the absence of *OH adsorption energies in the net distribution. Applying the extending blocking rules for unary Pd(111) or Ag(111) caps the total coverage of O* at 1/3 ML (Fig. 5b), but the adsorption energy of the O* species is unaltered.

Fig. 5 Left column shows gross (grey background) and net (full colour) *OH/O* adsorption energy distributions and illustration of the adsorbate coverage when the adsorbates are blocking according to the mutual blocking rules for (a) unary Pd(111) or Ag(111), (c) intermetallic PdAg(111), and (e) equimolar solid solution PdAg(111). Right column shows gross (grey background) and net (full colour) *OH/O* adsorption energy distributions and illustration of the adsorbate coverage when the adsorbates are blocking in accordance with the extended blocking scheme for (b) unary Pd(111) or Ag(111), (d) intermetallic PdAg(111), and (f) equimolar solid solution PdAg(111).

Obtaining the net distribution on the intermetallic PdAg surface is also quite simple. Here, *OH on-top adsorption on Pd is the most stable species, and with both mutual blocking rules (Fig. 5c) and extended blocking rules (Fig. 5d), every Pd atom has an adsorbed *OH. The only difference between the two sets of rules is whether *OH can adsorb on Ag sites; these are, however, much less reactive compared to Pd sites and not active for ORR until a large overpotential is reached.

The net *OH/O* adsorption energy distributions for the solid solution surface are more interesting because, unlike unary Pd(111) and the intermetallic, they have both O* and *OH species at the optimal adsorption energies for ORR. With the mutual blocking rules (Fig. 5e), there is a dip in O* species at the optimal ΔΔG^Pd_O* ∼0.1 eV, but this is compensated by a large presence of *OH species within the same energy range. With extended blocking (Fig. 5f), these trends are further enhanced and the surface now has a majority of *OH species, of which many are located at the optimal energy for ORR.

3.4 ORR activity of different PdAg(111) alloy surface manifestations

We now convert the net *OH/O* adsorption energy distribution to the ORR activity at U = 0.86 V vs. RHE. We follow the approach of ref. 13 and estimate the kinetic limited current from each surface site, j_ki, as:

Furthermore, we cap the maximum current from each surface site to a diffusion-limited current (j_D) and obtain the total current per number of sites (j_tot) divided by the diffusion limit:

We note that imposing a diffusion limit at each site is a choice we make, and that we could instead have imposed a diffusion limit on the total current. At U = 0.86 V vs. RHE, optimal sites for ORR have the maximum j_ki current equal to 1. If we also set j_D equal to 1, j_tot/j_D has a maximum value of ½, and a lower value reflects how far a given surface is from the optimum.

In the first case, we assume that the only Pd(111) or Ag(111) surfaces exist on Pd_xAg_(1−x) alloys and that the amount of each surface correspond to the bulk composition. On Pd(111), ΔΔG^Pd_O* is almost 0 eV. This means that Pd(111) has a similar limiting potential to Pt(111) (at 0.74 and 0.76 V vs. RHE, respectively), and very similar ORR activity, which is also observed experimentally.⁶³ It also means that the Pd(111) activity is low at 0.86 V vs. RHE, with mutual blocking giving 3 times higher activity than extended blocking, because the O* coverage is 1 monolayer with mutual blocking and 1/3 with extended blocking. Still, Ag(111) has no discernible ORR activity, so the ORR activity increases linearly with Pd molar fraction (Fig. 6a).

Fig. 6 Predicted ORR activity at U = 0.86 V vs. RHE for (a) unary surfaces, (b) intermetallic, and (c) solid solution as a function of Pd content in the alloy.

In the second case, we also include the PdAg intermetallic (111) surface, such that at x = 0.5, there is only the PdAg intermetallic surface, while at x > 0.5, both the intermetallic and Pd(111) surfaces exist, and at x < 0.5, both the intermetallic and Ag(111) surfaces exist. This situation would occur if the only stable structures in the convex hull were Pd(111), the PdAg intermetallic, and Ag(111). The intermetallic has close to optimal *OH adsorption energies with *OH coverage up to 0.5 monolayer, so the predicted activity of 0.21 (both for mutual blocking and extended blocking) is many times higher than Pd(111) at U = 0.86 V vs. RHE (Fig. 6b). The activity reduction, from the 0.5 optimum to 0.21, is mainly due to half the surface being inert Ag atoms. The activity falls off linearly away from x = 0.5 towards both Ag(111) and Pd(111).

In the last case, we consider solid solutions, where the surface composition is equal to the bulk composition. The advantage of a solid solution is that some sites have the optimal adsorption energy for ORR, while the drawback is that these sites are fewer in numbers compared to the PdAg intermetallic. This means that the ORR activity is high over a large range of composition space, but also that it is never as high as the PdAg intermetallic. Another feature is that the solid solution can have abrupt changes in the limits of unary metals, as seen at x → 1. This is also seen in modelling of CO₂ hydrogenation to CO on CuAg solid solution (111) surfaces.⁶⁴

In Fig. 5, it is worth noting that the relative activities are very similar, irrespective of whether we apply mutual blocking or extended blocking rules. This is comforting, since we like our conclusions to be robust no matter our exact choice of rules. It means that the unoccupied sites with the extended blocking rules are not significantly different to the occupied sites; it mostly just means that we have fewer sites than with mutual blocking. Likewise, we imagine that the analysis would also yield similar results if we allowed for O* adsorption in the HCP hollow sites, since these sites adsorb O* similarly to the FCC hollow sites. It also indicates that the surface structures we end up with do not have to look exactly like the experimentally realized surfaces to make good models for the activity. The largest discrepancy between mutual blocking and extended blocking is at the x → 1 limit, so this seems a good experimental test case to improve our understanding of how adsorbate–adsorbate blocking influences catalysis.

The last part in the PdAg alloy analysis is to qualify if the rules we have imposed are resulting in predictions and consequences that are consistent with experimental observations. Our main findings are: (1) On unary Pd(111), O* is adsorbing at a lower potential than *OH. (2) *OH on top of Pd would be active for ORR; however, the surface is covered with O* leaving no room for adsorption on the active site. (3) By alloying with Ag, the oxygen coordinated directly with Ag is destabilized, and in some cases *OH on top of Pd becomes the most stable, which increases the activity. We have earlier investigated *CO stripping and oxidation on Pd and Pt experimentally in ref. 65. The observations are that on the Pt and Pd blank CV in water, the oxidation happens around the same potential, which is consistent with *OH adsorption on Pt, and O* adsorption on Pd. The CO stripping becomes possible at potentials where *OH can adsorb on the CO-covered surface. After the initial CO stripping, CO oxidation continues due to adsorption of additional CO from the solution. On Pt, it stops at the potential where *OH adsorption is no longer possible, as measured with the blank CV. On Pd, it stops at a potential 0.1 V higher than the oxidation on the blank CV and there is a small region with negative current after the CO oxidation has stopped and at higher potentials than the oxidation on the blank CV, consistent with a reduction of oxygen on the surface of Pd. Our interpretation is that only *OH can oxidize *CO and, on Pd, O* is more stable than *OH, which requires a 0.1 V higher potential. An adsorption 0.1 eV weaker than Pt is known to be optimal for oxygen reduction. Experiments with PdAg have shown improved ORR activity in an alkaline electrolyte compared to Pd.¹² For this to be consistent with calculated adsorption energies, we attribute this increased activity to *OH on Pd sites. These two experimental studies are indicating that the effects described here are real and important for catalytic activity.

It is interesting to note that nobody has described this when studying uniform surfaces; it is a phenomenon we first encountered in the simulations on HEA surfaces, as these complex and heterogeneous surfaces forced us to think of adsorbate–adsorbate interactions and gross and net distributions of adsorption energies. We now have several examples where the HEA surfaces reveal fundamental effects in adsorption, which have been overlooked for decades because we mostly have investigated uniform surfaces.

3.5 Net distributions and ORR activity for Ag₁₄Ir₁₆Pd₃₀Pt₁₄Ru₂₆ high-entropy alloy

In the last section, we extend the net distribution and ORR activity analysis to a Ag₁₄Ir₁₆Pd₃₀Pt₁₄Ru₂₆ HEA (111) surface and to the full potential range relevant for ORR rather than just 0.86 V vs. RHE. We note that the Ag₁₄Ir₁₆Pd₃₀Pt₁₄Ru₂₆ HEA composition is arbitrarily chosen, and we therefore focus on aspects that we believe are more generally valid for HEA catalysts.

The gross *OH and O* adsorption energy distributions for the Ag₁₄Ir₁₆Pd₃₀Pt₁₄Ru₂₆ (111) surface (Fig. 7a) show that the most stable species are *OH species and that many of the distributions are situated at ΔΔG^alloy_*OH/O* of −0.25 eV or below. The latter lowers the ratio of sites with ΔΔG^alloy_*OH/O* = ∼0.1 eV and negatively impacts the overall ORR activity; however, the ratio of optimal sites could be increased by adjusting the HEA composition.⁶⁶

Fig. 7 Ag₁₄Ir₁₆Pd₃₀Pt₁₄Ru₂₆ HEA (111) surface (a) gross *OH/O* adsorption energy distributions and net *OH/O* adsorption energy distributions obtained with (b) mutual blocking or (c) extended blocking rules. Potential-dependent ORR current of the PdAg intermetallic, the equimolar PdAg solid solution, and the Ag₁₄Ir₁₆Pd₃₀Pt₁₄Ru₂₆ HEA given that either (d) the gross distributions, (e) mutual blocking net distributions, or (f) extended blocking net distributions determine the ORR activity.

The net *OH and O* adsorption energy distributions are shown in Fig. 7b for mutual blocking and in Fig. 7c for extended blocking. The distributions show that *OH, to a large extent, outcompetes O* species on the surface, even though *OH species only seem to be slightly more favoured than O* species in the gross distributions. We also see that the net *OH distribution from extended blocking is significantly reduced at higher adsorption energies compared to mutual blocking, so it seems that the adsorption sites below ΔΔG^alloy_*OH/O* = −0.25 eV, are more damaging to the ORR activity when extended blocking rules apply.

Lastly, we estimate the current of the PdAg intermetallic, the equimolar PdAg solid solution, and the Ag₁₄Ir₁₆Pd₃₀Pt₁₄Ru₂₆ HEA, for the full potential range relevant for ORR, given that either the gross distributions (Fig. 7d), mutual blocking net distributions (Fig. 7e), or extended blocking net distributions (Fig. 7f) determine the ORR activity. The ORR activity is lowest for the Ag₁₄Ir₁₆Pd₃₀Pt₁₄Ru₂₆ HEA, no matter what rules apply. This shows that if we already have an active bimetallic alloy, changing to a HEA will generally just dilute the number of active sites. The gross distributions contain both one monolayer of *OH and one monolayer of O* so the combined current reaches j_tot/j_D = 2, which is unphysical. Interestingly, both the PdAg intermetallic and the PdAg solid solution currents contain features that differ between the different applied rules, e.g. the presence (mutual blocking) or absence (extended blocking) of a secondary uptake in the PdAg intermetallic current at 0.5 V vs. RHE related to the presence or absence of *OH on the Ag atoms in the surface. Experiments could validate certain rules based on such features. In contrast, except for a scaling factor, the current from Ag₁₄Ir₁₆Pd₃₀Pt₁₄Ru₂₆ has the same features everywhere in the potential range, independent of the applied rules. This makes it difficult to validate specific rules by experiments, but it also means that there is little reason to construct net distributions, since the gross distributions are likely adequate descriptors for the ORR current. This latter finding is in line with previous findings from our group.⁶⁷

4 Conclusions

Including the adsorbate–adsorbate interaction is a frontier in modelling adsorption and catalysis on complex surfaces. It seems impossible to include it in a straightforward manner; therefore, we need to assume some rules for the adsorbate–adsorbate interaction that are very similar to the rules for a Langmuir isotherm. We impose two slightly different sets of rules resulting in different representations of the covered surface, but yielding very similar trends in activity. We identify that the adsorbate–adsorbate interaction can affect the activity positively by allowing for adsorption on-top Pd atoms, which on a unary Pd surface would be blocked by O* adsorption in the FCC hollow site. There are experimental indications that this can be a real effect that is important for catalysis. This is one of a few phenomena we have found by investigating HEA surfaces. In this study, it is very easy to challenge the assumptions; however, we see that as a strength compared to having the assumptions enter in a less explicit manner in the models.

Data availability

All DFT data for this article are from earlier published work.¹⁴ The scripts needed to construct the machine learning regression model are freely accessible at our Github repository: https://github.com/catalyticmaterials/CHEAT. The scripts necessary to reproduce the results and figures presented herein are available at: https://erda.ku.dk/archives/a4aaf5f6badf4aada8daebcd1896bdde/published-archive.html.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We acknowledge the Danish National Research Foundation Centers of Excellence, the Center for High Entropy Alloy Catalysis (Project DNRF149) and financial support from the European Union under ERC Synergy grant DEMI, GA no. 101118768.

References

1. I. Langmuir, J. Am. Chem. Soc., 1918, 40, 1361–1403.
2. G. S. Karlberg, T. F. Jaramillo, E. Skúlason, J. Rossmeisl, T. Bligaard and J. K. Nørskov, Phys. Rev. Lett., 2007, 99, 126101.
3. J. Rossmeisl, G. S. Karlberg, T. Jaramillo and J. K. Nørskov, Faraday Discuss., 2009, 140, 337–346.
4. C. Liang, R. R. Rao, K. L. Svane, J. H. L. Hadden, B. Moss, S. B. Scott, M. Sachs, J. Murawski, A. M. Frandsen, D. J. Riley, M. P. Ryan, J. Rossmeisl, J. R. Durrant and I. E. L. Stephens, Nat. Catal., 2024, 7, 763–775.
5. J. Rossmeisl, K. D. Jensen, A. S. Petersen, L. Arnarson, A. Bagger and M. Escudero-Escribano, J. Phys. Chem. C, 2020, 124, 20055–20065.
6. N. M. Marković, T. J. Schmidt, V. Stamenković and P. N. Ross, Fuel Cells, 2001, 1, 105–116.
7. K. Kodama, A. Shinohara, N. Hasegawa, K. Shinozaki, R. Jinnouchi, T. Suzuki, T. Hatanaka and Y. Morimoto, J. Electrochem. Soc., 2014, 161, F649.
8. Q. He, B. Shyam, M. Nishijima, D. Ramaker and S. Mukerjee, J. Phys. Chem. C, 2013, 117, 4877–4887.
9. M. K. Plenge, J. K. Pedersen, A. Bagger and J. Rossmeisl, J. Catal., 2024, 430, 115322.
10. J. K. Nørskov, T. Bligaard, A. Logadottir, S. Bahn, L. B. Hansen, M. Bollinger, H. Bengaard, B. Hammer, Z. Sljivancanin, M. Mavrikakis, Y. Xu, S. Dahl and C. J. H. Jacobsen, J. Catal., 2002, 209, 275–278.
11. J. K. Nørskov, T. Bligaard, J. Rossmeisl and C. H. Christensen, Nat. Chem., 2009, 1, 37–46.
12. J. A. Zamora Zeledón, M. B. Stevens, G. T. K. K. Gunasooriya, A. Gallo, A. T. Landers, M. E. Kreider, C. Hahn, J. K. Nørskov and T. F. Jaramillo, Nat. Commun., 2021, 12, 620.
13. T. A. A. Batchelor, T. Löffler, B. Xiao, O. A. Krysiak, V. Strotkötter, J. K. Pedersen, C. M. Clausen, A. Savan, Y. Li, W. Schuhmann, J. Rossmeisl and A. Ludwig, Angew. Chem., Int. Ed., 2021, 60, 6932–6937.
14. C. M. Clausen, M. L. S. Nielsen, J. K. Pedersen and J. Rossmeisl, High Entropy Alloys & Materials, 2023, 1, 120–133.
15. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864–B871.
16. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133–A1138.
17. B. Hammer, L. B. Hansen and J. K. Nørskov, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 7413–7421.
18. J. J. Mortensen, L. B. Hansen and K. W. Jacobsen, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 035109.
19. J. Enkovaara, C. Rostgaard, J. J. Mortensen, J. Chen, M. Dułak, L. Ferrighi, J. Gavnholt, C. Glinsvad, V. Haikola, H. A. Hansen, H. H. Kristoffersen, M. Kuisma, A. H. Larsen, L. Lehtovaara, M. Ljungberg, O. Lopez-Acevedo, P. G. Moses, J. Ojanen, T. Olsen, V. Petzold, N. A. Romero, J. Stausholm-Møller, M. Strange, G. A. Tritsaris, M. Vanin, M. Walter, B. Hammer, H. Häkkinen, G. K. H. Madsen, R. M. Nieminen, J. K. Nørskov, M. Puska, T. T. Rantala, J. Schiøtz, K. S. Thygesen and K. W. Jacobsen, J. Phys.: Condens. Matter, 2010, 22, 253202.
20. H. J. Monkhorst and J. D. Pack, Phys. Rev. B, 1976, 13, 5188–5192.
21. J. D. Pack and H. J. Monkhorst, Phys. Rev. B, 1977, 16, 1748–1749.
22. S. R. Bahn and K. W. Jacobsen, Comput. Sci. Eng., 2002, 4, 56–66.
23. A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng and K. W. Jacobsen, J. Phys.: Condens. Matter, 2017, 29, 273002.
24. L. Vegard, Z. Phys., 1921, 5, 17–26.
25. C. M. Clausen, J. K. Pedersen, T. A. A. Batchelor and J. Rossmeisl, Nano Res., 2022, 15, 4775–4779.
26. N. Garcia-Araez, J. Phys. Chem. C, 2011, 115, 3075–3082.
27. A. S. Bandarenka, H. A. Hansen, J. Rossmeisl and I. E. L. Stephens, Phys. Chem. Chem. Phys., 2014, 16, 13625–13629.
28. I. C. Man, H.-Y. Su, F. Calle-Vallejo, H. A. Hansen, J. I. Martínez, N. G. Inoglu, J. Kitchin, T. F. Jaramillo, J. K. Nørskov and J. Rossmeisl, ChemCatChem, 2011, 3, 1159–1165.
29. J. K. Pedersen, C. M. Clausen, L. E. J. Skjegstad and J. Rossmeisl, ChemCatChem, 2022, 14, e202200699.
30. J. K. Nørskov, J. Rossmeisl, A. Logadottir, L. Lindqvist, J. R. Kitchin, T. Bligaard and H. Jónsson, J. Phys. Chem. B, 2004, 108, 17886–17892.
31. A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai and S. Chintala, in Advances in Neural Information Processing Systems, Curran Associates, Inc., 2019, vol. 32.
32. I. Loshchilov and F. Hutter, arXiv, 2019, preprint, arXiv:1711.05101,  DOI:10.48550/arXiv.1711.05101.
33. T. Löffler, A. Ludwig, J. Rossmeisl and W. Schuhmann, Angew. Chem., Int. Ed., 2021, 60, 26894–26903.
34. J. R. Kitchin, K. Reuter and M. Scheffler, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 075437.
35. G. P. Schwartz, Surf. Sci., 1978, 76, 113–129.
36. F. J. Kuijers and V. Ponec, J. Catal., 1979, 60, 100–109.
37. B. I. Podlovchenko, Y. M. Maksimov and A. G. Utkin, Russ. J. Electrochem., 2015, 51, 891–898.
38. L. A. Cipriano, H. H. Kristoffersen, R. L. Munhos, R. Pittkowski, M. Arenz and J. Rossmeisl, Nanoscale, 2023, 15, 16697–16705.
39. A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder and K. A. Persson, APL Mater., 2013, 1, 011002.
40. S. Hayat, A. B. Ziya, N. Ahmad and F. Bashir, Phys. Solid State, 2020, 62, 54–58.
41. I. Karakaya and W. T. Thompson, Bull. Alloy Phase Diagrams, 1988, 9, 237–243.
42. J. Sopoušek, A. Zemanová, J. Vřešťál and P. Brož, J. Alloys Compd., 2010, 504, 431–434.
43. A. S. Pavlenko, E. A. Ptashkina, E. G. Kabanova, G. P. Zhmurko, M. A. Kareva and V. N. Kuznetsov, Calphad, 2023, 81, 102533.
44. I. E. L. Stephens, A. S. Bondarenko, F. J. Perez-Alonso, F. Calle-Vallejo, L. Bech, T. P. Johansson, A. K. Jepsen, R. Frydendal, B. P. Knudsen, J. Rossmeisl and I. Chorkendorff, J. Am. Chem. Soc., 2011, 133, 5485–5491.
45. C. M. Clausen, T. A. A. Batchelor, J. K. Pedersen and J. Rossmeisl, Adv. Sci., 2021, 8, 2003357.
46. M. Mavrikakis, B. Hammer and J. K. Nørskov, Phys. Rev. Lett., 1998, 81, 2819–2822.
47. P. Strasser, S. Koh, T. Anniyev, J. Greeley, K. More, C. Yu, Z. Liu, S. Kaya, D. Nordlund, H. Ogasawara, M. F. Toney and A. Nilsson, Nat. Chem., 2010, 2, 454–460.
48. M. Escudero-Escribano, P. Malacrida, M. H. Hansen, U. G. Vej-Hansen, A. Velázquez-Palenzuela, V. Tripkovic, J. Schiøtz, J. Rossmeisl, I. E. L. Stephens and I. Chorkendorff, Science, 2016, 352, 73–76.
49. H. H. Kristoffersen and J. Rossmeisl, J. Phys. Chem. C, 2022, 126, 6782–6790.
50. M. Polak and L. Rubinovich, Surf. Sci. Rep., 2000, 38, 127–194.
51. H. A. Hansen, J. Rossmeisl and J. K. Nørskov, Phys. Chem. Chem. Phys., 2008, 10, 3722–3730.
52. V. Tripkovic and T. Vegge, J. Phys. Chem. C, 2017, 121, 26785–26793.
53. A. Hitotsuyanagi, M. Nakamura and N. Hoshi, Electrochim. Acta, 2012, 82, 512–516.
54. M. F. Li, L. W. Liao, D. F. Yuan, D. Mei and Y.-X. Chen, Electrochim. Acta, 2013, 110, 780–789.
55. V. R. Stamenkovic, B. Fowler, B. S. Mun, G. Wang, P. N. Ross, C. A. Lucas and N. M. Marković, Science, 2007, 315, 493–497.
56. V. Climent, R. Gómez, J. M. Orts and J. M. Feliu, J. Phys. Chem. B, 2006, 110, 11344–11351.
57. M. Wakisaka, H. Suzuki, S. Mitsui, H. Uchida and M. Watanabe, Langmuir, 2009, 25, 1897–1900.
58. V. Tripković, E. Skúlason, S. Siahrostami, J. K. Nørskov and J. Rossmeisl, Electrochim. Acta, 2010, 55, 7975–7981.
59. M. F. Nygaard, M. L. S. Nielsen and J. Rossmeisl, Angew. Chem., Int. Ed., 2025, 64, e202417308.
60. J. Rossmeisl, J. Greeley and G. s. Karlberg, in Fuel Cell Catalysis, John Wiley & Sons, Ltd, 2009, pp. 57–92.
61. L. C. Grabow, B. Hvolbæk and J. K. Nørskov, Top. Catal., 2010, 53, 298–310.
62. C. T. Campbell, G. Ertl, H. Kuipers and J. Segner, Surf. Sci., 1981, 107, 220–236.
63. S. Kondo, M. Nakamura, N. Maki and N. Hoshi, J. Phys. Chem. C, 2009, 113, 12625–12628.
64. D. Yan, H. H. Kristoffersen, I. E. Castelli and J. Rossmeisl, J. Phys. Chem. C, 2022, 126, 19756–19765.
65. M. Paula Salinas-Quezada, J. K. Pedersen, P. Sebastián-Pascual, I. Chorkendorff, K. Biswas, J. Rossmeisl and M. Escudero-Escribano, EES Catal., 2024, 2, 941–952.
66. T. A. A. Batchelor, J. K. Pedersen, S. H. Winther, I. E. Castelli, K. W. Jacobsen and J. Rossmeisl, Joule, 2019, 3, 834–845.
67. C. M. Clausen, O. A. Krysiak, L. Banko, J. K. Pedersen, W. Schuhmann, A. Ludwig and J. Rossmeisl, Angew. Chem., Int. Ed., 2023, 62, e202307187.

---