Electron-hole pairs during the adsorption dynamics of O2 on Pd(100) - Exciting or not?

During the exothermic adsorption of molecules at solid surfaces dissipation of the released energy occurs via the excitation of electronic and phononic degrees of freedom. For metallic substrates the role of the nonadiabatic electronic excitation channel has been controversially discussed, as the absence of a band gap could favour an easy coupling to a manifold of electronhole pairs of arbitrarily low energies. We analyse this situation for the highly exothermic showcase system of molecular oxygen dissociating at Pd(100), using time-dependent perturbation theory applied to first-principles electronic-structure calculations. For a range of different trajectories of impinging O2 molecules we compute largely varying electron-hole pair spectra, which underlines the necessity to consider the high-dimensionality of the surface dynamical process when assessing the total energy loss into this dissipation channel. Despite the high Pd density of states at the Fermi level, the concomitant non-adiabatic energy losses nevertheless never exceed about 5% of the available chemisorption energy. While this supports an electronically adiabatic description of the predominant heat dissipation into the phononic system, we critically discuss the non-adiabatic excitations in the context of the O2 spin transition during the dissociation process.


Introduction
Energy conversion at interfaces is at the centre of the rapidly growing field of basic energy science. This concerns desired conversions like solar to chemical energy, but also unavoidable by-products like the dissipation of chemical energy into heat. An atomistic understanding of the involved elementary processes is in all cases only just emerging, but is likely to question established views and macro-scale concepts. With respect to the dissipation of heat freed during exothermic surface chemical reactions such a prevailing view is that of a rapid equilibration with the local heat bath provided by the solid surface. The view is nurtured by the rare event dynamics resulting from the typically sizeably activated nature of surface chemical processes: While the actual elementary processes themselves take place on a picosecond time scale, times between such rare events are orders of magnitude longer. The understanding is then that in these long inter-process time spans any released chemical energy is rapidly distributed over sufficiently many surface phononic degrees of freedom to warrant a description in terms of a mere heat bath with defined local temperature. At the atomic scale this equilibration with the surface heat bath leads to an efficient loss of memory of the adsorbates about their history on the solid surface between subsequent rare events. This motivates the description in terms of a Markovian state dynamics that is e.g. underlying all present-day microkinetic formulations in heterogeneous catalysis [1,2,3]. In turn, the freed reaction energy only enters the determination of the local heat bath temperature, commonly achieved through a continuum heat balancing equation [4,5,6]. Figure 1. Schematic illustration of electron-hole pair excitation during the impingement of a gas-phase molecule on a metal surface. An electron is excited from an occupied state below the Fermi level ε F to an unoccupied state above, resulting in an excited state commonly referred to as electronhole (e-h) pair.
In many cases this prevalent framework seems to allow for an accurate account e.g. of catalytic conversions. However, particularly for nanostructured surfaces and highly exothermic reactions it is presently unclear whether this effective description is sufficient. With respect to the exothermicity this suspicion comes from recalling that freed enthalpies can well be of the order of several eVs. This is e.g. generally the case for the dissociative adsorption of molecular oxygen at catalytically relevant transition metal surfaces, which is the specific surface chemical reaction we want to focus on here. Several eVs is an enormous amount of energy on the scale of phononic degrees of freedom, which calls for efficient dissipation channels to achieve the assumed quasi-instantaneous local equilibration. At metallic substrates electronic excitations could hereby potentially represent an important additional energy sink. Unfortunately, no consensus has hitherto been reached concerning the role of this additional channel for the gas-surface dynamics and subsequent energy dissipation. Among others Tully and coworkers have frequently stressed that the lack of a band gap in metallic systems should in principle allow for electronic excitations, namely electron-hole (e-h) pairs, of arbitrarily low energies to easily couple with the nuclear motion of a particle impinging from the gas-phase [7,8,9,10,11,12,13,14]. This is depicted schematically in figure 1. On the other hand, recent comparisons to experimental data for hydrogen and nitrogen molecules on metallic surfaces indicate that electronically adiabatic descriptions provided within the Born-Oppenheimer approximation (BOA) seem to describe the initial interaction dynamics governing the dissociation extremely well [15,16,17]. As e.g. recently stressed by Juaristi et al [16] a proper high-dimensional description of the nuclear motion appears hereby significantly more important than the effect of e-h pair excitations.
For the particular case of oxygen dissociation the spin-triplet ground-state of gas-phase O 2 nevertheless brings in another aspect that could require to go beyond an adiabatic BOA treatment. With the spin of adsorbed O atoms at metal surfaces quenched, a spin transition needs to occur at some stage during adsorption process and could then well affect the gas-surface dynamics. The well-known paradigm system for this is O 2 at Al(111), where only an explicit account of a hindered triplet-singlet transition could reconcile first-principles dynamical simulations with the experimentally measured low sticking coefficient for thermal molecules [19,20,21,22]. This non-adiabatic hindrance was traced back to the inefficiency of both coupling mechanisms generally discussed to relax the spin selection rules that suppress reactions like oxygen dissociation [19,20,21,22,23]: The low mass number of Al leads to a small spin-orbit coupling and the low Al density of states (DOS) at the Fermi level prevents efficient spin quenching through the preferred tunnelling of minority spin electrons between substrate and adsorbate.
In this context the aim of the present study is to specifically address the electronic non-adiabaticity of the O 2 dissociation process at Pd(100). Apart from the obvious catalytic relevance, the specific motivation for the choice of the Pd substrate is its extraordinary high DOS at the Fermi level, which is among the highest known among the low-index transition metal surfaces [24,25,26,27]. Much in contrast to Al(111) this and the much higher mass number of Pd should favor an adiabatic spin transition, an expectation which we find confirmed by good agreement of a completely adiabatically computed initial sticking coefficient with existing experimental data. We are going to publish these results in detail elsewhere [28], thereby also addressing the most intriguing "different dissociation properties" of cartwheeling and helicoptering O 2 molecules on Pd(100). These have been observed only very recently [29,30], most strikingly exemplifying the present limits of the aforementioned atomistic understanding of an elementary process like oxygen dissociation. Here, however, the focus is on the excitation of e-h pairs instead, which according to the general arguments of Tully and coworkers [7,9,10,11,12,13,14] should be particularly facilitated for such a high Fermi-level DOS system. The central objective is thus to assess how much of the total amount of ∼ 2.6 eV chemisorption energy (vide infra) in this system is dissipated into this electronic channel. In addition the computed spin-resolved e-h pair spectra will be critically analysed to look for features that relate to the adiabatic spin transition.

