Theoretical Simulations of Reactive and Nonreactive Scattering of Light Diatomic Molecules from Metal Surfaces : Past , Present , and Future

In everyday life we are surrounded by surfaces and, therefore, by phenomena involvingmolecule-surface interactions. Furthermore, the processes of heterogeneous catalysis, which are governed by molecule-surface interactions, are of huge practical importance, because the production ofmost synthetic compounds involves catalytic processes, which explains the tremendous effort that surface science scientists have invested to understand the basic principles underlying elementary interactions between light molecules and surfaces. This effort was recognized in 2007 with the Nobel prize in chemistry awarded to Gerhard Ertl. Here we revise some of the most relevant studies performed so far in this field. We also point out the major challenges that the surface science community may face in this field in the years to come.


Introduction
In everyday life we are surrounded by phenomena involving molecule-surface interactions.For example, the corrosion of a coin is due to the interaction between the oxygen molecules in the air and the metal surface atoms, which cause a structural damage leaving a layer of oxidized material (rust) on the coin.Another example is the green appearance of the domes of some buildings, which is due to the oxidation of copper, material from which domes are made.Atomic and molecular interactions on surfaces also play a key role in many industrial processes, such as corrosion, friction, lubrication, oxidation, hydrogen storage, and heterogeneous catalysis.Heterogeneous catalysis is of tremendous practical importance.The production of most synthetic compounds involves catalytic processes because most of the chemical reactions relevant to chemical industries are too slow in the absence of a catalyst.Therefore, understanding the basic principles that govern the geometry and electronic structure of metal surfaces and the elemental processes occurring on them, such as molecular reactivity and molecular scattering, has been and still is one major scope in surface science.At this point, it should be noticed that the importance of this research field was recognized in 2007 with the Nobel prize in chemistry awarded to Professor Ertl for the detailed description of the sequence of elementary molecule-surface reactions by which vast quantities of ammonia are produced [1].Ammonia production is basic for the fertilizers industry.
Metal surfaces serve as catalysts for many chemical reactions.Some of these reactions, more efficient on a metal surface than in the gas phase, are, for example, hydrogenation of O, C, N, and S to obtain H 2 O, CH 4 , NH 3 , and H 2 S; oxidation of ammonia to nitric acid, which is a basic reaction in the production of fertilizers; methanol synthesis from CO and H 2 ; oxidation of ethylene to ethylene oxide, a basic reaction in the production of antifreezes; and dehydrogenation of butane to butadiene, a reaction of primary importance in the production of synthetic rubber.
A detailed knowledge, at atomic scale, of the dynamic processes that govern the molecular reactions on surfaces is essential to design and develop new and improved catalysts.In this regard, detailed theoretical studies of these kinds of processes are of most fundamental interest.Theoretical simulations are used not only to provide accurate insights into experimental results, but also to predict new trends.They can be used to decrease the number of trial and error experiments conventionally used in the design of new devices.For example, five chemical compounds could be combined in thousands of different ways to build a catalyst; the use of theoretical simulations can reduce the number of experiments needed to optimize the structure of this catalyst.
A wealth of chemical and physical processes, involving light molecules, are possible on surfaces.These processes can be divided into three main groups.
(1) Molecular adsorption (see Figure 1): it is a physical process in which a molecule coming from the vacuum (gas phase) hits a surface somehow and sticks on it.This process may induce a relatively high concentration of molecules (or atoms) at the place of contact, and, therefore, to the formation of a molecular or atomic film on the surface.Four primary physical mechanisms are associated with the molecular adsorption mechanism: (a) molecular physisorption: a molecule, coming from the vacuum, gets adsorbed on the surface and weakly binds to it, due to van der Waals forces; (b) molecular chemisorption: a molecule, coming from the vacuum, gets adsorbed on the surface and strongly binds to it, due to the formation of new chemical bonds; (c) dissociative chemisorption: a molecule, coming from the vacuum, breaks its bond and its atoms get adsorbed on the surface thanks to the formation of new chemical bonds; (d) abstraction mechanism: a molecule, coming from the gas phase, breaks its bond; one of its atoms gets absorbed on the surface and the other one escapes back to the vacuum.
(2) Molecular desorption (see Figure 2): it is the opposite process to the molecular adsorption one.In this case a molecule previously absorbed on the surface is released.Physical mechanisms associated with molecular desorption are (a) Langmuir-Hinshelwood: two atoms or molecules adsorbed on a surface, in thermal equilibrium with it, meet each other and react; as a result a new molecule is formed and leaves the surface; (b) Eley-Rideal: an atom coming from the vacuum reacts with an atom adsorbed on the surface (in thermal equilibrium with it) forming a new molecule, which desorbs from the surface; (c) Hot-atom mechanism: a reaction takes place between an adsorbed atom, in thermal equilibrium with the surface, and an atom that has recently arrived from the vacuum, which is not in thermal equilibrium yet.

