# Adsorption dynamics of molecular nitrogen at an Fe(111) surface

M.A. Nosir, L. Martin-Gondre, G.A. Bocan, R. Díez Muiño

► To cite this version:

M.A. Nosir, L. Martin-Gondre, G.A. Bocan, R. Díez Muiño. Adsorption dynamics of molecular nitrogen at an Fe(111) surface. Physical Chemistry Chemical Physics, 2017, 19, pp.7370-7379. 10.1039/c6cp07174k. hal-01514077

HAL Id: hal-01514077

https://hal.science/hal-01514077v1

Submitted on 2 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Adsorption Dynamics of Molecular Nitrogen
at a Fe(111) Surface

M. A. Nosir,*,†,‡ L. Martin-Gondre,¶ G. A. Bocan,§ and R. Díez Muiño†,‡

†Centro de Física de Materiales CFM/MPC (CSIC-UPV/EHU), Donostia-San Sebastián,
Spain

‡Donostia International Physics Center (DIPC), Donostia-San Sebastián, Spain

¶ Institut UTINAM - Université Bourgogne Franche-Comté, Besançon, France

§CONICET and Centro Atómico Bariloche (CNEA), Argentina

E-mail: mohamed_ahmed@ehu.es

Abstract

We present an extensive theoretical study of N₂ adsorption mechanisms on a Fe(111) surface. We combine the static analysis of a six-dimensional potential energy surface (6D-PES), based on ab-initio density functional theory (DFT) calculations for the system, with quasi-classical trajectory (QCT) calculations to simulate the adsorption dynamics. There are four molecular adsorption states, usually called γ, δ, α, and ε, arising from our DFT calculations. We find that N₂ adsorption into the γ-state is non-activated, while there is a threshold energy associated with the entrance channel for the other three adsorption states. Our QCT calculations confirm that there are activated and nonactivated paths for the approach of N₂ to the Fe(111) surface, in agreement with previous experimental investigations. Molecular dynamics for a surface temperature Tₛ = 300 K and impact energies Eᵢ in the 0 - 5 eV range show the relative occupancy of the γ, δ, α, and ε states. The δ-state, however, is only marginally populated despite its adsorption energy being very similar to that of the γ-state. Our QCT calculations
trace the dependence of the molecular trapping on the surface temperature $T_s$ and ini-
tial impact energy $E_i$ and quantify the rates of the different competitive channels that
eventually lead to molecular adsorption.

Introduction

Since Haber and Bosch designed the first industrial reactor of iron-catalyzed ammonia synthesis from $H_2$ and $N_2$, great efforts have been made to improve the design of the reactor and understand the nature of the catalytic process as well as the associated reactions. Concerning the latter, the difficulty is that elementary reactions at surfaces are complex phenomena, which can be highly dependent on minor details of the electronic interaction.¹,² A generally accepted conclusion extracted from previous studies is that the dissociative adsorption of $N_2$ is the rate-determining step in the catalytic formation of ammonia at surfaces.³⁻⁶ A large body of research has thus been devoted to the study of fundamental processes of nitrogen atoms and molecules at metal surfaces.

Processes that take place at the interface between the gas and solid phases often show strong sensitivity to the crystal surface structure.⁷ In particular, metal surfaces consisting of identical atomic species can give rise to completely different energy landscapes when switching between different arrangements of surface atoms. Therefore, the chemical reactivity can be enhanced or hindered depending on the particular crystal face under study.

A standard experimental technique for studying the interaction between gases and metal surfaces are molecular beam experiments. Existing molecular beam measurements show dramatic changes, for instance, in the initial sticking probability or sticking coefficient $S$ for $N_2$ molecules at two different tungsten faces, W(100) and W(110). At relatively low surface temperature and thermal energies of the beam, the sticking coefficient is roughly $S \approx 1$ on W(100).⁸⁻¹⁰ However, the dissociation probability is decreased by two orders of magnitude on W(110) for the same range of surface temperature and kinetic energies.¹¹,¹² The variation
in the reactivity between the two tungsten faces is due to the contribution of dynamical effects at distances relatively far from the surface and not to specific features of the surface active sites in which the atoms are eventually adsorbed after dissociation.<sup>13-15</sup>

Fe surfaces are also reported to undergo tremendous variations in the chemical activity when going from one crystal face to another. Under low pressures, the maximum value of the nitrogen molecule dissociation probability is obtained at the most open iron face, Fe(111). A less reactive face is the Fe(100), while the lowest reactivity was observed for the close-packed Fe(110) surface.⁷ Although Fe(111) is the most reactive iron face, the observed dissociation probability is still extremely low $S \approx 10^{-6}$ at incident energies of $E_i \approx 0.1$ eV. The sticking is increased by five orders of magnitude and approaches a saturation value of $S \approx 10^{-1}$ when the initial kinetic energy is increased up to $E_i \approx 4.3$ eV.<sup>16,17</sup>

Density Functional Theory (DFT) has been widely used for quantitative studies of the atomic and molecular nitrogen adsorption features on the three low-index iron faces.<sup>18-24</sup> A detailed study of the dynamics of the adsorption process using a DFT-based potential energy surface was recently performed for the $N_2$ on Fe(110) system.<sup>20</sup> A comparison between the adsorption kinetics of nitrogen on Fe(100) and Fe(111) has been also addressed in Ref.<sup>25</sup> To the best of our knowledge, there are no full-dimensional DFT-based dynamical studies for the adsorption dynamics of $N_2$ on the Fe(111) surface, in spite of this being the most reactive face. The general purpose of this work is to fill, at least partially, this gap.

