Stability of heterogeneous single-atom catalysts: a scaling law mapping thermodynamics to kinetics

Heterogeneous single-atom catalysts (SACs) hold the promise of combining high catalytic performance with maximum utilization of often precious metals. We extend the current thermodynamic view of SAC stability in terms of the binding energy (Ebind) of single-metal atoms on a support to a kinetic (transport) one by considering the activation barrier for metal atom diffusion. A rapid computational screening approach allows predicting diffusion barriers for metal–support pairs based on Ebind of a metal atom to the support and the cohesive energy of the bulk metal (Ec). Metal–support combinations relevant to contemporary catalysis are explored by density functional theory. Assisted by machine-learning methods, we find that the diffusion activation barrier correlates with (Ebind)2/Ec in the physical descriptor space. This diffusion scaling-law provides a simple model for screening thermodynamics to kinetics of metal adatom on a support. A rapid screening approach is developed to access the stability of single-atom catalysts based on the binding energy of the metal atom to the support and the cohesive energy of the bulk metal.


INTRODUCTION
Supported metals constitute an important class of heterogeneous catalysts widely used by the chemical and automotive industry, in controlling environmental pollution, and in energy conversion technology [1][2][3][4] . The importance of strong metal-support interactions has long been recognized 5 . Of special interest is the stabilization of single-metal atoms on a support, which is rapidly becoming a new frontier in materials chemistry and catalysis [6][7][8][9][10][11] . The single-atom nature of the catalysts can overcome the scarcity of particular elements, such as precious group metals, ubiquitous in catalysis 1,8,12 .
A general requirement for successful supported metal catalysts is high stability of the active phase in terms of the exposed metal surface. Sintering through Ostwald ripening leads to a reduction of the number of exposed metal atoms and is a common deactivation pathway of heterogeneous catalysts [13][14][15] . Diffusion into lattice sites of a strongly interacting support with defects can also lead to catalyst deactivation 16 . Obtaining stable singleatom catalysts (SACs) is a particular challenge and requires thorough knowledge about metal-support interactions. Several studies attempted to describe SAC stability in terms of the binding energy (E bind ) of a single-metal atom to the support [17][18][19] . All these investigations presume that stronger binding of a metal on a support will make it less prone to sintering. However, the diffusion activation barrier (E a ) of a metal atom on a support is the most relevant factor to its stability. Although several reports dealt with diffusion pathways of metal atoms on a support 20,21 , the thermodynamic and kinetic aspects of stability of SACs have not been systematically explored yet.
Stable SACs should exhibit diffusion barriers substantially higher than the thermal energy to avoid rapid sintering under realistic reaction conditions 22 . Despite significant advances in computing power, the calculation of activation barriers remains more demanding than that of binding energies of reaction intermediates. Moreover, some supports such as oxides with specific magnetic properties (mostly Cr, Mn, Fe, and Co-based systems, e.g., LaFeO 3 perovskite) pose significant challenges in convergence [23][24][25][26] , making investigations at the density functional theory (DFT) level often very expensive. Thus, it would be desirable to identify correlations between the diffusion activation barrier and the binding strength to assess the stability of a wide range of SACs. A universal scaling relation that holds for diverse metal-support combinations would be preferred.
In the present study, we aim at identifying a correlation between the kinetic stability (E a ) and the thermodynamic stability (E bind ) of a range of single transition metal atoms on various supports. As there is usually no (or little) change in the energy of the initial and final states for diffusion events between adjacent similar adsorption sites, the Brønsted-Evans-Polanyi principle is not applicable. We turn to machine-learning (ML) methods to identify relevant physical descriptors and draw new correlations. We choose two reducible metal oxides (CeO 2 , TiO 2 ), two stable metal oxides (MgO, ZnO), a perovskite (SrTiO 3 ), and two 2-dimensional materials (MoS 2 and graphene) as shown in Supplementary Fig. 1, relevant to many topics in contemporary heterogeneous catalysis. Stepped CeO 2 (111) was included as a support model with a corrugated and more reactive surface. A total of 11 transition metals (Cu, Ag, Au, Ni, Pd, Pt, Co, Rh, Ir, Fe, and Ru) were included as catalytically active single-atom centers, resulting in a dataset of 99 points. We identify a universal correlation between E a with the cohesive energy (E c ) of bulk metal and E bind of metal atoms to a support. This scaling law provides an accurate description of the stability of SACs and allows rapid screening of other systems solely on the basis of easily computable physical descriptors.

