Surface carbon species formation from ethylene decomposition on Pd(100): a first-principles-based kinetic Monte Carlo study

Abstract

Based on the activation barriers and reaction energies from periodic density functional calculations, we conducted kinetic Monte Carlo (kMC) simulations of surface carbon species formation from ethylene decomposition on a Pd(100) surface. A comprehensive reaction network of ethylene decomposition involving such intermediates as CH₂CH, CHCH, CH₂C, CHC, CC, CH₂ and CH was proposed. Our kMC simulations show that the most probable pathway of ethylene decomposition on Pd(100) is CH₂CH₂ → CH₂CH → CH₂C → CHC → CC, among which the dehydrogenation of CH₂CH₂ to CH₂CH is the rate-limiting step with the activation barrier of 1.51 eV, followed by CH₂CH₂ → CH₂CH → CHCH → CHC → CC, whose rate-limiting step is the dehydrogenation of CH₂CH to CHCH with the activation barrier of 1.59 eV. The two most probable pathways produce a carbon dimer as the final product, since the activation barrier of the C–C bond cleavage reaction is so high (2.32 eV) that it is almost impossible for it to occur before the metal surface is totally poisoned by surface carbon species. Another three feasible pathways are: (i) CH₂CH₂ → CH₂CH → CHCH → CH → C, (ii) CH₂CH₂ → CH₂CH → CHCH → CHC → CH + C → C and (iii) CH₂CH₂ → CH₂CH → CH₂C → CHC → CH + C → C, whose final products contain surface carbon monomers. And the reactions involving C–C bond cleavage are the rate-limiting step of the three pathways. Simple as the reaction network of ethylene decomposition looks, it is still difficult to analyze the decomposition mechanism merely according to the activation barriers from DFT calculations. Our work here demonstrates that kMC simulations can nicely tackle the problem on competitive reaction pathways, each of which involves some reactions with relatively low activation barriers (e.g. the dehydrogenation reactions involved in ethylene decomposition) and some other reactions with relatively high activation barriers (e.g. the C–C bond cleavage reactions involved in ethylene decomposition).

---

1. Introduction

Pd catalysts, especially when alloying with Au, are highly effective in promoting various reactions such as CO oxidation,^1,2 H₂ oxidation to H₂O₂,³ alcohol oxidation,⁴ and vinyl acetate synthesis.⁵ The surface carbon species formation over supported Pd catalysts was first reported in 1970,⁶ which can lead to catalytic deactivation in reactions involving hydrocarbons, such as ethylene. Ethylene can be transformed into such important industrial chemicals as vinyl acetate, vinyl chloride and vinyl alcohol. To better understand the surface carbon species formation from olefin hydrocarbons, ethylene as the simplest olefin deserves a comprehensive and thorough investigation. Ethylene decomposition to surface carbon species can be considered to be fundamental and important for understanding deactivation of transition metal catalysts.

The reaction network of ethylene decomposition is shown in Fig. 1, which involves eight dehydrogenation reactions and six C–C bond cleavage reactions. Chen et al.^7,8 conducted investigations on the transformation and total decomposition of ethylene on several different metal (111) surfaces and carried out DFT-based analysis on each reaction step. They found that for kinetic reasons, most species tend to undergo C–H bond scission rather than C–C bond cleavage. However, since each possible pathway of ethylene decomposition involves some dehydrogenation reactions with relatively low activation barriers and some C–C bond cleavage reactions with relatively high activation barriers, no explicit overall preferred pathway was proposed. Therefore, further efforts need to be made to pick out the preferred pathways from all the possible ones of ethylene decomposition. Although many efforts have been devoted to ethylene transformation on Pd(111) surface,^9–11 there is scarce work conducted on Pd(100) surface which is thought to be more active than Pd(111) surface as for the vinyl acetate synthesis from ethylene acetoxylation.^5,12,13 Thus, we carried out the work on Pd(100) surface.

Fig. 1 Reaction pathways of ethylene decomposition to surface carbon species on Pd(100) surface. The black arrows and the red arrows represent dehydrogenation reactions and C–C bond cleavage reactions, respectively. In the following text, Rx will stand for the corresponding reaction in this figure with the letter x near the arrow.