Existing measurements provide a general picture of the nitrogen adsorption dynamical process on Fe(111).<sup>26-33</sup> For a wide range of substrate temperatures and for thermal energies of the incoming $N_2$ molecules, two molecular adsorption states have been identified, namely, a weakly bound state $\gamma$-$N_2$, with the molecular axis perpendicular to the surface plane, and a more strongly bound adsorption state $\alpha$-$N_2$, with the molecule adsorbed on the hollow site, parallel to the surface plane, and with stretched internuclear distance.<sup>26-33</sup> At low surface temperatures and when the surface is already saturated with $\gamma$-$N_2$, an additional adsorption state, $\delta$-$N_2$, has been experimentally observed with the molecule on top of a second
layer atom and perpendicular to the surface.³² Furthermore, a fourth adsorption state (ε-state) was reported in Ref. ³¹ Available theoretical studies predict similar adsorption states, in agreement with the experimental measurements.³⁴

In this context, the main goal of our study is to analyze in detail the *molecular adsorption dynamics* of N₂ on Fe(111). For low kinetic energies of the incident molecules, this is the dominant adsorption channel, as dissociative chemisorption was shown experimentally to be negligible.¹⁶ A preliminary short account of the static interaction between N₂ and Fe(111) surfaces was recently reported.³⁵ We found that the most favorable site for N atomic adsorption was the hollow site, with the atom placed at a distance of Z ≈ 0.5 Å from the surface and with an adsorption energy of E<sub>ads</sub> ≈ -5.82 eV. Concerning the N₂ molecular adsorption process and in accordance with the experimental observations, we find four molecular adsorption states, namely the γ, δ, α, and ϵ states. The four molecular adsorption sites play a competitive role in the dynamics of the N₂ molecular adsorption.

The article is arranged as follows: Sec.2 describes the computational details of the electronic structure calculations, the interpolation of the DFT data points that are used to build the continuous six dimensional (6D) potential energy surface (PES), by using the corrugation reducing procedure (CRP), and the methods that have been used to perform classical trajectory calculations. In Sec. 3, results from the DFT calculations as well as the dynamical study of the system are presented. The conclusion and outlook are included in Sec. 4.

# Methodology

## Electronic Structure Calculations

All ab-initio data are obtained using the DFT³⁶,³⁷ “Vienna Ab initio Simulation Program” (VASP) code³⁸⁻⁴⁰ that operates with a plane-wave basis set. The exchange-correlation energy is calculated with the generalized gradient approximation (GGA)⁴¹ using the Perdew-Wang energy functional (PW91).⁴² The (valence) electron-core interaction is described with
Projector Augmented Wave (PAW) pseudopotentials.³⁹,⁴⁰ Due to the magnetic nature of the Fe(111) surface, spin-polarized DFT calculations are performed. The energy cutoff in the plane-wave expansion is 400 eV; the fractional occupancies are determined through the broadening approach of Methfessel and Paxton with N=1 and σ = 0.2.

<figure><img src="image_2.png" /><figcaption>Figure 1: (a) Top view of the 2 × 2 unit cell of Fe(111) and (b) representation of the coordinates system used to describe the N₂/Fe(111) interaction.</figcaption></figure>

A 2 × 2 unit cell, as shown in Fig. 1.a, is used to model the periodicity of the Fe(111) surface. In our preliminary study of some features of the system,³⁵ the structural relaxation of the Fe(111) surface layers was investigated. Our results showed an independence between the interlayer relaxation and the number of layers in the slab, provided a minimum number of layers was used. We perform the present DFT calculations within the frozen surface approximation and using a slab consisting of 9-layers and separated by a vacuum layer of 15 Å width, which is considered to be large enough to avoid the interaction between neighbouring slabs. The two-dimensional Brillouin-zone integrations are performed by sampling the *k*-points with 6 × 6 × 1 Monkhorst-Pack grids,⁴³ which provide an adequate balance between numerical accuracy and computational burden.

The lattice constant for the Fe bulk crystal as well as the relaxation of the Fe(111) surface and the N₂ molecule equilibrium bond length were reported in a previous work by the authors.³⁵ Briefly, we obtained a theoretical lattice constant $a = 2.855$ Å, which is in good agreement with the experimental value of 2.867.⁴⁴ Our value of the theoretical lattice constant leads to an interlayer spacing of which leads to an interlayer spacing of $d = a/2\sqrt{3} \approx 0.824$ Å in
bulk. For a slab consisting of 9 layers, the net interlayer spacing values are $d_{12} \approx 0.74$ Å, $d_{23} \approx 0.64$ Å and $d_{34} \approx 0.90$ Å, where $d_{ij}$ is the distance between the i- and j- layers. The equilibrium internuclear distance for a nitrogen molecule in the gas phase is found to be $r_{eq} = 1.1125$ Å.

Here, we obtain the full set of DFT calculations for the six degrees of freedom of the nitrogen molecule on Fe(111). The coordinates of the $N_2$ molecule on the Fe(111) surface are illustrated in Fig 1.b, and defined in the following subsection.

## Geometry

The interaction energy of the $N_2$ molecule with the Fe(111) surface is described with a six-dimensional (6D) potential energy surface (PES) that depends on the molecular center of mass position $R_{CM} = (X_{CM}, Y_{CM}, Z_{CM})$ and the internal coordinates $(r, \theta, \phi)$, that respectively stand for the internuclear N-N distance, the polar angle relative to the surface normal, and the azimuthal angle relative to the X-axis on the surface plane (shown in Fig. 1.b). The PES is constructed from a grid of ab-initio energies, built out of 51 equidistant $Z_{CM}$ positions in the range $0 \text{ Å} \le Z_{CM} \le 5 \text{ Å}$ ($Z = 0$ Å is the topmost surface plane), 10 non-equidistant interatomic distances in the $0.7125 \text{ Å} \le r \le 2.125 \text{ Å}$ range, and 42 $(X_{CM}, Y_{CM}, \theta, \phi)$ configurations chosen following symmetry arguments. Thus, a total amount of 21420 electronic structure calculations were performed.

