Stabilizing mechanism of single-atom catalysts on a defective carbon surface

Single-atom (SA) catalysts represent the ultimate limit of atom use efficiency for catalysis. Promising experimental progress in synthesizing SA catalysts aside, the atomic-scale transformation mechanism from metal nanoparticles (NPs) to metal SAs and the stabilization mechanism of SA catalysts at high temperature remain elusive. Through systematic molecular dynamics simulations, for the first time, we reveal the atomic-scale mechanisms associated with the transformation of a metal NP into an array of stable SAs on a defective carbon surface at a high temperature, using Au as a model material. Simulations reveal the pivotal role of defects in the carbon surface in trapping and stabilizing the Au-SAs at high temperatures, which well explain previous experimental observations. Furthermore, reactive simulations demonstrate that the thermally stable Au-SAs exhibit much better catalyst activity than Au-NPs for the methane oxidation at high temperatures, in which the substantially reduced energy barriers for oxidation reaction steps are the key. Findings in this study offer mechanistic and quantitative guidance for material selection and optimal synthesis conditions to stabilize metal SA catalysts at high temperatures.


INTRODUCTION
Single-atom (SA) catalysts contain isolated single-metal atoms dispersed on a support surface, and represent the ultimate limit of atom use efficiency for catalysis 1,2 . Other desirable catalytic properties of SA catalysts include high selectivity, tunable high activity, and unique metal coordination environments [3][4][5][6] . As a result, SA catalysts have attracted tremendous attention as a new frontier in heterogeneous catalysis for many reactions, including oxidation 7,8 , hydrogenation 9,10 , electrocatalysis [11][12][13][14] , and so on 15,16 . However, SAs have the tendency to aggregate into particles at an elevated temperature due to their excess surface-free energy, which is known as thermal deactivation or sintering [17][18][19] . Thus, recent studies have sought to improve the thermal stability of SAs by forming strong metal-support interactions at 800-1200°C on different supports (e.g., TiO 2 10,20 , CeO 2 7 , FeO x 21 , nitrogen-doped carbon 22 , and carbon nanofibers 23 ). In some of these experiments, a counterintuitive transformation from nanocluster to SAs is observed at a high temperature 22 . However, it is challenging to uncover the underlying mechanism of this reverse dispersion solely using experiments, given that SAs are formed at atomic scale, and the transformation occurs in a period of milliseconds 23 . Therefore, encouraging experimental evidence of high-temperature-assisted conversion of metal nanoparticles (NPs) into SAs aside, the underlying SA formation mechanisms remain elusive. To this end, it is desirable to decipher such underlying mechanisms through comprehensive atomistic simulations of the NP-to-SA process to capture key atomicscale events that occur at a short period of time.
Recently, carbon nanomaterials have received enormous interest in the field of electrochemistry as a highly promising catalytic support because of its many advantages, including lightweight, low cost, adjustable porosity, high chemical and thermal stability, and controllable chemical properties by heteroatom doping [24][25][26][27][28] . In particular, carbon nanomaterial surface has shown high capacity for stabilizing metal SAs, resulting from the existence of defects in the surface 22 .
Motivated by recent experiments 22,23 , here we carry out comprehensive molecular dynamics (MD) simulations to decipher the underlying mechanism of the dispersion of metal NPs into metal SAs on a carbon surface under high-temperature shock heating, using Au as a model material. We reveal that defects in the carbon surface (e.g., various types of vacancies) play a pivotal role in dispersing Au-NPs and further stabilizing Au-SAs. We find that thermodynamically stable metal SA formation is governed by the interplay between the density/type of defects on the carbon surface, the cohesion energy between metal atoms, and the metal-carbon-binding energy. To demonstrate the high catalytic efficiency of metal SA catalysts, we perform reactive simulations to show that Au-SA catalysts can enhance the methane oxidation efficiency at 2000 K up to 94%, nearly four times higher than that of Au-NPs, a highly desirable performance in fuel and propellant combustion. The mechanistic findings from this study offer fundamental insight into facile and robust approaches to synthesize thermodynamically stable metal SA catalysts.