Kinetic Monte Carlo (kMC) simulation method is a stochastic approach to examine the events occurring on the catalyst surface.¹⁴ In principle, kMC allows us to keep track of the exact dynamic evolution of all surface species as a function of time and reaction conditions. Therefore, kMC can directly and explicitly exhibit the possible pathways of a complex reaction network. It has been widely applied in heterogeneous catalysis.^15–20 Thus, in this work, kMC simulation method was employed to explore the possible pathways of ethylene decomposition on Pd(100). The most likely pathway of ethylene decomposition will be given to provide theoretical reference and guidance for further application and development of efficient commercial catalysts, and the control of carbon deposit.

2. Computational methods

All the DFT calculations were conducted with the program package Dmol3 in Materials Studio,²¹ employing the generalized gradient approximation (GGA) with revised Perdew–Burke–Ernzerhof functional (RPBE).²² The detailed setups and results of DFT calculations are available in our previous work.²³ The intrinsic barriers obtained from the DFT calculations are shown in Table 1.

Table 1 Activation barrier (ΔE/eV mol⁻¹) and reaction energy (ΔH/eV mol⁻¹) of the elementary reaction step involved in ethylene decomposition on Pd(100)

Reaction	Pd(100)
	ΔE	ΔH
CH2CH2* → CH2CH* + H*	1.51	0.79
CH2CH* → CH2C* + H*	1.49	0.74
CH2CH* → CHCH* + H*	1.59	0.45
CH2C* → CHC* + H*	1.48	0.23
CHCH* → CHC* + H*	1.49	0.55
CHC* → CC* + H*	1.15	0.13
CH2* → CH* + H*	1.09	−1.06
CH* → C* + H*	1.45	−0.11
CH2CH2* → 2CH2*	2.63	1.07
CH2CH* → CH2* + CH*	2.13	0.15
CH2C* → CH2* + C*	2.35	−0.69
CHCH* → 2CH*	1.71	−0.43
CHC* → CH* + C*	1.75	−1.16
CC* → C* + C*	2.32	−1.36
2H* → H2(g)	0.75	0.75

---

The detailed DFT-based kMC approach has been reported in some previous publications.^15,17–20 The primary information can be provided as follows.

Variable step size method was adopted for kinetic Monte Carlo algorithm. In the kMC models, the Pd(100) surface is represented by a two-dimensional square periodic grid of 64 × 64 lattices containing three types of surface sites: top sites, bridge sites, and hollow sites. The preferred adsorption site of each species was determined by the DFT calculations. Periodic boundary conditions were employed to provide an adequate representation of the periodicity exhibited by the Pd(100) surface. In this way, surface species reappear at the opposite side of the lattice after jumping across the boundary. The initial state of present kMC simulations corresponded to an ethylene atmosphere with relatively low pressure which can continuously impinge on the Pd(100) surface with a reasonable rate. The rates of such elementary events as surface reactions and desorption are calculated through transition state theory (TST), where the rates can be calculated from eqn (1).

where v_i is the pre-exponential factor, R is the gas constant, T is the reaction temperature and ΔE_i is the activation barrier of the elementary event i. A theoretical value of 10¹³ s⁻¹ was chosen for the pre-exponential factors of all the surface reactions and desorption steps.²⁴ Since ethylene is the only source considered in the adsorption step on Pd(100) surface, the adsorption rate of ethylene is directly given with a moderate value that is not too large nor too small instead of calculating it from a given temperature and pressure. The temperature (450 K) of vinyl acetate synthesis from ethylene acetoxylation was adopted. Given the fact that the main effect of diffusion is to bring the adlayer to steady state or equilibrium, and the diffusion rate calculated from eqn (1) does this at a much shorter time scale than the time scale of all the other elementary events. Thus, in principle, we can reduce the diffusion rate to the extent that the adlayer is still brought to steady state or equilibrium without changing the kinetics. Several different diffusion rates were applied in different simulations.

3. Results and discussion

