Experimental and theoretical study of rotationally inelastic diffraction of H₂(D₂) from methyl-terminated Si(111)

This study details the first rotationally inelastic diffraction (RID) of molecular hydrogen from a hybrid organic-semiconductor interface, bolstering the understanding of technologically relevant systems, such as fuel cells and biosensing electronics, where the interaction of H₂ with the surface charge density distribution of these materials is of paramount interest.^1–7 The nature of gas interactions at solid surfaces has been thoroughly examined through characterizations of chemisorption for a variety of interfaces, providing a strong basis for understanding the processes involved in surface chemical reactions.⁸ While a traditional route for understanding molecular chemisorption is monitoring the fraction of molecules that stick to a given surface,⁹ it has been theoretically demonstrated and experimentally proven that diffraction of molecules from a surface can provide complementary and precise information regarding the gas-surface interaction potential.^10,11 Diffraction patterns of diatomic molecules in particular have highlighted the role of rotational degrees of freedom in the chemisorption process, indicating a direct effect on dissociative probabilities and revealing fundamental details concerning the interaction of gases with surface charge densities.^8,12

Methyl-terminated Si(111) features a complete (1 × 1) methyl termination of its underlying lattice, endowing this surface with improved interfacial electronic properties and surface passivation, and establishing it as a model system for the understanding of organic-functionalized systems. This interface has been thoroughly characterized by many experimental and theoretical techniques: the surface structure and extent of methyl termination have been surveyed by elastic helium atom scattering,^13–15 scanning tunneling microscopy,^16,17 and synchrotron-based x-ray photoelectron spectroscopy (XPS),¹⁸ and the vibrational dynamics were studied via high-resolution electron energy loss spectroscopy,¹⁹ transmission infrared spectroscopy,^20,21 and inelastic helium atom scattering in conjunction with density functional perturbation theory.^14,15 This study employs rotationally inelastic diffraction and accompanying scattering calculations to further explore the surface charge density of CH₃–Si(111) and probe the interaction potential with H₂, broadening the understanding of the anisotropic features of this system. This potential anisotropy can serve a fundamental role in dissociative chemisorption for thin-film systems.

Rotationally inelastic diffraction (RID), whereby diatomic molecules impinging upon a surface convert translational energy into rotational energy (or vice versa) and are scattered into unique angular channels, was first reported in experiments involving diffraction of H₂ from MgO and LiF in the 1930s.^22–24 Since then, this technique has seen improvements through gains in angular resolution that have allowed a more precise investigation of the RID peaks.^25,26 High angular resolution and a wide range of incident energies have enabled the widespread use of RID, with the goal of investigating the gas-surface potential for a variety of systems.^{22,25,27–29} These studies have provided a wealth of information not only on the nature of interfacial dynamics but also on the theory that has been developed to study them.^22,23,30,31

The dependence of diatomic diffraction on molecular orientation processes makes theoretical modeling of the interaction a distinctly more complicated task than for monatomic systems such as He,³² because the rotational transitions that occur at the surface are highly sensitive not only to the anisotropy of the interactions but also to the corrugation of the gas-surface interaction potential and the coupling of these two factors.²³ However, a number of combined theoretical and experimental studies have already attempted to understand better the effect of parallel momentum transfer on elastic^33–35 and rotationally inelastic diffraction.^22,36 State-of-the-art theoretical models have shown some limitations, for example, in accurately reproducing the intensities of RID peaks relative to their elastic counterparts. However, they have proven very useful in reproducing general trends,^{11,30,31,37,38} which are the result of the main features characterizing the underlying potential energy surfaces (PESs), such as the corrugation and the anisotropy.

This paper presents the first rotationally inelastic diffraction measurements on an organic-functionalized semiconductor. In contrast to the various metal and alloy surfaces that have been characterized via RID, CH₃–Si(111) represents a new soft-film, i.e., low Debye temperature system for experimental studies. High-resolution rotationally inelastic diffraction of H₂ and D₂ has been employed to study the anisotropy of the CH₃–Si(111) surface via comparison with quantum dynamics simulations. The low-energy molecular diffraction measurements performed herein are primarily surface-sensitive, revealing information on the structure, charge density, and interfacial properties of the surface and its interactions with impinging molecules. This low incident energy permits a valid use of the rigid rotor assumption and energetically forbids vibrational excitations at the energies used.³² Rotationally inelastic diffraction spectra for H₂ and D₂ are compared, demonstrating much greater rotational excitation probabilities for D₂. Measurements of experimental RID excitations relative to elastic scattering events are examined as a function of incident angle and beam energy. The precision of these rotational probabilities is improved by employing the Debye-Waller model to account for the attenuation of diffraction peak intensities as a function of increasing surface temperature. Experimentally measured RID probabilities indicate a greater likelihood, not unexpectedly, of rotational excitation for higher beam energies, but show no significant dependence on incident angle. The potential energy surface (PES) of the H₂(D₂)/CH₃–Si(111) system has been modeled by interpolation of a density functional theory (DFT) energies data set, and used to study rotationally inelastic scattering and to simulate RID probabilities by means of quasi-classical and quantum dynamics. Quantum dynamics simulations have been used to assess the accuracy of the PES through a direct comparison with experimental results, whereas quasi-classical trajectories have been used to track down the aspects of the molecule-surface PES that lead to rotational excitation. These classical turning points show both large corrugation and high anisotropy for the interaction, as revealed by a thorough survey of the polar angular dependence of the PES landscape on these regions.

