Dynamical Study of the Dissociative Chemisorption of CHD3 on Pd(111).

The specific reaction parameter (SRP) approach to density functional theory has been shown to model reactions of polyatomic molecules with metal surfaces important for heterogeneous catalysis in the industry with chemical accuracy. However, transferability of the SRP functional among systems in which methane interacts with group 10 metals remains unclear for methane + Pd(111). Therefore, in this work, predictions have been made for the reaction of CHD3 on Pd(111) using Born-Oppenheimer molecular dynamics while also performing a rough comparison with experimental data for CH4 + Pd(111) obtained for lower incidence energies. Hopefully, future experiments can test the transferability of the SRP functional among group 10 metals also for Pd(111). We found that the reactivity of CHD3 on Pd(111) is intermediate between and similar to either Pt(111) or Ni(111), depending on the incidence energy and the initial vibrational state distribution. This is surprising because the barrier height and experiments performed at lower incidence energies than investigated here suggest that the reactivity of Pd(111) should be similar to that of Pt(111) only. The relative decrease in the reactivity of Pd(111) at high incidence energies is attributed to site specificity of the reaction and to dynamical effects such as the bobsled effect and energy transfer from methane to the surface.

Introduction

An important heterogeneously catalyzed industrial process is steam reforming, where methane and steam react over a metal catalyst (typically Ni[1]) and subsequently form carbon monoxide and hydrogen. At high temperature, the dissociation of methane, i.e., breaking the first CH bond, is a rate-controlling step in steam reforming on a wide variety of metals.[2,3] Therefore, a detailed study of the CH bond breaking is warranted to improve catalysts. However, the reaction of molecules on metal surfaces remains difficult to simulate due to the complexity of molecule–metal surface interactions.[4−8] The so-called specific reaction parameter (SRP) approach to density functional theory (DFT), though, has been shown to provide chemically accurate results, i.e., with errors smaller than 1 kcal/mol (4.2 kJ/mol), for a number of molecule–metal surface reactions.[9−14]

Within the SRP-DFT approach, two density functionals are mixed, of which one overestimates and one underestimates the reaction probability, according to an empirically determined parameter to create an SRP functional. Recently, an SRP functional was developed (the SRP32-vdW functional) that gave chemically accurate results not only for the molecule–surface reaction it was developed for (CHD3 + Ni(111)[12]) but also for methane interacting with a metal from the same periodic table group (CHD3 + Pt(111)[13]) and with a stepped surface of Pt (CHD3 + Pt(211)[13−15]). However, it remains unclear whether this transferability is common among all group 10 metals. Therefore, in this work, we perform predictive Born–Oppenheimer molecular dynamics (BOMD) calculations for CHD3 + Pd(111) with the SRP32-vdW functional in the hope that future experiments will test the transferability of the SRP functional describing methane interacting with all group 10 metal surfaces. Although in our previous work, we usually referred to our direct dynamics calculations with SRP functionals as “ab initio molecular dynamics” (AIMD) calculations, we have changed our wording of the method from AIMD to BOMD as “ab initio” can be misleading in the context of calculations based on a semiempirical density functional.

To ensure the validity of the BOMD method, we address conditions for which the total energy of the molecule (translational + vibrational) exceeds the minimum zero-point energy corrected barrier height of the system addressed. This ensures that the accuracy of the quasi-classical trajectory (QCT) method used in the BOMD dynamics is not much affected by quantum effects like tunneling and classical artifacts like zero-point energy violation.[16,17] Second, for laser-off conditions, we only address conditions where at least 60% of the incident CHD3 is in its initial vibrational ground state, and in predictions for initial-state selective reaction, we only address CH stretch excited CHD3, to avoid problems with artificial intramolecular vibrational relaxation that might otherwise affect QCT calculations.[18,19] Third, we address surface temperatures (here, 500 K) well above the surface Debye temperature (140 ± 10 K for Pd(111)),[20] thereby ensuring that the energy transfer between the molecule and surface can be well described with quasi-classical dynamics.[13,21,22]

Also, we will perform a rough comparison with existing experimental data for CH4 + Pd(111),[23] although a direct comparison is not possible due to the low experimental reaction probabilities making BOMD calculations untractable and the employed high-nozzle temperatures for which BOMD performs badly due to intramolecular vibrational-energy redistribution among excited vibrational states.[12]

Alloys are of special interest for catalysts[24] as they can increase both reactivity and selectivity.[25] For example, by combining a highly active metal like Pt with a less reactive metal such as Cu, a catalyst with high activity and selectivity can be produced, without the typical issues such as catalyst poisoning.[26] Recent work has predicted that the Pt–Cu(111) single-atom alloy is considerably more reactive than Pd–Cu(111), even though the barrier height difference is only 8.4 kJ/mol.[27] It was suggested that dynamical effects such as the “bobsled effect”[28,29] played a major role in the relatively lower reactivity of Pd–Cu(111) compared to that of Pt–Cu(111).[27] The so-called bobsled effect causes molecules with a high incidence energy to slide off the minimum energy path (MEP) for late barrier systems as the molecule is not able to make the turn before the barrier on the potential energy surface (PES) and therefore needs to overcome a higher barrier height than the minimum barrier height.[28,29] Since it was shown that the barrier geometries and potential energy surfaces (PES) above the doped atoms were similar to those found for the pure (111) surfaces of the respective doped elements, these dynamical effects can also be investigated by comparing methane interacting with Pd(111) and Pt(111).

