Quantum Dynamics of Dissociative Chemisorption of H2 on the Stepped Cu(211) Surface.

Reactions on stepped surfaces are relevant to heterogeneous catalysis, in which a reaction often takes place at the edges of nanoparticles where the edges resemble steps on single-crystal stepped surfaces. Previous results on H2 + Cu(211) showed that, in this system, steps do not enhance the reactivity and raised the question of whether this effect could be, in any way, related to the neglect of quantum dynamical effects in the theory. To investigate this, we present full quantum dynamical molecular beam simulations of sticking of H2 on Cu(211), in which all important rovibrational states populated in a molecular beam experiment are taken into account. We find that the reaction of H2 with Cu(211) is very well described with quasi-classical dynamics when simulating molecular beam sticking experiments, in which averaging takes place over a large number of rovibrational states and over translational energy distributions. Our results show that the stepped Cu(211) surface is distinct from its component Cu(111) terraces and Cu(100) steps and cannot be described as a combination of its component parts with respect to the reaction dynamics when considering the orientational dependence. Specifically, we present evidence that, at translational energies close to the reaction threshold, vibrationally excited molecules show a negative rotational quadrupole alignment parameter on Cu(211), which is not found on Cu(111) and Cu(100). The effect arises because these molecules react with a site-specific reaction mechanism at the step, that is, inelastic rotational enhancement, which is only effective for molecules with a small absolute value of the magnetic rotation quantum number. From a comparison to recent associative desorption experiments as well as Born-Oppenheimer molecular dynamics calculations, it follows that the effects of surface atom motion and electron-hole pair excitation on the reactivity fall within chemical accuracy, that is, modeling these effect shifts extracted reaction probability curves by less than 1 kcal/mol translational energy. We found no evidence in our fully state-resolved calculations for the "slow" reaction channel that was recently reported for associative desorption of H2 from Cu(111) and Cu(211), but our results for the fast channel are in good agreement with the experiments on H2 + Cu(211).

Introduction

The rate-limiting step in heterogeneous catalysis is often a dissociative chemisorption reaction.[1,2] Hydrogen (H2) dissociation is important to the heterogeneously catalyzed production of syngas and ammonia[3] and has recently gained industrial importance with the production of methanol from CO2 over a Cu/ZnO/Al2O3 catalysts, in which the rate-limiting step is considered to be the dissociation of H2.[4−6] Stepped, kinked, or otherwise defective surfaces more closely resemble real catalytic surfaces, as catalyzed reactions tend to proceed at the corners or edges of nanoparticles.[7,8] A better theoretical understanding of the reaction dynamics of H2 dissociation on stepped surfaces could well be a first step to the design of new catalysts from first principles.[9]

H2 reacting on copper surfaces is a prototypical example of a highly activated late barrier system.[10−13] For the flat Cu(111), Cu(110), and Cu(100) surfaces, a plethora of experimental[13−24] and theoretical[12,15−43] results have been reported, which are generally in good agreement with each other. This large body of work has allowed for the development of a chemically accurate description of molecular beam experiments using the semiempirical specific reaction parameter approach to density functional theory (SRP-DFT).[35] Recently, molecular beam adsorption experiments[44] and associative desorption experiments[45] for H2 reacting on Cu(211) have been reported, allowing for a more stringent comparison between theory and experiment for this system that more closely resembles a catalytic particle. Theoretical reaction dynamics of H2 reacting on stepped or defective surfaces have only been reported sparingly, most notably for D2 on Cu(211),[46] H2 on Pt(211),[47−51] and H2 on defective Pd(111).[52]

In our previous work, we and others have shown that the Cu(211) surface is less reactive than the Cu(111) surface,[46] which indicates that predictions based on the d-band model of No̷rskov and Hammer[53,54] are not always reliable. In the d-band model, increased reactivity at steps, defects, or otherwise less coordinated surface atoms is ascribed to a reduced width of the d band[53,55] and a shift of the center of the d band toward the Fermi level at these sites. In the case of Cu(211), the breakdown of the d-band model is due to the geometric effect of the lowest barrier to the reaction of H2 on Cu(111) not being situated at a top site.[46]

Due to the corrugated nature of the molecule surface interaction and the denser distribution of barriers to the reaction, it is unclear whether quantum effects can have a significant effect on the reaction dynamics of H2 reacting on Cu(211). Our main goal is to investigate if including quantum effects during the dynamics significantly affects observables such as the macroscopic molecular reaction probability and rotational quadrupole alignment parameters. To this end, we will mainly focus on a comparison of fully state-resolved quantum dynamical (QD) and quasi-classical trajectory (QCT) reaction probabilities for H2 incident on Cu(211), and the effect of Boltzmann averaging over all rovibrational states populated in a molecular beam experiment. Employing the time-dependent wave packet (TDWP) method,[56,57] we have carried out QD calculations mainly for H2. Due to the low mass of H2, quantum effects are presumed to be most prevalent for H2, and energy transfer to the surface during collision is expected to be small. Performing this large body of calculations for D2 would have been much more expensive because its larger mass necessitates the use of denser numerical grids and longer propagation times.

Another aim will be to investigate if the reaction dynamics of H2 dissociation on the stepped Cu(211) differs from the reaction dynamics at low Miller index copper surfaces, for which the reaction dynamics is reasonably similar.[30−32,43] This is relevant because the Cu(211) surface has Cu(111) terraces and Cu(100) steps, and considering this question might thus provide more insight in how a stepped surface can alter reaction mechanisms. Rotational quadrupole alignment parameters for vibrationally excited molecules are similar in behavior for Cu(111)[16,32] and Cu(100).[30,31] We will investigate whether the same holds for H2 + Cu(211).

Recent associative desorption experiments on Cu(111) and Cu(211),[45] which were in good agreement with earlier theoretical and experimental works,[16,30,34,35,46,58] have shown a never before reported “slow” reaction channel to be active for both Cu(111) and Cu(211). In this channel, the reaction could be facilitated by trapping on the surface and distortion of the surface due to thermal motion forming a reactive site.[45] Our calculations on sticking of H2 are carried out using the static surface approximation, which suggest that we might not be able to model this slow channel. We do however make a direct comparison to the experimental effective barrier heights obtained by applying the principle of detailed balance and direct inversion of time-of-flight measurements reported by Kaufmann et al.[45] for the fast channel.

The highly accurate potential energy surface (PES) used in our calculations and our previous work[46] has been constructed using the corrugation reducing procedure (CRP)[59] together with the SRP48 SRP density functional,[32] which was proven to be chemically accurate for H2 dissociating on Cu(111).[35] It has also been shown previously that the SRP functional for H2 + Cu(111) is transferable to H2 + Cu(100).[30] All our calculations have been carried out using the BOSS model, which works well for activated H2 dissociation on metals at low surface temperatures.[26,33−36,60]

