Applied Physics Report 95-4

Isotope effects in the Eley-Rideal dynamics of the recombinative desorption of hydrogen on a metal surface

Mats Persson

Department of Applied Physics, Chalmers University of Technology and Göteborg University, S-412 96, Göteborg, Sweden

Bret Jackson

Department of Chemistry, University of Massachusetts, Amherst, Massachusetts 01003, USA (February 9, 1995)

Abstract

The dynamics of the direct reaction between a gas phase H(D) atom and a D(H) atom adsorbed on Cu(111) is explored for the first time, using a fully three-dimensional flat-surface model. The reactive scattering problem is treated quantum mechanically by solving the time-dependent Schrödinger equation. Quasi-classical trajectory studies have also been performed. We find strong isotope effects in the product distributions that are not simply due to differences in the zero-point motion of the initially adsorbed atom.

In many cases the detailed mechanisms behind gas-surface reactions are most convincingly revealed by dynamical studies based on state specific investigations. An interesting surface reaction mechanism in this respect is the so-called Eley-Rideal (ER) mechanism, in which a gas phase species is assumed to react directly with an adsorbed species to form a molecule that immediately desorbs from the surface. The reaction dynamics can be probed both by preparing the initial state of the reactants and by interrogating the product state. This has been done in recent gas cell [1,2] and molecular beam experiments [3–7], providing strong evidence for the ER mechanism. Much of this evidence has been supported by theoretical dynamics studies based on fully quantum mechanical models [8–11]. An elegant way to probe the ER dynamics is to study the effects of isotopic substitution on the product state. This has been partly explored in our theoretical study [9] of the direct reaction between a gas phase H(D) atom and a (D)H atom adsorbed on a few model metal surfaces, and also in the molecular beam study by Rettner [5] of this reaction on a Cu(111) surface. Hydrogen on Cu(111) is an interesting system since it has been widely studied, and has become a prototypical system for the understanding of dissociation dynamics at surfaces [12]. Isotope effects have also been studied by kinetic measurements of the recombination of gas phase H(D) with D(H) adsorbed on Si(111) and Si(100), which have been suggested to follow an ER mechanism [13,14].

Isotope effects in ER dynamics can arise from changes in both the relative time-scales for motion of the reactants and in the zero-point energy of the initially adsorbed atom. In our first theoretical study, the isotope effects were identified as resulting from the former change, and were rationalized in terms of a simple impulsive collision model [9]. These results were in agreement with later measurements by Rettner [5] of isotope effects on the product translational energy distribution. However, these theoretical studies were based on a restricted collinear geometry involving two degrees of freedoom. A subsequent study [11] based on a fully three-dimensional flat-surface shows that near-collinear configurations give in many cases a vanishingly small contribution to the overall behaviour and stressed the key role played by the zero-point motion of the initially adsorbed atom on the reaction cross section and the internal state distributions. This latter study did not address the effects of isotopic substitution.

In this letter we explore for the first time the ER dynamics of the recombination reaction between a gas phase H(D) atom and a D(H) atom adsorbed on Cu(111), using a fully three-dimensional model based on a flat-surface approximation. The reactive scattering problem is treated quantum mechanically by solving the time-dependent Schrödinger equation. Quasi-classical trajectory calculations are also performed. We use a model potential energy surface (PES) that is consistent with available knowledge about the interaction of H and H_{2} with Cu(111). We find strong isotope effects on the HD product translational energy and rovibrational distributions that do no simply arise from the differences in the zero-point motion of the initially adsorbed atom.

Our previous paper, henceforth referred to as PJ, provided a detailed description of the model and method of implementation and the reader is referred to that work [11]. Only enough detail is provided here to make the letter reasonably self contained. For easy reference the initially adsorbed atom and incident atom will be referred to as the targon and the incidon, respectively. Moreover we will use the terminology H-on-D, D-on-H and so on to describe the different isotopic combinations of the incidon and targon.