The reaction of methane on metal surfaces remains fundamentally important due to many dynamically interesting effects. For example, in partially deuterated methane, the CH bond can selectively be broken by exciting the CH stretch mode.[12,13,30−33] Also, the dissociative chemisorption of methane is vibrational-mode-specific,[34,35] and the mode specificity is dependent on the metal surface.[35−37] Moreover, steric effects play a significant role.[38] Finally, the reaction of methane is site-specific.[2,13,15,39] For all of these reasons, we will present in this work a detailed analysis of the results from the BOMD calculations and compare them to the results obtained on Pt(111) and Ni(111).

This paper is structured as follows. Section will present the methods used in this work. In Section , the results are presented and discussed. Finally, in Section , a brief summary will be given.

Method

For the BOMD and electronic structure (DFT) calculations, the Vienna ab initio simulation package (VASP version 5.3.5)[40−44] is used. The first Brillouin zone is sampled by a Γ-centered 6 x 6 x 1 k-point grid, and the plane wave basis set kinetic energy cutoff is 400 eV. Moreover, the core electrons have been represented with the projector augmented wave method.[44,45] The surface is modeled using a four-layer (3 x 3) supercell, where the top three layers have been relaxed in the z direction and a vacuum distance of 13 Å is used between the slabs. Due to the computational cost, a small vacuum distance (i.e., 13 Å) is required, which effectively raises the barrier height by 4.9 kJ/mol. Therefore, 4.9 kJ/mol is added to the translational energy to counteract this shift, as done earlier for, e.g., CHD3 + Ni(111).[12] To speed up convergence, first-order Methfessel–Paxton smearing[46] with a width parameter of 0.2 eV has been applied. The employed computational setup is confirmed to be converged within chemical accuracy (1 kcal/mol or 4.2 kJ/mol), as shown by convergence tests provided in the Supporting Information.

The transition state is obtained with the dimer method[47−50] as implemented in the VASP transition-state tool package and is confirmed to be a first-order saddle point. Forces are converged within 5 meV/Å, where only the methane is relaxed.

We use the SRP32-vdW functional previously used for CHD3 + Ni(111), Pt(111), Pt(211), Pt(110), Pt(210), Cu(111), and Cu(211).[12−14,27,51−54] The exchange functional is defined aswhere EPBE and ERPBE are the exchange parts of the Perdew, Burke, and Ernzerhof (PBE)[55] and revised PBE (RPBE)[56] exchange–correlation functionals, respectively, and x = 0.32. Since it has been shown that modeling van der Waals interactions is vital for describing the reaction of methane on a metal surface,[13,14] the vdW correlation functional of Dion and co-workers (vdW-DF1)[57] is used.

A surface temperature of 500 K is simulated in the BOMD calculations, where the atoms in the top three layers are allowed to move and the expansion of the bulk due to the surface temperature is simulated by expanding the ideal lattice constant[58] (3.99 Å) by a factor of 1.0049.[59] The parameters used to simulate the molecular beams are taken from ref (13), which describes experiments performed for CHD3 + Pt(111). For every BOMD data point between 500 and 1000 trajectories were run, with a time step of 0.4 fs, for a maximum total time of 1 ps. A trajectory is considered to result in a reaction if a bond is 3 Å long, or longer than 2 Å for 100 fs, and in scattering if the molecule–surface distance is 6.5 Å and the velocity vector is pointing away from the surface. If neither has occurred after 1 ps, the trajectory is considered trapped. Other technical details of the BOMD calculations and the sampling of the initial conditions can be found in recent work.[12,13,16,27]

Results

Activation Barriers

The barrier heights and geometries of CHD3 on Pd(111) are compared to the barrier data on Ni(111) and Pt(111) in Table . θ is the angle between the dissociating bond and the surface normal, β is the angle between the surface normal and the umbrella axis, which is defined as the vector going from the geometric center of the three nondissociating hydrogen atoms to the carbon atom, and γ indicates the angle between the umbrella axis and the dissociating bond (see Figure ). The minimum barrier geometry on Pd(111) is similar to the minimum barrier geometry on Ni(111), with the main difference being that the barrier on Pd is at a larger distance from the surface than on Ni. However, the barrier height on Pd is much closer to that on Pt(111), being only 5.4 kJ/mol higher than on Pt(111). Based on the minimum barrier heights reported in Table , we would expect the reactivity of Pd(111) to be closest to that of Pt(111). Furthermore, the minimum barrier is located on the top site, which is typical for methane on a metal surface.[12,13,27,61]

Table 1

Minimum Barrier Geometries of Methane on Ni(111),[12] Pd(111), and Pt(111)[13]a