The Fe(111) surface geometry is illustrated in Fig. 1.a. *Top*, *hcp* and *fcc* sites correspond respectively to Fe atoms at the first, second and third plane; while *hollow*, *top-hcp* and *top-fcc* sites are respectively located at mid-distance between adjacent *top* sites, between *top* and *hcp* sites, and between *top* and *fcc* sites. Before continuing to the angular grid it is important to observe that, for example, a molecule with a given $\theta - \phi$ orientation does not necessarily find equivalent environments when positioned on different *hollow* sites. An additional +120° or -120° rotation may be required. This is also the case for *top-hcp* and *top-fcc* sites. Therefore it is important, for each XY site, to define a reference position on which $\phi$ values are applied.
These reference ($X_{CM}$, $Y_{CM}$) positions are i) top: (0,0), ii) hcp: ($a\frac{\sqrt{2}}{2}$, $a\frac{1}{\sqrt{6}}$), iii) fcc: ($a\sqrt{2}$, $a\frac{2}{\sqrt{6}}$), iv) hollow: ($a\frac{\sqrt{2}}{2}$,0), v) top-hcp: ($a\frac{\sqrt{2}}{2}$, $a\frac{2}{\sqrt{6}}$), vi) top-fcc: ($a\sqrt{2}$, $a\frac{1}{\sqrt{6}}$); with a the lattice constant defined in Fig. 1.a.

Regarding the 2D angular grid, it is made of

* Vertical $\theta = 0^{\circ}$ configurations for all six XY sites.

* Horizontal $\theta = 90^{\circ}$ configurations with azimuthal orientations $\phi = 0^{\circ}$ - $90^{\circ}$ with equivalent periodic symmetry of $60^{\circ}$ (i.e. $0^{\circ}$ is equivalent to $60^{\circ}$ and $90^{\circ}$ equal to $30^{\circ}$) for top, hcp and fcc sites, and $\phi = 0^{\circ}$ - $30^{\circ}$ - $60^{\circ}$ - $90^{\circ}$ for hollow, top-hcp and top-fcc sites.

* Tilted $\theta = 45^{\circ}$ configurations with azimuthal orientations $\phi = 30^{\circ}$ - $60^{\circ}$ - $90^{\circ}$ for top, hcp and fcc sites, and $\phi = 0^{\circ}$ - $90^{\circ}$ - $270^{\circ}$ for hollow, top-hcp and top-fcc sites.

## Interpolation Procedure

Based on the DFT calculations presented above, we build a global six-dimensional potential energy surface 6D-PES by means of a combination of interpolation techniques. A key interpolation tool is the corrugation reducing procedure (CRP).$^{45}$ CRP has been shown to provide very accurate results for the numerical interpolation of the interaction between diatomic gas-phase molecules and metal surfaces.$^{46}$ Briefly, the major problem for interpolating numerically the DFT data arises from the large variation of energy values due to the strong repulsion for some configurations at distances close to the surface. The clever idea behind the CRP method is that, for these very repulsive regions, the energy is mainly determined by the interaction between the individual atoms of the molecule and the surface. Then, by subtracting the atom-surface interaction, an auxiliary smoother function is produced which is much easier to interpolate than the highly corrugated 6D-PES. This interpolation function $I^{6D}(X, Y, Z, r, \theta, \phi)$ is given by

$$I^{6D}(X, Y, Z, r, \theta, \phi) = V^{6D}(X, Y, Z, r, \theta, \phi) - \sum_{i=1}^{2} V_{i}^{3D}(X_{i}, Y_{i}, Z_{i}),$$
where $(X, Y, Z, r, \theta, \phi)$ represent the six-degrees of freedom of the molecule on the surface as shown in Fig. 1.b, $V^{6D}(X, Y, Z, r, \theta, \phi)$ is the DFT interaction energy between the N₂ molecule at $(X, Y, Z, r, \theta, \phi)$ and the surface, and $V_i^{3D}(X_i, Y_i, Z_i)$ is the interaction between atom i-th of the molecule at $(X_i, Y_i, Z_i)$ and the surface. The calculations of the 3D-PES for the atom-surface interaction should be therefore addressed as well. The 3D-PES is highly corrugated and a similar CRP method can be used. The 3D interpolation function can be written in this case as

$$I^{3D}(X_i, Y_i, Z_i) = V_i^{3D}(X_i, Y_i, Z_i) - \sum_{j=1}^{n} V^{1D}(|R_i - R_j|),$$