The recombination of two atoms involves six degrees of freedom, assuming both electronic adiabaticity and that the substrate is rigid, as motivated by the relatively small masses of H and D. The positions of the targon with mass m_{t} and the incidon with mass m_{i} are given by the vectors, r_{t} and r_{i}, respectively. We usually transform these to a center-of-mass vector R = (m_{t}r_{t} + m_{i}r_{i})/M and a relative vector r = r_{i} − r_{t} with the associated masses: the total mass M = m_{t} + m_{i} and the reduced mass μ = m_{t}m_{i}/M. We note in passing that Z is dependent on the considered isotopic combination while M and μ are independent. The components of R are Z, the distance above the surface plane, and R_{∥}, the position over the surface plane. The surface geometry makes it useful to describe r in cylindrical coordinates ρ, ϕ and z, where r = ρ(cos ϕx̂ + sin ϕŷ) + zẑ.

The flat-surface approximation reduces the dimensionality of the problem by introducing three constants of motion. If the metal surface is flat, which is a reasonable approximation in view of the small diffusion barriers of H atoms on a metal surface, then the potential V is a function of only three variables: ρ, z and Z. The center of mass motion parallel to the surface is described by the momentum P_{∥}, which is conserved. The other constant of motion is J_{z}, the angular momentum conjugate to ϕ. We thus need only evolve ρ, z and Z (or ρ, z_{i} and z_{t}). The actual geometry of the reactive scattering problem is depicted in the inset of Fig. 1.

In this study we only consider the case where the incidon is at normal incidence. Thus, P_{∥} = J_{z} = 0, and we need only evolve a single wave function Ψ(ρ, z, Z). The general case of non-normal incidence is a bit more complex, and was fully explored in PJ for H-on-H. In cylindrical coordinates, the kinetic energy takes a simple form when one works with the weighted wave function Φ(ρ, z, Z) = ρ^{\frac{1}{2}} Ψ(ρ, z, Z), which evolves according to the effective Hamiltonian Ĥ^{eff} = K̂^{eff} + V̂^{eff} where,

K̂^{eff} = −\frac{h̄^{2}}{2μ} (\frac{∂^{2}}{∂ρ^{2}} + \frac{∂^{2}}{∂z^{2}} ) − \frac{h̄^{2}}{2M} \frac{∂^{2}}{∂Z^{2}}
V̂^{eff} = V (ρ, z_{i}, z_{t}) − \frac{h̄^{2}}{8μρ^{2}} ,  (1)

where we have expressed V in mass independent coordinates. The usual kinetic energy terms due to J_{z} and P_{∥} are zero at normal incidence. We would like to stress that this model gives a fully three-dimensional description of the free molecule.

The model is solved by propagating the initial wave function in time using well known pseudospectral techniques. The wavefunction Φ(ρ, z, Z) is represented by its values on a three-dimensional grid of points [15]. Initially, Φ(ρ, z, Z; t = 0) = ρ^{\frac{1}{2}} u_{0}(z_{t})G(z_{i})F (ρ) where u_{0}(z_{t}) is the ground state ν_{t} = 0 wave function for the targon-metal vibration, and G(z_{i}) is Gaussian wave packet centered at a position far above the surface, with a momentum determined by the incident energy E_{i} and an energy uncertainity ΔE. In this work we will only consider a single E_{i} = 75 meV. In order to satisfy proper boundary conditions, a smooth cut-off function F (ρ) has also been introduced in the initial wave function. The wavefunction is evolved in time numerically using a short time propagator algorithm based on the split-operator method. The scattered wave function is analyzed using a projection operator-correlation function approach.

In order to elucidate the role of quantum effects on the dynamics we have made a direct comparison with results from standard quasi-classical calculations. The coordinates ρ, z_{i}, z_{t}, and their conjugate momenta were evolved according to the classical Hamiltonian using the velocity Verlet algorithm. The reactive range of impact parameters, b = ρ(t → −∞), and the initial phase of the targon vibration, whose energy was set equal to its quantum-mechanical zero-point value, were sampled independently with a total of ∼ 10^{4} trajectories. ”Quantum numbers” j and ν were assigned to the product HD molecule by setting h̄j and h̄(ν+ \frac{1}{2}) equal to the classical angular momentum and the classical vibrational action variable, respectively. Average quantities were computed by averaging over all trajectories with proper weighting by the circular area element bdb. For instance, the quasi-classical reaction cross section is given by,