During the kMC simulations, the times each elementary event happens are recorded. The times each elementary event happens for different diffusion rates are listed in Table 2, and the final coverage of each species for different diffusion rates is shown in Table 3. Considering the stochastic nature of kMC simulation method, the convergence of the results is satisfactory. Every time the diffusion rate increases by one order of magnitude, the corresponding total simulation time will also increase by almost one order of magnitude. In order to balance the calculation efficiency and the validation of the computation results, in the following work, the diffusion rate of 0.1 s⁻¹ will be used for further calculations and the results will be used to analyze the reaction network of ethylene decomposition.

Table 2 The times each elementary event happens for different diffusion rates

Reaction	Diffusion rate (s^−1)
	10	1	0.1	0.01
CH2CH2(g) → CH2CH2*	1 118 025	1 058 203	897 993	689 189
CH2CH2* → CH2CH* + H*	2389	2423	2550	2684
CH2CH* → CH2C* + H*	1635	1628	1649	1697
CH2CH* → CHCH* + H*	126	128	141	115
CH2C* → CHC* + H*	1227	1215	1258	1357
CHCH* → CHC* + H*	80	82	102	84
CHC* → CC* + H*	1281	1268	1325	1401
CH2* → CH* + H*	0	0	0	0
CH* → C* + H*	18	29	12	8
CH2CH* → CH2* + CH*	0	0	0	0
CH2C* → CH2* + C	0	0	0	0
CHCH* → 2CH*	23	29	15	6
CHC* → CH* + C*	6	4	6	5
CC* → 2C*	0	0	0	0
CH2CH2* → CH2CH2(g)	1 115 396	1 055 554	895 274	686 353
2H* → H2(g)	3138	3126	3229	3372

---

Table 3 The final coverage of each species for different diffusion rates^a

Species	Diffusion rate (s^−1)
	10	1	0.1	0.01
a V denotes the vacancy.
CH2CH2	0.059	0.058	0.044	0.037
CH2CH	0.15	0.16	0.18	0.21
CHCH	0.0056	0.0042	0.0061	0.0061
CH2C	0.10	0.10	0.095	0.083
CHC	0.0049	0.0061	0.0071	0.0085
CC	0.31	0.31	0.32	0.34
CH2	0	0	0	0
CH	0.0083	0.0081	0.0056	0.0022
C	0.0059	0.0081	0.0042	0.0032
H	0.12	0.13	0.14	0.15
V	3.2	3.2	3.2	3.16

---

During the decomposition of ethylene to surface carbon species on Pd(100), the metal surface becomes increasingly less active over time because the final products, surface carbon species, poison the metal surface. Besides, with the elimination of surface H in the form of H₂ desorption, the entire process can be regarded as irreversible. Thus, the rates of most of the elementary reactions increase in a short time, quickly hit the peak and then start to decline with the formation of surface carbon species. Fig. 2 shows a series of snapshots of the species on Pd(100) surface at different simulation time. The final state of our kMC simulation is a surface carbon species poisoning state as is shown in Fig. 2f. Although the coverage of vacancies looks large compared with the coverage of the adsorbates, the vacancies are almost not available for further adsorption because the nearest adsorbed species will repel those species that attempt to be adsorbed in those vacancies. At this time, ethylene adsorption and ethylene desorption achieve a dynamic equilibrium almost without other elementary events occurring. The diffusion rate of H on the metal surface and the desorption rate of H₂ are rather faster than the rates of the dehydrogenation reactions and the C–C bond cleavage reactions. For the lack of a continuing availability of surface hydrogen, the rate of the hydrogenation reactions is also very slow, although the activation barriers of some of them are lower than those of the corresponding dehydrogenation reactions. The coverage of each species over different simulation time is shown in Fig. 3. As can be seen from Fig. 3, in the first 1.0 × 10⁵ s, the coverage of CH₂CH₂, CH₂CH, CH₂C, CC and H changes greatly, after 1.0 × 10⁶ s, the corresponding coverage changes slightly and after 1.0 × 10⁷ s, the fluctuation of the coverage of each species is negligible and keeps steady. The only fluctuation of ethylene is caused by its adsorption and desorption. Therefore, the results calculated with the simulation time of 1.0 × 10⁷ s are employed in the analysis below. The traffic of the network of ethylene decomposition on Pd(100) is shown in Fig. 4.

