Non-adiabatic effects during the dissociative adsorption of O2 at Ag(111)? A first-principles divide and conquer study

We study the gas-surface dynamics of O2 at Ag(111) with the particular objective to unravel whether electronic non-adiabatic effects are contributing to the experimentally established inertness of the surface with respect to oxygen uptake. We employ a first-principles divide and conquer approach based on an extensive density-functional theory mapping of the adiabatic potential energy surface (PES) along the six O2 molecular degrees of freedom. Neural networks are subsequently used to interpolate this grid data to a continuous representation. The low computational cost with which forces are available from this PES representation allows then for a sufficiently large number of molecular dynamics trajectories to quantitatively determine the very low initial dissociative sticking coefficient at this surface. Already these adiabatic calculations yield dissociation probabilities close to the scattered experimental data. Our analysis shows that this low reactivity is governed by large energy barriers in excess of 1.1 eV very close to the surface. Unfortunately, these adiabatic PES characteristics render the dissociative sticking a rather insensitive quantity with respect to a potential spin or charge non-adiabaticity in the O2-Ag(111) interaction. We correspondingly attribute the remaining deviations between the computed and measured dissociation probabilities primarily to unresolved experimental issues with respect to surface imperfections.


Introduction
A predictive materials science modelling based on microscopic understanding requires a thorough knowledge of all underlying elementary processes at the atomic scale. The (dissociative) adsorption of individual oxygen molecules at metal surfaces is such an elementary process that is of crucial relevance for a wide range of applications, with heterogeneous catalysis forming just one prominent example. Recent work on the dissociative adsorption of O 2 at Al(111) has severely challenged the prevalent understanding of this process [1,2,3,4]. Significant spin non-adiabatic effects, i.e. a coupled electronic-nuclear motion beyond the standard Born-Oppenheimer approximation, have been proposed as rationalization for the experimentally observed low initial adsorption probability. This behaviour has been associated to the low density-of-states (DOS) at the Fermi-level of the simple metal surface, which prevents an efficient spin quenching through the tunnelling of electrons between substrate and adsorbate. This hinders the transition from the spin-triplet ground state of gas-phase O 2 into the singlet state upon adsorption at the metal surface. The low sticking coefficient can then be seen as a result of similar spin-selection rules that are well known from gas-phase chemistry.
In particular from the point of view of heterogeneous catalysis, interest naturally shifts to low-index noble metal surfaces. Only recently, a spin-asymmetry in electronhole pair excitation spectra for O 2 at Pd(100) have further fortified the previously proposed DOS-related mechanism [5]. On the other hand, in contrast to palladium, the DOS at the Fermi-level is equally low for aluminium and coinage metals. This raises the question whether similar spin non-adiabatic effects significantly contribute to their established low reactivity with respect to O 2 dissociation [6,7,8,9,10,11,12,13]. Here, we aim to contribute to this with a first-principles based modelling of the O 2 gas-surface dynamics specifically at the Ag(111) surface, for which particularly low adsorption probabilities have been reported experimentally [10,14]. Earlier simplified models based on low dimensional analytic potentials have already tried to understand the nature of such low reactivity in terms of charge transfer limitations from the Ag(111) surface to the impinging O 2 molecule [15,16]. Notwithstanding, on Ag(100) recent molecular dynamics (MD) calculations based on an ab initio six-dimensional potential energy surface (PES) showed that the comparable absence of O 2 dissociation on that surface can be fully explained within an adiabatic picture [17]. There, the system is simply characterized by large dissociation barriers of about 1.1 eV that appear when the molecule is close to the surface. This motivates a similar ab initio investigation also for Ag(111), with the objective to settle the question up to what extent the low reactivity on the Ag(111) surface is really a signature of electronic non-adiabaticity or only the result of large dissociation barriers as found for Ag(100).
The overall structure of the paper is as follows. The next section provides a detailed description of the different steps employed to model the dynamics of O 2 on Ag(111) within a divide-and-conquer approach. Section 3 commences with a brief review of the existing experimental data, emphasizing controversies with respect to the absolute order of magnitude of dissociative sticking for incidence energies E i below 1 eV. We then proceed with a discussion of the quality and relevant static properties of the computed adiabatic PES, revealing as a central characteristic barriers to dissociation larger than 1.1 eV at close distances to the surface. The ensuing classical trajectory calculations demonstrate that these barriers lead to dissociation probabilities that already fall within the large error bars set by the scattered experimental data. With heights that are larger than the gas-phase O 2 spin singlettriplet gap (0.98 eV [18]) and the O 2 electron affinity (0.44 eV [18]) these barriers furthermore completely mask any signs of a potential spin or charge non-adiabaticity in the O 2 -Ag(111) interaction in the dissociative sticking coefficient. The remaining deviations of the computed sticking curve to the measurements for E i < 1 eV are thus more likely the result of the unresolved issues in the experimental data, e.g. with respect to surface imperfections as discussed by the experimentalists [14,19].