Theory
Several approaches of varying accuracy have been suggested to compute electronic excitation upon surface adsorption. At the high accuracy end is certainly the "direct" first-principles simulation of e-h pair spectra within time-dependent density-functional theory (DFT) and an Ehrenfest dynamics for the nuclei [31,32,33]. However, for the intended semi-quantitative estimate of electronic dissipation, and in view of the high computational demands imposed by the O 2 at Pd(100) system, more approximate treatments are appealing. Notwithstanding, the rather rough description of the substrate electronic structure in the time-dependent Newns-Anderson model developed by Mizielinski et al [34,35,36,37] raises some concerns. The same holds for approaches based on electronic friction theory [8,38,39,40], as their results can depend sensitively on the way how friction coefficients are calculated [16,41,18,42,43,44,45]. More severely, already first applications within a forced oscillator model for the substrate electrons identify the proper description of spin transitions as an intrinsic shortcoming of electronic friction theory [46,47]. For the present purposes we therefore opt for an approach introduced recently by Timmer and Kratzer (TK) [48], which relies on perturbation theory applied to a time-dependent DFT framework. As the authors have demonstrated convincingly, it does not suffer from any troubles with the description of spin, and for the frequently studied reference system of hydrogen atoms at Al(111) very good agreement could be reached with the accurate timedependent DFT/Ehrenfest benchmark data [33]. At the same time the approach is computationally very efficient, as it only requires less demanding ground-state DFT calculations as input.
In the following we recapitulate the TK approach. This is done self-contained and in considerable detail to highlight several aspects that deserve specific attention for the aspired application to the O 2 at Pd(100) system. Among those is the spin transition, which is why also a comparison to electronic friction theory is included. Furthermore, several improvements of the original TK approach are described. This concerns a numerically more efficient implementation to make the application to computationally more demanding systems like the present one tractable. As previous work has mainly been focused on purely one-dimensional model trajectories (of monoatomic adsorbates), it also comprises the generalization to curved trajectories (of multi-atomic adsorbates).