DFT investigation
We computed E bind and corresponding E a , at the generalized gradient approximation (GGA)-DFT level (Perdew-Burke-Ernzerhof (PBE) functional). In this study, we focused on the diffusion of metal adatoms between two adjacent same sites on idealized support surfaces without vacancies, decoration by functional groups such as hydroxyls or other chemical impurities. The corresponding values are collected in the Supplementary Information (Supplementary Table 1). The adopted binding energies correspond to the most stable adsorption configurations determined by considering different adsorption sites of metal SAs on the supports and possible spin states where necessary (see Supplementary Table 2). Previously, it has been suggested that E bind is a suitable descriptor for the stability of SACs 17,27 . Obviously, a stronger metal-support interaction reflected by a higher binding energy implies a lower mobility on the support. However, Fig. 1 shows that E a correlates only moderately (R 2 = 0.83) with E bind . Graphene and the basal planes of MoS 2 only weakly interact with metal atoms and, therefore, will not likely yield stable SACs. The correlation is especially poor for supports that exhibit relatively strong interaction with single-metal atoms. Therefore, such a linear correlation is not suitable for making useful predictions for a wider range of single-metal atom/support pairs. Our goal is to identify a more accurate correlation to improve predictability and general applicability for the diffusion barrier of metal-support systems.
Machine-learning study We applied three different ML algorithms commonly used for regularization and descriptor selection 28,29 , namely ridge 30 , least absolute shrinkage and selection operator (LASSO) 31 , and elastic net 32 regression. As an additional descriptor, we selected the (computed) cohesive energy (E c ) of the metal, because this parameter is an intrinsic metal property, representing the reactivity of the metal. Moreover, E c can be looked up for all transition metals and is in principle also easily computed accurately at the DFT level. Previous work explored the description of E bind in terms of parameters such as E c 17,27,33 . We assessed the dependence of E bind on E c , as shown in Supplementary Fig. 2. Although on a specific surface, E bind generally increases with E c , there is no strong correlation between these two parameters when other surfaces and materials are also considered. This lack of correlation is because E bind strongly depends on other properties of the metal-support systems, such as the ionization energy of the bulk metal and the work function and oxygen vacancy formation energy of support surfaces. We found that linear correlations between the two primary descriptors and E a were not satisfactory either (see Supplementary Fig. 3). Therefore, we populated the hypothesis space by introducing polynomial and natural logarithm terms of the primary descriptors (E c and E bind ) as secondary descriptors. The mathematical form of the model contains 87 unique descriptors x j for j = 1,…, p = 87 in the hypothesis space (see the complete list in Supplementary Tables 3 and 4), and is written as, The magnitude of the regression coefficients weights across models given the standardized dataset are visualized in the heat map shown in Fig. 2 and Supplementary Fig. 4. The ridge, LASSO, and elastic net algorithms select descriptors containing (E bind ) 2 as the most weighted and informative descriptors. In this work, all models are trained using 80% of the dataset (the training set). For the models requiring hyperparameter tuning including LASSO, ridge, and elastic net, a tenfold cross validation with ten repeats on the training set is used to prevent overfitting. The quality of a predictive model is usually evaluated through the prediction accuracy of future data and the interpretability of the model. Therefore, the root mean square error (RMSE) of the 20% dataset withheld (the testing set, not used in training of the models) is calculated as the best discriminant of model performance as  shown in Supplementary Fig. 5. In addition, we use R 2 values of all data to describe the general goodness-of-fit of the model. Instead of populating the descriptors manually, genetic programming (GP) based on symbolic regression can also be used to efficiently explore the hypothesis space 34 . By running the GP program with different random seeds, we found that the fittest model with the lowest testing RMSE of 0.262 eV takes a particularly simple form: Thus, both the feature-selection algorithms and GP identify (E bind ) 2 /E c as the most significant descriptor for the diffusion barrier of a single-metal atom. R 2 given by Eq. (2) is as high as 0.93, indicating a strong correlation. Next, we employ (E bind ) 2 /E c as the sole descriptor in an ordinary least square fitting procedure to further reduce the error (see Supplementary Tables 5 and 6). The resulting correlation with a testing RMSE of 0.220 eV takes the following form Diffusion scaling-law for single-atom catalysts (DSL-SAC) The strong correlation in the underlying computational data is visualized in Fig. 3a, and the strong quadratic dependency of E a on E bind is shown in Fig. 3b. While Eq. 3 does not necessarily imply an underlying physical concept, it is useful to consider the correlation in the form The ratio σ (E bind /E c ), which lies between 0.16 and 1.06 for all metal-support combinations explored in this study, measures the binding energy of the metal as an isolated atom to the support surface with respect to the intrinsic binding strength of the metal atom in bulk metal (E c ), and can be considered as a correction factor to E bind . For a reactive surface like CeO 2 (100), σ lies between 0.73 and 1.04 and thus the correction is small. On the other hand, the correction is much larger, i.e., σ lies between 0.16 and 0.34, when the same metals are on a more inert support, e.g., graphene.
Low σ values mean that E a is lowered with respect to its intuitively expected proportionality with the E bind . The physical meanings of σ still remains to be explored, but we anticipate that this term is related to the degree of the metal-support interaction. Figure 3c shows the performance of various ML algorithms. The LASSO, elastic net, and ridge algorithms have RMSE values of 0.198, 0.198, and 0.264 eV, respectively. R 2 values do not vary much across the models. A comparison between the LASSO algorithm, which has the lowest RMSE, and our diffusion scalinglaw for single-atom catalysts (DSL-SAC) is presented in Fig. 3d. Both models predict satisfactorily the diffusion barrier E a with the majority of points scattered around the parity line. Therefore, we suggest the simple form of DSL-SAC to estimate diffusion behaviors of isolated metal atoms on a support. Overall, this ML analysis validates statistically the intrinsic correlation between E a and (E bind ) 2 /E c for SACs. Additionally, some efforts were made to estimate binding energies of metal SAs on supports using several physical descriptors 17,27 . Therefore, we expect that the prediction of binding energies (E bind ) can further accelerate the screening of diffusion barriers (E a ) of metal SAs on supports and save computational expenses.
Validation of DSL-SAC model To further validate our DSL-SAC model, we determined the stability of isolated Pd and Pt atoms on the (100) surface of LaFeO 3 perovskite. Given its complex magnetic properties and the Jahn-Teller effect, it is computationally expensive to optimize the geometry of a given starting configuration of LaFeO 3 (see Supplementary Table 7) Table 8)