All chemicals were used as received. Water (≥18.2 MΩ cm) was obtained from a Barnstead Nanopure system. Czochralski-grown n-Si wafers (Virginia Semiconductor, Fredericksburg, VA), 381 ± 25 μm thick, were double-side polished, doped with phosphorus to a resistivity of 1 Ω cm, and oriented to within 0.1° of the (111) crystal plane. CH₃–Si(111) surfaces were prepared according to a published procedure.³⁹ The wafers were cut into 1 cm × 3 cm pieces and rinsed sequentially with water, methanol (≥99.8%, EMD), acetone (≥99.5%, BDH), methanol, and water. Organic contaminants were removed and the surfaces were oxidized by immersing the wafers in a freshly prepared piranha solution (1:3 v/v of 30% H₂O₂(aq) (EMD):18M H₂SO₄ (EMD)) at 90–95 °C for 10 min. The piranha solution was drained and the wafers were rinsed with copious amounts of water. Atomically flat H-terminated surfaces were prepared⁴⁰ by immersing the cleaned wafers in buffered HF(aq) (Transene Co., Inc., Danvers, MA) for 18 s, rinsing with water, and immediately placing the wafers in an Ar-purged solution of NH₄F(aq) (40%, Transene Co., Inc.) for 9 min. The wafers were agitated periodically to remove bubbles that formed on the surface. The Si samples were removed from the NH₄F(aq) solution, rinsed with water, and dried under a stream of N₂.

The Si wafers were chlorinated inside a N₂-purged glove box with <10 ppm O₂. A saturated solution of PCl₅ (99.998% metal basis, Alfa Aesar) in chlorobenzene (anhydrous, ≥99.8%, Sigma-Aldrich) was preheated with an initiating amount (<1 mg/ml) of benzoyl peroxide (≥98%, Sigma-Aldrich). The wafers were rinsed with chlorobenzene and reacted in the PCl₅ solution at 90 ± 2 °C for 45 min. The reaction solution was drained and the wafers were rinsed with copious amounts of chlorobenzene, followed by tetrahydrofuran (THF, anhydrous, inhibitor-free, ≥99.9%, Sigma-Aldrich).

The Cl-terminated surfaces were alkylated in a 1.0 M–3.0 M solution of CH₃MgCl (Sigma-Aldrich or Fisher Scientific) in THF. The reaction solution was heated to 50 ± 2 °C for >12 h. The reaction solution was drained and the wafers were rinsed with THF, submerged in THF, and removed from the N₂-purged glove box. The samples were sonicated sequentially for 10 min in each of THF, methanol, and water. The wafers were broken into 1 cm² pieces, dried under a stream of N₂, and sealed under Ar for shipment from Pasadena, CA to Chicago, IL. X-ray photoelectron spectroscopy (XPS) indicated that the surfaces were fully terminated with Si–C bonds and that there was no detectable surface oxidation.