RESULTS AND DISCUSSION
Decipher the effect of defects in carbon surface on transformation of a Au-NP into Au-SAs We carry out MD simulations (see the "Methods" section for details) to reveal the effect of defects in a carbon surface on the transformation of a Au-NP into Au-SAs at high temperature. To this end, we compare the transformation of a Au-NP containing 123 Au atoms (~1.5 nm in diameter, referred to as Au 123 hereinafter) on three different types of carbon surface: defect free (i.e., a 10-by-10-nm pristine carbon surface), low defect density (186 missing carbon atoms in a 10-by-10-nm carbon surface), and high defect density (370 missing carbon atoms in a 10-by-10-nm carbon surface), as shown in the first panel in Fig. 1b, c. After initial equilibration at room temperature (300 K), the temperature is rapidly increased to and maintained at 1500 K, and then the system gradually cools back down to room temperature, mimicking the thermal loading in experiments 22 . Figure 1a depicts a representative temperature profile of the simulation system. The simulations capture the process of the evolution of the Au-NP on the three types of carbon surface at high temperatures, and the formation of Au-SAs on the carbon surface (see Supplementary Videos 1-3). In all three scenarios, the Au-NP is shown to experience Brownian-like random walk on the carbon surface, owing to thermal fluctuations at 1500 K. On the defect-free carbon surface, some Au atoms separate from the Au-NP during the random walk of the NP. These isolated Au atoms are not stable and keep moving disorderly on the defect-free carbon surface. Some isolated Au atoms can move close to and collide intensively with the remaining of the Au-NP, and merge back to be part of the NP again. Some isolated Au atoms may also form small clusters (e.g., three-atom bundle) on the defect-free carbon surface. At the end of the simulation, the remaining Au-NP still contains 90 Au atoms, and there are 30 Au-SAs and 1 tri-Au-atom bundle on the defect-free carbon surface at room temperature (Fig. 1b). Note that these Au-SAs are highly mobile and unstable, and thus could merge with Au-NP or other Au atom clusters.
By contrast, the random walking of the Au-NP on a defective carbon surface keeps leaving Au-SAs along its trace (Fig. 1c, d).
There are two features of such Au-SAs distinct from those on a defect-free carbon surface: (1) Au-SAs separated from the Au-NP remain bonded to and thus become immobile at the defect sites on the carbon surface, and (2) even the remaining Au-NP later collides with or passes through these Au-SAs; the Au-SAs do not merge back into the Au-NP. As a result, the Au-NP size keeps decreasing, while the Au-SAs are well dispersed on the defective carbon surface. The higher the defect density in the carbon surface, the better efficiency of Au-SA dispersion on the carbon surface. At the end of the simulations, the remaining Au-NP on a carbon surface with less defects (Fig. 1c) has only 23 atoms left, with the other 100 Au atoms well dispersed at the vacancy sites on the carbon surface; on a carbon surface with sufficient defects (Fig. 1d), all Au atoms in the NP can be dispersed and remain immobile at the vacancy sites during the high-temperature stage and after cooling down to room temperature. In other words, an Au-NP can be completely converted into stable Au-SAs on a carbon surface with sufficient defects.
Mechanism of Au-NP to Au-SA conversion The evolution of Au-NP on a carbon surface as described in Fig. 1 shows that the formation of Au-SAs on a defective carbon surface should satisfy three conditions: (1) Au-NP can move randomly on the carbon surface at high temperature; (2) Au-SAs can separate from the Au-NP and then be trapped by defect sites; (3) there are sufficient defect sites on the carbon surface. We use molecular statistics method and nudged elastic band (NEB) calculations to Fig. 1 Characterization of Au-NP and Au-SA formed on a defective carbon surface. a The system temperature change during the evolution process of Au-NP on a carbon surface. b-d Snapshots of the evolution of Au-NP (red) on a defect-free carbon surface (gray), a carbon surface with less defects, and a carbon surface with more defects, respectively.
investigate the underlying mechanism of the NP-to-SA conversion (see the "Methods" section for details). Figure 2 plots the energy variation of the system during the process of a Au-SA detaching from a Au-NP (with 123 Au atoms) and then attaching to a defect-free carbon surface (Fig. 2a) and on a defective carbon surface (Fig. 2b). As shown in Fig. 2a, as the Au-SA starts to detach from the Au-NP, the energy increases and reaches a peak value of 1.71 eV when the Au-SA is fully detached. Then as the detached Au-SA approaches and lands on the underlying defect-free carbon surface, the energy decreases modestly due to the formation of weak Van der Waals interaction between the Au-SA and the carbon surface. The corresponding energy in this final state is 1.19 eV higher than that of the initial state of Au 123 on the defect-free carbon surface. Therefore, the final state is thermodynamically unstable.
By contrast, as shown in Fig. 2b, as the Au-SA starts to detach from the Au-NP, the energy increases and reaches a peak value of 1.71 eV when the Au-SA is fully detached, the same as that in Fig. 2a. Then as the detached Au-SA approaches and lands on the underlying defective carbon surface, energy decreases drastically, which is attributed to the formation of Au-C chemical bond. As a result, the corresponding energy in this final state is 3.54 eV lower than that of the initial state of Au 123 on the defective carbon surface. Therefore, the final state is thermodynamically stable. Note that in Fig. 2b, we use a carbon surface with a bivacancy as a demonstration of the role of defects on stabilizing Au-SA. As to be further elaborated in the next section, such a stabilizing effect is universal for 18 different Au-SA/vacancy-defect bonding scenarios on a carbon surface (Fig. 4).
To further understand the continuous formation of Au-SAs from a Au-NP on a defective carbon surface as shown in Fig. 1c, d, we calculate the average cohesive energy per Au atom in a freestanding Au-NP (E cohensive ), which is plotted in Fig. 3 as a function of the number of Au atoms (n) in the NP (i.e., the Au-NP size, d). As the Au-NP size decreases, the average cohesive energy per Au atom decreases. Such a dependence can be fit by the following relation: where 3.91 eV represents the cohesive energy per Au atom in bulk Au, which agrees well with the experimental measurement of 3.93 eV 29 .
The trend in Fig. 3 can be understood by the increasing number of Au atoms on the NP surface as the NP size decreases, which in turn results in a lower average cohesive energy due to less saturated metallic bonds for the Au atoms on the surface. Therefore, once a Au-NP starts to dissociate and spread Au-SAs on a defective carbon surface, its size continuously decreases, leading to an ever-decreasing cohesive energy of the remaining Au atoms in the NP. As a result, the more Au-SAs dissociated and bonded to the defective carbon surface, the more readily the remaining Au-NP can be further dissociated as long as there are available defect sites on the carbon surface.
Thermal stability of Au-SAs on a defective carbon surface Defects in a carbon surface could have various sizes and types, as evident from the MD simulations shown in Fig. 1. Next, we demonstrate the thermal stability of Au-SAs on the most representative types of vacancy defects on a carbon surface. For comparison, we first calculate the binding energy of a Au-SA with a defect-free carbon surface (e.g., as shown in the last inset of Fig. 2a), which is 0.45 eV.    i The trace of a Au-SA on a defect-free carbon surface, and j the trace of a Au-SA trapped in various types of defect sites over the time at 3000 K. vacancy increases, it is possible that each defect site could trap more than one Au-SA, as evident from our MD simulation. For example, on a 3-V defect, it is possible to form three different Au-C bonding structures with 1, 2, or 3 Au-SAs trapped (Fig. 4c). We define the bonding energy per Au-SA in each bonding structure, E b , as the difference between the energy of the bonded structure and the sum of the energies of the carbon surface with the same defect and the freestanding Au-SAs to be trapped onto the defect site, normalized by the number of Au-SAs in such a bonding structure. The corresponding values of E b are listed below each bonding structure and further plotted in Fig. 4h in comparison with those of a Au-SA with a defect-free carbon surface. It is evident that the binding energy of Au-SA on a defect site is 6-18 times (i.e., 3.02-9.56 eV) higher than that on a defect-free carbon surface (0.45 eV). In other words, once a Au-SA is trapped on a defect site on the carbon surface, it is thermodynamically stable.
To verify that Au-SAs trapped on a defect site on the carbon surface are immobile and thus immune from merging and clustering at high temperatures, we further simulate the trace of a Au-SA on a defect-free carbon surface (Fig. 5i), and the trace of a Au-SA trapped in various types of defect sites on a carbon surface (Fig. 5j) at 3000 K. With a much lower binding energy on a defectfree carbon surface, the Au-SA travels randomly on the carbon surface driven by the thermal fluctuation. By contrast, the Au-SA trapped in the defect site remains nearly immobile even at such a high temperature (see Supplementary video 4). To further understand the trapping effect of defects on Au-SAs, we also compare the diffusion coefficients (see Supplementary Note 1 for calculation of diffusion coefficients) of a Au-SA on a defect-free and various defective carbon surfaces at different temperatures (from 500 to 3000 K) (Fig. S1). It is shown that the Au-SA on the defect-free carbon surface has a substantially increased mobility as the temperature increases (e.g., 42 times increase in the diffusion coefficient from 500 to 3000 K). By contrast, the Au-SA trapped in the defect site has a significantly lower mobility even at high temperature. For example, the diffusion coefficients of the Au-SA on various defect sites at 3000 K are about 10 times lower than those of the Au-SA on a defect-free carbon surface.
The above studies offer mechanistic understanding of the dispersion of Au-NPs into Au-SAs on a defective carbon surface, and the remarkable stability of such Au-SAs against high temperature due to the strong chemical bonds between the Au-SAs and the surrounding carbon atoms.
Catalytic reactivity: Au-SAs versus Au-NPs Methane-fueled liquid rocket engines are being developed (e.g., the Raptor by SpaceX) as a possible replacement for liquid hydrogen rocket engines, given the better stability of methane over long periods of time and a higher density over liquid hydrogen 30 . Therefore, how to enhance the oxidation efficiency of methane becomes crucial for improving the combustion performance of methane-fueled liquid rocket engines. Nanoparticulate gold is shown to be catalytically active for oxidation reactions.
Here we show that the thermally stable Au-SAs can significantly enhance the catalytic reactivity in the oxidation of methane in comparison with Au-NPs.
To investigate the catalytic efficiency of Au-SA, we simulate the methane oxidation in a system with 50 CH 4 molecules and 100 O 2 molecules catalyzed by Au-SAs and Au-NP, respectively, at 2000 K. During the methane oxidation process, CH 4 and O 2 first dissociate and then react to form many types of intermediate radicals and molecules (e.g., CH 3 , CH 2 O, CH 2 , CH, CH 3 O, CO, H 2 O, and OH) (Fig.  S3). Intermediate radicals (e.g., CH 3 , CH 2 O, CH 2 , CH, CH 3 O, and OH) can further dissociate and form oxidation products H 2 O and CO (Fig. 5a, b). The evolution of methane oxidation catalyzed by Au-NP and Au-SA is shown in Supplementary Videos 5 and 6. We define the methane conversion rate as the ratio of the amount of decomposed methane after the reaction to the total amount of methane before the reaction. As shown in Fig. 5c, the consumption of CH 4 catalyzed by Au-SA is much faster than that catalyzed by Au-NP. As shown in Fig. 5d, the methane conversion ratio at 1500-ps reaction under the catalyst of Au-NP is 24%, while the methane conversion ratio at 1500-ps reaction under the catalyst of Au-SAs reaches a value of 94%. The catalytic efficiency of Au-SAs is about four times higher than that of Au-NP, which can be attributed to the fact that only the Au atoms on the surface of the Au-NP participate in the catalyzation, while every Au-SA can contribute to the catalyzation process.
To further reveal the enhancing effect of Au-SA catalyst on methane oxidation quantitatively, we calculated and compared the energy barriers associated with different oxidation steps of the oxidation reaction for two scenarios: with Au-SA catalyst and without Au-SA catalyst, using NEB method. By tracing the evolution path of the elements belonging to CH 4 during the whole methane oxidation process, two major catalytic mechanisms concerning interatomic exchange are proposed as in Eqs. (2) and (3).
To compare the associated energy barriers in each step between two scenarios: without catalyst and with Au as catalyst for the two catalytic mechanisms in Eqs. (2) and (3), two catalytic cycles, each of which contains four elementary steps, are schematically shown in Fig. 5e, f. For the reaction path in Eq. (2), C-H bonds successively cleave at the Au surface, leading to the formation of C atom, which fills in the vacancy site on carbon surface. The dehydrogenation of methane requires to overcome four energy barriers (corresponding to the four steps in Eq. (2)) of 0.75, 2.46, 0.19, and 0.24 eV, respectively (3.64 eV in total). Without the catalyst, the corresponding four energy barriers for each reaction step are significantly higher, i.e., 4.39, 6.12, 4.33, and 4.44 eV, respectively (19.28 eV in total). Similarly, for the reaction path shown in Eq. (3) whose final products are CO 2 and H 2 O, without the catalyst, the associated energy barriers in each step (Fig. 5f) are 8.29, 2.91, 4.32, and 5.16 eV, respectively (20.68 eV in total). By contrast, with Au as catalyst, the corresponding four energy barriers decrease significantly to 0.93, 0.23, 2.87, and 0.82 eV, respectively (4.85 eV in total). Much lower activation energy for the methane oxidation prompts the occurrence of a reaction under the catalyst of Au atom.
The above comparison clearly epitomizes the significant effect of Au-SA catalyst on enhancing both the conversion and energy efficiency of methane oxidation (due to a drastically lower overall energy barrier).
In sum, for the first time, we reveal the atomic-scale mechanisms associated with the transformation of a Au-NP into an array of stable Au-SAs on a defective carbon surface at a high temperature, through systematic MD simulations. Such atomicscale transformation mechanisms feature a low-energy barrier for a Au-SA to decompose from a Au-NP and a high Au-C-binding energy to stabilize the resulting Au-SAs. This transformation mechanism sheds crucial light on understanding the experimental observation of stable Au-SA formation on a defective carbon nanofiber. Modeling study also reveals the pivotal role of defects in carbon surface in trapping and stabilizing the Au-SAs. Furthermore, reactive simulations demonstrate that the thermally stable Au-SAs exhibit much better catalyst activity than Au-NPs for the methane oxidation at high temperatures, in which the substantially reduced energy barriers for oxidation reaction steps are the key.