Methodology
Methodologically, the very low dissociation that characterizes this system makes the use of on-the-fly ab initio molecular dynamics (AIMD) to determine the sticking coefficient [20] computationally unaffordable. Instead, our simulations are based on a divide-and-conquer approach [21,22,23], in which an accurate PES is first constructed from first-principles, and next, classical MD calculations are performed using a continuous representation of this PES to describe the molecule-surface interaction. Once the former is obtained, the latter come at negligible computational cost. In contrast to AIMD, the divide-and-conquer approach thus enables the calculation of a large number of trajectories. Even for very small sticking probabilities and various initial conditions of molecules in the impinging molecular beam (E i , incidence angle etc.), sufficient statistical sampling can thus be performed. Additionally, detailed information on the PES properties, which allows for a better rationalization of the factors ruling the gas-surface dynamics, provide another advantage. Simulations based on this methodology have been able to reproduce not only reactive probabilities, but also more subtle quantities like the rovibrational state distributions of scattered molecules [24,25].
The key ingredient of this impressive success that ensures reliable quantitative results is the use of multi-dimensional ab initio-based PESs. In this respect, several interpolation methods to obtain continuous representations in the context of gassurface dynamics have been developed in the last years, reaching energy precisions on the order of meV in parts most relevant for the former [26,27,28,29,30]. In the present work, the PES interpolation is based on neural networks (NN) [28,29,30] that we adapt to the specific symmetry of the O 2 /Ag(111) system as explained in section 2.2. In order to keep the multi-dimensionality of the O 2 /Ag(111) PES in a tractable limit, we follow common practice and adopt the frozen surface approximation. This approximation has been successfully applied to simulate the reactivity of various diatomic molecules (H 2 , N 2 , O 2 , . . . ) on metal surfaces [1,24,31,32,33]. One rationalization for this approximation has been that the critical bond activation, that will irreversibly lead to dissociation or the reflection back into the gas-phase, often occurs at short interaction times, for which substrate mobility is not yet a critical issue. Obviously, the actual dynamics that proceeds after such a critical activation and closer to the surface is then described wrongly -the more, the heavier the adsorbate. But this has less consequences for the sticking probability as long as only the distinction "direct dissociative adsorption" or "immediate back-scattering" matters. The situation is less clear in case of molecular trapping that typically involves long interaction times. We return to this problem in section 2.3 below.

Adiabatic PES from DFT calculations
Within the frozen-surface approximation the remaining six-dimensional PES only depends on the degrees of freedom of the O 2 molecule, given by the cartesian coordinates R A = (X A , Y A , Z A ) and R B = (X B , Y B , Z B ) of the two constituent oxygen atoms, A and B, respectively. In the following, another very common equivalent coordinate system as depicted in figure 1 is made use of: Molecular configurations are conveniently described by the centre of mass of the molecule (X, Y, Z), the interatomic distance d, and the orientation relative to the surface, defined by the polar and azimuth angles ϑ and ϕ, respectively. The origin is placed on top of a silver atom within the surface plane, and (ϑ, ϕ) = (90 • , 0 • ) corresponds to the molecular axis oriented parallel to the [110] direction. To construct the 6D PES, we interpolate with the NN a large set of energies that have been first calculated with density-functional theory (DFT) for different positions and orientations of the molecule over the surface. Note that the role of subsurface oxygen [34,35] is not considered in this work because we are interested in the measured initial sticking coefficient, i.e. in the dissociation probabilities that are measured at very low coverages (about 0.1 ML or less). We have performed more than 7000 spin-polarized DFT total energy calculations, starting with so-called elbow plots, i.e. (Z, d) cuts through the 6D PES, at high symmetry sites of the surface. Additional irregularly spaced molecular configurations, typically in form of "scans" both in lateral (X, Y ) and angular (ϑ, ϕ) dimensions, were added iteratively as needed during the neural network interpolation (vide infra).
An energy cut-off of 400 eV in the plane-wave basis set together with ultra-soft pseudopotentials as contained in the default old Cambridge library of the CASTEP code [36] have been used to describe the interaction between electrons and nuclei, relying on the exchange-correlation functional due to Perdew, Burke and Ernzerhof (PBE) [37]. Based on these computational settings, the computed bulk Ag lattice constant and bulk modulus compare well to experimental values [38] ( a Ag 0,PBE = 4.14Å and B Ag 0,PBE = 97.1 GPa vs. a Ag 0,exp = 4.07Å and B Ag 0,exp = 109 GPa, respectively), and also the free O 2 bond-length is with d O2 0,PBE = 1.24Å well described (d O2 0,exp = 1.21,Å [18]). Notwithstanding, similar to other frozen-core calculations, the O 2 gas-phase binding energy is underestimated compared to accurate full-potential PBE results [39]. With 5.64 eV this fortuitously yields a value close to the experimental gas-phase data [40] (5.12 eV). For the binding energies addressed below this is not to be of concern, as this underestimation largely cancels in the total energy differences constituting the latter, and a potentially introduced constant shift does not affect the dynamics on the interpolated PES.
The Ag(111) surface is modelled with a periodic five-layer slab separated by 12Å of vacuum. In order to sufficiently suppress spurious interactions of the impinging O 2 molecule with its periodic images in the supercell geometry we are forced to use a (3 × 3) surface unit cell. The irreducible wedge of the first Brillouin zone corresponding to this large unit cell is then sampled with a (4 × 4 × 1) Monkhorst-Pack grid of special k-points. Systematic convergence tests focusing on the adsorption of atomic and molecular oxygen at high-symmetry sites confirm that the binding energies are numerically converged to within 50 meV at these settings. For the clean Ag(111) surface we reproduce the established marginal 1% outward relaxation of the first interlayer distance [41]. Within the frozen-surface approximation we neglect this and compute the O 2 /Ag(111) PES above a bulk-truncated slab. As discussed in detail by Behler et al. for O 2 on Al(111) [1,2], also in this adiabatic PES the spin of the O 2 molecule is gradually quenched from the triplet ground-state of O 2 at large Z distances to the spin unpolarized state close to the surface. We use the configuration where the O 2 molecule is located in the middle of the vacuum region as zero reference. This corresponds to the full decoupling limit between O 2 and Ag(111), since at these distances, unlike at Al(111) [1,2], there is no spurious charge transfer to the molecule. In the employed sign convention negative PES values indicate exothermicity with respect to this limit, while positive values indicate endothermicity.

