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

Authors: C. Díaz
Journal: Advances in Chemistry
 (2014)

---

## Abstract

In everyday life we are surrounded by surfaces and, therefore, by phenomena involving molecule-surface interactions. Furthermore, the processes of heterogeneous catalysis, which are governed by molecule-surface interactions, are of huge practical importance, because the production of most 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.


---

## Body

## 1. 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 H2O, CH4, NH3, and H2S; oxidation of ammonia to nitric acid, which is a basic reaction in the production of fertilizers; methanol synthesis from CO and H2; 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 Figure1): 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 Figure2): 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.
(3)
Molecular scattering (see Figure3): a molecule, coming from the vacuum, collides with the surface and is reflected back to it. During this physical process the molecule can transfer (gain) energy to (from) the surface. Based on this energy transfer we can distinguish the following phenomena:
(a)
elastic scattering: there is not transfer of energy between the molecule and the surface during the collision;
(b)
vibrationally inelastic scattering: the molecular vibrational energy increases or decreases during the collision, due to an energy transfer towards (from) the molecular vibrational motion from (toward) the surface;
(c)
rotationally inelastic scattering: the molecular rotational energy increases or decreases during the collision, due to an energy transfer towards (from) the molecular rotational motion from (towards) the surface;
(d)
diffraction: the parallel (to the surface) momentum of the molecule changes by discrete quantities, due to the periodicity of the surface; as a result the angular distribution of the scattered molecule presents a discrete peaks distribution.(a) dissociative adsortion mechanism; (b) molecular adsorption mechanism; and (c) abstraction mechanism.
(a)
(b)
(c)Schematic representation of the (a) Langmuir-Hinshelwood mechanism; (b) Eley-Rideal mechanism; and (c) hot-atom mechanism.
(a)
(b)
(c)Schematic representation of the (a) elastic scattering mechanism; (b) rotational excitation mechanism; (c) vibrational excitation mechanism; and (d) diffraction mechanism.
(a)
(b)
(c)
(d)In this review we will focus mainly on three of these mechanisms: (i) dissociative chemisorption [2, 3]; (ii) molecular scattering; and (iii) diffraction [4–6].


## 2. 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 P
0. Once the shutter is turned down, the beam enters the main chamber increasing the signal measured by the QMS to P
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 P
(
t
). Thus, the relative decrease of the QMS signal gives us the absolute dissociation probability
(1)
R
(
t
)
=
P
1
-
P
(
t
)
P
1
-
P
0
.The dissociative adsorption probability can be also measured using temperature-programmed desorption (TPD) 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, Θ
i
+
Θ
f
=
const, and the angular distributions of the diffracted particles are measured by rotating continuously the crystal and thus the incident angle Θ
i 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.


## 3. 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.


### 3.1. Potential Energy Surfaces

The first PES describing the electronic structure of a molecule-surface system was published in 1932 by Lennard-Jones [19]. This PES was built based on an analytical expression. Similar PESs based on London-Eyring-Polanyi-Sato (LEPS) potentials have been quite popular since then (see, e.g., [20–26]). Other analytical PESs are based on symmetry adapted functions [27, 28].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 Figure4) are selected. For each configuration, defined by the position of the molecule over the surface (X and Y) and its orientation (θ and ϕ), a set of r (atom-atom distance) and Z (molecule-surface distance) values are computed—see Figure 5 for coordinates definition. In a second step the interpolation is performed over this data set.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.
(a)
(b)Figure 5
Schematic representation of the coordinate system.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 atom-surface 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(2)
V
6
D
CRP
=
S
6
D
(
X
,
Y
,
Z
,
z
,
θ
,
φ
)
+
∑
i
=
1
2
‍
R
i
(
X
i
,
Y
i
,
Z
i
)
,
where S
6
D is the smoother function (two-body term) and R
i 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 H2/Pd(111) [31], H2/Ni(100), and H2/Ni(110) [35], H2/Pt(111) and H2/Cu(100) [36], H2/Pd(110) [37], N2/W(110) [38], H2/NiAl(110) [39], H2/Pt(211) [40], H2/Ru(0001) [41], H2/Cu(110) [42], N2/W(110) [43], H2/W(100), and H2/W(110) [44], O2/Ag(100) [45], H2/Cu/Ru(0001), and H2/Pd/Ru(0001) [46], H2/Pd(100) [47], H2/Cu(111) [48], and H2/C(1
×
1)-Ti/Al(100) [49]. In Figure 6 we show, as an example, a typical 2D(r
,
Z) cut of the 6D PES obtained for the system H2/Cu(111).Figure 6
2D cut through the H2/Cu(111) PES obtained by applying the CRP interpolation method to set of DFT-PW91 data. This 2D cut represents a H2 molecule dissociating on a bridge site (see Figure 4).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:
(3)
V
6
D
MS
=
∑
j
=
1
N
data
‍
w
j
(
X
i
)
T
j
(
X
i
)
,
with X
i representing the coordinates of the molecule, w
j a normalized weigh function, and T
j 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 ofab 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 to describe the energy landscape of a number of systems, such as H2/Pt(111) [52], N2/stepped-Ru(0001) [53], N2/Ru(0001) [54], H2/CO/Ru(0001) [55], H2/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(4)
V
6
D
NN
=
f
2
(
W
01
2
+
∑
j
‍
W
j
i
2
f
1
(
W
o
j
1
+
∑
i
=
1
6
‍
W
i
j
1
X
i
)
)
,
where X
i represent the six coordinates of the H2 molecule, f
2 and f
1 are nonlinear functions, and W
i
j are the parameters of the representation, the so-called weights. This method has been successfully applied to study, for example, the dissociative chemisorption of H2 on K(2
×
2)/Pd(100) [34], O2 on Al(111) [58], and O2 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 H2/Pt(111).


### 3.2. Dynamics

