A Rodríguez-Fernández
email: alberto.rodriguez-fernandez@u-bordeaux.fr

L Bonnet

¶ C Crespos

¶ P Larrégaray

R Díez Muiño

Accuracy: The Scattering of H 2 on Pd(111)

Keywords: Réactions gaz-surface, mécanique semi-classique, dynamique moléculaire ab-initio, binning gaussien, correction d'adiabaticité, surfaces métalliques et oxydées Gas-surface reactions, semi-classical mechanics, ab-initio molecular dynamics, gaussian binning, adiabaticity correction, metallic and oxidized surfaces

L'objet de cette thèse est l'étude théorique de processus réactifs et non réactifs se produisant à l'interface gaz-solide. Deux méthodes de trajectoires classiques, différentes et complémentaires, ont été utilisées afin de simuler la dynamique de ces processus. La première met en jeu un ensemble conséquent de trajectoires classiques obtenues en résolvant numériquement les équations de Hamilton sur une surface d'énergie potentielle (SEP) construite au préalable. Ces trajectoires sont pondérées par des poids statistiques conformément à deux contraintes semiclassiques: la pondération gaussienne et la correction d'adiabaticité. Cette approche, dans un esprit quantique, a été appliquée à la collision entre H 2 et la surface de Pd(111). Dans un premier temps, nous nous sommes limités au cas où H 2 se trouve dans son état rovibrationnel fondamental. Nous avons par la suite considéré ses états rotationnels excités. Il nous est alors apparu nécessaire de modifier la correction d'adiabaticité sur la base d'arguments semiclassiques rigoureux. Dans les deux cas, les prédictions des probabilités de collage et de réflexion résolues en états se sont avérées être en accord remarquable avec celles obtenues par des calculs quantiques exacts, contrairement aux prédictions classiques standards. L'approche classique dans un esprit quantique pourrait ainsi s'avérer d'une grande utilité pour les études à venir. La seconde méthode utilisée dans ce travail, connu sous le vocable de Ab-Initio Molecular Dynamics (AIMD), permet de calculer les forces entre noyaux à partir de la théorie de la fonctionnelle densité et d'en déduire classiquement leurs déplacements. Contrairement à l'approche précédente, l'AIMD n'exige pas la construction généralement difficile d'une SEP (le prix à payer, toutefois, est que le coût numérique de chaque trajectoire est nettement plus élevé qu'avec la méthode précédente). L'AIMD nous a permis d'étudier le processus de dissociation de H 2 sur la surface de W(110). La fonctionnelle utilisée inclut un terme de van der Waals, iii iv qui provoque un acroissement de l'attraction à longue distance, compensé par une augmentation de la répulsion à courte distance. La combinaison des deux effets diminue de façon appréciable la probabilité de dissociation, alors en meilleur accord avec le résultat expérimental obtenu à l'aide d'une surface propre. Lorsque des atomes d'oxygène sont adsorbés au préalable sur la surface, la probabilité de dissociation chute considérablement. Cet effet est d'autant plus fort que la quantité d'oxygène adsorbé est forte. Un modèle de phase ordonnée a été utilisé afin d'expliquer l'absence de collage pour le taux de courverture Θ = 0.35 ML observé expérimentalement. Les atomes d'oxygène dévient les molécules H 2 des étroits passages conduisant au collage en l'absence des atomes d'oxygène. Ceci élimine toute chance de collage pour de forts taux de recouvrement. En revanche, pour de faible taux, on s'attend à ce qu'une dynamique similaire à celle sur la surface propre se produise au dessus des atomes de tungstène, et à une distance suffisamment grande des atomes d'oxygène.

General Introduction

Gas-surface interactions are key to understand many processes of industrial importance. This is one of the reasons for which they attract the attention of worldwide theoretical and experimental scientists from multiple fields [1]. Advances in these studies have contributed greatly to the chemical industry through the understanding of heterogeneous catalysis [2, 3]. Heterogeneous catalysis is widely used to speed up reactions of industrial interest, such as the synthesis of ammonia and nitric and sulphuric acids or the desulphurization of petroleum. Gas-surface interactions are relevant, as well, to the preparation and performance of semi-conductor materials [4], key components of the manifold electronic devices used in the contemporary world. A good understanding of these interactions is also key to important environmental studies related to the acid rain [5] and the ozone hole. The study of gas interactions with metallic surfaces is paving the way to a future hydrogen-based economy by contributing to the study of hydrogen purification, storage, detection, and fuel cells [6]. Not less important are the contributions to the description of gas dynamics on the surface of plasma-facing components in the future fusion reactors.

In the last decades, our understanding of gas-surface interactions has been growing at gigantic steps. The advances in high-vacuum techniques have allowed researchers in the field to obtain cleaner surfaces to be used in multiple experiments.

Also, the refinement of molecular beam experiments has greatly improved the measurement of the kinetics of reactions on surfaces in a systematic way [7]. From the theoretical point of view, the development of ab initio electronic structure calculations based on density functional theory (DFT) have very much improved the description of elementary processes at metal surfaces. The rapid evolution of computational systems and the generation of optimized codes to implement new and old methods have also contributed greatly to the growth of the field. Important advances have been made in the generation of potential energy surfaces (PES) [8, 9] Chapter 1

The Scattering of H 2 on Pd(111) The scattering of H 2 on Pd(111) is investigated using the quasi-classical trajectory (QCT) method. Gaussian binning (GB) and the adiabaticity correction (AC) are employed in order to account for product vibrational quantization and classical overweighting of adiabatic nonreactive collisions, respectively. Some antecedents to the study are introduced in Section 1.1. Details of the calculations are given in Section 1.2. In Section 1.3, publications detailing the implementation of both GB and the AC as well as the results for the study system are presented. Section 1.4

summarize the main results.

Introduction

Palladium surfaces are excellent catalysts for chemical reactions involving hydrogen, such as the hydrogenation of unsaturated organic compounds. Also, its hydrogen-absorbing properties elevate it to a prominent position in a future hydrogen-based economy [6]. Due to these characteristics, the study of the interaction of H 2 with Pd surfaces is of great importance in the industry. In the past three decades this system have been the focus of studies to analyze dissociative adsorption, inelastic scattering and surface diffraction among other processes [12-21].

The abundance of theoretical and experimental results on the interaction of H 2 with Pd surfaces and the small number of degrees of freedom (DOF) of the problem, makes it a very good candidate to be used as a benchmark system to test semiclassical corrections. Among them, the Gaussian binning (GB) and the Adiabaticity correction (AC) are of special interest when studying collisions at low incidence energies on a frozen Pd surface.

The (GB) procedure was introduced at the end of the last century [22] to take into account Bohr's quantization principle in the analysis of the final trajectory results.

After its introduction, the GB was justified in the scope of the classical S-matrix theory [23, 24]. This binning procedure came to replace the Standard binning (SB) [23, 25] in situations when there is only a few (less than ∼ 5) number of vibrational states available to the products. While in the SB all trajectories are equally contributing to the process; in the GB a gaussian weight is assigned to each trajectory based on how close the final classical actions are to their quantum counterparts. The closer they are, the higher the weight and the contribution to the process. Rigorously, Bohr's condition of quantization is exactly satisfied when the Gaussian tends to the delta function. In practice, however, Gaussians with a full width at half maximum (FWHM) of ∼ 10 % are used to obtain a statistically relevant number of trajectories in each final state. This compromise leads to calculations under reasonable limits while taking into account Bohr's quantization rule [23].

For the past two decades, the GB procedure has been applied successfully to many gas-phase processes [26-32]. The fact that Gaussians with a FWHM of ∼ 10 % are used, means that the amount of trajectories needed to obtain the same level of statistical convergence as with SB is multiplicated tenfold. Then, to quantize multi-CHAPTER 1. THE SCATTERING OF H 2 ON PD (111) ple DOF it would be necessary to run 10 n (n being the number of DOF to be quantized) more trajectories than if SB was used. In practice, this is only possible when n is small due to computational limitations. To circumvent this issue a single Gaussian can be used to weight the energy deposited on the final quantized DOF. This procedure is usually called 1GB [25, 33].

Around ten years after the introduction of the GB procedure, the AC was proposed [23, 34]. Its domain of use has been restricted so far to processes involving a well along the reaction path and a significant number of non-reactive trajectories along which the internal motion of the diatomic molecule evolves adiabatically. In these cases the reaction probability calculated using QCT-GB is significantly lower than its counterpart obtained using quantum mechanics (QM) [23, 34]. As stated in

Ref. [23] the quantum-mechanical counterpart of adiabatic non-reactive paths is a diffracted wave that mainly gets trapped within the well. Is expected that at some point this wave will contribute to both reactive and non-reactive paths. This fact, of course, is not contemplated in classical calculations and leads to an underestimation of the reaction probabilities. In the frame of the AC, we assume that the redistribution of the previously mentioned paths is done in the same proportion as reactive and non-adiabatic non-reactive paths. This equals, in the simplest approach, to ignore vibrationally adiabatic non-reactive trajectories in the statistics. This idea has been developed more thoroughly in Ref. [21] and a variation is analyzed in Subsection 1.3.2.

In the case of normal incidence of H 2 on a Pd(111) surface, the standard quasiclassical trajectory (QCT) method fails to describe the dissociative adsorption probability at low energies. The addition of GB and the AC to the standard QCT procedure resulted in a good agreement [21] with quantum time-dependent wave packet (Q-TDWP) [15] results. In this chapter we extend the work carried in Ref. [21] through:

• The application of GB to all product quantized motions [35] in Subsection 1.3.1.

• The study of initial rotationally excited H 2 [36] in Subsection 1.3.2.

• The application of a slightly different AC based on a more rigorous application of semiclassical scattering theory [36] in Subsection 1.3.2.

1.2. CALCULATION DETAILS

Calculation details

Traditionally the QCT methodology is comprised of three basic steps:

• Initial conditions are generated.

• Classic equations of motion are integrated until certain conditions are met (the atom(s) or molecule(s) reach one of the exit channels or the maximum time of simulation is reached).

• The observables are calculated through statistical analysis.

The problem at hand is the scattering of H 2 molecules on a Pd(111) surface. The coordinate system employed is depicted in

Initial conditions

The generation of most of the initial spatial coordinates for each trajectory is done using a Monte-Carlo sampling method. While the sampling over X and Y is straightforward, the process to sample r, θ, φ and their conjugate momenta p r , p θ and p φ require extra steps. In Ref. [37] the interested reader can find a comprehensive analysis of the sampling method used in this work. The initial distance Z i from the center of mass of H 2 to the surface is fixed at 7 Å. This ensures that initially, minimal or no interaction exists between the molecule and the surface. In the XY plane, all the initial positions lay into a cell of dimensions ∆x∆ * sin 60. In addition to the spatial coordinates of the H 2 molecule, the initial collision energy E i , initial rovibrational energy E i rovib (ν i , j i ) and the initial rotational state j i need to be provided as well. In this work, initial collision energies are in the range of 2-300 meV. If the molecules are initially in the ground vibrational state this energy is not enough to provoke vibrational excitations. Moreover, only a few rotational 1.2. CALCULATION DETAILS j i = 0 j i = 1 j i = 2 j i = 3 j i = 4 E i rovib (ν i = 0, j i ) (eV) 0.272 0.286 0.315 0.357 0.482 Table 1.1: Initial rovibrational energies for ground vibrational state and different rotational states.

states will be available. As discussed before this is the region where GB replaces usually SB. In all cases normal incidence is considered giving:

P i X = 0; P i Y = 0; P i Z = -2E i M (1.1)
with P i X , P i Y , and P i Z being the components of the initial linear momentum and M the total mass of the H 2 molecule. The initial rovibrational energies have been calculated based on the initial values of the vibrational (ν i ) and rotational (j i ) states. In this work, only trajectories starting at the ground vibrational state are considered.

In subsection 1.3.1 molecules start from the ground rovibrational state while in subsection 1.3.2 the initial rotational states j i = 1 -4 are studied. In Table 1.1 the values of E i rovib (ν i , j i ) for these cases are summarized. The method to obtain these values is the same used in subsection 1.3.1 through equations ( 1) and (10). The fact that E i rovib (ν i = 0, j i = 0) = 0 is the main characteristic of the QCT when comparing with classical trajectory calculations.

Evolution in time

The evolution in time of H 2 molecules colliding with the Pd(111) surface has been carried through the integration of the classical equations of motion within the Born-Oppenheimer Static Surface (BOSS) approach. In practice, this means that only the H 2 molecule atoms are allowed to move while the atoms of the surface are frozen in their equilibrium positions. The movement of the electrons is supposed instantaneous in comparison with the movement of the nucleus according to the Born-Oppenheimer approximation (BOA). It is then enough to propagate the movement of the nucleus in the potential determined by the electronic distribution.

The equations to integrate are the set of Hamilton equations given by:

qk = ∂H ∂p k ; ṗk = - ∂H ∂q k (1.2)
where q k are the coordinates and p k the conjugated momenta associated to them. The Hamiltonian H in spherical coordinates is written as:

H = 1 2M (P 2 X + P 2 Y + P 2 Z ) + 1 2µ (p 2 r + p 2 θ r 2 + p 2 φ r 2 sin 2 θ ) + V(X, Y, Z, r, θ, φ) (1.3)
with V(X, Y, Z, r, θ, φ) as the interaction potential between the H 2 molecule and the surface. In our calculations, this is obtained through the same potential energy surface (PES) used in Ref. [15, 38]. This will allow later in the chapter to compare our results with the quantum time-dependent wave packet (Q-TDWP) ones in these references.

The method used to integrate the set of equations 1.2 is the predictor-corrector of Bulirsch and Stoer [39]. This allows for very accurate calculations without requiring too much computational effort. The initial integration time step has been set to 0.01 fs and the conservation of energy is enforced with an error under 0.1 meV.

A classical trajectory will be propagated until one of the following conditions are met:

• If Z > Z i and P Z > 0 (reflection).

• If r > r dis1 and p r > 0 (adsorption).

-If the total energy of the two atoms is negative both atoms are considered adsorbed.

-If only one of the atoms has negative total energy then we have only one atom adsorption.

• If at least one of the atoms has a positive velocity in the vertical direction and its z component is under the minimum distance to the surface defined in the

PES. (absorption)

-If only one atom complies with the previous conditions one atom absorption is considered.

-If the absorption condition is valid for both atoms then two-atom absorption is considered.

CALCULATION DETAILS

• If the propagation time exceeds the time integration limit of 15 ps (dynamic trapping)

Statistical analysis

After a considerable amount of trajectories have been propagated we can calculate the distribution of the final states of the products using statistical analysis. Usually, a discretization of the classical variables is needed to compare with their quantum counterparts. For the past half-century, this has mainly been done using the SB procedure [23, 25]. In this traditional approach, a final action y of the products will always contribute to the state corresponding to its nearest integer ȳ. All trajectories will contribute then to the final states with the same weight equals to one. The probability of a molecule ending in a determined state will be:

P ȳ = N ȳ N (1.4)
with N being the total number of trajectories and N ȳ the number of trajectories ending with an action y ∈ [ ȳ -0.5, ȳ + 0.5]. As mentioned previously, for small collision energies in the present system this method fail to describe the sticking probability.

If weights in the form of:

w GB (y) = 1 √ π exp - y -ȳ 2 .
(1.5) are used instead, the qualitative behavior of the sticking curve is recovered as demonstrated in Ref. [21]. This is how the GB procedure is applied in practice.

Of course, the denominator in Eq. 1.4 needs to be changed to renormalize the probabilities. The parameter is generally chosen to ensure a FWHM of ∼ 10%. In Fig. 1.3 we can observe a graphical representation of these two binning procedures.

As can be seen, only a small part of all trajectories will contribute with meaningful weights to the final state ȳ when GB is used. This is the reason why is necessary to propagate at least ten times more trajectories when using this procedure to obtain the same statistical efficiency as when using SB. When more than one final action needs to be quantized, each trajectory can be weighted by a product of the weights in Eq. 1.5. This leads to a geometric progression of the number of trajectories needed to obtain statistically accurate results.

In most cases, this solution is not efficient or even possible. Another option is to use only one Gaussian to weight the energy deposited on the final quantized DOF [25, 33]. The implementation of both of these methods and the AC is explained in more detail in the next section.

Publications

Introduction. In the last five decades, synergetic advances in theory and molecular beam experiments associated with a tremendous increase of computing power have allowed a detailed understanding of the dynamics of chemical elementary processes, particularly through the study of state-to-state reaction and/or inelastic probabilities. [1][2][3][4][5][6][7][8][9][10][11] For chemical reactions involving a small number of degrees-of-freedom (DOF) (typically lower or equal to 6), quantum dynamics calculations might be now regarded as exact since they generally provide results in excellent agreement with experiments. 12,13 However, reaction dynamics studies focus on systems of ever increasing size 9,[11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29] for which at this time the quasi-classical trajectory (QCT) method is the approach predominantly used to calculate state-to-state probabilities. A crucial issue is thus to develop semiclassical tools which, coupled with the QCT method, make its predictions more realistic when major quantum effects come into play in the dynamics. 30 Obviously, these tools must be validated on small systems for which benchmark quantum dynamics calculations are available. The development and test of such corrections is the central focus of this work (only processes taking place in a single electronic state and involving negligible tunneling are considered in the following).

Gaussian binning (GB) 17,[30][31][32] and the adiabaticity correction (AC) 33 have been proposed some years ago in order to account for product vibrational quantization and classical overweighting of adiabatic nonreactive collisions, respectively. These corrections, both supported by semiclassical arguments, 30,34 were recently shown 30,34,35 to significantly improve the capability of the QCT method to describe gas-phase and gas-surface barrierless reactions in the quantum regime, where only a few internal states are available to the reactants and/or products.

To date, GB has been applied only to the vibrational motions of final products. The originality of the present work is that GB is applied for the first time to all product quantized motions of a relatively complex collisional problem, the non activated reaction of H 2 on Pd(111). 1,[36][37][38] These full GB calculations, combined with the AC, are shown to reproduce surprizingly well benchmark quantum time-dependent wave packet (Q-TDWP) results at low 2 CHAPTER 1. THE SCATTERING OF H 2 ON PD(111) collision energies for which the standard QCT method fails. 1,[36][37][38] This specific system is here chosen since for it, dissociative adsorption, inelastic scattering as well as surface diffraction have been intensively scrutinized both theoretically and experimentally over the last two decades. 1,[36][37][38][39][40] The Letter is laid out as follows: we first detail the implementation of our methodology, then discuss the comparison of its results with Q-TDWP calculations, and finally conclude.

Theoretical Methodology. In this work, we assume that H 2 impinges the surface at normal incidence in the rovibrational ground state (v = 0, j = 0) with the collision energy E i . Since the surface was frozen in the Q-TDWP calculations, 1 we make the same standard assumption throughout this work. H 2 may either stick on the palladium surface or bounce back to the vacuum. E i will be small enough for the first excited vibrational state of H 2 to be unavailable after rebound.

The quantum probability of the sticking event, 41 obtained from a Density-Functional-Theory-based ground-state potential energy surface, 42,43 is displayed in Fig. 1 (black squares).

The amplitude fluctuations, particularly strong for small values of E i , are Feshbach resonances (as shown further in this work, some sudden drops, for instance around 40 and 150 meV, are not due to Feshbach resonances, but to the openings of quantum states of the reflected H 2 molecule). Feshbach resonances are usually seen when a significant amount of classical orbits are trapped for some time within the interaction region before either sticking or going back to the vacuum. The smaller the collision energy, the larger the average trapping lifetime, and the higher the amplitude and frequency of the quantum fluctuations.

Within the classical picture, rebound may lead to a final vibrational energy of H 2 lower than its zero point energy (ZPE), and this is all the more probable when the initial collision energy is much smaller than the ZPE. In such a case, the classical sticking probability artificially falls with respect to the quantum value. This is clearly seen in Fig. 1 (blue circles) for collision energies lower than 50 meV (some numerical details on the QCT calculations are given below). Recently, however, we remedied this issue by taking into account Bohr's condition of quantization in the analysis of trajectory results. 34 Transposed to molecular systems, Bohr's idea is that a quantum transition can be classically described provided that one only retains those classical paths connecting integer actions of the molecular system before the transition, with integer actions after the transition (actions are the classical analogs of quantum numbers). 44 So far, we have considered only the vibrational action of H 2 , which is enough to qualitatively recover the non-monotonous behaviour of the sticking probability. 34 In this work, however, we show that accounting for all the pertinent actions of the problem (detailed further below) leads to a drastic improvement of the predictions, not only for the sticking probability, but also for the distribution of the final quantum states of the system after reflection.

However, action quantization is not the only effect which must be taken into account in the classical trajectory approach. In fact, the title reaction involves in part trapped trajectories 42 for which H 2 wanders about the paladium surface before returning back to the vacuum or fragmenting (then leading to sticking), and in part trajectories involving a single rebound against the inner repulsive wall before reflecting back. Among these paths, several bounce back against isotropic parts of the wall in such a way that the internal motion of H 2 evolves adiabatically, i.e. in a quasi unperturbed way, throughout the whole collision.

4

CHAPTER 1. THE SCATTERING OF H 2 ON PD(111)