This paper is organized as follows. Section will outline the computational methods used, with Section detailing the coordinate system used and Section describing the AIMD calculations. Sections and 2.D describe the QCT and QD methods, respectively, while Section describes the calculations of observables. Section is the results and discussion section. Section is a comparison between QD and QCT reaction probabilities at the fully state-resolved level. Sections presents calculated rotational quadrupole alignment parameters for both H2 and D2. In Section , we compare theory to the experimental effective barrier heights reported by Kaufmann et al.[45]Section presents a comparison of BOMD, QCT, and molecular dynamics with electronic friction (MDEF) calculations for D2 in order to highlight the extent to which surface atom motion and electron–hole pair (ehp) excitation can be expected to affect the reaction probability in molecular beam experiments. In Section , we present fully quantum dynamical molecular beam simulations for H2 reacting on Cu(211), comparing to QCT calculations. Section presents conclusions.

Computational Methods and Simulations

In the following, we present details about the different simulations we have performed to describe the dynamics of H2(D2) incident on Cu(211). In our six-dimensional QD, QCT, and MDEF simulations, we used the static surface approximation. They are carried out on a six-dimensional PES that was previously developed by us[46] on the basis of the corrugation reducing procedure[49] and ∼116,000 DFT energy points computed with the SRP48 functional.[32] The SRP48 functional contains 48% RPBE[61] and 52% PBE[62] exchange correlation and was fitted to quantitatively reproduce experimental sticking probabilities for the reaction of H2(D2) on a flat Cu(111) surface.[32] The very similar SRP functional[35] performed well at describing the H2 + Cu(100) reaction.[30]

Coordinate System

The six-dimensional dynamics calculations account only for the motion along the six molecular degrees of freedom (DOF) of H2(D2), while the surface atoms are kept frozen at their ideal 0 K configuration as computed with DFT. The molecular coordinates include the center of mass (COM) position given by the coordinates X, Y, and Z, where Z is the molecule–surface distance, and X and Y are the lateral positions measured relative to a Cu reference atom at the step edge. Also included are the H–H bond distance r and the angular orientation of H2 given by the polar angle θ defined with respect to the surface normal and the azimuthal angle ϕ. The coordinate system is drawn in Figure a, and the Cu(211) surface unit cell is drawn in Figure b; additional details about the dimensions of the (1 × 1)Cu(211) unit cell are specified in the corresponding caption.

Figure 1

Coordinate system for H2(D2) on Cu(211). H atoms are drawn in blue, and Cu atoms are drawn in brown. Shown are (a) a side view on a (1 × 2)Cu(211) supercell and (b) a top view on a (1 × 1) unit cell. The six-dynamical molecular DOF are indicated, that is, the COM coordinates given by X, Y, and Z, where X and Y are the lateral coordinates and Z is the molecule surface distance. Further, the H–H bond distance is represented by r, and the angular orientation is represented by the polar angle θ and the azimuthal angle ϕ. The latter is defined with respect to the X axis, and the former is defined with respect to the macroscopic surface normal. The computed lengths of the lattice vectors of the (1 × 1) unit cell are L = 6.373 Å and L = 2.602 Å along X and Y.

Ab Initio Molecular Dynamics Simulations

To describe the reaction of D2 on Cu(211) at normal incidence with the BOMD technique, we employ a modified version of the Vienna ab initio simulation package[63−66] (VASP). Note that, in previous publications, we referred to the direct dynamics technique using SRP-DFT as the ab initio molecular dynamics (AIMD) method. Because this might be taken to imply that the SRP functional is not semiempirical, we abandoned this name and now refer to it as Born–Oppenheimer molecular dynamics (BOMD). The modifications of the computer package concern the propagation algorithm and were first introduced in previous works[67,68] on electronically nonadiabatic effects in gas–surface systems using VASP. To be consistent with our previous work on the system,[46] we adopt the same computational setup for the electronic structure calculations specified in the Supporting Information of ref (46). Here, we briefly recall only the most important details. The Cu(211) surface is represented using a five-layer slab model periodically repeated over a (1 × 2) supercell with a vacuum spacing of 15 Å. Ultrasoft pseudopotentials are used as well as plane waves corresponding to energies of up to 370 eV. The k-points are sampled using the Monkhorst grid scheme and an 8 × 8 × 1 mesh centered at the Γ point. Fermi smearing is used with a width of 0.1 eV.

BOMD simulations are performed at different average incidence energies and mimic corresponding molecular beam conditions at which Michelsen et al.[15] originally performed experiments on the dissociation of D2 on flat Cu(111). The inclusion of beam parameters in the simulations is explained below in Section . For each incidence energy point, we perform 500 trajectory calculations. This allows us to achieve an absolute standard error of smaller than 0.02 in the computed initial sticking coefficient. All BOMD trajectories start at a molecule–surface distance of Z = 7 Å and are propagated until dissociation or scattering of D2 has occurred. Here, we count trajectories to be dissociatively adsorbed if the D–D bond distance r is larger than 2.45 Å. A nonreactive scattering event is counted when trajectories return to the gas phase and have reached a molecule–surface distance of Z ≥ 7.1 Å. We use a time-step discretization Δt of 1 fs in the dynamics propagation and a maximum propagation time tf of 2 ps. Geometries between consecutive time steps are updated if the electronic structure energy is converged to 10–5 eV. The setup allows, on average, for an energy conservation error of typically ∼10 meV.

BOMD simulations performed within the static surface approximation employ the same slab model described in our earlier work.[46] Therein, the first four layers of the slab are relaxed through energy minimization (the positions of the fifth layer atoms are fixed during relaxation). The resulting optimized Cu(211) surface conserves the P1m1 space group and remains unchanged during the BOMD simulations. This prevents energy transfer to take place between the molecule and the surface due to excitation of surface atom motion upon scattering. To model a thermalized Cu(211) surface at a temperature Ts of 120 K according to experiments, we follow the NVE/NVT procedure explained in refs (32). and (69) and generate 10,000 slab configurations resembling the phase space. The initial condition of a BOMD trajectory at Ts = 120 K is set up by randomly mixing thermalized slab models with random configurations of D2 generated according to the molecular beam conditions.

Quasi-Classical Simulations

The MD(EF) simulations presented in this work use the 6D PES of ref (46) and assume quasi-classical conditions,[70] that is, initial conditions of the classical trajectories reflect the quantum mechanical energies of incident H2 (D2) in their initial rovibrational state(s). To do so, we use the method described in ref (69). The dynamics is studied by integrating a Langevin equation[71] numerically using the stochastic Ermak–Buckholz algorithm,[72] and the methodology is outlined in refs (69) and (73). Note that, in the nondissipative limit, that is, the MD case, the Langevin equation obeys Newton’s equation of motion for which the propagation algorithm is also suitable. In the MDEF case, energy dissipation between the molecule and surface is mediated through electronic friction as computed from the local density friction approximation within the independent atom approximation (LDFA-IAA) model.[74] Specifically, friction coefficients of the hydrogen atoms are represented as a function of the electron density of the ideal bare Cu(211) surface. The latter is extracted from a single DFT calculation (see ref (69) for details).