Within the BOA, once the electronic potential (the PES) is known, we can solve the time-dependent nuclear Schrödinger equation:(5)
H
^
Ψ
(
R
,
r
;
t
)
=
i
∂
Ψ
(
R
,
r
;
t
)
∂
t
.
In this equation, Ψ
(
R
,
r
;
t
) represents the wave function describing the system, H
^ 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
(6)
H
^
=
-
1
2
M
∂
∂
Z
2
-
1
2
M
si
n
2
γ
[
∂
∂
X
2
-
2
cos
⁡
γ
∂
∂
X
∂
∂
Y
+
∂
∂
Y
2
]
-
1
2
μ
∂
∂
r
2
+
j
^
2
2
μ
r
2
+
V
6
D
(
X
,
Y
,
Z
,
r
,
θ
,
ϕ
)
.
In this equation γ is the angle between the X and Y coordinate axis [61] (γ
=
90
° if Cartesian coordinates are used). M and μ are the total and reduced mass of the molecule, respectively; j
^ is the rotational operator and V
6
D 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.So that, at the end of the calculation, we obtain the monoenergetic state-resolved scattering probabilitiesS
(
ν
,
J
,
n
,
m
;
E
), with ν and j being the vibrational and the rotational quantum numbers, E the incidence energy of the molecule, and (n
,
m) the diffraction state. From these S
(
ν
,
J
,
n
,
m
;
E
) probabilities, the monoenergetic dissociative adsorption probabilities (S
(
ν
,
J
;
E
)) are computed as
(7)
S
(
ν
,
J
;
E
)
=
1
-
∑
n
,
m
‍
S
(
ν
,
J
,
n
,
m
;
E
)
.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(8)
q
˙
i
=
∂
H
∂
p
i
,
p
˙
i
=
-
∂
H
∂
q
i
,
q
i and p
i being the coordinates and the conjugated momenta of the system, respectively, or the Newton equations of motion
(9)
M
i
∂
2
∂
t
2
R
i
=
-
∂
∂
R
i
V
6
D
(
R
i
)
.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 a decrease of the potential well curvature and, therefore, a decrease of the ZPE value (see Figure7). This adiabatic energy transfer from vibration to translation is used by the molecule to overcome the reaction barrier and dissociate.Figure 7
Schematic representation of the vibrational softening.Classical dynamics can also be used to study quantum observables, such as rotational and vibrational quantum numbers [65, 66] and even diffraction probabilities [67–69]. Simulating rotational excitation using classical dynamics is possible by evaluating the closest integer that satisfied the quantum rigid rotor formula [
-
1
+
(
1
+
4
L
2
/
ħ
2
)
(
1
/
2
)
]
/
2, L 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 S
/
π
-
1
/
2, S 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 (
n
,
m
) (see Figure 8). Then, theclassical diffraction probability, for a given diffraction peak (
n
,
m
), is given by the number of trajectories leading to a parallel momentum change contained in the Wigner-Seitz cell of the (
n
,
m
) peak, divided by the total number of trajectories.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.
(a)
(b)


## 3.1. Potential Energy Surfaces

The first PES describing the electronic structure of a molecule-surface system was published in 1932 by Lennard-Jones [19]. This PES was built based on an analytical expression. Similar PESs based on London-Eyring-Polanyi-Sato (LEPS) potentials have been quite popular since then (see, e.g., [20–26]). Other analytical PESs are based on symmetry adapted functions [27, 28].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 Figure4) are selected. For each configuration, defined by the position of the molecule over the surface (X and Y) and its orientation (θ and ϕ), a set of r (atom-atom distance) and Z (molecule-surface distance) values are computed—see Figure 5 for coordinates definition. In a second step the interpolation is performed over this data set.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.
(a)
(b)Figure 5
Schematic representation of the coordinate system.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 atom-surface 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(2)
V
6
D
CRP
=
S
6
D
(
X
,
Y
,
Z
,
z
,
θ
,
φ
)
+
∑
i
=
1
2
‍
R
i
(
X
i
,
Y
i
,
Z
i
)
,
where S
6
D is the smoother function (two-body term) and R
i 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 H2/Pd(111) [31], H2/Ni(100), and H2/Ni(110) [35], H2/Pt(111) and H2/Cu(100) [36], H2/Pd(110) [37], N2/W(110) [38], H2/NiAl(110) [39], H2/Pt(211) [40], H2/Ru(0001) [41], H2/Cu(110) [42], N2/W(110) [43], H2/W(100), and H2/W(110) [44], O2/Ag(100) [45], H2/Cu/Ru(0001), and H2/Pd/Ru(0001) [46], H2/Pd(100) [47], H2/Cu(111) [48], and H2/C(1
×
1)-Ti/Al(100) [49]. In Figure 6 we show, as an example, a typical 2D(r
,
Z) cut of the 6D PES obtained for the system H2/Cu(111).Figure 6
2D cut through the H2/Cu(111) PES obtained by applying the CRP interpolation method to set of DFT-PW91 data. This 2D cut represents a H2 molecule dissociating on a bridge site (see Figure 4).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:
(3)
V
6
D
MS
=
∑
j
=
1
N
data
‍
w
j
(
X
i
)
T
j
(
X
i
)
,
with X
i representing the coordinates of the molecule, w
j a normalized weigh function, and T
j 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 ofab 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 to describe the energy landscape of a number of systems, such as H2/Pt(111) [52], N2/stepped-Ru(0001) [53], N2/Ru(0001) [54], H2/CO/Ru(0001) [55], H2/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(4)
V
6
D
NN
=
f
2
(
W
01
2
+
∑
j
‍
W
j
i
2
f
1
(
W
o
j
1
+
∑
i
=
1
6
‍
W
i
j
1
X
i
)
)
,
where X
i represent the six coordinates of the H2 molecule, f
2 and f
1 are nonlinear functions, and W
i
j are the parameters of the representation, the so-called weights. This method has been successfully applied to study, for example, the dissociative chemisorption of H2 on K(2
×
2)/Pd(100) [34], O2 on Al(111) [58], and O2 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 H2/Pt(111).


## 3.2. Dynamics