surface	site	ZC‡ (Å)	r‡ (Å)	θ‡ (deg)	β‡ (deg)	γ‡ (deg)	Eb (kJ/mol)
Ni(111)	top[12]	2.18	1.61	135.7	164.7	29.1	97.9 (85.3)
Ni(111)	fcc	2.09	1.63	128.5	157.3	30.7	121.1 (105.5)
Ni(111)	hcp	2.16	1.74	132.9	167.8	35.6	134.6 (120.7)
Ni(111)	bridge	2.06	1.65	126.3	154.8	29.5	135.1 (120.5)
Ni(111)	t2f	2.07	1.90	126.5	171.1	45.3	99.1 (88.8)
Ni(111)	t2b	2.12	1.63	130.4	160.0	31.0	113.9 (99.1)
Pd(111)	top	2.23	1.61	135.9	165.0	29.1	84.1 (70.1)
Pd(111)	fcc	2.14	1.73	133.0	160.8	27.8	132.6 (116.9)
Pd(111)	hcp	2.18	1.75	133.8	161.5	27.7	133.6 (118.1)
Pd(111)	bridge	2.14	1.76	130.8	161.9	31.1	125.6 (110.9)
Pd(111)	t2f	2.17	1.82	137.5	178.0	40.6	108.4 (96.1)
Pd(111)	t2b	2.18	1.76	132.8	165.8	33.0	132.5 (118.3)
Pt(111)	top[13]	2.28	1.56	133	168	35	78.7 (66.5)
Pt(111)	fcc	2.47	1.91	139.7	166.9	27.2	163.5 (145.8)
Pt(111)	hcp	2.59	1.90	122.1	161.2	39.1	158.0 (144.7)
Pt(111)	bridge	2.36	1.77	136.2	164.3	29.0	146.2 (128.1)
Pt(111)	t2f	2.31	1.64	149.5	179.2	29.7	117.7 (101.6)
Pt(111)	t2b	2.45	1.81	140.5	172.6	32.0	152.9 (136.5)

The zero-point energy corrected barriers are given in the brackets.

Figure 1

Transition state of methane on Pd(111), indicating the geometry angles as used in Table . θ is the angle between the CH bond and the surface normal, β is the angle between the umbrella axis and the surface normal, and γ is the angle between θ and β.

Moreover, barriers are also obtained above the fcc, hcp, bridge, top-2-fcc (t2f), and top-2-bridge (t2b) sites, by fixing the carbon in the X and Y directions above the aforementioned sites. For these barrier geometries, the angles are similar, but the dissociating bond length does increase, making the barrier even later. The distance of the carbon atom to the surface is smaller for Pd(111) and Ni(111) than at the top site, whereas in most cases, it is larger for Pt(111). The difference between the barrier heights obtained on Pd(111) and Ni(111) above these sites is considerably smaller than between the barrier heights at the top site. For Pt(111), the obtained barrier heights at the sites other than the top site are considerably higher than those of Pd(111) and Ni(111). The general trend observed here is that when going from Ni(111) to Pt(111), the difference between the barrier heights at the top site and at the other sites increases. We also note that, among the sites other than the top sites, the lowest barrier occurs on the t2f site for all metals. For Ni(111), this barrier is almost as low as the top site so that it may play an important role in the dynamics.

Finally, the adsorption energies of CH3 and H on Pd(111) are compared to those on Ni(111) and Pt(111) in Tables and 3. For CH3, Pd(111) is an intermediate of Ni(111) and Pt(111). The difference between the adsorption energies at the hollow and top sites is smaller for Pd(111) than for Pt(111), but for both, the preferred site is the top site, as opposed to Ni(111) where the preferred sites are the hollow sites. This may also explain why the barrier for dissociation on the t2f site is so low on Ni(111). However, Pd(111) is very similar to Ni(111) concerning the adsorption of hydrogen, where the binding of hydrogen to the top site is considerably weaker than to the other sites.

Table 2

Adsorption Energy of CH3 on Ni(111),[60] Pd(111), and Pt(111)[60]a

surface	site	ZC (Å)	adsorption energy (kJ/mol)
Ni(111)[60]	bridge	1.69	–155.2
Ni(111)[60]	fcc	1.55	–175.2
Ni(111)[60]	hcp	1.56	–172.5
Ni(111)[60]	top	1.98	–143.9
Pd(111)	bridge	1.85	–158.2
Pd(111)	fcc	1.75	–160.5
Pd(111)	hcp	1.77	–152.9
Pd(111)	top	2.09	–188.4
Pt(111)[60]	bridge	1.86	–120.2
Pt(111)[60]	fcc	1.78	–115.2
Pt(111)[60]	hcp	1.82	–105.4
Pt(111)[60]	top	2.10	–180.8

Note that the adsorption energies on Ni(111) and Pt(111) were calculated with the PBE functional.

Table 3

Adsorption Energy of H on Ni(111),[60] Pd(111), and Pt(111)[60]a

surface	site	ZH (Å)	adsorption energy (kJ/mol)
Ni(111)[60]	bridge	1.04	–256.4
Ni(111)[60]	fcc	0.91	–270.2
Ni(111)[60]	hcp	0.91	–269.3
Ni(111)[60]	top	1.47	–212.8
Pd(111)	bridge	0.98	–255.2
Pd(111)	fcc	0.81	–268.0
Pd(111)	hcp	0.81	–262.7
Pd(111)	top	1.56	–223.9
Pt(111)[60]	bridge	1.06	–256.5
Pt(111)[60]	fcc	0.92	–261.3
Pt(111)[60]	hcp	0.91	–256.5
Pt(111)[60]	top	1.56	–257.2

Note that the adsorption energies on Ni(111) and Pt(111) were calculated with the PBE functional.

Sticking Probability