QCT calculations are used here (i) to model fictitious molecular beam experiments using realistic beam parameters and (ii) to perform initial state-resolved calculations. In the former case, 100,000 QCT calculations per energy point are computed, whereas state-resolved sticking coefficients are evaluated per energy point over 50,000 trajectories. As with BOMD, all MD(EF) trajectories start at a molecule-surface distance of Z = 7 Å. A time step of Δt = 0.5 ℏ/Eh (≈0.012 fs) is used for the propagation resulting in an energy conservation error for the MD simulations of smaller than 1 meV. To determine dissociative adsorption and nonreactive scattering, we impose the same conditions used for the BOMD simulations (see above).

Quantum Dynamics Simulations

To perform 6D quantum dynamics simulations, we solve the time-dependent Schrödinger equationusing the time-dependent wave packet (TDWP) approach as implemented in our in-house computer package.[56,57] In eq , Q̲=(X,Y,Z,r,θ,ϕ) is a six-dimensional position vector, Ψ(Q̲;t) is the time-dependent nuclear wave function of the system, and Ĥ(Q̲) is the time-independent Hamiltonian, which readsHere, M and μ are the mass and the reduced mass of H2, respectively, and ∇ and Ĵ are the nabla and the angular momentum operators. The 6D PES, V(Q̲)=V(X,Y,Z,r,θ,ϕ), is taken from ref (46) and was computed with the SRP48 functional.[32] The initial wave function is represented as a product of a Gaussian wave packet u(Z0, k0) centered around Z0, a two-dimensional plane wave function ϕ(k0, k0) along X and Y, and the rovibrational wave function ψν, (r, θ, ϕ) of incident H2where the two-dimensional plain wave function and the Gaussian wave packet are defined asHere, σ is the width of the wave packet centered around the wave vector , and k0 are the initial wave vectors of the COM. The width σ is chosen in such a way that 90% of the Gaussian wave packet is placed in an energy range E ∈ [Emin, Emax]. Equation is solved numerically using the split operator method with a time step Δt. We apply a quadratic form of optical potentials[75] in the scattering (large values of Z) and adsorption regions (large values of r). The scattered fraction of the wave function is analyzed through the scattering matrix formalism,[76] and the scattering probability Psc is computed accordingly. Substracting Psc from 1 then yields the sticking probability S0.

Parameters for the wave packet calculations defining the initial wave packet, grid representation, time step, and the optical potentials are compiled in Table . The final propagation time can vary since we stop simulations if the remaining norm on the grid is below 0.01.

Table 1

Input Parameters for the 6D Quantum Simulations on the Reactive Scattering of H2 on Cu(211)a

	0.05–0.22 eV		0.2–0.6 eV		0.57–1.4 eV			D2
parameters	ν0	ν1	ν0	ν1	ν0J ∈ [0, 7]	ν0J ∈ [8, 11]	ν1	ν1J6
Zstart (bohr)	–2.0	–2.0	–2.0	–2.0	–2.0	–2.0	–2.0	–2.0
NZspec	280	280	280	280	280	280	280	280
NZ	180	180	176	176	176	176	176	176
ΔZ (bohr)	0.1	0.1	0.08	0.08	0.08	0.08	0.08	0.08
Rstart (bohr)	0.8	0.8	0.8	0.8	0.8	0.8	0.8	0.8
NR	60	60	56	56	56	56	56	56
ΔR (bohr)	0.15	0.15	0.15	0.15	0.15	0.15	0.15	0.15
NX	36	36	36	36	36	36	36	42
NY	12	12	12	12	12	12	12	16
NJ	26 / 25	30 / 29	26 / 25	32 / 31	38 / 37	42 / 41	36 / 35	42
NmJ	26 / 25	30 / 29	26 / 25	32 / 31	30 / 29	42 / 41	28 / 27	40
Complex absorbing potentials
ZCAP start (a0)	8.9	8.9	8.88	8.88	8.88	8.88	8.88	8.88
ZCAP end (a0)	15.9	15.9	12.0	12.0	12.0	12.0	12.0	12.0
ZCAP optimum (eV)	0.16	0.16	0.3	0.3	0.95	0.95	0.95	0.3
ZspecCAP start (a0)	18.1	18.1	16.8	16.8	18.16	18.16	18.16	16.8
ZspecCAP end (a0)	25.9	25.9	20.32	20.32	20.32	20.32	20.32	20.32
ZspecCAP optimum (eV)	0.16	0.16	0.3	0.3	1.2	1.2	1.2	0.3
RCAP start (a0)	4.55	4.55	4.55	4.55	4.55	4.55	4.55	4.55
RCAP end (a0)	9.65	9.65	9.05	9.05	9.05	9.05	9.05	9.05
RCAP optimum (eV)	0.12	0.12	0.3	0.3	1.0	1.0	1.0	0.3
Propagation
Δt (ℏ/Eh)	2	2	2	2	2	2	2	2
tf (ℏ/Eh)	44000	44000	14000	14000	10000	10000	10000	20000
Initial wave packet
Emin (eV)	0.05	0.05	0.2	0.2	0.57	0.57	0.57	0.2
Emax (eV)	0.22	0.22	0.6	0.6	1.4	1.4	1.4	0.6
Z0 (a0)	13.50	13.5	11.44	11.44	11.44	11.44	11.44	11.44

All wave packets were propagated until the remaining norm was less than 1%.

Computation of Observables

To incorporate the effect of a molecular beam on the computed sticking coefficient, we need to take into account the distributions of translational energies and rovibrational state population due to a nozzle temperature Tn. The probability to find a molecule with velocity v + dv and a rovibrational state described by the vibrational quantum number ν and the angular momentum quantum number j here is given bywhere the flux-weighted velocity distribution Pflux is a parameterized function of Tn and determined by the width parameter α and the stream velocity v0 according to[77]where C is a normalization constant. The ensemble representation of the rovibrational state population distribution readswithHere, kB is the Boltzmann constant, and Eν,  is the energy of the quantum state characterized by ν and j. The first and second Boltzmann factors describe the vibrational and rotational state populations, respectively. Note that the rotational temperature is Trot = 0.8Tn,[18] whereas the vibrational temperature applies Tvib = Tn. This setting is in agreement with the observation that rotational but no vibrational cooling occurs during gas expansion in the nozzle. The factor w(J) in eq is due to ortho- and para-hydrogen molecules present in the beam. For H2, w(J) is 1/4 (3/4) for even (odd) values of J, and for D2, w(J) = 2/3 (1/3) for even (odd) values of J.

In the case of classical dynamics calculations (MD, MDEF, and BOMD), the probability distributions P(v, ν, J, Tn) is randomly sampled as described in ref (69) using the different beam parameters on H2 and D2 listed in Table . The sticking coefficient per energy point is given by the ratio of the number of adsorbed trajectories Nads and the total number of computed trajectories N, that is, S0 = Nads/N. To extract quantum mechanical results on H2 beam simulations, a direct sampling of P(v, ν, J, Tn) is not feasible. Instead, initial state-resolved reaction probabilities Rmono(E, ν, J) are first computed as functions of the monochromatic incidence energy E by degeneracy averaging fully initial state-resolved reaction probabilities PR(E, ν, J, m) over the magnetic rotational quantum number m, that is,

Table 2

Molecular Beam Parameters Taken from Experiments Performed on the H2(D2) + Cu(111) Systema