Within the BOA, once the electronic potential (the PES) is known, we can solve the time-dependent nuclear Schrödinger equation:(5)
H
^
Ψ
(
R
,
r
;
t
)
=
i
∂
Ψ
(
R
,
r
;
t
)
∂
t
.
In this equation, Ψ
(
R
,
r
;
t
) represents the wave function describing the system, H
^ 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
(6)
H
^
=
-
1
2
M
∂
∂
Z
2
-
1
2
M
si
n
2
γ
[
∂
∂
X
2
-
2
cos
⁡
γ
∂
∂
X
∂
∂
Y
+
∂
∂
Y
2
]
-
1
2
μ
∂
∂
r
2
+
j
^
2
2
μ
r
2
+
V
6
D
(
X
,
Y
,
Z
,
r
,
θ
,
ϕ
)
.
In this equation γ is the angle between the X and Y coordinate axis [61] (γ
=
90
° if Cartesian coordinates are used). M and μ are the total and reduced mass of the molecule, respectively; j
^ is the rotational operator and V
6
D 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.So that, at the end of the calculation, we obtain the monoenergetic state-resolved scattering probabilitiesS
(
ν
,
J
,
n
,
m
;
E
), with ν and j being the vibrational and the rotational quantum numbers, E the incidence energy of the molecule, and (n
,
m) the diffraction state. From these S
(
ν
,
J
,
n
,
m
;
E
) probabilities, the monoenergetic dissociative adsorption probabilities (S
(
ν
,
J
;
E
)) are computed as
(7)
S
(
ν
,
J
;
E
)
=
1
-
∑
n
,
m
‍
S
(
ν
,
J
,
n
,
m
;
E
)
.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(8)
q
˙
i
=
∂
H
∂
p
i
,
p
˙
i
=
-
∂
H
∂
q
i
,
q
i and p
i being the coordinates and the conjugated momenta of the system, respectively, or the Newton equations of motion
(9)
M
i
∂
2
∂
t
2
R
i
=
-
∂
∂
R
i
V
6
D
(
R
i
)
.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 a decrease of the potential well curvature and, therefore, a decrease of the ZPE value (see Figure7). This adiabatic energy transfer from vibration to translation is used by the molecule to overcome the reaction barrier and dissociate.Figure 7
Schematic representation of the vibrational softening.Classical dynamics can also be used to study quantum observables, such as rotational and vibrational quantum numbers [65, 66] and even diffraction probabilities [67–69]. Simulating rotational excitation using classical dynamics is possible by evaluating the closest integer that satisfied the quantum rigid rotor formula [
-
1
+
(
1
+
4
L
2
/
ħ
2
)
(
1
/
2
)
]
/
2, L 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 S
/
π
-
1
/
2, S 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 (
n
,
m
) (see Figure 8). Then, theclassical diffraction probability, for a given diffraction peak (
n
,
m
), is given by the number of trajectories leading to a parallel momentum change contained in the Wigner-Seitz cell of the (
n
,
m
) peak, divided by the total number of trajectories.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.
(a)
(b)