Results for the reaction of methane on Pd(111) using BOMD are compared to those on Ni(111) and Pt(111) in Figure for laser-off conditions and ν1 = 1 (exciting the CH stretch mode with one quantum). Note that three additional points for Ni(111) have been calculated at ⟨Ei⟩ = 71.4, 89.2, and 101.1 kJ/mol for ν1 = 1 using the same computational setup as in ref (12). Additionally, results for ⟨EI⟩ = 146.6 kJ/mol were obtained in the original work of ref (12) but have not been reported before because there were no experimental data for this incidence energy. Contrary to expectations based on the minimum barrier heights only (see Table ), for laser-off conditions, the reaction probability on Pd(111) is similar to that on Ni(111). It should be noted that for Ni(111), a slightly higher surface temperature is used (550 K) than for Pd(111) and Pt(111) (500 K). However, this should not affect the results considerably as the surface temperature does not play a large role at high incidence energies, which will be discussed more in-depth in Section . For ν1 = 1 at lower incidence energy, the reaction probability is similar on all three systems investigated. Interestingly, on Pd(111), the reaction probability does not increase from 102 to 112 kJ/mol. It is possible that this is related to the site dependence of the reaction, which will be discussed later in Section . The generally much lower laser-off reactivity of Pd(111) compared to that of Pt(111) at high incidence energy is also consistent with the prediction that Pt–Cu(111) is much more reactive than Pd–Cu(111) at high incidence energies.[27] Finally, the trapping probabilities are not included in the reaction probability, as they are smaller than 0.5%.

Figure 2

Reaction probability of CHD3 on Ni(111) (blue), Pd(111) (black), and Pt(111) (red) for laser-off (a) and ν1 = 1 (b) using BOMD simulations. Results for Ni(111) and Pt(111) are taken from refs (12) and (13), respectively. The error bars represent 68% confidence intervals.

The bond selectivity is shown in Figure , where the fraction of CH bond cleavage under laser-off and state-resolved ν1 = 1 conditions is compared. When the CH stretch mode is excited, the dissociation of CHD3 is very selective toward CH cleavage, whereas under laser-off conditions, CH cleavage is close to statistical (25%). This is similar to what has been observed for CHD3 + Ni(111)[12,30] and CHD3 + Pt(111).[13] However, it remains unclear why on Pd(111) for laser-off conditions the fraction of CH cleavage is considerably lower for 112 kJ/mol compared to the other incidence energies under laser-off conditions. This may well be a statistical anomaly since a statistical analysis using Fisher’s exact test[62] cannot reject the null hypothesis that the CH dissociation probability is the same for all incidence energies. We also note that at higher incidence energies and laser-off conditions the CH cleavage ratio is somewhat lower than 0.25, which we attribute to the presence of CD-excited vibrational states in the beam,[12] also noting that there may be some artificial energy flow between these modes in classical dynamics calculations.

Figure 3

Fraction of reactions that occurred through CH bond cleavage for CHD3 on Ni(111) (blue), Pd(111) (black), and Pt(111) (red). Laser-off and ν1 = 1 results are indicated by solid and dashed lines, respectively. The error bars represent 68% confidence intervals.

Dynamics during the Reaction

The angles as indicated in Figure during the BOMD trajectories are shown in Figure and Table for the reacted trajectories. It is observed that both the initial θ and β angles, i.e., the angles that describe the orientations of the dissociating bond and umbrella axis, are close to the transition-state geometry. Moreover, during the dynamics, a considerable amount of bending between the dissociating bond and umbrella axis (γ angle) is observed. Finally, for all of the angles considered, some steering is observed, in the sense that at the time of reaction, the distributions describing the reacting molecules have moved somewhat toward the transition-state value of the angle described. However, the dynamics is not rotationally adiabatic (at the initial time step the orientational distribution of the reacting molecule is not statistical), as observed before for Ni(111)[12] and Pt(111).[13] This has consequences for how the rotations should be treated[5] in the reaction path Hamiltonian (RPH) approach of Jackson and co-workers.[63] Furthermore, the aforementioned dynamical behavior of the angles is not only typical for methane reacting on a group 10 metal surface (as can be seen in Figure ) but also for methane reacting on Cu(111).[27]

Figure 4

θ, β, and γ angles of methane during BOMD for all laser-off and ν1 = 1 reacted trajectories at the initial time step (dashed line) and when a dissociating bond reaches the transition-state value (solid line). The dotted lines indicate the transition-state values. Blue is Ni(111),[12] black is Pd(111), and red is Pt(111).[13]

Table 4

Average Value of the θ, β, and γ Angles with the Standard Error (σm) and Standard Deviation (σ) for All Laser-Off and ν1 = 1 Reacted Trajectories When a Dissociating Bond Reaches the Transition-State Value

surface	θ (deg) ± σm(σ)	β (deg) ± σm(σ)	γ (deg) ± σm(σ)
Ni(111)	117.0 ± 0.3 (11.3)	142.1 ± 0.4 (13.6)	31.3 ± 0.3 (12.4)
Pd(111)	123.5 ± 0.5 (11.0)	143.9 ± 0.6 (14.1)	27.9 ± 0.5 (11.4)
Pt(111)	123.5 ± 0.5 (10.1)	150.0 ± 0.6 (12.2)	34.1 ± 0.6 (12.8)