Tn (K)	⟨Ei⟩ (kJ/mol)	v0 (m/s)	E0 (eV)	α (m/s)
Seeded molecular H2 beams (Ts = 120 K)
1740	19.9	3923	0.160	1105
1740	28.1	4892	0.250	1105
1740	38.0	5906	0.364	945
2000	18.2	3857	0.155	995
2000	25.1	4625	0.223	1032
2000	44.1	6431	0.432	886
Seeded molecular D2 beams (Ts = 120 K)
2100	62.6	5377	0.829	649
2100	69.2	5658	0.860	717
2100	80.1	6132	0.849	830
Pure molecular H2 beam (Ts = 120 K)
1435	31.7	5417	0.307	826
1465	32.0	5446	0.310	830
1740	38.0	5906	0.364	945
1855	40.5	6139	0.394	899
2000	44.1	6431	0.432	886
2100	47.4	6674	0.465	913
2300	49.7	6590	0.454	1351
Pure molecular H2 beam (Rendulic et al.)
1118.07	25.1	3500	0.12794	1996
1331.89	29.9	3555	0.13200	2342
1438.82	32.3	3380	0.11932	2611
1501.19	35.7	3151	0.10371	2819
1581.35	35.5	3219	0.10816	2903

The parameters are used in this work to simulate the reaction of molecular hydrogen on Cu(211) as it would occur in experiments analogous to those performed on Cu(111). The parameters v0, α, Tn represent the stream velocity of the beam, the width of the beam, and the nozzle temperature at an average translational incidence energy ⟨E⟩, respectively. Parameters were taken from refs (15) and (34).

The initial sticking probability S0(⟨E⟩) is then calculated as a function of average incidence energy ⟨E⟩ by averaging over the rovibrational (ν, J) states populated in the beam (see eq ) and the flux-weighted distribution of the incidence translational energies of the beam according to

We note that, although S0(⟨E⟩) is written and plotted in studies as a function of average incidence only, it also implicitly depends on Tn through the distribution P′(E, ν, J, Tn) of incidence energies and the rovibrational state populations

P′(E, ν, J, Tn) makes the initial sticking also depend implicitly on incident beam conditions other than just Tn due to the occurrence of the flux-weighted distribution of incidence energies Pflux′(E; Tn), which depends on a number of factors including the molecular beam geometry, backing pressure, and whether or not a seeding gas is used and can be described by the parameters E0 and ΔE0 according to

Instead of averaging over incidence energies using Pflux′(E; Tn) as done in eq , it is also possible to average over the flux-weighted velocity distribution of the molecules in the beam, Pflux(v; Tn), and the derivation Pflux′(E; Tn) from Pflux(v; Tn) is discussed in ref (77). For a particle of mass m, the parameters are defined as E0 = mv02/2 and ΔE0 = 2E0α/v0.

To obtain sticking coefficients S0, we perform 114 state-resolved calculations (corresponding to 342 wave packet calculations) for an energy range of E ∈ [0.05,1.4] eV. The initial states of incident H2 considered here to evaluate eq are characterized by the quantum numbers J ∈ [0,11] for ν = 0 and J ∈ [0,7] for ν = 1 and m ∈ [0, J].

The rotational quadrupole alignment parameter as a function of ν and j is a measure of the extent to which the reaction depends on the orientation of the molecule. The rotational quadrupole alignment parameter is calculated from the fully state-resolved reaction probability as follows[78]

Results and Discussion

Fully State-Resolved Reaction Probabilities

In order to highlight the difference between a QD and QCT treatment of the H2 + Cu(211) system, we first present initial state-resolved reaction probabilities in Figure a–c. QD calculations have been performed for a large number of rovibrational states. All input parameters can be found in Table . The biggest differences between QD and QCT calculations at the fully state-resolved level are observed for the lowest rovibrational states, as shown in Figure b,c. The differences get increasingly smaller with increasing J for J > 1. From QCT data at higher translational energies that are not shown in this figure, it is clear that all states converge toward an asymptotic maximum reaction probability, which depends slightly on the rovibrational state with respect to the maximum reaction probability. We note that for a very high J, J > 10 (not shown here), QD predicts a marginally smaller (less than 2%) asymptotic maximum reaction probability, while Figure c suggests that the opposite is true for the vibrational ground state and the first vibrationally excited state.

Figure 2

Reaction probability computed with QD calculations (solid lines) and QCT calculations (dashed lines) for normal incidence. Panel (a) shows degeneracy averaged reaction probabilities for J = 1 for both the ground state and the first vibrationally excited state. Panel (b) shows the m = 0, 1 states belonging to J = 1 for both the ground state and the first vibrationally excited state. Panel (c) shows the J = 0 state for both the ground state and the first vibrationally excited state as well.

Figure b shows the largest discrepancy between the QCT and QD calculations observed. Here, |m| = J pertains to a “helicoptering” H2 molecule, and m = 0 pertains to a “cartwheeling” H2 molecule rotating in a plane perpendicular to the surface normal. The preference for reacting parallel to the surface (i.e., m = J having a higher reaction probability than m = 0) is bigger for QD calculations than QCT calculations. This difference is negligible, however, when looking at degeneracy averaged reaction probabilities, which are shown in Figure a. This also holds for the states not shown here. When looking at degeneracy averaged reaction probabilities, the agreement between the QCT and QD method is excellent.

In our calculations, we see no evidence of the “slow channel” reactivity reported by Kaufmann et al.[45] in their very recent paper, that is, reaction at low translational energies. We can now rule out quantum effects during the dynamics as the source of this slow channel reactivity, in which the reaction supposedly is inhibited by translational and promoted by vibrational energy.[45] When looking at the individual rovibrational states that exhibit the biggest difference in reactivity between QD and QCT calculations, no evidence of the slow reaction channel is present in our results. The translational energy range sampled in our calculations should overlap with the translational energy range where the slow channel is reported to be active by Kaufmann et al.[45] We therefore propose that the observed slow reaction channel must originate from surface motion at a very high surface temperature (923 K), which has not been incorporated into our QD calculations and is challenging to incorporate in QCT calculations.[79]

Rotational Quadrupole Alignment Parameters

As might be suspected from Figure b from the larger preference for a parallel reaction orientation for J = 1, calculated rotational quadrupole alignment parameters show a large difference between QCT and QD calculations for the J = 1 states shown there. However, here, we will now focus on rotational quadrupole alignment parameters for two particular rovibrational states of H2, (ν = 0, J = 7) and (ν = 1, J = 4) (Figure a), and those of D2, (ν = 0, J = 11) and (ν = 1, J = 4) (Figure b). These two sets of states were selected because they are very similar in rotational energy to the two rovibrational states for which rotational quadrupole alignment parameters for D2 desorbing from Cu(111) have been measured experimentally[16] and studied theoretically using the BOMD method.[32] Results for both states of D2 reacting on Cu(111) have been included in Figure b. Note that a positive A0(2)(ν, J) indicates a preference for a parallel reaction orientation, a negative value indicates a preference for a perpendicular orientation, and zero means the reaction proceeds independent of the orientation.

Figure 3