To measure the elastic and rotationally inelastic diffraction peaks, these experiments required the use of an energy- and momentum-resolved ultra-high vacuum (UHV) scattering apparatus. This instrument has been described in full detail elsewhere;⁴¹ in brief, this apparatus consists of three primary sections: a differentially pumped beam-source manifold, an UHV crystal chamber, and a rotatable mass spectrometer detector. For beam generation, high pressures (800–2000 psi) of ultra-high purity gases (He, H₂, D₂) are expanded through a 15 μm diameter nozzle source that is cooled by a closed-cycle helium refrigerator to generate an intense and nearly monoenergetic (Δv/v ≤ 1% for He, 10% for H₂, D₂) supersonic neutral atomic beam. The beam energy is dependent on the nozzle temperature which can be adjusted for beam energies in the range of 35–90 meV. For diffraction and time-of-flight measurements, a mechanical chopper was used to modulate the beam prior to collision with a duty cycle of 50%. For collection of time-of-flight data, the beam is modulated with a pseudorandom chopping sequence for cross-correlation analysis.⁴² Collimation of the beam occurs through a series of apertures, resulting in a 4 mm spot size on the crystal (chopper-to-crystal distance = 0.4996 m) in the UHV scattering chamber (base pressure 3 × 10⁻¹⁰ Torr). The crystal was mounted onto a six-axis manipulator which can be positioned precisely to control the incident angle, θ_i, the azimuth, φ, and the tilt, χ, with respect to the scattering plane. Post-collision, the atoms travel along a 1.0234 m (crystal-ionizer distance) triply differentially pumped rotatable detector arm, with an angular resolution of 0.29° FWHM, after which they are ionized by electron bombardment. The ions then pass through a quadrupole filter where they are mass-selected before being collected by an electron multiplier. Angular distributions were obtained by scanning the detector at 0.1° computer-controlled increments while holding the incident angle at a fixed value.

The rotational experiments are carried out with n-H₂ or n-D₂ molecular beams, with varied stagnation pressures. The rotational populations of H₂ and D₂ were not measured directly, but instead were calculated based on previous theoretical and experimental results. Generally, the ratio of ortho- to para-H₂(D₂) is determined by the source temperature, and subsequent ortho–para conversion does not occur during the expansion.⁴³ The occupation of rotational states follows a near Boltzmann distribution within the respective ortho/para families, and can be characterized by an effective rotational temperature, T_R, which can be expressed with an empirical fit⁴⁴ for n-D₂ as

where the reference temperature, T_ref, is 298 K, T₀ is the beam temperature, and P₀d is the stagnation pressure times the nozzle diameter, given in units of Torr cm. The beam energies, rotational temperatures, and corresponding rotational populations for the experiments performed herein are given in Table I. In cases where high-translational-energy D₂ beams were created by seeding in a 1:1 mixture of H₂ and D₂, the rotational populations of the D₂ molecules were calculated using the measured nozzle temperature instead of the temperature derived from time-of-flight analysis, allowing for a more precise determination of these data points. All experiments reported herein used incident energies below 120 meV for D₂ such that over 99% of the impinging D₂ molecules are in their vibrational ground states. Additionally, vibrational excitations can be ignored due to the large inter-level spacing for D₂ molecules (∼380 meV).⁴⁵ 

Theoretical analysis of the H₂(D₂)/CH₃–Si(111) system has been performed within the Born-Oppenheimer static surface approximation (BOSSA). The Born-Oppenheimer approximation is justified by the different time scales associated with the motion of nuclei and electrons. The static surface approximation (SSA) is justified by the mass mismatch between the surface-terminating CH₃ layer and the H₂ and D₂ projectiles (although low recoil effects could be expected). Furthermore, experimental results have been extrapolated to a surface temperature of 0 K via the Debye-Waller correction, which allows a direct comparison between experimental measurements and SSA theoretical results. Working within this framework, a six-dimensional (6D) PES is first computed, and then is included in the nuclear Hamiltonian to perform dynamics simulations.

The 6D PES, for which the degrees of freedom are illustrated in Figure 1, has been computed by interpolation of a density functional theory energy data set. To perform the DFT periodic calculations, the plane-wave-based code VASP (Vienna ab initio Simulation Package)^46,47 has been used. In these calculations, the exchange-correlation energy of the electrons has been described using the generalized gradient approximation (GGA); in applying the GGA, the Perdew-Burke-Ernzerhof (PBE) functional⁴⁸ has been used. Additionally, the ion cores have been described using the PAW (projector augmented-wave) method.⁴⁹ 

The H₂(D₂)/CH₃–Si(111) system has been modeled using a five-layer slab and a 2 × 2 hexagonal surface unit cell, as shown in Figure 1. The size of the unit cell has been chosen to mitigate interaction between molecules in adjacent cells, which are present in the calculations due to the use of periodic boundary conditions (PBCs). To avoid artifacts caused by the use of PBC in the z direction, a vacuum layer of 21 Å has been set. A 7 × 7 × 1 Monkhorst-Pack grid of k-points was used to sample the Brillouin zone.⁵⁰ The cutoff energy for the plane-wave basis has been set to 650 eV. The lattice bulk parameter has been optimized, finding a value of a = 5.48 Å, which agrees well with the known experimental value (5.431 Å).⁵¹ The surface interlayer distances have also been relaxed until the forces were below a 1 meV/Å threshold.