Although the barrier height on Pd(111) is considerably lower than on Ni(111), the barrier geometries are similar and thus dynamical effects such as the bobsled effect[28,29] would be expected to play similar roles. That the bobsled effect plays a role in the reaction of CHD3 on group 10 metal surfaces can be seen in Figure , where the distance of the carbon atom to the surface is shown when a bond dissociates. Both laser-off and ν1 = 1 trajectories that go on to react tend to slide off the MEP due to the bobsled effect and thus react over higher barriers. This deviation from the MEP increases with incidence energy, which is observed above all high-symmetry sites and thus is not related to the site over which the methane reacts. Furthermore, the bobsled effect is considerably smaller for Pt(111) than for Pd(111) and Ni(111), which leads to methane having to react over relatively higher barriers on Pd(111) and Ni(111) than on Pt(111).

Figure 5

Distance of the carbon atom to the surface when a bond dissociates, i.e., r = r‡, under laser-off conditions (solid lines) and ν1 = 1 (dashed lines). The blue squares, black circles, and red triangles indicate Ni(111), Pd(111), and Pt(111), respectively. The horizontal dashed lines indicate the transition-state value. The error bars represent 68% confidence intervals.

For similar values of the reaction probability, the bobsled effect on the reaction dynamics of CHD3 under laser-off conditions (predominantly ν1 = 0) is larger than for ν1 = 1. The reason is that larger incidence energy is required for ν1 = 0 to react than for ν1 = 1 so that ν1 = 0 CHD3 tends to slide further of the minimum energy path than ν1 = 1 CHD3. To observe this, see, e.g., Figure for Ni(111), observing the differences between laser-off conditions and ν1 = 1 for the lowest incidence energy for which a laser-off result is available on the one hand and for the lowest incidence energy for which a ν1 = 1 result is available on the other hand, and Figure to confirm that these conditions correspond to similar reaction probabilities. This has consequences for the vibrational efficacy, which is defined as the energy shift between the ν1 = 1 and 0 (≈laser-off) reaction probability curves divided by the energy difference between ν1 = 1 and 0, and defines how efficiently vibrational excitation promotes the reaction relative to increasing the translational energy. The larger bobsled effect on Ni(111) and Pd(111) than on Pt(111) partly explains why the vibrational efficacies for these systems (0.9–1.3 for Ni(111) and 0.7–0.9 for Pd(111)) exceed that obtained for Pt(111) (0.3–0.8, see Table S3, and also ref (12) for Ni(111) and ref (13) for Pt(111)). Furthermore, the large bobsled effect we find for CHD3 on Ni(111) is in line with one of the explanations Smith et al.[36] provided for the high vibrational efficacy of the asymmetric stretch mode of CH4 reacting on Ni(111), i.e., that ν3 = 1 CH4 reacts at the transition state, while ν3 = 0 CH4 slides off the MEP and has to pass over a higher barrier. We note that in the modeling of the reaction the molecule should be allowed to slide off the MEP to account for the bobsled effect on the vibrational efficacy. One reason that a too low vibrational efficacy was obtained for ν3 = 1 CH4 on Ni(111) in ref (64) may have been that the RPH calculations used a harmonic approximation for motion orthogonal to the MEP and an expansion in harmonic vibrational eigenstates with up to one quantum only in all modes combined. It is possible that such a limited expansion is not capable of describing the effect that the molecule may slide off the reaction path, as perhaps indicated by the reaction probability of methane in its vibrational ground state becoming smaller for particular incidence energies if the expansion is enlarged to also contain states with up to two vibrational quanta.[5]

Previously, it was observed that the minimum energy path (MEP) on Pd(111) is less favorable from a dynamical point of view than on Pt(111) due to the fact that the MEP makes a sharper turn on Pd(111) than on Pt(111).[27] Therefore, it is expected that at low incidence energies and ν1 = 1 where dynamical effects such as the bobsled effect are less important, the reactivity on Pd(111) is similar to that on Pt(111), whereas at higher incidence energies and laser-off conditions, dynamical effects cause the reactivity on Pd(111) to be similar to that on Ni(111) for the reaction of CHD3 in its vibrational ground state (to which laser-off reaction bears a close resemblance at a low nozzle temperature).

Another important dynamical aspect of the reaction of methane is the energy transfer from the molecule to the surface.[54]Figure compares for scattered trajectories this energy transfer from CHD3 to Cu(111),[54] Pt(111),[65] Ni(111),[12] and Pd(111), where the energy transfer ET is defined aswhere V and K are the potential free and kinetic energy of methane at the initial (i) and final (f) steps of the scattered trajectories, respectively. In general, it is observed that the lower the surface atom mass is, the higher the energy transfer is from methane to the surface. This is also predicted by the hard-sphere Baule model,[66] where the mass ratio between the molecule and the surface atom plays a large role in the energy transfer. This energy transfer is described bywhere μ = m/M (m is the mass of the projectile and M is the mass of a surface atom) and γ is the angle between the velocity vector of the molecule and the line connecting the centers of the hard spheres of the molecule and surface atom at impact. Surprisingly, the relatively simple Baule model does not only qualitatively but also semiquantitatively predict the energy transfer from methane to the metal surfaces, except for Ni(111), in contrast to what was previously predicted.[65] Typically, the Baule model is actually taken as an upper limit by treating the collision as a head-on collision (γ = 0), from which one obtains the well-known Baule limitHowever, when an empirical average for the γ angle distribution is used, in what we call the refined Baule model, the following average energy transfer (used in Figure ) is obtained[67]Considering the close to the spherical shape of methane, it is probable that the hard-sphere approximation made by the Baule model will typically hold. This is also suggested by Figure , which shows remarkably good agreement of the computed energy transfer with that predicted by the refined Baule model for Pt, Pd, and Cu. Additional work will be required to test the validity of the refined Baule model for other systems and investigate the considerably lower energy transfer we computed to Ni(111). Since the energy transfer from methane to Pd is higher than to Pt, less energy will be available for the reaction on Pd and thus the reaction probability should be further diminished on Pd compared to that on Pt. This effect will be larger at higher incidence energies as the difference in energy transfer between Pd and Pt will increase (see Figure ). Moreover, as the energy transfer to Pd(111) and Ni(111) is expected to be equal, differences in reaction probabilities on Pd(111) and Ni(111) are most likely not caused by the energy transfer from methane to the metal surface.