Panel (a) shows rotational quadrupole alignment parameters, A0(2)(ν, J), for two rovibrational states of H2: (ν = 0, J = 7) and (ν = 1, J = 4). Panel (b) shows rotational quadrupole alignment parameters for two rovibrational states of D2: (ν = 0, J = 11) and (ν = 1, J = 6). Solid lines correspond to QD calculations, and dashed lines correspond to QCT calculations. Panel (b) also shows experimental results for D2 on Cu(111) (black line)[16] and BOMD results for D2 (ν = 1, J = 6) on Cu(111) (green line).[32]

We observe that the predicted rotational quadrupole alignment parameters eventually tend to zero with increasing translational energy, as all molecules irrespective of the orientation will have enough energy to traverse the barrier. It is also clear that for H2 (ν = 0, J = 7), the agreement between QCT and QD calculations is excellent. The slight deviations at the lowest translational energies can be attributed to noise in the very low reaction probabilities of the underlying individual states.

The increase in the rotational quadrupole alignment parameter with decreasing translational energy, for the H2 (ν = 0, J = 7) and D2 (ν = 0, J = 11) states, is comparable to what is reported in the literature for H2 and D2 associatively desorbing from Cu(111) and Cu(100).[16,30−32] This monotonic increase in the rotational quadrupole alignment parameter with decreasing translational energy can be explained by a static effect of orientational hindering, in which slow- or nonrotating molecules will scatter when their initial orientation does not conform to the lowest barrier geometry.[31] Specifically, the molecule must be in favorable orientation to begin with in order to react, especially with the energy available to reaction being close to the threshold energy.

The blue lines in Figure a correspond to the (ν = 1, J = 4) rovibrational state of H2 and those in Figure b to the (ν = 1, J = 6) rovibrational state of D2. In contrast to the previously described states (ν = 0, high J), the rotational quadrupole alignment parameter now first increases with increasing translational energy until reaching a maximum around 0.43 eV for H2 and 0.52 eV for D2 before decreasing toward zero with increasing translational energy. From Figure , it is clear that, around the maximum, the agreement between the QD and QCT calculations is not as excellent for H2 than D2, although the agreement is still good.

The downturn of the rotational quadrupole alignment parameter with decreasing translational energy seen here for D2 and H2 in their (ν = 1) states colliding with Cu(211) was not observed for D2 desorbing from Cu(111) for which, as can be seen in Figure b, only a monotonous increase with decreasing translational energy has been reported.[16,32] A slight downturn of the rotational quadrupole alignment parameters has been predicted for vibrationally excited H2 reacting on Cu(100),[30,31] although the downturn was too small to lead to a negative rotational quadrupole alignment parameter. Because the behavior predicted for (ν = 1) hydrogen colliding with Cu(211) qualitatively differs from that observed previously for Cu(111) and Cu(100), we will now first attempt to explain the dependence of the rotational quadrupole alignment parameter on incidence energy that we predict for D2 (ν = 1, J = 6) and then discuss the case of D2 (ν = 0, J = 11).

From the literature, it is known that the behavior of the rotational quadrupole alignment parameter as a function of incidence energy can be related to features of the molecule–surface interaction at the preferred reaction site of the molecule for the initial rovibrational state considered.[31] For example, vibrationally excited H2 with a translational energy close to the threshold to reaction was found to prefer to react on a top site of Cu(100) due to features in the PES being more favorable; for instance, the increased lateness of the barrier at this site allowed more efficient conversion of energy from vibration to motion along the reaction path.[31,37,38] Next, the dependence of the rotational quadrupole alignment parameter on incidence energy of vibrationally excited H2 on Cu(100) could be explained on the basis of the anisotropy of the molecule–surface interaction energy at the top site. In our explanation of the behavior seen for H2 and D2 on Cu(211), we will therefore proceed in a similar manner.

Figure a shows the reaction density of D2 (ν = 1, J = 6) extracted from QCT calculations projected onto the Cu(211) unit cell. Here, we focus specifically on the D2 (ν = 1, J = 6) rovibrational state because it has been experimentally measured on Cu(111),[16] but the same mechanism appears to be present in our data for H2 (ν = 1, J > 2). All reacted trajectories up to a translational energy of 0.35 eV have been included. It is immediately clear that molecules in this particular state prefer to react on the t1 top site,[46] which, in the case of Cu(211), is at the step, with small outliers in reactivity pointing toward the bottom of the step. The t1 barrier is an extremely late barrier (r = 1.44 Å ), as can be seen in Table 3 of ref (46). The very late barrier allows for efficient conversion of vibrational energy to motion along the reaction coordinate.[10,11,40]

Figure 4

Three plots of the reaction density of D2 projected onto the Cu(211) unit cell. Panel (a) shows the reaction density of D2 (ν = 1, J = 6); all reacted trajectories up to a translational energy of 0.35 eV are included. Panel (b) shows the reaction density of D2 (ν = 0, J = 11); all reacted trajectories up to a translational energy of 0.35 eV are included. Panel (c) shows the reaction density of D2 (ν = 0, J = 2); all reacted trajectories up to a translational energy of 0.65 eV are included.

Figure shows a representative reactive trajectory of D2 (ν = 1, J = 6, m = 0) with a translational energy of 0.3 eV and plots the classical angular momentum, JC, as a function of the propagation time. JC is decreased before reaching the barrier, and a minimum in JC is reached at the transition state, where r becomes equal to 1.44 Å corresponding to the t1 barrier.[46] In the majority of the reacted trajectories, the minimum of JC is reached when r reaches the value of the t1 transition state, even when the molecule would make one or more bounces on the surface. This is a clear indication that rotational de-excitation takes place before the molecule reaches the transition state. This suggests that the reaction proceeds through rotational inelastic enhancement,[31] that is, the reaction is promoted by rotational energy flowing to the reaction coordinate. The bump in JC (i.e., its increase) still relatively far away from the surface is a feature that is also present in the majority of reactive trajectories. It is not completely clear to us what the cause is of this increase in JC still relatively far away from the surface before proceeding toward the transition state. We speculate that the increasing vicinity to the surface turns on the anisotropy of the molecule–surface interaction, thereby coupling rotational motion and stretching motion and providing a mechanism for the rotational energy to remain more constant while the bond extends and compresses due to the molecular vibration. This mechanism could consist in the classical angular momentum increasing when the bond extends to offset the effect of the bond extension on the rotational constant (upon bond extension the rotational constant decreases and, if not compensated, this would decrease the rotational energy). This could possibly explain the hump observed in JC at t ≈ 4500 atomic units of time in Figure .

Figure 5

Single representative reactive trajectory of D2 (ν = 1, J = 6, m = 0) with a translational energy of 0.3 eV. Blue shows the angular momentum (JC), red shows the bond length (r), and green shows the center of mass distance to the surface (Z).

There is also indirect evidence for rotationally enhanced reaction of D2 (ν = 1, J = 6, m = 0) in our QD calculations. Figure a shows inelastic scattering probabilities for D2 (ν = 1, J = 6, m = 0), and Figure b shows inelastic scattering probabilities for D2 (ν = 1, J = 6, m = 6). From a pairwise comparison of data with the same color between Figure a and Figure b, it is clear that D2 (ν = 1, J = 6, m = 0) has a considerably higher probability to rotationally de-excite in the scattering process compared to D2 (ν = 1, J = 6, m = 6). This suggests that the reaction of (ν = 1, J = 6, m = 0) is also rotationally enhanced in the quantum dynamics if the de-excitation occurs before the barrier is reached and the released rotational energy is transferred to motion along the reaction coordinate.

