Stereodynamics effects in grazing-incidence fast-molecule diffraction

Stereodynamics effects in grazing-incidence fastmolecule diffraction M. del Cueto,aA. S. Muzas,bF. Martín,cde and C. Díazf Grazing-incidence fast-projectile diffraction has been proposed both as a complement and an alternative to thermal-energy projectile scattering, which explains the interest that this technique has received in recent years, especially in the case of atomic projectiles. On the other hand, despite the richer physics involved, molecular projectiles have received much less attention. In this work, we present a theoretical study of grazing-incidence fast-molecule diffraction of H2from KCl(001) using a six-dimensional density functional theory based potential energy surface and a time-dependent wavepacket propagation method. The analysis of the computed diffraction patterns as a function of the molecular alignment, and their comparison with the available experimental data, where the initial distribution of rotational states in the molecule is not known, reveals a puzzling stereodynamics effect of the diffracted projectiles: Diffracted molecules aligned perpendicular, or quasi perpendicular, to the surface reproduce rather well the experimental diffraction pattern, whereas those molecules aligned parallel to or titled with respect to the surface do not behave as in the experiments. These results call for more detailed investigations of the molecular beam generation process. 1 Introduction Grazing-incidence fast-atom (GIFAD) and fast-molecule (GIFMD) diffraction, measured experimentally for the first time in 20071,2, can be considered as a good alternative to thermal-energy atom (TEAS) and molecule scattering (TEMS), proposed and widely developed between the 70’s and 90’s (see Refs.3–6 and references therein). The physical mechanism behind GIFAD and GIFMD, already described in pioneering theoretical works7–9, is the effective strong decoupling between the fast motion parallel to the surface, which does not induce diffraction, and the slow motion perpendicular to it, responsible for the observed diffraction peaks. As already discussed in Refs.7,10, under GIFAD conditions the projectile feels a periodic potential along the incidence direction (the xaxis in our case) because the potential barely varies over the interval δz=dtanΩ,dbeing the lattice parameter along the incidence direction and Ωthe grazing angle. Consequently, the projectile feels a periodic potential as long as d<< aDepartment of Chemistry, University of Liverpool, Liverpool, L69 3BK, UK. bCentro de Física de Materiales CFM/MPC (CSIC-UPV/EHU), 20018 Donotia-San Sebastián, Spain. cDepartamento de Química Módulo 13, Universidad Autónoma de Madrid,28049 Madrid, Spain. dCondensed Matter Physics Center (IFIMAC), Universidad Autónoma de Madrid, 28049 Madrid, Spain. eInstituto Madrileño de Estudios Avanzado en Nanociencia (IMDEA-Nanociencia), Cantoblanco 28049, Madrid, Spain. fDepartamento de Química Física, Facultad de CC. Químicas, Universidad Complutense de Madrid, 28040 Madrid, Spain; E-mail: [email protected] E[tanΩ(∂V/∂z)−1]. A periodic potential will result in a parallel momentum change ∆K|| x=0, because ∆K|| x=−1/vxRdx(∂V/∂z) = −1/vx[V(d,y)−V(0,y)] = 0. Therefore, the only diffraction peaks observed under GIFAD and GIFMD conditions are the ones perpendicular to the incidence direction. Thus, in GIFAD, the incidence energy inducing diffraction can be easily varied by changing the total energy (further accelerating the molecular beams) and/or changing the incidence angle -theoretically it has been shown that, if phonons are neglected, GIFAD conditions are preserved up to 10◦11–13, although experimentally the maximum angle has been limited to approximately 2◦so far. At this point, it is also important to remark that even if the decoupling between the parallel and the perpendicular motions was not fully met, GIFAD measurements could still be performed14. Thus, GIFAD allows one to record diffraction for a much wider energy range than in TEAS, i.e., to access regions of the potential energy surfaces (PESs) unreachable to TEAS15. Furthermore, quantum decoherence effects,mainly due electronic excitations16, are not a more important problem than in TEAS17, although phonons may also have some influence in diffraction patterns18,19. GIFAD has been widely used since 2007 (see15,20 and Refs. therein). Using this technique, it has been possible to analyze from relative simple surfaces such as clean insulators1,2,21–26, semiconductors27–33, and metal surfaces34,35, to more complex systems such as superstructures absorbed on a metal substrate36,37, reconstructed oxide surfaces14,38, graphene grown on 6H-SiC(0001)39, a monolayer of silica adsorbed on Mo40,41, or organic monolayers adsorbed on metal sufarces42,43. GIJournal Name, [year], [vol.], 1–11 | 1 FAD has also shown its ability to follow phase transitions in real time at an organic-inorganic interface44. This vast experimental effort has been accompanied by the subsequent development and implementation of theoretical tools, allowing for a more precise interpretation of the experimental measurements. All these tools are based on the Born-Oppenheimer approximation (BOA), where the PESs can be computed using pairwiselike analytical functions18,19,30,39,45–48 or density functional theory (DFT)21,25,26,33,49–55. PESs computed at multiconfiguration self-consistent field level, considering a cluster model, have also been used1,17. To perform the dynamics calculations, several approaches have been adopted, from classical49,51,52 to quantum1,12,17,26,30,47,55 methods, as well as semi-classical28,36,41,56 and semi-quantum18,21,25,50,53,57–59 methods. H H r q f (XCM YCM ZCM) (z1y1z1) (z2y2z2) Top K Top Cl Bridge K-K Bridge K-Cl Z X Y Fig. 1 In the left panel, we represent the irreducible KCl(0001) unit cell (grey triangle) and the high symmetry geometries site used in the CRP interpolation (black dots). In the right panel, we represent the molecular and atomic DOFs. GIFMD, on the other hand, has received much less attention, despite the fact that, as in the case of TEAS60,61, the H2molecule: (i) is as easy to generate as atomic H (a widely used projectile in GIFAD experiments); (ii) is lighter than He (another widely used projectile in GIFAD), which would further reduce surfacephonons inelastic processes; and (iii) can reveal aspects of the surface landscape that may be relevant in other contexts due to the internal degrees of freedom (DOFs) and, in the case of the ionic surfaces, to the interaction of its quadrupole moment with the electric field created by the ionic crystal, which is very sensitive to the surface details. Experimentally, GIFMD measurements have been performed, for example, for LiF(001)2,62, silica/Mo(112)41, and alanine/Cu(110)43, but the subsequent analysis has been much less thorough than in the case of GIFAD. The first GIFMD theoretical studies already pointed out interesting potentialities of this technique. For example in Refs.10,63, it was shown that 1reflectivity probabilities resulting from GIFMD, for several H2/metal surface systems, were surprisingly similar to dissociative adsorption probabilities obtained in TEAS, which implies that GIFMD could be used to estimate dissociative adsorption probabilities at thermal and quasi-thermal energies, provided that the dissociative adsorption probability is approximately equal to 1 minus the reflectivity. Thus, GIFMD would make the dissociative adsorption saturation limit experimentally accessible. Here, it is worth noticing that the atoms generated upon molecular dissociation in GIFAD do not get adsorbed on the surface due to their large parallel energy, but they are scattered. Thus experimentally by measuring the atomic reflectivity probabilities one could infer the dissociative adsorption probabilities. Focusing on diffracted molecules, our previous studies for H2/LiF(001)64,65 have revealed a strong dependence of the diffraction patterns, and diffractograms, on both the initial rotational state of the molecules (Jand MJ) and the crystallographic incidence direction. A similar dependence has also been found for normal incidence at low impact energy66. On the other hand, the initial vibrational state of the molecule was shown to play a negligible role, as well as, rotationally inelastic processes. These findings have been rationalized in terms of the interaction between the quadrupole moment of the molecule and the electric field associated with the ionic crystal, which differs from one crystallographic direction to another65. To verify that the conclusions drawn for H2/LiF(001) are valid for other surfaces with a marked ionic character, in this work, we have studied GIFMD for H2/KCl(001) focusing on the role played by the initial rotational state (J, MJ) of the molecule. The comparison between our simulated diffraction patterns and the experimental ones available in the literature67 reveals a puzzling stereodynamics effect. Diffracted molecules aligned perpendicular (or quasi perpendicular) to the surface reproduce fairly well the experimental pattern, whereas molecules with other alignments do not. 2 Methodology Our theoretical study is based on two main pillars: 1. The inclusion of the six nuclear degrees of freedom of the molecules, freezing the surface ones. 2. The validity of the BOA, which allows us to divide our calculations in two steps: (a) Electronic structure calculations based on DFT. (b) Six-dimensional quantum dynamics simulations. 2.1 Potential energy surface To build the required continuous six-dimensional (6D) PES, we have applied a modified version of the corrugation reducing procedure (CRP)68 to a set of energies computed at DFT theory level for 4116 differents configurations of H2over KCl(001). The main idea behind the CRP is that the 6D-PES can be divided in a 6D smooth function, which is relatively easy to interpolate with high accuracy, and a corrugated function that represents most of the corrugation coming from the atom-surface corrugation. In this way, by removing the corrugated function from the full 6D PES, one can accurately interpolate the remaining 6D smooth function, and then add back the corrugated function to obtain the whole PES. Using the modified-CRP version described in detail in Ref.65, the 6D PES is written as: V6D(R) = I6D(R)+V3D(r1)Lz0,δz(r1)+V3D(r2)Lz0,δz(r2),(1) 2 | 1–11 Journal Name, [year], [vol.], Site Top K Top Cl Bridge K-K Bridge K-Cl (θ,φ)1(0,0) (0,0) (0,0) (0,0) (θ,φ)2(π 2,0) (π 2,0) (π 2,0) (π 2,0) (θ,φ)3(π 2,π 4) (π 2,π 4) (π 2,π 4) (θ,φ)4(π 2,π 2) (θ,φ)5(π 2,3π 10 ) (θ,φ)6(π 2,4π 5) (θ,φ)7(π 4,0) (π 4,0) (θ,φ)8(π 4,π 4) (π 4,π 4) (θ,φ)9(3π 4,0) (θ,φ)10 (3π 4,)(3π 10 ) Table 1 (XCM,YCM,θ,φ) molecular configuration used in the interpolation process. See Fig. 1 for (XCM ,YCM) site description. For each of this configuration 196 DFT single point energies are computed by varying ZCM (from 0.25 to 7.0 Å) and r (from 0.4 to 2.3 Å). where I6D(R)is a 6D smooth function (Rrepresenting the molecular DOFs, XCM,YCM,ZCM, r, θ,ϕ-see Fig. 1 ) that can be easily interpolated. To evaluate this function, we use a cubic spline interpolation for each of the 2D cuts (ZCM,r) of the DFT data set, and subsequently, we perform a symmetry adapted Fourier interpolation along θand ϕfor each high symmetry site (see Fig. 1). Eventually, we carry out a Fourier interpolation in (XCM ,YCM) and ϕbased on the previous ones. V3D(r)represents the interaction between one atom of the molecule and the surface, a highly repulsive part of the 6D potential. Here rirepresent the atomic coordinates (see Fig. 1). Finally Lz0,δz, written as: Lz0,δz(ri) = h1+exphzi−z0 δzii−1,(2) is a logistic control function that allows one to control the amount of CRP correction applied to the obtained I6D. We have taken the H/KCl(001) 3D-PES from Ref.55 as a basis for our corrugated function (see eq. 1). But one may also use any other mathematical function whose subtraction from the 6D PES leads to a smooth function easy to interpolate. The 4116 6D-DFT single point energies, grouped in 21 (XCM,YCM,θ,φ) molecular configurations (see Fig. 1 and Tab.1), have been computed using the plane-wave based code VASP69–71. In performing DFT calculations, we have relied on the generalized gradient approximation (GGA) through the PBE functional72. The cutoff energy for the plane wave expansion was set to 550 eV, the projector augmented wave (PAW) method73,74 has been used to take into account the interaction of the core electrons with nuclei, and a 3×3×1 kpoints grid was used to sample the Brilloiun zone. To describe the molecule/surface system with periodic boundary conditions, we have used a 5-layer slab and a 2×2 surface unit cell. To avoid interactions of the projectile with its periodic images and with the top periodic image of the surface, we have placed a vacuum layer of 20 Å between the slabs in the Z direction. (see Ref.55 for further details). Fig. 2 displays several 2D cuts of the computed 6D PES. One can see that the characteristics of the PES are consistent with the nature of the interaction between H2and an ionic surface. Thus, the cartwheel configuration is energetically more stable than the helicopter one over Cl−top, because the negative charge located on the Cl−ions interacts favorably with the excess of positive charge located in the hydrogen atoms, whereas the opposite occurs over K+top, because the positive charge located on the K ions interacts more favorably with the excess of negative charge located in the H-H bond. We also observe that K+top exhibits the largest anisotropy, i.e., the highest variation of energy with respect to molecular alignment, whereas the smaller anisotropy corresponds to K-Cl bridge. As expected, these findings are, generally speaking, similar to those observed for H2/LiF(001)65,66,75,76. 2.2 Time-dependent wave packet propagation (TDWP) method Once we have a good representation of the PES, we solve the time-dependent Schrödinger equation, ˆ HΦ(R,r;t) = i∂Φ(R,r;t) ∂t,(3) to obtain diffraction probabilities as a function of the incident conditions (Ei,θi,φi) and the initial rovibrational state of the molecule (vi,J,MJ). In eq. 3, ˆ His written (in atomic units) as: ˆ H=−1 2Mh∂2 ∂X2 CM +∂2 ∂Y2 CM +∂2 ∂Z2 CM i−1 2µ ∂2 ∂r2+ˆ J2 2µr2+V6D,(4) where Mand µrepresent the total and reduced mass of the molecule, respectively, and ˆ Jthe rotational operator, whose eigenfunctions are the spherical harmonics YJMJ(θ,φ). To solve eq. 3, we have used a TDWP method77 as implemented in78. In this implementation, widely used to study molecule-surface interactions at thermal and quasi-thermal energies79–81, the wave packet is propagated according to eq. 3 using the split-operator method82. A direct-product discrete-variable representation (DVR), with constant grid spacing, is used to represent the dependence of the wave function on XCM,YCM,ZCM, and r83. To represent the dependence on θand φ, we used a nondirect-product finite-basis representation (FBR) of spherical harmonics. Gauss-associated-Legendre and Fourier transforms are employed to transform the wave function from FBR to DVR, and vice versa84. The initial wave function is written as: Φ0(XCM,YCM,ZCM,r,θ,φ) = ϕv,J(r)YJ,mJ(θ,φ) ×1 √AeiK0RZdkZb(kZ)1 √2πeikZZ, i.e., as a product of a rovibrational wave function describing the molecular initial state, a plane wave describing the motion parallel to the surface, and a Gaussian wave packet describing the motion perpendicular to it. The scattered wave function is absorbed by a second-order complex absorbing potential (CAP) placed in the asymptotic region, and it is analyzed using the Balint-Kurti formalism66,85. The relevant parameters used in the calculations are summarized in Tab. 2. Note that, as a result of this analysis, the obtained diffraction peaks are delta functions. Thus, to take into account the Journal Name, [year], [vol.], 1–11 | 3 1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0 0.6 0.7 0.8 0.9 1 0.6 0.7 0.8 0.9 1 K_top q=0 Cl_top q=0 0.6 0.7 0.8 0.9 1 1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0 0.6 0.7 0.8 0.9 1 0.6 0.7 0.8 0.9 1 eVeV 0.6 0.7 0.8 0.9 1 6.0 5.5 5.0 4.5 4.0 3.5 3.0 2.5 2.0 0.6 0.7 0.8 0.9 1 0.6 0.7 0.8 0.9 1 6.0 5.5 5.0 4.5 4.0 3.5 3.0 2.5 2.0 r(A) o Z (A) o abcd efgh q=0º q=0º q=0º q=0º q=90º q=90º q=90º q=90º Fig. 2 2D cuts through the H2/KCl(001) 6D PES. Top panels: cartwheel configurations; bottom panels: helicopter configurations. From left to right: top K, top Cl, K-K bridge, and K-Cl bridge. The spacing between the contour levels is 0.1 eV. The dotted lines show the 0.1 eV isovalue. Blue balls represent Cl atoms, green balls K atoms, and black balls H atoms. experimental resolution, our theoretical results have been convoluted using a 2D Gaussian function with typical values of the widths σΘ=0.4 deg. and σλ⊥=0.01 Å. At this point, it should be noticed that although the experimental diffraction chart corresponds to a collection of data recorded on the Laue circle, at different normal energies (different perpendicular wavelengths), phonons and electronic excitations inelastic processes, as well as beam collimating conditions may induce a spread of the recorded data. None of these effects are considered in our calculations. It is also important to clarify that when comparing with experimental measurements, we are just comparing rotationally elastic results. Relative to the theory-experiment comparison, it is also worth mentioning that our simulations reveal that: (i) as in the case of H2/LiF(001)64, rotational excitation upon molecular scattering is a minority process, totally negligible for excitation into J=9 (see Appendix); (ii) the only populated diffraction peaks (elastic and inelastic) are those perpendicular to the incidence direction, in the reciprocal space, as one would expect if the decoupling between the normal and the parallel motions is fulfilled. 3 Results and discussion We have carried out dynamics calculations for a wide range of initial energies aiming to match the range considered in the experiments67. For that purpose, we have run up to 11 dynamic calculations for each set (v,J,MJ) defining the initial internal state of the molecule. All these calculations comprise normal energy values, E⊥, between 50-600 meV. The corresponding parallel energy, E||, is chosen so that the incidence angle with respect to the surface is kept constant at Ω=2◦, i.e., E|| =E⊥ tan2Ω. Although Initial wave packet Initial position, Z020 Width, ∆Z00.23-0.64 Grid parameters Z minimum value 1.5 Grid points, NZ224 Specular grid points, Nsp Z480 Grid spacing, ∆Z0.10 r minimum value 0.6 Grid points, Nr40 Grid spacing, ∆r0.15 Grid points, NX, NY24, 24 Maximum Jin rotational basis 14-15 Time propagation Time step, ∆t 2 Total propagation in Z(r) 15000-65000 CAP in Z(Zsp) Initial value, Zmin (Zsp min) 13.50 (26.50) CAP length 10.30 (22.90) CAP in r Initial value, rmin 3.60 CAP length 2.85 Other parameters Analysis values, Z∞13.50 Table 2 Parameters used in the TDWP calculations. All values are given in atomic units. Details about the grid parameters can be found in Ref.78. CAP stands for complex absorbing potential 4 | 1–11 Journal Name, [year], [vol.], Q Y Fig. 3 Schematic representation of the diffraction process of H2from a KCl(001) surface, with the projectile path represented in red. We indicate in the detection plane both the defection angle Θ, and the diffraction order, n. The inset shows a top view of the surface, indicating the correspondence between Φand the crystallographic direction. we are using a grazing angle larger than the one used in the experiments (Ωexp =0.7◦-see Ref64), this is not relevant here because our results reveal that for H2GIFMD from KCl(001) only E⊥matters. At this point, it worth noting that in our simulations we are not imposing any decoupling between fast parallel and slow normal motions. This decoupling, which the reason why only diffraction peaks perpendicular to the incidence direction are observed, comes out naturally from our full dimensional calculations. Thus, we can make (λ⊥,Θ) plots for each set of values (v, J,MJ). Here, λ⊥=h/√2mE⊥,mbeing the mass of the molecule, hthe Planck constant, and Θthe deflection angle, defined as the angle between the components parallel and perpendicular of the molecule wave vector. Θis related to the diffraction order, n, by the expression: dsinΘ=nλ⊥, where dis the channel periodicity in incidence direction. A schematic representation of the diffraction process is shown in Fig 3. From this figure we can see that tanΘ=tanΨ/tanΩ Due to the computational effort required to perform this GIFMD simulations, we have only performed calculations along the [100] direction (Φ=45◦), because, in this case, the experimentally measured diffraction chart shows a characteristic structure with peaks that appear and vanish as a function of the perpendicular energy (perpendicular λ)67, thus facilitating, in the absence of experimental measured diffractograms, a detailed comparison with theoretical results. This is not the case for other directions considered in existing experiments, where all possible diffraction peaks are presented in the charts, and no modulation is observed. Note that the lack of available experimental diffractograms prevents us from a detailed comparison based on relative intensities. It is also important to point out here that, to carry out a direct comparison with experimental diffraction patterns, one should know the rovibrational distribution in the molecular beam. However, contrary to TEMS experiments86–89, determining such rovibrational distributions is a real challenge in GIFMD Deflection angle (deg) Wavelength (Å) 0.9 0.8 0.7 0.6 0.5 0.4 -40 -20 0 20 40 (a) 164 Chapter 7. Surface analysis with diffraction techniques: H2/KCl(001) effect we have studied is that of the initial vibrational state vi. To that end, we have run simulations for (vi=0, Ji=0)and (vi=1, Ji=0), at the [100] and [110] incidence directions. In Figures 7.6A and 7.6B, we represent results obtained with a PBE-PES, along the [100] and [110] crystallographic incidence directions, respectively. From this Figure, we can observe that the diffraction results for the vibrationally excited projectiles (bottom panels) are equivalent to the one obtained for the vibrational ground state (top panels), with just minor changes in the relative intensities and minor wavelength shifts. This result agrees with our previous study of H2diffraction on LiF(001)[54], where we observed a negligible dependence on the initial vibrational state. Thus, we have studied the effect of the initial rotational (Ji,mJi) states on the diffraction patterns, keeping always vi=0. FIGURE 7.6: Theoretical H2/KCl(001)diffraction patterns as a(Q,l? i)function, for the [100] (Panel A) and [110] (Panel B) incidence directions. Top-row shows (vi=0, Ji=0)results, and bottom-row shows (vi=0, Ji=0)results. 164 Chapter 7. Surface analysis with diffraction techniques: H2/KCl(001) effect we have studied is that of the initial vibrational state vi. To that end, we have run simulations for (vi=0, Ji=0)and (vi=1, Ji=0), at the [100] and [110] incidence directions. In Figures 7.6A and 7.6B, we represent results obtained with a PBE-PES, along the [100] and [110] crystallographic incidence directions, respectively. From this Figure, we can observe that the diffraction results for the vibrationally excited projectiles (bottom panels) are equivalent to the one obtained for the vibrational ground state (top panels), with just minor changes in the relative intensities and minor wavelength shifts. This result agrees with our previous study of H2diffraction on LiF(001)[54], where we observed a negligible dependence on the initial vibrational state. Thus, we have studied the effect of the initial rotational (Ji,mJi) states on the diffraction patterns, keeping always vi=0. FIGURE 7.6: Theoretical H2/KCl(001)diffraction patterns as a(Q,l? i)function, for the [100] (Panel A) and [110] (Panel B) incidence directions. Top-row shows (vi=0, Ji=0)results, and bottom-row shows (vi=0, Ji=0)results. -40 -20 0 20 40 Fig. 4 Simulated H2/KCl(001) GIFMD patterns along the crystallographic direction [100] for the initial rovibratonal state: left panel (v=0,J=0); right panel (v=1,J=0). experiments and, in fact, it has not been achieved yet. Therefore, any comparison with existing experiments can only be done at a qualitative level. Fig. 4 shows the simulated diffraction patterns (λ⊥,Θ) for (v=0,J=0) and (v=1,J=0). As can be seen, the diffraction patterns barely depend on the initial vvalue, so that one can neglect the vibrational DOF in our ensuing analysis. Nevertheless a closer look at the diffraction patterns reveals that vibrationally excited molecules lead in general to slightly larger intensities for lower diffraction orders than molecules in the vibrational ground state. A similar phenomenon was observed for H2/LiF(001)64,65. The role of the molecule rotational DOFs is more complex. In Fig. 5, we compare the experimental diffraction spectra available in the literature67 with our simulated ones for several initial rotational states, J, of the molecule. To obtain these results, we have computed diffraction probabilities for all allowed MJvalues and subsequently obtained the MJ-averaged results. Fig 5 shows that the simulated diffraction spectrum for J=0 resembles the experimental one, whereas the simulated spectra for J>0 are rather different. One may argue that this result suggests that the molecules in the experimental beam are mainly in their rotational ground state, i.e., they have an extremely low rotational temperature. However, in view of the beam generation technique used in experiment15,20, neutralization of fast H+ 2molecular ions by alkali atoms, which is likely to produce hot H2molecules, this seems to be unlikely. To further analyze the role of the molecular rotation in GIFMD, we have simulated the diffraction spectra (λ⊥,Θ) as a function of Jand MJ(see Fig. 6). Fig. 6 shows that our simulated diffraction patterns reproduce rather well the experimental ones, for molecules in a cartwheel or quasi-cartwheel configuration, i.e., with MJclose to 0, whereas molecules with MJclose to J (helicopter) show diffraction patterns clearly different from the experimental ones. At this point, it should be noticed that the high order diffraction peaks present in the diffraction charts are a consequence of the strong corrugation felt by both helicopter and cartwheel molecules. This strong corrugation is the result of the local interaction between the H atoms (slightly positively Journal Name, [year], [vol.], 1–11 | 5 Fig. 5 H2/KCl(001) GIFMD patterns along the crystallographic direction [100]. (a) experimental pattern (data taken from ref.67); (b), (c), (d), and (e) simulated pattern for J=0,J=1,J=2, and J=9, respectively. charged) and the molecular bond (slightly negatively charged) with the surface anions and cations (see fig.7). On the other hand, as we discuss below, the substantial differences found in diffraction charts for MJ=0 and MJ=1 are due to the dissimilar corrugation features felt by the molecule aligned parallel and perpendicular to the surface. In summary, the results shown in Fig. 6 confirm the main role that molecular alignment may play in GIFMD as already suggested in Ref.67. To understand the stereodynamics effect revealed in Fig. 6, we have analyzed the surface corrugation felt by cartwheel and helicopter molecules at their theoretical classical turning points distribution throughout the KCl(001) unit cell in the limit of a sudden collision with fixed molecule alignment and normal energy scaling. The 2D and 3D representations of the classical turning points, shown in Fig. 7, give us a measure of the corrugation felt by the molecule as a function of its alignment. One can see that cartwheel and helicopter molecules feel a quite different corrugation. Cartwheel (MJ=0) molecules see little differences between K+and Cl−top sites, i.e., the classical turning points over these sites are quite similar, so that the corrugation observed is the result of a double periodicity (two types of diffracted trajectories). In contrast, for helicopter molecules the classical turning point over the K+top site is much higher than over the Cl−top site, so that the observed corrugation is the result of a single periodicity (one type of diffracted trajectories). Using corrugation arguments, we can qualitatively interpret the results shown in Fig. 6 in terms of supernumerary rainbows, which appear as modulations in the relative intensity of the Bragg peaks (see Fig. 2-11 of Ref.20). The interference between trajectories sharing the same phase A and A’ produces diffraction patterns with all diffraction Bragg peaks present, whereas the interference between trajectories with different phases, A and B, but same deflection angle produces quantum supernumerary rainbows. The convolution of supernumerary rainbows and Bragg peaks patterns produces the modulation of the intensity of the Bragg peaks. The resulting modulation is closely related to the corrugation felt by the projectile and may lead to the disappearance of some diffraction peaks, as shown in Fig.6 for cartwheel molecules. We can also rationalize the results shown in Fig. 6 in terms of the geometrical structure factor (Sg), which explains the modulation of the intensity of the diffraction peaks due to the presence of two more types of atoms in a crystal90. In the case of KCl, Sg=fK(G)eiGdK+fCleiGdCl ,(5) where Gand dK/Cl represent the reciprocal lattice and the atomic basis vectors, respectively. The form factors fidepend on G, but also on the detailed characteristics of the surface-projectile electronic interaction. As shown in Fig. 7. the corrugation felt by cartwheel and helicopter molecules is quite different, and therefore, the form factors and the geometrical structure factor are quite different, which explains the remarkable differences observed in the diffraction spectra as a function of the MJvalue. One important question at this point is where this apparent stereodynamics selectivity comes from. Several phenomena could be behind this effect: • The H2generation mechanism15,20 might impose the creation of molecules in a cartwheel configuration as a results of selection rules governing the formation of H+ 2. Confirming or rejecting this hypothesis would require detailed experimental and theoretical studies of the neutralization reaction, which so far has not been sufficiently investigated. • Beyond Debye-Waller effects, which alone are not able to reconcile theoretical results obtained at T=0K with experimental TEMS results at T>0K (see91), one could speculate that the coupling between the vibrational modes along the [100] crystallographic direction (where K and Cl atoms alternate) and the rotational motion of the molecule may induce a change in the molecule alignment preferentially toward the cartwheel configuration, which cannot be described within the frozen surface approximation. To assess this speculation further experimental and theoretical analysis would be required. From the theoretical point of view, to include the effect of phonons in the dynamics is not an easy task. In the case of H2scattering at low energy, sevendimensional quantum dynamics, including the perpendicu6 | 1–11 Journal Name, [year], [vol.], 0.9 0.8 0.7 0.6 0.5 0.4 0.9 0.8 0.7 0.6 0.5 0.4 Wavelength (Å) -50 -30 0 30 50 -50 -30 0 30 50 -50 -30 0 30 50 -50 -30 0 30 50 -50 -30 0 30 50 Deflection angle (deg.) Fig. 6 H2/KCl(001) GIFMD patterns along the crystallographic direction [100]. (a) Experimental pattern (data taken from ref.67); (b), (c), (d), (e), (f), (g), (h), (i) and (j) simulated pattern for (J=0,MJ=0), (J=2), MJ=0), (J=2,MJ=1), (J=2,MJ=2, (J=9,MJ=0), (J=9,MJ=1), (J=9, MJ=3), (J=9,MJ=6) and (J=9,MJ=9), respectively. lar motion of the subsurface layer atoms, and the phonon sudden approximation have been tested to describe vibrational excitation and state-to-state scattering probabilities of H2interacting with a metal surface92. However, the performance of these approaches to deal with insulating surfaces and diffraction phenomena is still unclear. In the case of fast grazing incidence, phonons have been included in the dynamics to study GIFAD from insulating surfaces using a semiquantum phonon-surface initial value representation18,19. Results for He and Ne GIFAD from LiF(001) reveal that thermal lattice vibrations can affect the relative intensity in the diffraction pattern and even the interference maxima. However, at this point it should be reminded, on the one hand, that H2is lighter than He and Ne atoms and, on the other hand, that 6D quantum dynamics simulations, based on the surface frozen approximation, were able to describe fairly well diffraction patterns of H2GIFMD from LiF(001) in comparison with experiment. In summary, further theoretical and experimental studies are required to elucidate the role of phonons on H2GIFMD from KCl(001). Finally, one might suggest that shortcomings of the DFT functional used to built the 6D-PES used in this study could also be blamed for the disagreement between theory and experiment. One may wonder if the GGA-PBE functional could yield better results for cartwheel than for helicopter aligned molecules, but such an effect has never been observed for molecule-surface interactions at low energy, where the effect of DFT-functional on the dynamics has been widely investigated (see93 and reference therein). Same argument may hold for van der Waals (vdW) effects. The inclusion of vdW effects in GIFAD has been shown to modify the relative intensity53 and even the position of the minima and maxima of interference in the diffraction patterns30,53,55, which could explain the difference between the experimental diffraction pattern shown in Fig. 6 (a) and our simulated patterns for cartwheel, and quasi-cartwheel, aligned molecules (see Figs. 6 (b), (c), (f) ,and (g)). However, vdW can hardly explain the disappearance of the structure in the diffraction patterns of helicopter or quasi-helicopter aligned molecules. It should also be remembered that similar DFT parameters were used to build the H2/LiF(001) PES used to accurately reproduce experimental diffraction pattern64,65. 4 Conclusion We have studied diffraction of H2from KCl(001) under fast grazing incidence as a function of the internal molecular degrees of freedom. The comparison between our results and the available experimental diffraction patterns reveals a striking stereodynamics effect. Diffracted molecules aligned perpendicular, or quasi perpendicular, to the surface reproduce rather well the experimental diffraction pattern, whereas those molecules aligned parallel or tilted to the surface fail completely in reproducing the experimental observations. This stereodynamics effect has not been observed previously in the few studies found in the literature dealing with molecular projectiles. The source of this phenomenon is not totally clear, although possible causes are the Journal Name, [year], [vol.], 1–11 | 7 [100] (a) MJ=0 E⊥=300 meV [100] (c) MJ=0 E⊥=64 meV [100] (b) MJ=J E⊥=300 meV [100] (d) MJ=J E⊥=64 meV [100] [100] [100] [100] 8 6 4 2 0 YCM(! 𝐴) 8 6 4 2 0 8 6 4 2 0 8 6 4 2 0 YCM(! 𝐴) YCM(! 𝐴)YCM(! 𝐴) 0 2 4 6 8 XCM(! 𝐴) XCM(! 𝐴)XCM(! 𝐴) XCM(! 𝐴) 0 2 4 6 8 0 2 4 6 8 0 2 4 6 8 [100] ZCM(! 𝐴) Fig. 7 2D and 3D (insets) representations of the classical turning points (Z(Å)) for cartwheel (MJ=0) and helicopter (MJ=J) molecules for two perpendicular kinetic energies (E⊥). The squares delimited by white dashed lines represent the KCl(001) unit cell -see Fig. 1. Note that E⊥= 64 and 300 meV correspond to λ⊥=0.8 and 0.37 Å, respectively. experimental molecular beam generation process and moleculephonons interactions. To elucidate which of these sources, if any, is the main responsible for the observed behavior, further experimental and theoretical studies are required. We hope that the results presented here will further motivate the experimental groups working in the field to improve the current techniques of generation of molecular beams with the aim of controlling the rotational state (J,MJ) of the molecules, at the same level that is already achieved in TEAS86–88. Conflicts of interest There are no conflicts to declare. Acknowledgements The authors are grateful to Prof. L. Mendez for enlightening discussions about the beam generation process, to Prof. H. Winter and Dr. E. Meyer for useful discussions about their experimental results, and to Prof. G.-J. Kroes, Dr. E. Pijper and Dr. M. F. Somers for allowing us to used their quantum dynamics code. This work has been supported by the MICINN projects PID2019105458RB-I00 and PID2019-106732GB-I00, ’Severo Ochoa’ Programme for Center of Excelence in R&D (CEX2020-001039-S), ’María de Maeztu’ Programme for Units of Excellence in R&D (CEX2018-000805-M), and ANPCyT project PICT-2016 2750. We acknowledge the allocation of computer time by the Red Española de Supercomputación and the Centro de Computación Científica at the Universidad Autónoma de Madrid (CCC-UAM). M. del Cueto and A. S. Muzas acknowledge the FPI program of the MICINN co-financed by the European Social Fund. Appendix: Rovibrational excitation In figures 8 and 9, we compare rotationally elastic diffraction probabilities with the most populated rotational excited or deexcited diffraction channels for two diffraction peaks (0,0) and (¯ 2,2) and four initial rovibrational states of the molecule (v=0, J2, mJ=0,2) (Fig. 8) and (v=0, J=9, mJ=0,9) (Fig. 9). From Fig. 8, we can see that for J=2, rotationally inelastic probabilities are one order of magnitude smaller than the rotationally elastic peaks, except for the higher normal energies and lower elastic probabilities, where elastic and inelastic probabilities are of the same order of magnitude. At this point, it is worthy to remark that the presence in the beams of the molecules with low rotational states should be testimonial. On the other hand, for J=9 (see Fig. 9) rotational deexcited probabilities J=9 →J=7 (the most populated rotationally inelastic channel) are between two and three orders of magnitude smaller than rotationally elastic peaks. These results are similar to those obtained for H2/LiF(001)64. Relative to rotational excitation is also worth pointing out that only rotationally inelastic diffraction peaks (RID’s) perpendicular to the incidence direction are significantly populated. In the case of the crystallographic direction analyzed here, h100i, only (¯n,n) and (n,¯n) RID’s are populated. This can be observed from the data in Tab. 3, where we show raw elastic and inelastic diffraction probabilities obtained for the second order diffraction peaks, for a normal energy equal to 400 meV. Similar results are obtain for other diffraction orders. 8 | 1–11 Journal Name, [year], [vol.], (v=0, J=2, mJ=0) →(v=0, J=2, mJ=0) (v=0, J=2, mJ=0) →(v=0, J=0, mJ=0) (v=0, J=2, mJ=0) →(v0, J=4, mJ=0) Peak Probability (¯ 2, 0) 0.363220 10−10 (¯ 1,¯ 1) 0.171869 10−10 (¯ 1, 1) 0.548651 10−01 (0, ¯ 2) 0.294560 10−10 (0, 2) 0.000000 10+00 (1, ¯ 1) 0.548639 10−01 (1, 1) 0.000000 10+00 (2, 0) 0.000000 10+00 Peak Probability (¯ 2, 0) 0.264077 10−10 (¯ 1,¯ 1) 0.161417 10−11 (¯ 1, 1) 0.488435 10−03 (0, ¯ 2) 0.227826 10−10 (0, 2) 0.000000 10+00 (1, ¯ 1) 0.488627 10−03 (1, 1) 0.000000 10+00 (2, 0) 0.000000 10+00 Peak Probability (¯ 2, 0) 0.223183 10−10 (¯ 1,¯ 1) 0.119949 10−10 (¯ 1, 1) 0.215639 10−02 (0, ¯ 2) 0.216098 10−10 (0, 2) 0.000000 10+00 (1, ¯ 1) 0.215635 10−02 (1, 1) 0.000000 10+00 (2, 0) 0.000000 10+00 Table 3 Diffraction peaks probabilities for a initial rovibrational state (v=0, J=2, mJ=0) and E⊥=400eV. 0 0.5 1 1.5 2 2.5 Prob. x 10 0.2 0.4 0.6 Normal energy (eV) 0 0.5 1 1.5 Prob. x 103 0.2 0.4 0.6 Normal energy (eV) 0 0.5 1 1.5 Prob. x 103 0.2 0.4 0.6 Normal energy (eV) 0 0.5 1 1.5 Prob. x 10 Jf=9 Jf=7 0.2 0.4 0.6 0.2 0.4 0.6 Normal energy (eV) 0 1 2 3 Prob. x 104 0.2 0.4 0.6 Normal energy (eV) 0 1 2 3 Prob. x 104 Ji=9 mJi=0 (0,0) Ji=9 mJi=0 (2,2) Ji=9 mJi=9 (0,0) Ji=9 mJi=9 (2,2) AB D C Jf=7 Jf=7 Jf=7Jf=7 _ _ Fig. 9 Elastic (black line) and rotationally deexcited probabilities (red lines) as a function of the normal incidence energy. The diffraction peak and the initial rotational state is shown in the legends. In all four case, the molecules are in the vibrational ground state. The inset show a zoom of the rotational excitation probabilities. 0 0.5 1 1.5 Prob. x 10 Jf=2 Jf=0 Jf=4 0.2 0.4 0.6 Normal energy (eV) 0 0.4 0.8 1.2 Prob. x 10 0.2 0.4 0.6 Ji=2 mJi=0 (0,0) Ji=2 mJi=0 (2,2) Ji=2 mJi=2 (0,0) Ji=2 mJi=2 (2,2) AB D C _ _ Fig. 8 Elastic (black line) and rotationally excited and deexcited probabilities (green and red lines) as a function of the normal incidence energy. The diffraction peak and the initial rotational state is shown in the legends. In all four case, the molecules are in the vibrational ground state. Finally, it is worthy to point out that none of the possible vibrational inelastic diffraction (VID’s) peaks are appreciably populated. Taken into account that the vibrational excited energy for H2(v=0) →H2(v=1) is 515.8 meV, the absence of VID’s further confirms that the energy transfer between the normal and the molecule internal motions is solely responsible for exciting the normal vibrational and rotational modes. Thus, for H2/KCl(001) under fast grazing incidence conditions the internal and the normal motions are coupled, whereas the internal and parallel motions seems to be decoupled. Notes and references 1 P. Rousseau, H. Klemliche, A. G. Borisov and P. Roncin, Phys. Rev. Lett., 2007, 98, 016104. 2 A. Schuller, S. Wethekam and H. Winter, Phys. Rev. Lett., 2007, 98, 016103. 3 H. Hoinkes, Rev. Mod. Phys., 1980, 52, 933. 4 D. Frankl, Prog. Surf. Sci., 1983, 13, 285. 5 J. Barker and D. Auerbach, Surf. Sci. Rep., 1984, 4, 1. 6 D. Farías and K. Rieder, Rep. Prog. Phys., 1998, 61, 1575. 7 D. Farías, C. Díaz, P. Nieto, A. Salin and F. Martín, Chem. Phys. Lett., 2004, 390, 250. 8 E. A. Andreev, Russ. J. Chem. Phys., 2002, 76, 5164. 9 D. Danailov, J. H. Rechtien and K. J. Snowdon, Surf. Sci., 1991, 259, 359. 10 C. Díaz, P. Rivière and F. Martín, Phys. Rev. Lett., 2009, 103, 013201. 11 A. Zugarramurdi and A. Borisov, Phys. Rev. A, 2012, 86, 062903. 12 A. Zugarramurdi and A. Borisov, Nucl. Instrum. Meth. Phys. Res. B, 2013, 317, 83. 13 A. S. Muzas, F. Gatti, F. Martín and C. Díaz, Nucl. Instr. Meth. B, 2016, 382, 49. 14 M. Busch, J. Seifert, E. Meyer and H. Winter, Phys. Rev. B, 2012, 86, 241402. 15 M. Debiossac, P. Pang and P. Roncin, Phys. Chem. Chem. Phys., 2021, 23, 7615. 16 J. Lienemann, A. Schuller, D. Blauth, J. Seifert, S. Wethekam, Journal Name, [year], [vol.], 1–11 | 9