Figure 6

Energy transfer from methane to Ni(111) (blue squares), Pd(111) (black circles), Pt(111)[65] (red triangles), and Cu(111)[54] (green diamonds) compared to the refined Baule model. The solid lines without symbols indicate results predicted by the refined Baule model, whereas the dashed and dotted lines with solid and open symbols indicate laser-off and ν1 = 1 results, respectively.

As can be seen from Figure , at high incidence energy, the site over which CHD3 reacts on Pd(111) is close to statistical for both laser-off reaction and ν1 = 1. However, at lower incidence energy, it is observed that the top site is the most reactive site, followed by the bridge site. This means that at lower incidence energy mostly only the minimum barrier is accessed, since it is located at the top site as discussed in Section . Therefore, at lower incidence energies, a large portion of the surface would be catalytically inactive. This corresponds with the lack of increase in the reactivity of ν1 = 1 on Pd(111) from 102 to 112 kJ/mol, as it is also observed that the distribution of reaction sites shifts toward the less reactive sites (i.e., the bridge and hollow sites). Moreover, the reaction of CHD3 on Pt(111) shows a similar site-specific behavior as CHD3 reacting on Pd(111). At lower incidence energy, the reaction on Ni(111) again occurs predominantly over the top site; however, the second most reactive site is now the hollow site instead of the bridge site. In general, all of the considered metal surfaces show nonstatistical behavior, where the top site is usually favored, with the main difference being the ordering of the sites according to their reactivity. This behavior is also predicted by the site-specific barriers discussed in Section .

Figure 7

Fraction of closest high-symmetry site (i.e., top, hollow, and bridge) to the impact site of reacting methane on Ni(111) (blue), Pd(111) (black), and Pt(111) (red) for laser-off and ν1 = 1 as a function of the incidence energy when a bond dissociates, i.e., r = r‡. The dotted line indicates the statistical average for the high-symmetry site. The error bars represent 68% confidence intervals.

Figure shows the site-specific reaction probability of CHD3, of which the reaction probabilities add up to the total reaction probability. Again, we see that Ni(111), Pd(111), and Pt(111) exhibit similar site-specific reaction probabilities. Most of the reactivity is observed above the top site, whereas the hollow and bridge sites play a considerably smaller role. Here, the difference in reaction probability between Pd(111) and Pt(111) under laser-off conditions can be seen more clearly. The difference in reaction probability for the top site is large, whereas the difference for the hollow and bridge sites is generally much smaller. Therefore, the considerably lower reactivity of CHD3 on Pd(111) than on Pt(111) under laser-off conditions is mostly due to the difference in the top-site reactivity. However, this difference is not caused by the difference in minimum barrier heights; probably, it is caused by the difference in barrier heights that can be dynamically accessed due to the bobsled effect. Furthermore, it remains unclear whether the large variation in reaction probability for Pd(111) and Ni(111) at the top site for ν1 = 1 is a statistical anomaly or a systematic feature. Also, the partial contribution of each site is compared to the total reaction probability for each surface in Figure S2, which again shows the aforementioned differences in site-specific reactivity.

Figure 8

Reaction probability of CHD3 on the high-symmetry sites (i.e., top, hollow, and bridge) on Ni(111) (blue), Pd(111) (black), and Pt(111) (red) for laser-off and ν1 = 1 as a function of the incidence energy when a bond dissociates, i.e., r = r‡. The error bars represent 68% confidence intervals.

While the difference between the low vibrational efficacy computed for CHD3 + Pt(111) on the one hand and the higher vibrational efficacies on Pd(111) and Ni(111) on the other hand could be explained on the basis of the bobsled effect (see above), the reason for the higher vibrational efficacy on Ni(111) (0.9–1.3) than on Pd(111) (0.7–0.9, see Table S3) could not be explained in this way. On the basis of the minimum barrier heights and geometries collected in Table , it is tempting to speculate that the t2f site could play a role in this, as it has a much lower barrier on Ni(111) than on Pd(111), and a later barrier on Ni(111) than on Pd(111). The plot of the impact sites for the reactive trajectories with ⟨Ei⟩ = 89 kJ/mol for ν1 = 1 on Ni(111) (Figure S3) can be construed to offer some support for this idea, as quite a few reactive impacts are seen near the corners of the triangles making up the t2f and t2h sites. However, to gather further support for this idea, better statistics are needed, which could perhaps be obtained on the basis of QCT dynamics on a PES also incorporating the effect of surface atom motion.[54]

In the reaction of CHD3 on Pd(111), no steering in the xy plane is observed (on average, a movement of 0.06 Å in the xy plane), as is typical for the reaction of CHD3 on a metal surface.[5,12,27,51,52] As a result, it should be a good approximation to treat the reaction with a sudden approximation for motion in X and Y, as done, for instance, with the RPH model of Jackson and co-workers,[5] and firmly established to be valid for CH4 + Ni(111),[68] and also for H2O + Ni(111).[69]