Experimental Techniques in Brief
Molecule-surface interactions experiments are, in general, carried out in ultrahigh vacuum systems (see [7] and references therein).These systems are formed by a number of connected chambers, which can be individually pumped.They contain a number of skimmers and orifices that are used to create a supersonic molecular beam (SMB) from an initial pulsed nozzle source, at high stagnation pressure.These SMBs are well defined by its stream velocity and the width of its velocity distribution.The average translational energy of the molecular beam can be varied by heating or cooling the nozzle.
SMB experiments can be used to measure the dissociative adsorption probability of a molecule-surface system using the so-called King and Wells' method [8].In this method, the partial pressure of the molecular gas is monitored as a function of time using a quadrupole mass spectrometer (QMS).At the beginning of the measurement the supersonic beam of molecules is produced in a secondary chamber equipped with a shutter, which can intercept the beam.While the shutter intercepts the molecular beam, the QMS, located in the main chamber, measures the background signal  0 .Once the shutter is turned down, the beam enters the main chamber increasing the signal measured by the QMS to  1 .In the main chamber a second shutter intercepts the beam preventing it from striking the surface target.When this second shutter is removed from the path of the beam, the molecules can hit and stick on the surface.Whenever adsorption takes place the QMS signal decreases as ().Thus, the relative decrease of the QMS signal gives us the absolute dissociation probability The dissociative adsorption probability can be also measured using temperature-programmed desorption (TPD) Advances in Chemistry techniques (see, e.g., [9]).In this case, the molecules, previously adsorbed on a surface, are desorbed by heating the surface and simultaneously are detected by a mass spectrometer.In these kinds of experiments, the adsorption probabilities are obtained from coverage versus exposure measurements.Thus, the adsorption probability is determined directly as the ratio of the number of molecules that adsorb to the number of molecules that strike the surface.The number of the adsorbed molecules is proportional to the area under the TPD trace.The reliability of the TPD technique lies in the validity of the detailed balance principle [10], which assumes that the dissociation probability can be measured by looking at the reverse process, the molecular desorption.
Associative desorption measurements are also performed using time-of-flight (TOF) techniques [11], which allow obtaining molecular-state-specific probabilities.In general, TOF experiments are performed on a two-chamber apparatus.One of the chambers contains the crystal into which the molecules are adsorbed after permeation through the bulk.The molecules desorbing from the surface are then probed by laser ionization detection in the second chamber.And the TOF distributions are obtained by recording the time of flight of the photoions to a multichannel plate detector.To determine quantum state distributions, the relative signal intensity for each quantum state desorbed from the surface is compared with the corresponding signal obtained with an effusive Knudsen source of molecules.
TOF experiments are also used to measure rotational and vibrational inelastic scattering probabilities, in combination with stimulated Raman pumping (SRP) [12] and resonance enhanced multiphoton ionization (REMPI) techniques.SRP can be used to excite vibrationally the initial molecular beam, thus selecting the initial vibrational state.For this aim, the initial molecular beam is crossed with two focused laser beams, which are chosen in such a way that the frequency difference between them matches the vibration of the molecule, allowing an efficient excitation of the molecules from the ground state to an excited state.The rotational and vibrational populations as well as the quadrupole alignment of the scattered molecules can be determined using REMPI [13,14].In applying this technique, the probe laser beam is focused on the molecular beam few mm in front of the surface where it ionizes the molecules, and the ions are collected and detected with a microchannel electron multiplier plate.
SMB techniques can also be used to measure diffractive scattering probabilities [15].In these kinds of experiments the monochromatic character of the molecular beams plays a crucial role; spread-velocity beams do not allow observing diffraction peaks.To measure diffraction from SMB experiments two types of apparatus are commonly used: (i) in the fixed-angle setups, the angle between the incident and the reflected beam is fixed; that is, Θ  + Θ  = const, and the angular distributions of the diffracted particles are measured by rotating continuously the crystal and thus the incident angle Θ  varies during data acquisition; (ii) in the rotary setups, the detector is able to rotate around the crystal from a given incidence angle.Therefore, the diffraction beam can be measured in a large region of the reciprocal space.The advantage of this latter setup is the possibility of determining absolute diffraction probabilities.