Semiclassical developments show that the quantum counterpart of these adiabatic paths is the diffraction of the incident wave by reflection on the isotropic parts of the wall. 34 Now, diffraction excites the rotational motion of H 2 to the detriment of its recoil motion with respect to the surface (see Fig. 3 in Ref. 34 ). Moreover, H 2 is on average attracted by the Pd(111) surface when moving away from it. Since there is less recoil energy to break free of the attractive forces and reach the vacuum, one expects that part of the diffracted wave gets trapped within the interaction region. As in previous works, [33][34][35]45 we assume for simplicity's sake that all the diffracted wave is trapped. In such a case, the diffracted wave eventually contributes to both sticking and inelastic bouncing, just as the rest of the incident wave. The classical transposition of this scenario of diffraction-mediated-trapping is that the probability carried by adiabatic paths is redistributed to the remaining paths. More details on this expected mechanism can be found in Ref. 34 For simplicity's sake, we shall assume that the redistribution is democratic, which amounts to ignoring adiabatic paths and renormalize to unity the probability carried by the remaining trajectories. The results of the present work show that the combination of this adiabaticity correction and Gaussian binning leads to a remarkable description of the dynamics at low energies. [33][34][35]45 In addition, we have found no other quantum effect playing a significant role in the title reaction. In the remainder of this section we consider the practical implementation of the previous ideas.

Gaussian binning as a practical application of Bohr's condition of quantization: Actions. The vibrational action x, 44 rotational action J and diffractional actions (a n , a m ) that we shall be considering at the end of each non reactive trajectory are given by

x = 1 2π pdr - 1 2 , (1)
J = -1 + 1 + 4L 2 f 2 , (2) 5 1.3.
PUBLICATIONS

a n = ∆ 2π P f X (3)
and

a m = ∆ 4π P f X + √ 3P f Y . (4)
We work in h unit throughout this work. r is the H 2 bond length, and p its conjugate momentum. The cyclic integral in Eq. 1 is calculated by running each trajectory over one vibrational period once H 2 has reached the vacuum, i.e., no longer interacts with the surface.

L f is the modulus of the final classical rotational angular momentum L f of H 2 . Eq. 2 is deduced by equating L 2 f and J(J + 1). P f X and P f Y are the projections of the final linear momentum vector of H 2 on the X and Y axis of Pd(111), chosen in such a way that (X = 0, Y = 0) and (X = ∆ = 2.75 Å, Y = 0) are the positions of two nearest neighbour surface atoms. The Z axis is normal to the surface. a n and a m are the classical analogs of Miller indices n and m. More details on these diffractional actions can be found in the supporting information.

The fifth action of the problem is the projection J Z of J on the Z axis orthogonal to the surface, where J is the vector of modulus J parallel to L f . In contrast to (x, J, a n , a m ), however, this action is not involved in the quantization of energy partitioning (see further below). Since GB aims at best selecting those trajectories complying with this quantization, and J Z is of no use for this purpose, we ignore it from now on.

We call a and a the quadruplets (x, J, a n , a m ) and (x, J, a n , a m ), respectively, where y is the nearest integer of y. In the following, we assume that a contributes to no other quantum state than a.

Energies. At the end of a given non reactive trajectory, the total energy E satisfies

E = P 2 Z 2M + E(a) (5)
6

CHAPTER 1. THE SCATTERING OF H 2 ON PD(111)
with

E(a) = E rovib (x, J) + E (a n , a m ), (6)
E rovib (x, J) = p 2 2µ + V (r) + L 2 f 2µr 2 (7)
and

E (a n , a m ) = 8π 2 3M ∆ 2 a 2 n + a 2 m -a n a m . (8)
M and µ are the total and reduced masses of H 2 , respectively. P 2 Z /(2M ) is the energy of recoil along the Z-axis. E rovib (x, J) is the rovibrational energy of H 2 (V (r) being its potential energy in the vacuum). E (a n , a m ) is the energy in the translational motion parallel to the surface, obtained from its definition in Cartesian coordinates and Eqs. 3 and 4. E(a) is thus the energy in the degrees-of-freedom submitted to spatial constraints (bound space, surface periodicity), thus having a discrete spectrum in quantum mechanics. The energy of quantum state a, closest to a, is given by

E(a) = E rovib (x, J) + E (a n , a m ). (9)
E rovib (x, J) is calculated as follows: p, in Eq. 1, is replaced by

p = ± 2µ E trial - J( J + 1) 2µr 2 -V (r) 1/2 (10)
where E trial is initially taken at the bottom of the well of the effective potential V (r) + J( J + 1)/(2µr 2 ). This makes p equal to 0. From Eq. 1, the corresponding value of x is thus -1/2.

E trial is then increased until x = x (to the required accuracy), and E rovib (x, J) is identified as the last value of E trial . The dependence of E rovib (x, J) on x and J clearly appears in the previous reasoning. Analogously, E rovib (x, J) is univocally determined by x and J (although this is not explicit on the right-hand side of Eq. 7).

Standard Binning: SB. In this traditional approach, 17,30 trajectories are usually assigned unit statistical weight. This is what we shall assume for sticking trajectories. On the other hand, for non reactive trajectories, the final value of j can only be even since the initial value of j is 0 (this is due to parity conservation for homonuclear diatomic molecules).

Therefore, these trajectories are assigned the weight

w SB (a) = 2(1 + δ J0 ) (11)
if J is even, 0 if J is odd. The Kronecker delta takes into account the fact that for J = 0, J belongs to the range [0,0.5] whereas for J ≥ 2, J belongs to the range [ J -0.5, J + 0.5].

The blue circles in Fig. 1 were obtained by means of this approach.

Action-Based Gaussian Binning: GB. Within the Gaussian Binning procedure, 31 the statistical weight differs for each trajectory depending on the values of the final actions.

Consider the statistical weight

w GB (a) = G(x)G(J)G(a n )G(a m ) (12)
with

G(y) = 1 √ π exp - y -y 2 . ( 13
)
Bohr's condition of quantization would exactly be realized if would tend to an infinitely small (positive) value. In such a case, one would indeed have G(y) = δ(yy), 30 and only those trajectories leading to integer actions would contribute to the process. This is, however, only of formal interest. In practice, one must increase the value of until w GB (a) takes a significant value for at least ∼ 10 2 trajectories. The normalization of w GB (a) is identical to that of w SB (a), i.e., for small enough, integrating w GB (a) over the unit cube of the action space centered at a given quantum state q = (v, j, n, m) leads to unity. by means of four Gaussians, the alternative is to pseudoquantize E(a) by means of a single Gaussian. 17,46 Formally, trajectories are thus assigned the statistical weight

Energy-based Gaussian

w 1GB (a) = 1 √ π e - E(a)-E(a) 2 E(0) 2 w SB (a), (14)
where E(0) is the harmonic zero-point vibrational energy. Eq. 14 is a straightforward application of Eq. ( 13) of Ref. 17 . w 1GB (a) clearly puts strong emphasis on those trajectories leading to E(a) close to E(a). As mentioned above, is increased until w 1GB (a) takes a significant value for at least ∼ 10 2 trajectories. For a statistical (i.e., uniform) distribution of a, it was analytically proved that GB and 1GB lead to rigorously identical results. 46 Moreover, recent simulations confirmed the similarity of GB and 1GB predictions for several non statistical processes. 9 However, we shall see further below that this is only true for total energies sufficiently distant from product energy levels.

We have numerically checked that w 1GB (a), as given by Eq. 14, is normalized to unity over the unit cubes of the action space centered at the available quantum states (again, for small enough). However, we wish to emphasize that this is generally not so, and numerically renormalizing w 1GB (a) to unity is a prerequisite for a safe use of the 1GB procedure.

Practical Detection of Adiabatic Paths. As seen before, adiabatic paths involve a single rebound against weakly anisotropic repulsive walls before reflecting back to the vacuum with nearly no perturbation of the internal motion of H 2 . In practice, we have found that trajectories involving a single rebound spend less than 100 fs at a distance from the surface less than 3 Å. We have thus defined adiabatic paths by those trajectories complying with the previous criterion, and leading to the final vibrational action x within the narrow range [-0.02,0.02], centered at the initial value of x. We have checked that the width of this window does not affect the results as long as it is kept much smaller than 1, the action spacing between two neighboring vibrational states.

Sticking and Rebound Probabilities. When the AC correction is not taken into account, the sticking probability is given by:

P X S = N S N S + Σ X R (15) with Σ X R = N -N S k=1 w X (a k ). (16)
X stands for SB, GB or 1GB, according to the procedure employed to calculate the weights.

N is the total number of trajectories, and N S the number of trajectories leading to sticking.

a k is the value of a for the k th trajectory, and Σ X R is the total weight of the N -N S nonreactive trajectories.

Analogously, the probability of reflection in quantum state q = (v, j, n, m) reads

P X R (q) = Σ X R (q) N S + Σ X R (17) with Σ X R (q) = N -N S k=1 w X (a k ) δ a k q , (18)
δ a k q being equal to 1 if a k = q, 0 otherwise. When the AC correction is applied, the expressions of the probabilities are unchanged, but Σ X R and Σ X R (q) are calculated only from nonadiabatic paths.

A total of 200 energies were considered, spaced 2 meV apart on average, and 2 million trajectories were run at each energy. The Bulirsch-Stoer integrator was used, with an initial time step of 0.01 fs. was taken at 0.06 and 0.003 for the calculation of GB and 1GB weights, respectively. This ensures that the percentage of trajectories whose weights are larger than half their maximum possible value is close to 0.5% for both GB and 1GB-QCT calculations.

Their statistical efficiencies, i.e., the rates of convergence of their predictions in terms of number of trajectories run, are thus comparable. As shown further below, the above choice of parameters make the calculations well converged. 10 CHAPTER 1. THE SCATTERING OF H 2 ON PD(111) For convenience's sake, we will call the quantum-corrected QCT calculations by the acronyms of the corrections involved (GB or 1GB, and AC whenever used). In the upper panel of Fig. 2, the GB and 1GB sticking probabilities (green triangles and red circles, respectively) are compared with the Q-TDWP one (black squares). In contrast with the structureless SB probability (blue circles in Fig. 1), the GB and 1GB

probabilities involve structures related to the opening of inelastic channels (see further below) that resemble those present in the quantum probability, except that they are shifted downwards all the stronger as the collision energy is small. However, as shown in the lower panel of Fig. 2, the adiabaticity correction shifts the GB and 1GB probabilities upwards in such a way that the gap between the resulting GB-AC and 1GB-AC probabilities on one hand, and the quantum probability on the other hand, is dramatically reduced. This finding 11

1.3. PUBLICATIONS
clearly illustrates the fact that the statistical weight carried by adiabatic paths is overestimated by the Gaussian binning procedures and that its redistribution among the remaining trajectories improves the predictions.

The standard deviation (SD) for the 1GB-AC probability was obtained from five batches of one million trajectories at 40, 150, and 300 meV, corresponding to the openings of the rovibrational states of the reflected H 2 molecule. We arrived at 0.74%, 0.14% and 0.51%, respectively, i.e. 0.46% on average. Because the red curve in the lower panel of Fig. 2 was obtained from batches of 2 million trajectories, 0.46% is an upper bound of the SD corresponding to this curve. Moreover, we have recalculated the 1GB-AC probability at the previous energies for taken at 0.0015 and 0.006 instead of 0.003. The variation of the probability was found equal to 0.3% on average. These results illustrate the good level of convergence of our calculations. Note that a comparable accuracy is expected for GB probabilities, owing to the fact that the statistical efficiencies of the GB and 1GB procedures are similar (in this work).

In Fig. 3, the same comparison as in the lower panel of Fig. 2 is done for rovibrational state-resolved reflection probabilities. These are obtained by summing the full state-resolved reflection probabilities over Miller indices. The accuracy of 1GB-AC predictions is remarkable in the vicinity of thresholds (see panels (b) and (c)) and very satisfying elsewhere. Let us emphasize here that the probability of rovibrational excitation is a magnitude difficult to theoretically describe due to the high sensitivity of the final result to fine details of the interaction and dynamics.

A zoomed-in view of the Q-TDWP and 1GB-AC probabilities around their respective thresholds is given in Fig. 4 for states (0,2) and (0,4). Moreover, the contributions of the five first diffraction orders to the 1GB-AC probabilities are also shown. These contributions are obtained by summing the full state-resolved probabilities over the Miller indices corresponding to a given value of E (n, m) (see Eq. 8). For the zero order, the sum involves the single term (n, m) = (0, 0), and E (n, m) = 0; for the first order, the sum involves the six terms (n, m) = (1, 0), (-1, 0), (0, 1), (0, -1), (1, 1) and (-1, -1), and E (n, m) = 8π 2 /(3M ∆ 2 ).

We do not enter into the details for higher orders. Their thresholds are indicated by vertical dash-dotted lines in Fig. 4. Note that each abrupt increase of these fully resolved reflection probabilities is concomitent with the abrupt decrease of the sticking probability, as shown by a careful comparison between Fig. 4 and the black squares and red circles in Fig. 2. In contrast, GB-AC reflection probabilities are only satisfying far from thresholds (see Fig. 3). In order to explain why, we ignore the quantization of the diffractional actions and use the rigid rotor harmonic oscillator approximation so as to simplify the system. The intrinsically quantized part of the energy after reflection, then, reduces to 14 CHAPTER 1. THE SCATTERING OF H 2 ON PD(111)

30
E(a) = ω(x + 1 2 ) + J(J + 1) 2µr 2 e . (19)
ω is the vibrational frequency of H 2 and r e its equilibrium bond length. Since H 2 initially heats the surface with the collision energy E i while being in the rovibrational ground state of energy ω/2, the maximum possible value of E(a) is given by

E(a) = E i + ω 2 . ( 20
)
The points satisfying this identity form the upper bound of the available area in the (x, J) plane. In panels (a) and (b) of Fig. 5, this frontier is represented within the unit square associated with quantum state (0,2) for three values of E i . The available area is located on the left side of the frontier. At E a i , (0,2) is not available. E b i corresponds to its opening.

E c i is such that E c i -E b i is equal E b i -E a i .
Contour plot representations of the GB and 1GB weights, given by Eqs. 12 and 14 respectively (neglecting the diffractional actions), are also displayed. Integrating them over the region where they are larger than half their maximum possible value leads to 0.5 % in both cases (as for the calculations previously presented).

Hence, their statistical efficiencies are similar. In panel (c), the integrals of the GB and 1GB weights obtained assuming a uniform distribution of the actions limited to the left side of the frontier are represented in terms of E i . In fact, the sudden jump previously evoked appears to be described far better by the 1GB procedure than by the GB one, as expected from the contour levels in the upper panels.

Conclusion.

For the first time in quasi-classical trajectory (QCT) calculations, all the relevant final actions of a reactive molecular system have been pseudo-quantized, in the spirit of Bohr, by the Gaussian binning procedure. 30 The process under scrutiny is the scattering of H 2 on Pd(111) down to low collision energies (∼ 10 meV). 1,[36][37][38] This full Gaussian binning, combined with the adiabaticity correction successfully used in recent years for several gasphase reactions, [33][34][35]45 leads to an unprecedented agreement with state-of-the-art quantum Blue lines correspond to three different collision energies E a i , E b i and E c i whose values can be guessed from panel (c). E b i is the energy of state (0,2). The domains of the (x, J) space available at these energies are located on the left side of the blue lines in panels (a) and (b). Integrating the GB and 1GB weights assuming a uniform distribution in the available domain at E i leads to the green and red curves in panel (c), respectively. The sudden jump observed with 1GB shows the suitability of this procedure to realistically describe probabilities around energy thresholds.

16

CHAPTER 1. THE SCATTERING OF H 2 ON PD(111) dynamics for both sticking and state-resolved reflection probabilities. This is particularly true when energy-based Gaussian binning, known as 1GB, is used. 17,46 This agreement is lost if all the final actions are not pseudoquantized, especially at low collision energies for which only a few quantum states are involved. These results highlight more than ever the ability of the QCT method to accurately describe reaction dynamics in the quantum regime provided that final trajectory results are analyzed in a quantum spirit. Note that, so far, we have focused on electronically adiabatic processes invoving no tunneling. However, Gaussian binning and the adiabaticity correction could easily be coupled with surface hopping and tunneling corrections whenever necessary. When H 2 is reflected by the Pd(111) surface, the variation of the component of its velocity parallel to the surface is quantized, due to the surface periodicity defined by the primitive vectors (a 1 , a 2 ) (see panel (a) of Fig. S1). This is mathematically expressed through Bragg condition for diffraction:

K i + G nm = K f (S.1)
where K i and K f are the parallel components of the initial and final wave vectors k i and k f , respectively, and G nm is a reciprocal lattice vector defined by:

G nm = nb 1 + mb 2 (S.
a i • b j = 2πδ ij (S.3)
From Eq. S.3 and the fact that |a 1 | = |a 2 | = ∆, one finds after some steps of algebra: For a normal incidence, the X and Y components of K i are both equal to zero. From Eq.

b 1 = 2π ∆ ; - 2π ∆ √ 3 (S.4) b 2 = 0; 4π ∆ √ 3 (S.5) a 1 a 2 a 1 = a 2 = ∆ = 2.79 Å ∆ (a) b 2 b 1 b 1 = b 2 = 2.60 Å-1 (b)
S.1 and the fact that in h unit, P f = K f , we have:

P f X = G nm X = nb 1X + mb 2X = n 2π ∆ (S.6) P f Y = G nm Y = nb 1Y + mb 2Y = -n 2π ∆ √ 3 + m 4π ∆ √ 3 (S.7)
Inverting these expressions leads to:

n = ∆ 2π P f X (S.8) m = ∆ 4π P f X + √ 3P f Y (S.

Introduction

Constant progress of molecular beam experiments allow to measure with an ever increasing accuracy state-resolved integral and differential cross sections for gas-phase reactions, and sticking or state-resolved reflection probabilities for gas-surface reactions. 1,2 Similar progress are thus required as far as theoretical predictions of these measurements are concerned.

For systems of low dimensionality (typically lower than ∼ 6), exact quantum scattering calculations generally provide predictions in excellent agreement with experimental measurements. 3,4 Nowadays, however, more and more reactions of larger dimensionality are under scrutiny. [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24] For such processes, the quasi-classical trajectory (QCT) method is the only general tool applicable in practice. Besides, the QCT method has a strong interpretative power. [25][26][27] Nevertheless, quantum corrections must be introduced in these calculations in order to make them realistic when major quantum effects strongly influence the dynamics.

Since the early 2000s, two corrections have been introduced in the QCT method, namely, Gaussian binning [28][29][30][31][32][33] (GB) and the Adiabaticity correction (AC). 32,[34][35][36] GB introduces some pseudo-quantization in the QCT method by taking into account Bohr's condition of quantization for the final quantized degrees-of-freedom. If one strictly applies Bohr ideas, only those trajectories starting from the reagents with integer actions and reaching the products with integer actions should be taken into account. These trajectories, however, form a set of zero measure compared to the whole set of possible classical paths. To go round this difficulty, trajectories are assigned Gaussian statistical weights such that the closer the final actions to integer values, the larger the weight. For more than 2 or 3 quantized degrees-of-freedom (DOFs) of the final fragments, however, GB lacks precision when applied as such to a quantum state having an energy close to the total energy available. 37 In such 2 CHAPTER 1. THE SCATTERING OF H 2 ON PD(111)

a case, a more accurate variant of GB consists in Gaussian weighting the energy deposited on the final quantized DOFs, rather than the actions. 33 Since only one Gaussian weight is involved in this procedure, it is usually called 1GB. 38 We shall use both 1GB and GB in the following. The second correction, AC, takes into account the fact that the QCT method tends to overestimate the statistical weight of those nonreactive trajectories along which the internal motion of the reagents is weakly perturbed and evolves adiabatically.

Such trajectories seem to be quite common in complex-forming processes with no entrance barrier. This correction is supported by the semiclassical initial value representation (SCIVR) of molecular collisions. 35,36 In a recent work, 37 called Part I in the following, the 1GB-AC-QCT approach was applied to the gas-surface reaction H 2 + Pd(111), with H 2 in its rovibrational ground state, and nearly quantitative agreement was found with exact time-dependent quantum scattering calculations. The goal of the present work, called Part II further on, is to check whether this agreement is maintained when H 2 is rotationally excited. We have chosen this specific process because dissociative adsorption, inelastic scattering as well as surface diffraction have been intensively scrutinized both theoretically and experimentally over the last two decades. [39][40][41][42][43][44] Note that this process takes place on a single electronic state and involves negligible tunneling. Therefore, no other quantum corrrection than 1GB and AC was found to be necessary.

The paper is laid out as follows. Section 2 details our methodology, which is slightly modified on the basis of SCIVR with respect to Part I in order to improve its consistency with first principles. Briefly, the modification consists in assigning new Gaussian weights to part of the reflected paths implicated in the adiabaticity correction. Our QCT calculations in a quantum spirit, obtained using the methods of both Part I and Part II, are compared with time-dependent quantum scattering calculations in Section 3. Section 4 concludes.

Dynamical Method

We assume that H 2 in the internal state (v i = 0, j i ) impinges the surface at normal incidence with the collision energy E i . Eventually, the molecule either sticks on the palladium surface, or bounces back to the vacuum. Whatever the initial rotational excitation, E i is small enough to make vibrational excitation unavailable after reflection. Our primary goal is to simulate the dynamics of the previous process in the case of the excited rotational states j i = 1 -4, though we will also revisit the case j i = 0.

We briefly summarize here the approach described in detail in Part I. 37 For nonreactive trajectories, the four final actions of the system are the vibrational action x f , the rotational action J f and the diffractional actions (a n f , a m f ). These are given by:

x f = 1 2π pdr - 1 2 , (1)
J f = -1 + 1 + 4L 2 f 2 , (2)
a n f = ∆ 2π P f X (3)
and

a m f = ∆ 4π P f X + √ 3P f Y . (4)
r is the H 2 bond length and p its conjugate momentum. The cyclic integral in Eq. 1 is calculated by running each trajectory over one vibrational period once H 2 has reached the vacuum, i.e., no longer interacts with the surface. L f is the modulus of the final classical rotational angular momentum L f of H 2 . Eq. 2 is deduced by equating L 2 f and J f (J f + 1). P f X and P f Y are the projections of the final linear momentum vector of H 2 on the X and Y axis of the Pd(111) surface, chosen in such a way that (X = 0, Y = 0) and (X = ∆ = 2.75 Å, Y = 0) 4 CHAPTER 1. THE SCATTERING OF H 2 ON PD(111) are the positions of two nearest neighbour surface atoms. The Z axis is normal to the surface.

a n f and a m f are the classical analogs of Miller indices n and m. Actions are in h unit. We call a and ā the quadruplets (x f , J f , a n f , a m f ) and (x f , Jf , ān f , ām f ), respectively, where ū is the nearest integer of u. Finally, x i and J i are the analogs of x f and J f for the initial H 2 molecule. x i and J i are taken at v i = 0 and j i , respectively.

In the 1GB approach, 33,38 the statistical contribution to final quantum state q = (v f , j f , n f , m f ) of a non reactive trajectory ending with a is given by:

w 1GB (a, q) = 1 √ π e -[ E(a)-E(ā) 2 E(0) ] 2 w SB (a, q). (5)
The right-most term is the statistical weight assigned within the standard binning (SB)

procedure. Here, it reads

w SB (a, q) = 2(1 + δ Jf 0 )δ āq (6)
if Jf has the same parity as j i , 0 in the contrary case. The Kronecker δ āq is the usual SB weight corresponding to the boxing procedure. The factor 2 doubles the contribution of trajectories complying with parity conservation, thus, counterbalancing the fact that the remaining trajectories are ignored. The Kronecker delta δ Jf 0 takes into account the fact that for Jf = 0, J f belongs to the range [0, 0.5] whereas for Jf larger than 0, J f belongs to the range [ Jf -0.5, Jf + 0.5]. E(a) is the final energy in all the degrees-of-freedom minus the Z-axis orthogonal to the surface. E(ā) is the analogous energy for quantum state ā. The mathematical expressions of these energies are given in the section Energies of Part I (see Eqs.

(5-10)). 37 E(0) is the harmonic zero-point vibrational energy. w 1GB (a, q) clearly puts strong emphasis on those trajectories leading to E(a) close to E(q). is selected in such a way that w 1GB (a, q) takes a significant value for at least ∼ 10 2 trajectories. In practice, adiabatic non reactive paths are defined by those reflected trajectories involving a single rebound, spending less than 100 fs at a distance from the surface less than 3 Å;and leading to the final vibrational action x f within the narrow range [-0.02, 0.02]. In Part I, these 5

1.3. PUBLICATIONS paths were ignored in the statistics. Finally, reactive trajectories are assigned unit statistical weight (sticking is supposed to occur whenever r reaches 2.25 Å;and p is positive).

Calling N the total number of trajectories run, excluding adiabatic reflected paths, and N S the number of trajectories leading to sticking, the sticking probability P I S and the probability P I R (q) of reflection in quantum state q used in Part I read:

P I S = N S N S + Σ 1GB R (7) with Σ 1GB R = q Σ 1GB R (q) (8)
and

Σ 1GB R (q) = N -N S k=1 w 1GB (a k , q), (9)
and

P I R (q) = Σ 1GB R (q) N S + Σ 1GB R . (10)
In Eq. 9, the sum is over the non adiabatic reflected trajectories.

We now introduce the modification evoked in the introduction. The upper panel of Fig. 1 shows the rotational state distribution of reflected H 2 molecules for initial j i equal 0 and an initial collision energy E i of 65 meV. Two contributions are displayed, corresponding to adiabatic trajectories (blue area) and nonadiabatic trajectories (brown area). Since the discrimination between both sets of paths used in Part I involves the vibrational action but not the rotational one (see the practical definition of adiabatic paths a bit before Eq. 7), there can be no assurance that the paths which we called adiabatic in Part I are rotationally adiabatic. In other words, our criterion does not allow to distinguish paths being fully adiabatic (FA), i.e. both vibrationally and rotationally, and paths being partially adiabatic (PA), i.e. vibrationally but not rotationally (note that for simplicity's sake, we ignore diffractional actions in our reasoning). However, we have a visual way of detecting the existence of type-6 CHAPTER 1. THE SCATTERING OF H 2 ON PD(111)

FA paths: these imply in the rotational state distribution a single narrow peak centered at j i (here, 0), such as the one visible in the blue distribution of the upper panel of Fig. 1.

We will define it by Jf = 0 further on. Such peaks are indeed due to compact volumes of initial conditions which all lead to J f ∼ j i . For the process at hand, however, the peak is only visible for j i = 0 and E i lower than ∼90 meV. In the following, we will assume that the peak is entirely due to fully adiabatic paths. In all other cases, no peak has been found, as illustrated in the lower panel of Fig. 1. The rotational state distribution for j i = 1 and the same collision energy as previously is indeed flat around J f = 1. Consequently, the adiabaticity is only partial in this case. For the sake of rigor, we now modify our formulation of the statistical weights in order to take into account the previous observations. As a first step, it is necessary to examine how type-FA trajectories look like. We have thus randomly selected a few tens of paths among the ones contributing to the narrow peak previously considered. It is found that they all rebound against the Pd(111) surface in the neighborhood of a bridge site. Moreover, the distribution of the polar angle θ with respect to the Z-axis at the instant of the rebound is found from Fig. 2 to be strongly polarized, while the distribution of the azimuth angle φ is uniform. The time-dependences of the distance Z between H 2 and the surface on one hand, and the potential energy on the other hand, are displayed in Fig. 3. At time 0, the incoming H 2 is at ∼ 7 Å; from the Pd surface with a collision energy of 65 meV and a normal incidence. Until ∼ 100 fs, Z linearly varies in the course of time, thus implying that H 2 does not interact with the surface. However, the process is barrierless in the entrance channel, so at some point, the with respect to 200 fs, illustrating thereby the full adiabaticity of the internal motion of H 2 ;

back to the vacuum, nothing has changed, neither its vibrational energy, nor its rotational energy, nor the direction of its velocity vector, still orthogonal to the Pd surface. As for the potential energy (see Fig. 3), its time dependence is consistent with the previous curve.

Since the incoming H 2 is vibrating, the potential energy strongly oscillates. The minima, fitted by the red line, correspond to those instants where the bond length of H 2 takes its equilibrium value. They thus define the interaction potential energy between H 2 and the surface. The whole potential is thus given by the sum of this interaction potential and the vibrational potential of H 2 . In the region of the configuration space (X, Y, Z, r, θ, φ) visited by type-FA trajectories, the interaction potential is expected to be nearly independent on (X, Y, θ, φ). Otherwise, some energy would be transfered to these coordinates during the collision, perturbating thereby the rotational motion of H 2 as well as its motion parallel to the (X, Y ) plane, which is not the case. The interaction potential along type-FA paths can thus be approximated by a function of Z alone. Since the time dependence of Z is roughly linear except during the rebound, the Z-dependence of the previous function has nearly the same shape as the fit of the minima (red curve in Fig. 3) deep into the well. 9

PUBLICATIONS

We now modify the method of Part I within a reduced dimensional semiclassical treatment taking into account the key topographical features of the interaction potential discussed above. The dimensional reduction is aimed at simplifying the following mathematical developments while preserving the essence of the dynamics involved in the title process. The treatment is as follows: we freeze the X, Y and φ coordinates, being thus left with Z, r and θ. H 2 is thus a vibrator and planar rotator (supposed to rotate within the (X, Z) plane)

whose center-of-mass is constrained to lie on the Z-axis. In three dimensions, θ belongs to the range [0, π], but here, it belongs to the range [0, 2π]. The phase space state of the system is thus (Z, r, θ, P, p, J), where P is the linear momentum conjugate to Z and J is the rotational action (or angular momentum in h unit) conjugate to θ (p is defined right below Eq. 4). However, the semiclassical treatment used in the following makes use of action angle coordinates for the internal DOFs, so we will alternately specify the phase space state by (Z, q, θ, P, x, J), where q is the vibrational phase of H 2 and x is the vibrational action conjugate to q. We assume that H 2 is in the initial state (v i = 0, j i = 0). We denote the mass of H 2 by µ, its vibrational frequency by ω and its moment of inertia by I. We call S v f j f the probability amplitude that after reflection, H 2 is in the final state (v f , j f ). Following Miller, 45,46 this scattering matrix element can be expressed in terms of classical paths as:

S v f j f = 1 4π 2 dq s i dθ s i ∂q s f ∂q s i ∂θ s f ∂θ s i - ∂q s f ∂θ s i ∂θ s f ∂q s i 1/2 exp i (x f -v f )q s f + (J f -j f )θ s f + Ω f i (11)
with

Ω f i = - t 0 dτ Z Ṗ + q ẋ + θ J . ( 12
)
q s and θ s are the shifted angles defined by

q s = q - µZω P (13) 10
CHAPTER 1. THE SCATTERING OF H 2 ON PD(111)

and

θ s = θ - µZJ P I . ( 14
)
Contrary to q and θ, q s and θ s have the property to be constant in the asymptotic channel, i.e., they only evolve within the interaction region. Like for quantum numbers, the labels i and f assigned in Eq. ( 11) to q s and θ s refer to their initial and final values, respectively. The trajectories used to calculate S v f j f are started at Z i and stopped at Z f , distances beyond which H 2 is uncoupled with the surface. Time τ is equal to 0 at Z i and t at Z f . The initial conditions are the negative value of P corresponding to the collision energy E i , x i = v i = 0,

J i = j i = 0, q s i and θ s i . q s f , θ s f , x f , J f
and Ω f i are functions of q s i and θ s i . We note that Eq. ( 14) and the fact that J i = 0 make θ s i and θ i equal. For simplicity's sake, we ignore the phase index necessary for an accurate description of interference effects. 47 Historically, Eq. ( 11) is the starting point of the SCIVR, 45,46,48,49 remarkably accurate for direct inelastic collisions, 45,47,50 direct collinear reactions 51,52 and molecular spectroscopy. 53,54 The interested reader will find detailed developments on SCIVR in the previously mentioned references and references therein.

In line with our previous observations, we assume that among the trajectories defined above, type-FA trajectories are those for which H 2 is nearly perpendicular to the Pd surface at the instant of the rebound. For simplicity's sake, we also assume that along these paths, there is no stiring effect reorienting H 2 on way to the surface; the approaching H 2 diatom is thus supposed to be already oriented nearly perpendicular to the surface before experiencing the attractive forces. In the following, we define type-FA paths as those for which θ s i belongs to the two ranges of width ∆ respectively centered at 0 and π (since θ s i belongs to the range

[0, 2π], the range of width ∆ "centered" at 0 is [0, ∆/2] ∪ [2π -∆/2, 2π]
). The interaction potential met by type FA-paths, assumed to have the shape of the red curve in Fig. 3, is represented in terms of Z in Fig. 4 (also in red). In a preliminary step, however, we will consider the simple model of interaction potential corresponding to the blue dashed curve.

This potential involves the same strongly repulsive wall but no well. We emphasize that 11 1.3. PUBLICATIONS the interaction potential felt by type-FA trajectories is the same for all of them, i.e., the potential does not depend on θ s within the two previously defined ranges of width ∆. The rest of the trajectories (type-PA plus type-NA, NA meaning non adiabatic) are defined by angular ranges which are not necessary to specify. Each type of path is assumed to remain of the same type when modifying its initial vibrational phase q s i within [0, 2π]. This assumption was validated from trajectory calculations. At last, since our purpose here is to reweight reflected trajectories, we assume that no reactive trajectory leading to sticking is involved in the present process model.

We now come back to Eq. ( 11) and focus on the contribution S F A v f j f to S v f j f due to type-FA paths. As previously stated, we start with the blue dashed interaction potential in Fig. 4. Along type-FA paths, the internal motion of H 2 is conserved, i.e., ẋ = J = 0 (it was numerically checked in a previous study 35 that ẋ = 0 along both type-FA and type-PA paths). Consequently, x f = x i = 0 and J f = J i = 0. Second, the time dependence of Z along type-FA paths is identical for all these paths (they are supposed to bounce off the 12 CHAPTER 1. THE SCATTERING OF H 2 ON PD (111) wall in the same manner since the latter is isotropic within the previously defined ranges of width ∆); from Eq. ( 12), we thus see that Ω f i takes the same value Ω F A f i for all fully adiabatic paths. Third, one may show that 17) to (19) in Ref. 47 ). 11) over the initial angles corresponding to type-FA paths (we recall that q s i belongs to [0, 2π] while

∂q s f ∂q s i = ∂θ s f ∂θ s i = 1 (see Eqs. (
θ s i belongs to [0, ∆/2] ∪ [π -∆/2, π + ∆/2] ∪ [2π -∆/2, 2π]
). After some mathematical steps, one finds:

S F A v f j f = δ 0v f sin( ∆j f 2 ) πj f [1 + e -iπj f ]e iΩ F A f i (15)
(δ 0v f is Kronecker delta). One immediately sees that according to SCIVR, type-FA trajectories do not contribute to the vibrational excitation of H 2 , in agreement with the classical description. On the other hand, they can contribute to its rotational excitation, against common sense. In the present case, rotational excitation is nothing else than the consequence of a diffraction process (see Fig. (3) in Ref. 35 ). One notes that semiclassical theory respects the fact that for homonuclear diatomic molecules, odd transitions are prohibited.

Along type-PA trajectories, the interaction potential is supposed to depend on θ. The rotational motion is thus perturbed throughout the rebound against the strongly repulsive wall. On the other hand, the vibrational motion is nearly uncoupled with the rotational motion so we can neglect its perturbation. Consequently, ẋ = 0,

x f = x i = 0, ∂q s f ∂q s i = 1 and ∂q s f ∂θ s i = 0.
Integrating the RHS of Eq. ( 11) over q s i as before and using the stationary phase approximation 45,47,55,56 to integrate over θ s i leads to the contribution to S v f j f due to type-PA paths:

S P A v f j f = δ 0v f 1 √ 2π P aths ∂J f ∂θ s i -1/2 e iΩ P A f i (16)
(it is unnecessary to specify Ω P A f i which will disappear in the following). The sum is over those paths complying with the boundary conditions J f = j f . 13 1.3. PUBLICATIONS

52

For type-NA trajectories, one still uses the stationary phase approximation 45,55,56 to integrate the RHS of Eq. ( 11) over both angles corresponding to type-NA paths. The final result provides the contribution to S v f j f due to these paths:

S N A v f j f = 1 2π P aths ∂x f ∂q s i ∂J f ∂θ s i - ∂x f ∂θ s i ∂J f ∂q s i -1/2 e iΩ f i . ( 17
)
The sum is over those paths complying with the boundary conditions

x f = v f and J f = j f .
The rotational term in the RHS of Eq. ( 16) and the RHS of Eq. ( 17) are standard expressions of classical S-matrix theory. 45,56 When (i) squaring the modulous of the full S-matrix element

S F A v f j f +S P A v f j f +S N A v f j f , (
ii) using the results of Refs. 57,58 to make delta distributions appear, and (iii) neglecting interference terms generally quenched in realistic situations by summing over non measurable quantum numbers (see Sec. 2.2 in Ref. 59 ), the whole population of state (v f , j f ) appears to be given by P

F A v f j f + P P A v f j f + P N A v f j f with P F A v f j f = 2δ 0v f [1 + (-1) j f ] sin( ∆j f 2 ) πj f 2 , (18)
P P A v f j f = δ 0v f 2π P A paths dθ s i δ(J f -j f ) (19)
and

P N A v f j f = 1 4π 2 N A paths dq s i dθ s i δ(x f -v f )δ(J f -j f ). (20)
We now consider the realistic red interaction potential in Fig. 4. Unlike the previous case, the diffraction created by the rebound against the isotropic repulsive wall may rotationally excite H 2 within the well to such an extent that H 2 has not enough translational energy along the Z coordinate to come back to the vacuum. If so, H 2 is trapped within the well.

This is what we previously called diffraction-mediated trapping. 35 If the collision energy E i is very small, j f can only be 0. On the other hand, diffraction may excite many rotational 14 CHAPTER 1. THE SCATTERING OF H 2 ON PD(111)

states within the well where the available energy is much larger than in the vaccum. Hence, the probability that the j f = 0 state is excited is small. Since this probability is also the probability that H 2 turns back to the vacuum, diffraction-mediated trapping (DMT) is expected to be important. As in previous works, we assume it to be total, i.e., all the diffracted wave is trapped. If so, the latter eventually contributes to both sticking and inelastic bouncing, just as the rest of the incident wave. The classical transposition of this DMT scenario is that the probability carried by adiabatic paths is redistributed to the remaining paths.

For simplicity's sake, we assume that the redistribution is democratic, so one may ignore adiabatic paths and renormalize to unity the probability carried by the remaining trajectories. [35][36][37] Moreover, Eq. ( 19) suggests that type-PA paths should be assigned thin Gaussian

weights G(J f -j f ) = exp[-(J f -j f ) 2 / 2 ]/( √ π
) in order to mimic the delta function (see Eq. ( 11) in Ref.; 29 is generally taken at 0.06, thus making the Gaussian 10% wide). Finally, Eq. ( 20) suggests that type-NA paths should be assigned the product of Gaussian weights

G(x f -v f )G(J f -j f )
. Now, if all the dimensions of the process are taken into account, the weights for type-PA paths will depend on the final diffractional actions a n f and a m f in addition to J f whereas the weights for type-NA paths will depend on all the final actions, as in our previous calculations (see Part I). 37 For type-PA trajectories, we chose the product

G(J f -j f )G(a n f -n f )G(a m f -m f )
of three Gaussian weights, each one depending on a single action, while for type-NA paths, we chose the same 1GB weights as in Part I, 37 given by Eqs. ( 5) and ( 6) of this work. These choices have been made so as to maximize the accuracy of our calculations at the thresholds of new quantum states given the number of trajectory results available.

To recapitulate, the sticking probability obtained when dividing in four sets the whole bunch of trajectories run [(i) reactive, (ii) type-FA, (iii) type-PA and (iv) type-NA)] reads:

P II S = N S N S + Σ N A R + Σ P A R (21) 15 1.3.

PUBLICATIONS

where

Σ N A R = q Σ N A R (q) (22)
with

Σ N A R (q) = N N A k=1 w 1GB (a k , q), (23)
and

Σ P A R = q Σ P A R (q) (24)
with

Σ P A R (q) = N P A k=1 w GB (a k , q), (25)
w GB (a, q) = G(J f -j f )G(a n f -n f )G(a m f -m f )w SB (a, q) (26)
and

G(u) = exp[-u 2 / 2 ]/( √ π ). (27)
As previously, N S is the number of trajectories leading to sticking and type-FA paths are ignored. Moreover, N N A and N P A are the numbers of type-NA and type-PA trajectories.

In Eqs. ( 23) and ( 25), the sums are obviously over the type-NA and type-PA trajectories, respectively. Finally, the state-resolved reflection probability reads:

P II R (q) = Σ N A R (q) + Σ P A R (q) N S + Σ N A R + Σ P A R . (28)
The type of calculations published in Part I (or using Eqs. ( 7)- (10)

Results and discussion

The classical trajectory calculations were performed by means of the Density-Functional-

Theory-based ground state potential energy surface of Busnengo et al. 60,61 The Bulirsch-Stoer integrator was used with an initial time step of 0.01 fs. Trajectories were started at Z = 5 Å with normal incidence. Reflected paths were stopped at the same value of Z and reactive trajectories at r = 2.5 Å. For H 2 in the initial state (v i = 0, j i = 1), a total of 133 energies spaced 3 meV aparts were considered within the range 2 -398 meV, and 800 000 trajectories were run at each energy. For the rest of the initial states, 50 energies were considered within the range 2 -149 meV with the same energy spacing and number of trajectories per energy as previously. A few calculations beyond 149 meV showed that the differences between quantum and classical results are no longer significant. 62 In any cases, the Gaussian width parameter was taken at 0.003 and 0.06 for the calculation of 1GB and GB weights, respectively. Fig. 5 shows the quantum time-dependent wave-packet (Q-TDWP) 40,62 and purely classical (QCT) sticking probabilities for H 2 in the initial state (v i = 0, j i = 1). Clearly the QCT method fails to correctly reproduce the uneven structure of the quantum probability below ∼270 meV. This behavior is similar to that observed when H 2 is in the initial state (0, 0). 63 On the other hand, it is shown in panel (a) of Fig. 6 that the QCT-SC-I and QCT-SC-II Sticking probability curve in panel (a)) is slightly lower than the QCT-SC-I probability (red curve). This is due to the fact that type-PA reflected trajectories are not ignored in QCT-SC-II calculations while they are in QCT-SC-I calculations. Since in the present case, the disagreement with quantum mechanics is increased by the account of type-PA trajectories, these are likely to be more rotationally adiabatic and, thus, more fully adiabatic than expected from our criterion to distinguish between type-FA and type-PA trajectories. This suggests that this criterion should be refined. We note that the larger the collision energy, the stronger the couplings between translational, vibrational and rotational motions within the interaction region, the smaller the amount of adiabatic trajectories and the closer the green and red curves (whenever there are no adiabatic paths, QCT-SC-I (Eqs. ( 7)-( 10)) and QCT-SC-II (Eqs. ( 21)-( 28))

j i = 1 Q-TDWP QCT
methods lead to strictly the same predictions).

For higher initial rotational excitations (j i = 2 -4), sticking probabilities are displayed in Fig. 7 while rotationally elastic and inelastic reflection probabilities are represented in Fig. 8. The agreement between the QCT predictions in a quantum spirit and the quantum results appears to be even better than for j i = 0 (see Part I 37 and below) and j i = 1. On the other hand, standard QCT calculations make disappear any structure in the probabilities, as for j i = 0 and 1. The drops in the sticking probabilities at 98 meV for j i = 2, and 125 meV for j i = 3 (see panels (a) and (b) of Fig. 7), correspond to the openings of the rotational excitation channels 2 → 4 and 3 → 5, as clearly seen in panels (a) and (b) of Fig. 8. In both cases the thresholds are very well described.

Finally, we applied our new QCT-SC-II approach to the case (v i = 0, j i = 0), using the same energy spacing and trajectory amount per energy as in Part I. We recall that among adiabatic paths, type-FA paths were identified as those leading to Jf = 0 for E i lower than ∼90 meV (see the discussion around Fig. 1). The resulting sticking probability is represented by the green triangles in Fig. 9. The red circles in the same figure were obtained by means of the QCT-SC-I approach. Above ∼ 40 meV, the green sticking probability is slightly lower than the red one, as for j i = 1, since the former takes into account more reflected trajectories than the latter. We note that above ∼ 200 meV, QCT-SC-II probabilities are in slightly better agreement with QM probabilities than QCT-SC-I ones. On the other hand, it is quite probable that below 40 meV, most adiabatic paths are type-FA paths. In such case, the two formulations lead to the same predictions (compare Eqs. 7-10 with Eqs. 21-28), and the red and green curves overlap, as observed in Fig. 9.

j i = 3 j f = 3 j f = 1 j f = 5 0 50 100 150 (c) j i = 4 j f = 2 j f = 0 j f = 4 Q-TDWP QCT-SC-I QCT-SC-II
Although in the case of the H 2 /Pd(111) system, the QCT-SC-II method leads to predictions very similar to the QCT-SC-I method, the former is based on a more rigorous application of semiclassical scattering theory than the latter. Moreover, the calculations are 20 CHAPTER 1. THE SCATTERING OF H 2 ON PD(111)

of the same difficulty (or simplicity) for both approaches. We might thus be more inclined to use the QCT-SC-II method in future studies.

Conclusions

The classical approach in a quantum spirit, initially devised for gas-phase reactions, differs from the standard quasiclassical trajectory (QCT) method by the statistical weights finally assigned to the trajectories. While these weights are identical for all the trajectories in the QCT method, they are either Gaussian, 28,29,32,33,38 or simply zero 32,[34][35][36] for the previous approach. Gaussian weights take into account Bohr's condition of quantization for the final bounded degrees-of-freedom. Zero weights take into account the fact that, within the quantum picture, some parts of the incident wave may be diffracted at very small collision energies. The first type of weighting is called Gaussian binning (GB 28,29,32 or 1GB 33,38 ) while the second is called adiabaticity correction (AC). 32,[34][35][36] In a recent work (Part I), 37 we applied the 1GB-AC-QCT method to the heterogeneous gas-surface reaction between H 2 in its internal ground state and Pd(111) surface. The sticking probability and the state-resolved reflection probabilities obtained by means of this approach were found to be in close agreement with the corresponding probabilities obtained from exact 21 1.3. PUBLICATIONS quantum scattering calculations. In contrast, the standard QCT method fails to reproduce the previous quantities, except in the high energy limit (as it should be).

The purpose of this work was to check whether the previous agreement was maintained when H 2 is rotationally excited. We not only applied the method of Part I, but also a variation of it regarding the adiabaticity correction, put on firmer semiclassical grounds than the method of Part I. The results displayed in Figs. 6789, of similar quality for both methods, provide a very positive answer. The rugged collision energy dependences of quantum sticking probabilities are accurately reproduced for all the rotational excitations considered. This comment holds for the state-resolved reflection probabilities, particularly around the opening of inelastic channels. As a matter of fact, for non activated reactions at relatively low energies, our two-part study shows that the classical approach in a quantum spirit offers a remarkable description of gas-surface reaction dynamics features that the QCT method in its standard implementation is unable to provide.

Summary

For the first time, the quantization in Bohr's spirit of all the relevant final actions of a classically described reactive molecular system has been performed. The process studied is the scattering of H 2 on Pd(111). Simulations have been carried out for the H 2 in its internal ground state as well as for initial rotational excitations. The action-based Gaussian binning (GB) and energy-based Gaussian binning (1GB) have been employed to take into account Bohr's quantization principle in the analysis of the final trajectory results. The final actions to be quantized are the vibrational and rotational actions of the reflected H 2 diatom, as well as the diffractive actions (classical analogs of Miller indices). On the other hand, the adiabaticity correction (AC) takes into account the fact that, within the quantum picture, some parts of the incident wave may be diffracted at very small collision energies.

For this system, the quasi-classical trajectory method fails to describe the sticking probability curve at low energies. Both combinations of binning procedures (GB and 1GB) with the AC lead to results in a better agreement with quantum timedependent wave packet (Q-TDWP) results. In particular, the use of 1GB shows a remarkably accurate description of the zones near the energy thresholds of the opening of new channels in the state-resolved reflection probability curves. At these energies, the suitability of the use of 1GB compared to GB is demonstrated.

Also, a variation of the AC based on firmer semiclassical grounds has been used

for the same study system. The results are very similar and the quantum sticking and state-resolved reflection probabilities are accurately reproduced. Computationally this variation is not heavier than the traditional version. In general, the prospects look bright but more information is needed before making the decision of which AC should be used.

It is planned to repeat this study for different gas-phase and gas-surface reactions to figure out whether the quality of the present agreement between quantum corrected QCT and exact quantum probabilities are systematically maintained or not. If it is indeed confirmed that the accuracy is equally good for a large variety of systems and dynamical situations, the current methodology can be a reliable alternative to computationally heavy quantum dynamics calculations. It may be also applicable to systems of increasing complexity and/or energy regimes not currently

Introduction

One of the biggest challenges for the human race in the XXI century is to make the transition to clean and renewable energy sources. Atomic fusion is among the main candidates to accomplish this task. Nowadays, the largest and most expensive initiative to demonstrate the feasibility of fusion energy is ITER, the facility located in southern France as part of an international collaboration. ITER is a 16-billioneuros project that currently dictates most of the material specifications and loading scenarios that are considered to be relevant for fusion facilities. The main reason for the persistent delay in the fusion calendar is the gigantic technical effort involved.

The design, construction, and operation of a fusion power plant still pose plenty of problems that require coordinated advances in many different disciplines and technologies, including materials research [40].

Amidst the multiple problematics under investigation, our work focus on the modifications of the interaction of H 2 with a W surface when oxygen has been previously adsorbed. Questions about changes in the properties of the materials after processes such as adsorption, absorption, diffusion, or recombination still do not have clear answers. This specific system is of interest in the fusion community due to the use of W in several plasma-facing components.

The first step to tackle this problem is to understand the dynamics of the interaction of H 2 with a clean W surface. Experimental studies using molecular beams [41] show that at normal incidence, the measured values of the initial sticking probability S 0 increase monotonically. Of particular interest is the fact that this magnitude has a strong dependence on the crystallographic face of the surface. In general, values measured for the W(110) surface are smaller than those measured in the W(100)

and W(111) surfaces. Similar results have been observed for incident HD [42] and N 2 [43, 44] molecules. Other processes like hot atom and Eley-Rideal recombination of hydrogen and its isotopes on W surfaces have also been studied [45-51].

Theoretical studies carried by Busnengo and Martínez [52] reproduce the higher reactivity of W(100) as compared with W(110), as well as the general trends of the sticking coefficient in both faces. This study uses classical dynamics on a sixdimensional potential energy surface (PES) built from density functional theory (DFT) calculations.

Oxygen atoms adsorbed in a W(110) surface form an over-layer on top of the threefold sites as shown by experimental studies [53-55]. This is also the preferential adsorption site for hydrogen atoms [56, 57]. In their experimental studies, Whitten and Gomer observed that the amount of H adsorbed in the W(110) plane decreased with the amount of O preadsorbed [58]. For ratios O/W 0.35, hydrogen adsorption drops to zero.

For the past thirty years, ab initio molecular dynamics (AIMD) has been gaining traction in the field of molecule-surface interactions [59, 60]. This method has proven to be effective and adaptable in the treatment of the dynamics of atoms and molecules at surfaces. In full AIMD calculations, the adiabatic forces between the atoms are calculated on the fly at every step of the numerical integration. The most common implementations use plane-wave-based density functional theory (DFT) [61-63] techniques to evaluate the energy and gradient. Even with the accelerated development in computer systems, one of the major drawbacks of AIMD is its computational cost. This makes it hard to employ binning procedures like GB and 1GB that require a big number of trajectories to achieve statistical accuracy.

Ever since AIMD calculations have been used in surface chemistry, a main focus has been on the dissociative adsorption (mostly on metal surfaces) of H 2 and other diatomic molecules [20, 64, 65]. In this chapter, we employ AIMD calculations to study the following problems:

• The adsorption mechanisms of H 2 on a clean W(110) surface [66] in Section 2.3.1. Van der Waals dispersion effects in the interaction between molecule and surface are included and discussed in this section as well.

• The adsorption mechanisms of H 2 on oxidized W(110) surfaces in Section 2.3.2.

Calculation details

In a similar fashion to Chapter 1, the AIMD can be divided into three steps:

• Initial conditions are generated.

• Classic equations of motion are integrated until certain conditions are met (the atom(s) or molecule(s) reach one of the exit channels or the maximum time of simulation is reached).

• The observables are calculated through statistical analysis.

In this chapter, we will study the dynamics of H 2 molecules when approaching clean and oxidized W(110) surfaces. All calculations are based on density functional theory (DFT) using the Vienna Ab initio Simulation Package (VASP) [9, 67-69].

Within this code, ionic relaxations are initially performed to obtain the equilibrium positions of the surface atoms as well as the equilibrium interatomic distance of the H 2 molecule in the vacuum. These values will be used as part of the initial coordinates for each trajectory (Subsection 2.2.1). Once the dynamics start, DFT is also used inside the code to calculate the electronic ground state of the molecule-surface system at every integration step. This result is then used to calculate the forces employed to integrate the equations of motion (Subsection 2.2.2).

The current section provides general information about the methodology that we use and approaches that are common to all systems studied. Specific details about VASP calculations are given in Subsections 2.3.1 and 2.3.2 for each particular surface.

The coordinate system employed through this work is depicted in Fig. 2.1 (a).

Initial conditions

In Chapter 1 all our calculations were done considering a frozen surface. This approximation places the surface atoms in their equilibrium positions and does not allow them to move during the dynamics. In our AIMD calculations, on the other hand, surface atoms are still placed initially at their equilibrium positions but the movement is allowed during the dynamics. The latter means that some energy transfer between the incident molecule and the surface is in principle included in our calculation. The first step to find the initial equilibrium positions is to perform ionic relaxation calculations for the bulk of the material. In our case, the W crystal structure presents a body-centered cubic (bcc) arrangement of atoms that will be infinitely reproduced in three dimensions (3D). This permits VASP calculations that use periodic boundary conditions. In this step, we obtain the lattice constant a.

When performing surface calculations the super-cell approach is used to ensure the 3D periodicity of the system. First, a 3D cell (supercell) composed of finite atomic layers representing the surface(crystal slab) and vacuum is constructed. This is the base unit that will be periodically repeated in all directions as shown in Fig. 2.2.

The shape and dimensions of the supercell depend on each particular system and the type of calculations. In general, for collisions of a molecule with a surface, the dimensions of the supercell need to be large enough to avoid interactions between the periodic images of the incident molecule. The number of atomic layers of the slab needs to be large enough to reproduce the electronic properties of the surface.

On the other hand, extra layers mean more atoms in the slab and this entails longer calculation times. We usually conduct different tests to choose an adequate size of the supercell while keeping the calculation times under reasonable limits. In Sub- The initial position of the center of mass in the z-direction is chosen as the middle point between two consecutive slabs to minimize interaction with the surface.

The influence of the surface temperature has been ignored and therefore the initial momentum of surface atoms is zero. Initially, the H 2 molecule is in its ground rovibrational state. This fact together with ignoring the ZPE makes the initial values of the momentum associated with r, θ, and φ equal to zero as well. Only normal incidence has been studied in our calculations making:

P i X = 0; P i Y = 0; P i Z = -2E i M (2.1)
with P i X , P i Y , and P i Z being the components of the initial linear momentum and M the total mass of the H 2 molecule. E i is the initial kinetic energy of the incident molecule.

Evolution in time

In the scope of the BOA, only the movement of the nuclei is propagated in the potential determined by the electronic distribution. In our simulations, the incoming molecule atoms as well as the surface atoms are all allowed to move. The equations to integrate are Newton's equations of motion:

F i (t) = m i d 2 r i (t) dt 2 (2.2)
where m i and r i are the mass and position of the ith atom (surface and incident molecule atoms). The force F i (t) acting on the ith atom is calculated from the electronic density as mentioned before. These equations are integrated within VASP using a Verlet algorithm [70]. The integration time step used in all our simulations is 2 fs. AIMD calculations using VASP require as well as input the number of steps in the time integration, NSW. An auxiliary script can be used to run calculations for a fixed amount of steps NSW and check at the stopping point if molecules have reached one of the final channels mentioned in Subsection 1.2.2. If that is the case, the simulation for that specific trajectory is stopped. If not, a new calculation is started for another NSW steps from the final states of the previous calculation. This is done to save computational resources. In this work, the maximum time of integration for a given trajectory is typically set at 2 ps.

Statistical analysis

Due to the high computational cost of the methods used in this chapter, only 400 trajectories have been calculated for each collision energy per surface. This makes it impossible to use the GB procedures from Chapter 1. The sticking coefficient S 0 is obtained using SB through:

S 0 = N S N (2.3)
where N S is the number of trajectories ending in dissociative adsorption and N is the total number of trajectories. Further details about the exit channels in our AIMD calculations are similar to those explained in Chapter 1.

Publications

increase in the far-distance attraction that is compensated by a stronger repulsion at short distances. The combination of both effects appreciably decreases the value of the dissociation probability. The successful comparison of our results with experimental information confirms that the methodology employed can be considered as a rich and accurate instrument to study the dissociation of hydrogen on surfaces.

Introduction

In general, scientific and technological progress has very much improved life conditions for humankind. Scientific advance and the subsequent social and economical development have succeeded in making human life longer, healthier, and significantly dignified. The current challenge is to make this development sustainable and truly global. Among many other tasks, there is a continuous need of finding alternative energy sources that are efficient, sustainable, and environment friendly. Population growth and economic activity will very probably increase the energy demand worlwide in the decades to come.

For a long time, fusion power has been a promising candidate to become a carbon-free, sustainable, long-term energy source. Fusion power is conceptually based on the generation of electricity from the heat released in nuclear fusion reactions. The most common reaction used for this purpose is the fusion of hydrogen isotopes to obtain helium and neutrons.

One of the most challenging engineering problems in fusion reactors is the design and maintenance of the vacuum vessel in which the plasma moves. 1 The fusion process requires to confine the hydrogen isotopes, which act as fuel, under extreme temperature and pressure conditions in order to create a hot plasma where fusion reactions can occur. The enormous quantity of heat generated during the reaction makes of tungsten, a robust metal with the highest melting point, the standard material of choice to build many of the plasma-facing components.

Among plasma-facing components, the most vulnerable one is probably the divertor.

The divertor is instrumental to extract heat, impurities, and ash from the plasma, as well 2

CHAPTER 2. AB-INITIO MOLECULAR DYNAMICS OF HYDROGEN ON TUNGSTEN SURFACES

as to protect parts of the plasma-facing walls. It receives a high flux (∼ 10 20 cm -2 s -1 ) of relatively low-energy (∼ 100 eV) hydrogen and helium particles. The huge flux of hydrogen isotopes and helium towards the divertor can deteriorate the W properties. Surface blistering and exfoliation, among other problems, can appear due to the trapping and diffusion of these atoms within the material. Clustering of helium and hydrogen can lead to bubble formation and migration of these bubbles to the surface induces the blistering. In this scenario, hydrogen recombination at the surface and molecular desorption is possible as well. The eventual fate of these desorbed molecules, i.e., whether they dissociate again at the surface or they return to the plasma, is relevant for the full process: the injection of atomic and molecular hydrogen from the walls again into the plasma can significantly perturb its ideal behavior. Therefore, for a correct operation of the global system, it is essential to understand and accurately model the interaction of hydrogen and its isotopes with W. The role of the W surface as an agent that either promotes or inhibits some of the elementary reactive processes leading to adsorption, recombination, diffusion, or desorption is crucial in this respect. Current theoretical and experimental activity on this topic is therefore extensive.

The traditional surface science approach, in which reactions take place under wellcontrolled conditions, shows however that the interaction between hydrogen and tungsten surfaces is a complex problem. Small differences in temperature, pressure, or crystalline face of the surface can dramatically change the output of a reaction. A remarkable feature in the adsorption of H 2 on W is, for instance, the acute dependence on the crystalographic face of the surface, as proven in molecular beam experiments. For initial energies below 1 eV and normal incidence, adsorption of H 2 can be as much as a factor of two higher in W(100) than in W(110). 2 This effect is even more drastic when the incident molecule is N 2

where the difference between the sticking coefficient in W(100) and W(110) is two orders of magnitude. 3,4 The disparity can be understood from the differences in the long-distance interaction between the incident molecules and the surface, as well as from the associated 3

PUBLICATIONS

phase space to access the precursor well where dynamic trapping takes place. 5,6 A similar explanation in terms of long-distance interaction and dynamic trapping can be invoked to explain the differences in the chemisorption of H 2 on W(100) and W(110), as proven by Busnengo et al. 7 In the particular case of H 2 on W(110), which is the one treated in this work, molecular beams experiments show that the initial sticking probability S 0 at normal incidence increases monotonically from S 0 ∼ 0.1 to S 0 ∼ 0.3 for kinetic energies of the beam below 350 meV. From the theoretical point of view, Busnengo and Martínez performed classical dynamics for H 2 on W(110) and W(100) using a six-dimensional potential energy surface (PES) built from density functional theory (DFT) calculations. 7 Their results reproduce the higher reactivity of W(100) as compared with W(110), as well as the general trends of the sticking coefficient in both faces. They show that the difference in reactivity between W(100) and W(110) can be understood in terms of an indirect channel, in which dynamic trapping plays a role, that is significantly reduced in W(110). There is also a large theoretical activity on the study of hot atom and Eley-Rideal recombination processes of hydrogen and its isotopes on W surfaces. [9][10][11][12][13][14][15][16] Absorption and diffusion of hydrogen on W surfaces have been analyzed as well. 17 Our goal in this work is to provide additional insights on the adsorption mechanisms of H 2 on W(110) and improve the accuracy of the currently existing theoretical results.

For this purpose, we use ab initio molecular dynamics (AIMD), which is an effective and adaptable method to treat the dynamics of atoms and molecules at surfaces. In AIMD, the adiabatic forces between the atoms are calculated on the fly, at every step of the numerical integration. The value of the electronic energy is obtained in our case from DFT. In some of our calculations, we include van der Waals (vdW) contributions to the exchange correlation 4 The improvement in the theoretical description of the H 2 /W adsorption is an objective per se, but it will also set the ground for further analysis of other systems of interest. One of the questions that is under study in the fusion community is the modification that the presence of oxygen at the W surface introduces in processes such as adsorption, absorption, diffusion, or recombination. The results obtained in this work will serve as a reference in a future comparison between adsorption features in clean and oxidized surfaces. The methodology employed, proven successful here, can be equally applied without much additional computational effort to the case of oxidized W surfaces.

The outline of the paper is as follows. The details of the theoretical procedure and the numerical implementation are described in Sec. 2. Our results and the comparison with available experimental data are shown in Sec. 3. Our main conclusions are summarized in Sec. 4.

Theory and Methods

AIMD calculations within the Born-Oppenheimer approximation (BOA) have been performed using the Vienna Ab initio Simulation Package (VASP). 18 These calculations are computationally expensive. Parameters have been chosen looking for an adequate balance between calculation times and accuracy. A statistically significant number of trajectories has been simulated using either the Perdew-Burke-Ernzerhof xc-functional (PBE) or the vdW xc-functional developed by Dion et al. (vdW-DF). 19 We choose this specific form of the exchange-correlation term because it has proven to provide satisfactory results in previous calculations of the sticking coefficient for systems with, in principle, similar features in the energy landscape. It was shown to be among the best choices to describe the dissociation 5

PUBLICATIONS

probabilities of H 2 on Ru(0001). 20 It also describes reasonably well the dissociation probabilities of N 2 on W(110) at low energies and normal incidence. 21 In all calculations, the electron-core interaction is treated within the plane augmented wave (PAW) approximation through the PBE pseudopotentials provided with VASP. 22,23 AIMD calculations in VASP using the vdW-DF xc-functional are carried using the routine written and provided by Jiri

Klimeš 24 based on the Román-Pérez and Soler algorithm. 25 The energy cutoff used in all of the calculations was 250 eV, the recommended maximum value for the pseudopotentials we are using. Tests using higher cutoff energies show that the difference in the interaction energy between the H 2 molecule and the surface is below 2%. In Fig. 1 the coordinate system employed is depicted. X cm , Y cm and Z cm represent the coordinates of the center of mass of the H 2 molecule while r, θ and φ account for the internuclear distance, the polar angle, and the azimuthal angle, respectively. Different lattice constant a = 3.17 meV in the case of the ones using the vdW-DF xc-functional. In this energy range theoretical and experimental results are available 2,7 and a comparison with these works will be presented in the next section. The integration step for all trajectories is 2 fs and the energy convergence error for electronic relaxation is 10 -6 eV. The dissociation probabilities are calculated from 400 trajectories per energy. The criteria for the dissociation of the H 2 molecule is established as a separation in the internuclear distance between H atoms of r > 1.5 Å and dr/dt > 0.

When the H 2 molecule is moving away from the surface and the distance between molecule and surface is Z cm > 4.5 Å, a reflection event is considered.

Diverse theoretical data about the chemisorption of H 2 on W(100) and W(110) surfaces have been obtained previously through molecular dynamics (MD) calculations. 7 The six dimensional PES employed in this calculations by Busnengo et al. 7 was built using the corrugation reducing procedure (CRP) 28 to interpolate energy values previously calculated with DFT. The PW91 xc-functional 29 was used in the scope of the generalized gradient approximation (GGA) to get these energies. Classical (C) and quasi-classical (QC) calculations were performed within the MD approach. The initial zero point energy (ZPE) of the former is not taken into account while the quantum ZPE of H 2 in vacuum is used for the latter. In the following, MD-C and MD-QC will be employed to refer to the results of these calculations. There are however some features absent in our calculation that may explain the still existing quantitative difference with the experimental results. First of all, we do include the motion of the W atoms in the dynamics but we do not fix any initial surface temperature, while the experimental data shown in Fig. 2 possible improvements for our calculation include quantum effects in the dynamics 30,31 and electronic non-adiabatic effects. 32,33 In surface, at roughly Z cm = 2 -2.5 Å. Discussions about this behavior can be also found in Ref.

Results and discussion

7 and is similar to what happens with N 2 on the same surface. 5,6 For both xc-functionals this effect looks the same, although the dissociation process is more direct for the PBE functional.

However, the distance at which adsorption is effectively determined depends on the initial collision energy, as can be seen in Fig. 5. The latter point can be explained if we consider that molecules with higher initial energy can get closer to the surface as it allows them to surpass to certain extent the potential energy of the surface on the repulsive zones of the PES. Below 150 meV all molecules that reach 2.5 Å will get dissociated while for higher energies that distance is reduced to approximately 2 Å.

In Figs. 2 and5 it is noticeable that the probability of sticking and the probability of H 2 molecules reaching distances close to the surface is always higher when using the PBE xc-functional than when using the vdW-DF xc-functional. Even if the adsorption is determined far away from the surface, the subsequent dynamics of each trajectory may vary. In the case of the PBE functional, most of the dissociation takes place through a fast and direct process. In the case of the vdW functional, however, most of the molecules that are able to cross the bottleneck of Z cm = 2 -2.5 Å spend some time hovering and bouncing on the surface before becoming eventually dissociated. This dynamical process includes the motion of the molecules looking for a more suitable position and/or angular configuration that allow them to approach the surface and finally dissociate.

The time delay in the dissociation introduced by this dynamic trapping is longer for the vdW-DF xc-functional than for the PBE functional. In average, trajectories leading to sticking experience 0.9 rebounds before dissociation when using PBE while the number of rebounds raise to 1.8 when the vdW-DF xc-functional is employed.

Our calculations show that the preferential path for the H 2 molecules to approach the W(110) surface is above the top sites. This is shown in Fig. Z cm coordinate is closer to the one displayed. Through cuts at different distances from the surface it can be observed a movement toward the top positions. This trend is independent of the xc-functional employed, although it seems that the higher concentration near these spots happens around 2.1 Å when vdW-DF is used and 1.7 Å when using PBE. This general behavior is consistent with previous theoretical works. 7 Apart from the positions in the XY plane, the angular configuration also plays an important role when deciding which molecules are undergoing adsorption on the surface. Figure 8 show the distribution per molecule of the polar angle θ at the same distances from the surface than in Fig. 7. For both xc-functionals is easy to spot the predominance of angles close to 90 • . The distribution of the azimuthal angle φ on the contrary does not present any remarkable features. It needs to be mentioned that for trajectories that spend a considerable time near the surface (most of them when using vdW-DF) the data from Figs. 7 and8 are collected when molecules approach the surface for the first time. It seems then that molecules 13 with axis parallel to the surface and above the top positions have the greatest chances to get closer to the surface and undergo a future dissociative adsorption. This behavior helps to explain why the dissociation is effectively decided far away from the surface. At roughly 2.1 Å for the initial kinetic energy of 200 meV, only the H 2 molecules with suitable conditions (see Figs. 7 and8) will not be reflected back to the vacuum, creating a bottleneck above top sites. In our simulations almost the totality of molecules passing this point will be dissociated hinting that, once the bottleneck is crossed, there is a very limited phase space for the molecules to be reflected back to the vacuum.

2.3. PUBLICATIONS Z cm = 9 Å PBE Z cm = 2.5 Å Z cm = 2.1 Å Z cm = 1.7 Å vdW-DF
0 • 90 • 180 • Z cm = 9 Å 0 • 90 • 180 • PBE Z cm = 2.5 Å 0 • 90 • 180 • Z cm = 2.1 Å 0 • 90 • 180 • Z cm = 1.7 Å 0 30 60 90 120 150 Z cm = 9 Å vdW-DF Z cm = 2.5 Å Z cm = 2.1 Å Z cm = 1.7 Å

Conclusions

In summary, we have used AIMD to understand the dissociation dynamics of H 2 on W(110).

We have shown that the use of the vdW-DF xc-functional improves the comparison of the theoretical description with available experimental information. Our analysis of the process shows that the dissociation is decided at a distance relatively far from the surface, in a way similar to what happens in the dissociation of H 2 on W(100) or in the dissociation of N 2 on W surfaces. Almost all hydrogen molecules able to reach a given distance from the surface (≈ 2 -2.5 Å, depending on the initial kinetic energy) are eventually dissociated. When using the vdW-DF xc-functional, we obtain a short dynamic trapping at the surface before the dissociation takes place. This is a mechanism that nevertheless has no influence in the final values of the dissociation probability.

The use of AIMD is an alternative to other theoretical approaches based on the construction of interpolated PESs. A weakness of AIMD is that it can be computationally expensive and therefore the number of trajectories that are launched can be significantly smaller than those launched using interpolated PESs, with the subsequent damage in statistical accuracy for low probability processes. It has other advantages, however, and one of them is that it can address more complex systems, often without much additional computational effort. In the case of hydrogen interacting with tungsten surfaces, a problem of current interest is the effect that the oxidation of the surface may have in the hydrogen adsorption process. This is of particular interest for the fusion community, but not only. The methodology used in this work can be easily extended to the case of oxidized W surfaces and the results presented here can be used as a benchmark to quantify the effect that the presence of oxygen atoms adsorbed on the surface may have on the dissociation dynamics of H 2 on these surfaces.

Work along these lines is in progress. show that the influence of O adsorbates on the H 2 dissociation dynamics on W(110) goes much beyond the occupancy of possible H adsorption sites.

Introduction

For decades, the generation of energy in nuclear fusion processes has been promoted as a promising enviromental-friendly alternative to current power sources. Unfortunately, the expectations have not been fulfilled yet. One of the main reasons for the permanent delay is the enormous engineering challenge associated to the construction of a fusion power plant.

A complex and integrated approach is required, including contributions from many different scientific fields.

Fusion is based on the confinement of hydrogen isotopes under extreme temperature and pressure conditions. Fusion reactions take place in the hot plasma thus created. In the development of viable fusion power plants, a crucial issue is an adequate selection of materials. 1 Tungsten has become the standard constituent of the plasma-facing walls because of its high melting point and thermal conductivity. During operation of the reactor, there is a huge flux of hydrogen isotopes towards the tungsten walls, which can become eventually deteriorated.

The role of oxygen contamination in the chemistry of W surfaces is also a frequent topic of investigacion in the context of fusion devices. Low levels of oxygen impurities are usually found in fusion devices and may trigger the oxidation of the W surface. 2 Tungsten oxides can eventually become volatile 3 and affect the plasma properties. For all these reasons, there is a large amount of research currently devoted to characterize how W systems behave when in contact with hydrogen gas. 4 Interest in the interaction between hydrogen and metal surfaces has been historically associated to heterogenous catalysis, the basic mechanism for many large scale chemical industry processes. 5,6 Hydrogen adsorption on W is one of the simplest chemical reactions that one may envision on a surface. It has been studied for more than one hundred years, 2

PUBLICATIONS

since Irving Langmuir placed a filament of tungsten in a vessel containing hydrogen gas. 7 In spite of the extensive coverage of the problem, it is only recently that a basic atomiclevel understanding of the mechanisms ruling different reactive processes of hydrogen on W surfaces has been reached. [8][9][10][11][12] Under well-controled conditions, the interaction between hydrogen molecules and W surfaces is quite subtle. Small differences in temperature, pressure, or surface termination can drastically change the chemical products. For initial energies below 1 eV and normal incidence, the adsorption of H 2 is much higher in W(100) than in W(110). 13 This effect is not exclusive of H 2 . For N 2 , the difference between the sticking coefficient in W(100) and W(110) reaches up to two orders of magnitude. 8,14 The explanation of the different reactivity between the two faces has been explained, both for H 2 and N 2 , in terms of long-distance interaction and dynamic trapping. 10,11,15 We recently published a detailed theoretical analysis of the dissociative adsorption of H 2 on W(110) based on ab-initio molecular dynamics (AIMD). 12 Our AIMD methodology relies on density functional theory (DFT) and classical molecular dynamics. For normal incidence, molecular beams experiments show that there is a monotonical increase of the initial sticking probability S 0 on clean W(110), from S 0 ∼ 0.1 to S 0 ∼ 0.3 when tuning the kinetic energy of the H 2 beam up to 350 meV. 13 Our theoretical results reproduce reasonably well the experimental data and show that the fate of each trajectory is determined at distances relatively far from the surface. 12 An accurate description of the long-range interaction between the molecule and the surface therefore requires the inclusion of van der Waals (vdW) terms.

One of the advantages of AIMD is that the complexity of the system under study can be often increased without much additional computational effort. In the particular case of the current work, we extend the methodology previously used for H 2 on W(110) to the case of H 2 impinging on oxidized W surfaces, an involved problem for which previous studies are scarce. The reliable results we obtained for the adsorption of H 2 on W(110) can be used as a reference to quantify the effect induced by the presence of adsorbed oxygen in the H 2 dissociation dynamics.

CHAPTER 2. AB-INITIO MOLECULAR DYNAMICS OF HYDROGEN ON TUNGSTEN SURFACES

The adsorption of oxygen on W(110) is experimentally well characterized. Oxygen molecules dissociate on W(110). 9,16,17 The close packed nature of the W(110) does not favor however the formation of bulk oxides, which often appear in other W faces. No experimental evidence has been found for oxygen penetration into the bulk through the W(110) face. Isolated oxygen atoms generally adsorb on three-fold coordinated sites and can form well ordered structures for coverages below one monolayer (1ML). 9,18,19 Ordered oxygen islands and disordered domains coexist for some oxygen coverages. Phase transitions can arise as well. 16,20,21 At room temperature and for 0.5 ML, 0.75 ML, and 1 ML coverages, the structures p(2x1), p(2x2), and p(1x1) are respectively observed. 16,19,22 There is not much experimental information on the effect that adsorbed oxygen atoms may have in a posterior adsorption of hydrogen on W. Using different experimental techniques, Whitten et al. showed that preadsorbed oxygen decreases the amount of hydrogen that can be adsorbed on W(110). 23 For coverages of oxygen above 0.35 ML and temperatures of 90 K, they found no significant evidence of H adsorption. 23 The latter data suggest that the role of adsorbed oxygen goes beyond the possible blocking of adsorption sites for the incoming H atoms. For O coverages below 0.35 ML, the adsorption of H is also clearly below the one found in the clean W(110) surface.

Our purpose in this work is to investigate which are the mechanisms that inhibit the dissociative adsorption of hydrogen when oxygen atoms are present at the W(110) surface.

Our approach is based on dynamical grounds. We focus into three different situations, representative of three different regimes: a clean W(110) surface, a W(110) surface with 0.25 ML O coverage, and a W(110) surface with 0.5 ML O coverage. In this way we can follow the evolution of the system from what can be considered as normal rates of H adsorption up to no H adsorption at all. We will show that O atoms adsorbed at the W surface act as repulsive centers that modify the dynamics of the approaching H 2 molecules. As a result, the surface active sites in which H 2 dissociate are shielded and become inaccessible.

The outline of the paper is as follows. The details of the theoretical procedure and the

Theory and Methods

A statistically significant number of calculations have been carried out for hydrogen molecules impinging on clean and oxidized W(110) surfaces. The coordinate system employed is depicted in Fig. 1 (a). X cm , Y cm and Z cm represent the coordinates of the center of mass of the H 2 molecule while r, θ and φ account for the internuclear distance, the polar angle of the internuclear axis, and the azimuthal angle of the internuclear axis, respectively. Oxygen atoms, when present, are not represented in this panel for clarity purposes. The Z cm distance to the surface is always measured relative to the upper layer of W atoms. The unitary cell of the W(110) surface employed in our simulations is shown in Fig. 1 (b). The black circle inside the cell represents an oxygen atom in the threefold hollow site. This is the chemisorption position found in various experiments 9,18,22 for oxygen in the W(110) surface.

The chemisorption position is independent of coverage. 9 All distances are given in function of the W lattice constant a.

Two oxidized surfaces with different coverages have been used in our simulations and are represented in Fig. 2. In panel (a) the ordered phase (2x1) at coverage Θ = 0.5 monolayers (ML) is shown. All over this article, Θ stands for coverage in place of the commonly used θ to avoid notation conflicts with the polar angle. The meaning of its value remains the same: the ratio between the number of oxygen atoms and the number of tungsten atoms in the first layer of the surface. The (2x1) phase has been observed and studied thoroughly experimentally. 9,18,22 Gray and black circles represent W and O atoms respectively. Dashed gray lines show the boundaries of the cell in the XY plane that we use in our calculations.

Experimental results 23 show that the hydrogen adsorption falls to zero for coverages Tungsten atoms in the corners of the depicted cell do not have any oxygen atoms in the nearby threefold sites. Therefore, our model for 0.25 ML coverage helps to understand the situation in which H 2 molecules collide with surface regions in which the O atoms are located at distances relatively far from the W atoms, but they are still affecting the dynamics. At the same time, the size of the cell is small enough to permit a sufficiently large number of AIMD calculations, which are computationally expensive. From this point, the notation Θ = 0.5 ML and Θ = 0.25 ML will be used to refer to the surfaces corresponding to Fig. 2.

Θ ≥ 0.
All our calculations have been performed using the Vienna Ab initio Simulation Package 6 2.3. PUBLICATIONS (VASP). 24 The electron-core interaction is treated with the plane augmented wave (PAW) approximation. 25,26 The exchange-correlation energy is calculated with the generalized gradient approximation (GGA) and the vdW xc-functional developed by Dion et al. (vdW-DF). 27 The implementation in VASP of this functional is performed through the routine written and provided by Jiri Klimeš 28 based on the Román-Pérez and Soler algorithm. 29 The construction of the slabs used to represent the different surfaces in Fig. 2 goes as follows. First, bulk calculations for W are performed obtaining the lattice constant a = 3.20

Å. Second, a 5-layer periodic slab of W atoms is relaxed. Each layer in the supercell contains four W atoms as represented with dashed lines in Fig. 2. Last, oxygen atom(s) are placed on the threefold site(s) at 1.5 Å over the surface and the system is relaxed again.

Using this method we obtain bond lengths between the O atoms and the nearby W atoms the sides of the rhombus showed in Fig. 2 and the last value representing the vertical distance between the first layer of one slab and the first layer of the next one. The vertical distance between two consecutive slabs is approximately 18 Å. These dimensions ensure that the selfinteraction of hydrogen is maintained to a minimum while calculation times are kept under reasonable bounds. The total force acting on each atom is kept under 10 meV/Å when the system reaches equilibrium. The energy cutoff used in all the calculations is 400 eV. The fractional occupancies are determined through the broadening approach of Methfessel and Paxton 30 with N = 1 and σ = 0.4. A mesh of 5 x 5 x 1 k-points within the Monkhorst-Pack method 31 is used in all cases.

In addition to these parameters, the integration step for all trajectories in AIMD calculations is 2fs and the energy convergence error for electronic relaxation is 10 -6 eV. Initially, the H 2 molecules are located above the surface at a distance of 9 Å with hydrogen atoms separated by the equilibrium distance of 0.741 Å. We consider a grid of initial kinetic energies in the range 50-300 meV, with a 50 meV step. Previous calculations in this range on a clean W(110) surface are available 12 and will be used for comparison purposes. For each energy, 400 trajectories have been simulated. The initial coordinates of the H 2 molecule in the XY plane, as well as the polar θ and azimuthal φ angles, are obtained using a conventional Monte-Carlo sampling over the whole cell limited by dashed gray lines in Fig. 2. The W and O atoms are initially fixed to their equilibrium positions but they are allowed to move during the dynamics. The latter means that some energy transfer between the incident molecule and the surface is included in our calculations. The zero-point energy of the incident molecule is ignored thus all calculations are purely classical. The criteria for the dissociation of the H 2 molecule is established as a separation in the internuclear distance r between H atoms of r > 1.5 Å and dr/dt > 0. When the H 2 molecule is moving away from the surface and the distance between molecule and surface is Z cm > 4.5 Å, a reflection event is considered.

As mentioned, previous calculations 12 of H 2 impinging on a clean W(110) surface using the same methodology and in the same energy range are available. Similar parameters have been used in both cases in order to compare the results (same supercell dimensions, xcfunctional, smearing method, etc). In particular, the same set of initial coordinates have been used for H 2 molecules over the clean surface and the surfaces with coverages Θ = 0.25

ML and Θ = 0.5 ML. This means that for each energy we generate 400 different sets of initial coordinates as described previously. Then, exactly the same sets are used in AIMD calculations for different coverages. This approach allows us to compare specific trajectories of interest at different coverages with the same initial conditions. to what happens when a clean W(110) surface is used. 12 Although not directly comparable, the absolute values of the sticking coefficient seem consistent with the experimental findings:

Results and discussion

the amount of H adsorbed on W(110) surfaces with 0.25 ML O coverage is roughly 1/4 of the amount of H adsorbed on clean W(110) surfaces. 23 The time evolution of the Z cm coordinate for all trajectories with an initial collision energy of 200 meV is displayed in Fig. 4. Different lines represent trajectories with different initial coordinates for a given coverage. The ending of the lines before 400 fs corresponds to a dissociation event. Several differences can be spotted for different coverages. First of all, reflection events appear at higher Z cm for the oxidized surfaces. The spread of these events also appears to be wider in Z cm and in time. there is no reflection event with turning point above 2.8 Å in the clean surface. However, for Θ = 0.5 ML there are trajectories that get reflected even at Z c m = 3.4 Å. This is a conspicuous effect of the presence of O repulsive centers above the W atoms.

Apart from that, in panels (a) and (b) the number of trajectories in which molecules are trapped near the surface at 400 fs is much higher in the latter. In a clean surface, only a small percentage of molecules were still trapped at 400 fs and none at 600 fs for all energies studied. After that time, all trapping ended up in dissociation. For Θ = 0.25 ML coverage, the percentage of trajectories trapped for a time longer than 400 fs near the surface increases greatly. These molecules are hovering over the whole cell and a big part of them end up dissociating and getting adsorbed in the surface. Nonetheless, around ∼20-25% of the trapped molecules go back to the vacuum. This behavior was not observed when a clean surface was used. 12 It seems that the adsorbed oxygen atoms open new paths in the phase space for H 2 to return to the vacuum. Finally, a small percentage of the total number of molecules are still trapped after 2 ps. The amount depends on the initial kinetic energy, Top panels in Fig. 5 show the distribution of the Z cm coordinate at the rebound point for all reflected trajectories with initial kinetic energy of 200meV (same as in Fig. 4). The bars under Z cm = 2.8 Å have been colored in black. This value corresponds to the maximum Z cm coordinate at the rebound point for molecules colliding with a clean surface. Reflection above this point only appears for the oxidized surfaces and the corresponding bars are colored in red. Bottom panels in Fig. 5 show the positions in the XY plane of the rebound points for the same trajectories as the top panels. It seems apparent that H 2 molecules being reflected at higher distances (in red) correspond to reflections near O atoms. Reflection at lower distances (black) appears to be due to the interaction with W atoms when there are no O atoms nearby. The amount of H 2 molecules getting reflected above 2.8 Å for Θ = 0.25 ML is roughly half of that for Θ = 0.5 ML. This fact strengthens the conclusion that these 11 O and W atoms is the main cause of this apparent inconsistency. On the other hand, the average distance to the surface at the rebound points change with the initial kinetic energy.

The higher the energy the more penetrative power the molecules will have resulting in lower averages of the vertical rebound point. For molecules getting reflected on W atoms (black group) these values are in the range ∼ 3-2.5 Å corresponding to energies in the range 50-300 meV respectively. Other than that, for different initial kinetic energies, the general behavior is similar to the one showed in Fig. 5.

The reflection process seems very similar for trajectories in the black group no matter the presence of O atoms or not. It might seem strange then that, for coverages of Θ = 0.5 ML, there are no trajectories leading to dissociation of H 2 in the surface. More than that, molecules do not even get closer to the surface according to Fig. 4 (a). Results for H 2 colliding with a clean W(110) surface 12 indicate that the H 2 path to approach the W atoms and undergo dissociation is narrow. Specifically, only molecules approaching on top of W atoms and with the orientation of the molecular axis roughly parallel to the surface will find favorable conditions. Figure 6 shows the position of molecules in the XY plane for all simulated trajectories with an initial kinetic energy of 200 meV. The first column of the panels shows the initial position (t = 0). The other three panels show the positions at times closer to the first interactions with the surface. The meaning of the gray and black circles, as well as the black and red dots, is the same as in the lower panels of Fig. 5. The new blue dots represent molecules that will undergo dissociation near the surface.

In the panels corresponding to the oxidized surfaces, we can observe the strong repulsive effect of O atoms on H 2 molecules. For Θ = 0.5 ML coverage, most molecules are pushed away from the top sites towards an imaginary line connecting the threefold sites that are not occupied by oxygen. This line would be the location further away from the adjacent rows of oxygen. Energetically, further approach to the surface would require the H 2 molecules to be even closer to the W top positions and not in the zone that this imaginary line crosses.

This helps to explain the lack of molecules approaching shorter distances to the surface. All previous analysis leads to the conclusion that adsorbed oxygen atoms effectively act as repulsive centers at the surface. This idea is confirmed by a parallel analysis in energy terms.

0 • 90 • 180 • Z cm = 9 Å 0 • 90 • 180 • Clean surface Z cm = 2.5 Å 0 • 90 • 180 • Z cm = 2.1 Å 0 • 90 • 180 • Z cm = 1.7 Å

Conclusions

In summary, we have used AIMD to study the interaction of H 2 molecules and a W(110)

surface with a Θ = 0.5 ML coverage of oxygen atoms. Calculations have been performed

as well in a second surface with an oxygen coverage of Θ = 0.25 ML. The model we use for the distribution of O atoms at the W surface helps to explain the differences in the dissociative adsorption probability between a clean W(110) surface 12 and a surface with an oxygen coverage of Θ = 0.5 ML. In the former, only H 2 molecules on top of W atoms and parallel to the surface will reach distances close enough to undergo dissociation. In the latter, the oxygen atoms push the incident H 2 molecules away from the top positions, effectively preventing dissociative adsorption at the surface. Calculations on the model surface with oxygen coverage of Θ = 0.25 ML show that the effective influence of the oxygen atoms seems to be restricted to only the closest W atoms. For the rest of W atoms, the dynamics are similar to those in a clean surface. The main difference is an increase in the time that it will take for molecules close to the surface to dissociate and the presence of molecules returning to vacuum from these close distances. Differences in the interaction energy between the H 2 molecule and the surface help to explain this behavior.

Our AIMD simulations show that the preadsorbed O atoms efficiently modulate the dissociation dynamics of the H 2 molecules moving in their vicinity. Quantitatively, our results for the sticking coefficient are consistent with available experimental information. ?

A direct comparison with experimental data for the case of 0M L ≤ Θ ≤ 0. at the surface and form well-ordered oxygen islands. Statistical kinetic simulations would be necessary on top of our calculations to describe the full process. We consider that the AIMD calculations presented here provide nevertheless a very useful ingredient to understand the complex dynamics of H 2 when dissociating over such a complex system, as the W(110)

surface partially covered with O atoms is.

2.4. SUMMARY two surfaces is observed.

The use of AIMD is an alternative to other theoretical approaches based on the construction of interpolated PESs. A weakness of AIMD is that it can be computationally expensive and therefore the number of trajectories that are launched can be significantly smaller than those launched using interpolated PESs, with the subsequent damage in statistical accuracy for low probability processes. It has other advantages, however, and one of them is that it can address more complex systems, often without much additional computational effort.

Conclusions

In the first chapter of this manuscript, we quantize, in Bohr's spirit, all the relevant final actions of a classically described reactive molecular system. The process at hand is the scattering of H 2 on Pd(111). In the beginning, the study focused on collisions where H 2 is initially in the rovibrational ground state. Then, rotationally excited states were considered. Two semiclassical corrections are used to assign statistical weights to simulated trajectories: Gaussian binning (GB) and the adiabaticity correction (AC).

The implementation of the quasi-classical trajectory method in its standard form fails to describe the sticking probability curve at low energies for this system. GB (action and energy based) combined with the AC lead to results in a much better agreement with quantum time-dependent wave packet (Q-TDWP) ones. In particular, the use of the energy-based Gaussian binning shows a remarkably accurate description of the rovibrational scattering probabilities near the thresholds.

A variation of the adiabaticity correction based on firmer semiclassical grounds was also introduced and applied to the study problem. The results do not vary much, reproducing accurately the quantum sticking and state-resolved reflection probabilities with the same computational cost. on the initial kinetic energy) are eventually dissociated. For the surface with the higher oxygen coverage, there is no dissociation independently of the initial collision energy. For the surface with the lower oxygen coverage, the sticking probability curve has qualitative similarities with the clean surface one, but the presence of preadsorbed oxygen provokes a considerable drop. In most cases, the dissociation is also decided at distances relatively far from the surface. There is, however, a small percentage of molecules returning to the vacuum after being trapped on the surface.