Time-dependent perturbation theory
The basis for the following derivation is a classical trajectory Q(t) of a particle (a single atom or a molecule consisting of N atoms) from the gas phase impinging on the metal surface. In general, Q is a 3N dimensional vector of Cartesian coordinates of the atoms, which might be replaced by a suitable reaction coordinate description. It is assumed that the particle collides with the surface and is reflected, which is commonly referred to as a "full round trip". Consequently, for very small and very large times t the particle is assumed to be infinitely far away from the surface. This is meant in a physical sense such that the particle-surface interaction vanishes. In the philosophy of a perturbative approach, the effect of energy loss on the actual trajectory is neglected. Accordingly, a rigid surface that does not allow for any phonon excitation is considered, and energy dissipation into the electronic degrees of freedom of the surface is assumed to be sufficiently small not to significantly change the trajectory itself.
Thanks to the Runge-Gross theorem [49], the time dependent electronic manybody problem can be mapped onto an effective single particle Hamiltonian where t e is the kinetic energy contribution and the time-dependent effective potential is given by σ denotes collinear spin in two different channels within the framework of spin-DFT for collinear spin [50]. v σ ext describes the electron-ion interaction and explicitly depends on time due to the motion of the impinging particle along its trajectory Q(t). The dependence of the Hartree potential v σ H and exchange-correlation potential v σ xc on the time-dependent density ρ(r, t) is dropped in order to simplify readability.
At the "beginning" and the "end" of the full round trip, the (stationary) potential describes vanishing particle-surface interaction as mentioned before, yielding a (timeindependent) Hamiltonian So far, no approximations have been made. To obtain a scheme that is more computationally tractable than the direct time-dependent DFT approach used e.g. by Lindenblatt and coworkers [31,32,33], the time-dependent effective potential is now approximated by its counterparts in a "series of snapshots" of respective separate non-time-dependent ground state problems along the considered trajectory Q(t) as described by stationary DFT, , with the interaction part of the approximated effective potential v σ (t) vanishing at the boundaries of the considered trajectory Following the arguments of TK [48], the concomitant neglect of deviations from the instantaneous ground state charge density can only be justified for systems where the densities associated with excited electrons and holes are small compared to the overall charge density. From a theoretical point of view, this can, of course, be confirmed a posteriori if the calculated electron-hole pair spectra are sufficiently small (vide infra). Equations (2.4) and (2.5) define an approximation for (2.1), which is a quantum mechanical textbook example for a starting point of timedependent perturbation theory: The unperturbed Hamiltonian h (0) gets acted on by a (implicitly) time-(and spin-) dependent perturbation potential v σ (Q(t)). This potential might not be small, but should only be non-zero, cf. (2.6), during a short time interval, given by the time of "intense interaction" (collision) of the impinging particle with the surface. Up to first order, the transition rate p ij for an excitation within one spin channel σ from an occupied state |ε σ i into an unoccupied state |ε σ j (i.e. i = j), which belong to the part of h (0) that describes the (clean) surface, is then given by States associated with the impinging particle (spanning the other part of the Hilbert space which h (0) acts on) are not considered here, as the e-h pair excitations to be described are supposed to be located within the substrate. For closer contact with actual calculations, discrete sets of initial and final substrate states rather than continuous band structures are denoted by indices i and j, respectively, in the following. In practice, theses indices hence encapsulate k-point and band indices of the Kohn-Sham states of the clean surface.

Electron-hole pair spectra
The excitation spectrum for a complete round trip is obtained from (2.8) by integrating over time and summing over allowed transitions: Here ω are positive and non-zero excitation energies. Furthermore, here and in the following, for these kind of sums over transitions appropriate weight factors (depending on symmetry) for the k-point part of the indices i and j are implicitly considered. Integrating (2.9) by parts leads to (2.10) Here, the boundary term vanishes due to the properties of the interaction potential v σ as given by (2.6). Reinserting the result of (2.10) back into (2.9) gives where the "dressed" (transition) matrix elements are ultimately the key ingredients to be calculated as outlined in section 2.6 below.
Up to now, only occupations of substrate states corresponding to zero temperature have been considered in (2.8), (2.9) and (2.11). Generalization to finite electronic temperatures is straightforward by introducing Fermi occupation factors Not allowing for de-excitations of "thermally smeared" electrons, (2.11) then gets extended to the following final working expression for the electron-hole pair spectrum: (2.14) Individual spectra for electrons and holes can be obtained consistently by considering respective transitions only relative to the Fermi level, where, as before, the excitation energies ω are positive and non-zero. As pointed out by TK [48] (2.14), (2.15a) and (2.15b) stay mathematically well defined even for |ε j − ε i | → 0. Nevertheless, numerical evaluations have difficulties in these regions, which is why they are omitted in plots of calculated spectra below.