Fig. 2 The snapshots of the species on Pd(100) surface (light yellow-green ball: CH₂CH₂; dark blue ball: CH₂CH; dark green ball: CHCH; purple ball: CH₂C; red ball: CHC; light blue ball: CC; light green ball: CH; dark yellow-green ball: C; pink ball: H; gray ball: vacancy).

Fig. 3 The coverage of each species over different simulation time. (a) The coverage of each species over simulation time of 1.0 × 10⁵ s. (b) The coverage of each species over simulation time of 1.0 × 10⁶ s. (c) The coverage of each species over simulation time of 1.0 × 10⁷ s.

Fig. 4 The traffic of the network of ethylene decomposition on Pd(100).

On Pd(100) surface, two types of adsorption modes have been identified for ethylene adsorbed on transition metal surfaces: π adsorbed mode (ethylene sits at the top site with two carbon atoms of ethylene binding with one Pd atom) and 2σ adsorbed mode (ethylene sits at the bridge site with two carbon atoms of ethylene binding with two Pd atoms, respectively). π adsorbed ethylene is thought to form CH₂ through C=C bond cleavage (Rb), while 2σ adsorbed ethylene is regarded as the sources of dehydrogenation of ethylene (Ra).^9–11 Given that the activation barrier of C=C bond cleavage is too high with the value of 2.63 eV, and on Pd(100) 2σ adsorbed ethylene whose adsorption energy is 1.46 eV is more stable than π adsorbed ethylene whose adsorption energy is 1.36 eV, only 2σ adsorbed ethylene is considered in this work. As a result, the possible pathway starting from Rb is ruled out in the following work. Actually our calculation results show that such elementary events with relatively high activation barriers as most of the C–C bond cleavage reactions almost don't happen, which can justify our assumption in some extent. When ethylene is adsorbed on the Pd(100) surface, most ethylene molecules adsorbed on the surface desorb back into the gas phase because the rates of ethylene desorption and ethylene dehydrogenation are of the same magnitude in this work. Only a small portion of ethylene molecules (2550/897 993) convert to vinyl through dehydrogenation.

As a product of ethylene dehydrogenation, CH₂H can participate in the following reactions as a reactant on Pd(100) surface: dehydrogenation to CH₂C (Rc) or to CHCH (Rd), and C–C bond cleavage to CH₂ + CH (Re). The two dehydrogenation reactions are almost faster by six orders of magnitude than the C–C bond cleavage reaction. This is because the activation barrier, 2.63 eV, of the C–C bond cleavage is 0.54 eV and 0.64 eV higher than those of the former two dehydrogenation reactions. Actually in the present kMC calculations, no C–C bond cleavage reaction takes place for CH₂CH. The difference of 0.1 eV in activation barrier leads to an order of magnitude difference in the rates of the two dehydrogenation reactions. As can be seen from Fig. 4, most of CH₂CH (1649/2550) converts to CH₂C, while only a small part of CH₂CH (141/2550) transforms to CHCH.

There are two possible reactions beginning from CH₂C on Pd(100) surface: dehydrogenation to CHC (Rg), and C–C bond cleavage to CH₂ + C (Rf). The activation barrier of dehydrogenation to CHC, 1.48 eV, is 0.87 eV lower than that of C–C bond cleavage to CH₂ + C, 2.35 eV. Therefore, Rg is almost ten orders of magnitude faster than Rf. Virtually no Rf occurs in this work. Most CH₂C (1258/1649) turns into CHC through dehydrogenation.

It is likely for CHCH to be dehydrogenated to CHC (Rh), or generate CH + CH via C–C bond cleavage (Ri). The rate of (Ri) is slower than that of (Rh), since the activation barrier of the former, 1.71 eV, is higher than that of the latter, 1.49 eV. The majority of CHCH (1258/1649) undergoes dehydrogenation to yield CHC, while the minority of CHCH (15/141) proceeds by C–C bond cleavage to form CH, which is the precursor of carbon monomer.