METHODS
The MD simulations described in the main text are performed using ReaxFF potential, as implemented in the LAMMPS 31 . The ReaxFF reactive force field, based on a bond-order methodology 32 , provides a good compromise between accuracy and computational efficiency. The ReaxFF parameters used in this work are optimized using density functional theory (DFT) calculations with the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional and norm-conserving pseudopotentials 33 . The Au/C/S/H force field parameter 33 is used for the pure Au/C system to simulate the evolution of Au-NP on carbon surface. This potential has demonstrated the capability of accurately representing the interaction between carbon and gold 34,35 . To study the accuracy of the Au/C/S/H force field parameters, the binding energy of Au atom with various carbon vacancies shown in Fig. 4a is calculated using DFT. The comparison between ReaxFF simulations and DFT calculations in Fig. S4 shows that the Au/C/S/H ReaxFF potential reproduces the same qualitative trends in binding energies as what is obtained with DFT. To compare the catalytic efficiency of Au-SAs and Au-NPs on methane oxidation, a merged ReaxFF force field is used with Au/O/ C/H parameters (see Tables S1-S6 for ReaxFF parameters in tabular form, and Supplementary Note 2 for ReaxFF parameters in ReaxFF format) obtained by combining the Au/O parameters from Joshi et al. 36 with Au/C/ S/H parameters from Jarvi et al. 33 , and with C/O/H/S parameters from Rahaman et al. 37 . To investigate the accuracy of the Au/O/C/H force field parameters, we further calculate the activation energy in each step for the catalytic mechanisms in Eqs. (2) and (3) using DFT as compared in Table S7. It shows that the agreement between the activation energies obtained from ReaxFF simulations and DFT calculation is reasonably good for all steps in both catalytic mechanisms. In this paper, all the microstructures and deformation mechanisms are analyzed in the software OVITO 38 .