Figure 6

Rotational inelastic scattering probabilities for D2 for two different initial rovibrational states as a function of translational energy. Panel (a) shows rotational inelastic scattering probabilities for D2 (ν = 1, J = 6, m = 0). Panel (b) shows rotational inelastic scattering probabilities for D2 (ν = 1, J = 6, m = 6). Colors correspond to the final rovibrational state of the molecule, with blue being (ν = 1, J = 0), red being (ν = 1, J = 2), and green being (ν = 1, J = 4).

There are four possible mechanisms that affect the reaction probability and may affect the rotational quadrupole alignment parameters: two enhancing mechanisms and two steric hindering mechanisms.[31] Here, we have focused on one enhancement mechanism, inelastic rotational enhancement, since the evidence presented in Figures and 6a,b is consistent with this mechanism. Inelastic rotational enhancement requires the reaction to take place on a site with a low anisotropy in ϕ and a large anisotropy in θ at the barrier.[31] The main reasons for proposing the presence of this mechanism are the sharp downturn of the quadrupole alignment parameters for (ν = 1, J > 2) rovibrational states in Figure a,b and the rotational de-excitation seen in Figures and 6a,b. We note that inelastic rotational enhancement is the only mechanism that predicts a lowering of the rotational quadrupole alignment parameters.[31] A complete overview of the four mechanisms and what features of the PES they depend on can be found in Table III of ref (31).

A feature of the t1 site that facilitates the conversion of rotational energy to motion along the reaction coordinate is a low anisotropy of the potential in ϕ combined with a large anisotropy in θ. Figure shows the anisotropy at the t1 barrier[46] (r and Z are kept constant here), the top panel shows the anisotropy in ϕ, and the bottom panel shows the anisotropy in θ. It is clear that the anisotropy in θ is substantial, while the anisotropy in ϕ is very small compared to the anisotropy in θ. Somers et al.[31] have shown that the high anisotropy in θ may facilitate inelastic rotational enhancement. Inelastic rotational enhancement is expected to be most effective for low |m| states with J > 2, and the mechanism would lead to decreased rotational quadrupole alignment parameters.[31] The reason for the decrease in the rotational quadrupole alignment parameters is that m is approximately conserved, so that a decrease in J is possible only for a low |m|.

Figure 7

Anisotropy of the interaction potential at the t1 top site barrier[46] for θ (top panel) and ϕ (bottom panel).

It is also clear from Figure b that from the point of view of the orientational dependence of reaction, Cu(211) cannot be described as a combination of (100) steps and (111) terraces. The monotonic increase in the rotational quadrupole alignment parameter for D2 reacting on Cu(111)[16,32] is very similar to the behavior reported for Cu(100).[30,31] A slight downturn at translational energies close to the threshold to reaction has been reported in the case of Cu(100), indicating that the inelastic rotational enhancement mechanism is taking place. The downturn is however small and does not lead to negative quadrupole alignment parameters as we show here for H2 and D2 reacting on Cu(211). This is a clear indication that the reaction dynamics of the Cu(211) surface are distinct from the reaction dynamics of its component Cu(111) terraces and Cu(100) steps when looked at individually. This is most likely because the energetic corrugation of the Cu(211) surface is much lower compared to Cu(111) and Cu(100), a feature that favors the reaction of vibrationally excited molecules if sites with late barriers are present.

We now turn to an explanation for the monotonic decrease in the rotational quadrupole alignment parameter predicted for the (ν = 0, high J) states of H2 and D2 colliding with Cu(211) in Figure a,b. No downturn of the rotational quadrupole alignment parameter is observed for the (ν = 0) states even though D2 (ν = 0, J = 11, m = 11) reacts at the step as well as D2 (ν = 1, J = 6), as can be seen in Figure b. The lack of a downturn in the rotational quadrupole alignment parameter arises because the D2 (ν = 0) states react using a different mechanism. Figure shows a representative reactive trajectory of D2 (ν = 0, J = 11, m = 11), and it is clear that the angular momentum only drops after the transition state has been reached. This is a clear combination of elastic rotational enhancement for the helicopter molecules together with orientational hindering for the cartwheeling molecules, which causes the increase in the rotational quadrupole alignment parameters of the (ν = 0) molecules.[31] We note that D2 (ν = 0, J = 11) reacting on the step at the t1 site is due to the high initial rotational quantum number. The t1 barrier is slightly higher in energy than the lowest barrier to the reaction, but at this site, the reaction is less rotationally hindered if the molecule rotates in a plane parallel to the surface, and the barrier is much later than at the lowest b2 barrier on the terrace. This allows molecules in the vibrational ground state that are rotating fast in helicopter fashion and have incidence energies close to the threshold to reaction to react there by converting rotational energy to motion along the reaction path as the bond extends and the rotational constant of the molecule drops, while j remains roughly the same.

Figure 8

Single representative reactive trajectory of D2 (ν = 0, J = 11, m = 11) with a translational energy of 0.3 eV. Blue shows the angular momentum (JC), red shows the bond length (r), and green shows the center of mass distance to the surface (Z).

Above, we have shown that D2 in its (ν = 0, J = 11) and (ν = 1, J = 6) states prefers to react near the t1 site, that is, on or near the steps (see Figure a,b). This might seem to contradict an earlier conclusion that, at low incidence energies, D2 prefers to react on the terrace.[46] However, this conclusion was based on molecular beam experiments and simulations of those experiments, and under the conditions addressed,[46] the (ν = 0, J = 11) and (ν = 1, J = 6) states would hardly have population in them. A more appropriate picture of the reaction density for molecules under the conditions of ref (46) is shown in Figure c. There, it can be seen that D2 (ν = 0, J = 2) (this state would be highly populated in the beams used and simulated in ref (46)) prefers to react at the terrace b2 site, which has the lowest barrier to the reaction. The reaction density for D2 (ν = 0, J = 2) is in line with earlier findings that molecules in the vibrational ground state with low j react at the lowest barrier to the reaction[31,37,38] as well as with the findings for D2 + Cu(211) of ref (46).

Kaufmann et al.[45] did not measure rotational quadrupole alignment parameters in their recent study. We believe that the downturn of the rotational quadrupole alignment parameter at low incidence energies, which has not been observed before with this large downward shift for both H2 and D2 reacting on copper, may well be experimentally verified for both isotopes on Cu(211). Specifically, the reaction probability of H2 and D2 is large enough, and the (ν = 1, J = 4) rovibrational state of H2 and the (ν = 1, J = 6) state of D2 have large enough Boltzmann weights at reasonable surface temperatures (923 K) to make the downturn measurable. Comparing experimental rotational quadrupole alignment parameters to theoretical ones will provide a very stringent and detailed way of testing the accuracy of the electronic structure calculations used in the construction of the PES.