To build the continuous 6D PES, V_6D(x, y, z, r, Θ, φ), representing the ground electronic state structure of the H₂(D₂)/CH₃–Si(111) system, this work makes use of the modified Shepard (MS) interpolation method, originally developed by Collins et al.^52,53 to study gas-phase reactions, and later adapted to study reactive scattering of molecules from surfaces.^54–57 Specifically, a recent implementation of the MS method,⁵⁸ which includes strict plane group symmetry and translational periodicity, is utilized herein. In the MS method, the interpolated PES is described by a weighted series of Taylor expansions centered on a number of DFT energy data points, N. These data points are sampled throughout the configuration space of the system. Thus, the global PES at any configuration $Q→$ is given by

where w is a weighting function, G_CNP and G_PG are the molecular permutation and plane symmetry groups, and S is the Taylor series expansion of the PES in the vicinity of the geometry of data point (g∘i), which denotes that the quantity for data point i has been transformed according to the symmetry operation g ∈ G_CNP × G_PG. Note that $Q→$ in Equation (2) does not represent the set of Cartesian coordinates used to compute the DFT energy points, but a set of redundant internal coordinates relative to them through the Wilson B matrix.⁵⁸ The Taylor series expansion used here is expressed as

where E(i) is the energy at the data point geometry i, ΔE(i) is the vector of first derivatives at data point i with respect to elements of ς(i), which are local coordinates resulting from a linear combination of the redundant internal coordinates $Q→$⁠, F(i) is the matrix of second derivatives at data point i with respect to elements of ς(i), and Δς(i) is the displacement of the point $Q→$ from the data point geometry $Q→$ in ς(i) coordinates.

It is important to note that the MS method uses a non-homogeneous sampling of the configuration space, so that more DFT energy points are used in the dynamically relevant regions. These regions are selected by using classical dynamics through a feedback process, hereafter called the GROW process. The first step begins with an initial basic version of the PES, defined in this case by only 50 DFT energy points. Then a small batch of classical trajectories is run on this basic PES, and from these trajectories new geometries are selected and added to the PES, thereby augmenting it. The new geometries are chosen according to two different criteria:⁵⁹ either new energy data points are added to the region most frequently visited by the trajectories or they are added to regions suspected to be the most inaccurate ones. Periodically, a larger batch of trajectories is run and used to compute some observables, which in this case are the rotational excitation probabilities. If the probabilities change significantly with the number of DFT data points added to the PES, the procedure goes back to the second step of running a small batch of trajectories. If the probabilities do not change significantly, the PES is considered converged, as illustrated in Figure 2. At this point, it should be pointed out that in order to properly sample the dynamical regions relevant to this analysis, the incident conditions of the classical trajectories are selected to correspond to the experimental conditions.

Quasi-classical dynamics—i.e., classical dynamics including the zero-point energy of the molecule—have been used to both grow and scrutinize the PES. To compute quasi-classical trajectories, the classical equations of motion are solved using the velocity-Verlet algorithm.⁶⁰ For each initial energy (E_i) and incident angle (θ_i), the classical scattering probability is calculated as an average over the molecular initial conditions, i.e., over the internal coordinates and conjugated momenta. The initial molecular conditions are sampled using a Monte Carlo method. To ensure low statistical error, approximately 2.5 × 10⁴ trajectories are computed for each set of initial conditions (E_i, θ_i). In these calculations, a molecule is considered reflected (and the integration ends) whenever the final distance between the molecule and the surface, z_f, becomes equal to initial distance, z_i, with the molecular velocity vector pointing towards the vacuum. To analyze rotational excitation upon scattering using classical dynamics, it must be taken into account that the classical angular momentum (L) follows a continuous distribution. Therefore, to analyze rotational excitations, the continuous representation is transformed into a discrete one.^61,62 This transformation is performed, in general, by evaluating the closest integer that satisfies

However, quantum selection rules for homonuclear diatomic molecules only allow rotational transitions for which Δ j = ± 2. Therefore, to obtain D₂ rotational excitation probabilities, the initial rotational state of the molecule is considered, and only the closest even or odd integers that satisfy Equation (4) are evaluated.

Diffraction probabilities are computed herein by solving the time-dependent Schrödinger (TDS) equation of the nuclear Hamiltonian using the multi-configuration time-dependent Hartree (MCTDH) method^63,64 which has already been successfully used to study diffraction of atomic projectiles,⁶⁵ as well as reactive scattering of molecular projectiles.^66–68 In the MCTDH method, the nuclear wavefunction is written as a sum of products of single-particle functions (SPFs). In this particular case, each SPF (χ) combines up to two degrees of freedom, and so the nuclear function of this system can be written as follows:

where Q represents the set of nuclear coordinates, and (N_x,y, N_z, N_Θ,φ) represents the number of SPFs used to describe each mode. The SPFs are in turn represented by linear combinations of time-independent primitive basis functions. In order to reduce the computational effort, and taking advantage of the lack of reactivity of D₂/CH₃–Si(111) in the energy range considered here, five-dimensional (5D) calculations have been performed, in which the atom-atom distance has been kept frozen at the equilibrium D₂ gas-phase distance (see Table II for calculation parameters). Within the MCTDH framework, the equations of motion for both the expansion coefficients and the SPFs are derived from the Dirac-Frenkel variational principle which leads to a set of coupled equations. In general, solving this coupled-equation system requires less computational effort than standard time-dependent wave packet propagation (TDWP) methods,⁶⁹ because the number of SPFs needed is smaller than the number of time-independent basis functions used in the standard TDWP methods. To obtain elastic and inelastic diffraction probabilities, a flux analysis is carried out with the aid of a complex absorbing potential located in the non-interaction z region.⁷⁰ 

Finally, it should be noted that to take full advantage of the MCTDH formalism, the multidimensional non-separable PES has to be rewritten as a linear combination of products of one- or two-dimensional functions. This transformation can be performed using the POTFIT algorithm, which is based on the approximation theorem of Schmidt,⁷¹ provided with the Heidelberg MCTDH package. Thus, the multi-dimensional PES is rewritten as follows:

where $ζj(u)$ is the jth one- or two-dimensional function used to describe the single particle mode u (the so-called natural potential). These functions are the ones used to expand the SPFs, which are represented in a primitive grid of points. The expansion coefficients c_jkl are determined by the overlap between the multi-dimensional PES and the natural potentials (see Table III for representative parameters for PES refits).

Angular distributions of H₂ scattered from CH₃–Si(111) are shown in Figure 3 along the high-symmetry $1̄21̄$ (⁠$Γ̄−M̄$⁠) and $011̄$ (⁠$Γ̄−K̄$⁠) azimuthal alignments; both diffraction spectra were taken under identical conditions, with a room temperature beam (E_i = 81.4 meV) and cold surface temperature (T_s = 140 K). The zeroth-order specular (θ_i = θ_f) diffraction peak and elastic first-order diffraction peaks are clearly resolved here; the relatively large first-order diffraction peak intensities compared to specular indicate a significant corrugation of the gas-surface potential. Elastic diffraction peaks arise when the kinematic condition for Bragg diffraction for in-plane scattering is satisfied such that

where $ΔK→$ is the change in the surface-parallel component of the H₂ wavevector $k→i$⁠, θ_i and θ_f are the initial and final scattering angles relative to the surface normal, and $G→mn$ is the surface reciprocal lattice vector. Information on the surface geometry and lattice constant can be determined by analyzing the spacing between diffraction peaks; the spacings between specular (⁠$ΔK→=0$⁠) and first-order diffraction peaks are $ΔK→=1.90Å−1$ and 3.29 Å⁻¹ along $1̄21̄$ and $011̄$⁠, respectively, corresponding to a hexagonally packed methyl adlayer with a real-space lattice constant of 3.82 Å. This value is in excellent agreement with previous helium diffraction measurements from this methyl-terminated Si(111) surface,¹³ with both identifying a (1 × 1) commensurate monolayer of methyl groups on Si(111).

Unlike with helium diffraction, when a H₂ or D₂ molecule scatters from the surface, it is capable of exchanging energy between its internal rotational and translational degrees of freedom; this exchange must still conserve total energy such that

Here, ħ is the reduced Planck constant, m is the mass of H₂ or D₂, E_i and E_f are the initial and final kinetic energies of the projectile, respectively, and $ΔEintrot$ is the energy exchanged between rotational and translational modes, which equals a difference between the molecule’s rotational energy levels as determined by the rigid rotor model; note that symmetry constraints within the hydrogen molecule impose a rotational selection rule of Δj = ± 2. This phenomenon of internal energy exchange results in RID peaks scattered at distinct angles from their parent elastic diffraction peaks such that

