Single atom catalysts on amorphous supports: A quenched disorder perspective

Catalysts are critical in the efficient transformation of conventional and alternative hydrocarbon sources to fuels and chemical products, as well as in new energy technologies such as solar fuels and fuel cells. Many important and long-studied catalysts involve isolated metal centers dispersed on amorphous supports such as silica and silica-alumina. Examples include Mo(VI), W(VI), and Re(VII) sites supported on a variety of oxides for olefin metathesis,^1–3 Cr(VI)/silica and Cr(VI)/AlPO₄ catalysts for olefin polymerization,^4,5 Ti(IV)/silica catalysts for olefin epoxidation,^6,7 and oxide-supported V(V) catalysts for methanol partial oxidation.^8,9 The active sites of these catalysts are often represented with skeletal formulas like those in Figure 1.

For a homogeneous catalyst, the skeletal formula corresponds to the same molecular geometry for each catalyst site. For sites comprised of a metal atom grafted to an amorphous support, each site is slightly different from the others. Many catalysts made in this way are known to have very low fractions of active sites. Nevertheless, some closely related materials involving amorphous supports have recently been described as single-site catalysts.^10–14 While this is appropriate in the sense that each site is independent of the others, amorphous supports present many slightly different grafting environments for metal ions, so the sites are unlikely to have identical reactivities. Taylor and others recognized very early that even small structural differences between sites can lead to important activity differences.^15–17 In addition to their activity, the differences between sites may also influence adsorption, selectivity, and even cause rate-limiting steps to differ between sites.¹⁸ Similarly, Langmuir anticipated non-uniformity among adsorption sites and suggested that isotherms should be summed empirically over individual sites (i) with varying adsorption constants K_i to obtain a generalized coverage θ = Σ_iK_iP/(1 + K_iP).¹⁹ 

Under vacuum, desorption from a surface occurs once per site as the temperature increases, so the number of first-desorption events can be counted systematically at each temperature. Consequently, temperature-programmed desorption (TPD) can directly probe the distribution of desorption energies for an ensemble of non-uniform sites.^20–23 Using a similar approach, van Bokhoven et al. characterized the activation energies for decomposition of pre-adsorbed sarin on γ-alumina.^24,25 The study by van Bokhoven et al. was also made possible by the once-per-site nature of the decomposition process.

One might anticipate that similar approaches could characterize the distribution of steady-state rates on an ensemble of sites. However, characterizing the distribution of reactivity is fundamentally more difficult because, at each temperature, the steady-state rate is a superposition of the activities for all sites. New imaging techniques can identify highly active sites (or patches of active sites) as sources of products, and these tools may become important for single atom catalysts on amorphous supports.^26–28 Traditionally, the only direct probes of non-uniform activity were active-site counting techniques: selective poisoning, radiolabeling, stopped-flow techniques, and reaction-desorption studies like that of van Bokhoven et al.^24,25,29 Even these techniques rely on uncertain (and sometimes questionable) assumptions, e.g., that poisons bind selectively to the most active sites. Nevertheless, carefully designed active-site counting experiments often suggest that low fractions, ca. 1%-10%, of sites are responsible for most of the catalytic activity.^30–32 Catalysts with low fractions of active sites pose challenges for structural characterization methods, which report primarily on the structures of common sites instead of the rare and highly active sites.³³ In the context of the Phillips catalyst, McDaniel noted that such spectra may be highly misleading in efforts to identify mechanisms.⁴ 

Ab initio computational studies of amorphous catalysts are also prohibitively difficult. In contrast to crystalline catalyst materials like zeolites, metals, and ordered oxides where computation has dramatically advanced our understanding of structure-activity relationships,^34–41 computational advances for amorphous catalysts have been limited by a lack of detailed structural information. We and others have used small cluster models containing specific ring sizes to model amorphous silica and silica-alumina.^3,42–47 Other investigators have melted, quenched, and cut silicas and then capped them with hydroxyl groups to model the hydroxyl-terminated surfaces of amorphous silica.^48,49 Some investigators have used crystalline structures like zeolites⁴⁴ or β-cristobalite^50,51 to model amorphous silica-alumina and silica, respectively.

To some degree, these models will always differ from real amorphous supports and from the processes by which they are made. Typically, one is not certain which reactions will be affected by such structural differences and which reactions will not. Recently, Goldsmith et al. developed an algorithm that quantifies the sensitivity of activation energies to structural differences between active sites without relying on any specific structural model.^52,53 Rigid constraints on the peripheral atoms of a small cluster model were included to mimic connections to a solid matrix. Previous approaches also included rigid constraints, but the peripheral atom positions were based on the structures of crystalline materials or on ad hoc considerations. Real amorphous supports present a quenched distribution of peripheral atom positions corresponding to a multitude of possible site environments. Since exhaustively sampling the quenched distribution is not yet possible, the algorithm of Goldsmith et al. instead generates the structures of minimum energy active sites for each activation energy across the range of activation energies.^52,53