In all cases, H 2 needs to go through a narrow path in the phase-space before getting to distances close enough to the surface to get dissociated. Specifically, only molecules close to the top positions of W atoms and with an axis parallel to the surface will satisfy the required conditions. The pre-adsorption of oxygen on the surface does not change this behavior in zones where the W atoms are relatively separated from the O atoms. For an oxygen coverage of Θ = 0.5 ML, these zones do not exist due to the closely packed rows of oxygen atoms.

Resumen de la tesis

Las interacciones gas-superficie son clave para comprender muchos procesos de importancia industrial, atrayendo la atención de científicos teóricos y experimentales de todo el mundo [1]. Los avances en estos estudios han contribuido enormemente a la industria química a través de la comprensión de la catálisis heterogénea [2, 3]. Este proceso se usa ampliamente para acelerar reacciones a niveles industriales como la síntesis de amoníaco, ácidos nítrico y sulfúrico o la desulfurización del petróleo. Las interacciones gas-superficie son relevantes para la preparación de materiales semiconductores [4], componentes clave en los múltiples dispositivos electrónicos usados en el mundo contemporáneo. Una buena comprensión de estas interacciones tambien contribuye a la realización de estudios ambientales como los relacionados con la lluvia ácida [5] y el agujero de la capa de ozono. El estudio de las interacciones de los gases con superficies metálicas es clave en el trayecto hacia una futura economía basada en el hidrógeno [6], especialmente, en procesos como la purificación, el almacenamiento y la detección hidrógeno, así como la contrucción de celdas de combustible. No menos importantes son las contribuciones a la descripción de la dinámica de los gases en la superficie de los componentes en los futuros reactores de fusión.