Au-NP evolution on carbon surface
To investigate the transformation of Au-NP to Au-SA, geometry minimization of the Au/carbon surface system is performed prior to MD simulations. At the beginning of the simulation, we generate an ensemble of velocities using a random number generator at 300 K, and equilibrate the system under constant energy and fixed volume for 2.5 ps, with a time step of 0.25 fs in the microcanonical ensemble (NVE). Meanwhile, we reset the temperature to 300 K of the system by explicitly rescaling its velocities every step. During the equilibrating process, the temperature fluctuation of system is effectively controlled using this method. The temperature is then increased to 1500 K rapidly in 1.5 fs, and the system is maintained at 1500 K for 622.5 ps in the canonical ensemble (NVT) and Nosé-Hoover thermostat. After that, we decrease the temperature to 300 K slowly in 600 ps. In the simulation, the carbon surface is 10 × 10 nm in size, and the defects are created by removing carbon atoms randomly. The carbon surface with more defects (370 vacancies), less defects (186 vacancies), and no defects is considered. The initial Au-NP contains 123 Au atoms with 1.5 nm.

Methane oxidation assisted by Au catalyst at high temperature
We study the reactions of CH 4 /O 2 mixture and dissociation of CH 4 on Au-NP and Au-SA-based catalysts. The defective carbon surface with only a Au-NP and with only Au-SAs, as shown in Fig. 1d (I and V), is used as a catalyst for the methane oxidation, respectively. The size of the box in every simulation is 10 × 10 × 6 nm with 50 CH 4 and 100 O 2 gas-phase molecules, and the periodic boundary condition is implemented in all three directions. Before MD simulations, every system is first energy minimized via a conjugate gradient algorithm to eliminate simulation artifacts that can arise from initial high-energy contacts. A time step of 0.25 fs is selected to ensure energy conservation at the temperatures used in this study. The reaction temperature is initially set at 300 K, and is increased at a rate of 10 K/ps, yielding a final temperature of 2000 K after 1 × 10 6 iterations.

NEB calculation
To characterize in detail the observed mechanism from an energetic point of view, and to give an estimate of the activation energy barriers, we reproduce the dispersion pathway of Au-SA during the Au-NP evolution and the methane oxidation pathway by means of the NEB methodology and the climbing image (CI) method 39,40 implemented in the ReaxFF code and provided to us by Stijn Huigh and Erik Neyts 41 . For the Au-SA dispersion pathway, the well-relaxed system with a Au 123 NP on carbon surface is the initial structure; the final well-relaxed structure is one Au atom leaving the Au 123 NP and landing on carbon surface. CI-NEB calculations are performed to search six transition structures. According to the calculated energy of the six transition structures, the energy barrier for a Au-SA to dissociate from the Au 123 NP and to bond to the carbon surface can be obtained. For the methane oxidation pathway, we first trace the reaction pathway of carbon atom in dynamic simulations to see which products the atom ever formed. Then, we determine the final and initial structures to calculate the reaction energy barrier using NEB methodology.