Interpolation by symmetry-adapted neural networks
Neural networks represent a very flexible class of composed functions that in principle allows to approximate (physical) potential energy surfaces to arbitrary accuracy [42,43]. However, no physical properties of the latter are built therein a priori. The use of NNs in the context of gas-surface dynamics, pioneered by Lorenz et al. [28], has been described in detail before [29,30]. Therefore, we only give a brief account of the novel aspects of how symmetry is taken into account for the present system. Instead of using the "straightforward physical" coordinates depicted in figure 1, the neural network interpolation is performed on a properly transformed set, which has been termed symmetry functions before [30]. This way, two molecular configurations that are equivalent by symmetry are enforced to have same energy on the NN-interpolated PES by construction. By the rigid incorporation of these important physical properties of a specific system, the interpolation problem is thus simplified considerably (i.e. better quality with less input data). The construction of such a coordinate transformation is a rather complicated task, as it has to capture both periodic and point group symmetry on high-symmetry sites of the surface correctly. This implies certain demands for the behaviour under a certain (discrete but infinite) set of transformations including the molecular centre of mass position on the surface (X, Y ) and the molecular orientation (ϑ, ϕ), but not the distance from the surface Z and the internuclear distance d. Artificial symmetries have been introduced unintentionally in past attempts and only recently been found to have a significant effect on the calculated sticking probability [44]. For the present system, symmetries resulting from the combination of the threefold symmetry of the (111) surface of the closed-packed fcc metal silver and the exchange symmetry of the homonuclear O 2 molecule can be further increased by exploiting the near degeneracy of the O 2 interaction with the fcc and hcp hollow sites: For equivalent molecular configurations at these two sites, we observe a statistical relation in our DFT data set with a root mean square error (RMSE) of less than 5 meV and a mean absolute deviation of less than 3 meV. Similar results have been obtained before for H 2 on Pt(111) [45]. For the present work we will neglect this small difference and thus treat the fcc and hcp sites as equivalent. This increases the symmetry of the surface from three-to sixfold, as illustrated by the high-symmetry points and irreducible wedge in figure 2, together with the primitive lattice vectors a 1 and a 2 . Inspired by the symmetry functions obtained for a fcc(111) surface with threefold symmetry [30], the following set of symmetry-adapted coordinates Q i have been designed and carefully verified not to yield any artificial symmetries as "side effects". A detailed derivation and their simple generalization to other common lowindex surfaces of both close-packed fcc and hcp metals will be detailed elsewhere [46]: (2.1i) g 1 and g 2 are functions R 2 → R with the periodicity given by the surface lattice vectors a 1 and a 2 . They have been constructed such that points r = (x, y) in the surface plane which are equivalent by the sixfold symmetry depicted in figure 2 are mapped onto the same pair of function values (g 1 (x, y), g 2 (x, y)) [46]: a is the surface lattice constant given by 1/ √ 2 a Ag 0,PBE for the present system. The building blocks of g 1 and g 2 are low order terms of the corresponding Fourier series as indicated by the respective reciprocal lattice vectors Here b 1 and b 2 are the primitive vectors of the reciprocal lattice given by the defining relation and the crystallographic convention to indicate a negative direction by a bar over the corresponding number is adopted. Along the lines of previous work [30], neural network training has been simplified by limiting highly repulsive energies to ≤ 5 eV, which is still much larger than the energy range of interest for the ensuing molecular dynamics trajectories discussed below. Using the adaptive extended Kalman filter algorithm including the modifications described in [29] for training, the sensitivity of the NN is further focused on the most relevant parts of the PES by assigning training weights for E DFT ≥ 2 eV (2.5) to the input data. More than forty different NN topologies have been tried. In each case, the interpolation quality has been carefully monitored by plotting numerous twodimensional (d, Z), (X, Y ) and (ϑ, ϕ) cuts through the such obtained PESs, including information about the interpolation errors of individual data points. An example for a corresponding elbow is shown in figure 3 (a). With increasing knowledge gained about the PES and dynamics thereon, we decided for a {9-16-16-16-1 tttl} topology of the NN. In this commonly used notation [2,28,29,30], the first number refers to the nine symmetry-adapted coordinates (Q 1 , . . . , Q 9 ) in the input layer. The three numbers in the middle denote the number of nodes in the hidden layers, and the last number indicates that there is a single node yielding the PES value in the output layer. Like in previous work [2,29], a hyperbolic tangent and a linear function serve as activation functions in the hidden and output layer and are denoted by t and l, respectively. Figure 3 (b) summarizes the high quality of the finally obtained interpolation. For statistically significant subsets of the data points in the dynamically relevant entrance channel and chemisorption well (vide infra), the RMSE is smaller than 26 meV. Considering points from the training and test sets separately, RMSEs for the latter are larger by only a few meV -thus indicating good predictive properties for arbitrary molecular configurations to be encountered during the dynamics.