Comparing to Experimental E0(ν, J) Parameters

Next, we will make a direct comparison with the state-specific, or degeneracy averaged, reaction probabilities reported by Kaufmann et al.[45] From their experiments, they could derive dissociative adsorption probabilities by applying the principle of detailed balance to the measured time-of-flight distributions. However, comparing the relative saturation value of the reaction probability obtained from associative desorption experiments to the zero coverage absolute saturation values predicted by theory is not straightforward. The authors of the experimental paper pose several ways of scaling the experimental data in order to make a comparison to theoretical work possible. Scaling the experimental data to experimental molecular adsorption results introduces the uncertainties related to the direct molecular adsorption experiment used as a reference in this process. Theory calculates sticking probabilities in the zero coverage limit. When scaling the experimental desorption data to the experimental adsorption data, the zero coverage limit will only be a lower bound, especially when a molecular beam experiment with a very broad translational energy distribution is chosen as a reference.

We opt for the simplest and most direct method to scale to the relative experimental associative desorption data. In order to compare to the experimental E0(ν, J) parameters, where E0(ν, J) is the translational energy for which the reaction probability of the (ν, J) state is half of the maximum reaction probability measured for the (ν, J) state, we use the maximum translational energy sensitivity reported in Tables S7 and S9 of ref (45). Theoretical E0(ν, J) are taken to be the translational energy to which the reaction probability is half that of the reaction probability at the maximum translational energy for which the experiment is sensitive. This method also corresponds to what is showcased in Figure 13a of ref (45).

Figure shows E0(ν, J) parameters for H2 and D2 reacting on Cu(211). The agreement between the theory and experiment is excellent for H2. We calculated mean absolute and mean signed deviations between the experimental and theoretical E(ν, J) parameters (see Table ). It is clear from Figure and Table that the agreement between the theory and experiment is excellent in the case of H2, for which the total mean absolute deviation (MAD) (n–1∑ | E0, exp,n – E0,n | ) and mean signed deviation (MSD) (n–1∑0, exp,n – E0,n) values for QD and QCT calculations fall within chemical accuracy. We note that, for H2, the agreement is best for vibrationally excited molecules, while the reverse is true with respect to D2. For D2, the agreement is not yet within the chemically accuracy, mainly due to the slightly bigger discrepancies between the theory and experiment for the first vibrationally excited state. The theory, however, does not reproduce the rotational hindering that can be seen in the experimental data, that is, E0(ν, J) does not first increase with J until a maximum before falling off with increasing J. The theory shows no such behavior; here, the E0(ν, J) parameter falls off with increasing J for all methods investigated here.

Figure 9

E0(ν, J) parameters as a function of J for H2 and D2 reacting on Cu(211). Blue dots represent QCT results, red dots represent QD results, green dots represent MDEF results, and black dots represent experimental results.[45]

Table 3

Mean Absolute and Mean Signed Deviations for the Theoretical E0(ν, J) Parameters Compared to Experimental Values Shown in Figure

	MAD (eV) H2			MSD (eV) H2
parameters	total	ν = 0	ν = 1	total	ν = 0	ν = 1
QCT	0.0362	0.0384	0.0289	0.0209	0.0384	–0.0044
QD	0.0362	0.0449	0.0235	0.0241	0.0449	–0.006
MDEF	0.0509	0.0531	0.0272	0.0342	0.0532	0.0069
MDEF*	0.0239	0.0237	0.0306	0.0076	0.0236	–0.0157

Experiments on associative desorption of H2 from Cu(111)[18,45] and that of D2 from Cu(111)[15,32,45] likewise found the rotational hindering effect on reaction for low j. As for H2 and D2 interacting with Cu(211), we have not been able to reproduce this subtle effect in calculations on H2 + Cu(111)[34,35] and D2 + Cu(111)[32,34] in electronically adiabatic dynamics calculations. Here, we find that MDEF calculations on H2 and D2 + Cu(211) do not reproduce the trend either, suggesting that, in the previous calculations, neglecting electron–hole pair excitation was not the cause of the discrepancy between the theory and experiment. However, it is possible that calculations modeling electron–hole pair excitation with orbital-dependent friction (ODF) will succeed in recovering the subtle trend observed in experiments. For this, it may well be necessary that the ODF coefficients explicitly model the dependence of the tensor friction coefficients on the molecule’s orientation angles; earlier MDEF calculations on H2 + Cu(111) using ODF coefficients have not yet done this.[25]

According to Figure , the reactivity measured experimentally in the associative desorption experiments is, for most (ν, J) states, larger than that predicted theoretically, with the experimental E0(ν, J) being lower. With the use of the same scaling method to relate theory to experiment, Kaufmann et al.[45] obtained the same result for H2 and D2 reacting on Cu(111), and also in their case, they compared theory on the SRP48 functional.[32] To some extent, these results are odd, as calculations for H2 and D2 + Cu(111) using the original SRP functional showed that the theory overestimated the experimentally measured sticking coefficients.[35] However, also in this work, the theory generally underestimated the reactivity measured in associative desorption experiments.[35]

The paradox noted above may be explained on the basis of the BOSS model used in the calculations. This model neglects the effect of ehp excitation. Modeling this effect on sticking experiments should lower the theoretical reactivity, with computed sticking curves shifting to higher energies. Modeling the effect on associative desorption experiments should show the opposite effect if the modeling is done correctly, that is, starting with molecules being formed at the transition state and then desorbing.[35,80] The effect of ehp excitation in such calculations should lead to translational energy distributions of desorbed molecules being shifted to lower translational energies. The reaction probability curves obtainable from these distributions by assuming detailed balance (which, strictly speaking, is not applicable if ehp excitation is active) should then lead to computed reaction probability curves (E0(ν, J) values) shifted toward lower energies, in better agreement with the experiment (see Figure ).

The above also explains why our present MDEF calculations led to decreased agreement with the experiment: In these calculations, we modeled the associative desorption experiment as an initial-state selected dissociative chemisorption experiment, in which ehp excitation should have the opposite effect. If we assume the ehp excitation to have an effect that is similar in magnitude but opposite in sign with respect to the QCT calculations, then the net effect of modeling ehp excitation is to increase the agreement with the experiment to the extent that chemical accuracy is obtained for both (ν = 0) and (ν = 1) H2 on Cu(211). This is illustrated by the MDEF* mean absolute and mean signed deviations in Table . The MDEF* values have been calculated by subtracting the difference between the MDEF and QCT values from the QCT values. We finally note that we have assumed that the surface temperature does not much affect the measured E0(ν, J) through surface atom vibrational motion, which is in line with experiments,[23,24] as discussed in the Supporting information of Díaz et al.[35]

Classical Molecular Beam Simulations

One of the goals of this project was to carry out a molecular beam simulation using the QD method. Since surface atom motion and ehp excitations cannot be incorporated in QD calculations, we have also performed molecular beam simulations using the BOMD, QCT, and MDEF methods for D2 impinging on Cu(211) in order to quantify their effects on the reactivity measured in a molecular beam experiment. As discussed together with the comparison between our state-resolved reaction probabilities and the associative desorption experiments of Kaufmann et al.,[45] there are some effects on the reactivity from surface atom motion and ehp excitations though the effect falls within chemical accuracy. The molecular beam experiments we treat here were carried for a surface temperature of 120 K.[15,34]