Molecule-Surface Dynamics Simulations
Since the early 90s the interactions of light molecules with surfaces have received more and more attention from theoretical surface scientists, essentially due to the development of multidimensional quantum dynamics methods [2,3,16,17] and the development of interpolation methods able to build flexible potential energy surfaces (see [18] and references therein).
Molecule-surface interactions have been usually described within the static surface Born-Oppenheimer approximation (SS-BOA).The validity of the SS approximation is strongly supported by the mass mismatch between the atoms of the molecule and the metal atoms of the surface whereas the BOA is supported by the velocity mismatch between nuclei and electrons, the latter ones being faster by few orders of magnitude than the former ones.Within the SS-BOA we describe the motion of the nuclei on a continuum potential energy surface (PES).Thus, the first step in any adiabatic simulation is to determine the electronic structure of the system, that is, the electronic landscape on which the nuclei move.
Although, for a number of molecule-surface systems, these analytical PESs have been able to describe successfully their electronic structure, the lack of flexibility of these PESs has impelled the development of methods based on interpolation of density functional theory (DFT) data.Generally speaking, the idea behind all these interpolation methods is quite similar.The first step is always to build a DFT data set; to this aim a number of high-and low-symmetry configurations (see Figure 4) are selected.For each configuration, defined by the position of the molecule over the surface ( and ) and its orientation ( and ), a set of  (atom-atom distance) and  (molecule-surface distance) values are computed-see Figure 5 for coordinates definition.In a second step the interpolation is performed over this data set.
First interpolation methods were applied to the development of low dimensional PESs [29,30].But, nowadays, to construct a PES based on the interpolation of a DFT data set, including all the molecular degrees of freedom (DOFs), has become a routine work thanks to the development of methods such as the corrugation reducing procedure (CRP) [31], the modified Shepard (MS) method [32,33], and the neural networks (NN) method [34].
The key idea behind the CRP method is that most of the corrugation of a molecule-surface PES is due to the atomsurface interactions and, therefore, the subtraction of this atomic contribution from the PES leads to a much smoother function, which can be interpolated more accurately.Thus, the 6D PES can be written as where  6D is the smoother function (two-body term) and   are the 3D PES representing the atom-surface interactions (one-body term).In the CRP method the interpolation is performed by combining analytical functions and numerical techniques.The major advantage of the CRP method is that the precision of the PES can be systematically improved by adding extra DFT data, whereas the major disadvantage is that it cannot be extended straightforwardly to describe polyatomic molecules (more than two atoms) interacting with surfaces.To date, the CRP method has been successfully used to build PESs for a wide variety of molecule-surface systems, for example H  [46], H 2 /Pd(100) [47], H 2 /Cu(111) [48], and H 2 /C(1 × 1)-Ti/Al(100) [49].In Figure 6 we show, as an example, a typical 2D(, ) cut of the 6D PES obtained for the system H 2 / Cu (111).
In the case of the MS method, initially developed to study gas phase reactions [50,51], the 6D PES is written as a weighted series of second-order Taylor expansions: with   representing the coordinates of the molecule,   a normalized weigh function, and   the second-order Taylor expansion of the potential centered on each data point.The main advantages of this method with respect to the CRP one are that (1) it can be extended straightforwardly to describe the interaction of polyatomic molecules with surfaces and (2) it requires a smaller number of ab initio points, because it focuses the interpolation on the dynamically relevant regions of the PES, which are found by means of classical dynamics calculations.The MS method has been already used Advances in Chemistry to describe the energy landscape of a number of systems, such as H 2 /Pt(111) [52], N 2 /stepped-Ru(0001) [53], N 2 /Ru(0001) [54], H 2 /CO/Ru(0001) [55], H 2 /Pd(111) [56], and H/H-Si( 111) [57].
The NN method takes into account the symmetries underlying the system using nonfitting functions, which do not require any assumption about the functional form of the underlying problem.In this method the 6D PES is written as where   represent the six coordinates of the H 2 molecule,  2 and  1 are nonlinear functions, and   are the parameters of the representation, the so-called weights.This method has been successfully applied to study, for example, the dissociative chemisorption of H 2 on K(2 × 2)/Pd(100) [34], O 2 on Al(111) [58], and O 2 on Ag(111) [59].Furthermore, this method can also be combined with the CRP method.For example, Ludwig and Vlachos [60] used a combination of these two methods to construct a PES representing the energy landscape for H 2 /Pt(111).