Comparison to electronic friction theory
When electronic friction theory [38,39,40] is combined with the forced oscillator model [46,47] and introducing further approximations to obtain excitation spectra [51,52], the equations directly comparable to (2.15a) and (2.15b) read as follows when adapted to the present notation: based on the electronic friction coefficient For the sake of simplicity and better readability, but without loss of generality, a one-dimensional description via a reaction coordinate Q(t) has been employed here. Comparing (2.12) and (2.17) the key ingredient to be calculated for both theories is very similar: Matrix elements of the type with v σ eff (Q(t)) as defined in (2.5). Obviously, it does not matter for the derivatives in (2.12) and (2.18) that the constant part v σ eff; non−int , cf. (2.5), is not subtracted here. In fact, due to the energy differences entering λ σ fric in (2.16a) and (2.16b), this can also be read as matrix elements for electronic transitions -similar to λ σ ij . Their dependence on the electronic structure is less direct though, i.e. "only" due to the forced oscillator model connected to the electronic friction coefficient η σ . Additionally, it is important to note, that instantaneous Eigenstates |ε σ i (t) , |ε σ j (t) are used to obtain the latter. Obviously, the present theory can in principle be modified in a straightforward way to use instantaneous Eigenstates as well in (2.12). However, this leads to conceptual issues with shifts of Eigenvalues in the DFT calculations at different points of a trajectory and technical problems with (relative) phase shifts of the instantaneous Eigenstates.
Apart from this, the main difference between the present and electronic friction theory comes from the fact that squared moduli of the (respective) matrix elements appear before the time integration is carried out in (2.17) and (2.18). As discussed extensively by TK [48], this is the reason why (2.12), (2.14), (2.15a) and (2.15b) remain mathematically well defined even in case of an adiabatic spin transition: The spin population assigned to the impinging adsorbate typically (at least for hydrogen atoms impinging on different substrates) [31,32,33,46,47,48] is proportional to singularity of dv σ eff (Q(t))/dQ, for which the aforementioned main difference is of crucial (mathematical) importance.

Dissipated energy
Finally, the energy dissipated into electron-hole pair excitations, which is the key objective here, is obtained as energy-weighted integral over the excitation spectrum, (2.14): Again, a comparable expression is also obtained within electronic friction theory, which additionally would allow to "measure" the dissipated energy for only a part of the considered trajectory Q(t). However, with friction theory not being (directly) applicable here, (2.19) should serve well for providing the desired estimate of energy dissipated into electron-hole pairs.

Realistic first-principles trajectories
In previous work addressing e-h pair excitation typically strictly one-dimensional model trajectories Q(t) along the surface normal have been considered when studying adsorption of mono-atomic adsorbates with perpendicular impingement over high symmetry sites [31,32,33,46,47,48,51,52]. The motion of incoming poly-atomic particles (like the diatomic molecule O 2 ) is in contrast less likely to be reducible to a single spatial dimension. To keep the input data entirely scalar, an effectively one-dimensional description of coordinates and velocities, i.e. a suitable reaction coordinate Q, must therefore be employed. A canonical choice is the arc length associated with a given trajectory Q(t): However, we note that in case of non-negligible changes of relative coordinates within the adsorbate, (2.21) does not provide immediately obvious physically meaningful information about the trajectory of an adsorbate. For a stiff diatomic molecule like O 2 this might not be much of a concern, unless considering large vibrational excitations e.g. during scattering at high kinetic energies. Still, it is interesting to note that (2.12) and all its descendants remain invariant, if trajectories are described by scaled and shifted coordinates and consistent velocities (2.22b) When α = 1/ √ N ads for an adsorbate consisting of N ads atoms, the corresponding arc lengthQ(t) then conveniently corresponds to the absolute distance travelled by the centre of mass of the adsorbate (or likewise any of its constituents).
The trajectories employed in the present work have been obtained from molecular dynamics simulations on a six-dimensional (Born-Oppenheimer) potential energy surface (PES) V 6D , describing the molecular degrees of freedom of an isolated O 2 molecule on a frozen Pd(100) surface. This continuous PES representation is constructed by a highly accurate interpolation of discrete DFT energies computed within the Perdew, Burke and Ernzerhof (PBE) generalized gradient approximation to electronic exchange and correlation [53,54]. Further details and a statistical analysis of the resulting dissociation dynamics will be published elsewhere [28]. Briefly, the interpolation is based on neural networks, adapted to properly take symmetry into account [55]. 4000 single-point DFT calculations in supercell geometries, carefully sampling molecular positions and orientations over the Pd(100) substrate, serve as database for the interpolation. The supercell geometry is characterized by (3 × 3) supercells, a five layer slab, and a vacuum distance of 15Å. A (4 × 4 × 1) Monkhorst-Pack grid is used for k-point sampling. Electronic states are described by a plane wave basis using a cut-off energy of 400 eV together with ultrasoft pseudopotentials as implemented in the Castep code and supplied within the default Castep pseudopotential library [56].