CHC may be dehydrogenated to CC (Rj), or form CH + C through C–C bond cleavage (Rk). The rate of (Rj) is larger because of its lower activation barrier with the value of 1.15 eV. A large fraction of CHC (1325/1360) forms CC via dehydrogenation, while a small percentage of CHC (6/1360) undergoes C–C bond cleavage to generate CH and carbon monomer.

CC may convert to carbon monomer through C–C bond cleavage (Rl). Indeed, CC, or carbon dimer, stays on the surface as one of the final products, poisoning the Pd(100) surface, because the activation barrier of Rl is so high with the value of 2.32 eV that it is almost impossible for Rl to occur.

The results show that no CH₂ is generated during the simulation, suggesting that no corresponding C–C bond cleavage reactions involving CH₂ as the product occur. And the following CH is not formed through CH₂ dehydrogenation, but through the C–C bond cleavage of CHCH and CHC as mentioned above.

CH can be obtained from the C–C bond cleavage of CHCH (Ri) or the C–C bond cleavage of CHC (Rk). Besides, some CH (12/36) produces C monomer through dehydrogenation, while the rest remains on the surface for lack of available hollow sites for further dehydrogenation.

Carbon monomer, one of the final products, can derive from CH through dehydrogenation or form as one product of C–C bond cleavage reaction of CHC, poisoning the Pd(100) surface. The amount of the carbon monomer is relatively small compared with that of the carbon dimer.

The black arrows and the red arrows represent dehydrogenation reactions and C–C bond cleavage reactions, respectively. Those elementary events that are not exhibited here actually don't happen before the Pd(100) is totally poisoned. The number near the letter denotes the times the corresponding reaction happens during the simulation.

Hence, it can be directly and explicitly seen from Fig. 4 that the most probable pathway of ethylene decomposition on Pd(100) is CH₂CH₂ → CH₂CH → CH₂C → CHC → CC, among which the dehydrogenation of CH₂CH₂ to CH₂CH is the rate-limiting step with the activation barrier of 1.51 eV, followed by CH₂CH₂ → CH₂CH → CHCH → CHC → CC, whose rate-limiting step is the dehydrogenation of CH₂CH to CHCH with the activation barrier of 1.59 eV. The two most probable pathways produce carbon dimer as the final product, since the activation barrier of the C–C bond cleavage reaction is so high (2.32 eV) that it is almost impossible for it to occur before the metal surface is totally poisoned by surface carbon species. Another three feasible pathways are: (i) CH₂CH₂ → CH₂CH → CHCH → CH → C, (ii) CH₂CH₂ → CH₂CH → CHCH → CHC → CH + C → C and (iii) CH₂CH₂ → CH₂CH → CH₂C → CHC → CH + C → C, whose final product contain surface carbon monomer. And the reactions involving C–C bond cleavage are the rate-limiting step of the three pathways.

4. Conclusions

Based on the DFT calculations, we presented a comprehensive investigation on ethylene decomposition on Pd(100) with kMC simulation method. Our kMC simulations directly and explicitly show that the most probable pathway of ethylene decomposition on Pd(100) is CH₂CH₂ → CH₂CH → CH₂C → CHC → CC, among which the dehydrogenation of CH₂CH₂ to CH₂CH is the rate-limiting step with the activation barrier of 1.51 eV, followed by CH₂CH₂ → CH₂CH → CHCH → CHC → CC, whose rate-limiting step is the dehydrogenation of CH₂CH to CHCH with the activation barrier of 1.59 eV. The two most probable pathways produce carbon dimer as the final product, since the activation barrier of the C–C bond cleavage reaction is so high (2.32 eV) that it is almost impossible for it to occur before the metal surface is totally poisoned by surface carbon species. Another three feasible pathways are: (i) CH₂CH₂ → CH₂CH → CHCH → CH → C, (ii) CH₂CH₂ → CH₂CH → CHCH → CHC → CH + C → C and (iii) CH₂CH₂ → CH₂CH → CH₂C → CHC → CH + C → C, whose final product contain surface carbon monomer. And the reactions involving C–C bond cleavage are the rate-limiting step of the three pathways. Simple as the reaction network of ethylene decomposition looks, it is still difficult to analyze the decomposition mechanism merely according to the activation barriers from DFT calculations. Our work here demonstrates that kMC can nicely tackle the problem on competitive reaction pathways, each of which involves some reactions with relatively low activation barriers (e.g. the dehydrogenation reactions involved in ethylene decomposition) and some other reactions with relatively high activation barriers (e.g. the C–C bond cleavage reactions involved in ethylene decomposition).

