Understanding High-Temperature Chemical Reactions on Metal Surfaces: A Case Study on Equilibrium Concentration and Diffusivity of CxHy on a Cu(111) Surface

Chemical reactions on metal surfaces are important in various processes such as heterogeneous catalysis and nanostructure growth. At moderate or lower temperatures, these reactions generally follow the minimum energy path, and temperature effects can be reasonably described by a harmonic oscillator model. At a high temperature approaching the melting point of the substrate, general behaviors of surface reactions remain elusive. In this study, by taking hydrocarbon species adsorbed on Cu(111) as a model system and performing extensive molecular dynamics simulations powered by machine learning potentials, we identify several important high-temperature effects, including local chemical environment, surface atom mobility, and substrate thermal expansion. They affect different aspects of a high-temperature surface reaction in different ways. These results deepen our understanding of high-temperature reactions.


More computational details 1.1 Vibrational contributions on energy barriers
Three vibrational terms on energy barrier are zero-point-energy (ZPE), vibrational enthalpy, and vibrational entropy contributions: where is the vibrational frequency of the i-th mode, h is the Planck's constant, is the Boltzmann constant. TS and IS denote the transition and initial state, respectively.

Widom insertion to calculate the excess chemical potential
Since all hydrocarbon species considered in this study have a favored adsorption site on the top two Cu atom layers (surface hollow sites for H/CH3, surface bridge sites for C2, and subsurface interstitial site for C), we choose an insertion space containing surface/subsurface Cu atoms. In the vertical direction, we set the upper limit of the space (the upper limit of the lowest atom for multi-atom species) to be the average height of the top Cu atom layer plus one and half of the Cu-Cu layer spacing, the lower limit (the lower limit of the highest atom for multi-atom species) to be the average height of the second Cu layer minus half of the Cu-Cu layer spacing.
If we define the concentration of a species in this insertion space as the number ratio between species and Cu atoms in the thin volume, since the number of surface/subsurface Cu atoms equals the number of surface hollow sites, the 3D species S2 concentration in the insertion space is identical to the 2D species coverage for H/CH3. For other species, such a mapping can be obtained by simply applying a constant factor.
Given a specific species and temperature, we do the following steps: 1. Run GAP-MD for 110,000 steps at the given temperature with a pure-Cu surface model.

2.
Abandon the first 10,000 steps and extract one structure from every 10 steps in the remaining 100,000 steps.
3. For each extracted substrate structure, randomly insert the given species into the insertion space，calculate and record the energy change due to insertion. 4. Repeat step 3 for 100 times. Finally, we have 1,000,000 samples to calculate the average value of exp (− / ). The excess chemical potential can thus be acquired according to its definition = − × 〈exp (− / )〉 For CH/CH3/C2, we fixed the C-H/C-C bond length to the optimized value in the adsorption state. The energy difference caused by the bond-length constraints can be estimated via MD simulations of the system with and without bond length constraints at each temperature. It turns out to be small (about 0.14 eV for CH at 1000 K) and it is added as an energy correction to the excess chemical potential.
For monoatomic species (C and H), the insertion is performed by simply giving three evenly distributed random numbers to determine the insertion position. For diatomic species (CH and C2), the position of one C atom is randomly given. After that, a random direction is given by providing three random numbers according to the normal distribution. The position of the other atom is determined by the bond length and the random direction. For CH3, a random direction is given to determine the rotational symmetry axis and a random number is given to determine the rotation angle of H atoms around the symmetry axis.
The ideal-gas chemical potential at a reference density is defined as follows. For monomer, it contains only the translational part: where is the density of substrate atoms which slightly changed due to the volume expansion at a higher temperature, and Λ is the thermal de Broglie wavelength: For multi-atom species, it additionally contains the vibrational and rotational parts: where vi is the i-th vibrational frequency of a multi-atom species, I is the degenerate moment of inertia for a linear species, IA/B/C are the three principal moments of inertia for a non-linear species and is the symmetry number.