Dynamics.
Within the BOA, once the electronic potential (the PES) is known, we can solve the time-dependent nuclear Schrödinger equation: In this equation, Ψ(R, r; ) represents the wave function describing the system, Ĥ its Hamiltonian, and R and r represent the molecular center of mass and the molecular internal coordinates, respectively.In the case of a diatomic molecule, the Hamiltonian (in atomic units) can be written as In this equation  is the angle between the  and  coordinate axis [61] ( = 90 ∘ if Cartesian coordinates are used). and  are the total and reduced mass of the molecule, respectively; ĵ is the rotational operator and  6D the PES.
There are several methods to solve the time-dependent Schrödinger (TDS) equation.In the time-dependent wave packet (TDWP) method as implemented by Kroes et al. [17,62] the TDS is solved by numerically exact propagation of the wave packet using a time-independent basis-set.This method is divided into three steps: (1) the choice of the initial wave packet; (2) the propagation of the wave packet; (3) the asymptotic analysis.
For further details on this method see [3,17] and references therein.
A promising alternative method that can be also used to solve ( 5) is the so-called multiconfiguration time-dependent Hartree method (MCTDH) [63,64].The main idea behind this method is to expand the wave function in a basis-set where not only the coefficient of the expansion, but also the basis functions themselves are time-dependent.The use of time-dependent basis-sets reduces their size, decreasing the computational cost with respect to the TDWP method.
Although, in general, the nuclear motion has to be described using quantum dynamics, classical dynamics is a very useful tool to get simple physical interpretations of quantum results and experimental measurements.In the classical dynamics method a classical trajectory is obtained by integration of the classical equations of motion.
We can integrate either the Hamilton equations of motion and   being the coordinates and the conjugated momenta of the system, respectively, or the Newton equations of motion Disregarding the set of equations used, a classical dissociation probability is computed by averaging over internal coordinates and conjugated momenta of the molecules, which can be sampled using a standard Monte Carlo method.
In performing classical dynamics we can distinguish between pure classical and quasiclassical dynamics.In the latter, the zero point energy (ZPE) of the molecule is included in the calculation, whereas in the former the ZPE is assumed to be zero.In general, pure classical dynamics yields better results for nonactivated systems (see below for definition), and quasiclassical dynamics yields better results for activated systems (see below).For activated systems the ZPE may play a significant role due to the so-called vibrational softening, which happens whenever a molecule approaches an attractive surface.In this case, the attractive force between the atoms of the molecule and the surface becomes larger than the intramolecular force, and as a result the force constant associated with the vibrational motion is reduced, which induced Molecule isolated Molecule + surface   a decrease of the potential well curvature and, therefore, a decrease of the ZPE value (see Figure 7).This adiabatic energy transfer from vibration to translation is used by the molecule to overcome the reaction barrier and dissociate.Classical dynamics can also be used to study quantum observables, such as rotational and vibrational quantum numbers [65,66] and even diffraction probabilities [67][68][69].Simulating rotational excitation using classical dynamics is possible by evaluating the closest integer that satisfied the quantum rigid rotor formula [−1+(1+4 2 /ℏ 2 ) (1/2) ]/2,  being the classical angular momentum of the molecule.In the case of vibrational excitation, vibrational quantum numbers can be simulated by evaluating the closest integer that satisfied / − 1/2,  being the action variable.To simulate diffraction probabilities, the parallel momentum change is discretized by dividing the reciprocal space using as pattern the Wigner-Seitz cell associated with each lattice point (, ) (see Figure 8).Then, the classical diffraction probability, for a given diffraction peak (, ), is given by the number of trajectories leading to a parallel momentum change contained in the Wigner-Seitz cell of the (, ) peak, divided by the total number of trajectories.

A Little of History: From Low to High Dimensional Simulations
First quantum dynamics calculations of dissociative adsorption of a diatomic molecule on a metal surface were performed by Jackson and Metiu [70], for H 2 /Ni (100).In this study only the molecular bond () and the distance moleculesurface () (see Figure 5) were taken into account; that is, only the vibrational and the translational motions towards the surface were properly described.Two years later a similar study was performed for H 2 /Cu(100) [71].In 1992 Sheng and Zhang published the first 3-dimensional (3D) quantum study for H 2 /Ni(100) [72], including, in addition to  and , the polar rotational motion  (see Figure 5).Similar 3D studies were performed by Mowrey [73] for H 2 (D 2 )/Ni (100) and by Darling and Holloway [74] and Dai and Zhang [75] for dissociative adsorption of H 2 on several Cu surfaces.Other 3D simulations were performed by taking into account ,  (translational motion along the line joining two surface atoms neighbors-see Figure 5) [76].More complex simulations including a forth degree of freedom were carried out since 1994 [77,78].These simulations included , , , and  [24,79,80], with  representing the molecular azimuthal motion (see Figure 5) or , , , and - being the axis perpendicular to  (see Figure 5).Full dimensional quantum calculations, including the 6 DOFs of the molecule, were performed for H 2 /Pd(100) in 1995 [28] and later on for H 2 /Cu(111) [61] and H 2 /Cu(100) [81,82].See Section 5 for more examples.