En las últimas décadas, el interés por las interacciones gas-superficie ha crecido a pasos agigantados. Los avances en las técnicas de alto vacío han permitido obtener superficies más limpias que son utilizadas luego en múltiples experimentos. El refinamiento de los experimentos con haces moleculares ha mejorado enormemente la medición de forma sistemática de la cinética de reacciones en superficies [7]. Desde el punto de vista teórico, el desarrollo de cálculos de estructura electrónica ab initio basados en la teoría del funcional de la densidad ha mejorado la descripción teórica de procesos elementales en superficies metálicas. Además, se han logrado grandes avances en la generación de superficies de energía potencial [8, 9] y en el desarrollo de métodos de dinámica cuántica [10, 11]. La rápida evolución de los sistemas computacionales y la generación de códigos optimizados para implementar nuevas y antiguas metodologías, también ha contribuido en gran medida al crecimiento de este campo.

Incluso con todos los avances realizados en los últimos años, el evento más simple entre una molécula y una superficie sigue siendo un proceso complejo de estudiar. Es necesario comprender correctamente la dinámica de la molécula que se aproxima a la superficie y vuelve al vacío. Además, en la región cercana a la superficie, la interacción dependerá de muchos factores como el estado rovibracional de la molécula y la estructura electrónica, geometría, temperatura y naturaleza de la superficie, entre otros. La información sobre estos eventos "simples" es crucial para la comprensión y extrapolación a sistemas químicos más complejos.