The algorithm further assumes that site populations decrease with increasing site energy.^52,53 The assumption seems reasonable at first inspection but raises an important question: which energies govern these site populations? Candidates include energies of the model sites before the catalyst was grafted, energies of the bare catalyst sites after grafting, energies of the sites in the most abundant surface intermediate (MASI) state, and energies of the sites at some critical point in the catalyst preparation history.

This troublesome question could be avoided by beginning with an atomistically detailed model of the amorphous support. Fortunately, structural models are advancing in their complexity, in their fidelity to the actual preparation history, in the quality of the model chemistry, and in the diversity of structural features represented. Ugliengo,⁵⁴ Johnson,⁵⁵ and others⁵⁶ are working toward realistic ab initio models of amorphous silica with variable hydroxyl contents. Adaptive kinetic Monte Carlo⁵⁷ and event-cataloging techniques^58,59 have the potential to simulate hydroxyl condensations and structural rearrangements during slow temperature quenches. On the experimental side, Freund and others^60–63 have imaged thin layers of amorphous silica and aluminosilicate with atomic-scale resolution. These new imaging capabilities may soon enable direct comparisons between theory and catalysis experiments on surfaces with precisely known structural features.

This letter provides a statistical framework for efficiently computing quenched averages⁶⁴ of kinetic properties such as the overall rate, the activation energy, and the fraction of active sites. The formalism also provides an interpretation of the kinetics in terms of an effective synthesis temperature and a separate catalyst operation temperature. Finally, the statistical framework quantitatively prioritizes those aspects of future calculations that will most strongly influence the predictions.

We consider only the simple case in which each active site on a catalyst has the same rate-limiting step and the same most abundant intermediate state, and does not deactivate. Suppose that each site in a quenched structural distribution of sites has a local geometry x and a rate law with power-law dependence on concentrations C_i.^65,66 The turnover frequency at each site is $k ( x ) ∏ i C i α$⁠, i.e., a product of concentrations and a site-dependent (non-elementary) rate constant of Arrhenius form⁶⁷ 

where β = 1/k_BT, and where the prefactor A(x) may depend on geometry but only weakly depends on temperature. ΔE(x) is the activation energy for the site with structure x and not necessarily the site-averaged activation energy. All of the analysis below could be done instead for a transition-state theory rate constant,⁶⁸ but since catalytic reactions are never elementary, the overall rate does not precisely follow transition-state theory. The rate constant averaged over all sites is

where ρ(x) is the normalized probability density for sites of specific geometry x, and where the integration with unspecified bounds is an average $$ over all site geometries. Wu et al. used Eq. (2) for equilibrium site averages in which different site activities emerge from different patterns of neighboring adsorbates on a Pt(111) surface.⁶⁹ In this work, the catalyst structure ρ(x) is a quenched-in non-equilibrium distribution. It cannot be obtained with currently existing computational procedures. We assume that ρ(x) is irreversibly set by the temperature and other conditions during catalyst synthesis. We assume that, for high-melting materials like silica and silica-alumina, the distribution ρ(x) does not depend on the temperature T used for catalyst operation.

The distribution ρ(x) can be formally projected onto ΔE to obtain

We can also define k(ΔE) as a conditional average of the rate constant over all sites with activation energy ΔE,

From $ρ ̃ ( Δ E )$ and k(ΔE), all remaining kinetic properties are easily determined. For example, the site-averaged rate constant is

In this statistical framework, $ρ ̃ ( Δ E )$ is a density of sites at activation energy ΔE, analogous to the density of states in equilibrium statistical mechanics. Not surprisingly then, the calculations of certain properties resemble familiar calculations from equilibrium statistical mechanics.⁷⁰ For example, the Arrhenius activation energy is

where the new, k-weighted average for any property Y is

The observed (k-weighted) activation energy in (6) may be very different from the activation energy of a typical site. Additionally, the illustrative example below shows that the statistical framework can also quantify the fraction of active sites.

Suppose that a catalyst has a Gaussian (G) distribution of activation energies ρ_G(ΔE) with mean ΔE₀ and standard deviation σ. Equation (5) gives 〈k〉_G = Aexp[ − βΔE₀ + β²σ²/2]. The average rate is not the rate at the average activation energy. Instead, the apparent activation energy is