Six-Dimensional Molecule-Surface Simulations
When a light diatomic molecule approaches a metal surface, at low energy (up to few eV's), the molecule can either dissociate or get reflected.Relative to the behavior of the dissociative adsorption probability as a function of the molecular incidence energy, diatomic molecule-metal surface systems are classified as activated and nonactivated systems (see Figure 9).Activated systems show a monotonous increase of the dissociative adsorption probability as a function of the incidence energy.This behavior is due to the presence of a minimum reaction barrier (MRB) in PES.At this point, it should be pointed out that the MRB is only one of the characteristics of the PES that influence the interaction between a molecule and a surface, but it is not the only one, as it is shown in the following examples.If the PES does not present a MRB (although other barriers are present in the PES), the dissociative adsorption probability exhibits a nonmonotonous behavior.The dissociation probability first decreases when the incidence energy increases, reaching a minimum, and then it increases with the incidence energy.This nonmonotonous behavior is associated with the socalled dynamic trapping [83,84].At low incidence energies the molecules are attracted by the attractive regions of the PES and get trapped on the surface, rebounding several times until they find a reactive path, that is, a path without a barrier, and dissociate.The trapping probability decreases when the incidence energy increases, decreasing the dissociation probability.On the contrary, the direct dissociation probability, which only depends on the dissociation barriers, increases monotonously with the incidence energy.The combination of both phenomena yields this characteristic nonmonotonous behavior (see Figure 10).In the following, we show few significant examples of both activated and nonactivated systems.they are the simplest molecule-surface systems and, therefore, theoretical models and hypothesis can be more easily tested on them.
In Figure 12 we show some examples of dissociative adsorption of H 2 (V = 0,  = 0) on nonactivated surfaces.Once again we can see that the subtle behaviour of the dissociative adsorption probability as a function of the incidence energy is surface-dependent and that the classical trajectory method gives pretty good results in comparison with the more sophisticate, and computational demanding time, quantum calculations.
(2) The monoenergetic dissociation probabilities (  ; ) are used to simulate experimental dissociation probabilities ((  )) performing a convolution over the distribution of the molecular beam using the expression where the flux weighted velocity distribution (V;   ) is given by and  (a constant),  (the width of the velocity distribution), and V 0 (the stream velocity) define the experimental molecular beam.
The parameters describing the molecular beam (, V 0 ) may vary from one experiment to another one, which explains to a large extent the differences on the dissociation probabilities as a function of the incidence energy found experimentally for some systems.For example, experimental results on dissociative adsorption probabilities for H 2 /Cu (111) measured by Rettner et al. [92] are significantly smaller than that obtained by Berger et al. [86] (see Figure 13).At this point, it should be pointed out that for long time there Figure 13: Dissociation probabilities as a function of the incidence energy for a pure H 2 molecular beam on Cu (111).Solid lines: quantum theoretical simulations from [111].Black circles: experimental data from [92].Red triangles: experimental data from [86].
was not a clear explanation about why, at similar translational energies, these two sets of experimental dissociation probabilities differed that much.The explanation based on the different molecular beam parameters was offered in [111] through a detailed theoretical analysis.In [111] it was shown, on one hand, that an appropriate convolution of the theoretical probabilities, using the correct experimental molecular parameters, allows simulating both experimental data sets.And, on the other hand, an appropriate convolution of the theoretical probabilities is required in order to perform adequate analysis of the experimental data.This kind of convolution is, since then, performed routinely [49,55].These theoretical simulations can also be used to analyze the associative desorption [114] probability, and other observables [104], such as the following.
(i) the vibrational ( ] ) and rotational (  ) efficacy, which are used to evaluate the effectiveness of the vibrational and rotational energy in promoting reaction. ] can be computed as where  0 (]) is the incidence energy needed to obtain a reaction probability  equal to half of the saturation value when the molecule is initially in the vibrational state ] and  ] is the vibrational energy of the molecule in the gas phase.And   can be computed by the slope of the line where   represents the rotational energy.
(ii) The associative desorption probabilities (  (], ; )) can be computed, invoking detailed balance, from the state-resolved dissociative adsorption probabilities as where  is the surface temperature and  the Boltzmann constant.
(iii) The average desorption energies can be computed from state-resolved dissociative adsorption probabilities as (iv) The quadrupole alignment parameter ( (2)  0 ()) gives us a measurement of the reactivity of the molecule as function of its orientation, that is, as a function of   . (2)  0 () can be computed using the following equation: Thus, a positive value indicates a preference for reactivity when the molecule is oriented parallel (helicopter) to the surface (  = ), whereas a negative value indicates, on the contrary, a preference for reactivity when the molecule is oriented perpendicular (cartwheel) to the surface (  = 0).
(v) The vibrational and rotational excitation probabilities of molecules scattered from surfaces.

Scattering of H
2 from Metal Surfaces.The scattering of H 2 from metal surfaces at low energy, below 1 eV, can be considered as the complementary process to dissociative adsorption-although, strictly speaking, it is the complementary process to sticking, which includes dissociative adsorption and molecular adsorption.By measuring the molecular distribution as a function of the scattering angle (Θ), and the molecular rotational and vibrational excitation upon scattering [115,116], we can obtain extra information about the corrugation of the PES and, therefore, about the reactivity of the system.From a theoretical point of view, molecular scattering is a very useful phenomenon to test our tools.The behavior of a scattered molecule depends very much on the subtle characteristics of the PES; therefore the study of the molecular scattering procces allows us to evaluate the accuracy of PES and the accuracy of the dynamics methods.For example, a detailed comparison between the experimental measurements and the theoretical 6D simulations, for vibrational and rotational survival probabilities and rotational excitation probabilities [65], had revealed the minor role played by nonadiabatic effects on the scattering of H 2 and D 2 from Cu (111) and accordingly on reactivity.On the other hand, TOF simulations for H 2 (V = 1,  = 3)/Cu (111) [117] have shown that the state-of-the-art adiabatic theoretical models Figure 14: TOF spectrum of H 2 (] = 0) scattered from Cu (111).Black dashed line: experimental data (  = 400 k) taken from [91]; red solid line: theoretical data from [117].
underestimate the vibrational excitation.This phenomenon can be observed in Figure 14, where it can be seen that the experimental gain peak (at short times), due to vibrational excitation from H 2 (V = 0) to H 2 (V = 1,  = 3), is three times higher than the simulated one, whereas the agreement between theory and experiment for the specular peak is pretty good.In [117] it has been suggested that the discrepancy may be due to the use of the SS-BOA.As we discuss below (see Section 6), one the challenges that surface science theorists will have to face in the next future will be to go beyond this approximation.
The theoretical TOF probabilities shown in Figure 14 were computed by using the following equation: where   (  ) represents the distance travelled by the molecule from the source (surface) to the surface (detector), V  is the velocity of the scattered molecule,  ], represents the Boltzmann population of the initial molecule state (], ) in the molecular beam divided by the population of the final state (]  ,   ) (at the nozzle temperature used in the experiment), and (V,  → V  ,   ) is the transition probability from the initial state (], ) to the final state (]  ,   ).The angular distributions of H 2 , upon scattering from metal surfaces, have been studied, for example, for Pd(111) [118] and Pd (110) [119].These studies revealed a quite different behavior of the angular distributions (see Figure 15).A detailed analysis of these distributions revealed that the H 2 molecules are reflected closer to Pd (110) than to the Pd (111), and therefore they are more sensitive to the surface corrugation in the case of Pd (100).This analysis also allowed establishing the signature of the dynamic trapping, a cosinelike angular distribution of the reflection probability per solid angle (see inset Figure 15).Here, it is important to remark that the same cosine-like distribution is observed in trappingdesorption experiments [120].But the origin of trappingdesorption is totally different from the origin of dynamic trapping.The former phenomenon is due to the thermal equilibrium, that is, to the energy exchange between the surface and the molecule, which requires long interaction times (millisecond) in contrast with dynamic trapping that is a faster process (picoseconds).