Los cálculos de dinámica cuántica proporcionan una descripción muy precisa de los procesos de dinámica molecular. Esto los convierte en un método ideal para el estudio de interacciones en sistemas gas-superficie. Lamentablemente, todavía estamos lejos de poder emplear cálculos cuánticos completos para reacciones que involucren más de unos pocos (∼ 6) grados de libertad. En estos sistemas, los cálculos cuánticos tienen una complejidad cada vez mayor y los costos computacionales son prohibitivos. Por otro lado, los métodos dinámicos clásicos proporcionan una alternativa eficiente y, por lo general, ofrecen una comprensión más intuitiva del problema en cuestión. Por supuesto, estos tratamientos clásicos a menudo deben combinarse con correcciones semiclásicas cuando entran en juego efectos cuánticos en la dinámica.

Este manuscrito está dedicado al estudio teórico de varios procesos reactivos y no reactivos que tienen lugar en la interfaz entre gases y sólidos. Se utilizaron dos métodos clásicos de trayectorias, diferentes y complementarios, para simular la dinámica de estos procesos.

En un primer capítulo, usamos grandes conjuntos de caminos clásicos obtenidos mediante la resolución numérica de las ecuaciones de Hamilton en una superficie de energía potencial previamente construida. Luego, a estos caminos clásicos se les asignan pesos estadísticos basados en dos correcciones semiclásicas: el binning gaussiano y la corrección de adiabaticidad. Este enfoque en un espíritu cuántico se aplica a la dispersión de H 2 en una superficie de Pd(111).

En este capítulo, por primera vez, se ha realizado la cuantificación en el espíritu de las reglas de Bohr de todas las acciones finales relevantes en un sistema molecular reactivo descrito clásicamente. Tanto el binning gaussiano basado en las acciones (GB), como el binning gaussiano basado en las energías (1GB), se han empleado para analizar los resultados finales de las trayectorias teniendo en cuenta el principio de cuantificación de Bohr. Las acciones finales a cuantificar son las acciones vibracionales y rotacionales de la molécula de H 2 reflejada, así como las acciones difraccionales (análogos clásicos de los índices de Miller). Por otro lado, la corrección de adiabaticidad (AC) tiene en cuenta el hecho de que, en mecánica cuántica, algunas partes de la onda incidente pueden difractarse a energías de colisión muy pequeñas. Se han realizado dos tipos de simulaciones: cuando el H 2 está inicialmente en su estado fundamental interno y cuando se encuentra inicialmente en un estado rotacionalmente excitado.

En este sistema, el método de trayectorias cuasi-clásico no describe correctamente la curva de probabilidad de adsorción a bajas energías. La combinación de procedimientos de binning (GB y 1GB) con la AC conducen a resultados en una mejor concordancia con los resultados cuánticos de cálculos de paquetes de ondas dependientes del tiempo. En particular, el uso del 1GB muestra una descripción notablemente precisa de las zonas cercanas a los umbrales energéticos de la apertura de nuevos canales en las curvas de probabilidad de reflexión. En estas zonas, se demuestra la idoneidad del uso del 1GB en comparación con el GB.

Adicionalmente, en este mismo sistema de estudio, se ha utilizado una variación de la AC basada más firmemente en argumentos semiclásicos. Los resultados son muy similares a los obtenidos utilizando la implementación tradicional y las probabilidades de reflexión y adsorción obtenidas cuánticamente se reproducen con precisión. La aplicación de esta alternativa no incrementa los recursos computacionales a emplear.

En el futuro se planea repetir este tipo de estudios para diferentes reacciones en fase gaseosa y de gas-superficie para determinar si la calidad del actual acuerdo entre la metodología empleada en este trabajo y las probabilidades cuánticas exactas se mantiene con sistematicidad. Si se confirma que la precisión es igualmente buena para una gran variedad de sistemas y situaciones dinámicas, la metodología actual puede ser una alternativa fiable a los cálculos de dinámica cuántica que re-quieren una utilización excesiva de los recursos computacionales. Además, esta metodología puede ser aplicable a sistemas de complejidad creciente y/o regímenes de energía no accesibles actualmente en dinámica cuántica, como incidencias a bajas energías.

En un segundo capítulo hemos empleado el método conocido como dinámica molecular ab initio (AIMD) para simular la evolución temporal del sistema, a partir de las fuerzas internucleares obtenidas mediante la teoría del funcional de la densidad. Estas fuerzas se calculan sobre la marcha, para cada una de las posiciones de los átomos durante la dinámica, y se utilizan luego para mover clásicamente los núcleos. A diferencia del enfoque anterior, AIMD no requiere la construcción de una superficie de energía potencial. El precio a pagar, sin embargo, es que el costo numérico de cada trayectoria es mucho más alto que con el método anterior.

En este capítulo se han empleado cálculos completos de AIMD para estudiar los mecanismos de adsorción de H 2 en superficies de W(110). Se han estudiado tres superficies, una limpia y dos con cobertura de oxígeno igual a Θ = 0.25 ML y Θ = 0.5 ML. Hemos demostrado que, en la superficie limpia, el uso de un funcional que incluye correcciones de van der Waals (vdW-DF) mejora la comparación entre la descripción teórica y la información experimental disponible. Este es el funcional empleado en el resto de los cálculos realizados en este trabajo.

En una superficie limpia, la disociación de la molécula H 2 se decide a una distancia relativamente alejada de la superficie, de forma similar a lo que ocurre en la disociación de H 2 sobre W(100) o en la disociación de N 2 en superficies de W. Casi todas las moléculas de hidrógeno capaces de alcanzar una distancia determinada de la superficie (≈ 2 -2.5 Å, dependiendo de la energía cinética inicial) finalmente se disocian. Cuando se utiliza el funcional vdW-DF vemos que existe un breve atrapamiento dinámico en la superficie antes de que se produzca la disociación. Este es un mecanismo que, sin embargo, no influye en los valores finales de la probabilidad de disociación. Para la superficie de tungsteno con una cobertura de oxígeno Θ = 0.5 ML no existe disociación, independientemente de la energía de colisión inicial. Para una cobertura Θ = 0.25 ML, la curva de probabilidad de adsorción tiene similitudes cualitativas con la obtenida para una superficie limpia; pero la presencia de oxígeno provoca una disminución considerable de la probabilidad. En la mayoría de los casos, la disociación también se decide a distancias relativamente alejadas de la superficie. Sin embargo, existe un pequeño porcentaje de moléculas que regresan al vacío después de quedar atrapadas en la superficie.

En todos los casos estudiados, la molécula de H 2 necesita atravesar un camino estrecho en el espacio de fases antes de poder alcanzar distancias lo suficientemente cercanas a la superficie como para disociarse. Específicamente, solo las moléculas cercanas a las posiciones sobre los centros de los átomos de W y con una orientación del eje paralela a la superficie cumplen con las condiciones requeridas. La preadsorción de oxígeno en la superficie no cambia este comportamiento en zonas donde los átomos de W están relativamente separados de los átomos de O. Para una cobertura de oxígeno de Θ = 0.5 ML, estas zonas no existen debido a las filas compactas de átomos de oxígeno en la superficie. Estos átomos empujan las moléculas de H 2 incidentes lejos de las posiciones sobre los centros de los átomos de W evitando efectivamente la disociación cerca de la superficie. Para Θ = 0.25 ML, se observa una mezcla de lo que sucede en las otras dos superficies.

El uso de AIMD es una alternativa a otros enfoques teóricos basados en la construcción de superficies de energía potencial. Una debilidad de AIMD es que puede ser computacionalmente costoso y, por lo tanto, el número de trayectorias que se lanzan puede ser significativamente menor que las lanzadas utilizando superficies de energía potencial interpoladas, con el consiguiente perjuicio en la precisión estadística para procesos de baja probabilidad. Sin embargo, tiene otras ventajas. Una de ellas es que puede abordar sistemas más complejos, a menudo sin mucho esfuerzo computacional adicional.

Con la culminación de este trabajo hemos dado respuesta a varias interrogantes relacionadas con la interacción de hidrógeno molecular con superficies de paladio y tungsteno. Por supuesto, es necesaria la realización de otros estudios si queremos comprender completamente la dinámica en estos sistemas. En general, los resultados de este trabajo de tesis constituyen un avance importante en la comprensión de la dinámica en procesos químicos elementales en superficies metálicas.

TITRE: Dynamique classique des réactions gaz-surface: des fondements aux applications.

RÉSUMÉ: L'objet de cette thèse est l'étude théorique de processus réactifs et non réactifs se produisant à l'interface gaz-solide. Deux méthodes de trajectoires classiques, différentes et complémentaires, ont été utilisées afin de simuler la dynamique de ces processus. La première met en jeu un ensemble conséquent de trajectoires classiques obtenues en résolvant numériquement les équations de Hamilton sur une surface d'énergie potentielle (SEP) construite au préalable. Ces trajectoires sont pondérées par des poids statistiques conformément à deux contraintes semiclassiques: la pondération gaussienne et la correction d'adiabaticité. Cette approche, dans un esprit quantique, a été appliquée à la collision entre H 2 et la surface de Pd(111). Dans un premier temps, nous nous sommes limités au cas où H 2 se trouve dans son état rovibrationnel fondamental. Nous avons par la suite considéré ses états rotationnels excités. Il nous est alors apparu nécessaire de modifier la correction d'adiabaticité sur la base d'arguments semiclassiques rigoureux. Dans les deux cas, les prédictions des probabilités de collage et de reflections résolues en états se sont avérées être en accord remarquable avec celles obtenues par des calculs quantiques exacts, contrairement aux prédictions classiques standards. L'approche classique dans un esprit quantique pourrait ainsi s'avérer d'une grande utilité pour les études à venir. La seconde méthode utilisée dans ce travail, connu sous le vocable de Ab-Initio Molecular Dynamics (AIMD), permet de calculer les forces entre noyaux à partir de la théorie de la fonctionnelle densité et d'en déduire classiquement leurs déplacements. Contrairement à l'approche précédente, l'AIMD n'exige pas la construction généralement difficile d'une SEP (le prix à payer, toutefois, est que le coût numérique de chaque trajectoire est nettement plus élevé qu'avec la méthode précédente). L'AIMD nous a permis d'étudier le processus de dissociation de H 2 sur la surface de W(110). La fonctionnelle utilisée inclut un terme de van der Waals, qui provoque un acroissement de l'attraction à longue distance, compensé par une augmentation de la répulsion à courte distance. La combinaison des deux effets diminue de façon appréciable la probabilité de dissociation, alors en meilleur accord avec le résultat expérimental obtenu à l'aide d'une surface propre. Lorsque des atomes d'oxygène sont adsorbés au préalable sur la surface, la probabilité de dissociation chute considérablement. Cet effet est d'autant plus fort que la quantité d'oxygène adsorbé est forte. Un modèle de phase ordonnée a été utilisé afin d'expliquer l'absence de collage pour le taux de courverture Θ > 0.35 ML observé expérimentalement. Les atomes d'oxygène dévient les molécules H 2 des étroits passages conduisant au collage en l'absence des atomes d'oxygène. Ceci élimine toute chance de collage pour de forts taux de recouvrement. En revanche, pour de faible taux, on s'attend à ce qu'une dynamique similaire à celle sur la surface propre se produise au dessus des atomes de tungstène, et à une distance suffisamment grande des atomes d'oxygène.

TITLE: Classical dynamics of gas-surface scattering: fundamentals and applications.

ABSTRACT: This thesis manuscript is devoted to the theoretical study of several reactive and non-reactive processes that take place at the gas-solid interface. Two classical trajectory methods, different and complementary, were used to simulate the dynamics of these processes. The first one relies on large sets of classical paths obtained by numerically solving Hamilton equations on a previously constructed potential energy surface (PES). Classical paths are then assigned statistical weights based on two semiclassical corrections: Gaussian binning and the adiabaticity correction. This approach, in a quantum spirit, was applied to the scattering of H 2 on a Pd(111) surface. First, the study focused on collisions where H 2 is initially in the rovibrational ground state. Then, rotationally excited states were considered. On this occasion, a variation of the adiabaticity correction based on firmer semiclassical grounds was introduced. In both cases, the predictions of the sticking and state-resolved reflection probabilities were found to be in remarkably good agreement with those obtained through exact quantum timedependent calculations, contrary to standard quasi-classical trajectory predictions. The classical approach in a quantum spirit could thus be very useful for future studies.

The second method used in this work, known as Ab-Initio Molecular Dynamics (AIMD), calculates the inter-nuclear forces from density functional theory and uses them to classically move the nuclei. Contrary to the previous approach, AIMD does not require the very demanding construction of a PES (the price to pay, however, is that the numerical cost of each trajectory is much higher than with the previous method). AIMD allowed us to study the dissociation process of H 2 on W(110) surfaces. The functional we use includes a van der Waals term which provokes an increase of the far distance attraction that is compensated by a stronger repulsion at short distances. The combination of both effects appreciably decreases the value of the dissociation probability, bringing it closer to the experimental result when a clean surface is used. When oxygen atoms are previously adsorbed on the surface, the dissociation probability drops considerably. This effect increases with the amount of oxygen on the surface. An ordered phase of O adsorbates on the W surface is used to explain the nonexistent sticking probability for coverages Θ > 0.35 ML observed experimentally. We show that the oxygen atoms push the H 2 molecules away from the narrow bottlenecks that open the paths to dissociation in the absence of oxygen atoms. This effectively eliminates any chance of dissociation in the surface for high coverages. At lower coverages, our calculations demonstrate that the dissociation dynamics resemble those in the clean surface just in very specific surface regions.