Finally, we will summarize the general trends observed and how they affect the reaction probability, which are also shown in Table . First, the bobsled effect is considerably more important for Pd(111) and Ni(111) than for Pt(111), making Pt(111) considerably more reactive than the other surfaces, especially for laser-off conditions. Moreover, the energy transfer of methane to Pt(111) is smaller than to Pd(111) and Ni(111), again making Pt(111) relatively more reactive. However, the site-specific reactivity is increasingly more important when going from Ni(111) to Pt(111), reducing the reaction probability on Pt(111) the most. The vibrational efficacy plays an increasingly more important role when going from Pt(111) to Ni(111), increasing the reaction probability for ν1 = 1 on Ni(111) the most. Furthermore, the initial angular distribution of the molecule and concomitant steering are equally important on all surfaces considered here. These dynamical effects combined cause the reaction probability on Ni(111) and Pd(111) to be similar and on Pt(111) comparatively higher, for laser-off conditions. Additionally, they explain why the reactivity is rather similar on all of these surfaces for ν1 = 1. In this, we suspect that the site specificity plays the most important role in almost equalizing laser-off reaction on Pd(111) and Ni(111), while the vibrational efficacy should also be important to making the ν1 = 1 reaction probabilities almost equal on these two surfaces.

Table 5

Dynamical Features and How They Affect Qualitatively the Reaction Probability of CHD3 on Ni(111), Pd(111), and Pt(111)a

dynamical feature	Ni(111)	Pd(111)	Pt(111)	largest effect on
bobsled effect	–––	–––	–	laser-off
energy transfer	––	––	–	laser-off
site specificity	–	––	–––	laser-off
vibrational efficacy	+++	++	+	ν1 = 1
angular distribution	–	–	–	both

The number of pluses and minuses indicates how much the effect increases or reduces the reaction probability, respectively, when the aforementioned surfaces are compared.

Due to the combined effects of decreased site specificity and increased vibrational efficacy, it is conceivable that Ni(111) becomes more reactive than Pd(111), and/or Pd(111) becomes more reactive than Pt(111) toward ν1 = 1 CHD3 at higher incidence energies than results are shown for in Figure b. It would be a considerable challenge to explore this experimentally, for two reasons:[70,71] (i) at higher incidence energies, the extraction of the reactivity of ν1 = 1 CHD3 requires a subtraction of an increasingly large “laser-off” signal from a “laser-on” signal that might actually decrease, because laser excitation takes place from a rotational level that is less populated at the higher associated Tn, and (ii) the extraction requires the approximation that the reactivity of the vibrational ground state equals that averaged over the vibrational states populated in the beam under laser-off conditions, of which the validity decreases with incidence energy.

Discussion of Reactivity of Pd(111) vs Ni(111) and Pt(111); Comparison with Experiment

Experimentally, at low incidence energies (<70 kJ/mol) (see Figure ), the reactivity of Pd(111) toward CH4 is similar to that of Pt(111), whereas Ni(111) is about 3 orders of magnitude less reactive than Pt(111).[23,72−75] It should be noted that the experiments at low incidence energies were performed with CH4 using various nozzle and surface temperatures (see Table ), making a direct quantitative comparison between the experiments on CH4 + Pt(111) and CH4 + Pd(111) and with the BOMD results for CHD3 difficult. Therefore, we will discuss the general trends observed for the reaction of methane on Pt(111) and try to extrapolate this to Pd(111).

Figure 9

(a) Experimental reaction probability of CH4 on Ni(111) (blue), Pd(111) (black), and Pt(111) (blue) under laser-off conditions. Results for Ni(111) and Pt(111) are taken from ref (72) and refs (72−74), respectively. The Pd(111) results (black circles and triangles) are taken from ref (23), where the circles and triangles indicate incidence angles of 0 and 28°, respectively, and the black line is a linear regression fit those points. (b) Reaction probability of CH4 and CHD3 on Pd(111) and Pt(111) obtained with experiment (closed symbols) and BOMD (open symbols) under laser-off conditions. For CH4 + Pt(111), only the results from ref (74) are shown. The red squares and diamonds indicate results for CHD3 + Pt(111) taken from refs (16) and (13), respectively.

Table 6

Seeding Gas, Surface Temperature, and Nozzle Temperature Employed in the Experiments Shown in Figure

system	refs	seeding gas	surface temperature (K)	nozzle temperature (K)
CH4 + Ni(111)	Bisson et al.[72]	H2	475	323–373
CH4 + Pd(111)	Tait et al.[23]	He	550	470–885
CH4 + Pt(111)	Luntz et al.[73]	H2, He, Ar	800	300
CH4 + Pt(111)	Oakes et al.[74]	He	550	500–1000
CH4 + Pt(111)	Bisson et al.[72]	H2	600	323–373
CHD3 + Pt(111)	Nattino et al.[16]	He	120	500–850
CHD3 + Pt(111)	Migliorini et al.[13]	H2	500	400–650