Efficient matrix element calculation
The electronic structure data obtained concomitantly with the single-point energetics for the PES then also forms the basis for the evaluation of the matrix elements entering (2.12). In several applications of electronic friction theory, the corresponding matrix elements given by (2.17) are related to the nuclear kinetic energy terms, which are neglected within the Born-Oppenheimer approximation. Following ideas from Head-Gordon and Tully [7,8], the commonly used plane wave implementation by Lorente and Persson [57,58,59] makes use of first order changes in Kohn-Sham states and Eigenvalues to obtain the aforementioned matrix elements. Obviously, this cannot work in the present context when dealing with unperturbed states of the substrate.
In the original TK work [48] and in an equivalent strategy previously used within electronic friction theory [42], the matrix elements are calculated directly through the derivative of the effective Kohn-Sham potential. Quite in contrast, here we follow a different strategy instead, namely to compute at all points Q(t n ) for which electronic structure data is available. Afterwards, each of these matrix elements is interpolated along Q and the analytical derivative of the interpolation gives ∂M σ ij (Q)/∂Q. Through this then also yields the matrix elements required to determine the e-h spectrum via (2.12). Following common practise to rely on sufficient decoupling between periodic images [31,32,33,48,57,58,59], only intra-k transitions are considered in this procedure, which becomes exact for a periodic overlayer [42].
The crucial advantage of this new strategy is that it allows for efficient use of existing infrastructure in most plane wave codes, because the derivative of the potential does not need to be constructed in a form which can be used to act on the states directly. This is particularly useful (if not essential) when dealing with non-local and/or ultrasoft pseudopotentials. Notwithstanding, the matrix elements calculated for all Q(t n ) have to be kept in memory to perform the interpolation and the ensuing derivative. This does lead to higher memory demands than for other schemes, which by virtue of (memory) parallelisation nevertheless does no longer constitute an invincible problem. We have carefully validated and benchmarked our implementation strategy by recomputing the H at Al(111) data reported in the original TK work [48]. Quantitative agreement is reached at a minute fraction of the computational cost. Similar observations have also been made by Timmer and Kratzer in their more recent work [51] after also adapting to our implementation strategy.
In theory, the loss of orthonormality between substrate states when using ultrasoft pseudopotentials [60,61] could pose a serious problem for (2.8), hence shaking the foundations of the perturbative approach. However, in practice, the aforementioned tests for H at Al(111) have shown that this is not a general concern. Selected comparisons of the results for the present system to results obtained with normconserving pseudopotentials indicate some quantitative dependence of the ultimately deduced dissipated energy, yet, without touching on any of the physical conclusions discussed below. Much more problematic is the well-known spurious ferromagnetism of Pd within semi-local DFT [62,63]. Test calculations indicate that this leads to a dramatic increase of the dissipated energy that can be of the order of several 100 % in relative terms. For all results presented below, matrix elements have therefore been evaluated using the Kohn-Sham states of a non-spin polarized clean Pd(100) surface that have been "cloned" into both spin channels of the (spin-polarized) adsorption calculations including the O 2 molecule. In this respect the fact that only the derivative of the potential enters (2.12) can be seen as particular (additional) virtue of the employed scheme for this specific system: It allows to exploit a cancellation of errors in differences (as already in case of the energetics) -quite in contrast to a conventional treatment e.g. within time-dependent DFT. and hole (at negative excitation energies ω) spectra P σ ex,el ( ω) and P σ ex,ho ( ω) according to (2.15b) and (2.15a). (d) Total e-h pair spectrum P σ ex ( ω) given by (2.14) together with resulting dissipated energies according to (2.19). All spectra are for a half round trip with energies ω relative to the Fermi energy. Both majority (↑, violet) and minority (↓, blue) spin channels are shown.