Contents 1 . 1 6 1. 2 8 1. 2 . 1 9 1. 2 . 2 10 1. 2 . 3 12 1. 3 13 1. 3 . 1 40 1. 4

 116282192210231231331404 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Calculation details . . . . . . . . . . . . . . . . . . . . . . . . . . . . Initial conditions . . . . . . . . . . . . . . . . . . . . . . . . . Evolution in time . . . . . . . . . . . . . . . . . . . . . . . . . Statistical analysis . . . . . . . . . . . . . . . . . . . . . . . . Publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . When Classical Trajectories Get to Quantum Accuracy: The Scattering of H 2 on Pd(111) . . . . . . . . . . . . . . . . . . . 14 1.3.1.1 Supporting Information: Diffractional actions . . . 37 1.3.2 When Classical Trajectories Get to Quantum Accuracy: II. The Scattering of Rotationally Excited H 2 on Pd(111) . . . . Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

Fig. 1 .

 1 1 (a). Blue and black spheres represent Pd and H atoms respectively. X, Y, Z are the coordinates of the center of mass and r is the interatomic distance of the H 2 molecule. θ and φ are the polar and azimuthal angles respectively. The Palladium crystal structure presents a facecentered cubic (fcc) arrangement of atoms. In 1.1 (b) the unit cell of the Pd(111) face is depicted. The distances are expressed in function of the lattice constant a. ∆ is the distance between two nearest neighbour surface atoms.

Figure 1 . 1 :

 11 Figure 1.1: Coordinate system employed to describe the interaction of H 2 with the Pd(111) surface. Blue and black spheres represent Pd and H atoms respectively.

Figure 1 .

 1 2 shows this cell and a randomly generated set of N positions for the center of mass of H 2 in the plane. As in Fig. 1.1, blue circles represent Pd atoms and black dots the center of mass of the H 2 molecule. The values of N in panels (b) and (c) correspond to the number of trajectories per collision energy simulated in the calculations presented in Subsections 1.3.2 and 1.3.1 respectively. Panel (a) is only shown for comparison purposes. In all cases, the initial positions in the XY plane cover the whole cell evenly.

Figure 1 . 2 :

 12 Figure 1.2: Two-dimensional working cell. Blue circles represent Pd atoms and black dots randomly generated positions for the center of mass of H 2 in the XY plane. The values of N in panels (b) and (c) correspond to the number of trajectories per collision energy simulated in the calculations presented in Subsections 1.3.2 and 1.3.1 respectively. Panel (a) is only shown for comparison purposes.

CHAPTER 1 .

 1 THE SCATTERING OF H 2 ON PD(111)

CHAPTER 1 .Figure 1 . 3 :

 113 Figure 1.3: Graphical comparison between SB and GB procedures.

3 1Figure 1 :

 31 Figure 1: Quantum (black squares) and classical (blue circles) sticking probabilities for H 2 (0,0) + Pd(111) as a function of the collision energy.

8 CHAPTER 1 .

 81 Binning: 1GB. Instead of pseudoquantizing the four actions THE SCATTERING OF H 2 ON PD(111)

Figure 2 :

 2 Figure 2: Collision energy dependence of the sticking probability for H 2 (0,0) + Pd(111) according to quantum calculations (black squares) and various quantum-corrected classical calculations involving GB, 1GB and AC (green triangles and red circles).

Figure 3 :

 3 Figure 3: Collision energy dependence of rovibrationally elastic and inelastic scattering probabilities for H 2 (0,0) + Pd(111) according to quantum calculations (black squares), GB-AC-QCT calculations (green triangles) and 1GB-AC-QCT calculations (red circles); Final state (v, j) = (0, 0) (a), (0, 2) (b), and (0, 4) (c).

Figure 4 :

 4 Figure 4: Quantum (black squares) and 1GB-AC (red circles) reflection probabilities for (v, j)=(0,2) (panel(a)) and (0,4) (panel (b)). The contributions of the five first orders of diffraction to the 1GB-AC probability are shown (see text for the definition of an order).

15 1Figure 5 :

 155 Figure5: GB and 1GB weights are represented in panels (a) and (b), respectively. These weights have been rescaled in such a way that their maximum values are both equal to 1. Blue lines correspond to three different collision energies E a i , E b i and E c i whose values can be guessed from panel (c). E b i is the energy of state (0,2). The domains of the (x, J) space available at these energies are located on the left side of the blue lines in panels (a) and (b). Integrating the GB and 1GB weights assuming a uniform distribution in the available domain at E i leads to the green and red curves in panel (c), respectively. The sudden jump observed with 1GB shows the suitability of this procedure to realistically describe probabilities around energy thresholds.

  Acknowledgement A.R.F. acknowledges financial support by the University of Bordeaux. The authors acknowledge the support of France Grilles for providing computing resources on the French National Grid Infrastructure. Computer time was provided as well by the Pôle Modélisation HPC facilities of the Institut des Sciences Moléculaires UMR 5255 CNRS Université de Bordeaux, co-funded by the Nouvelle Aquitaine region as well as the MCIA (Mésocentre de Calcul Intensif Aquitain) resources of the Université de Bordeaux and of the Université de Pau et des Pays de l'Adour. The authors acknowledge financial support by the Spanish Ministerio de Economia, Industria y Competitividad Grant No. FIS2016-76471-P. This work was conducted in the scope of the transborder joint Laboratory "QuantumChemPhys: Theoretical Chemistry and Physics at the Quantum Scale (ANR-10-IDEX-03-02)." Supporting Information Available: This document contains detailed information about the diffractional actions.

( 45 )

 45 Bonnet, L.; Larrégaray, P.; Lara, M.; Launay, J.-M. Theoretical Study of Barrierless Chemical Reactions Involving Nearly Elastic Rebound: The Case of S( 1 D) + X 2 , X = H, D. J. Phys. Chem. A 2019, 123, 6439. (46) Bonnet, L.; Espinosa-García, J. The method of Gaussian weighted trajectories. V. On the 1GB procedure for polyatomic processes. J. Chem. Phys. 2010, 133, 164108. 23 1.3. PUBLICATIONS Supporting Information: Diffractional actions

Fig. S 1 :

 1 Fig. S 1: (a) Primitive basis vectors and (b) reciprocal lattice basis vectors for Pd(111).

3 1

 3

Figure 1 :

 1 Figure1: Rotational state distribution of reflected H 2 molecules for E i = 65 meV, j i = 0 (panel a) and j i = 1 (panel b). The blue and brown contributions are due to adiabatic and nonadiabatic paths, respectively. In ordinate, the trajectory number is for each bar of the histogram. In panel a, the peak in the blue distribution is due to paths being both vibrationally and rotationally adiabatic (see text for further details).

Figure 2 :

 2 Figure 2: Distribution of the polar angle θ at the instant of the rebound against the Pd(111) surface for type-FA trajectories.

H 2 8 CHAPTER 1 .

 281 molecule experiences an attractive force which accelerates its motion towards the Pd surface. The slope of the blue curve, thus, becomes slightly stronger from ∼ 120 fs up to ∼ 200 fs, where the rebound takes place. The rest of the blue curve is nearly symmetric THE SCATTERING OF H 2 ON PD(111)

Figure 3 :

 3 Figure 3: Time-dependence of the distance Z between the center-of-mass of H 2 and the surface for a typical type-FA trajectory (blue curve); same dependence for the potential energy (black curve).