## 4. 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 H2/Ni(100). In this study only the molecular bond (r) and the distance molecule-surface (Z) (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 H2/Cu(100) [71]. In 1992 Sheng and Zhang published the first 3-dimensional (3D) quantum study for H2/Ni(100) [72], including, in addition to r and Z, the polar rotational motion θ (see Figure 5). Similar 3D studies were performed by Mowrey [73] for H2(D
2)/Ni(100) and by Darling and Holloway [74] and Dai and Zhang [75] for dissociative adsorption of H2 on several Cu surfaces. Other 3D simulations were performed by taking into account r, X (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 r, Z, θ, and ϕ [24, 79, 80], with ϕ representing the molecular azimuthal motion (see Figure 5) or r, Z, X, and Y-Y being the axis perpendicular to X (see Figure 5).Full dimensional quantum calculations, including the 6 DOFs of the molecule, were performed for H2/Pd(100) in 1995 [28] and later on for H2/Cu(111) [61] and H2/Cu(100) [81, 82]. See Section 5 for more examples.


## 5. 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 Figure9). 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 so-called 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.Figure 9
Dissociation probability of H2 on Ni(111) (black circles) and Ni(110) (red triangles). These systems are nonactivated and activated, respectively. Experimental data are taken from [176].Figure 10
Theoretical dissociation probability versus incidence energy for H2/Pd(111). Black line with circles: total probability; red line with triangles: direct probability; blue line with squares: trapping probability.


### 5.1. Dissociative Adsorption of H2 on Metal Surfaces

The dissociative adsorption probability of H2 on different metal surfaces has been widely studied during the last 20 years. This phenomenon is the first step in hydrogenation and dehydrogenation processes, which are of most importance for many industrial processes. From a purely theoretical point of view, the interest in these systems resides in the fact that they are the simplest molecule-surface systems and, therefore, theoretical models and hypothesis can be more easily tested on them.Since the development of efficient methods to build 6D PES, the number of theoretical simulations performed to study the dissociative adsorption of H2 on metal surfaces, pure metal, alloys, and precovered-metal surfaces has experienced a significant increase. For example, in Figure 11 we show the dissociative adsorption probability of H2, in its rovibrational ground state (
v
=
0
,
J
=
0
), for a number of activated surfaces—this is merely a sample of systems, by no means a complete list. From this figure two conclusions can be extracted: (i) although, qualitatively speaking, for all the systems the dissociative adsorption probability increases monotonously when the incidence energy increases, the minimum energy leading to reaction and the slope of the curve vary from one system to another; (ii) by comparing quantum (Figure 11(a)) and classical (Figure 11(b)) calculations we can conclude that classical (quasiclassical, in fact) dynamics simulations give qualitative good results, with a smaller computational cost.Dissociative adsorption probability of H2, in its rovibrational ground state (
v
=
0
,
J
=
0
), as a function of the incidence energy. (a) Quantum calculations for Cu(110) (black line) [105]; Pt(111) (red line) [177]; Ru(0001) (green line) [41]; PdRu(0001) (blue line) [178]; NiAl(110) (magenta line) [179]. (b) Quasiclassical calculations for Cu(110) (black line) [42]; W(110) (green line) [44]; Ni(111) (brown line) [35]; Ni(100) (orange line) [35]; PdRu(0001) (blue line) [46]; NiAl(110) (magenta line) [39].
(a)
(b)In Figure12 we show some examples of dissociative adsorption of H2
(
v
=
0
,
J
=
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.Dissociative adsorption probability of H2, in its rovibrational ground state (
v
=
0
,
J
=
0
), as a function of the incidence energy. (a) Quantum calculations for Pd(111) (black line) [62]; Pd(100) (green line) [28]. (b) Classical calculations for Pd(111) (black line) [83]; Ni(110) (red line) [35].
(a)
(b)Among the systems mentioned above, H2/Cu has been considered for long time the prototypical activated system, which explains the huge number of experimental [85–98] and theoretical [21, 24, 42, 61, 65, 75, 78, 79, 99–113] works devoted to this system. As a result of this important theoretical effort, H2/Cu(111) has been the first system in which theoretical simulations have been able to reproduce experimental observables with chemical accuracy [111]. To achieve this aim, beyond using a well chosen DFT functional, an appropriate comparison with experimental measurements was needed. Experimentally, dissociative adsorption probabilities are measured using SMB techniques. To know the characteristics of this beam, that is, its molecular rovibrational distribution and its energy width, is crucial to perform a meaningful comparison between theory and experiment. Thus, in order to compare with experiments, the theoretical monoenergetic state-resolved probabilities have to be convoluted using the SMB parameters. This convolution is performed through the following steps:(1)
The monoenergetic state-resolved dissociation probabilitiesS
(
ν
,
J
;
E
), obtained from the simulations, are used to compute monoenergetic probabilities S
(
T
n
;
E
), which only depend on nozzle temperature and on the incidence energy as
(10)
S
(
T
n
;
E
)
=
∑
ν
,
J
‍
F
B
(
ν
,
J
;
T
n
)
S
(
ν
,
J
;
E
)
,
with F
B
(
ν
,
J
;
T
n
) being the Boltzmann factor given by
(11)
F
B
(
ν
,
J
;
T
n
)
=
(
2
J
+
1
)
e
-
(
E
v
/
k
T
n
)
×
w
(
J
)
e
-
(
E
r
/
0.8
k
T
n
)
N
,
where N is the normalization factor and w
(
J
) is the factor characterizing the nuclear spin statistics—for example, in the case of H2 (D
2) w
(
J
=
0
)
=
(
1
/
4
)
(
2
/
3
) and w
(
J
=
1
)
=
(
3
/
4
)
(
1
/
3
).
(2)
The monoenergetic dissociation probabilitiesS
(
T
n
;
E
) are used to simulate experimental dissociation probabilities (S
(
T
n
)) performing a convolution over the distribution of the molecular beam using the expression
(12)
S
(
T
n
)
=
∫
v
=
0
v
=
∞
‍
f
(
v
;
T
n
)
S
(
T
n
;
E
)
d
v
∫
v
=
0
v
=
∞
‍
f
(
v
;
T
n
)
,
where the flux weighted velocity distribution f
(
v
;
T
n
) is given by
(13)
f
(
v
;
T
n
)
d
v
=
C
v
3
e
[
-
(
v
-
v
0
)
2
]
/
α
2
d
v
,
and C (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 H2/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 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].Figure 13
Dissociation probabilities as a function of the incidence energy for a pure H2 molecular beam on Cu(111). Solid lines: quantum theoretical simulations from [111]. Black circles: experimental data from [92]. Red triangles: experimental data from [86].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 (χ
J) efficacy, which are used to evaluate the effectiveness of the vibrational and rotational energy in promoting reaction. χ
ν can be computed as
(14)
χ
ν
=
E
0
(
ν
=
0
,
J
=
0
)
-
E
0
(
ν
=
1
,
J
=
0
)
E
ν
(
v
=
1
,
J
=
0
)
-
E
ν
(
v
=
0
,
J
=
0
)
,
where E
0(ν) is the incidence energy needed to obtain a reaction probability S equal to half of the saturation value when the molecule is initially in the vibrational state ν and E
ν is the vibrational energy of the molecule in the gas phase. And χ
J can be computed by the slope of the line
(15)
E
0
(
ν
=
0
,
J
)
=
C
-
χ
J
×
E
J
(
J
)
,
where E
J represents the rotational energy.
(ii)
The associative desorption probabilities (P
D
(
ν
,
J
;
E
)) can be computed, invoking detailed balance, from the state-resolved dissociative adsorption probabilities as
(16)
P
D
(
ν
,
J
;
E
)
=
E
e
[
-
E
/
k
T
]
S
(
ν
,
J
;
E
)
,
where T is the surface temperature and k the Boltzmann constant.
(iii)
The average desorption energies can be computed from state-resolved dissociative adsorption probabilities as(17)
〈
E
(
ν
,
J
)
〉
=
∫
‍
S
(
ν
,
J
;
E
)
E
2
e
-
E
/
k
T
d
E
∫
‍
S
(
ν
,
J
;
E
)
E
e
-
E
/
k
T
.
(iv)
The quadrupole alignment parameter (A
0
(
2
)
(
J
)) gives us a measurement of the reactivity of the molecule as function of its orientation, that is, as a function of m
j. A
0
(
2
)
(
J
) can be computed using the following equation:
(18)
A
0
(
2
)
(
J
)
=
∑
m
j
=
-
J
J
‍
S
(
ν
,
J
,
m
J
;
E
)
(
3
m
J
2
-
J
(
J
+
1
)
)
/
J
(
J
+
1
)
∑
m
j
=
-
J
J
‍
S
(
ν
,
J
,
m
j
;
E
)
.
Thus, a positive value indicates a preference for reactivity when the molecule is oriented parallel (helicopter) to the surface (m
J
=
J), whereas a negative value indicates, on the contrary, a preference for reactivity when the molecule is oriented perpendicular (cartwheel) to the surface (m
J
=
0).
(v)
The vibrational and rotational excitation probabilities of molecules scattered from surfaces.


### 5.2. Scattering of H2 from Metal Surfaces

The scattering of H2 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 H2 and D
2 from Cu(111) and accordingly on reactivity. On the other hand, TOF simulations for H2
(
v
=
1
,
J
=
3
)/Cu(111) [117] have shown that the state-of-the-art adiabatic theoretical models underestimate the vibrational excitation. This phenomenon can be observed in Figure 14, where it can be seen that the experimentalgain peak (at short times), due to vibrational excitation from H2
(
v
=
0
) to H2
(
v
=
1
,
J
=
3
), is three times higher than the simulated one, whereas the agreement between theory and experiment for thespecular 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.Figure 14
TOF spectrum of H2
(
ν
=
0
) scattered from Cu(111). Black dashed line: experimental data (T
s
=
400 k) taken from [91]; red solid line: theoretical data from [117].The theoretical TOF probabilities shown in Figure14 were computed by using the following equation:
(19)
f
(
t
,
T
n
)
=
c
+
N
×
[
∑
v
J
,
v
′
J
′
≠
‍
(
v
i
3
v
s
v
0
4
)
(
v
i
v
0
)
4
×
[
-
(
v
i
-
v
0
α
)
2
]
hhhhhhhhhhh
×
P
(
v
=
v
′
,
J
=
J
′
⟶
v
′
,
J
′
)
hhhhhhhhhhh
+
x
t
×
(
∑
v
J
,
v
′
J
′
≠
‍
(
v
i
3
v
s
v
0
4
)
hhhhhhhhhhhhhhhhhh
×
(
1
x
i
v
i
-
2
+
x
s
v
i
v
s
-
3
)
hhhhhhhhhhhhhhhhhh
×
exp
⁡
[
-
(
v
i
-
v
0
α
)
2
]
hhhhhhhhhhhhhhhhh
×
w
v
J
×
P
(
v
,
J
⟶
v
′
,
J
′
)
∑
v
J
,
v
′
J
′
≠
‍
(
v
i
3
v
s
v
0
4
)
)
]
,
where x
i
(
x
s
) represents the distance travelled by the molecule from the source (surface) to the surface (detector), v
s is the velocity of the scattered molecule, w
ν
,
J represents the Boltzmann population of the initial molecule state (ν
,
J) in the molecular beam divided by the population of the final state (ν
′
,
J
′) (at the nozzle temperature used in the experiment), and P
(
v
,
J
→
v
′
,
J
′
) is the transition probability from the initial state (ν
,
J) to the final state (ν
′
,
J
′).The angular distributions of H2, 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 H2 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 cosine-like 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 trapping-desorption experiments [120]. But the origin of trapping-desorption 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).Figure 15
Reflection probabilities per degree as a function of the scattering angle (Θ) for H2/Pd(110) (solid lines) and H2/Pd(111) (dashed lines). Data are taken from [118, 119], respectively. Black lines: incidence angle, Θ
i
=
15
°; red lines: Θ
i
=
45
°; blue lines: Θ
i
=
75
°. Inset: reflection probability per solid angle as a function of the scattering angle for H2/Pd(110) (Θ
i
=
45
∘).


### 5.3. Diffraction of H2 from Metal Surfaces

The diffraction of hydrogen molecules from a metal surface [6] can be considered as a unique technique to gauge the molecule-surface PES and the dynamics. First detailed comparison between experimental diffraction data and 6D adiabatic quantum calculations was performed on H2/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 getsdirty (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 H2/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 H2 from Cu(111) have revealed a very low out-of-plane diffraction probability [69]. This result, together with previous results for H2 diffraction from Pt(111), NiAl(110), Pd(111), and Pd(110), was induced to suggest that the presence of intense out-of-plane diffraction peaks in the diffraction spectrum could be considered a signature of a notable dissociative adsorption probability. However, detailed analyses of the systems H2/Ru(0001), H2/Cu(111), and H2/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 H2/Ru(0001) is the most reactive system, whereas the most intense out-of-plane diffraction is found for H2/CuRu(0001).(a) dissociative adsorption probability of H2 as a function of the incidence energy. (b) in-plane and out-of-plane (see Figure 8) diffraction spectra for an incidence energy = 0.075 eV and Θ
i
=
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 H2/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 Figure 8). Solid black bars: CuRu(0001); dashed red bar: Ru(0001); dotted bars: Cu(111).
(a)
(b)
(c)


## 5.1. Dissociative Adsorption of H2 on Metal Surfaces

The dissociative adsorption probability of H2 on different metal surfaces has been widely studied during the last 20 years. This phenomenon is the first step in hydrogenation and dehydrogenation processes, which are of most importance for many industrial processes. From a purely theoretical point of view, the interest in these systems resides in the fact that they are the simplest molecule-surface systems and, therefore, theoretical models and hypothesis can be more easily tested on them.Since the development of efficient methods to build 6D PES, the number of theoretical simulations performed to study the dissociative adsorption of H2 on metal surfaces, pure metal, alloys, and precovered-metal surfaces has experienced a significant increase. For example, in Figure 11 we show the dissociative adsorption probability of H2, in its rovibrational ground state (
v
=
0
,
J
=
0
), for a number of activated surfaces—this is merely a sample of systems, by no means a complete list. From this figure two conclusions can be extracted: (i) although, qualitatively speaking, for all the systems the dissociative adsorption probability increases monotonously when the incidence energy increases, the minimum energy leading to reaction and the slope of the curve vary from one system to another; (ii) by comparing quantum (Figure 11(a)) and classical (Figure 11(b)) calculations we can conclude that classical (quasiclassical, in fact) dynamics simulations give qualitative good results, with a smaller computational cost.Dissociative adsorption probability of H2, in its rovibrational ground state (
v
=
0
,
J
=
0
), as a function of the incidence energy. (a) Quantum calculations for Cu(110) (black line) [105]; Pt(111) (red line) [177]; Ru(0001) (green line) [41]; PdRu(0001) (blue line) [178]; NiAl(110) (magenta line) [179]. (b) Quasiclassical calculations for Cu(110) (black line) [42]; W(110) (green line) [44]; Ni(111) (brown line) [35]; Ni(100) (orange line) [35]; PdRu(0001) (blue line) [46]; NiAl(110) (magenta line) [39].
(a)
(b)In Figure12 we show some examples of dissociative adsorption of H2
(
v
=
0
,
J
=
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.Dissociative adsorption probability of H2, in its rovibrational ground state (
v
=
0
,
J
=
0
), as a function of the incidence energy. (a) Quantum calculations for Pd(111) (black line) [62]; Pd(100) (green line) [28]. (b) Classical calculations for Pd(111) (black line) [83]; Ni(110) (red line) [35].
(a)
(b)Among the systems mentioned above, H2/Cu has been considered for long time the prototypical activated system, which explains the huge number of experimental [85–98] and theoretical [21, 24, 42, 61, 65, 75, 78, 79, 99–113] works devoted to this system. As a result of this important theoretical effort, H2/Cu(111) has been the first system in which theoretical simulations have been able to reproduce experimental observables with chemical accuracy [111]. To achieve this aim, beyond using a well chosen DFT functional, an appropriate comparison with experimental measurements was needed. Experimentally, dissociative adsorption probabilities are measured using SMB techniques. To know the characteristics of this beam, that is, its molecular rovibrational distribution and its energy width, is crucial to perform a meaningful comparison between theory and experiment. Thus, in order to compare with experiments, the theoretical monoenergetic state-resolved probabilities have to be convoluted using the SMB parameters. This convolution is performed through the following steps:(1)
The monoenergetic state-resolved dissociation probabilitiesS
(
ν
,
J
;
E
), obtained from the simulations, are used to compute monoenergetic probabilities S
(
T
n
;
E
), which only depend on nozzle temperature and on the incidence energy as
(10)
S
(
T
n
;
E
)
=
∑
ν
,
J
‍
F
B
(
ν
,
J
;
T
n
)
S
(
ν
,
J
;
E
)
,
with F
B
(
ν
,
J
;
T
n
) being the Boltzmann factor given by
(11)
F
B
(
ν
,
J
;
T
n
)
=
(
2
J
+
1
)
e
-
(
E
v
/
k
T
n
)
×
w
(
J
)
e
-
(
E
r
/
0.8
k
T
n
)
N
,
where N is the normalization factor and w
(
J
) is the factor characterizing the nuclear spin statistics—for example, in the case of H2 (D
2) w
(
J
=
0
)
=
(
1
/
4
)
(
2
/
3
) and w
(
J
=
1
)
=
(
3
/
4
)
(
1
/
3
).
(2)
The monoenergetic dissociation probabilitiesS
(
T
n
;
E
) are used to simulate experimental dissociation probabilities (S
(
T
n
)) performing a convolution over the distribution of the molecular beam using the expression
(12)
S
(
T
n
)
=
∫
v
=
0
v
=
∞
‍
f
(
v
;
T
n
)
S
(
T
n
;
E
)
d
v
∫
v
=
0
v
=
∞
‍
f
(
v
;
T
n
)
,
where the flux weighted velocity distribution f
(
v
;
T
n
) is given by
(13)
f
(
v
;
T
n
)
d
v
=
C
v
3
e
[
-
(
v
-
v
0
)
2
]
/
α
2
d
v
,
and C (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 H2/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 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].Figure 13
Dissociation probabilities as a function of the incidence energy for a pure H2 molecular beam on Cu(111). Solid lines: quantum theoretical simulations from [111]. Black circles: experimental data from [92]. Red triangles: experimental data from [86].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 (χ
J) efficacy, which are used to evaluate the effectiveness of the vibrational and rotational energy in promoting reaction. χ
ν can be computed as
(14)
χ
ν
=
E
0
(
ν
=
0
,
J
=
0
)
-
E
0
(
ν
=
1
,
J
=
0
)
E
ν
(
v
=
1
,
J
=
0
)
-
E
ν
(
v
=
0
,
J
=
0
)
,
where E
0(ν) is the incidence energy needed to obtain a reaction probability S equal to half of the saturation value when the molecule is initially in the vibrational state ν and E
ν is the vibrational energy of the molecule in the gas phase. And χ
J can be computed by the slope of the line
(15)
E
0
(
ν
=
0
,
J
)
=
C
-
χ
J
×
E
J
(
J
)
,
where E
J represents the rotational energy.
(ii)
The associative desorption probabilities (P
D
(
ν
,
J
;
E
)) can be computed, invoking detailed balance, from the state-resolved dissociative adsorption probabilities as
(16)
P
D
(
ν
,
J
;
E
)
=
E
e
[
-
E
/
k
T
]
S
(
ν
,
J
;
E
)
,
where T is the surface temperature and k the Boltzmann constant.
(iii)
The average desorption energies can be computed from state-resolved dissociative adsorption probabilities as(17)
〈
E
(
ν
,
J
)
〉
=
∫
‍
S
(
ν
,
J
;
E
)
E
2
e
-
E
/
k
T
d
E
∫
‍
S
(
ν
,
J
;
E
)
E
e
-
E
/
k
T
.
(iv)
The quadrupole alignment parameter (A
0
(
2
)
(
J
)) gives us a measurement of the reactivity of the molecule as function of its orientation, that is, as a function of m
j. A
0
(
2
)
(
J
) can be computed using the following equation:
(18)
A
0
(
2
)
(
J
)
=
∑
m
j
=
-
J
J
‍
S
(
ν
,
J
,
m
J
;
E
)
(
3
m
J
2
-
J
(
J
+
1
)
)
/
J
(
J
+
1
)
∑
m
j
=
-
J
J
‍
S
(
ν
,
J
,
m
j
;
E
)
.
Thus, a positive value indicates a preference for reactivity when the molecule is oriented parallel (helicopter) to the surface (m
J
=
J), whereas a negative value indicates, on the contrary, a preference for reactivity when the molecule is oriented perpendicular (cartwheel) to the surface (m
J
=
0).
(v)
The vibrational and rotational excitation probabilities of molecules scattered from surfaces.


## 5.2. Scattering of H2 from Metal Surfaces

The scattering of H2 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 H2 and D
2 from Cu(111) and accordingly on reactivity. On the other hand, TOF simulations for H2
(
v
=
1
,
J
=
3
)/Cu(111) [117] have shown that the state-of-the-art adiabatic theoretical models underestimate the vibrational excitation. This phenomenon can be observed in Figure 14, where it can be seen that the experimentalgain peak (at short times), due to vibrational excitation from H2
(
v
=
0
) to H2
(
v
=
1
,
J
=
3
), is three times higher than the simulated one, whereas the agreement between theory and experiment for thespecular 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.Figure 14
TOF spectrum of H2
(
ν
=
0
) scattered from Cu(111). Black dashed line: experimental data (T
s
=
400 k) taken from [91]; red solid line: theoretical data from [117].The theoretical TOF probabilities shown in Figure14 were computed by using the following equation:
(19)
f
(
t
,
T
n
)
=
c
+
N
×
[
∑
v
J
,
v
′
J
′
≠
‍
(
v
i
3
v
s
v
0
4
)
(
v
i
v
0
)
4
×
[
-
(
v
i
-
v
0
α
)
2
]
hhhhhhhhhhh
×
P
(
v
=
v
′
,
J
=
J
′
⟶
v
′
,
J
′
)
hhhhhhhhhhh
+
x
t
×
(
∑
v
J
,
v
′
J
′
≠
‍
(
v
i
3
v
s
v
0
4
)
hhhhhhhhhhhhhhhhhh
×
(
1
x
i
v
i
-
2
+
x
s
v
i
v
s
-
3
)
hhhhhhhhhhhhhhhhhh
×
exp
⁡
[
-
(
v
i
-
v
0
α
)
2
]
hhhhhhhhhhhhhhhhh
×
w
v
J
×
P
(
v
,
J
⟶
v
′
,
J
′
)
∑
v
J
,
v
′
J
′
≠
‍
(
v
i
3
v
s
v
0
4
)
)
]
,
where x
i
(
x
s
) represents the distance travelled by the molecule from the source (surface) to the surface (detector), v
s is the velocity of the scattered molecule, w
ν
,
J represents the Boltzmann population of the initial molecule state (ν
,
J) in the molecular beam divided by the population of the final state (ν
′
,
J
′) (at the nozzle temperature used in the experiment), and P
(
v
,
J
→
v
′
,
J
′
) is the transition probability from the initial state (ν
,
J) to the final state (ν
′
,
J
′).The angular distributions of H2, 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 H2 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 cosine-like 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 trapping-desorption experiments [120]. But the origin of trapping-desorption 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).Figure 15
Reflection probabilities per degree as a function of the scattering angle (Θ) for H2/Pd(110) (solid lines) and H2/Pd(111) (dashed lines). Data are taken from [118, 119], respectively. Black lines: incidence angle, Θ
i
=
15
°; red lines: Θ
i
=
45
°; blue lines: Θ
i
=
75
°. Inset: reflection probability per solid angle as a function of the scattering angle for H2/Pd(110) (Θ
i
=
45
∘).


## 5.3. Diffraction of H2 from Metal Surfaces

The diffraction of hydrogen molecules from a metal surface [6] can be considered as a unique technique to gauge the molecule-surface PES and the dynamics. First detailed comparison between experimental diffraction data and 6D adiabatic quantum calculations was performed on H2/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 getsdirty (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 H2/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 H2 from Cu(111) have revealed a very low out-of-plane diffraction probability [69]. This result, together with previous results for H2 diffraction from Pt(111), NiAl(110), Pd(111), and Pd(110), was induced to suggest that the presence of intense out-of-plane diffraction peaks in the diffraction spectrum could be considered a signature of a notable dissociative adsorption probability. However, detailed analyses of the systems H2/Ru(0001), H2/Cu(111), and H2/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 H2/Ru(0001) is the most reactive system, whereas the most intense out-of-plane diffraction is found for H2/CuRu(0001).(a) dissociative adsorption probability of H2 as a function of the incidence energy. (b) in-plane and out-of-plane (see Figure 8) diffraction spectra for an incidence energy = 0.075 eV and Θ
i
=
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 H2/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 Figure 8). Solid black bars: CuRu(0001); dashed red bar: Ru(0001); dotted bars: Cu(111).
(a)
(b)
(c)


## 6. 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–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].The need for such models is becoming more and more visible as surface scientists are facing more and more complex systems, for example, heavy and polyatomic molecules interacting with surfaces, such as NO [128, 131, 132], O2 [45, 58, 59, 133], N2 [38, 43, 134–136], and CH4 [137–143].


### 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 N2/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 H2 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–149]. Most recently, several DFT-based 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
(20)
m
i
d
2
r
i
d
t
2
=
-
∇
i
V
6
D
(
r
i
,
r
j
)
-
η
(
r
i
)
d
r
i
d
t
,
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 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.Figure 17
Vibrational state population distribution of NO(ν
=
15)/Au(111) as a function of the vibrational quantum number (ν). Black circles: experimental results from [131]; blue triangles: adiabatic theoretical results from [132]; red squares: IESM results from [132].


### 6.2. 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 H2 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
(21)
V
9
D
(
R
A
,
R
B
,
R
S
)
=
V
6
D
(
R
A
-
R
S
,
R
B
-
R
S
)
+
M
2
(
w
x
2
X
S
2
+
w
Y
2
Y
S
2
+
w
Z
2
Z
S
2
)
,
where w
X
,
Y
,
Z are the surface oscillation frequencies, M is the surface atoms mass, and R
A
,
B
,
S are the vectors (X
,
Y
,
Z) defining the position of the molecule and surface atoms.Figure 18
Schematic representation of the surface oscillator model.This method has been successfully used to study, for instance, the rotational excitation of H2 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 usingab 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 molecule-surface 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.Figure 19
Quadrupole alignment parameter as a function of the incidence energy for H2 on Cu(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].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 H2 on several Pd surfaces. This method includes dissipation and thermal fluctuations, thanks toghost 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
(22)
d
2
R
A
,
B
d
t
2
=
-
1
m
A
,
B
∇
R
A
,
B
V
6
D
(
R
A
-
R
S
,
R
B
-
R
s
)
,
d
2
R
s
d
t
2
=
-
1
m
s
∇
R
S
V
6
D
(
R
A
-
R
s
,
R
B
-
R
S
)
-
ω
S
2
R
S
+
Λ
g
s
U
,
d
2
U
d
t
2
=
-
ω
g
2
U
+
Λ
g
s
R
S
+
γ
g
d
U
d
t
+
F
r
(
Δ
t
)
.
Here U
(
U
x
,
U
y
,
U
z
) and ω
g are the coordinates and frequency matrix associated with the ghost oscillator, Λ
g
s is the coupling matrix, which couples the ghost and the surface oscillators, γ
g represents the damping matrix associated with the friction force, and F
r
(
Δ
t
) is the random force.Figure 20
Schematic representation of the generalized Langevin oscillator model.GLO model has been used by Goikoetxea et al. [135] to study the molecular adsorption of N2 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 N2 from W(110) [136].Figure 21
Maximum of the molecular adsorption of N2 on Fe(110) as a function of the surface temperature. Data taken from [135].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.


### 6.3. 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 H2 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 O2 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).Figure 22
Dissociation probability as a function of the incidence energy for O2 on Al(111). Black solid circles: experimental results from [161]; red solid squares: spin-constrained DFT results from [58]; blue triangles: adiabatic results.This failure of the standard DFT functionals to accurately describe the interaction of O2 with Al(111), also observed for adsorption of O2 on Si(111) [164], is due to the triplet-to-singlet spin conversion, which is not properly described by standard DFT. Contrary to most diatomic molecules, O2 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 O2 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 crucial mostly for systems with a low density of states at the Fermi level. For example, standard DFT adiabatic dynamics for O2/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–169] using an empirical correction scheme and by Dion et al. [170–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].