Diffraction of H
2 from Metal Surfaces.The diffraction of hydrogen molecules from a metal surface [6] can be considered as a unique technique to gauge the moleculesurface PES and the dynamics.First detailed comparison between experimental diffraction data and 6D adiabatic quantum calculations was performed on H 2 /Pt(111) [121].The agreement between the experimental diffraction peaks probabilities and the theoretical ones, for several incidence energies, was found to be remarkable-experimental results were extrapolated to 0 K surface temperature using the Debey-Wallet attenuation model (see [122] and references therein).This excellent agreement, together with the very good agreement obtained for dissociative adsorption [62], allowed discarding any significant role of nonadiabatic effects on this system.A good agreement with experimental diffraction peaks was also obtained for D 2 /NiAl(110) [123,124].Although, in this latter case, the agreement is not that good for rotationally inelastic diffraction (RID).In this case, it has been suggested that these discrepancies may be due to inaccuracies in the PES [124], discarding surface phonons and significant nonadiabatic effects.
The systems mentioned above are activated systems, i.e., under experimental conditions most of the molecules (if not all of them) are reflected, which simplify the uptake of experimental measurements.But diffraction of nonactivated systems represents a major challenge, because under experimental conditions most of the molecules dissociate; therefore, the number of the scattered molecules is very low, and the surface gets dirty (covered with hydrogen) very fast.Thus, the surface temperature needs to be kept high enough to avoid hydrogen coverage, increasing phonons background and making more difficult the extrapolation of the measurements to 0 K. First measurements and subsequent simulations on these kinds of systems were carried out on H 2 /Pd(111) [123].The most remarkable result obtained from this study was a pronounced out-of-plane diffraction, showing a highly corrugated PES.This pronounced out-of-plane diffraction was considered to be associated with high reactive systems.This study revealed the importance of scanning the whole space (in-plane and out-of-plane) in order to infer trustworthy properties of the PES.Diffraction was also studied for Pd (110) [125].In this case, the theoretical simulations showed a large number of diffraction peaks, with very low probability, which prevent the experimental observation.The presence of many diffraction peaks with low intensity was considered a signature of the dynamic trapping.
On the other hand, diffraction measurements of H 2 from Cu (111) have revealed a very low out-of-plane diffraction probability [69].This result, together with previous results for H 2 diffraction from Pt (111), NiAl (110), Pd (111), and Pd (110), was induced to suggest that the presence of intense out-ofplane diffraction peaks in the diffraction spectrum could be considered a signature of a notable dissociative adsorption probability.However, detailed analyses of the systems H 2 /Ru(0001), H 2 /Cu (111), and H 2 /CuRu(0001) have shown that there is not a lineal relationship between dissociative adsorption and out-of-plane diffraction [69].In Figure 16 we can observe that H 2 /Ru(0001) is the most reactive system, whereas the most intense out-of-plane diffraction is found for H 2 /CuRu(0001).