Figure 4 :

 4 Figure 4: Realistic Z-dependence of the interaction potential energy for H 2 /Pd(111) (red curve); simple model of interaction potential involving the same repulsive wall as the red curve followed by a flat potential up to the vacuum (blue dashed curve).

  the vibrational and rotational motions are uncoupled along type-FA trajectories). Using these relations, one can exactly integrate the right-hand-side (RHS) of Eq. (

Figure 5 :

 5 Figure 5: Quantum (black squares) and standard QCT (blue circles) sticking probabilities for H 2 (0,1) + Pd(111) as a function of the collision energy.

Figure 6 :

 6 Figure 6: Panel (a): collision energy dependence of the sticking probability for H 2 (0,1)+ Pd(111) according to quantum (black squares), QCT-SC-I (red circles) and QCT-SC-II (green triangles) calculations. The standard QCT sticking probability is also plotted for comparison purposes (blue circles). Panels (b-d): same as before for the rovibrationally elastic and inelastic scattering probabilities. The final rotational states j f = 1, 3 and 5 are specified in each panel (the final vibrational quantum number is necessarily zero).

Fig. 6

 6 Fig.6show the collision energy dependence of the rotationally elastic and inelastic reflection probabilities. Panels (b), (c) and (d) correspond to final rotational states j f = 1, 3 and 5, respectively. The vertical dashed-dotted lines in panels (c) and (d) correspond to those of panel (a) and the QCT-SC-I and QCT-SC-II curves accurately describe the thresholds corresponding to the opening of the (0,3) and (0,5) states, in agreement with the quantum probabilities. We note that below ∼120 meV, the QCT-SC-II sticking probability (green 18

Figure 7 :

 7 Figure7: Collision energy dependence of the sticking probabilities for H 2 (0,j i = 2 -4) + Pd(111) according to quantum (black squares), QCT-SC-I (red circles), QCT-SC-II (green triangles) and standard QCT (blue circles) calculations.

Figure 8 :

 8 Figure 8: Collision energy dependence of rovibrationally elastic and inelastic scattering probabilities for H 2 (0, j i = 2 -4) + Pd(111) according to quantum (black symbols), QCT-SC-I (red symbols) and QCT-SC-II (green symbols) calculations. Final rotational states are indicated close to their corresponding curves.

Figure 9 :

 9 Figure 9: Collision energy dependence of the sticking probability for H 2 (0, 0) + Pd(111) according to quantum (black symbols), QCT-SC-I (red symbols) and QCT-SC-II (green symbols) calculations.

AcknowledgementA 22 CHAPTER 1 .

 221 .R.F. acknowledges financial support by the Université de Bordeaux. The authors acknowledge the support of France Grilles for providing computing resources on the French National Grid Infrastructure. Computer time was provided as well by the Pôle Modélisation HPC facilities of the Institut des Sciences Moléculaires, UMR 5255, CNRS, Université de Bordeaux, co-funded by the Nouvelle Aquitaine region as well as the MCIA (Mésocentre de Calcul Intensif Aquitain) ressources of the Université de Bordeaux and of the Université de Pau et des Pays de l'Adour. This work was conducted in the scope of the transborder joint Laboratory "QuantumChemPhys: Theoretical Chemistry and Physics at the Quantum Scale (ANR-10-IDEX-03-02)." RDM acknowledges financial support by the Gobierno Vasco-UPV/EHU Grant No. IT-1246-19. Finally, we wish to thank Dr. Cristina Diaz, of the Universidad Autónoma de Madrid, for providing us with the Q-TDWP results. THE SCATTERING OF H 2 ON PD(111)

  Gray and black spheres represent W and H atoms respectively. X, Y, Z are the coordinates of the center of mass and r is the inter-atomic distance of the H 2 molecule. θ and φ are the polar and azimuthal angles respectively. In panel (b) the unit cell of the W(110) face is depicted. The distances are expressed in function of the lattice constant a. We use the system of coordinates shown in Fig.2.1 (a) for generating the initial conditions in the dynamics and for checking if the H 2 molecule has reached one of the exit channels. For VASP calculations the cartesian system of coordinates is employed instead.

Figure 2 . 1 :

 21 Figure 2.1: (a) Coordinate system employed to describe the system H 2 /W(110). Gray and black spheres represent W and H atoms respectively. (b) Surface unit cell of W(110). Gray circles represent W atoms.

sections 2 .

 2 3.1 and 2.3.2, the characteristics of the supercells employed and the tests conducted are discussed.The crystal slab is constructed based on the specific plane of the surface to study, in our case, the (110) plane as shown in Fig.2.1 (b). The lattice constant obtained in the previous step is used to calculate the initial separation of the atoms. Once the slab is constructed a second ionic relaxation takes place to ensure that all atoms are in the equilibrium positions. When studying oxidized surfaces, a third ionic relaxation is enforced after the placement of oxygen atom(s) in the experimentally and theoretically predicted positions over the previously relaxed slab.

Figure 2 . 2 :

 22 Figure 2.2: Periodic repetition of a supercell containing a (2x2) surface unit. In each supercell the crystal slab is composed of three layers. Crystal slabs are separated by vacuum in the z direction. The supercell is indicated by the dashed blue lines.

Figure 2 . 3 :

 23 Figure 2.3: Initial sample distribution of 400 molecules in the XY plane for a clean W(110) surface. Grey circles show the positions of W atoms on the first layer of the slab. Each black dot represents the center of mass of an incident molecule. Grey dashed lines are the boundaries of the supercell in the XY plane.

CHAPTER 2 .

 2 AB-INITIO MOLECULAR DYNAMICS OF HYDROGEN ON TUNGSTEN SURFACESterm. Comparison of our results with previous theoretical work on the same system allows us to evaluate the relevance of the vdW dispersion forces in the dynamics, as well as the effect of making a full AIMD description of the problem.

  Å and a = 3.20 Å are used depending if the PBE xc-functional or the vdW-DF xc-functional are employed in the simulations.

Figure 1 : 6 CHAPTER 2 .

 162 Figure 1: (a) Coordinate system employed to describe the system H 2 /W(110). Gray and black spheres represent W and H atoms respectively. (b) Surface unit cell of W(110). Gray circles represent W atoms.

Figure 2

 2 Figure 2 shows the sticking probabilities of H 2 on W(110) as a function of the collision energy when using the PBE and vdW-DF xc-functionals. For comparison purposes, MD-C and MD-QC results from Ref. 7 are displayed as well. MD-C results using the PW-91 xcfunctional and AIMD results using the PBE xc-functional are very close quantitatively and qualitatively, despite the use of different functionals and methodology. On the other hand, once the vdW term is added, the sticking probability drops considerably (between 0.2-0.3) although the trend of the curve is still similar. The latter results are closer quantitatively to what was obtained in molecular beam experiments 2 as can be seen in Fig. 3. This suggests that the use of the vdW-DF xcfunctional put us one step closer to accurately describe dissociative adsorption in this system.

Figure 2 :

 2 Figure 2: Sticking probabilities for H 2 + W(110) as a function of the collision energy. Comparison between molecular dynamic classical and quasiclassical calculations (using a PES obtained with the PW91 xc-functional) from Ref. 7 and our classical AIMD calculations using the PBE and vdW-DF xc-functionals.

Fig 4 ,

 4 the time evolution of the distance between the H 2 molecule's center of mass and the surface is displayed. Each line represents one trajectory starting with random initial values for X cm , Y cm , θ and φ. The initial distance to the surface is fixed at Z cm = 9 Å and the separation of hydrogen atoms is set at the equilibrium distance r = 0.741 Å. Lines ending before 400 fs symbolize a dissociation event at that point. The initial collision energy is 200 meV in all cases. For other initial collision energies tested in our calculations, the general behavior is very similar. Differences arise, of course, in the number of reflected and dissociated molecules that varies according to Fig.2.Both panels of Fig.4exhibit a similar character for the whole process. After approaching the surface, either the molecule bounces back quickly or finds its way to get close to the W atoms. If the latter happens, the trajectories end up in a dissociation event in most of the cases. This effectively means that dissociative adsorption is determined far away from the

Figure 3 :

 3 Figure 3: Sticking probabilities for H 2 + W(110) as a function of the collision energy. Comparison between initial sticking probabilities obtained in molecular beam experiments as reported in Ref. 2 and classical ab initio molecular dynamics using the PBE and vdW-DF xc-functionals.

Fig. 6 10 CHAPTER 2 .Figure 4 :

 61024 Figure 4: Time evolution of the distance in the z-axis between the H 2 molecule's center of mass and the surface. Each line represents one trajectory starting with random initial values for X cm , Y cm , θ and φ. The initial distance to the surface is fixed at Z cm = 9 Å and the separation of hydrogen atoms is set at the equilibrium distance r = 0.741 Å. Lines ending before 400 fs symbolize a dissociation event at that point. Trajectories simulated using the PBE xc-functional are displayed in panel (a) while the vdW-DF xc-functional was employed for trajectories in panel (b).

Figure 5 :

 5 Figure 5: Probability to reach the distance Z cm = 2.5 Å and Z cm = 2 Å for H 2 molecules approaching the W(110) surface as a function of the initial kinetic energy. Black lines correspond to the use of the PBE xc-functional while red lines are obtained using the vdW-DF functional. The sticking curves from Fig. 2 and 3 are shown here as well to facilitate the comparison.

7 .

 7 Black and red dots (for calculations using PBE and vdW-DF respectively) represent the position of the center of mass of H 2 on trajectories leading to dissociation over the 2x2 cell used in this work (gray dashed lines). Grey circles show the position of W atoms on the first layer to help with the interpretation at different distances from the surface. In each of the sub-panels the position of the black and red dots has been obtained from the step of the calculations whose

Figure 6 :

 6 Figure 6: Interaction energy between the W(110) surface and the H 2 molecule for two trajectories leading to dissociative adsorption with initial collision energy of 200 meV. The energy is plotted as a function of the vertical distance Z cm between the molecule center of mass and the W surface. Both trajectories (one simulated by means of the PBE xc-functional and the other by means of the vdW-DF xc-functional) have the same set of initial coordinates.

Figure 7 :

 7 Figure 7: Evolution of the molecules that eventually dissociate. Different vertical distances Z cm are considered. Black and red dots represent the center of mass of molecules over the 2x2 cell employed (gray dashed lines). Grey circles correspond to the W atoms of the surface. Upper panel corresponds to trajectories simulated using the PBE xc-functional. For the lower panels, vdW-DF xc-functional has been used. The initial collision energy is 200 meV for all trajectories.

14 CHAPTER 2 .

 142 AB-INITIO MOLECULAR DYNAMICS OF HYDROGEN ON TUNGSTEN SURFACES

Figure 8 :

 8 Figure 8: Evolution of the molecules that eventually dissociate. Different vertical distances Z cm are considered. Distribution over polar angles θ for a binning of 10 • at different distances from the surface. Upper panel corresponds to trajectories simulated using the PBE xcfunctional. For the lower panels, vdW-DF xc-functional has been used. In both cases the total number of trajectories simulated is 400. The initial collision energy is 200 meV for all trajectories.

  Acknowledgement A.R.F. acknowledges financial support by the University of Bordeaux. This work was conducted in the scope of the transborder joint Laboratory "QuantumChemPhys: Theoretical Chemistry and Physics at the Quantum Scale" (ANR-10-IDEX-03-02). This work has been supported in part by the Basque Departamento de Educación, Universidades e Investigación, the University of the Basque Country UPV/EHU (Grant No. IT1246-19) and the Spanish Ministerio de Ciencia e Innovación (PID2019-107396GB-I00/AEI /10.13039/501100011033).

4 2. 3 .

 43 PUBLICATIONS numerical implementation are described in Sec. 2. Our results are shown in Sec. 3. The main conclusions are summarized in Sec. 4.

5 CHAPTER 2 .Figure 1 :

 521 Figure 1: (a) Coordinate system employed to describe the system. Gray and black spheres represent W and H atoms respectively. (b) Surface unit cell of W(110) with an O atom in the stable chemisorption position. Gray circles represent W atoms. The small black circle represents an oxygen atom in the threefold hollow site following results from Ref. 9,18,22.

Figure 2 :

 2 Figure 2: Oxidized surfaces with different coverages. (a) Ordered phase (2x1) corresponding to Θ = 0.5 ML (b) Theoretical model corresponding to Θ = 0.25 ML. Gray and black circles represent W and O atoms respectively. Dashed gray lines show the boundaries of the cell employed in the calculations.

of 2 . 7 CHAPTER 2 .

 272 16 Å, close to the experimental value of 2.08 Å found in Ref. 18. This corresponds to a perpendicular distance of 1.21 Å measured from the topmost layer of W atoms in the surface. The supercell dimensions are 5.55 x 5.55 x 27.22 Å, with the first two values being AB-INITIO MOLECULAR DYNAMICS OF HYDROGEN ON TUNGSTEN SURFACES

Figure 3

 3 Figure3shows the sticking probabilities of H 2 on W(110) at various oxygen coverages as a function of the collision energy. Experimental results23 show that for coverages Θ > 0.35 ML hydrogen adsorption gets inhibited by the presence of oxygen. Our simulation results are in agreement with this conclusion: at a coverage Θ = 0.5 ML no adsorption is observed in the whole energy range studied. On the other hand, at Θ = 0.25 ML coverage, the sticking coefficient is non zero and it increases with the collision energy in a similar fashion

For example, for this initial collision energy 9 CHAPTER 2 .

 92 AB-INITIO MOLECULAR DYNAMICS OF HYDROGEN ON TUNGSTEN SURFACES

Figure 3 :

 3 Figure 3: Sticking probabilities for the studied system as a function of the collision energy. The probabilities drop as the coverage increases for all the collision energies. Similar behavior has been observed experimentally in Ref. 23.

10 2. 3 .

 103 PUBLICATIONSgoing from 4 % at 50 meV to 2 % at 150 meV and dissapearing at 200 meV. For Θ = 0.5 ML coverage, molecules can not reach distances near the surface where dissociation occurs for lower coverages. This effectively prevents these events from appearing. The main features of this figure remain unchanged for different collision energies.

Figure 4 :

 4 Figure 4: Time evolution of the distance in the z-axis between the H 2 molecule's center of mass and the surface. Different lines represent trajectories with different initial coordinates for a given coverage. The ending of the lines before 400 fs represents a dissociation event. Initial collision energy is 200 meV in all cases.

CHAPTER 2 .

 2 AB-INITIO MOLECULAR DYNAMICS OF HYDROGEN ON TUNGSTEN SURFACES rebounds occur due to O atoms.

Figure 5 :

 5 Figure 5: Top panels: distribution of the Z cm coordinate at the rebound point for all reflected trajectories with collision energy 200meV. Bottom panels: positions of the rebound points for the same trajectories as top panels in the XY plane. Black and red bars (top) and dots (bottom) correspond to molecules that get reflected on W and O atoms respectively. Grey and black circles stand for W and O atoms in the surface as in Fig. 2. Dashed gray lines show the edges of the supercell in the plane XY.

13

 13

CHAPTER 2 .

 2 AB-INITIO MOLECULAR DYNAMICS OF HYDROGEN ON TUNGSTEN SURFACES

Figure 6 :

 6 Figure 6: Evolution of the molecules in the XY plane. Different times are considered. Dots represent the center of mass of molecules over the 2x2 cell employed (gray dashed lines).Black and red dots correspond to molecules reflected at different distances from the surface as in Fig.5. Blue dots represent molecules that will undergo dissociation near the surface. Grey and black circles symbolize the W and O atoms of the surface. The initial collision energy is 200 meV for all trajectories. Each row of panels corresponds to different oxygen coverages on the W(110) surface.

Figure 7

 7 Figure7shows the distribution over polar angles θ for trajectories that end in dissociation

Z

  cm = 2.5 Å Z cm = 2.1 Å Z cm = 1.7 Å

Figure 7 :

 7 Figure 7: Evolution of the molecules that eventually dissociate. Different vertical distances Z cm are considered. Distribution over polar angles θ for a binning of 10 • at different distances from the surface. Top panels correspond to results on a clean surface and bottom panels to an oxygen coverage of Θ = 0.25 ML. The initial collision energy is 200 meV for all trajectories.

Figure 8 Figure 8 :

 88 Figure 8 shows the interaction energy between H 2 and the studied surfaces for molecules with initial kinetic energy equal to 200meV. Each panel shows three lines corresponding to a molecule approaching to a clean surface and surfaces with O coverages of Θ = 0.25 ML and Θ = 0.5 ML. The initial coordinates of the molecules represented in every panel are the same. We have checked that the trajectories in Fig. 8 are representative of the general behavior for all kinetic energies considered.

  In the second chapter, full ab initio molecular dynamics calculations are used to study the adsorption mechanisms of H 2 on W(110) surfaces. At first, the dynamics of H 2 on a clean surface were studied. The changes in the dissociation mechanisms of H 2 in presence of preadsorbed oxygen on the surface were investigated next. On the clean surface, the sticking probability curve calculated using the vdW-DF xcfunctional is compared with theoretical and experimental results. The inclusion of the van der Waals term improves the comparison of the theoretical description with available experimental information. This is the xc-functional used for the rest of this work.In a clean surface, the dissociation of the H 2 molecule is decided at a distance relatively far from the surface, in a way similar to what happens in the dissociation of H 2 on W(100) or in the dissociation of N 2 on W surfaces. Almost all hydrogen molecules able to reach a given distance from the surface (≈ 2 -2.5 Å, depending

	Testing in other systems is still necessary to figure out whether the quality of the
	present agreement between quantum corrected QCT and exact quantum probabili-
	ties are systematically maintained or not. If it is indeed confirmed that the accuracy
	is equally good for a large variety of systems and dynamical situations, the current
	methodology can be a reliable alternative to computationally heavy quantum dy-
	namics calculations. It may be also applicable to systems of increasing complexity
	and/or energy regimes not currently accessible in quantum dynamics, such as low
	incident energies.

  [8, 9] et le développement de méthodes de dynamique quantique[10, 11].L'évolution rapide des systèmes informatiques et la génération de codes optimisés pour mettre en oeuvre des méthodologies nouvelles et anciennes, a également largement contribué à l'essor de ce domaine.Malgré tous les progrès réalisés ces dernières années, la rencontre la plus simple entre une molécule et une surface reste un processus complexe à étudier. Il est nécessaire de bien comprendre la dynamique de la molécule lorsqu'elle s'approche de la surface puis retourne dans le vide. De plus, dans la région proche de la surface, l'interaction dépendra de nombreux facteurs tels que l'état rotationnel de la molécule et la structure électronique, la géométrie, la température et la nature de la surface, entre autres. Les informations sur ces événements «simples» sont cruciales pour pouvoir appréhender des systèmes chimiques plus complexes.Dans ce chapitre, la quantification dans l'esprit des règles de Bohr, a été réalisée pour la première fois pour toutes les actions finales pertinentes. Aussi bien le binning gaussien basé sur des actions (GB), que le binning gaussien basé sur les énergies (1GB), ont été utilisés pour analyser les résultats finaux des trajectoires en tenant compte du principe de quantification de Bohr. Les actions finales à quanti-

	sequent damage in statistical accuracy for low probability processes. It has other non actuellement accessible en dynamique quantique, comme les incidences à faible
	advantages, however, and one of them is that it can address more complex systems, énergies.
	often without much additional computational effort. Dans un deuxième chapitre, nous avons utilisé la méthode connue sous le voca-
	ble dynamique moléculaire ab initio (AIMD) pour simuler l'évolution temporelle du
	système, à partir des forces inter-nucléaires obtenues par la théorie de la fonction-Résumé de la thèse fier sont les actions de vibration et de rotation de la molécule H 2 réfléchie, ainsi que nelle de la densité. Ces forces sont calculées à la volée, pour chacune des positions
	les actions de diffraction (analogues classiques des indices de Miller). D'autre part, d'atomes pendant la dynamique, et sont ensuite utilisées pour déplacer classique-
	Les interactions gaz-surface sont essentielles pour comprendre de nombreux proces-la correction d'adiabaticité (AC) prend en compte le fait qu'en mécanique quan-ment les atomes. Contrairement à l'approche précédente, AIMD ne nécessite pas la
	sus d'intérêts industriels, attirant l'attention des scientifiques théoriciens et expéri-tique, certaines parties de l'onde incidente peuvent être diffractées à des énergies construction d'une surface d'énergie potentielle. Le prix à payer, cependant, est que
	mentateurs du monde entier [1]. Les progrès de ces études ont grandement con-de collision très basses. Deux types de simulations ont été réalisés: lorsque H 2 est le coût de chaque chemin est beaucoup plus élevé qu'avec la méthode précédente.
	tribué à l'industrie chimique à travers la compréhension de la catalyse hétérogène initialement dans son état fondamental interne et lorsqu'il est initialement dans un
	état rotationnel excité.
	Sur une surface propre, la dissociation de la molécule H 2 est décidée à une dis-
	tance relativement éloignée de la surface, similairement à ce qui se passe dans la
	dissociation de H 2 sur W(100) ou dans la dissociation de N 2 sur des surfaces W.
	Presque toutes les molécules d'hydrogène capables d'atteindre une distance don-
	née de la surface (≈ 2 -2.5 Å, en fonction de l'énergie cinétique initiale) sont finale-
	Ce manuscrit est consacré à l'étude théorique de divers processus réactifs et non-ment dissociées. En utilisant la fonctionnelle vdW-DF, nous voyons qu'il y a un bref
	Au cours des dernières décennies, l'intérêt pour les interactions gaz-surface réactifs qui ont lieu à l'interface entre gaz et solides. Deux méthodes de trajectoires L'application de cette alternative n'augmente pas les ressources de calculs utilisées. piégeage dynamique par la surface avant que la dissociation ne se produise. C'est
	These atoms push the in-s'est accru à pas de géant. Les progrès des techniques de vide poussé ont permis classiques, à la fois différentes et complémentaires, ont été utilisées pour simuler la À l'avenir, il est prévu de répéter ce type d'études pour différentes réactions un mécanisme qui, cependant, n'influence pas les valeurs finales de la probabilité de
	cident H 2 molecules away from the top positions effectively preventing dissociation d'obtenir des surfaces plus propres qui sont ensuite utilisées dans de multiples ex-dynamique de ces processus. en phase gazeuse et de type gaz-surface pour déterminer si la qualité de l'accord dissociation. Pour une surface en tungstène avec un revêtement d'oxygène Θ = 0.5
	near the surface. At Θ = 0.25 ML, a mix of what happens in the other two surfaces périences. Le raffinement des expériences de faisceaux moléculaires a grandement Dans un premier chapitre, nous utilisons de grands ensembles de chemins clas-actuel entre la méthodologie utilisée dans ce travail et les probabilités quantiques ML, il n'y a pas de dissociation, quelle que soit l'énergie de collision initiale. Pour
	is observed. amélioré la mesure systématique de la cinétique de réaction de surface [7]. D'un siques obtenus en résolvant numériquement les équations de Hamilton sur une exactes est maintenue de façon systématique. Si la précision est confirmée pour une une couverture Θ = 0.25 ML, la courbe de probabilité d'adsorption présente des
	The use of AIMD is an alternative to other theoretical approaches based on the point de vue théorique, le développement de calculs de structure électronique ab surface d'énergie potentielle précédemment construite. On attribue alors à ces grande variété de systèmes et de situations dynamiques, notre méthodologie peut similitudes qualitatives avec celle obtenue pour une surface propre; mais la présence
	construction of interpolated PESs. A weakness of AIMD is that it can be compu-chemins des poids statistiques basés sur deux corrections semi-classiques: le bin-devenir une alternative fiable aux calculs de dynamique quantique qui nécessitent d'oxygène entraîne une diminution considérable de la probabilité. Dans la plupart
	tationally expensive and therefore the number of trajectories that are launched can une utilisation excessive des ressources de calcul. De plus, notre méthodologie peut
	be significantly smaller than those launched using interpolated PESs, with the sub-être applicable à des systèmes de plus en plus complexes et à des régimes d'énergie
	123

[2, 3]

. Ces procédés, tels que la synthèse de l'ammoniac, des acides nitrique et sulfurique ou la désulfuration de l'huile, sont largement utilisés pour accélérer les réactions aux niveau industriel. Les interactions gaz-surface sont importantes pour la préparation des matériaux semi-conducteurs

[4]

, composants clés de plusieurs dispositifs électroniques utilisés dans le monde contemporain. Une bonne compréhension de ces interactions contribue également à la réalisation d'études environnementales telles que celles portant sur les pluies acides

[5]

et le trou dans la couche d'ozone. L'Etude des interactions des gaz avec les surfaces métalliques sont essentielles pour arriver un jour à une économie basée sur l'hydrogène

[6]

, en particulier via des processus tels que la purification, le stockage et la détection de l'hydrogène, ainsi que la construction de réservoirs de carburant. Les contributions à la description de la dynamique des gaz à la surface des composants dans les futurs réacteurs à fusion sont également d'une grande importance. initio basés sur la théorie de la fonctionnelle de la densité a amélioré la description théorique des processus élémentaires sur les surfaces métalliques. En outre, de grands progrès ont été réalisés dans la génération de surfaces d'énergie poten-tielle Les calculs quantiques permettent d'arriver à une compréhension précise de la structure électronique et de la dynamique moléculaire. Cela en fait une méthode idéale pour l'étude des interactions dans les systèmes gaz-surface. Malheureusement, on est encore loin de pouvoir réaliser des calculs quantiques exacts pour étudier des réactions impliquant plus de ~6 degrés de liberté. Pour ces systèmes problématiques, les calculs quantiques sont de plus en plus complexes et leurs coûts sont prohibitifs. D'autre part, les méthodes de dynamique classique sont une alternative efficace et offrent généralement une compréhension plus intuitive du problème étudié. Bien sûr, ces traitements classiques doivent souvent être associées à des corrections semi-classiques lorsqu'entrent en jeu des effets quantiques dans la dynamique. ning gaussien et la correction d'adiabaticité. Cette approche dans un esprit quantique est appliquée à la collision de H 2 avec une surface de Pd(111). Dans ce système, la méthode de trajectoire quasi-classique ne décrit pas correctement la courbe de probabilité d'adsorption aux basses énergies. Combiner les procédures de binning (GB et 1GB) avec l'AC conduisent à des résultats en meilleur accord avec les résultats quantiques des calculs de paquets d'ondes dépendant du temps. En particulier, l'utilisation de 1GB montre une description remarquablement précise des zones proches des seuils énergétiques d'ouverture de nouveaux canaux dans les courbes de probabilité de réflexion. Dans ces domaines, cela démontre la pertinence d'utiliser 1GB par rapport à GB. De plus, dans ce même système d'étude, une variante de l'AC a été utilisée, s'appuyant plus fermement sur des arguments semi-classiques. Les résultats sont très similaires à ceux obtenus en utilisant l'AC précédente et les probabilités de réflexion quantique et d'adsorption obtenues sont reproduites avec précision. Dans ce chapitre, des calculs AIMD complets ont été utilisés pour étudier les mécanismes d'adsorption de H 2 sur les surfaces de W (110). Trois surfaces ont été étudiées, une propre et deux avec une couverture d'oxygène égale à Θ = 0.25 ML et Θ = 0.5 ML. Nous avons montré que, sur la surface propre, l'utilisation d'une fonctionnelle incluant des corrections de van der Waals (vdW-DF) améliore la comparaison entre la description théorique et les informations expérimentales disponibles. C'est cette fonctionnelle qui est utilisée dans le reste des calculs effectués dans ce travail.

des cas, la dissociation est également décidée à des distances relativement éloignées de la surface. Cependant, il y a un petit pourcentage de molécules qui retournent

Dissociation distance, usually taken as a value between two and three times the equilibrium interatomic distance in the vaccum.

1.3. PUBLICATIONS

CHAPTER 1. THE SCATTERING OF H 2 ON PD(111)

2.2. CALCULATION DETAILS

2.3. PUBLICATIONS

les effets de dissociation active près de la surface. Pour Θ = 0.25 ML, on observe un mélange de ce qui se passe sur les deux surfaces précédentes. L'utilisation de l'AIMD est une alternative à d'autres approches théoriques basées sur la construction de surfaces d'énergie potentielle. L'une des faiblesses de l'AIMD est qu'elle peut être coûteuse en temps de calcul et donc le nombre de trajectoires lancées peut être nettement inférieur à celui de trajectoires calculées sur des surfaces d'énergie potentielle interpolées, avec des dommages conséquents sur la précision statistique pour les processus à faible probabilité. Cependant, l'AIMD présente d'autres avantages. Le principal d'entre eux est que l'AIMD permet de s'attaquer à des systèmes plus complexes, souvent sans trop d'efforts calculatoires supplémentaires.Grâce à cette approche, nous avons répondu à plusieurs questions liées à l'interaction de l'hydrogène moléculaire avec les surfaces de palladium et de tungstène. Bien sûr, d'autres études sont nécessaires si l'on veut comprendre pleinement la dynamique de ces systèmes. De façon générale, les résultats de ce travail de thèse constituent une avancée importante dans la compréhension de la dynamique des processus chimiques élémentaires sur les surfaces métalliques.

Acknowledgement A.R.F. acknowledges financial support by the University of Bordeaux. This work was conducted in the scope of the transborder joint Laboratory "QuantumChemPhys: Theoretical Chemistry and Physics at the Quantum Scale" (ANR-10-IDEX-03-02). This work has been supported in part by the Basque Departamento de Educación, Universidades e Investigación, the University of the Basque Country UPV/EHU (Grant No. IT1246-19) and the Spanish Ministerio de Ciencia e Innovación (PID2019-107396GB-I00/AEI /10.13039/501100011033). dans le vide après avoir été piégées en surface. Dans tous les cas étudiés, la molécule H 2 doit emprunter un chemin étroit dans l'espace des phases avant de suffisamment se rapprocher de la surface pour pouvoir se dissocier. Plus précisément, seules les molécules proches des positions sur les centres des atomes W et avec une orientation de l'axe parallèle à la surface répondent aux conditions requises. La preadsorption d'oxygène en surface ne modifie pas ce comportement dans les zones où les atomes W sont relativement séparés des atomes O. Pour une couverture d'oxygène de Θ = 0.5 ML, ces zones n'existent pas en raison des lignes compactes d'atomes d'oxygène sur la surface. Ces atomes poussent les molécules H 2 incidentes éloignés des positions sur les centres des atomes W, évitant

When Classical Trajectories Get to Quantum Accuracy: II. The Scattering of Rotationally Excited H 2 on Pd(111)

A. Rodríguez-Fernández, * , †, ‡ L. Bonnet, * , †, ¶ C. Crespos, * , †, ¶ P. Larrégaray, * , †, ¶ and R. Díez Muiño * , §, ‡ †Université de Bordeaux, ISM, UMR 5255, F-33400 Talence, France. ‡Centro de Física de Materiales CFM/MPC (CSIC-UPV/EHU), Paseo Manuel de Lardizabal 5, 20018 Donostia-San Sebastián, Spain. ¶CNRS, ISM, UMR 5255, F-33400 Talence, France. §Donostia International Physics Center (DIPC), Paseo Manuel de Lardizabal 4, 20018 Donostia-San Sebastián, Spain. E-mail: alberto.rodriguez-fernandez@u-bordeaux.fr; claude-laurent.bonnet@u-bordeaux.fr; cedric.crespos@u-bordeaux.fr; pascal.larregaray@u-bordeaux.fr; rdm@ehu.es

Abstract

The classical trajectory method in a quantum spirit assigns statistical weights to classical paths on the basis of two semiclassical corrections: Gaussian binning and the adiabaticity correction. This approach was recently applied to the heterogeneous gassurface reaction between H 2 in its internal ground state and Pd(111) surface e.g. E-mail: alberto.rodriguez-fernandez@u-bordeaux.fr; claude-laurent.bonnet@u-bordeaux.fr; pascal.larregaray@u-bordeaux.fr; rdm@ehu.es

Abstract

The dissociation process of hydrogen molecules on W(110) was studied using density functional theory and classical molecular dynamics. We have calculated the dissociation probability for molecules with energies below 300 meV and analyzed the dynamics of the adsorption process. Our results show that the fate of each trajectory is determined at distances relatively far from the surface, at roughly 2-2.5 Å. This distance varies slightly with the initial kinetic energy of the molecule. Part of our simulations include van der Waals dispersion effects in the interaction between molecule and surface. We E-mail: alberto.rodriguez-fernandez@u-bordeaux.fr; claude-laurent.bonnet@u-bordeaux.fr; pascal.larregaray@u-bordeaux.fr; rdm@ehu.es

Abstract

Hydrogen molecules dissociate on clean W(110) surfaces. The H 2 dissociation is progressively inhibited as the W surface is pre-covered with oxygen. We use density functional theory and ab-initio molecular dynamics to analyze and understand the effect of the adsorbed O atoms on the H 2 dissociation process. The dissociation probability for H 2 molecules is calculated for energies below 300 meV and different O nominal coverages. We show that the adsorbed O atoms are repulsive centers that modulate the dynamics of the approaching H 2 molecules. The H 2 dissociation path becomes harder to follow when O atoms are nearby. In agreement with existing experimental information, we find that H 2 dissociation is absent for an O coverage of half a monolayer. Our results 1

Summary

In this chapter full ab initio molecular dynamics (AIMD) calculations have been employed to study the adsorption mechanisms of H 2 on W(110) surfaces. Three surfaces have been studied, a clean one and two with oxygen coverage equal to Θ = 0.25 ML and Θ = 0.5 ML. We have shown that in the clean surface the use of the vdW-DF xc-functional improves the comparison of the theoretical description with available experimental information. This is the xc-functional employed in the rest of the calculations in this work.

In a clean surface, the dissociation of the H 2 molecule is decided at a distance relatively far from the surface, in a way similar to what happens in the dissociation of H 2 on W(100) or in the dissociation of N 2 on W surfaces. Almost all hydrogen molecules able to reach a given distance from the surface (≈ 2 -2.5 Å, depending on the initial kinetic energy) are eventually dissociated. When using the vdW-DF xc-functional, we obtain short-time dynamic trapping at the surface before the dissociation takes place. This is a mechanism that nevertheless does not influence the final values of the dissociation probability.

For the higher O coverage surface studied here (Θ = 0.5 ML) there is no dissociation independently of the initial collision energy. For the lower coverage one (Θ = 0.25 ML), the sticking probability curve has qualitatively similarities to the results in the clean surface but the presence of oxygen provokes a considerable drop in the sticking values. In most cases, the dissociation is also decided at distances relatively far from the surface. There is, however, a small percentage of molecules returning to the vacuum after being trapped on the surface. In all cases, H 2 needs to go through a narrow path in the phase-space before reaching distances close enough to the surface to get dissociated. Specifically, only molecules close to the top positions of W atoms and with an axis parallel to the surface will satisfy the required conditions. The pre-adsorption of oxygen on the surface does not change this behavior in zones where the W atoms are relatively separated from the O atoms. For an oxygen coverage of Θ = 0.5 ML, these zones do not exist due to the closely packed rows of oxygen atoms. These O atoms push the incident H 2 molecules away from the top W positions effectively preventing dissociation near the surface. At Θ = 0.25 ML, a mix of what happens in the other