σ_{r} = 2π ∫_{0}^{∞} P_{r}(b)bdb,  (2)

where the opacity function, P_{r}(b), is determined by the relative number of reactive trajectories for each b.

An appropriate model PES describing the rearrangement of chemical bonds between the H and D atoms and a metal surface like Cu(111) is provided by a London-Eyring-Polanyi-Sato (LEPS) potential, V, of the form given in PJ. The two-body interactions describing the HD molecule and the H and D atoms adsorbed on Cu(111) are given by Morse potentials. The potential in the interaction region is modified by two Sato parameters. In this study we use the same potential parameters originally introduced in Ref. [9] and referred to as the PES-I model in PJ. This model gives a good description of the H-D interaction, with well-depth D_{m}=4.745 eV, and gives an H-metal vibrational frequency h̄ω ≈130 meV and an hydrogen chemisorption energy D_{a}=2.45 eV (excluding the zero-point energy h̄ω/2) in agreement with experimental and theoretical studies of H on Cu(111) [16,17]. The relatively large difference between the bond energies D_{m} and D_{a} results in a large exothermicity Q = D_{m}−D_{a}+E_{i}+h̄ω/2 ≈ 2.44 eV, where the isotopic differences in h̄ω/2 are negligible. The Sato parameters are choosen to give values for the height and position of the activation barrier to dissociative adsorption that are consistent with values suggested by studies of this part of the PES [12]. In fact our choosen values are in agreement with recent density functional calculations of the multi-dimensional PES for the dissociative adsorption of hydrogen [18]. Our construction results in an attractive PES which means that a large part of Q is released in the entrance channel. As demonstrated for the H-on-H combination in PJ, this results in a weak dependence on E_{i} and a highly vibrationally and rotationally excited H_{2} product. In the case of the HD product, Q allows a maximum possible rotational excitation j_{max} = 21 and vibrational excitation ν_{max} = 5.

The isotope effects observed in the results from our collinear model study were explained by a simple impulsive collision model [9]. The incidon was taken to be normally incident on the targon with zero impact parameter in a collinear arrangement as indicated in the inset of Fig. 1. For D-on-H, after the initial acceleration and subsequent collision the heavy incidon (D) remains close to or moving towards the surface while the light targon (H) recoils. For H-on-D on the other hand, the light incidon (H) reflects aways from the heavy targon (D) while the latter recoils. Because of the larger separation between H and D in the latter case by the time the targon is moving away from the surface, one expects a smaller reaction probability for H-on-D than for D-on-H. However, the molecules which do form for H-on-D should be more stretched and more vibrationally excited. These trends were identified in the quantum-mechanical results for the 2D collinear model. This argument is also consistent with Rettners study [5], in which he observed the H-on-D product molecules to be translationally colder and thus internally hotter, than for the D-on-H case. However, kinetic experiments comparing H-on-D and D-on-H rates suggest a weaker isotope dependence [13,14]. Also, in a real system only a vanishingly small number of particles are incident in a near-collinear arrangement. In fact, quasi-classical studies of H-on-H for the PES of interest in this work showed that most of the reactive trajectories had impact parameters, b, close to the H_{2} equilibrium bond length r_{0} ≈ 0.7Å. Another shortcoming of the collinear model is that it ingores the rotational motion of the product.

In Fig. 1 we show the calculated opacity functions, P_{r}(b), for H-on-D and D-on-H. These results looks qualitatively similiar with some exceptions. For both combinations there is a relatively narrow range of reactive b’s close to r_{0}. In contrast to D-on-H there are almost no reactive collisions for H-on-D when b is close to zero. This is consistent with the above lightheavy collision argument, but due to the area element in Eq. (2) these collisions contribute very little to the overall behaviour. Because the detailed structure of P_{r}(b) is quite different for the two isotopic combinations, one would expect that the underlying dynamics behind this difference should also affect the HD product distributions.