Electron-Hole Pair Spectra and Dissipated Energies
With the objective to compute e-h pair spectra and to estimate the concomitant energy dissipated into electronic excitations we focus our investigation on four distinct trajectories, selected to span the space of possible O 2 impingement over the Pd(100) surface. Namely these are two trajectories in which the O 2 molecule initially approaches with its molecular axis parallel to the surface and oriented along a [100] direction (a so-called side-on configuration), and two trajectories in which the O 2 molecule initially approaches with its molecular axis perpendicular to the surface (so-called head-on configuration). In order to sample the influence of the lateral corrugation of the surface potential one of the two side-on trajectories starts out with the O 2 center of mass located above a fourfold hollow site, and one trajectory with this centre above a twofold bridge site, cf. insets in figures 2 and 3. Similarly the two head-on trajectories start out above a hollow and above a top site, cf. insets     in figures 4 and 5. Henceforth we will use the intuitive abbreviations h-para, b-para, h-perp and t-perp for the four cases, respectively. Each trajectory begins with the O 2 center of mass at a perpendicular distance Z 0 = 6Å above the surface, where the adsorbate-substrate interaction potential, i.e. the PES V 6D , has already well decayed to zero. As an initial kinetic energy we choose 50 meV to mimic thermal molecules in the high-energy Boltzmann tail. In general, non-adiabatic energy losses should increase with higher velocities, which is why we expect the deduced dissipated energy to represent a useful upper bound for a thermal gas. However, in order to explicitly test the dependence on the initial velocity quantitatively we additionally consider a h-para trajectory with an increased initial kinetic energy of 400 meV, corresponding to an almost tripled initial velocity fromQ(0) = 5.5Å/fs toQ(0) = 15Å/fs. The results for the four trajectories h-para, b-para, h-perp and t-perp with low initial kinetic energy are compiled in figures 2 to 5, respectively. In all four cases the molecule is reflected, such that a half round trip as assumed in the perturbative approach can be consistently defined. With the definition given in (2.22a) the reaction coordinate Q(t) corresponds for the considered trajectories essentially to the vertical distance Z(t) of the O 2 center of mass above the surface. As shown in parts (b) of figures 2 -5 the half round trip proceeds in all cases in about 600 fs, in which the molecule approaches from its initial height of Q(0) ≈ Z(0) = 6.0Å towards the surface. Intuitive in terms of the expected corrugation of the surface potential, the point of reflection (defining the end of the considered half round trip) is with Z ≈ 2.6Å highest for the t-perp trajectory with the molecule head-on above the top site, and with Z ≈ 1.0Å lowest for the h-para trajectory with the molecule side-on above the hollow site. In the other two trajectories the reflection occurs at intermediate Z heights.
These differences in the point of closest encounter have direct consequences on the adiabatic spin transition of the O 2 molecule. As illustrated by the projected O 2 spin density in parts (a) of figures 2 -5 the molecule approaches the surface only sufficiently closely in the h-para trajectory for the spin to be completely quenched at the point of reflection. In this case exactly the square-root like decay of the spin density results that has been generally observed for hydrogen adsorption before [31,32,33,46,47,48] and which as described in section 2.3 would prevent application of electronic friction theory. In contrast, in the t-perp trajectory the transition to the singlet state has not even really started when the molecule gets reflected. However, the adiabatic spin transition of the oxygen molecule depends not only on the penetration depth, but also strongly on lateral position and orientation. This is nicely illustrated by the trajectories b-para and h-perp in figures 3(a) and 4(a), respectively. In both cases the minimum Z-height reached is comparable, but the projected spin density on both oxygen atoms remains nearly twice as large for b-para.
Significant acceleration beyond the initial velocity ofQ(0) = 5.5Å/fs only takes place in areas of attractive potential at low Z heights just before the molecule hits the repulsive Pauli wall. As shown in figure 2(b) this is obviously largest for the most favourable h-para configuration, which is very close to the main entrance channel for dissociation in this system [28]. While the maximum speed reached is in this case about five times the initial velocity, essentially no acceleration occurs throughout the entire trajectory in the most repulsive t-perp case, cf. figure 5(b). These differences have to zeroth order exactly the bearings on the generated e-h pairs and concomitant energy dissipation one would expect e.g. from electronic friction theory. As shown both separately for electrons and holes in figures 2(c) -5(c) and combined into one e-h pair spectrum in figures 2(d) -5(d), the qualitative shape of the spectra is for all trajectories about the same. However, the absolute magnitude is much larger for the high velocity attractive case h-para, resulting into a total amount of dissipated energy that is about an order of magnitude larger than in the other three cases. Notwithstanding, the latter three cases illustrate that at closer look also other factors like the orientation of the molecule matter: Specifically in the trajectories b-para and h-perp in figures 3(a) and 4(a), respectively, the O 2 molecule experiences the mentioned similarly repulsive interaction potential and reaches comparable maximum speed before reflection. Nevertheless, the energy dissipation into e-h pairs is by a factor of four less for the latter. This constitutes a nice confirmation of one of the key results of the electron friction work of Juaristi et al [16]: The importance of the "high" dimensionality of the molecule-substrate-interaction also extends to e-h pair excitations, such that a proper assessment of the role of this dissipation channel needs necessarily to rely on a representative set of impingement scenarios, not just one model trajectory.
Further insight into the velocity dependence comes from the final h-para trajectory, in which the initial kinetic energy has been raised to 400 meV. Despite the nearly tripled initial velocity, the corresponding trajectory and spectra (not shown) do not change qualitatively with respect to the data shown in figure 2. The total energy dissipated into e-h pairs is increased by around a factor of 1.3, resulting in a total amount of about 100 meV. Compared to the electronic friction theory results for N 2 on W(110) [16], the relative increase of the total dissipated energy is thus much less. There, for perpendicular incidence, a similar increase of the incidence kinetic energy has led on average to an increase in the energy dissipation of more than an order of magnitude. At present it is unclear whether these differences are due to the two different systems studied, or whether they indicate potential shortcomings of either theory. Direct application of the present theory to the N 2 on W(110) system will therefore constitute an interesting next step in our work.
Among the trajectories considered, the total dissipated energy per half round trip is thus largest for exactly the h-para trajectory that is "most interesting" for thermal impinging molecules in the sense that it falls very close to the dominant entrance channel to dissociation into which almost all molecules get steered [28]. For this hpara trajectory the total dissipated energy furthermore seems not very sensitive to the specific initial velocity. Considering also the approximate nature of the employed perturbative approach we therefore take ∼ 100 meV as a reasonable estimate of the order of magnitude of the electronic dissipation channel. With the present DFT-PBE setup we compute a total adsorption energy released in the adsorption of O 2 , i.e. in the exothermic dissociation into two separate O atoms at the Pd(100) surface, of 2.6 eV. As such the e-h pair excitation channel does not even take up 5 % of this energy, which provides an a posteriori justification for the assumption of a weak perturbation underlying the present approach. On the one hand, the insignificance of this channel for this system is somewhat consistent with the inability to detect chemicurrents during the adsorption of O 2 on polycrystalline palladium [64]. On the other hand, it is remarkably low compared to H impingement over the threefold hollow site of Al(111). The latter leads to a comparable release of chemisorption energy, but the computed amount of energy taken up by e-h pairs was one order of magnitude more, i.e. 1 eV (albeit likely favoured by a penetration of the adsorbate into the first substrate layer) [32]. Neither is it thus possible to a priori dismiss electronic excitations in the adsorption process at metal surfaces, nor does the difference between the two systems follow the picture expected from the much higher Fermi-level DOS in the Pd(100) system.