In Figure b, a few results concerning Pt(111) and Pd(111) are shown to qualitatively compare the effect of nozzle and surface temperatures, and the isotopic effect of using CH4 or CHD3. Nattino et al.[16] used CHD3 seeded in a He beam with Ts = 120 K, whereas Migliorini et al.[13] used CHD3 seeded in a H2 beam with Ts = 500 K. Typically, at the high incidence energies and reaction probabilities involved here, the surface temperature does not play a large role for the reactivity of methane.[54,73,76] Moreover, at high surface temperature, the seeding gas influences the kinetic energy and thus also the required nozzle temperature. Therefore, the slightly higher reaction probability of Nattino and co-workers[16] in the overlapping regime is caused by the higher nozzle temperature (as needed by He-seeded molecular beam studies) as a larger fraction of CHD3 in the beam will be vibrationally excited.

However, the surface temperature can cause the reaction probability at lower incidence energy to vary by up to 2 orders of magnitude, depending on the surface temperature and incidence energy.[54,64,73,76,77] This surface temperature effect probably causes the reaction probabilities obtained by Luntz and Bethune[73] (Ts = 800 K) to be considerably higher than those obtained by Oakes and co-workers[74] (Ts = 550 K) and Bisson and co-workers[72] (Ts = 600 K), whom all used CH4. On the other hand, the higher reaction probability obtained by Oakes et al. (Tn = 500–1000 K) compared to that by Bisson et al. (Tn = 323–373 K) is probably due to the higher employed nozzle temperature used by Oakes et al.

Furthermore, the effect of partially deuterating methane can be seen by comparing the results of Nattino et al. and Oakes et al. For the incidence energy range where data are available for both sets, the difference in surface temperature (i.e., Ts = 120 and 550 K, respectively) should only play a role for the low incidence energies. Moreover, the nozzle temperature employed by Nattino et al. is similar to that by Oakes et al. and thus should not make a large difference either. It is expected that these differences should also (partially) cancel out at high incidence energies. It has also been shown previously that using CHD3 instead of CH4 lowers the reaction probability.[73,78−80] However, the reaction probabilities obtained by Nattino et al. and Oakes et al. at high incidence energy are similar, where it is expected that the reaction probabilities obtained by Oakes et al. should be slightly higher than those by Nattino et al. It remains unclear why no difference at high incidence energy is observed between the two data sets, although it is possible that the molecular beams are considerably different making direct comparison difficult.

Finally, the reaction probability of CH4 on Pd(111) obtained by Tait et al.[23] is similar to that of Oakes et al. for CH4 + Pt(111), except for the highest incidence energies where Pd(111) is measured to be more reactive than Pt(111) toward methane. Both used the same surface temperature and similar nozzle temperature range, but Tait et al. used relatively less seeding gas and thus a higher nozzle temperature is employed for given incidence energy compared to Oakes et al., which perhaps explains the higher reaction probability for Pd(111) at high incidence energy. At energies that are higher than those for which CH4 + Pd(111) experimental results are available, our BOMD calculations predict a substantially lower reactivity of Pd(111) toward CHD3 than that of Pt(111). While this may seem odd in light of the experimental results for CH4 on Pt(111) and Pd(111), one should keep in mind that due to the simulated use of H2 seeding the incidence energy is higher while the nozzle temperature is lower for the calculations on CHD3 + Pd(111) and Pt(111), which leads to a larger importance of the bobsled effect and to a smaller importance of the promotion of reaction by vibrational excitation. Both effects disfavor the reaction on Pd(111). Nevertheless, experiments are clearly needed to verify our predictions for the reaction of CHD3 on Pd(111). For all of these reasons, we conclude that, experimentally, it is expected that the reactivity of CHD3 + Pd(111) should be slightly lower than that of CHD3 + Pt(111) at lower incidence energies. Qualitatively, this is also what we obtain from the BOMD calculations at higher incidence energies, although there the difference in reactivity is larger (see Figure ).

Conclusions

In this work, predictive calculations using BOMD have been performed for CHD3 on Pd(111) with the SRP32-vdW functional. The reactivity of Pd(111) is compared to those of Pt(111) and Ni(111) and is found to be intermediate between these systems. Although this is to be expected from the minimum barrier heights and experiments at low incidence energy, the reaction probability is also found to be dependent on dynamical effects such as the bobsled effect and energy transfer from methane to the metal surface. In general, at the lowest incidence energy and laser-off conditions when these dynamical effects are smaller, the reaction probability on Pd(111) is comparable to that on Pt(111), which is also observed by experiment. However, at higher incidence energies, these dynamical effects play a larger role and the reaction probability is more comparable to Ni(111). Furthermore, for ν1 = 1, all three systems investigated show similar reaction probabilities. Moreover, barriers across the surface need to be considered as the reaction of methane on a group 10 metal surface is highly site-specific, with the minimum barrier height and geometry varying across the surface. This variation in barrier heights across the surface also explains the similarity of the reactivity of Ni(111) and Pd(111) toward methane at high incidence energy. Interestingly, methane on Pd(111) and Ni(111) exhibits typically quite similar dynamical behavior such as the bobsled effect, energy transfer from methane to the surface, and the site-specific reactivity, whereas the dynamical behavior of methane on Pt(111) tends to be different from the aforementioned metal surfaces. This again causes reactivity of Pd(111) toward methane to shift more to that of Ni(111) than that of Pt(111). Our results also suggest why Pt–Cu(111) is predicted to be much more reactive than Pd–Cu(111) at high incidence energy. We hope that these predictions will inspire new experiments that will test the transferability of the SRP32-vdW functional to CHD3 + Pd(111).