Perspectives and Future Challenges
The surface static Born-Oppenheimer approximation has been essential to advance in the field of surface dynamics, but recent experiments on reactive and nonreactive scattering have questioned its applicability [126][127][128].Therefore, some of the major challenges that theorists are facing in the next future are related to the development of theoretical models, which will allow including accurately electrons and phonons excitations [129,130].
6.1.Beyond the Born-Oppenheimer Approximation.The possible influence of nonadiabatic effects and concretely of electron-hole (e-h) pair excitations is currently under debate within the surface dynamics community.This debate comes from experiments showing unexpected results.For example, experimental results for N 2 /Ru(0001) showing very low dissociative adsorption probabilities, for incidence energies well above the minimum reaction barrier [144,145], or the isotope-dependence behavior found for vibrational deexcitation of H 2 and D 2 upon scattering from Cu(100) have been considered as a fingerprint of strong nonadiabatic effects [115,116].This conclusion was supported by low dimensional adiabatic calculations, which were unable to reproduce these experimental results [145,146].However, subsequent theoretical studies, using high (6D) dimensional adiabatic calculations, have refuted this conclusion [54,65].But stronger experimental indications of nonadiabatic effects, coming from experiments showing e-h pair excitation accompanying molecular chemisorption [126] and ejection of electrons accompanying scattering of highly vibrationally excited molecules, call for the development of theoretical methods able to account for these effects [128].
The first attempts to include electron-hole pair excitations in dynamics studies of molecular processes on metal surfaces date back to the 80s [147][148][149].Most recently, several DFTbased approximate methods have been proposed to include nonadiabatic effects in the dynamics.Luntz et al. [146] have proposed to compute friction coefficients (), describing electron-hole pair damping along the minimum energy barrier reaction path.This method allows performing quasiclassical nonadiabatic dynamics simulations on 2D PESs.But it cannot be extrapolated straightforwardly to multidimensional (6D) calculations.In view of the importance of taking into account the full dimensionality the system may have on the appropriate description of many molecule-surface phenomena, the applicability of this method is rather limited.Juaristi et al. [150] have proposed an approximate method to include nonadiabatic effects by keeping the multidimensionality of the problem.The main approximation of this method, from now on called LDFA (local density friction approximation), is to consider that the atoms of the molecule move independently.Using the LDFA, the nonadiabaticity is introduced in the classical equations of motion, for each atom of the molecule, through a dissipative force.Thus, the equation of motion for each atom can be written as where the second term in the right-hand side of this equation is the dissipative force experienced by each atom.Shenvi et al. [151] have proposed to take into account nonadiabatic effects by using an independent electron surface hopping (IESM), in which the energy transfer to e-h pair excitations is described by considering hops between electronic adiabatic states.Using the IESM model, these authors have ).Black circles: experimental results from [131]; blue triangles: adiabatic theoretical results from [132]; red squares: IESM results from [132].
been able to reproduce qualitatively the experimental vibrational deexcitation probabilities obtained for NO/Au(111) [131,132] (see Figure 17), a system for which nonadiabatic effects have been found to be quite important [128,152].However, a model to accurately describe nonadiabatic effects in reactive and nonreactive scattering of molecules from metal surface is still to be developed.

Surface Temperature:
The Effect of Phonons.Although the static surface approximation has yielded qualitatively good results for reactive and non-reactive scattering for many molecule-surface systems, including the effect of surface temperature may be crucial to accurately describe some systems and/or observables.For example, there are strong indications that vibrational excitation of H 2 upon scattering from Cu(111) can only be accurately described if phonons are taken into account [117].So far, several methods have been proposed to include surface temperature, all of them within the classical dynamics framework.The simple one is the so-called surface oscillator (SO) model [153], in which a 3D harmonic oscillator is used to describe the collective motion of the surface atoms (see Figure 18).In this model, the coupling between the molecule and the surface atoms motion is described by a rigid coordinate shift of the 6D PES.Thus, the 9D PES can be written as  (111).Black solid circles: experimental results from [93]; blue squares: static surface approximation results from [104]; green triangles: surface oscillator model results from [104]; red diamonds: AIMD results from [106].
where  ,, are the surface oscillation frequencies,  is the surface atoms mass, and  ,, are the vectors (, , ) defining the position of the molecule and surface atoms.
This method has been successfully used to study, for instance, the rotational excitation of H 2 upon scattering from Pd (111) [66] and the quadrupole alignment parameters for D 2 /Cu (111).In the latter case, it was found that though the SO model yields better results than the static surface one, it is still too simple to obtain chemical accuracy [104].Much better results have been obtained using ab initio molecular dynamics (AIMD) [106] (see Figure 19).This method allows the motion of the surface atoms in such a way that the intricate molecule-phonons coupling can be taken into account.AIMD, used for the first time to study moleculesurface reactions by Groß and Dianat [154], computes the forces on the fly which increases hugely the computational cost, limiting the number of classical trajectories that can be computed.
An improvement over the SO model is the generalized Langevin oscillator (GLO) model [155] adapted to study molecule-surface interactions by Busnengo et al. [156] to analyze the effect of surface temperature in the dynamic trapping mediated adsorption process of H 2 on several Pd surfaces.This method includes dissipation and thermal fluctuations, thanks to ghost 3D oscillator (see Figure 20), which is coupled to the surface oscillator.In this model, the ghost oscillator is subject to friction and random forces.Within this framework model, the equations of motion for the 12 coordinates of the system are given by Here (  ,   ,   ) and   are the coordinates and frequency matrix associated with the ghost oscillator, Λ  is the coupling matrix, which couples the ghost and the surface oscillators,   represents the damping matrix associated with the friction force, and   (Δ) is the random force.GLO model has been used by Goikoetxea et al. [135] to study the molecular adsorption of N 2 on Fe (110), showing a strong dependence of the maximum of the molecular adsorption on the surface temperature (see Figure 21).This method has been also used to study, for example, the scattering of N 2 from W(110) [136].
Some effort to include surface phonons effects in quantum dynamics has been recently done [157].In this case, a 7D model including the perpendicular motion of the second layer atom was proposed.In spite of the big effort already invested by surface scientists, a method to fully take into account the complex motion associated with surface phonons is still to be developed.