Molecular dynamics simulations
On the thus obtained continuous NN-representation of the adiabatic PES we perform trajectory calculations using a Bulirsh-Stoer integrator [47] for the classical equations of motion in Cartesian coordinates. Trajectories start with the O 2 molecule at its calculated gas-phase bond length and at Z = 7Å above the surface, neglecting its initial zero point energy. A conventional Monte Carlo procedure is used to sample all possible initial O 2 orientations (ϑ, ϕ) and lateral positions (X, Y ). For each incident energy E i , the sticking coefficient is obtained based on 10 7 different trajectories at perpendicular incidence. ‡ Such a large number of trajectories can only be integrated because the cost to compute forces from the NN PES-representation is negligible compared to the DFT calculations. As outcome of each individual trajectory we consider the following possibilities: (i) The trajectory is assumed to yield a dissociation event when the O 2 bond length d is stretched beyond 2.4Å, i.e. twice the gas-phase distance, and the associated velocityḋ is still positive. This is a rather conservative criterion and we have verified that some variation of the critical bond distance has no effect on the obtained sticking curves.
(ii) The molecule is assumed to be reflected when its centre of mass reaches the initial starting distance of 7Å with a positive Z-velocity. The results were equally found to be insensitive to reasonable variations of this criterion. (iii) Molecular trapping, when the molecule is neither dissociated nor reflected after 24 ps. A detailed analysis of the obtained trajectories shows that this time span is well separated from those of reflection or dissociation events, which happen in the order of 2-3 ps, such that the precise value of 24 ps has little influence on the number of determined trapping events. Although one may expect that the longer the molecule spends close to the surface, the more probable would be energy dissipation into either electronic or phononic excitations that in turn would prevent these trapped molecules to eventually escape from the surface, no final statement can be done in such cases within the frozen surface approximation (vide infra).