Equation (8) highlights the importance of at least two temperatures: the temperature during catalyst synthesis, which influences $ρ ̃ ( Δ E )$⁠, and the temperature during catalyst operation which determines the activation barriers that can be surmounted. To illustrate realistic sizes of the ΔE₀ and σ parameters, we extract $ρ ̃ ( Δ E )$ from the site energy vs. activation energy profile for the Mo-catalyzed metallacycle rotation examined by Goldsmith et al.⁵³ If the synthesis temperature is 700 K, then a Boltzmann distribution applied to the previous results gives an activation energy distribution with mean ΔE_μ = 53.1 kJ/mol and standard deviation σ = 11.5 kJ/mol. We stress that this is only an illustrative calculation. The synthesized catalyst need not correspond to an equilibrium structure at any temperature. Now, suppose the catalyst is operated near 450 K. Equation (8) predicts that the site-averaged activation energy is only 20 kJ/mol. This example illustrates the degree to which sites on the low-activation energy tail of the distribution dominate the reactivity: in fact, it suggests that sites three standard deviations below the mean activation energy dominate the observed kinetics. More broadly, it suggests that the effects of quenched disorder might be even larger than the more widely appreciated effects of different density functionals and basis sets.

The Gaussian example also illustrates the need to thoroughly characterize the low energy tail of the activation energy distribution. Goldsmith et al. did not actually identify any site with activation energies below 30 kJ/mol, yet the Gaussian approximation says they should be important. Do sites with activation energies three standard deviations below the average even exist? In general, the rate-weighted distribution $k − 1 k ( Δ E ) ρ ̃ ( Δ E )$ emphasizes low activation energy tails of the distribution, making this tail critical for accurate predictions. For example, consider a uniform block (B) distribution with the same mean and standard deviation as the Gaussian (G) distribution. The corresponding block distribution is $ρ ̃ B ( Δ E ) = ( 2 σ 3 ) − 1$ for $Δ E μ −σ 3 <ΔE<Δ E μ +σ 3$⁠, and $ρ ̃ B ( Δ E ) =0$ otherwise. The site-averaged rate constant is $k B =Aexp [ − β Δ E μ ] sinh [ β σ 3 ] ( β σ 3 ) − 1$⁠, which depends very differently on the operating temperature compared to the rate that emerged from the Gaussian distribution. Figure 2 shows that the difference in shape causes a large difference in apparent activation energies for the Gaussian- and block-shaped models.

The statistical framework can also be used to predict the results of active site counting experiments. The fraction of sites that are active can be expressed as the fraction of sites f_X% which are responsible for X% of the activity. This quantity can be computed in two steps,

which identifies the appropriate region ΔE < ΔE_X% of the activation energy tail as shown in Figure 3. A second integration gives the corresponding fraction of active sites in the tail of the distribution $ρ ̃ B ( Δ E )$⁠,

Figures 2 and 3 demonstrate the need for more powerful methods to probe the lowest reaches of the activation energy distribution computationally. Unless we can accurately predict how rapidly $ρ ̃ B ( Δ E )$ decays at low activation energies, we cannot predict overall kinetics, activation energies, or the results of active-site counting experiments. The framework developed here can potentially be generalized to predict averaged kinetics from rate expressions like those that emerge from Langmuir-Hinshelwood-Hougen-Watson analyses.⁶⁶ 

Finally, the statistical framework in this letter can help reveal which aspects of a computational model are most critical for obtaining accurate predictions. For example, relative errors in the averaged rate constant, $δln k = k − 1 δ k$⁠, can be expanded (to first order) in the errors δρ(x), δA(x), and δΔE(x). As might be expected, $k$ is most sensitive to errors in the mapping from x to ΔE(x): $δln k = δ ln ρ ( x ) k + δ ln A ( x ) k − δ β Δ E ( x ) k$⁠. Perhaps less expected, the formalism also shows that all factors are primarily important near the peak in the k-weighted distribution, with errors in other locations being inconsequential.

For appropriate systems, it may eventually be possible to relax the assumption that the distribution is quenched-in by the synthesis protocol. The actual distribution of site structures might evolve over time, likely affecting a type of catalyst deactivation. One could, in principle, augment the present framework with time-dependent population balance equations⁷¹ for time-dependent distributions ρ(x, t) or for ρ(ΔE, t). However, the more immediate challenge, as emphasized in this work, is to understand the initial distributions ρ(x, 0) and ρ(ΔE, 0).

Atomistically detailed models of amorphous materials and surfaces are becoming available, and quantitative ab initio studies of catalysis on the quenched-distribution of surface sites will follow. We have outlined (1) how the distribution of activation energies can be constructed from a distribution of active site structures, (2) how apparent activation energies can be extracted from the activation energy distribution, and (3) how active site counting results can be predicted computationally. The statistical framework shows that the low activation energy tail in the activation energy distribution is particularly important for predicting kinetic properties. This letter should both motivate and guide the development of new computational algorithms and the modeling of catalysts with isolated sites on amorphous supports.

The authors thank the Department of Energy, Basic Energy Sciences for support under Grant No. DE-FG02-03ER15467.