References

1. R. W. Scott, C. Sivadinarayana, O. M. Wilson, Z. Yan, D. W. Goodman and R. M. Crooks, J. Am. Chem. Soc., 2005, 127, 1380.
2. Y. L. Yang, K. M. Saoud, V. Abdelsayed, G. Glaspell, S. Deevi and M. S. El-Shall, Catal. Commun., 2006, 7, 281.
3. J. K. Edwards, B. Solsona, P. Landon, A. F. Carley, A. Herzing, M. Watanabe, C. J. Kiely and G. J. Hutchings, J. Mater. Chem., 2005, 15, 4595.
4. D. I. Enache, J. K. Edwards, P. Landon, B. Solsona-Espriu, A. F. Carley, A. A. Herzing, M. Watanabe, C. J. Kiely, D. W. Knight and G. J. Hutchings, Science, 2006, 311, 362.
5. M. S. Chen, D. Kumar, C. W. Yi and D. W. Goodman, Science, 2005, 310, 291.
6. S. Nakamura and T. Yasui, J. Catal., 1970, 17, 366.
7. Z. X. Chen, H. A. Aleksandrov, D. Basaran and N. Rösch, J. Phys. Chem. C, 2010, 114, 17683.
8. D. Basaran, H. A. Aleksandrov, Z. X. Chen, Z. J. Zhao and N. Rosch, J. Mol. Catal. A: Chem., 2011, 344, 37.
9. V. Pallassana, M. Neurock, V. S. Lusvardi, J. J. Lerou, D. D. Kragten and R. A. van Santen, J. Phys. Chem. B, 2002, 106, 1656.
10. J. Andersin, N. Lopez and K. Honkala, J. Phys. Chem. C, 2009, 113, 8278.
11. L. V. Moskaleva, Z. X. Chen, H. A. Aleksandrov, A. B. Mohammed, Q. Sun and N. Rosch, J. Phys. Chem. C, 2009, 113, 2512.
12. F. Calaza, Z. Li and W. T. Tysoe, Catal. Lett., 2011, 141, 266.
13. B. Xing, Z. Z. Wei and G. C. Wang, J. Energy Chem., 2013, 22, 671.
14. A. F. Voter, Introduction to the kinetic Monte Carlo method, Radiation Effects in Solids, K. E. Sickafus, E. A. Kotomin and B. P. Uberuaga, Springer Netherlands, Dordrecht, 2007, pp. 1–23.
15. J. Rogal, K. Reuter and M. Scheffler, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 155410.
16. H. A. Aleksandrov, L. V. Moskaleva, Z. J. Zhao, D. Basaran, Z. X. Chen, D. Mei and N. Roesch, J. Catal., 2012, 285, 187.
17. F. Hess, A. Farkas, A. P. Seitsonen and H. Over, J. Comput. Chem., 2012, 33, 757.
18. L. Yang, A. Karim and J. T. Muckerman, J. Phys. Chem. C, 2013, 117, 3414.
19. Y. Li, S. H. Chan and Q. Sun, Nanoscale, 2015, 7, 8663.
20. H. Prats, L. Alvarez, F. Illas and R. Sayos, J. Catal., 2016, 333, 217.
21. B. Delley, J. Chem. Phys., 2000, 113, 7756.
22. B. Hammer, L. B. Hansen and J. K. Nørskov, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 7413.
23. M. H. Zhang, K. W. Yang and Y. Z. Yu, Appl. Surf. Sci., 2015, 328, 583.
24. D. H. Mei, M. Neurock and C. M. Smith, J. Catal., 2009, 268, 181.

---