Kirkwood coupling parameter method
For C and C2, simulations based on the Kirkwood coupling parameter method were also performed. This method is based on a series of artificial systems, in which the interaction between the adsorbate and the Cu environment is scaled with a factor α. We first calculate the excess chemical potential for the α=0.1 system using Widom insertion: Then, the excess chemical potential for the real system α=1 is calculated with thermodynamic integration: where 〈 〉 is the ensemble average of for the artificial system with scaled factor α.
To get 〈 〉 , we performed a series of MD simulations for different α ranging from 0.1 to 1 (the α=1 system is the real system). To get the scaled interaction, we calculate forces for both the complete model ( ) and the model with the adsorbate removed ( ). Then, the artificial forces are given by = + × ( − ) For atoms in the adsorbate, is zero. In this way, the interaction between Cu is kept but that between Cu and adsorbate is scaled with the factor α. Meanwhile, the energy difference between two systems with and without the adsorbate ( ) is recorded in each MD step to get its ensemble average 〈 〉 .
This artificial-force scheme for MD simulation can be directly used for the C monomer. For C2, it scales also the C-C interaction which is incorrect. Here, we simply treat C2 as a rigid body by fixing the C-C bond length in all MD simulations. The excess chemical potential obtained is finally corrected by applying the bond constrain correction at each temperature.

GAP training
For each hydrocarbon species, we first run ab initio MD at 3000 K for 50 steps. All these frames are used as the first batch of the training set, from which an initial GAP is obtained. Then, GAP-MD simulations are performed to obtain other configurations of the training set. To extract N structures at a temperature between 2000 and 600 K, we run a GAP-MD simulation at this temperature for 1000+100*N steps. We uniformly extract one structure in every 100 frames for the last 100*N steps. At temperatures below 600 K, the simulation steps are increased to 1000+1000*N steps to extract a structure from every 1000 frames. DFT calculations for all these extracted structures are performed with VASP to get their energies and forces. The value N for each temperature is given in Figure S1a.
To train the less-accurate initial GAPs used for generating more structures in the training set, we use only SOAP descriptor 1 with an expansion of the neighbor density up to lmax=6, nmax=10, and cutoff=5 Å. For the final GAPs used for calculating diffusivity and coverage properties, we use distance_2b descriptor with cutoff=4 Å, angle_3b descriptor with cutoff=3 Å, and Soap descriptor with lmax=8, nmax=12, and cutoff=5.5 Å. The parameter atom_sigma is set to be 0.3 and 0.5 for systems with and without H involved, respectively. Figure S1. (a) The iterative training schemes. "Frames" in the table gives the number of structures extracted as a training set at each temperature. More structures were collected at temperatures close to the experimental temperature of graphene CVD growth. (b) Energy distribution in each batch of the training set, and in the whole training set (denoted as "All"). An ab initio MD trajectory at 1300 K denoted as "AIMD" with the same number of frames to the training set is presented as a comparison. (c) Principal component analysis on structures of the whole training set and that of the 1300 K AIMD trajectory. SOAP descriptor 1 is adopted to convert the structure into a vector for the principal component analysis. Figure S2. Energy and force error analysis for Cu substrate with adsorbate C/C2/H/CH/CH3 or without adsorbate (denoted as "None"). The test set obtained from an AIMD simulation at 1200 K is in red and that at 1400 K is in blue. All these results indicate high accuracy with very low root-mean-square errors (RMSE).   Table S1. Vibrational correction terms and the final free energy barriers for C/C2/H/CH/CH3 under different temperatures. T is the temperature in K. ΔZPE, ΔCp and Δ(-TS) are the difference of ZPE, vibrational enthalpy, and vibrational entropy part between transition and initial states. ΔGvib is the summation of these three terms. "Barrier" is calculated as the summation of the NEB barrier and ΔGvib. The substrate, except for the top layer, is fixed to calculate the vibrational frequencies.    To obtain converged diffusivity, it is suggested that the total mean square displacement (TMSD) in the simulation should be up to a few hundred Å 2 . 2 Such a criterion is fulfilled in all our simulations for diffusivity ( Figure S10), indicating the trajectories are long enough. For Widom insertion, the typical number of insertion attempts is 10 5~1 0 6 . 3,4 In this study, Widom insertion simulations are performed with 10 6 attempts each. At the same time, we divide each simulation into five parts with equal steps and calculate the diffusivity or excess chemical potential based on each part. The obtained standard deviations are drawn as the error bars in Figure 2 and Figure S11, S10 which are reasonably small. Figure S10. TMSD of adsorbates in MD simulations at different temperatures. Figure S11. The excess chemical potential of C/C2/H/CH/CH3 adsorbates on Cu surface at different temperatures. The error bar gives the standard deviation of five independent sampling groups. Figure S12. The adsorption energy of C monomer on expanded Cu substrate with local and global expansion models. In the global expansion model, the entire surface slab is expanded in the x/y/z directions. All Cu atoms are fixed and only the C monomer is optimized to get the interaction energy.