Insight from experiment
With silver as the almost unique industrial catalyst employed in the selective epoxidation of ethylene [48,49], the interaction of O 2 with silver surfaces has already been extensively investigated during the last decades in an attempt to understand the fundamentals behind such exceptional catalytic functionality (e.g. see [50] and references therein). The richness and complexity of the O 2 interaction with silver is reflected in a variety of adsorbed species that have been identified even at just the flat low-index surfaces Ag(100), Ag(110), and Ag(111), when varying the exposure conditions such as surface temperature, oxygen coverage, gas temperature and gas pressure [48,51,52,53,54]. In ultra-high vacuum physisorption states are observed at lowest temperatures, while molecularly chemisorbed states prevail up to T s ∼ 150 K. Dissociative adsorption appears for higher crystal temperatures; in this case subsurface oxygen can be also important for coverages where silver oxides are about to be formed [34]. All these studies showed that the abundance of each state also depends on the surface structure, with the close-packed Ag(111) surface being the most inert towards thermal oxygen adsorption [51]. All molecular beam experiments performed on the three low-index faces show that, if present, both the molecular and the dissociative initial adsorption probabilities ‡ Results of calculations for various non-perpendicular incidence angles focusing on scattering will be reported in a forthcoming publication.
are always activated [6,7,8,9,10]. Despite this consensus, there remain unresolved experimental discrepancies on the Ag(111) surface regarding the actual order of magnitude of the dissociation probability at normal incidence. Using a molecular beam, Raukema et al. reported initial dissociation probabilities that increase from 10 −5 to 10 −3 as the incidence energy E i increases in the range 0.8-1.7 eV [10]. In contrast, similar experiments performed with a molecular beam of 0.098 < E i < 0.8 eV, estimated that the dissociation probability on the clean Ag(111) surface should be lower than 9 × 10 −7 at least in that E i -range [14]. This has been related [14] to the particular propensity of the Ag(111) surface towards surface imperfections, notably steps, or towards residual gas contamination at lower surface temperatures.
In this situation the main purpose of our modelling is evidently not so much to reach full quantitative agreement with experiment, but rather to establish the order of magnitude for the initial dissociation probability within an adiabatic picture. As we will show below this is already sufficient to reach important conclusions concerning the role of possible non-adiabatic effects. In addition it contributes to the discussion of possible sources for the differences between the experimental data sets of [10,14] and will provide a valuable basis for future measurements that aim to overcome the present scatter in the experimental data.

Adiabatic PES properties
Overall, the obtained O 2 /Ag(111) PES is rather repulsive. A systematic analysis of (d, Z) cuts for different molecular configurations shows that at intermediate distances from the surface (Z ≈ 2Å) the potential energy is already above one eV in almost all cases. In this respect, the (d, Z) cuts of figure 3 (a) and figure 4 (a) are rather representative. As expected, the repulsion is strongest over the top sites. Already at distances Z > 3Å far from the surface the emerging corrugation gives rise to a very shallow well with a depth of the order of 20-30 meV. As the employed semi-local PBE functional does not account for van der Waals interactions, this is unlikely a proper representation of the experimentally reported physisorption state, but instead rather an artefact of the NN interpolation. As we will discuss below, this shallow (and likely spurious) well only plays a role for the O 2 /Ag(111) gas-surface dynamics at lowest incidence energies and its effects on the dissociation probability are easily separated.
Much more relevant for the dynamics is instead a second well that we determine at Z = 2.3Å, and for a molecular configuration where the O 2 centre of mass is situated above a bridge site and with the molecular axis pointing to the two closest Ag atoms. Figure 4 (a) shows the elbow plot corresponding to this configuration. The computed internuclear distance is with d = 1.28Å only slightly stretched compared to the gas-phase value (d = 1.24Å). This indicates a rather weak interaction that we find equally reflected in the computed low binding energy of −40 meV and the modest shift of the O 2 stretch vibration from the computed 192 meV gas-phase value down to 145 meV. These properties are only little affected by the frozen surface approximation and the NN interpolation. After a direct DFT calculation of this state that also allows for surface relaxation the most notable change is the increase of the binding energy to −70 meV. The latter value is somewhat lower than the −0.17 eV binding energy determined before by Xu et al. for the same O 2 molecular adsorption state [35]. Considering the highly similar computational setup employed in the latter study we attribute this 0.1 eV difference primarily to their use of a smaller (2 × 2) surface unit-cell. It is tempting to identify this molecular state with the chemisorption state reported in the experimental literature. For this, one has to recognize though that (regardless of the small scatter) both theoretical binding energies are at variance with the experimental estimate for the chemisorption state of ∼ 0.4 eV from thermal desorption spectroscopy (TDS) [51]. This difference is noteworthy as preceding equivalent DFT calculations performed for O 2 at Ag(100) [17] and Ag(110) [55,56] reproduce the approximately −0.4 eV chemisorption energies measured at these facets [48,57]. We also confirmed that the here employed computational setup fully reproduces the O 2 binding energy reported by Alducin et al. at Ag(100) [17]. This suggests that both the employed PW91 and PBE functionals, which are very closely related, are capable of describing the O 2 -Ag interaction rather well. Similarly problematic is the observation that none of the experimentally reported vibrational frequencies [58,59] comes close to the 145 meV that we compute for the molecular state. Instead much lower frequencies below 100 meV are measured, which, however, have already been assigned to hydroxyl, water groups or imperfections [58,59]. Furthermore, as shown by Xu et al. [35] also sub-surface oxygen can strongly stabilize chemisorbed O 2 at Ag(111). In this situation we believe that the low computed binding energy of the molecular state rather supports the interpretation that the chemisorption state prepared in the experiments by Raukema, Campbell and others [10,19,51,53] is not a property of the pristine Ag(111) surface, but reflects the propensity of this surface towards defects or residual gas contamination, and we will return to this point below when discussing the dissociative sticking probability. Within the theory-theory comparison of O 2 at the three low-index surfaces, the lower binding energy computed here at Ag(111) as compared to Ag(100) [17] and Ag(110) [55,56] is instead perfectly consistent with the expected chemical inertness of the close-packed surface.
As a final important characteristic of the adiabatic PES we also determined the minimum energy path (MEP) that leads out of the molecular well to dissociation. This path goes along the (d, Z) PES cut shown in figure 4 (b), in which the molecule is located over a bridge site with its axis oriented along the direction joining the closest fcc/hcp sites. The transition state corresponds to an activation barrier of 1.1 eV, which is in very good agreement to the value reported before by Xu et al. [35]. It is located at Z ≈ 1.6Å, i.e. at a distance closer to the surface (Z ∼ 2 − 3Å) than where one would expect barriers due to spin or charge non-adiabaticity (vide infra).