## 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 N2/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 H2 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–149]. Most recently, several DFT-based 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
(20)
m
i
d
2
r
i
d
t
2
=
-
∇
i
V
6
D
(
r
i
,
r
j
)
-
η
(
r
i
)
d
r
i
d
t
,
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 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.Figure 17
Vibrational state population distribution of NO(ν
=
15)/Au(111) as a function of the vibrational quantum number (ν). Black circles: experimental results from [131]; blue triangles: adiabatic theoretical results from [132]; red squares: IESM results from [132].


## 6.2. 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 H2 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
(21)
V
9
D
(
R
A
,
R
B
,
R
S
)
=
V
6
D
(
R
A
-
R
S
,
R
B
-
R
S
)
+
M
2
(
w
x
2
X
S
2
+
w
Y
2
Y
S
2
+
w
Z
2
Z
S
2
)
,
where w
X
,
Y
,
Z are the surface oscillation frequencies, M is the surface atoms mass, and R
A
,
B
,
S are the vectors (X
,
Y
,
Z) defining the position of the molecule and surface atoms.Figure 18
Schematic representation of the surface oscillator model.This method has been successfully used to study, for instance, the rotational excitation of H2 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 usingab 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 molecule-surface 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.Figure 19
Quadrupole alignment parameter as a function of the incidence energy for H2 on Cu(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].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 H2 on several Pd surfaces. This method includes dissipation and thermal fluctuations, thanks toghost 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
(22)
d
2
R
A
,
B
d
t
2
=
-
1
m
A
,
B
∇
R
A
,
B
V
6
D
(
R
A
-
R
S
,
R
B
-
R
s
)
,
d
2
R
s
d
t
2
=
-
1
m
s
∇
R
S
V
6
D
(
R
A
-
R
s
,
R
B
-
R
S
)
-
ω
S
2
R
S
+
Λ
g
s
U
,
d
2
U
d
t
2
=
-
ω
g
2
U
+
Λ
g
s
R
S
+
γ
g
d
U
d
t
+
F
r
(
Δ
t
)
.
Here U
(
U
x
,
U
y
,
U
z
) and ω
g are the coordinates and frequency matrix associated with the ghost oscillator, Λ
g
s is the coupling matrix, which couples the ghost and the surface oscillators, γ
g represents the damping matrix associated with the friction force, and F
r
(
Δ
t
) is the random force.Figure 20
Schematic representation of the generalized Langevin oscillator model.GLO model has been used by Goikoetxea et al. [135] to study the molecular adsorption of N2 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 N2 from W(110) [136].Figure 21
Maximum of the molecular adsorption of N2 on Fe(110) as a function of the surface temperature. Data taken from [135].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.


## 6.3. 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 H2 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 O2 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).Figure 22
Dissociation probability as a function of the incidence energy for O2 on Al(111). Black solid circles: experimental results from [161]; red solid squares: spin-constrained DFT results from [58]; blue triangles: adiabatic results.This failure of the standard DFT functionals to accurately describe the interaction of O2 with Al(111), also observed for adsorption of O2 on Si(111) [164], is due to the triplet-to-singlet spin conversion, which is not properly described by standard DFT. Contrary to most diatomic molecules, O2 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 O2 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 crucial mostly for systems with a low density of states at the Fermi level. For example, standard DFT adiabatic dynamics for O2/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–169] using an empirical correction scheme and by Dion et al. [170–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].



---

*Source: 175351-2014-09-10.xml*