Looking for Accurate Exchange-Correlation Functionals.
Also relative to future challenges, it is worthy to mention some of the shortcomings inherent to state-of-the-art DFT simulations.Although a number of exchange-correlation functionals, such as PW91 [158], PBE [159], RPBE [160], or a specific reaction parameter approach applied to them [111] have shown to give good qualitative (and even quantitative) results, in comparison with experiment, for H 2 reactive and nonreactive scattering, other molecules present a major challenge.For example, up to now, none of the proposed functionals have been able to yield a sufficiently accurate PES to describe the dissociative adsorption of O 2 on Al (111).For this system, experimental results [161] show a very low dissociation probability at thermal energies, increasing monotonously with the incidence energy, whereas adiabatic molecular dynamics simulations show very high reactivities [162,163] independent of the incidence energy (see Figure 22).This failure of the standard DFT functionals to accurately describe the interaction of O 2 with Al (111), also observed for adsorption of O 2 on Si (111) [164], is due to the triplet-tosinglet spin conversion, which is not properly described by standard DFT.Contrary to most diatomic molecules, O 2 in gas phase (far from the surface) is in its triplet ground state, and when approaching the surface a transition to two oxygen atoms in their spin-singlet state occurs.Aiming to overcome this shortcoming a spin-constrained DFT approach has been proposed [58,165].In this model, the spin of the O 2 molecule is constrained to the Hilbert subspace, which prevents spin quenching and charge transfer before the molecule starts interacting with the surface; that is, the molecule is forced to travel in a spin-triplet configuration up to distances close to the surface.It is worth mentioning that an accurate description of triplet-to-singlet spin conversion becomes  (111).Black solid circles: experimental results from [161]; red solid squares: spin-constrained DFT results from [58]; blue triangles: adiabatic results.
crucial mostly for systems with a low density of states at the Fermi level.For example, standard DFT adiabatic dynamics for O 2 /Ag(100) [45] agrees with experimental data in the absence of dissociative adsorption for energies below 1 eV.Finally, we should mention that standard DFT does not take into account the effect of van der Waals (vdW) interactions.However these forces may play a prominent role for polyatomic molecules interacting with surfaces [166].Attempts to include vdW forces have been made by Grimme et al. [167][168][169] using an empirical correction scheme and by Dion et al. [170][171][172] who proposed to include the vdW forces by expansion to second order in a specific quantity contained in the long range part of the correlation functional.Tkatchenko and Scheffler [173] have proposed a scheme which uses a parameter-free method to define accurately the long range vdW forces from mean-field electronic structure calculations.To conclude, it is also worth mentioning that several efficient implementations of these nonlocal density functionals are already available in commercial codes [174,175].

Figure 4 :Figure 5 :
Figure 4: Schematic representation for a FCC(111) surface of (a) three high-symmetry surface sites; (b) some representative configurations of a diatomic molecule on the surface.

Figure 6 :
Figure 6: 2D cut through the H 2 /Cu(111) PES obtained by applying the CRP interpolation method to set of DFT-PW91 data.This 2D cut represents a H 2 molecule dissociating on a bridge site (see Figure 4).

Figure 7 :
Figure 7: Schematic representation of the vibrational softening.

Figure 8 :
Figure 8: Real lattice (a) and reciprocal lattice (b) for a FCC surface.The numbers within parentheses indicate the diffraction peaks.The dashed hexagons show the Wigner-Seitz cells corresponding to some of the lattice points.

Figure 16 :
Figure 16: (a) dissociative adsorption probability of H 2 as a function of the incidence energy.(b) in-plane and out-of-plane (see Figure8) diffraction spectra for an incidence energy = 0.075 eV and Θ  = 30 ∘ .Raw data have been convoluted using a Gaussian function of width  = 0.7 ∘ , and the peaks have been normalized to the specular peak for H 2 /CuRu(0001).(c) diffraction probabilities as a function of the diffraction order-defined by concentric hexagons built around the specular peak; peaks lying on the same hexagon belong to the same order (see Figure8).Solid black bars: CuRu(0001); dashed red bar: Ru(0001); dotted bars: Cu(111).

Figure 18 :Figure 19 :
Figure 18: Schematic representation of the surface oscillator model.

Figure 20 :
Figure 20: Schematic representation of the generalized Langevin oscillator model.

Figure 21 :
Figure 21: Maximum of the molecular adsorption of N 2 on Fe(110) as a function of the surface temperature.Data taken from [135].

1 Figure 22 :
Figure 22: Dissociation probability as a function of the incidence energy for O 2 on Al(111).Black solid circles: experimental results from[161]; red solid squares: spin-constrained DFT results from[58]; blue triangles: adiabatic results.