Spin transition
At a closer look at the e-h pair spectrum for the h-para trajectory shown in figure 2(d) a slight asymmetry between electrons and holes deserves further attention. The fact that excitations are stronger in the spin (down) minority channel, and specifically in the high-energy wing of the spectrum, is particularly interesting in light of the spin transition of the O 2 molecule, which, according to figure 2(a), is completed for this trajectory before the point of reflection is reached. A similar asymmetry, albeit much weaker, also seems to be present for the h-perp trajectory, cf. figure 4, correlating with the aforementioned similar degree of spin quenching in this case. Timmer and Kratzer have obtained nothing comparable in the e-h pair spectra of H on Al(111) [48]. In contrast, for the same system Lindenblatt and Pehlke [31] did determine a similar abundance of excited electrons and holes in the spin majority and minority channel, respectively. However, since the time-dependent DFT based dynamics is intrinsically non-adiabatic and the spectra are obtained a posteriori by relaxation back to the electronic ground state, this is not surprising. It can be easily understood by propagation on excited PESs during the non-instantaneous spin transition. In contrast, the present description is adiabatic with respect to the spin transition, and e-h pair excitations are evaluated based on unperturbed substrate states.
In order to put the observed spin asymmetry into perspective with the spin transition, a detailed analysis of the O 2 molecular states interacting with the Pd(100) surface along the trajectory Q(t) can be illustrative. For this we evaluate the density of states projected onto a Kohn-Sham state |ϕ σ O2 (Q) of spin σ of the free O 2 molecule (molecular PDOS) [65] It is important to emphasize that here the Kohn-Sham states of the interacting system of substrate and molecule |ε σ n (Q) are used for this analysis. The most relevant state with respect to overlap with the Pd valence d-band is specifically the degenerate 2π * state [19], which in the free triplet O 2 molecule is occupied by two electrons in the spin majority channel and unoccupied in the spin minority channel. Upon approaching the surface, the interaction with the substrate potential starts to broaden the discrete molecular levels. Their centre of energy along the trajectory Q(t) is each given by As shown in figure 6, at the beginning of the calculated trajectory this centre is for the majority 2π * ↑ level located about 2 eV below the Fermi-level, while it is at about 0. above the Fermi-level for the minority 2π * ↓ level. If from inspection of figure 2(d) we roughly take the dominant fraction of e-h pair excitations to spread over a range of 0.4 eV around the Fermi-level, we see that it is primarily the spin minority 2π * ↓ level that falls close to this range. To quantify this as the width of the 2π * level increases along the trajectory, we evaluate the number of projected states that lie within this relevant energy range as Altogether we therefore conclude that the observed stronger e-h pair excitations in the minority spin channel coincide with a larger amount of substrate states in the corresponding energy range, which overlap with the 2π * ↓ orbital of the O 2 molecule. A tempting interpretation is therefore to see the excess excitations as a reflection of the preferred tunneling of excited minority spin electrons that leads to the quenching of the oxygen spin-density in the here assumed adiabatic spin transition picture. If the electronic structure during the "real" spin transition does not behave in a radically different fashion, this connection between the different energetic position of majority and minority O 2 2π * level and corresponding different overlap with the energy range spanned by e-h pair excitations would prevail. In that case, this then provides an intuitive mechanism how electron tunneling from the (here forcibly) non-spin polarized substrate electronic structure compensates the adsorbate spin.