As can be seen from Table I, there are indeed strong isotope effects in the calculated reaction cross sections and averages of the different HD product distributions. The quasiclassical (QC) results are in good agreement with the quantum mechanical (QM) results for the average internal and translational energies, but are less accurate for σ_{r} and the averages ⟨ν⟩ and ⟨j⟩ of the HD rovibrational distributions. However, the trend of ⟨ν⟩ and ⟨j⟩ with isotope is reproduced. The H-on-D combination leads to lower reactivity than D-on-H, in agreement with the above impulsive collision arguement, but this argument fails with respect to the internal state distribution. The D-on-H combination gives larger ⟨E_{int}⟩ and ⟨ν⟩ than H-on-D. With regard to the magnitude the values for ⟨E_{Z}⟩ are in rough agreement with Rettner [5], who measures ⟨E_{Z}⟩ ≈ 1.1 and 0.85 eV, but the trend with isotope is opposite.

Strong isotope effects are also seen in the detailed final state distributions. In Fig. 2 we compare the QM and QC HD product translational energy distributions, P (E_{Z}), for the H-on-D and D-on-H combinations. The QM-QC agreement is fairly good, although, of course, the QC distribution lacks the detailed structure which reflects the quantization of the final ro-vibrational states. The QC distribution for H-on-D is much narrower than the QM one which is also seen in the rovibrational distributions as discussed below. Note that by total energy conservation, E_{Z} + E_{int} = Q, the probability distributions for E_{Z} and E_{int} display complementary behaviour.

In Fig. 3 we present the final vibrational state distributions, P (ν), for the H-on-D and D-on-H combinations. The isotope effects on P (ν) are quite pronunced. The D-on-H combination not only gives rise to a much larger ⟨ν⟩ but also a much broader distribution than H-on-D. In fact, the D-on-H combination gives an appreciable population up to ν_{max} = 5. These isotope effects are not simply a result of the corresponding differences in the zeropoint motion of the targon. This is shown by the QC results for P (ν) in Fig. 3, where results are presented using both the proper zero point energy of the targon, as well as a zero point energy corresponding to the mass of the other isotope. This modification of the QC zero-point motion leads to only small changes in P (ν) for both H-on-D and D-on-H. The minor broadening of P (ν) for H-on-D shows that differences between the QC and QM results is partly due to the approximate description of the zero-point motion in the QC treatment.

As shown in Fig. 4, the rotational state probability distributions, P (j), are also quite different between the two isotopic combinations. Both combinations give rise to a broad range of rotational excitations from j = 0 to j_{max} = 21. The H-on-D combination tends to populate higher j states than for the D-on-H combination. This behaviour is consistent with the variation of the QC reaction probability with impact parameter shown in Fig. 1. Once again we see that the QC distribution is too narrow for H-on-D but fairly good for D-on-H, and also that the aforementioned modification of the QC zero-point motion give rise to only minor changes in the QC P (j).

In conclusion, we have presented the first fully three-dimensional and quantummechanical study of isotope effects on the Eley-Rideal dynamics of hydrogen recombination on a metal surface. These effects are found to be strong and are not simply due to differences between the zero-point motions of the initially adsorbed H or D atom. We hope that this work will stimulate measurements of the product internal state distributions.

ACKNOWLEDGMENTS

We are grateful to Dianne Lynch for assisting us with the computations on the CM-5 at Thinking Machines Inc. The work by MP has been supported by the Swedish Natural Science Research Council (NFR) and in part by USA’s National Science Foundation under grants DMR 91-03466 and DMR 94-07055. Bret Jackson acknowledges support from the Division of Chemical Sciences, Office of Basic Energy Sciences, Office of Energy Research, US Department of Energy, under Grant # DE-FG02-87ER13744. Allocation of computer time by the Center for Parallel Computers (PDC) in Stockholm is also gratefully acknowledged.

REFERENCES

[1] R.I. Hall, I. Cadez, M. Landau, F. Pichou, and C. Schermann, Phys. Rev. Lett. 60 (1988) 337.

[2] P.J. Eenshustra, J.H.M. Bonnie, J. Los, and H.J. Hopmann, Phys. Rev. Lett. 60 (1988) 341.

[3] K.R. Lykke and B.D. Kay, in: Laser Photoionization and Desorption Surface Analysis Techniques, vol. 1208, ed. N.S. Nogar, (SPIE, Bellingham, WA, 1990) p. 18.