Dissociative sticking coefficient
We compile the results of our adiabatic sticking coefficient simulations for normal incidence in figure 5. Two different regions can clearly be discerned. At incidence energies below 0.4 eV there is an appreciable amount of molecular trapping events. This component rapidly decreases with increasing E i , as characteristic for physisorption. Note that the non-monotonous scatter of this curve around E i ∼ 0.2−0.4 eV is already no longer statistically relevant. For the employed 10 7 trajectories the corresponding probabilities around 10 −7 mean that we have observed either one or two molecular trapping events. Analysis of the trajectories in this entire low-energy trapping branch reveals that these events are all due to dynamical trapping in the shallow PES well located at distances Z > 3Å. As we are rather sceptical about the meaning of this shallow well we would not put too much emphasis on this low-energy results, even though the order of magnitude of this trapping and the incidence energy range fit rather well to those of the physisorption trapping discussed by Raukema and Kleyn [60,61]. We note however that a proper description of any kind of molecular adsorption requires to go beyond the frozen surface approximation by including surface temperature effects and energy dissipation channels in the simulations. In any case, this lowest-energy trapping part is well separated from the sticking properties at incidence energies higher than 0.4 eV, and for the latter part the shallow PES well plays no further role.
This other part that we will exclusively discuss henceforth starts at incidence energies above 1.1 eV and is characterized by a steeply increasing dissociative sticking probability. In contrast, in the intermediate range 0.4 eV< E i <1.1 eV we have not observed a single dissociation event in the total of 10 7 trajectories performed for every incidence energy. From our calculations we can thus put an upper bound to the real dissociative sticking probability at these incidence energies of 10 −7 . Above 1.1 eV this probability increases steeply and seems to saturate above 10 −2 at the highest incidence energies considered. Analysis of the trajectories shows that the dissociation follows a rather direct mechanism ruled by the large energy barriers existing at distances below 2Å from the surface. Almost all of the dissociating molecules follow paths very close to the MEP shown in figure 4(b) above. This explains the coincidence of the onset of sticking above 1.1 eV with the activation energy along this path and rationalizes the overall very low sticking probability at this surface in terms of the narrow dissociation channel and thus reduced configurational space.
If we compare these sticking results with the experimental data from Raukema et al. [10] that is also included in figure 5 we first notice a gratifying overall agreement in the sense that the computed adsorption probabilities are in the right (low) order of magnitude and above threshold show the correct monotonous increase up to the highest incidence energies considered. However, quantitatively, there are notable deviations. This concerns foremost a much steeper slope of the theoretical curve especially for E i < 1.3 eV. As a consequence, there is still appreciable sticking of the order 10 −5 − 10 −4 at 0.8-0.9 eV incidence energies in the experiments, whereas the theoretical curve drops steeply to below 10 −7 already at E i = 1.1 eV.
One reason that immediately comes to mind to generally rationalize such discrepancies is the approximate nature of the employed xc functional. In the present case, this is unlikely to explain all of the deviations though. As described before, almost all of the direct dissociation goes along paths close to the MEP shown in figure 4(b). If the xc functional messed up the MEP barrier height, this would thusto zeroth order -predominantly have resulted in a rigidly shifted theoretical sticking curve to too high or too low E i , and less to a wrong slope. If the employed PBE functional e.g. overestimated the true MEP barrier height, a "better" functional would shift the curve to lower E i . This would then lead to a better agreement at the onset region. However, simultaneously it would increase the difference to the experimental curve at the highest incidence energies, with the theoretical curve then seriously overshooting the experimental data there.
From this reasoning we thus rather suspect other factors to stand behind the disagreement between the theoretical and experimental sticking curve. The influence of vibrationally excited O 2 in the beam can be neglected, according to the low populations (less than 10 −3 for ν = 3 and than 10 −2 for ν = 2) that have been estimated by Raukema et al. based on the nozzle temperature [10]. Another factor that could be important in particular in light of the abrupt initial increase of the theoretical sticking curve is the finite energetic width of the experimental beam [62].
To mimic the effect of the latter in the theoretical data we recompute the curve simply as a convolution with an energy distribution function f β,E0 (E) for the beam: This equation is evaluated numerically, andS 0 (E ′ ) describes the corresponding curve plotted in figure 5 (i.e. a piecewise linear interpolation based on the explicitly calculated points). f β,E0 (E) is taken to be a shifted Maxwell-Boltzmann distribution Formulated in the velocity domain, it is known to properly describe supersonic molecular beams [63] and has also been used by Raukema for the evaluation of time-offlight (TOF) experiments for the present system [64].
is a normalization factor, and the width parameter β is chosen to reproduce the experimental full width at half maximum (FWHM) of 0.2 eV in the energy range of interest [60]. As seen in figure 5, this correction leads to a smoother theoretical curve that has its onset with 0.9 eV now at slightly lower incidence energies. Nevertheless, it does not help to reconcile experiment and theory especially at these E i close to threshold, where the experiment by Raukema et al. [10] yields a dissociative sticking that is more than two orders of magnitude higher than the computed one. Our results are instead more compatible to those of [14], which estimated a S 0 (E i ) below 9 × 10 −7 for E i = 0.8 eV.
Another factor that should be kept in mind is the finite surface temperature T s in the experiments that is not accounted for at all in the present theoretical approach. For the normal incidence data shown in figure 5 Raukema et al. used T s = 400 K [10]. For larger incidence angles they also reported data measured at a lower T s = 220 K. There, the lowering of the surface temperature led to an increase of the slope of the sticking curve above E i = 1.0 eV [10]. Unfortunately, no such lower temperature data was presented for normal incidence. If we assume that the trend seen at larger incidence angles prevails, we would expect normal incidence measurements at lower T s to yield a larger slope than the one of the experimental T s = 400 K curve shown in figure 5 -in closer agreement to the calculated sticking data.
In this respect it is also important to note that from the observed non-monotonous dependence of S 0 (E i ) on T s and on the incidence conditions (E i , incidence angle) Raukema et al. proposed the existence of two distinct dissociation mechanisms [10]: (i) a direct one that governs dissociation at normal incident energies E i > 1 eV and (ii) a precursor-mediated one at work in the threshold energy regime (E i < 1 eV), in which a molecular chemisorption state acts as the precursor.
In the absence of substrate mobility and concomitant phononic dissipation channels in the model, the closest we can get to a precursor-mediated mechanism within the present theoretical approach are the molecularly trapped events we indeed observe for E i > 1.1 eV. The probability for such trapping is rather energy independent and about 10 −6 − 10 −5 . Even if we thus assume each of these trajectories to eventually lead to dissociation, at this magnitude this contribution only significantly affects the total sticking coefficient at the direct onset around E i = 1.1 eV as shown in figure 5.
The resulting kinked sticking curve then has a shape that is very reminiscent of the experimental curve measured by Raukema et al. [10] for an off-normal incidence angle of 30 • and low surface temperature T s = 220 K. Unfortunately, again for normal incidence the reported data does not extend to similarly low incidence energies and surface temperatures to resolve this kink. Vice versa, we can at present only speculate if an account of substrate mobility in the calculations would for T s = 400 K modify the present frozen surface total sticking curve in figure 5 (dissociation plus trapping) to such an extent to reach quantitative agreement with the shown experimental curve. Particularly because of our scepticism that the here computed molecularly bound state on the adiabatic PES corresponds to the experimentally reported molecular chemisorption state, we do not believe that this will be the case. From the dependence on surface temperature Raukema et al. estimate a trapping probability of their precursor molecular chemisorption state of the order of 10 −4 at incidence energies as low as E i = 0.5 eV [10]. At such low energies our computed molecular state can not be accessed at all.
Our current interpretation is therefore instead that the discrepancy in the initial sticking coefficient below E i = 1 eV is primarily caused by another molecular precursor state that is absent in the theoretical calculations and that is not a property of pristine Ag(111). This state would then also correspond to the strongly bound state seen in the TDS experiments [48,51,52,53]. This interpretation that a molecular precursor state associated to surface imperfections is predominantly responsible for the dissociation probabilities in the intermediate energy range below E i = 1.0 eV has already been advocated by Rocca et al. [14] to rationalize the discrepancy in their measured sticking data and the one of Raukema et al. precisely in this energy range. Our analysis fully supports this view and places an upper limit of 10 −7 to the dissociative sticking at ideal Ag(111) for incidence energies up to at least E i = 0.9 eV. In this situation, new experiments are self-evidently required to finally settle this cause. Particularly for normal incidence they should be performed for a wide range of incidence energies (as in Raukema's work) and for as low surface temperatures as possible (as to not getting riddled by residual gas contamination).
Finally, we come back to the original motivation and the possibility that the discrepancy between experiment and theory could also be due to non-adiabatic effects. However, such effects that either involve sudden charge transfer [16] or spin flip hindrance [1] tend to release the system in an excited state and would thus in general further lower the dissociative sticking, which is already too low in its present adiabatic form. Due to the specific characteristics of the adiabatic PES not even such a lowering seems likely though. In contrast to e.g. the O 2 /Al(111) system [1,2], the adiabatic dissociative sticking at Ag(111) is governed by very late energy barriers larger than 1.1 eV and at Z-distances very close to the surface. Potential spin-or charge transferinduced hindrances instead occur earlier in the entrance channel at Z-distances around 2-3Å [1,2,16]. Even if such hindrances would give rise to additional barriers, those molecules that are able to surmount them would thus still be confronted with the adiabatic barriers when they approach closer to the surface and into the PES region where the spin quenching or charge transfer has then already taken place. A noticeable spin-selection rule induced lowering of the dissociative sticking as compared to our present adiabatic curve would therefore only occur, if the additional restrictions in the entrance channel were larger than those due to the late adiabatic barriers. This is not possible for energetic reasons, because the present adiabatic MEP barrier is with 1.1 eV already higher than the gas-phase O 2 singlet-triplet gap of 0.98 eV [18]. Similar arguments hold for possible non-adiabaticity effects due to charge transfer limitations, considering that the electron affinity of O 2 is 0.44 eV [18]. As such we reach the same conclusions with respect to non-adiabatic effects as for O 2 dissociation at Ag(100) [17]. Also there, the dissociative sticking was found to be governed by large and late barriers above 1 eV on the adiabatic PES. Correspondingly, the computed sticking coefficients at Ag(100) and at Ag(111) are strikingly similar as shown in figure 6. The higher values at Ag(100) for lower incidence energies are due to the higher accessible configurational space as illustrated in figure 6 by the elbow plots along the MEP. Intriguingly, at Ag(100) experiments do confirm the absence of dissociative sticking at 0.9 eV [9], which (considering the analogy of the theoretical findings at the two Ag surfaces) indirectly gives further support to our assessment that surface imperfections cause the larger dissociative sticking around of 10 −4 at this energy in the experiments by Raukema et al. on Ag(111).

Summary and conclusions
In summary, we have employed a DFT-based divide and conquer approach to study the reactivity of O 2 molecules impinging onto the Ag(111) surface. The determined adiabatic dissociative adsorption probability is extremely low and only exceeds 10 −7 for incidence energies above 1.1 eV. Our analysis directly relates this low reactivity to large energy barriers in excess of 1.1 eV located very close to the surface. This excludes that electronically non-adiabatic effects, either due to spin-selection rules [1,2] or due to charge transfer limitations [15,16], contribute to the low sticking coefficient. The late adiabatic barriers exceed both the gas-phase O 2 singlet-triplet gap and the O 2 electron affinity. Even if non-adiabatic effects would give rise to additional barriers in the entrance channel, those molecules that are able to surmount them would thus still be confronted with higher barriers when they approach the adiabatic PES region closer to the surface. As such, the sticking coefficient is insensitive to a possible nonadiabaticity in the O 2 -Ag(111) interaction and other quantities will have to be found as useful signatures for such effects.
The highly activated nature of O 2 dissociation at Ag(111) is in good agreement with the existing experimental data. Still, quantitatively, large deviations of more than two orders of magnitude are obtained at intermediate incidence energies around 0.8-0.9 eV, when taking the particular data set from Raukema et al. [10] as reference. Instead, our data is well consistent with the estimate of 9 × 10 −7 for this energy range by Rocca et al. [14], who had already suggested surface imperfections as reason for the higher dissociation probability in Raukema's data. The analysis of possible uncertainties in the present approach fully supports this point of view. At best, the high surface temperature employed in the experiments by Raukema et al. could be an alternative reason, and in this respect it would be intriguing to repeat the present adiabatic calculations with a model that includes some account of substrate mobility. However, more important to finally settle the cause would be new experiments at lower surface temperatures and specifically investigating the role of surface imperfections (like the step density) for the dissociative sticking of O 2 at Ag(111).