Figure 4(a) shows a diffraction spectrum for H₂ scattering from CH₃–Si(111), which includes a small peak associated with the specular j = 0 → 2 rotational excitation; note that RID peaks are labeled using the (j_i, j_f, m, n) notation, whereas elastic diffraction peaks are labeled with the (mn) notation. By positioning the rotatable detector arm at the final scattering angle of this RID peak, time-of-flight measurements can be used to confirm the inelastic scattering of the hydrogen molecules. A cross-correlation time-of-flight spectrum taken at the same conditions as this peak is shown in Figure 4(b); a diffuse elastic and rotationally inelastic hydrogen peak are observed to be fully resolved from one another. The flight time separation between the elastic and inelastic peaks can be used to calculate the energy exchange as a result of inelastic scattering from the surface (Figure 4(b), inset). The RID peak shown in Figure 4(b) arrived at a longer time than the elastic peak, indicating a loss in energy and velocity of the H₂ molecules; specifically, the calculated energy agrees well with the expected value of 45.4 meV for a j = 0 → 2 rotational excitation for an incident H₂ molecule.

In order to obtain highly resolved elastic and inelastic diffraction peaks in this range of relatively low beam energies, the surface probe was switched from H₂ to D₂. The predominant advantage in using D₂ is that the energy required to transition between rotational states is half that of H₂ (D₂: j = 0 → 2, 22.7 meV), as predicted by the rigid rotor model, leading to a greater probability of rotational excitation. In addition, the combination of nuclear spin states for each of the nuclei in H₂ and D₂ leads to the formation of ortho (symmetric) and para (antisymmetric) spin isomers, which pair with a set of rotational states to maintain the antisymmetry of the molecule as a whole. The degeneracy of these spin states causes only ∼25% of n-H₂ molecules to be in the j = 0 state, whereas ∼66% of n-D₂ molecules are in the rotational ground state.

The utility of D₂ as a surface probe is demonstrated in Figure 5(a), in which a diffraction spectrum for a room temperature D₂ beam from CH₃–Si(111) shows several high-intensity elastic and RID peaks. The resolution of inelastic peaks in D₂ spectra allows for precise measurements of rotationally inelastic transition probabilities; specifically, the rotational excitations can be measured as a function of scattering angle or beam energy, as seen in Figure 5(b).

As the setup of this instrument does not allow for direct measurement of the absolute incident beam flux, the probabilities of rotationally inelastic diffraction are evaluated as ratios of inelastic to elastic diffraction intensities.^26,28,32 This ratio enables the number density measured by the detector to be effectively converted into flux by accounting for the change in beam velocity that results from transferring energy between translational and rotational degrees of freedom. In addition, this approach obtains accurate values by accounting for the influences of instrumental broadening, finite crystal temperature, energy spread of the beam, and the initial rotational distributions of the incident beam.

Specifically, to relate a RID peak to its parent elastic peak, the following expression is used:

In this expression, I(j_i, j_f, m, n) and I(mn) are the peak area integrated intensities of D₂ molecules scattered from a crystal of finite surface temperature T_s. The square root term corrects for the velocity difference between elastic and rotationally inelastic scattering events, as discussed above.²⁸ Additionally, the rotational distributions of the impinging atoms are accounted for with n(j_i). By using this ratio rather than a pure probability, the effects of surface defects and beam geometry are eliminated and experimental error associated with evaluation of the Debye-Waller factor, W(T_s), is mitigated.

Attenuation of diffraction intensity due to thermal motion at the CH₃–Si(111) surface is accounted for via the Debye-Waller factor, which can be approximated in the semi-classical limit of a quantum-mechanical description of inelastic scattering as

for the specular peak, where D is the attractive well depth for a gas-surface interaction potential, M_eff is the effective surface mass that a given H₂ or D₂ molecule interacts with (assumed here to be 15 amu), k_B is Boltzmann’s constant, and Θ_D is the surface Debye temperature.^43,72 The attenuated intensities of RID peaks are corrected using the Debye-Waller factor associated with their parent diffraction peaks, which is a reasonable assumption based on the small parallel momentum transfer associated with the (0,2,0,0) and (0,2,0,$1̄$⁠) peaks.^61,73

The Debye-Waller factor can be quantified by relating the intensity of a peak to its ideal intensity for a lattice at 0 K, I₀, such that

As such, the natural log of I/I₀versusT_s produces a linear decay, from which the Debye-Waller factor for a given system can be extracted. Figure 6(a) shows the thermal attenuation of the specular peak at a given set of incident conditions over surface temperatures ranging from 140 to 350 K; the inset of this figure exhibits the linear decay which provides the Debye-Waller factor for the H₂/CH₃–Si(111) system. This exponential factor is derived from the normal and parallel momentum transfers during the scattering process (Δk_z and ΔK, respectively) and the associated mean-square displacements (MSD) of the crystal atoms (⁠$uz2andu∥2$⁠) such that

For the specular peak,