such that $V^{1D}(|R_i - R_j|$ is a one-dimensional potential describing the interaction between the N atom at $\bar{R}_i = (X_i, Y_i, Z_i)$ and the j-th Fe surface atom at $\bar{R}_j = (X_j, Y_j, Z_j)$ and $V_i^{3D}(X_i, Y_i, Z_i)$ is the actual DFT calculated potential. In our case, we use as one-dimensional $V^{1D}$ potential the DFT calculation that describes the approach of an N atom to the top site of the Fe(111) surface.

In order to evaluate the $I^{6D}$ for an arbitrary molecular configuration $(R_0^{CM}, r_0, \theta_0, \phi_0)$, a spline interpolation is carried out along the $Z_{CM}$ and $r$ coordinates, followed by 2D Fourier interpolations, first on the angular internal coordinates $\theta$ and $\phi$, and afterwards within the $X_{CM} - Y_{CM}$ plane. The interpolated value $I^{6D}(R_0^{CM}, r_0, \theta_0, \phi_0)$ can easily be converted back into an interpolated energy $V^{6D}(R_0^{CM}, r_0, \theta_0, \phi_0)$. Note that the 2D Fourier expansion in the $X_{CM} - Y_{CM}$ plane is here chosen over a 2D Periodic spline because it analytically respects the surface hexagonal symmetry. The interpolation of the 3D atomic PES is performed in a similar way, using a spline method for interpolating the Z variable and a Fourier expansion to interpolate over the $(X, Y)$ plane.

In order to test the accuracy of the CRP interpolation method over the six-degrees of freedom of the N₂/Fe(111) system, we perform electronic structure calculations for different molecule configurations that are not included in the CRP construction procedure and compare these DFT results with the interpolated energy values. Figure 2 shows an example of the PES
accuracy when interpolating over the molecular cartesian coordinates (Z, X, and Y). The left side in Fig. 2 shows the dependence of the potential energy on the vertical distance from the surface, Z, for the equilibrium internuclear distance in the gas phase ($r_{eq}$), and for different molecular configurations that have been included in the interpolation procedure. In the right side of Fig. 2, we present a test over the surface coordinates X and Y, for various values of the remaining coordinates Z, r, $\theta$, and $\phi$. Another example of the accuracy obtained in the interpolation procedure for the polar angle $\theta$ and the azimuthal angle $\phi$ is shown in the left and right sides of Fig. 3, respectively. The DFT data shown in Fig. 3 are again, out of the grid used to build the 6D-PES and the internuclear distance is fixed to $r = r_{eq}$. Results are shown for several high symmetry sites. In general, at distances relatively far from the surface, our interpolated values compare well with the ab-initio calculations. Nevertheless, due to the large increase in the corrugation at distances close to the surface, some errors in the interpolation do exist. Fortunately, for the configurations of relevance for the adsorption dynamics, the errors are well controlled and lower than 100 meV in most cases.

<figure><img src="image_3.png" /><figcaption>Figure 2: Comparison between the interpolated (solid lines) and DFT (symbols) energy values for N$_2$/Fe(111). The left side illustrates the dependence of the energy on the distance of the molecular center of mass from the surface Z, for fixed values of X,Y, $\theta$, and $\phi$ (shown in the legend). The internuclear distance is set at the equilibrium distance $r_{eq} = 1.1125$ Å. The right side shows the dependence of the potential energy on the surface (X,Y) coordinates along different paths, as illustrated in the inset of the Figure for various Z,r, $\theta$ and $\phi$ values.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Comparison between the interpolated (solid lines) and DFT (symbols) energy values, for fixed X,Y, Z and equilibrium internuclear distance $r_{eq} = 1.1125$ Å. The left side illustrates the dependence on the polar angle $\theta$, while the left side shows the dependence on the azimuthal angle $\phi$.</figcaption></figure>

## Quasi-Classical Trajectory Calculations

Quasi-classical trajectory (QCT) calculations are performed on the continuous 6D-PES, in order to study the $N_2$ molecular adsorption dynamics on the clean Fe(111) surface. We include the initial zero point energy (ZPE) of the molecule in the calculations and set it to be ZPE = 142 meV. The path for all trajectories starts with the molecule at its internuclear equilibrium distance in the gas phase ($r_{eq} = 1.1125$ Å) and with an initial distance from the surface of $Z_{cm} = 5.0$ Å. For such initial conditions, the 6D-PES values tend to zero. A conventional Monte Carlo procedure is used for sampling the initial molecular coordinates (X,Y) and the molecular orientation ($\theta$ and $\phi$). For surface motion, the generalized Langevin oscillator (GLO)$^{47-49}$ model is applied to include temperature effects. In this model, surface motion is represented in terms of a three-dimensional harmonic oscillator. In order to consider the coupling and energy-exchange of the surface with the bulk, an additional three-dimensional ghost oscillator is coupled to the surface oscillator. The latter is subjected to friction and random forces related to each other through the second fluctuation-dissipation theorem. The ghost particles allow us to represent the bulk of the solid as a thermal bath at
the chosen temperature.⁵⁰ In our calculations, we distinguish among three different events: (i) reflection, when the molecule center of mass reaches for a second time the initial starting distance of 5.0 Å, with positive velocity pointing towards the vacuum; (ii) dissociation, when the molecule bond length reaches the value of $r_{diss} = 2.1125$ Å; and (iii) molecular trapping, when the molecule is neither reflected nor dissociated after 10 ps. The calculations are performed for normal incidence within a wide interval of initial kinetic energies for the molecule, ranging between 50 meV and 5.0 eV. For each value of the initial kinetic energy (Ei), 100 000 trajectories are used to obtain sufficiently good statistics, with a maximum trajectory integration time of 10 ps. The surface sampling area used for representing the initial impact (X,Y) positions on the surface is shown as the shaded area in Fig 4.

<figure><img src="image_5.png" /><figcaption>Figure 4: Representation of the Fe(111) surface used in the quasi-classical trajectory calculations. Shaded area illustrates the sampling area of initial (X,Y) positions, where 'a = 2.855 Å' is the lattice parameter.</figcaption></figure>

# Results and Discussion

## Static Adsorption States

Let us start this section by presenting the results of our electronic structure calculations for the static interaction between the N₂ molecule and the Fe(111) surface. Figure 6 illustrates the interaction between the N₂ molecule and the Fe(111) surface as a function of the molecule-surface distance Z. The reference energy is taken to be that of the N₂ molecule
at 5 Å from the surface and at the equilibrium bond length ($r_{eq} = 1.1125$ Å). Our DFT calculations yield four possible adsorption states for N₂ on the Fe(111) surface. All of them are equilibrium molecular adsorption states in the multidimensional PES. A graphical representation of their geometry at the surface is shown in Fig. 5, together with the adsorption geometry of atomic N on Fe(111).

<figure><img src="image_6.png" /><figcaption>Figure 5: (a) Top view and (b) side view of the five nitrogen adsorption states on a clean Fe(111) surface. From left to the right, the molecular adsorption states $\gamma$, $\delta$, $\alpha$, and $\varepsilon$, and the atomically adsorption state called $\beta$.</figcaption></figure>

On one hand, there are two molecular adsorption states with the molecular axis perpendicular to the surface plane and with $r = r_{eq}$, on top of the first ('top' site) and second ('hcp' site) layer atoms. These states are usually called the $\gamma$- and $\delta$-state, respectively. Molecules adsorbed into the $\gamma$-state (black symbols) are associated with an adsorption energy well of $E_{ads} \approx 0.46$ eV, at $Z \approx 2.5$ Å from the surface. Molecular adsorption into the $\gamma$-state is detected experimentally in all the previous studies,$^{26-32}$ as well as predicted theoretically,$^{34,51}$ with an adsorption energy ranging between 0.2 and 0.4 eV.$^{26,27,31,32}$ The molecular nitrogen adsorption into the $\delta$-state (dark-grey symbols) is characterized by an adsorption energy of $E_{ads} \approx 0.45$ eV, slightly weaker than the $\gamma$-state and with molecular center of mass distance to the surface of $Z \approx 1.7$ Å. Regarding the $\delta$-state it is worth mentioning that the perpendicular distance of N₂ center of mass to the second layer Fe atom right below the N₂ molecule is 2.46 Å, also close to the value obtained for the $\gamma$-state. The $\delta$-state has been observed experimentally at low surface temperature and when the surface is saturated with $\gamma$-N₂.$^{32}$
It was found to be less stable than the γ-state.

On the other hand, in Fig. 6, there are two other adsorption states arising from our present study, with molecular axis parallel to the surface plane ($\theta = 90^\circ$), and with stretched molecular bond length. The α-state, which is denoted by green symbols in Fig 6, represents the deepest adsorption well of our DFT calculations. The adsorption energy is $E_{ads} \approx 0.64$ eV and the molecular center of mass is located at $Z = 1.2$ Å, on the hollow site and with an internuclear distance increased by about 0.1 Å ($r = 1.2$ Å). The α-state is well known experimentally$^{26-32}$ and theoretically,$^{34,51}$ with an adsorption energy larger than the one of the γ-state by about 200 meV.$^{26,27,31,32}$ In our calculations, we find the same energy difference. Finally, the ε-state is an additional adsorption state found in our calculations, which has been reported before only in Ref.$^{31,34,51}$ It is denoted by light-grey symbols in Fig. 6. Characterized by an adsorption energy of $E_{ads} \approx 0.45$ eV, the ε-state is identified close to the surface at $Z \approx 0.6$ Å, on top of the third (fcc) layer atoms, with an internuclear distance increased from the equilibrium distance by about 0.2 Å ($r = 1.3$ Å). In the ε-state, the two N atoms in the molecule are oriented in the azimuthal plane in such a way that they avoid the strong repulsion due to the first layer atoms ($\phi = 60^\circ$) and are placed close to the atomic adsorption site (hollow site). Freund et al. reported that the adsorption energy of the linearly bonded γ-state is very close to that of the ε-state.$^{31}$ We reach the same conclusion from our calculations. Both the α- and the ε-states are potential candidates to be intermediate states for dissociative adsorption at very low impact energy.

In Fig 7, a selection of two-dimensional (r,Z) cuts for the four local minima of the 6D-PES, as identified before, are introduced. The molecular orientation and the surface coordinates are shown in the inset of each contour plot. An energy scale ranging from -0.7 eV (red color) to 2.5 eV (blue color) is used to show the different attractive and repulsive regions of the potential for each elbow plot, with the zero energy value indicated by dashed lines. The upper panel in Fig. 7 represents the adsorption states γ (Fig 7.a) and δ (Fig 7.b), for which the molecules are adsorbed vertically on top of the first and second layer atoms of the sur-
<figure><img src="image_7.png" /><figcaption>Figure 6: Representation of the four local minima for N₂/Fe(111) system. The figure shows the dependence of the ab-initio potential energy on the distance of the molecular center of mass from the surface Z for fixed X,Y, r, θ and φ. The values of r, θ and φ are shown in the inset.</figcaption></figure>

face, respectively. The lower panel shows the adsorption states α (Fig 7.c) and ε (Fig 7.d), for which the molecules are adsorbed parallel to the surface. Except for the molecular adsorption path toward the γ-state, Fig 7 shows an entrance energy barrier into all adsorption states δ (≈ 200 meV), α (≈ 100 meV) and ε (≈ 400 meV). Figure 7 corresponds, however, to two-dimensional contour plots, i.e., to molecules that are forced to reach a particular site with fixed configuration. Due to the interaction with the surface, the initial configurations of the molecule change continuously through the molecule dynamical pathway for adsorption. Therefore, a conclusion based only on the static analysis of the PES for different but fixed configurations can be misleading.¹³,⁵²,⁵³ In the next subsection, we show the results of dynamical calculations performed to characterize the molecular preferred adsorption states over the continuous 6D-PES.
<figure><img src="image_8.png" /><figcaption>Figure 7: Contour plots for interpolated 2D (r,Z) cuts of the possible adsorption states ($\gamma$, $\delta$, $\alpha$, and $\varepsilon$). The molecular orientation and position over the unit cell is depicted in each contour plot. The energy scale on the right side is given in eV, ranging from -0.7 eV (red) to 2.5 eV (blue). The white regions show potential energy values higher than 2.5 eV. The dashed white lines correspond to the zero potential energy values. The contour lines are separated by an interval of 0.2.</figcaption></figure>
Adsorption Dynamics

In the following we present our results for the six-dimensional quasi-classical trajectory calculations for the N₂ adsorption dynamics on the clean Fe(111) surface. Figure 8 shows the molecular adsorption probabilities as a function of the initial kinetic energy of the molecule. Results are shown for normal incidence and for surface temperature, Tₛ, varying from 0 to 520 K. The adsorption probability decreases with increasing Tₛ, in accordance with experimental observations.⁸,²⁷ Molecules can not accommodate on the surface unless they are able to dissipate energy into the lattice. The energy difference between the incoming molecule and the substrate plays a crucial role in the adsorption process. In most of the energy range shown in Fig 8, the molecules are hot compared to the surface (i.e. Eᵢ > k<sub>B</sub> Tₛ). As Tₛ increases, less energy is usually dissipated into the lattice in average, and, consequently, the probability of trapping into the adsorption well is decreased.⁴⁹ If we focus on the dependence on the initial molecular kinetic energy, Eᵢ, the adsorption probability tends to increase with decreasing Eᵢ. Later in this manuscript we will show that the reason for that is the increasing contribution of different adsorption channels. For Eᵢ below 200 meV, however, the adsorption probability starts decreasing with decreasing Eᵢ due to the complex energy landscape for some particular configurations.

Details of the adsorption dynamics process for N₂ on Fe(111) with Tₛ = 300 K, are shown in Figs. 9 and 10, for normal incidence and low initial kinetic energies, Eᵢ = 50 and 200 meV. Snapshot of the molecular center of mass coordinates on the surface plane as well as of the orientation of the adsorbed molecules for various Z values are represented in the upper panel of the figures. Only those trajectories that correspond to molecules eventually adsorbed are shown. The lower panel corresponds to the final adsorption sites, orientations, heights and internuclear distances of the adsorbed molecules. For Eᵢ = 50 meV (Fig. 9), all adsorbed molecules are spread over the unit cell at Z=2.6 Å. Closer to the surface, all trajectories tend to change their direction and the molecules feel a strong attraction to the first layer atoms at Z=2.2 Å. On their way further toward the surface, more than 40 % of
<figure><img src="image_9.png" /><figcaption>Figure 8: Molecular adsorption probability as a function of the molecular initial incidence energy, for surface temperatures 0, 300 and 520 K and normal incidence.</figcaption></figure>

the molecules are attracted to the adsorption well of the first layer atoms, i.e., the γ-state. Moreover, repulsion due, in part, to the interaction with the first layer atoms will force the molecule to escape toward other positions in the unit cell. However, the kinetic energy of most of these molecules is insufficient to overcome the energy barrier to approach further the surface. Therefore, only few trajectories can approach as close as $Z = 1.8$ Å from the surface and most of them have an orientation parallel to the surface. This dynamical evolution can explain the final position of the adsorbed molecules: the majority of the molecules are adsorbed on the *top* site of the first layer atoms and very few molecules are finally adsorbed on the *hollow* site. This conclusion is consistent with the histograms in the left panel of the lower graph, where most of the trajectories will be adsorbed above 2.0 Å from the surface, with molecular orientations close to the perpendicular adsorption state and with internuclear distances similar to that of the γ-state. These results are consistent with the fact that there is no energy barrier for the molecules to reach the adsorption well of the first layer atoms from the gas-phase, which agrees with the experimental studies.²⁷,²⁸,³¹

Figure 8 shows that the maximum value of the adsorption probability is found at $E_i = 200$
<figure><img src="image_10.png" /><figcaption>Figure 9: Evolution of the coordinates of the N₂ molecules that become eventually adsorbed. 100 000 trajectories, with initial impact energy *E*ᵢ = 50 meV and with surface temperature *T*ₛ = 300 K are considered. (a) XY distribution of the center of mass of N₂ for different Z values. Above each sampling area, the percentage of adsorbed molecules that reach this particular Z distance, N<sub>Z</sub><sup>ads</sup>, is written. The polar angle θ-distribution is presented below the sampling area. (b) The final XY positions of the adsorbed molecules are represented in the right panel. The left panel shows the distribution of heights (Z), internuclear distances (r) and polar angles (θ) for the adsorbed molecules. The total molecular adsorption probability under these conditions is P=0.81.</figcaption></figure>
meV, for which the molecules can reach deeper positions toward the surface and have access to potential wells other than the γ-state. Similarly as in Fig. 9, we follow the dynamics of the adsorption molecules with an initial kinetic energy of $E_i = 200$ meV in Fig 10. In this case, approximately 17% of the molecules fall in the adsorption well of the top site (γ-state) before reaching $Z \approx 2.0$ Å. Most of the remaining trajectories start to point toward the hollow site of the unit cell. In addition, the number of trajectories is roughly preserved from $Z \approx 1.6$ up to $Z \approx 1.2$. Further movement toward the surface is restricted by a large potential energy barrier, and only few molecules can reach distances below 1.0 Å. The positions of the center of mass for these molecules will be located on top of the third layer atoms, the 'fcc site' of the unit cell, and they will be oriented parallel to the surface. In addition, we can notice the difference between the final adsorption states for the two initial incidence kinetic energies. In the lower panel of Fig. 10, there are more adsorbed molecules in the so-called α-state, and even a few number of trajectories will end up adsorbed into the ε-state. Comparison between Fig. 9 and Fig. 10 shows that the molecules cross through the γ-state before being adsorbed into the α-state. Indeed, there is a separation energy barrier of 200 meV between γ and α states. A similar conclusion was derived from previous experimental studies.<sup>28-31</sup>

Figure 11 shows the variation of the N<sub>2</sub>/Fe(111) relative adsorption probability into top, hcp, fcc and hollow sites of the unit cell, corresponding to the possible adsorption states γ, δ, α, and ε, respectively, with the initial molecular incidence energy at $T_s = 300$ K. At very low energy, in agreement with experimental observations,<sup>26,27</sup> adsorption into the γ-state has a relative adsorption probability close to one ($\sigma_γ \approx 1$). Adsorption into the α-state has a relative adsorption probability $\sigma_α \approx 10^{-2}$ for the same range of energies. As the initial molecular kinetic energy increases, a sharp decrease is observed for the adsorption into the γ-state, which eventually vanishes at high initial impact energy. Nevertheless, $\sigma_α$ will increase with increasing molecule initial energy, until reaching a saturation value of $\sigma_α \approx 0.55$ at $E_i = 0.8$ eV. Moreover, the adsorption into the δ-state is negligible over the full energy range. Such behaviour was previously extracted from experiments and static DFT calcu-
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10: Same as Fig. 9 with an initial incidence energy $E_i = 200$ meV. The total molecular adsorption probability under these conditions is $P=0.86$.</figcaption></figure>
lations in Ref. ³²,³⁴ According to the PES characteristics, direct adsorption into the ε-state
is associated with an entrance barrier larger than roughly 400 meV. An increase in σε with
increasing initial molecular impact energy is therefore observed, starting at 500 meV until
reaching a saturation value of σε ≈ 0.40 at Ei = 2.25 eV.

<figure><img src="image_13.png" /><figcaption>Figure 11: Relative molecular adsorption probability to *top*, *hcp*, *fcc* and *hollow* sites of the unit cell (γ, δ, ε, and α states, respectively) as a function of the initial incidence energy, *E*<sub>i</sub>, of the molecule, with surface temperature *T*<sub>s</sub>= 300 K.</figcaption></figure>

Let us remark at this point that we have restricted our analysis of the adsorption dynamics
to non-dissociative molecular adsorption. Dissociative adsorption, however, is a competitive
(and even dominant) adsorption channel for high temperatures and/or high incident kinetic
energies of the nitrogen molecule. For 3 eV initial kinetic energy, for instance, the dissocia-
tive adsorption probability is about 0.05.⁹ Our dynamical calculations include this additional
channel. We have preferred, however, to focus the current work only on the molecular ad-
sorption dynamics, which is a topic sufficiently rich, and leave the analysis of the dissociation
dynamics for a succeeding article.
Conclusion

In this work, we have presented a theoretical study for the interaction of N₂ on clean Fe(111) surfaces. First we performed ab-initio density functional theory calculations, taking into account the six-degrees of freedom (X, Y, Z, r, θ, φ) of the molecular system. The corrugation reducing procedure is used as a numerical interpolation method of the DFT data. The continuous 6D-PES shows a high accuracy. Four local adsorption minima have been obtained from our DFT calculations. Two adsorption states with the molecule perpendicularly oriented to the surface plane and on top of the first and the second layer atoms, have been found. Two additional states are reported, with the molecular axis oriented parallel to the surface plane, with a stretched bond length, and with the molecule located on the hollow site and on top of the third layer atoms.

We studied the molecular adsorption dynamics in the 6D-PES by performing quasi-classical trajectory calculations. The QCT calculations were performed for a wide range of initial impact energies, and using different surface temperatures. The molecule adsorption probability is strongly dependent on the surface temperature as well as on the initial impact energy. At very low energy, details of the dynamics indicate no entrance *energy* barrier for the molecules to reach the surface. Most of the molecules will be adsorbed into the first layer atoms. However, the adsorption probability into the molecular states, which are parallel to the surface and that are associated with an activation *energy* barrier, will be much higher for increasing initial impact energy. Our results show good agreement with previous experimental studies of N₂ on Fe(111), helping to understand the complex dynamics of the N₂ molecular adsorption on the Fe(111) surface. The analysis of the dynamics for N₂/Fe(111) will be continued in a forthcoming article where we will address the dissociation dynamics of the system.
## Acknowledgements

This work has been partially supported by the Basque Departamento de Educación, Universidades e Investigación, the University of the Basque Country UPV/EHU (Grant No. IT-756-13) and the Spanish Ministerio de Economía y Competitividad (Grant No. FIS2016-76471-P). M.A.N. acknowledges financial support by the Ministerio de Economía y Competitividad (Grant No. BES-2011-045536). Computational resources were provided by the DIPC Computing Center.

## References

(1) Somorjai, G. A. *Introduction to surface chemistry and catalysis*, 3rd ed.; John Wiley & Sons, Inc.: New York, 1984.

(2) Díez Muiño, R.; Busnengo, H. F. *Dynamics of gas/surface interactions: Atomic-level understanding of scattering processes at surfaces*; Springer Series of Surface Sciences: Berlin, 2013; Vol. 50.

(3) Hansen, F. Y.; Henriksen, N. E.; Billing, G. D.; Guldberg, A. *Surface Science* **1992**, *264*, 225.

(4) Hellman, A. et al. *The Journal of Physical Chemistry B* **2006**, *110*, 17719.

(5) Kitano, M.; Kanbara, S.; Inoue, Y.; Kuganathan, N.; Sushko, P. V.; Yokoyama, T.; Hara, M.; Hosono, H. *Nature Communications* **2015**, *6*, 6731.

(6) Iyngaran, P.; Madden, D. C.; Jenkins, S. J.; King, D. A. *Proceedings of the National Academy of Sciences* **2011**, *108*, 925.

(7) Spencer, N.; Schoonmaker, R.; Somorjai, G. *Journal of Catalysis* **1982**, *74*, 129.

(8) Rettner, C. T.; Schweizer, E. K.; Stein, H.; Auerbach, D. J. *Physical Review Letters* **1988**, *61*, 986.
(9) Rettner, C. T.; Stein, H.; Schweizer, E. K. *The Journal of Chemical Physics* **1988**, *89*, 3337.

(10) Pétuya, R.; Plötz, P.-A.; Crespos, C.; Larregaray, P. *The Journal of Physical Chemistry C* **2014**, *118*, 21904.

(11) Pfnür, H. E.; Rettner, C. T.; Lee, J.; Madix, R. J.; Auerbach, D. J. *The Journal of Chemical Physics* **1986**, *85*, 7452.

(12) Goikoetxea, I.; Juaristi, J. I.; Alducin, M.; Díez Muiño, R. *Journal of Physics: Condensed Matter* **2009**, *21*, 264007.

(13) Alducin, M.; Díez Muiño, R.; Busnengo, H. F.; Salin, A. *The Journal of Chemical Physics* **2006**, *125*, 144705.

(14) Alducin, M.; Díez Muiño, R.; Busnengo, H. F.; Salin, A. *Physical Review Letters* **2006**, *97*, 056102.

(15) Martin-Gondre, L.; Crespos, C.; Larregaray, P.; Rayez, J. C.; van Ootegem, B.; Conte, D. *The Journal of Chemical Physics* **2010**, *132*, 204501.

(16) Rettner, C. T.; Stein, H. *Physical Review Letters* **1987**, *59*, 2768.

(17) Rettner, C. T.; Stein, H. *The Journal of Chemical Physics* **1987**, *87*, 770.

(18) Mortensen, J.; Ganduglia-Pirovano, M.; Hansen, L.; Hammer, B.; Stoltze, P.; Nørskov, J. *Surface Science* **1999**, *422*, 8.

(19) Riikonen, S.; Foster, A. S.; Krasheninnikov, A. V.; Nieminen, R. M. *Physical Review B* **2009**, *80*, 155429.

(20) Goikoetxea, I.; Alducin, M.; Díez Muiño, R.; Juaristi, J. I. *Physical Chemistry Chemical Physics* **2012**, *14*, 7471.
(21) Goikoetxea, I.; Juaristi, J. I.; Díez Muiño, R.; Alducin, M. *Physical Review Letters* **2014**, *113*, 066103.

(22) Goikoetxea, I.; Alducin, M.; Díez Muiño, R.; Juaristi, J. I. *Physical Chemistry Chemical Physics* **2015**, *17*, 19432.

(23) Pick, S.; Légaré, P.; Demangeat, C. *Physical Review B* **2007**, *75*, 195446.

(24) Wang, T.; Tian, X.; Yang, Y.; Li, Y.-W.; Wang, J.; Beller, M.; Jiao, H. *The Journal of Physical Chemistry C* **2016**, *120*, 2846.

(25) Panczyk, T. *The Journal of Physical Chemistry C* **2007**, *111*, 3175.

(26) Bozso, F.; Ertl, G.; Grunze, M.; Weiss, M. *Journal of Catalysis* **1977**, *49*, 18.

(27) Ertl, G.; Lee, S.; Weiss, M. *Surface Science* **1982**, *114*, 515.

(28) Grunze, M.; Golze, M.; Hirschwald, W.; Freund, H. J.; Pulm, H.; Seip, U.; Tsai, M. C.; Ertl, G.; Küppers, J. *Physical Review Letters* **1984**, *53*, 850.

(29) Tsai, M.-C.; Ship, U.; Bassignana, I.; Küppers, J.; Ertl, G. *Surface Science* **1985**, *155*, 387.

(30) Whitman, L. J.; Bartosch, C. E.; Ho, W.; Strasser, G.; Grunze, M. *Physical Review Letters* **1986**, *56*, 1984.

(31) Freund, H.-J.; Bartos, B.; Messmer, R.; Grunze, H.; Kuhlenbeck, H.; Neumann, M. *Surface Science* **1987**, *185*, 187.

(32) Grunze, M.; Strasser, G.; Golze, M. *Applied Physics A* **1987**, *44*, 19.

(33) Alstrup, I.; Chorkendorff, I.; Ullmann, S. *Journal of Catalysis* **1997**, *168*, 217–234.

(34) Mortensen, J.; Hansen, L.; Hammer, B.; Nørskov, J. *Journal of Catalysis* **1999**, *182*, 479.
(35) Nosir, M. A.; Martin-Gondre, L.; Bocan, G.; Díez Muiño, R. *Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms* **2016**, *382*, 105.

(36) Hohenberg, P.; Kohn, W. *Physical Review* **1964**, *136*, B864.

(37) Kohn, W.; Sham, L. J. *Physical Review* **1965**, *140*, A1133.

(38) Kresse, G.; Furthmüller, J. *Physical Review B* **1996**, *54*, 11169.

(39) Kresse, G.; Joubert, D. *Physical Review B* **1999**, *59*, 1758.

(40) Blöchl, P. E. *Physical Review B* **1994**, *50*, 17953.

(41) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. *Physical Review B* **1992**, *46*, 6671.

(42) Perdew, J. P.; Wang, Y. *Phys. Rev. B* **1992**, *45*, 13244.

(43) Monkhorst, H. J.; Pack, J. D. *Physical Review B* **1976**, *13*, 5188.

(44) Kohlhaas, R.; Donner, P.; Schmitz-Pranghe, N. *Zeitschrift für Angewandte Physik* **1967**, *23*, 245.

(45) Busnengo, H. F.; Salin, A.; Dong, W. *The Journal of Chemical Physics* **2000**, *112*, 7641.

(46) Díaz, C.; Olsen, R. A.; Busnengo, H. F.; Kroes, G. J. *The Journal of Physical Chemistry C* **2010**, *114*, 11192.

(47) Adelman, S. A.; Doll, J. D. *The Journal of Chemical Physics* **1976**, *64*, 2375.

(48) Tully, J. C. *The Journal of Chemical Physics* **1980**, *73*, 1975.

(49) Busnengo, H. F.; Di Césare, M. A.; Dong, W.; Salin, A. *Phys. Rev. B* **2005**, *72*, 125411.
(50) Martin-Gondre, L.; Bocan, G.; Alducin, M.; Juaristi, J.; Díez Muiño, R. *Computational and Theoretical Chemistry* **2012**, *990*, 126.

(51) Tománek, D.; Bennemann, K. H. *Phys. Rev. B* **1985**, *31*, 2488.

(52) Díaz, C.; Vincent, J. K.; Krishnamohan, G. P.; Olsen, R. A.; Kroes, G. J.; Honkala, K.; Nørskov, J. K. *Phys. Rev. Lett.* **2006**, *96*, 096102.

(53) Gross, A.; Scheffler, M. *Phys. Rev. B* **1998**, *57*, 2493–2506.