[4] E.W. Kuipers, A. Vardi, A. Danon, and A. Amirav, Phys. Rev. Lett. 66 (1991) 116.

[5] C.T. Rettner, Phys. Rev. Lett. 69 (1992) 383.

[6] C.T. Rettner and D.J. Auerbach, Science 263 (1994) 365.

[7] C.T. Rettner, J. Chem. Phys. 101 (1994) 1529.

[8] P. Kratzer and W. Brenig, Surf. Sci. 254 (1991) 275.

[9] B. Jackson and M. Persson, J. Chem. Phys. 96 (1992) 2378; Surf. Sci. 269/270 (1992) 195.

[10] B. Jackson, M. Persson, and B.D. Kay, J. Chem. Phys. 100 (1994) 7687.

[11] M. Persson and B. Jackson, J. Chem. Phys. 102 (1995) 1078.

[12] See the contribution by H.A. Michelsen, C.T. Rettner, and D.J. Auerbach and also the separate contributions by B.E. Hayden and S. Holloway, respectively, in: Surface Reactions, ed. R.J. Madix, (Springer, Berlin 1993)

[13] D. D. Koleske, S. M. Gates, and J. A. Schultz, J. Chem. Phys. 99 (1993) 5619.

[14] D. D. Koleske, S. M. Gates, and B. Jackson, J. Chem. Phys. 101 (1994) 3301.

[15] Because H-on-D and D-on-H involve larger masses than H-on-H, smaller grid spacings and more gridpoints (256^{3}) are used than in PJ.

[16] E.M. McCash, S.F Parker, J. Pritchard, and M.A. Chesters, Surf. Sci. 215 (1989) 363.

[17] K. Gundersen, B. Hammer, K.W. Jacobsen, J.K. Nørskov, J.S. Lin, and V. Milman, Surf. Sci. 285 (1993) 27.

[18] B. Hammer, M. Scheffler, K.W. Jacobsen, and J.K. Nørskov, Phys. Rev. Lett. 73 (1994) 1400.

FIGURES

FIG. 1. Quasi-classical reaction probabilities, P_{r}(b), as a function of impact parameter, b, for H-on-D and D-on-H. Both panels correspond to E_{i}=75meV, ν_{t} = 0 and normal incidence. The inset in the upper panel depicts the reactive scattering geometry in the flat-surface model together with the coordinates of the targon (t) and the incidon (i). The dotted line, ρ = 0, indicates the geometry used in the collinear model.

FIG. 2. Quantum mechanical (solid lines) and quasi-classical (dashed ) translational energy distributions, P (E_{Z}), for the H-on-D and D-on-H combinations. The initial conditions are the same as in Fig. 1 and ΔE = 50meV.

FIG. 3. Quantum mechanical (open symbols) and quasi-classical (dashed and dotted lines) vibrational distributions, P (ν) for the H-on-D and D-on-H combinations. The dashed line is the result using the proper zero-point energy of the targon, while the dotted line is the result using a zero point energy corresponding to the other isotope. The initial conditions are the same as in Fig. 2.

FIG. 4. Quantum mechanical (open symbols) and quasi-classical (dashed and dotted lines) rotational distributions, P (j), for the H-on-D and D-on-H combinations. The dashed line is the result using the proper zero-point energy of the targon, while the dotted line is the result using a zero point energy corresponding to the other isotope. The initial conditions are the same as in Fig. 2.

TABLES

TABLE I. Quantum mechanical reaction cross sections σ_{r}, average translational and internal energies, ⟨E_{int}⟩ and ⟨E_{Z}⟩, and average vibrational quantum number, ⟨ν⟩, and rotational quantum number, ⟨j⟩, for H-on-D and D-on-H. Results are for normal incidence with E_{i} = 0.075 meV and ν_{t} = 0. The quasi-classical results are within the parenthesis.

H-on-D  D-on-H
σ_{r}(Å^{2})  0.27(0.31)  0.39(0.30)
⟨E_{Z}⟩(eV)  0.95(0.94)  0.77(0.75)
⟨E_{int}⟩(eV)  1.48(1.48)  1.67(1.69)
< ν >  1.05 (0.51)  2.13 (1.89)
< j >  11.2 (14.5)  8.8 (11.7)