Conclusions
In conclusion we have advanced and employed a first-principles perturbative approach to study the energy dissipation into e-h pair excitations during the adsorption of O 2 at Pd(100). Despite the unusually large Pd density of states at the Fermilevel, which in principle could render this system a prime candidate for electronic non-adiabaticity, we obtain energy losses into this electronic channel that are at maximum of the order of 5 % of the total released adsorption energy. In detail, this loss and the underlying electron-hole pair spectra differ considerably for different O 2 impingement configurations, i.e. they depend on the molecular orientation with respect to the surface and on the lateral position of impact. By far the largest amount of excitations is obtained for a side-on approach above the fourfold hollow site. This is also the statistically most relevant configuration as it corresponds closely to the entrance channel for dissociation into which the dominant fraction of impinging molecules gets steered. While these variations confirm the necessity to account for the "high" dimensionality of the molecule-substrate-interaction by considering appropriate statistical averages over ensembles of trajectories, we find only a weak dependence of the total dissipated energy into e-h pairs on the initial O 2 velocity.
All in all we therefore derive from our calculations a value of ∼ 100 meV as a reasonable estimate of the order of magnitude for the electronic dissipation channel in this system. Even when cautiously considering that the approximations underlying the employed perturbative approach lead to a similar underestimation of this quantity as was identified for H at Al(111) [48], e-h pair excitation thus still plays an insignificant role for the dissipation of the large amount of 2.6 eV adsorption energy released during dissociative adsorption. This assessment is very much in line with the trends observed in other studies beyond single atoms at metal surfaces, no matter if the impinging diatomic molecule carries a permanent dipole moment (HCl on Al(111) [66]) or not (H 2 on Cu(110), N 2 on W(110) [16]), as well as with the hitherto unsuccessful attempts to detect chemicurrents in experiments over polycrystalline palladium [64].
Notwithstanding, two aspects underlying the present state-of-the-art treatment need to be emphasized that could in principle radically question our conclusions on the role of e-h pair excitations. First, only non-reactive trajectories at a frozen surface are considered. While one trajectory is very close to the actual entrance channel for dissociation, this still can not capture possibly important non-adiabatic effects during the actual dissociation at a mobile substrate. Second, an adiabatic spin-transition from the triplet ground-state of the free O 2 molecule to the adsorbed singlet state is assumed. While a definite answer to the first point has yet to be found, there are good arguments that support the latter assumption. The high Pd density of states at the Fermi-level as well as the high Pd mass number favour both coupling mechanisms generally discussed to relax the spin selection rules that suppress reactions like oxygen dissociation, namely spin-orbit coupling and electron tunnelling between substrate and adsorbate. With respect to the latter it is noteworthy that an intriguing asymmetry of electron-hole pair excitations in the majority and minority spin channel is obtained in the present framework. The excess minority spin excitations span thereby an energy range around the Fermi-level in which there is a larger overlap of Pd substrate states with the minority 2π * ↓ orbital of the O 2 molecule. In the here implied adiabatic spin transition through efficient electron tunneling, this nicely illustrates the mechanism with which the non-spin polarized substrate electronic structure quenches the adsorbate spin-density.