Stability assessment
To assess the stability of single-metal atoms on a support, we estimate the characteristic time of their diffusion (τ diffusion ) using the following equations: where k diffusion is the rate constant, k B is Boltzmann's constant, T is the temperature, h is Planck's constant, and R is the gas constant. In this study, the discussed lifetime is the characteristic time of diffusion of metal adatoms. The actual lifetime of a single atom on a support strongly depends on this characteristic diffusion time but also on the metal loading. In the Supplementary Table 9, we estimate the SAC lifetime for 0.5 wt% Pt/CeO 2 (111), demonstrating that this time is typically 1-2 orders of magnitude higher than τ diffusion . Given the strong variation of the τ diffusion with temperature, we employ τ diffusion as a conservative metric to compare the lifetime of SACs. It is worthy of note that the high concentration of metal SAs on support surfaces may influence the lifetime prediction, especially for weakly bonded systems in which metal SAs easily diffuse and aggregate into particles via Ostwald ripening. Figure 4 shows the estimated SAC lifetime as a function of E bind and E c at room temperature and a typical high temperature of 1073 K at which SACs may be exposed to for several hours. This figure predicts the minimum required binding energies of single atoms (examples given for Pt, Pd, and Ni) on a support to obtain stable SACs for a day at the two indicated temperatures. In this study, we focused on idealized surfaces without defects (except for steps), hydroxyl groups or dopants, which may act as nucleation sites for sintering and in this way alter diffusion processes of metal adatoms [35][36][37] . Despite this simplification, the framework introduced is general and the correlations provide a guiding tool for initial materials selection. To put this into practice, we discuss the work of Datye and colleagues, who prepared Pt/ CeO 2 SACs by a vapor-phase synthesis at 1073 K for 12 h 1,38 . At this temperature, the binding energy to avoid Pt agglomeration is estimated to be 6.3 eV. We note that the highest binding energies on CeO 2 are provided by CeO 2 (100) (5.47 eV) and steps of CeO 2 (111) (5.54 eV), the binding energy on the most stable CeO 2 (111) surface being much lower (see Supplementary Fig. 7). Thus, we predict that the single Pt atoms would agglomerate under the given conditions. However, we note that the synthesis is carried out in oxygen, i.e., Pt is present as PtO 2 39 . The computed binding energy for a PtO 2 on a step of CeO 2 (111) is 7.5 eV, explaining the high stability of the SAC in the recent experimental work of Datye 40 . Recently, Lopez and co-workers found that Pt single atoms can be trapped on CeO 2 (100) in the form of Pt 2+ with four O ligands owing to the inherent surface oxygen mobility 21 . Dynamic charge transfer can evidently affect the binding strength and diffusion behavior of single Pt atoms on CeO 2 (100) 21 . Doping of a metal like Pt in the surface of ceria also leads to a much higher binding energy (>10 eV), representing a very stable catalyst under most conditions (see Supplementary Fig. 8 17 . In the present work, we extended this thermodynamic approach by determining scaling relations for the activation barrier of diffusion of the single atom, which is the key kinetic step in the sintering process. Modern ML approaches identify a correlation which includes, in addition to the binding strength, also the cohesive energy of the metal, which reflects the intrinsic chemical reactivity of the single-metal atom. Previously, Mavrikakis and co-workers found that the diffusion barrier of adsorbates such as atomic O and N and molecular CO species depends in a linear fashion on the corresponding adsorption energies on transition metal surfaces 42 . Our work emphasizes a more complex correlation for the diffusion activation barrier of metal adatoms on typical supports. We expect that the diffusion scaling-law defined for single transition metal atoms on supports will find general application and can be improved by taking also into account support defects, support dopants, and hydroxyl group terminations present on typical metal oxides. The model is validated by accurate predictions of diffusion activation barriers of single Cu, Ru, Pd, and Pt atoms on LaFeO 3 (100). The computational framework was also used to discuss the stability of Pt atoms trapped on a stepped CeO 2 (111) surface, highlighting the role of defects in stabilization of this relevant case study under experimental conditions.
In summary, we started from DFT calculations and determined a meaningful correlation of the activation barriers for diffusion of single-metal atoms on a support with two easily accessible parameters, E bind and E c . Various ML methods were used to assist the physical descriptor discovery without assuming the functional form of the model explicitly. Contrast to many complex ML or black-box deep learning models, the developed diffusion scaling-law (DSL)-SAC consisting of a single descriptor (E bind ) 2 /E c offers an interpretable and generalizable model providing a facile approach to screen large numbers of metal-support combinations. Our approach provides a step toward understanding the stability of SACs, by properly considering the activation barrier for diffusion rather than the simple thermodynamic metric of binding energies.
Our study also provides a powerful strategy to rationally design SACs with promising stability.

DFT calculations
Spin-polarized DFT calculations were carried out by the Vienna ab initio simulation package 43 . The ion-electron interactions are represented by the projector-augmented wave method and the electron exchange-correlation by the GGA with the PBE exchange-correlation functional 44 . For ceria systems, the DFT + U approach was used with U The climbing image nudged-elastic band algorithm was used to identify the transition states of the metal atom diffusion on a support 45,46 . The total energy difference was less than 10 −4 eV and the relaxation convergence criterion was set at 0.05 eV/Å.

Regularization and feature-selection algorithms
We consider the usual linear regression model with n observations and p descriptors: 31,32,47 x 1 , …, x p , where x j = (x 1j , … x nj ) T for j = 1,…, p. The response y = (y 1 , …, y n ) T is predicted by A model training procedure produces a vector of coefficientŝ β ¼β 0 ; ;β p by minimizing a loss function. For the ordinary least squares (OLS), the loss function is the residual sum of squares. Standardization of a dataset which rescales all descriptor values to a centered mean of 0 and have variance of 1 is a requirement for many ML algorithms. It ensures all descriptors vary within the same order of magnitude so that the selection algorithms can identify the dominate descriptors correctly. For descriptor p i , the scaled values z i are where μ i is the mean of training samples and σ i is the standard deviation.
The standardized form of model can be written as The coefficient weightsβ Ã ¼β Ã 0 ; ;β Ã p in the standardized form can be obtained from the original coefficients by, One should note that onlyβ Ã reflects the weightage of each descriptor in producing the responses. Here we introduce three regularization techniques including the elastic net, the LASSO, and ridge regression with the aim to select significant descriptors and to prevent overfitting in OLS. The coefficient weightsβ Ã are obtained through minimizing the loss function containing the residual sum of squares, L2 and L1 norm of the coefficients: The above equation is the general formula of the elastic net, which contains both L1 and L2 penalty functions. Rewriting Eq. (11) gives, The LASSO and ridge regression are special cases of the elastic net, where λ 1 = λ, λ 2 = 0, r 1 = 1 (L1 penalty only) or λ 1 ¼ 0; λ 2 ¼ λ; r 1 ¼ 0 (L2 penalty only), respectively (see Supplementary Fig. 9). The training of ML models requires the tuning of hyperparameters including the degree of shrinkage λ and the amount of L1 penalty r 1 .
Overfitting is prevented by shrinking the number of descriptors included in the model (L1) or the magnitude of the coefficients (L2). Removing descriptors from the model can be seen as setting their coefficients to zero. Ridge regression only shrinks the magnitude of the coefficients but keeps all descriptors in the model. The LASSO shrinks the number of descriptors and the nonzero coefficients indicate the significance of corresponding descriptors after selection. The elastic net performs both L1 and L2 regularization which shrinks both the number of descriptors included and the magnitude of remaining descriptors. The amount of L1 penalty r 1 determines the degree of shrinkage on the number of descriptors.
The primary descriptors are E c and E bind . We also included terms of the order −2, −1, −0.5, 0.5, and 2 as well as the natural logarithm of the primary descriptors (E c and E bind ) as secondary descriptors. As tertiary descriptors, we added pair interactions between primary and secondary descriptors (see Supplementary Table 2). We perform the training in the Scikit-learn ML package in Python 29 . The dataset is first randomly shuttled and split into the training (80% of the dataset) and testing set (20% of the dataset). We then performed tenfold cross validation on the training set. The loss function was evaluated by leaving 10% of the training set out and using the rest to fit the coefficients. This process was repeated for ten times. The best set of coefficients giving the minimal value of the loss function during the entire training procedure were selected. The final prediction errors were determined by the RMSE of the testing set.
Genetic programming (GP) GP with symbolic regression is a supervised learning technique to identify an underlying mathematic expression that best describes the relationship between input and output data. The search space consists of combinations of simple mathematical operators on the input descriptors. An evolutionary algorithm is used to evolve a population of randomly generated candidate models according to natural-selection rules (selection, crossover, and mutation). Each model is associated with a fitness value, which in our case is the RMSE value of diffusion barrier E a . The advantage of GP is that no manual combination of descriptors is needed. The population of models would evolve based on the rules towards the optimal, which can be seen as a stochastic optimization process. We include addition, subtraction, multiplication, division, square root, and natural logarithm as operators for two descriptors, E bind and E c . We perform GP analysis in the Python package gplearn 6 . The same set of training data (consisting 80% of the dataset) was used for training GP models and the rest 20% of the dataset was used as the testing set. We performed five parallel runs initialized with different random seeds (see Supplementary Fig. 10). In each run, the population size was 5000. The fittest models were selected after a sufficient number of generations (here we use 100) when the evolution process converged to one solution.