The perpendicular MSD and the well depth (D) of the gas-surface interaction potential can therefore be determined by plotting the derivative d2W(T_s)/dT_sversus cos²θ_i and extracting the slope and y-intercept, respectively.^13,74 Figure 6(b) shows diffraction decay rates at several angles, and the inset demonstrates the angular dependence of the Debye-Waller factor that provides a perpendicular MSD of (1.85 ± 0.30) × 10⁻⁵ Ų K⁻¹. Analogous data for the He/CH₃–Si(111) system are also shown in this figure, and its comparable slope provides a perpendicular MSD of (1.0 ± 0.1) × 10⁻⁵ Ų K⁻¹, just slightly below what is measured for the H₂ system.¹³ The difference in y-intercepts indicates a higher potential well depth for H₂ (D = 32 ± 9 meV) than for He (D = 7.5 ± 2.6 meV), as expected based on the higher degree of polarizability for H₂. Equation (11) uses these well depths to provide a surface Debye temperature of 723 K (503 cm⁻¹) for the H₂ system, which is considerably lower than 983 K (683 cm⁻¹) measured via He diffraction. While the surface Debye temperature measured for H₂/CH₃–Si(111) differs from He, both molecules seem to interact with a vibrational mode of the methyl adlayer: either Si–C bending (507 cm⁻¹) or Si–C stretching (683 cm⁻¹).⁷⁵ 

Figure 7 shows plots of r(0,2,0,0) and r(0,2,0,$1̄$⁠) for D₂ on CH₃–Si(111) as a function of beam energy, as calculated via Equation (10) with Debye-Waller corrections. There is a clear increase in rotational excitation probability with increasing beam energy, which is expected due to higher-energy incident molecules penetrating further into the surface charge density, thereby increasing the corrugation of the gas-surface interaction potential and the resultant torque on the non-spherical molecule. The dependence of rotational excitation probability on incident angle is weaker, with a more normal incident angle eliciting slightly more rotational probability for the first-order RID peak, and no apparent trend in angular dependence for the specular RID peak.

Quantum and classical dynamics simulations have been carried out with the goal of more accurately interpreting experimental measurements. First, to assess the accuracy of the theoretical tools employed herein, in particular the interpolated PES, quantum dynamics simulations have been compared with experimental measurements obtained at several representative sets of initial conditions (Figure 8). To perform this comparison, the quantum simulations consider the initial rotational distribution of the molecular beam (Table I). Overall, a good agreement is observed between both sets of data; in particular, the theoretical and experimental spectra exhibit an increase in the rotational excitation probability as the incident energy increases. These results indicate that the calculated PES is accurate enough to perform the required analysis.

In Figure 7, quantum ratios r(0,2,0,0) and r(0,2,0,$1̄$⁠) (dashed lines) are compared with the experimental ones; it should be noted that theoretical probabilities for a fixed incident energy were within 1° of the experimental scattering angle. The theoretical results generally agree with the increase in rotational excitation probability as a function of incident energy. Slight disagreement between theory and experiment is observed at higher beam energies, which is likely due to the strong corrugation of this system, as well as the use of a frozen-surface model for the PES, which becomes a less reasonable approximation when the energy of incident D₂ molecules nears the rotational barrier of the methyl group (∼100 meV for a surface temperature of 0 K).¹⁴ 

Having established the accuracy of the PES, a quasi-classical dynamics analysis can be performed on the regions of the PES that determine the characteristics of the diffraction spectra. To establish the validity of the classical analysis, quantum and classical total elastic and rotational excitation probabilities have been compared. Figure 9 displays the elastic and rotational excitation probabilities for H₂ (classical results only) and D₂ (quantum and classical results) as a function of the beam energy, along the $Γ̄−M̄$ direction, for three different incident angles. Several notable features are evident in this figure: (a) in the case of H₂, rotational excitation stays below 5% over the entire energy range investigated; (b) for D₂, quantum and classical simulations yield rather similar results—the classical rotational excitation probability fluctuates around 20%-25%, whereas the quantum probability is slightly higher (30%-35%). These results agree with the experimental results reported above, which show the presence of large RID in the diffraction spectra of D₂, whereas negligible RID has been observed for H₂ in diffraction spectra. This qualitative agreement between classical and quantum rotational excitation probabilities justifies further analysis of the systems using classical trajectories calculations, especially for incidence energies below 90 meV, beyond which point classic and quantum rotational excitation probabilities begin exhibiting different trends.