In Figure , we compare BOMD calculations performed for a surface temperature of 120 K (red) to QCT (black) and MDEF (green) calculations carried out on our six-dimensional PES. As an additional validation of the PES, we have also calculated one energy point using the BOMD method with a rigid surface (blue). Each BOMD point is based on 500 trajectories, each QCT and MDEF point on 100,000 trajectories. The molecular beam parameters were taken from refs (15) and (35) and can be found in Table . From the excellent agreement in Figure between the black and blue data points at 80.1 kJ/mol, it is clear that our PES was accurately fitted, as was previously demonstrated in Figure S2 of ref (46). There we showed that for the dynamically relevant region of the PES (VMAX < 2 eV), the PES has an RMSE of <0.035 eV. Therefore, results obtained from QD calculations performed on our PES should not be influenced much by any (small) lingering inaccuracies still present in the PES related to the fitting procedure. It can also be observed from Figure that the effect of surface motion is small and well within the limits of chemical accuracy with respect to incidence energy. Due to the fact that H2 has a lower mass, we expect that the effect of including surface motion during the dynamics will be even less pronounced for H2 than D2. We should also note here that when low surface temperature experiments are considered, as with the 120 K surface temperature here, it is known from the literature that the BOSS model works well for activated H2 dissociation on metals.[26,33,34,36,60]

Figure 10

Reaction probability as function of the average translational energy for D2 on Cu(211), with molecular beam parameters taken from Table . BOMD results with a surface temperature of 120 K are shown in red, MDEF results are shown in green, and QCT results are shown in black. The blue point is a BOMD result for D2 on Cu(211) with a rigid surface.

It can also be seen from Figure that including the effect of ehp’s as a classical friction force shifts the reaction probability curve slightly to higher energies and that the effect is rather small and linear with respect to the average translational energy. From the literature, it is also known that including ehp excitations in the dynamics of H2 reacting on Cu(111) has only a marginal effect on the reaction probability.[25,32,36,81]

Due to the very small contribution of surface atom motion and nonadiabatic effects incorporated in the MDEF calculations to the overall reaction probability, we pose that H2 impinging on Cu(211) is an excellent system to fully simulate a molecular beam experiment using quantum dynamics methods since large discrepancies between the theory and experiment can reasonably be attributed to quantum effects during the dynamics, as the BOSS model should be quite accurate.

Quantum Molecular Beam Simulations

Figure shows results of simulations for four sets of molecular beam experiments, with varying molecular beam conditions. The experiment of Rendulic et al.[14] has the broadest translational energy distributions. The molecular beam parameters are taken from (the supporting information of) refs (15)(34), and (35). Here, theoretical results obtained for the H2 + Cu(211) system are compared to theoretical results for the H2 + Cu(111) system, where, for all theoretical results, the SRP48 density functional was used. We only make a comparison to the theoretical work since, to the best of our knowledge, there exists no published experimental molecular beam dissociative adsorption data for H2 reacting on Cu(211).

Figure 11

Comparison between four sets of molecular beam simulations for H2 + Cu(111) and Cu(211), using the SRP48 functional, and for normal incidence. Reactivity is shown as a function of average translational energy. The red dots correspond to QCT calculations for H2 + Cu(111). The green and blue dots correspond to the QCT and QD calculations for H2 + Cu(211), respectively.

In order to make the best possible comparison between the QCT and QD results, both results are calculated from initial state-resolved reaction probabilities for the same set of initial states. The molecular beam reaction probabilities predicted by QCT and QD calculations are in excellent agreement (Figure ). The excellent agreement holds for the very broad molecular beams of Rendulic et al. in Figure a as well as for the translationally narrow molecular beams of Auerbach et al.[15] shown in Figure b–d. However, QCT predicts slightly higher reaction probabilities, especially for the lowest translational energies. The consistently higher QCT reaction probability can be attributed to zero-point energy (ZPE) leakage, which is not possible by design in the QD calculations wherein the ZPE is preserved.

The excellent agreement between the QCT and QD calculations implies that, on the scale of a molecular beam experiment, in which a large number of rovibrational states are populated, quantum effects during the dynamics affect the reaction probability only in a very limited manner for reaction probabilities (>0.1%). The similarity between the QCT and QD calculations also holds over a wide range of molecular beam conditions, ranging from high to low incidence energies and from high to low nozzle temperatures.

From Figure , it is also clear that for most incidence energies (>22 kJ/mol), Cu(211) is predicted to be less reactive than Cu(111), as was reported previously for D2 + Cu(211).[46] The lower reactivity of Cu(211) compared to Cu(111) cannot be explained by the d-band model.[54,55] In our previous paper, we and others showed that the d-band model does make accurate predictions of the reactivity of different facets when similar reaction geometries are considered but the breakdown of the predictive prowess of the d-band model is caused by the geometric effect of the lowest barrier to the reaction of H2 dissociation on the low index Cu(111) surface not being on a top site.

Based on the results in Figure , we can now say definitively that, on the scale of a molecular beam experiment, neglecting quantum effects during the dynamics cannot be invoked to explain the lower reactivity of Cu(211) than Cu(111). This corroborates the theoretical results obtained in previous works,[44,46] where QCT calculations were performed for D2 and H2, and S0 values were measured for D2 + Cu(111) for E > 27 kJ/mol. More generally, we can state that molecular beam sticking of H2 on cold Cu(211) is well described with quasi-classical dynamics, and this very probably also holds for H2 reacting on Cu(111) and Cu(100).

Conclusions

In this work, we have carried out a comprehensive study of the quantum reaction dynamics of H2 reacting on the Cu(211) surface. A large number of TDWP calculations have been performed for all important individual rovibrational states reasonably populated in a molecular beam experiment. Our main conclusion is that the reaction of H2 (D2) with Cu(211) is well described classically. This is especially true when simulating molecular beam experiments where one averages over a large number of rovibrational states and molecular beam energy distributions.

We have however found that the extent to which the reaction depends on the alignment of H2 is somewhat dependent on whether QD or the QCT method is used, requiring careful validation of the dynamical model depending on the type of experiment that is being simulated. The QD method predicts stronger alignment effects on the reactivity than the QCT method for low-lying rotational states.

A comparison to recent associative desorption experiments suggest and BOMD calculations appear to show that the effect of surface atom motion and ehp’s on the reactivity falls within chemical accuracy, even for the high surface temperature used in the associative desorption experiments. We saw no evidence in our fully state-resolved data for the recently reported “slow” reaction channel, even though we carried out calculations over a translational energy range where this reported reactivity should be manifested. We speculate that the “slow” reaction channel is related to surface atom motion and its modeling requires the description of this motion, which is why we did not see it here.

In contrast to the theoretical and experimental results for D2 reacting on Cu(111) and Cu(100), at low translational energy, we observe a sharp downturn of the rotational quadrupole alignment parameters for vibrationally excited molecules. This downturn can be attributed to a site-specific reaction mechanism of inelastic rotational enhancement.