Classical trajectories have revealed that most molecules are scattered at a classical turning point around 2 Å from the plane formed by the H atoms in the methyl groups, with a molecular bond length around 0.78 Å. Given that these values are almost independent of θ_i and E_i within the range of experimental incident conditions, the characteristics of the potential have been analyzed at these (z, r) values. Figure 10 illustrates the one-dimensional (1D) potential energy profile along both the (⁠$Γ̄−M̄$⁠) and (⁠$Γ̄−K̄$⁠) azimuthal directions. One important feature of the PES that can be observed in Figure 10 is that the corrugation of the potential due to H atoms is very small compared to that due to C atoms, as the only indication of the presence of H atoms is the small shoulder seen between the hill and the valley. Another interesting feature that can be extracted is that the projectiles are guided towards the hollow and bridge sites. It is clear that molecules with incident energies below 90 meV cannot adiabatically approach the top site or surrounding sites; therefore most of the molecules are scattered from the hollow and bridge sites. This behavior is also corroborated by analysis of classical trajectories, which indicates that the molecules are scattered far from the top site after being efficiently steered towards the hollow and bridge sites.

Finally, Figure 10 also reveals a marked anisotropy of the potential in the classical turning point regions. For a closer inspection of the potential anisotropy, Figure 11 displays the relative potential energy as a function of the molecular orientation angle, Θ, at several z-distances from the surface. The anisotropy of the potential is observed to increase rapidly when the molecule approaches the surface, except in the case of the top site, which molecules do not sample, as discussed above. This rapid increase in the corrugation around the classical turning points is responsible for the substantial rotational excitation found in this system. It is noted that insubstantial differences are obtained for the trajectories followed by H₂ and D₂ molecules, when similar incident molecular velocities are considered. It can thus be concluded that the anisotropy experienced by both isotopes is similar, and therefore that the different rotational excitation observed is only due to the differences in rotational level spacing (∼45 meV for H₂, ∼22 meV for D₂).

Further insight into the nature of the PES can be extracted from Figure 12, which shows a series of xy-cuts for several (Θ, φ) orientations. In these plots, z and r values have been chosen according to the average values at the classical turning points revealed by a classical trajectory analysis, as discussed above. This figure indicates strong corrugation in the PES and shows that the molecular projectile feels, although rather weakly, the H atoms that belong to the CH₃ groups.

The diffraction of H₂ and D₂ molecular beams from CH₃–Si(111) was complemented by electronic structure and scattering calculations to investigate the nature of the interaction potential and surface charge density for a technologically relevant organic-functionalized semiconductor interface. Time-of-flight measurements confirmed the presence of rotationally inelastic diffraction for both H₂ and D₂ on this surface, and diffraction measurements demonstrated a stronger probability of rotational excitation for D₂ as compared with H₂, as expected based on the lower energy required for an internal exchange of rotational energy. The probabilities of these j = 0 → 2 rotational transitions were quantitatively evaluated as a function of beam energy and scattering angle, accounting for the thermal attenuation caused by incoherent motion at the CH₃–Si(111) surface by implementing the Debye-Waller model. The interaction potential between H₂ (D₂) and the surface charge density was then examined via the combination of these experimental results and quantum and quasi-classical dynamics simulations carried out on a continuous potential energy surface for H₂(D₂)/CH₃–Si(111) constructed by interpolation of density functional theory energies. Both experimental and theoretical data show high rotational excitation probabilities that increase with incident energy and which are weakly dependent upon incident angle, for the angle range investigated here. Additionally, dynamics calculations have identified the classical turning point regions (2 Å over hollow and bridge sites), and close scrutiny of these regions reveals a large anisotropy in the potential as a function of molecular orientation, which increases rapidly when the molecule approaches these regions and is responsible for the large rotational excitation observed experimentally. Overall, this work has revealed important details regarding the interaction of molecular hydrogen with a model hydrocarbon-decorated interface, which are important for fuel cells and next-generation energy systems.

S.J.S. acknowledges support from the Air Force Office of Scientific Research Grant No. FA9550-15-1-0428, and the Material Research Science and Engineering Center at the University of Chicago, No. NSF-DMR-14-20709. F.M., C.D., A.S.M., and M.d.C. acknowledge the MICINN Project No. FIS2013-42002-R, CAM Project No. S2014/MIT-2850, NANOFRONTMAG, and the CCC-UAM and the University of Chicago Research Computing Center for allocation of computer time. C.D. acknowledges the MICINN Project No. CTQ2013-50150-EXP and the Ramón y Cajal program. A.S.M. and M.d.C. acknowledge the FPI program of the MICINN, co-financed by the E.S.F. N.T.P. acknowledges a National Science Foundation Graduate Research Fellowship. N.S.L. acknowledges support from the National Science Foundation (Grant No. CHE-1214152), and research was in part carried out at the Molecular Materials Research Center of the Beckman Institute of the California Institute of Technology.