We study the electronic coupling between an adsorbate and a metal surface by calculating tunneling matrix elements Had directly from first principles. For this, we employ a projection of the Kohn–Sham Hamiltonian upon a diabatic basis using a version of the popular projection-operator diabatization approach. An appropriate integration of couplings over the Brillouin zone allows the first calculation of a size-convergent Newns–Anderson chemisorption function, a coupling-weighted density of states measuring the line broadening of an adsorbate frontier state upon adsorption. This broadening corresponds to the experimentally observed lifetime of an electron in the state, which we confirm for core-excited Ar*(2p3/214s) atoms on a number of transition metal (TM) surfaces. Yet, beyond just lifetimes, the chemisorption function is highly interpretable and encodes rich information on orbital phase interactions on the surface. The model thus captures and elucidates key aspects of the electron transfer process. Finally, a decomposition into angular momentum components reveals the hitherto unresolved role of the hybridized d-character of the TM surface in the resonant electron transfer and elucidates the coupling of the adsorbate to the surface bands over the entire energy scale.

The electronic coupling Had between an adsorbate and a surface is a central quantity in many technologically important physical processes, including photochemistry, photovoltaics,1–3 and heterogeneous catalysis,4–7 not to mention its place in the very theory of the chemical bond.8,9 The coupling of an adsorbate to the d-band of transition metal surfaces has, for example, famously served as a descriptor to understand catalytic reactivity trends for a wide array of important chemical reactions.4–6,10–15 Electronic coupling likewise plays an important role in the interpretation of experimental studies of charge transfer (CT) between adsorbates and surfaces, where it is invoked to understand, e.g., non-adiabatic excited-state dynamics,16–18 desorption induced by (multiple) electronic transitions (DIET/DIMET),19,20 scanning-tunneling microscopy (STM),21–26 or pump–probe and core–hole clock spectroscopies.

Given the prevalence of electronic couplings in theoretical descriptions, there has been a high interest in determining them with computational methods,1,4,27 including, but not limited to, ab initio methods.4,28 Coupling matrix elements Had themselves are not physical observables though, and they thus evade direct comparison to experiment.29 Furthermore, their definition in the theory is in general non-unique, complicating efforts to calculate them and compare alternative approaches.2,30,31 A careful interplay between theory and experiment is therefore necessary for a meaningful quantitative determination of these coupling parameters, which are most often approximated4,32,33 and used only qualitatively.

Of the many experimental setups where electronic couplings play a pivotal role, CT experiments provide a most natural context in which to test and validate a quantitative theory of electronic coupling. For example, experiments based on core-hole clock spectroscopy have measured the ultrafast electron transfer (ET) dynamics of simple core-excited adsorbates, such as Ar*(2p3/214s), on transition metal surfaces to high precision.34–45 As illustrated in Fig. 1(a), the probe atoms are thereby physisorbed to ultra-cold surfaces and optically excited. The thus resulting excited donor state—labeled by the index d—then decays on femtosecond time scales into the unoccupied acceptor states—labeled with an a—of the surface. A schematic of the energetics of such a process exemplified here by adsorbed argon is given in Fig. 1(b). Interpretation of such experiments usually proceeds through theoretical models of various levels of complexity.46,47 From a theoretical perspective, systems like this provide an ideal test-bed for electronic coupling models. First, with their ET timescales in the regime of a few femtoseconds, nuclear dynamics plays little to no role in the mechanism. Together with the fact that Ar only physisorbs at such surfaces, this greatly reduces the number of possible influences an ET theory needs to consider.28 Yet, even such simple systems can show quite a complex ET behavior, strongly influenced by the type and density of the surface acceptor states and their coupling to the transferred electron.35,38 With the nowadays achievable experimental resolutions at the attosecond scale,39 they thus offer fundamental insights into the ET process and provide a valuable test case for theories of electron transfer.

FIG. 1.

(a) Diagram of the core-excited electron transfer system of an Ar atom adsorbed on a transition metal surface. Ar is thereby optically core-excited into a 4s donor state ψd, which can then decay into the spectrum of surface acceptor states {ψa}. (b) Schematic energy diagram of the process.

FIG. 1.

(a) Diagram of the core-excited electron transfer system of an Ar atom adsorbed on a transition metal surface. Ar is thereby optically core-excited into a 4s donor state ψd, which can then decay into the spectrum of surface acceptor states {ψa}. (b) Schematic energy diagram of the process.

Close modal

Early experiments38,48 and subsequent theoretical work46,47,49–51 have explained measured (and simulated) ET lifetimes primarily in terms of the surface electronic band structure, invoking the projected bandgap that is present in the dispersive sp-bands on many transition metal surfaces. These models thus centered solely on a coupling to the sp-type states, which is easily approximated using expressions common in STM simulations and which references their free-electron-like character.21–23,26 Even though neglecting the possible role of the d-band, these models have succeeded in explaining aspects of the ET process, including its spin-dependence on magnetic surfaces51 and the observed energy dependence of the ET lifetimes.47 The latter dependence results simply from the relative position of the resonance within the projected sp gap, which reduces the density of available acceptor states and hinders (resonant) ET with high k momentum (perpendicular to the surface). This simple sp-band model of ET is generally well motivated for Ar adsorbates on transition metal surfaces. Argon atoms physisorb more than 3 Å above the surface, and their states should thus overlap more strongly with diffuse and free-electron-like sp-states than with more localized d-states. Furthermore, the resonance energy lies some 3 eV above the substrate’s Fermi energy, and thus well above the typical metal d-band. Nevertheless, the limitations due to a complete neglect of possible d-band coupling were pointed out already by Menzel in 200038,48 and Gauyacq and Borisov in 2004.50 In these studies, a model based on sp-band arguments alone and neglecting a possible d-character of states was thought to be inadequate for a comprehensive understanding of differences in ET lifetimes observed on different surfaces. First-principles simulations were therefore called for explicitly to resolve the coupling to the d-type states, which is more difficult to describe.38,48 While recent ab initio methods, based on density-functional theory (DFT) combined with surface Green’s function techniques, have succeeded in accurate predictions of the ET lifetimes on a number of surfaces,39,47,51 the limited interpretability of these methods—which proceeds by the same sp-band model as above—has left the role of coupling to the d-band in ultrafast ET unresolved.

In this work, our objective is twofold. First, we provide a first-principles model of ultrafast ET on metal surfaces, which is based on the direct calculation of electronic couplings Had between the adsorbate and all surface bands.

For this, we employ a projection of the Kohn–Sham electronic structure upon a basis of suitably defined localized states or diabats. While this procedure is referred to in the charge transfer context as diabatization,30 projection upon localized states is of course commonplace in the interpretation of the DFT electronic structure. For example, the (angular-momentum) projected DOS,52,53 Mulliken partial charge analysis54 and Crystal Orbital Hamilton Population (COHP) bond order analyses,53,55–59 as well as methods based on maximally localized Wannier states,60–62 all employ local representations. Such projection schemes are however widely known to exhibit nontrivial convergence behaviors and depend on the choice of localized basis, which constrains them often to qualitative, though highly insightful use.52,53,55,56,63 In order to predict observable ET lifetimes, we therefore directly confront convergence aspects known to limit the quantitative application of electronic couplings, including the non-unique choice of localized (diabatic) basis,30,31,53,55,56,58,63 dependence of couplings on the underlying atomic basis set,31,55,63,64 Brillouin-zone sampling, and finite slab size effects.47,63,65,66 By calculating electronic couplings explicitly, we aim for enhanced interpretability compared to other highly accurate first-principles models,39,47,51 which nevertheless rely on the sp-band approximation for their mechanistic interpretation. Our key tool for this is the chemisorption function, from the celebrated model of Newns and Anderson.32,67–69 The chemisorption function derives the energetic broadening of an adsorbate resonance explicitly in terms of the couplings to the surface. This broadening in turn corresponds to the lifetime of an electron in the state and can be compared directly with experimental ET lifetimes, offering a crucial bridge between observable lifetimes and non-observable couplings. Second, by critically validating the scheme against resonant ET measurements of argon on five transition metal surfaces and by invoking rules for electronic couplings based on phase arguments from STM, we seek to provide a reliable and general protocol for calculating electronic couplings between adsorbate and substrate bands over the entire energy scale, noting again their potential for broader use.

In the following, we describe how we use the Projection-Operator Diabatization (POD)66 method to extract the argon/surface electronic couplings from Kohn–Sham DFT. We build here on insights from our earlier work, which offers a variant, POD2GS,31 with improved basis set convergence properties. We additionally build on earlier applications of POD, which calculated electronic coupling for photoexcited systems on surfaces.64,66,70–76 While these works limited themselves to nonperiodic cluster models, or periodic models within a Γ-point approximation, we demonstrate within our scheme that it is possible and indeed necessary to consider couplings throughout the full Brillouin zone to achieve a convergent, and thus physically interpretable, chemisorption function of the surface.

The chemisorption function67 describes the energetic broadening or linewidth of an adsorbate frontier state |ψd⟩ (donor) as it interacts with a continuum of states {|ψa⟩} (acceptors) on a surface, cf. Fig. 1,
(1)
where Ĥ is the Hamiltonian of the full system and ϵa the energy of the acceptor states. Its relation to CT processes is given, in a weak coupling regime,64 by the golden rule transition for a donor state of energy ϵd,
(2a)
with the lifetime then given by
(2b)

The chemisorption function thus emerges as the half-width at half-maximum (HWHM) of a Lorentzian broadening of the energy level ϵd. It is also referred to as the weighted density of states (WDOS),77 as well as a hybridization function,63 which emphasizes that it can in principle be evaluated at any energy ϵ. Of course, the relation to the lifetime of the donor state strictly holds only for Δ(ϵd), but approximately one can also evaluate the chemisorption function for other energies, for example, at the experimental resonance energy rather at a computed donor state energy ϵd.

Importantly, the chemisorption function thus allows us to approximately explore and predict the lifetime as a function of the energy of the donor state. To adapt its definition for practical use in a periodic DFT supercell model for the extended metal surface, two modifications are needed. First, due to the finite nature of periodic slab models and the employed basis set, only a discrete quasicontinuum of N levels is present. We therefore replace the Dirac δ-functions in Eq. (1) by artificial broadening functions, gσ(ϵϵa), to give Δ(ϵ) a smooth and continuous behavior. Specifically, we choose a Lorentzian distribution with HWHM σ = 0.2 eV and demonstrate in the supplementary material that the determined lifetimes do not depend on the choice of the latter parameter. Second, addressing the Bloch description of the electronic states in supercell models, the explicit state summation in Eq. (1) needs to be replaced by an appropriate integration over the system’s Brillouin zone (BZ).78 Already including the Lorentzian smoothing, we correspondingly arrive at k-resolved chemisorption functions
(3)
with the full chemisorption then resulting from BZ integration79,
(4)
In practice, this integral is performed by summing over the k-points with weights wk of a regular Monkhorst–Pack grid in the irreducible BZ.80 We note this uniform integration over the BZ is thus consistent with the conventional treatment of periodic boundary conditions within DFT,52 e.g., in population and bond order analysis55–57 and previous ET studies on surfaces.47,81 Furthermore, uniform BZ integration implies a uniform delocalization of the donor wavepacket in reciprocal space and thus a strongly localized donor state in real space. Due to the low-coverage model employed here—cf. Sec. III—we find the donor state indeed to be strongly localized as depicted in Fig. S26 of the supplementary material, but anticipate the need for a more involved k-space sampling scheme at higher coverages.
The donor |ψd⟩ and acceptor {|ψa⟩} states are the key inputs to compute the three quantities determining the chemisorption function
(5a)
(5b)
(5c)

Newns67 originally defined |ψd⟩ and {|ψa⟩} simply as the Bloch eigenfunctions of non-interacting (isolated) donor and acceptor fragments. We will term this the (fragd, fraga) diabatic basis, consisting in the present application case of the Bloch eigenstates of an infinitely periodic isolated argon layer and an isolated metal surface, both described in a periodic-boundary condition supercell.

For a more quantitative description of the ET process, we instead here rely on a diabatization procedure, which accounts for some weak hybridization between the donor and acceptor system. Specifically, we will draw on our previously described projection-operator diabatization approach, which is based on the partitioning and block-diagonalization of the Kohn–Sham DFT Hamiltonian of the full interacting system (POD2), and the subsequent orthogonalization of the diabats via a Gram–Schmidt (GS) procedure (POD2GS).31 

For ET from weakly physisorbed argon to a metal surface, we thus generate a POD2 diabatic spectrum {|ψd⟩} for the argon donor by extracting and diagonalizing its block in the Hamiltonian of the combined system [Fig. 2(b)]. This POD2 spectrum is found to contain a suitable argon 4s-like lowest unoccupied molecular orbital (LUMO) state, and we refer to this donor state as POD2d. Keeping to the original spirit of Newns, this is combined with the eigenstates of the isolated surface to yield the (POD2d, fraga) diabatic basis. Note that a full POD2-based approach where also the acceptor is described in the interacting picture (POD2d, POD2a) would not show appreciable differences to the (POD2d, fraga) representation due to the acceptor being only very weakly influenced by the donor.

FIG. 2.

Illustration of the diabatization procedure to obtain the chemisorption function Δk(ϵ) at the Γ-point for the case of the low-coverage argon model on ferromagnetic Fe(110) shown in (a). (b) Schematic of the quasi block-diagonal representation of the Kohn–Sham Hamiltonian achieved upon projection onto the (POD2d, fraga)GS diabatic basis. (c) Energy spectra of the majority-spin donor and acceptor diabatic blocks compared to the original adiabatic Kohn–Sham spectrum (occupied states in black, unoccupied states in blue). An Ar4s-like state (red solid line) emerges close to the experimental resonance (red dashed line) in the diabatic donor spectrum. Additionally, a ghost state resulting from basis set superposition appears at 1.6 eV below the Fermi level (black dashed line). (d) Real-space representation of the POD2D argon donor state, revealing its 4s4pz-type hybridization (isovalue 0.01 e1/2Å−3/2), as well as of four select acceptor states indicated in panel (e) (isovalues 0.05 (side views) and 0.02 (top view) e1/2Å−3/2). In each case, the electronic coupling Had to the donor and the (Mulliken) character are indicated. (e) Electronic couplings Had between the donor and all acceptor states are shown on the energy axis, colored according to the dominant character of the acceptor state in the metal slab (s = yellow, p = red, and d = blue). The resulting chemisorption function Δk(ϵ), cf. Eq. (3), is shown as a black line. The predicted Ar donor state energy ϵd is marked with a red vertical line. For comparison, the k = Γ DOS of the slab is shown light gray.

FIG. 2.

Illustration of the diabatization procedure to obtain the chemisorption function Δk(ϵ) at the Γ-point for the case of the low-coverage argon model on ferromagnetic Fe(110) shown in (a). (b) Schematic of the quasi block-diagonal representation of the Kohn–Sham Hamiltonian achieved upon projection onto the (POD2d, fraga)GS diabatic basis. (c) Energy spectra of the majority-spin donor and acceptor diabatic blocks compared to the original adiabatic Kohn–Sham spectrum (occupied states in black, unoccupied states in blue). An Ar4s-like state (red solid line) emerges close to the experimental resonance (red dashed line) in the diabatic donor spectrum. Additionally, a ghost state resulting from basis set superposition appears at 1.6 eV below the Fermi level (black dashed line). (d) Real-space representation of the POD2D argon donor state, revealing its 4s4pz-type hybridization (isovalue 0.01 e1/2Å−3/2), as well as of four select acceptor states indicated in panel (e) (isovalues 0.05 (side views) and 0.02 (top view) e1/2Å−3/2). In each case, the electronic coupling Had to the donor and the (Mulliken) character are indicated. (e) Electronic couplings Had between the donor and all acceptor states are shown on the energy axis, colored according to the dominant character of the acceptor state in the metal slab (s = yellow, p = red, and d = blue). The resulting chemisorption function Δk(ϵ), cf. Eq. (3), is shown as a black line. The predicted Ar donor state energy ϵd is marked with a red vertical line. For comparison, the k = Γ DOS of the slab is shown light gray.

Close modal
In neither the (fragd, fraga) nor the (POD2d, fraga) basis, the donor and acceptor states are necessarily orthogonal to each other, as they should be for a numerically stable and formally correct value of Had according to tight binding theory.31,82,83 To remove any finite overlap
(6)
we therefore GS orthogonalize each donor–acceptor state pair.31 As a consequence, the acceptor acquires a small orthogonalization tail, and the acceptor states are no longer orthogonal with respect to each other. For the present application, the effect of this orthogonalization is relatively small though. Consistent with earlier work,49 the POD2d donor state has a slightly hybridized 4s4pz nodal structure (cf. Fig. 2) and is already nearly orthogonal to the surface as a response to Pauli repulsion interactions. The orthogonalization leads thus only to a slight reduction of the electronic couplings, similar to the findings in earlier studies of molecular systems.31,82 In the following, we shall refer to the two orthogonalized diabatic basis sets as (fragd, fraga)GS and (POD2d, fraga)GS and compare the effect of the diabatization method on the electronic couplings and ET lifetimes below.
A key feature of the POD2 diabatization approach is that it can be performed at arbitrary points in k-space, thus giving full access to the k-resolved chemisorption functions Δk(ϵ) and allowing to properly perform the BZ integration in Eq. (4). Note, though, that this feature is not unique to the POD2 method, but could quite simply be included in other schemes as well. For example, other Hamiltonian fragmentation methods, such as fragment orbital DFT (FO-DFT)84,85 or frozen-density embedding,86,87 could be modified exactly analogously to the POD2 variants. Indeed, the (fragd, fraga)GS basis is here computed with such a modified FO-DFT.88 Similarly, though possibly with mildly more effort, other diabatization approaches, such as constrained DFT30 or the analytic overlap method,89 can be adapted to allow a proper integration over reciprocal space. Yet, in practice, k-point sampling has rarely been used in the context of adsorbate–substrate CT studies to date and, to the best of our knowledge, never been used to computing k-dependent couplings from first principles.22,79 We believe the reasons for this to be twofold. First, most of these methods were initially developed to study molecular systems,28 which are either non-periodic or which can be very well approximated through cluster models. The second reason is that even inherently periodic systems were in the past often either approximated through cluster models or simply treated at the Γ-point. In the present context and nomenclature, this neglect of appropriate BZ integration would correspond to approximating the chemisorption function as
(7)
This can be further simplified by assuming a constant effective coupling between the donor and the acceptor states to ultimately reduce Eq. (3) to
(8)
with DOS(ϵ) denoting the density of states of the acceptor. This is equivalent to the famous golden rule of Fermi66,74,90 and allows the interpretation of excited state lifetimes in terms of the DOS.38,48 Note, though, that this simple picture would not yield quantitative results for extended systems because the DOS does not converge with system size as it scales with the total number of electrons. Instead, one could replace the DOS with a localized version as in the Tersoff–Hamann approach.22,23

Similarly, the prevalent sp-band model of ET would correspond to simply setting to zero all couplings Had,k to d-like states in the acceptor spectrum (e.g., according to a Mulliken analysis). The resulting sp-WDOS differs only from the commonly applied sp-band models by computing the remaining couplings from first principles instead of approximating them, either trivially as in the popular Tersoff–Hamann scheme,21–23,26 or as the overlap of specifically constructed donor and acceptor states.47,51 Note that such an approach allows the interpretation of ET lifetimes in terms of the relative position of the resonance energy in the surface-projected bandgap in the sp-states, yielding a simple descriptor for ET processes.47,51 The rationale behind this descriptor is that the sp-bandgap may influence lifetimes both by reducing the total DOS of acceptor states at k = Γ and the overall coupling by forcing transfer to states with increasingly large k values.47,51 These, for a fixed energy, have less momenta in the k direction and hence less overlap and coupling with the adsorbate.

Focusing solely on CT, all of the above methods, including our WDOS approach, give direct access to excited state lifetimes. Similar results can, of course, be achieved by direct time propagation of the initial state.46,49,50 Most such prior approaches employed model Hamiltonians in the construction of the time propagator and indeed our diabatic Hamiltonian could be similarly employed.

Importantly, in this work, we compute the terms in Eq. (4) from a slab model of the acceptor system, cf. Sec. III. While slab models have found numerous applications in surface science their finite representation of, in principle, semi-infinite metals might lead to artificial gaps in the acceptor spectrum.47,51,65 In our approach, this is addressed by converging the results with respect to slab depth and by introducing a small broadening in Eq. (3). Given that only surface states are expected to show non-zero coupling to the donor, this approach seems justified. Indeed our results (cf. Fig. S3 in the supplementary material) show the WDOS to be well converged already at four-layered slabs. There is, however, another way to include the effects of the continuum on a material’s surface states. Using surface Green’s functions constructed from a Hamiltonian combining bulk and slab DFT calculations, Echenique and co-workers47,65 directly computed resonance lifetimes,51 with only the smallest of artificial broadenings for numerical reasons. A similar semi-infinite Hamiltonian could potentially be used to compute the couplings Had,k.

Finally, we note a variety of existing methods for calculating the chemisorption function in the correlated many-body physics community,60–63,91–93 where it is referred to as the imaginary component of the hybridization function. There, the chemisorption function for periodic systems is often calculated in a localized basis of Wannier functions, either by first calculating electronic couplings or directly through a Green’s function approach. While such a localization scheme is comparable to the diabatization schemes employed in charge transfer studies, our diabatic basis differs in that the diabats will remain fully periodic Bloch states, though localized on either the slab or adsorbate. A deeper comparison of approaches for calculating the chemisorption function will be relegated to future work.

We apply our chemisorption-function based scheme to low-coverage models of argon monolayers on five low-index transition metal surfaces, namely, the magnetic Fe(110), Co(0001) and Ni(111), as well as the non-magnetic Ru(0001) and Pt(111) surfaces. The extended surfaces are described in periodic boundary supercell models, comprising slabs consisting of eight layers and a vacuum region in excess of 40 Å. Ar is adsorbed at the metal surfaces’ top position on one side of the slab such that there is one adsorbate per surface unit-cell. At the employed large surface unit-cells detailed in Fig. 2, this then describes a dilute overlayer with lateral Ar–Ar distances exceeding 8 Å.

The electronic structure of all systems is described at the DFT level using the Perdew–Burke–Ernzerhof (PBE) functional94 with Tkatchenko–Scheffler95 dispersion correction. All calculations are performed with the all-electron DFT package FHI-aims.96,97 Tier 1 numeric atomic orbital basis sets are employed for all metals, a Pople 6-311+G** valence triple-ζ basis set from Basis Set Exchange98 is used for argon. Note that this split basis approach is necessary because while standard FHI-aims basis sets are able to very accurately represent the electronic structure of the metallic surfaces, they would yield erroneous unoccupied states of the Ar adsorbate. Especially, the crucial LUMO state would have shown a wrong level alignment and orbital geometry. The chosen triple-ζ Ar basis fully remedies this as depicted in Fig. 2. Real-space quantities are represented on FHI-Aims’ tight integration grids while BZ integration is performed on a (4 × 4 × 1) regular k-point grid.80 The Ar adatom and the topmost four metal layers are fully relaxed until residual forces fall below 5 meV/Å. A half core-hole is then included on Ar to simulate the effect of the core-excitation in the system’s initial state. Here, we use the excited transition potential (XTP) occupation constraint described in previous work, which is a charge-neutral and implicit approach.99,100 A half electron from the adsorbed argon core level is promoted to the Fermi level of the system and is mostly delocalized over the metal slab, leaving the argon with a Hirshfeld charge of +0.43e. Locally, on the adsorbed argon atom, this state thus corresponds to the charged transition potential (TP) method, which itself is also commonplace in the simulation of spectra of core-excited molecules.99 The constraint scheme was found to give superior alignment of the Ar4s* resonance on the surface with respect to the Fermi level (see Table S6), in comparison to, e.g., ground state (GS) and full core-hole (XCH) approaches. It also yields the suitably hybridized donor wavepacket upon diabatization in the (POD2d, fraga)GS scheme. To generate the pure 4s state of an isolated argon monolayer, for use in the (fragd, fraga)GS scheme, a full core hole (XCH) on the argon was used, with the 4s level occupied.

For the calculation of the more delicate electronic couplings, a single-point calculation is finally performed. The basis set setting of the topmost metal layer is increased to tier 4. Note that all other metal layers are left at tier 1 and the less crucial integration grids were used at light settings for reasons of computational efficiency. The SCF-converged Hamiltonian and overlap matrices are printed at each point on a (12 × 12 × 1) k-grid, and diabatization routines are performed in an external program based on SciPy.101 Consistent with the employed low-coverage model and the corresponding excitation of the quasi-isolated initial core electron to arbitrary momenta, the wavepacket of the excited Ar4s4pz state is assumed to be uniformly distributed in k-space. This uniform distribution was exploited in earlier work, which found the construction of non-trivial (Wannier) wavepackets, with increased weight near the k = Γ point, necessary only at higher coverages.47,51

As demonstrated in Figs. S1, S9, and S10 of the supplementary material, the resulting chemisorption function and the lifetimes derived from it are fully converged with respect to the employed finite k-grid, the Lorentzian smearing σ and the metal basis set. Unfortunately, convergence with respect to the Ar basis set is not that straightforward. Addition of further diffuse functions will increase the basis set superposition. Any such participation of argon basis functions in describing the slab density will lead to the appearance of a ghost state upon partitioning and diagonalization of the argon block in the POD2 scheme. This state at 1.6 eV below the Fermi level is already seen in the Γ-point in Fig. 2(c) at the triple-ζ basis set and fortifies for larger Ar basis sets. We find such larger basis sets to also decrease the 4s4pz hybridization of the true donor state. To avoid these problems, we therefore stick to the triple-ζ Ar basis set, which yields a comparable hybridization for the Ar@Fe(110) system as found in earlier work.47,49–51 An in-depth analysis of the influence of the Ar-basis can be found in the supplementary material.

We also note that, in direct opposition to earlier results,47 the k-grid integrated chemisorption function converges well with slab depth. As depicted in Figs. S3 and S4 of the supplementary material, already four layers of metal atoms show a mostly converged WDOS with no noticeable improvement beyond eight layers.

Finally, we also examined the convergence of the WDOS with donor coverage. As mentioned above, we target a dilute limit in our work, while in the experiment Ar actually forms a dense monolayer. Yet, at the small fluences of incidence photons only a small fraction of Ar atoms will be excited at any one time. Thus, the excited donor states will essentially be localized and dilute in real space, rather than forming an excited band in the Ar overlayer. As depicted in Fig. S6 of the supplementary material, we find a small variation of computed lifetimes going from the c(2 × 6) coverage used here to the dilute p(5 × 5) overlayer, similar to earlier studies,47 which however does not significantly influence resulting trends.

1. Γ-point results

To illustrate our approach, we first examine the properties of the chemisorption function at the Γ-point, Δk(ϵ), for the low-coverage c(2 × 6) Ar/Fe(110) system. These results are summarized in Fig. 2. Upon projection onto the (POD2d, fraga)GS diabatic basis, the Kohn–Sham Hamiltonian nearly block-diagonalizes, with the acceptor block only approximately diagonal as a result of the GS orthogonalization procedure. The projection yields diabatic energy levels, which are comparable to the levels of the original adiabatic Kohn–Sham spectrum as shown in Fig. 2(c). Most important for the ensuing discussion of adsorbate to surface ET, a donor state appears in the diabatic argon spectrum at ϵd = 3.13 eV above the Fermi level. This is close to the experimental value of 2.97 eV43 and we note that the half core-hole constraint is crucial to achieve this good alignment. As apparent from the real-space visualization in Fig. 2(d) this Ar donor state has a hybridized 4s4pz character and thus contains a nodal plane parallel to the surface, consistent with earlier theoretical49 and experimental45 work.

Figure 2(e) shows the computed electronic couplings Had,k to all acceptor states in the diabatic basis (for comparison, the couplings at k=S̄ are shown in Fig. S24 of the supplementary material. These couplings span a rather large range of strengths within each angular momentum channel. Most d-states show rather small couplings, seemingly supporting the prevalent sp-band model of ET. Yet, they are not zero, and in fact do show similar strengths as s- and p-like states in the energy region close to the donor level. In fact, the high symmetry of the Ar4s4pz donor state, as well as its position on the on-top site, allow to nicely understand these varying couplings in terms of textbook rules of orbital phase cancellation. Figure 2(d) illustrates this for four select acceptor states. Acceptor state (i) is a bulk d-state with hardly any presence in the surface layer and a correspondingly vanishing coupling. State (ii) is a p-like state, which is in-phase with the donor orbital in the k direction, resulting in strong overlap and coupling. State (iii) is the same state but out-of phase (antibonding) in the k direction, resulting in a cancellation of phase and a corresponding zero coupling. Finally, state (iv) is an intermediate case with a corresponding intermediate coupling.

In contrast to a mere (angular momentum projected) DOS the computed chemisorption function Δk(ϵ) shown in Fig. 2(e) appropriately accounts for these varying couplings, with its value Δk(ϵd) at the donor level energy then yielding the broadening for the Ar adsorbate and the concomitant ET lifetime. Coincidentally, the Γ-only WDOS displayed Fig. 2(e) already seems to provide, when evaluated at the donor resonance energy, lifetimes directly comparable to experimental results43,48,102 (see Table S5). Notwithstanding, we stress that the actual agreement here is in fact largely fortuitous. Primarily due to band folding, the WDOS merely evaluated in the Γ-point approximation is highly sensitive to the employed size of the periodic supercell. We illustrate this in Fig. S2 of the supplementary material for a range of cell sizes. It is found that the chemisorption function within a Γ-point approximation only slowly converges with cell size, making it insufficient for the prediction of lifetimes at reasonable slab sizes. To really achieve a quantitatively converged chemisorption function, a proper integration over the BZ is indispensable.

2. Brillouin-zone integrated results

Analogous to the Γ-point results in Fig. 2(e), we show in Fig. 3 now the Ar4s4pz chemisorption function on five transition metal surfaces after appropriate integration over the full Brillouin zone (full lines), cf. Eq. (4). To highlight the influence of BZ integration we depict the respective Γ-point results (dashed lines) alongside the k-converged chemisorption functions. Next to the noticeable differences between the two functions, we note that the BZ-integrated chemisorption function is actually largely independent of the slab size, unlike the Γ-point result. As depicted in Figs. S3 and S4 in the supplementary material, we find the BZ-integrated WDOS to converge very well with slab depth, showing very little to no change beyond four layers. Furthermore, we find significant differences between the integrated and Γ-only WDOS’ especially at the argon donor energies. These, in turn, lead to noticeably different predictions for the lifetimes [cf. Eq. (2)] of the core-excited Ar state (see Table S5).

FIG. 3.

Brillouin-zone integrated chemisorption function (solid line) for the Ar4s4pz compared to the chemisorption function obtained at Γ-point (dashed line). The (POD2d, fraga)GS method is used. The BZ-integration is performed on a (12 × 12 × 1) k-grid. Both functions contain a Lorentzian broadening of 0.2 eV. The experimental and calculated resonance energies are shown with vertical red and black lines, respectively. The DOS of the metal surface layer is shown for comparison.

FIG. 3.

Brillouin-zone integrated chemisorption function (solid line) for the Ar4s4pz compared to the chemisorption function obtained at Γ-point (dashed line). The (POD2d, fraga)GS method is used. The BZ-integration is performed on a (12 × 12 × 1) k-grid. Both functions contain a Lorentzian broadening of 0.2 eV. The experimental and calculated resonance energies are shown with vertical red and black lines, respectively. The DOS of the metal surface layer is shown for comparison.

Close modal

The k-point integrated and converged chemisorption function now allows us to compute ET lifetimes [cf. Eq. (2)] comparable to experiment as illustrated in Fig. 4(a) for the transfer from Ar*(2p3/214s) to the five metal surfaces. Additionally, the results are summarized in Table S6 of the supplementary material. Experimental lifetimes, lifetime errors (when reported), and resonance energies for ferromagnetic systems Fe(110), Co(0001) and Ni(111) were taken from the same publication,43 while those of Ru(0001)48 and Pt(111)102 are from separate, older works. For the core-hole measurement on Ar/Pt(111),102 no lifetime was reported. It can, however, be estimated with the same method as used in Ref. 48. The calculated lifetimes have an overall mean relative signed error of −6% and mean signed error of −0.14 fs compared to the experimental values. The qualitative and quantitative differences in lifetimes over the surfaces, and, for ferromagnetic systems, between spin channels, are captured excellently.

FIG. 4.

(a) Calculated electron transfer lifetimes (vertical bars) compared to experimental values (dots, error bars shown when reported). The lifetimes may be understood in terms of parameters describing the systems: (b) the value of the chemisorption function (WDOS) and (c) the DOS of the surface layer at the resonance energy, as well as their angular momentum decompositions, (d) the alignment of the surface bands relative to the resonance, and (e) the adsorption height of the argon above the on-top site on each surface. Correlations are analyzed in Fig. 5, and Figs. S13 and S14 of the supplementary material.

FIG. 4.

(a) Calculated electron transfer lifetimes (vertical bars) compared to experimental values (dots, error bars shown when reported). The lifetimes may be understood in terms of parameters describing the systems: (b) the value of the chemisorption function (WDOS) and (c) the DOS of the surface layer at the resonance energy, as well as their angular momentum decompositions, (d) the alignment of the surface bands relative to the resonance, and (e) the adsorption height of the argon above the on-top site on each surface. Correlations are analyzed in Fig. 5, and Figs. S13 and S14 of the supplementary material.

Close modal

As we demonstrate in the supplementary material (Table S5 and Fig. S11), the choice of diabatization scheme can have a significant effect on the accuracy of results. Overall, the (POD2d, fraga)GS diabatization scheme yields the most accurate resonance energies and lifetime predictions when compared with the alternative diabatization scheme (fragd, fraga)GS, and to the schemes applied without orthogonalization. The constrained 4s donor (fragd) generally couples more strongly to the surfaces than the hybridized POD2 donor state, leading to shorter and underestimated lifetimes. Gram–Schmidt orthogonalization of the acceptor diabats reduces the magnitude of couplings in both methods and leads to longer lifetimes. While for the magnetic systems, the (fragd, fraga)GS method seemingly yields lifetimes that agree well with experiment, the method fails for non-magnets (Fig. S17) and consistently underestimates the resonance energy. It is prudent, therefore, in such studies to consider the sensitivity of predictions to the employed diabatization scheme.

Having identified (POD2d, fraga)GS as the most suitable diabatization method, and established the agreement of our approach with experiment, we now examine the contributions of different acceptor states, specifically focusing on the angular momenta. For comparison, Figs. 4(b) and 4(c), respectively, depict the WDOS and surface-DOS evaluated at the resonance energies and both resolved by angular momentum contributions. To this end, we employ a Mulliken analysis of both the DOS and WDOS (cf. Eq. S16 in the supplementary material). The angular-momentum components of the WDOS (which we shall call the projected-WDOS or pWDOS) thus contain the weights of each state by character. The pWDOS can thus be seen as a projected DOS weighted by couplings, and thus shows the dominant character of the states that participate in the overall lifetime broadening (see Fig. S15 of the supplementary material). It does not necessarily show the character that is associated with strong coupling (indeed, coupling and character are nontrivially related, as was seen in Fig. 2).

Figure 4(b) shows that states that participate in the overall lifetime broadening have, on average, significant d as well as sp-character. This is of course consistent with the hybridized nature of states at the resonance energy. However, we see upon closer examination that it is the changing d-character of the states that is found to decisively determine the differences in lifetimes over the surfaces and spin channels, while the sp-character is relatively unchanged.

Considering other popular modes of interpreting ET lifetimes, the DOS depicted in Fig. 4(c) indeed bears some qualitative resemblances to the WDOS. Importantly, the WDOS, due to its inclusion of phase effects, captures more closely the variations among the surfaces, in closer agreement with the experimental lifetimes than the DOS. In particular, the DOS fails to capture the difference between the shortest (Ru) and longest (Pt) lifetimes. The WDOS, on the other hand, shows a much higher value on Ru than Pt, despite similar DOS values on the surfaces. As argon adsorbs approximately at the same height on these two systems—cf. Fig. 4(e), the large difference in lifetimes on Ru and Pt is strongly related to different electronic couplings, determined both by the spatial extent of the surface wavefunctions and specific phase effects.

Finally, in Fig. 4(d), we depict the onsets of the sp-gap and the d-band together with the resonance energy. The sp-gap onset is thereby evaluated from the primitive band structures of all five surfaces (cf. Fig. S27 in the supplementary material), while the d-band edge is extracted from the materials’ surface DOS at a cutoff of 0.2 (states/eV/atom). Both, the proximity of the resonance to the sp-bandgap onset and to the d-band appear to correlate well with the lifetimes. Yet, we note, for example, that for Fe, very different sp-bandgap onsets appear in the majority and minority spin channels, while the sp-DOS and sp-WDOS are nearly identical in both spin channels, indicating little influence of bandgap alignment on the couplings or DOS of the sp-channel (indeed, the parabolic gap should have a nearly constant DOS similar to the 2D free electron gas). In contrast, the difference of the d-band edge between the spin channels corresponds to a large difference in both the d-band DOS and WDOS.

For a clearer analysis of such relationships, we correlate the ET lifetimes against features of the electronic structure in Fig. 5. We find that the calculated lifetimes correlate strongly with (coupling-weighted) d-character of the acceptor states, and weakly with their sp-character. In Figs. 5(a) and 5(b), the (inverse) WDOS d-channel is seen to correlate strongly with the lifetimes, while the sp-channel of the WDOS correlates poorly. This correlation is even worse for the investigated magnetic systems. Furthermore, coming back to the simple Tersoff–Hamann picture, the ET rate should be proportional to the local DOS at the probe coordinate (approximated here as the DOS of the surface layer). In Figs. 5(c) and 5(d), we correlate the lifetimes against the inverted d and sp components of the DOS, respectively. Thereby, we observe a very slight correlation between the lifetimes and inverse DOS of the d-channel, but none at all for the sp-DOS, consistent with the lack of phase information in this simple picture. Additionally, while our predicted lifetimes correlate strongly with the proximity of the resonance to the sp-bandgap onset as depicted in Fig. 5(f), we find them to correlate poorly with the sp-component of the DOS [Fig. 5(d)] and the sp-WDOS [Fig. 5(b)]. Instead, there is a strong correlation of the lifetimes with the d-band edge, which may be explained by an increasing diffusivity of states with energy,47 or simply by the increase in the d-band DOS magnitude at the d-band edge. The energy of the d-states is captured in the d-band edge descriptor, and indeed strongly correlates with the d-band WDOS, cf. Fig. S14(c) in the supplementary material.

FIG. 5.

The predicted electron transfer lifetimes from Fig. 4 are tested for correlation with electronic structure properties of the surfaces. A linear regression is shown for all datapoints (black line) and separately for magnetic systems (violet line). The angular-momentum components of the WDOS are depicted: (a) contains the WDOS weighted by the d-character of states and (b) the WDOS weighted by the sp-character of states. The projected DOS of the surface is shown: (c) contains the DOS weighted by the d-character and (d) the DOS weighted by the sp-character of states in the surface layer. Finally, the proximity of the resonance energy to the edge of the d-band is depicted in (e), with the onset energy of the sp-bandgap depicted in (d).

FIG. 5.

The predicted electron transfer lifetimes from Fig. 4 are tested for correlation with electronic structure properties of the surfaces. A linear regression is shown for all datapoints (black line) and separately for magnetic systems (violet line). The angular-momentum components of the WDOS are depicted: (a) contains the WDOS weighted by the d-character of states and (b) the WDOS weighted by the sp-character of states. The projected DOS of the surface is shown: (c) contains the DOS weighted by the d-character and (d) the DOS weighted by the sp-character of states in the surface layer. Finally, the proximity of the resonance energy to the edge of the d-band is depicted in (e), with the onset energy of the sp-bandgap depicted in (d).

Close modal

We find that similar observations can be made when interpreting the energy dependence of the predicted electron transfer lifetimes. While the connection between the WDOS and excited-state lifetimes only strictly holds at the resonance energy cf. Fig. 2, assuming an unaffected nature of the donor state, one could use it to interpret experimental findings for other incidence energies. They show growing lifetimes for larger incidence energies beyond the resonance.35,48,103 Furthermore, local maxima in the ET rates below the resonance have been observed experimentally102 on Ar/Pt(111) and in accurate Green’s function-based simulations47 of Ar/Ru(0001). Both of these trends are reproduced in the WDOS as depicted in Figs. 3 and S19 of the supplementary material. Theoretically, the energy dependence of lifetimes has alternatively been attributed to the declining DOS of the d-states (already in 1998 by Menzel and co-workers34,103), or to the effect of the bandgap in the sp-states originally proposed by Gauyacq and Borisov.47,49,50 Upon examination of the WDOS components over the energy scale, as shown in Fig. 6 for Ar/Fe, we find that the energy dependence at resonance arises from significant contributions from both d- and sp components, combining the two explanations proposed earlier. The onset of the sp-bandgap, which is marked with an orange line in Fig. (6), is associated with a local maximum in the WDOS, followed by a sharp monotonic decay that differs from the oscillatory behavior observed at lower energies. The decay behavior is consistent with the decay of couplings to sp-type states with k in the gap region proposed earlier.47 While we relegate a more detailed study of such dependencies to future work (the present model is based on the folded band structure, cf. Fig. S26 of the supplementary material), we do observe this suppression of the p-channel WDOS in the bandgap region of all systems (Fig. S16 in the supplementary material)—against quite different behavior in the PDOS—indicating a systematic relationship between the WDOS and the band structure.

FIG. 6.

The chemisorption function of the Ar/Fe(110) c(2 × 6) system is compared to the surface-projected DOS, and to the band structure of the surface primitive cell. This is shown for the majority spin direction (left three panels) and minority spin (right three panels). The nontrivial relationship between the DOS and the WDOS arises due to the electronic couplings, which contain, e.g., phase cancellation effects. The sp-bandgap onset, marked with an orange horizontal line, is associated with a local maximum in the chemisorption function followed by a monotonic decay, consistent with the role of the bandgap identified in earlier work, and with the experimentally measured energy dependence of the ET process. The d-band edge is marked with a blue line, and the experimental and predicted Ar4s resonances are indicated. Other systems shown in Figs. S16 and S27 of the supplementary material.

FIG. 6.

The chemisorption function of the Ar/Fe(110) c(2 × 6) system is compared to the surface-projected DOS, and to the band structure of the surface primitive cell. This is shown for the majority spin direction (left three panels) and minority spin (right three panels). The nontrivial relationship between the DOS and the WDOS arises due to the electronic couplings, which contain, e.g., phase cancellation effects. The sp-bandgap onset, marked with an orange horizontal line, is associated with a local maximum in the chemisorption function followed by a monotonic decay, consistent with the role of the bandgap identified in earlier work, and with the experimentally measured energy dependence of the ET process. The d-band edge is marked with a blue line, and the experimental and predicted Ar4s resonances are indicated. Other systems shown in Figs. S16 and S27 of the supplementary material.

Close modal

Finally, observing the nontrivial relationship of the DOS and WDOS over the energy scale, we can furthermore better understand the (lack of) correlation of the ET lifetimes with the d- and sp-channels established above. For this, we may consider the nature of the couplings themselves. States with strong sp character, which have large spatial extent and presence on the surface, may nevertheless have vanishing couplings due to phase cancellation (cf. Fig. 2). This could be an additional explanation for the suppression of sp-character relative to d-character in the pWDOS in Figs. 6 and S16. The more localized d-bands, on the other hand, couple via their exponential tails and may thus be more consistently affected by phase cancellation than larger sp-like states. Chen related24–26,104 the coupling matrix elements to the gradients, rather than the squared-modulus (density), of participating states. As such, we may anticipate qualitatively different coupling behaviors from the localized (high-gradient) d-band states compared the more diffuse (and uniform) free-electron-like sp-bands above the Fermi level.

These effects are meaningfully elucidated by constraining the donor to have purely 4s character, within the (fragd, fraga)GS method. Upon removal of the 4pz hybridization, depicted in Figs. S16 and S17 of the supplementary material, we recover a much stronger relationship between the d-WDOS and d-DOS over the full energy scale. Nonlinearity between the sp-DOS and sp-WDOS is meanwhile enhanced, confirming the more complex nature of coupling to these states.

In this work, we examine the utility of the Newns–Anderson chemisorption function, also known as coupling-weighted DOS (WDOS), for the prediction and interpretation of ultrafast electron transfer lifetimes. Specifically, we revisit the case of core-excited Ar*(2p3/214s) atoms on five different transition metal surfaces, to find our WDOS results in very good agreement with experimental measurements.

We demonstrate that the WDOS converges well with system size, both lateral and with layer depth. A prerequisite for this convergence and indeed the agreement with experiment was found to be an appropriate sampling of each system’s Brillouin zone. Using the tried and tested k-point scheme by Monkhorst and Pack to integrate the WDOS, we find electron transfer rates from the excited argon to the metal surface to already converge for metal slabs of only four layers. This performance stands in contrast to the previously supposed inadequacy of finite slab models for the calculation of sensitive electron transfer properties, due to inherent finite size effects. Previously, such effects were overcome by combining bulk and slab calculations and using surface Green’s Functions techniques.47,65,66 Our model overcomes the finite size effects present in the Γ-point approximation simply with an integration over the Brillouin zone with a sufficiently dense k-grid and an energy broadening function, yielding a comparable accuracy to previous approaches for the present systems.

A critical component of the WDOS scheme is the calculation of electronic coupling matrix elements Had, for which a numerically stable diabatization procedure is necessary. To this end, we apply the DFT-based projection-operator diabatization scheme POD231 to generate the diabatic basis of the argon donor system, and use the eigenspectrum of the unperturbed surface as the diabatic basis of the acceptor (surface). The two systems are Gram–Schmidt orthogonalized, yielding the diabatic basis we term (POD2d, fraga)GS, and matrix elements are extracted by projecting the Kohn–Sham Hamiltonian onto this basis.

The resulting ab initio chemisorption function is found to be highly expressive, encoding rich information on orbital phase and symmetry, and yielding insight into the decay mechanism beyond a mere prediction of rates. A key advantage of this approach is its ability to explicitly capture complex phase cancellation effects over the full single-particle spectrum, hinting at an under-leveraged potential for diabatic representations of the Kohn–Sham Hamiltonian.30 Across the studied systems, highly nontrivial and nonlinear effects in the couplings are found to modulate the magnitude of the chemisorption function. This contrasts with the anticipated scaling between the electron transfer lifetimes and features of the sp-channel, which arises from a simplified free-electron-like view of couplings,47,49–51,105 as openly recognized in earlier works. Decomposing the WDOS by angular momentum contributions, we find the couplings to have a more consistent and linear effect upon the acceptor states with strong d-character, and it is the degree of d-hybridization of states that is then found to strongly correlate with the electron transfer lifetimes, explaining their surface and spin dependence and contributing to their energy dependence. By treating the DOS, the band structure, and the electronic couplings all in a consistent ab initio level of theory, the WDOS scheme thus resolves long-standing questions (formulated, e.g., by Menzel et al. in 2000,38,48 or later by Gauyacq and Borisov50) on the respective roles of the compact d-states vs the more spatially extensive sp-states in the resonant ET process.

Note that the role of orbital phase cancellation in electron transfer has long been recognized in scanning tunneling microscopy,24,25,104,106 where Chen’s derivative rules were found essential to understanding the contrast mechanism beyond the (phase-less) Tersoff–Hamann model.22,23 Our results for the ET lifetimes are thus perhaps unsurprising. For magnetic systems, the d-band is, after all, that which is responsible for ferromagnetism and the asymmetry of the electronic structure. Furthermore, the observation of a relatively constant contribution of sp-type states to the lifetimes on different surfaces is also well in line with the prevailing view (e.g., in heterogeneous catalysis11,12,33) that the coupling to the sp-band is a constant and largely independent of the substrate.

The present model does, however, also show some departure from prevailing assumptions about the electronic coupling. Noticeably, the phase modulation present in the sp-WDOS gives it a more rich and nontrivial dependence on the incidence energy, compared to the relatively structureless sp-projected DOS. This is in contrast with the common wide-band approximation, which assumes a constant (energy-independent) coupling to the sp-bands.7,11,32,68,107–109 Furthermore, the d-channel of the chemisorption function, which shows a stronger resemblance to the d-band DOS (and indeed is often approximated as such33), still shows nontrivial departures from the contour of the d-band DOS over the incidence energy scale.

Taken together, all of the above observations suggest a stronger role for nontrivial coupling effects in other approximative Newns–Anderson models, such as those forming the staple description of, e.g., many chemical processes on surfaces.4,7,11,15,32,33,110–112 While the present diabatization scheme benefits from the weakly interacting nature of the argon adsorbate, similar schemes could be constructed for chemisorbed systems. Additional possible applications include orbital-resolved scanning-probe microscopy22,24,104,113,114 and strongly correlated impurity models.60,61,63,91,92,115

See the supplementary material for the following: Equations (S1)–(S6) describe the Gram–Schmidt orthogonalization of the diabatic basis, (S7)–(S10) the Mulliken analysis of states, and (S11)–(S16) the angular momentum decomposition of the WDOS. Figures are provided to illustrate convergence of the chemisorption function with respect to: Brillouin zone sampling (S1, S2), slab depth (S3, S4), and basis sets and monolayer coverage [(S5)–(S10)]. Figure S11 and Table S5 compare diabatization schemes. Figures S12–S23 show k = Γ and BZ-averaged diabatic properties comprising the calculation of ET lifetimes on the five slab systems. Figures S25–S27 provide adiabatic DOS and band structures for comparison.

The authors thank Peter Feulner, Wolfgang Domcke, Alan Luntz, and Niri Govind for useful and enlightening discussions. S.G. was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Grant No. RE1509/24-1. H.O. recognizes support by the DFG under the Heisenberg program (Grant No. OB 425/9-1). This work was partially funded by the DFG under Germany’s Excellence Strategy (Grant No. EXC 852 2089/1-390776260).

The authors have no conflicts to disclose.

Simiam Ghan: Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Elias Diesen: Methodology (equal); Software (equal); Writing – review & editing (equal). Christian Kunkel: Methodology (equal); Software (equal); Writing – review & editing (equal). Karsten Reuter: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Resources (equal); Writing – review & editing (equal). Harald Oberhofer: Conceptualization (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
M. D.
Newton
,
Chem. Rev.
91
,
767
792
(
1991
).
2.
R. J.
Cave
and
M. D.
Newton
,
Chem. Phys. Lett.
249
,
15
19
(
1996
).
3.
R. J.
Cave
and
M. D.
Newton
,
J. Chem. Phys.
106
,
9213
9226
(
1997
).
4.
J. K.
Nørskov
,
J. Chem. Phys.
90
,
7461
7471
(
1989
).
5.
B.
Hammer
and
J. K.
Norskov
,
Nature
376
,
238
240
(
1995
).
6.
B.
Hammer
and
J. K.
Nørskov
,
Surf. Sci.
343
,
211
220
(
1995
).
7.
E.
Santos
,
P.
Quaino
, and
W.
Schmickler
,
Phys. Chem. Chem. Phys.
14
,
11224
11233
(
2012
).
8.
R.
Hoffmann
,
Solids and Surfaces: A Chemists View of Bonding in Extended Structures
(
Wiley-VCH, Inc.
,
1989
).
10.
F.
Abild-Pedersen
,
J.
Greeley
,
F.
Studt
,
J.
Rossmeisl
,
T. R.
Munter
,
P. G.
Moses
,
E.
Skúlason
,
T.
Bligaard
, and
J. K.
Nørskov
,
Phys. Rev. Lett.
99
,
016105
(
2007
).
11.
S.
Saini
,
J.
Halldin Stenlid
, and
F.
Abild-Pedersen
,
npj Comput. Mater.
8
,
163
(
2022
).
12.
T.
Bligaard
and
J. K.
Nørskov
,
Chemical Bonding at Surfaces and Interfaces
, edited by
A.
Nilsson
,
L. G.
Pettersson
, and
J. K.
Nørskov
(
Elsevier
,
Amsterdam
,
2008
), pp.
255
321
.
13.
A.
Nilsson
and
L. G. M.
Pettersson
,
Chemical Bonding at Surfaces and Interfaces
, edited by
A.
Nilsson
,
L. G.
Pettersson
, and
J. K.
Nørskov
(
Elsevier
,
Amsterdam
,
2008
), pp.
57
142
.
14.
A.
Gross
,
Theoretical Surface Science: A Microscopic Perspective
(
Springer
,
2003
).
15.
A.
Cao
and
J. K.
Nørskov
,
ACS Catal.
13
,
3456
3462
(
2023
).
16.
B. N. J.
Persson
and
M.
Persson
,
Solid State Commun.
36
,
175
179
(
1980
).
17.
P.
Saalfrank
,
Chem. Rev.
106
,
4116
4159
(
2006
).
18.
D.
Menzel
,
J. Chem. Phys.
137
,
091702
(
2012
).
19.
D.
Menzel
and
R.
Gomer
,
J. Chem. Phys.
41
,
3311
3328
(
1964
).
20.
J. A.
Misewich
,
T. F.
Heinz
, and
D. M.
Newns
,
Phys. Rev. Lett.
68
,
3737
3740
(
1992
).
21.
J.
Bardeen
,
Phys. Rev. Lett.
6
,
57
59
(
1961
).
22.
J.
Tersoff
and
D. R.
Hamann
,
Phys. Rev. B
31
,
805
813
(
1985
).
23.
J.
Tersoff
and
D. R.
Hamann
,
Phys. Rev. Lett.
50
,
1998
2001
(
1983
).
24.
C. J.
Chen
,
Phys. Rev. B
42
,
8841
8857
(
1990
).
25.
C. J.
Chen
,
Phys. Rev. Lett.
69
,
1656
1659
(
1992
).
26.
C. J.
Chen
,
Introduction to Scanning Tunneling Microscopy
(
Oxford University Press
,
2007
).
27.
L. S.
Cederbaum
,
H.
Köppel
, and
W.
Domcke
,
Int. J. Quantum Chem.
20
,
251
267
(
1981
).
28.
H.
Oberhofer
,
K.
Reuter
, and
J.
Blumberger
,
Chem. Rev.
117
,
10319
10357
(
2017
).
29.
C. R.
Moon
,
L. S.
Mattos
,
B. K.
Foster
,
G.
Zeltzer
,
W.
Ko
, and
H. C.
Manoharan
,
Science
319
,
782
787
(
2008
).
30.
T.
Van Voorhis
,
T.
Kowalczyk
,
B.
Kaduk
,
L.-P.
Wang
,
C.-L.
Cheng
, and
Q.
Wu
,
Annu. Rev. Phys. Chem.
61
,
149
170
(
2010
).
31.
S.
Ghan
,
C.
Kunkel
,
K.
Reuter
, and
H.
Oberhofer
,
J. Chem. Theory Comput.
16
,
7431
7443
(
2020
).
32.
D.-Q.
Liu
,
M.
Kang
,
D.
Perry
,
C.-H.
Chen
,
G.
West
,
X.
Xia
,
S.
Chaudhuri
,
Z. P. L.
Laker
,
N. R.
Wilson
,
G. N.
Meloni
,
M. M.
Melander
,
R. J.
Maurer
, and
P. R.
Unwin
,
Nat. Commun.
12
,
7110
(
2021
).
33.
A.
Vojvodic
,
J. K.
Nørskov
, and
F.
Abild-Pedersen
,
Top. Catal.
57
,
25
32
(
2014
).
34.
C.
Keller
,
M.
Stichler
,
G.
Comelli
,
F.
Esch
,
S.
Lizzit
,
D.
Menzel
, and
W.
Wurth
,
Phys. Rev. B
57
,
11951
11954
(
1998
).
35.
D.
Menzel
,
Chem. Soc. Rev.
37
,
2212
2223
(
2008
).
36.
S.
Lizzit
,
G.
Zampieri
,
K. L.
Kostov
,
G.
Tyuliev
,
R.
Larciprete
,
L.
Petaccia
,
B.
Naydenov
, and
D.
Menzel
,
New J. Phys.
11
,
053005
(
2009
).
37.
H.
Schlichting
,
D.
Menzel
,
T.
Brunner
,
W.
Brenig
, and
J. C.
Tully
,
Phys. Rev. Lett.
60
,
2515
2518
(
1988
).
38.
A.
Föhlisch
,
D.
Menzel
,
P.
Feulner
,
M.
Ecker
,
R.
Weimar
,
K. L.
Kostov
,
G.
Tyuliev
,
S.
Lizzit
,
R.
Larciprete
,
F.
Hennies
, and
W.
Wurth
,
Chem. Phys.
289
,
107
115
(
2003
).
39.
A.
Föhlisch
,
P.
Feulner
,
F.
Hennies
,
A.
Fink
,
D.
Menzel
,
D.
Sánchez-Portal
,
P. M.
Echenique
, and
W.
Wurth
,
Nature
436
,
373
376
(
2005
).
40.
S.
Neppl
,
U.
Bauer
,
D.
Menzel
,
P.
Feulner
,
A.
Shaporenko
,
M.
Zharnikov
,
P.
Kao
, and
D. L.
Allara
,
Chem. Phys. Lett.
447
,
227
231
(
2007
).
41.
P.
Kao
,
S.
Neppl
,
P.
Feulner
,
D. L.
Allara
, and
M.
Zharnikov
,
J. Phys. Chem. C
114
,
13766
13773
(
2010
).
42.
H.
Hamoudi
,
S.
Neppl
,
P.
Kao
,
B.
Schüpbach
,
P.
Feulner
,
A.
Terfort
,
D.
Allara
, and
M.
Zharnikov
,
Phys. Rev. Lett.
107
,
027801
(
2011
).
43.
F.
Blobner
,
R.
Han
,
A.
Kim
,
W.
Wurth
, and
P.
Feulner
,
Phys. Rev. Lett.
112
,
086801
(
2014
).
44.
P.
Feulner
,
F.
Blobner
,
J.
Bauer
,
R.
Han
,
A.
Kim
,
T.
Sundermann
,
N.
Müller
,
U.
Heinzmann
, and
W.
Wurth
,
e-J. Surf. Sci. Nanotechnol.
13
,
317
323
(
2015
).
45.
T.
Sundermann
,
N.
Müller
,
U.
Heinzmann
,
W.
Wurth
,
J.
Bauer
,
R.
Han
,
A.
Kim
,
D.
Menzel
, and
P.
Feulner
,
Surf. Sci.
643
,
190
196
(
2016
).
46.
J. P.
Gauyacq
,
A. G.
Borisov
, and
G.
Raşeev
,
Surf. Sci.
490
,
99
115
(
2001
).
47.
D.
Sánchez-Portal
,
D.
Menzel
, and
P. M.
Echenique
,
Phys. Rev. B
76
,
235406
(
2007
).
48.
W.
Wurth
and
D.
Menzel
,
Chem. Phys.
251
,
141
149
(
2000
).
49.
S.
Vijayalakshmi
,
A.
Föhlisch
,
F.
Hennies
,
A.
Pietzsch
,
M.
Nagasono
,
W.
Wurth
,
A. G.
Borisov
, and
J. P.
Gauyacq
,
Chem. Phys. Lett.
427
,
91
95
(
2006
).
50.
J. P.
Gauyacq
and
A. G.
Borisov
,
Phys. Rev. B
69
,
235408
(
2004
).
51.
M.
Müller
,
P. M.
Echenique
, and
D.
Sánchez-Portal
,
J. Phys. Chem. Lett.
11
,
7141
7145
(
2020
).
52.
P.
Kratzer
and
J.
Neugebauer
,
Front. Chem.
7
,
106
(
2019
).
53.
S.
Maintz
,
V. L.
Deringer
,
A. L.
Tchougréeff
, and
R.
Dronskowski
,
J. Comput. Chem.
37
,
1030
1035
(
2016
).
54.
R. S.
Mulliken
,
J. Chem. Phys.
23
,
1833
1840
(
1955
).
55.
R.
Dronskowski
and
P. E.
Bloechl
,
J. Phys. Chem.
97
,
8617
8624
(
1993
).
56.
S.
Steinberg
and
R.
Dronskowski
,
Crystals
8
,
225
(
2018
).
57.
V. L.
Deringer
,
A. L.
Tchougréeff
, and
R.
Dronskowski
,
J. Phys. Chem. A
115
,
5461
(
2011
).
58.
S.
Maintz
,
V. L.
Deringer
,
A. L.
Tchougréeff
, and
R.
Dronskowski
,
J. Comput. Chem.
34
,
2557
(
2013
).
59.
T.
Hughbanks
and
R.
Hoffmann
,
J. Am. Chem. Soc.
105
,
3528
(
1983
).
60.
J.
Kügel
,
M.
Karolak
,
J.
Senkpiel
,
P.-J.
Hsu
,
G.
Sangiovanni
, and
M.
Bode
,
Nano Lett.
14
,
3895
(
2014
).
61.
J.
Kügel
,
M.
Karolak
,
A.
Krönlein
,
J.
Senkpiel
,
P.-J.
Hsu
,
G.
Sangiovanni
, and
M.
Bode
,
Phys. Rev. B
91
,
235130
(
2015
).
62.
M.
Karolak
,
T. O.
Wehling
,
F.
Lechermann
, and
A. I.
Lichtenstein
,
J. Phys.: Condens. Matter
23
,
085601
(
2011
).
63.
M. P.
Bahlke
,
M.
Schneeberger
, and
C.
Herrmann
,
J. Chem. Phys.
154
,
144108
(
2021
).
64.
J.
Li
,
M.
Nilsing
,
I.
Kondov
,
H.
Wang
,
P.
Persson
,
S.
Lunell
, and
M.
Thoss
,
J. Phys. Chem. C
112
,
12326
12333
(
2008
).
65.
D.
Sánchez-Portal
,
Prog. Surf. Sci.
82
,
313
335
(
2007
).
66.
I.
Kondov
,
M.
Čížek
,
C.
Benesch
,
H.
Wang
, and
M.
Thoss
,
J. Phys. Chem. C
111
,
11970
11981
(
2007
).
67.
D. M.
Newns
,
Phys. Rev.
178
,
1123
1135
(
1969
).
68.
P. W.
Anderson
,
Phys. Rev.
124
,
41
53
(
1961
).
69.
T. B.
Grimley
,
Proc. Phys. Soc.
90
,
751
764
(
1967
).
70.
J.
Li
,
I.
Kondov
,
H.
Wang
, and
M.
Thoss
,
J. Phys. Chem. C
114
,
18481
18493
(
2010
).
71.
J.
Li
,
H.
Wang
,
P.
Persson
, and
M.
Thoss
,
J. Chem. Phys.
137
,
22A529
(
2012
).
72.
J.
Li
,
H.
Li
,
P.
Winget
, and
J.-L.
Brédas
,
J. Phys. Chem. C
119
,
18843
18858
(
2015
).
73.
Z.
Futera
and
J.
Blumberger
,
J. Phys. Chem. C
121
,
19677
19689
(
2017
).
74.
V.
Prucker
,
O.
Rubio-Pons
,
M.
Bockstedte
,
H.
Wang
,
P. B.
Coto
, and
M.
Thoss
,
J. Phys. Chem. C
117
,
25334
25342
(
2013
).
75.
V.
Prucker
,
M.
Bockstedte
,
M.
Thoss
, and
P. B.
Coto
,
J. Chem. Phys.
148
,
124705
(
2018
).
76.
M.
Pastore
and
F.
De Angelis
,
J. Am. Chem. Soc.
137
,
5798
5809
(
2015
).
77.
R. E.
Warburton
,
A. V.
Soudackov
, and
S.
Hammes-Schiffer
,
Chem. Rev.
122
,
10599
10650
(
2022
).
78.
C.
Kittel
and
C.-y.
Fong
,
Quantum Theory of Solids
(
Wiley
,
1987
).
79.
S.
Gosavi
and
R. A.
Marcus
,
J. Phys. Chem. B
104
,
2067
2072
(
2000
).
80.
H. J.
Monkhorst
and
J. D.
Pack
,
Phys. Rev. B
13
,
5188
5192
(
1976
).
81.
C.-P.
Hsu
and
R. A.
Marcus
,
J. Chem. Phys.
106
,
584
598
(
1997
).
82.
E. F.
Valeev
,
V.
Coropceanu
,
D. A.
da Silva Filho
,
S.
Salman
, and
J.-L.
Brédas
,
J. Am. Chem. Soc.
128
,
9882
9886
(
2006
).
83.
B.
Baumeier
,
J.
Kirkpatrick
, and
D.
Andrienko
,
Phys. Chem. Chem. Phys.
12
,
11103
11113
(
2010
).
84.
K.
Senthilkumar
,
F. C.
Grozema
,
F. M.
Bickelhaupt
, and
L. D. A.
Siebbeles
,
J. Chem. Phys.
119
,
9809
9817
(
2003
).
85.
H.
Li
,
J.-L.
Brédas
, and
C.
Lennartz
,
J. Chem. Phys.
126
,
164704
(
2007
).
86.
T. A.
Wesolowski
and
A.
Warshel
,
J. Phys. Chem.
97
,
8050
8053
(
1993
).
87.
M.
Pavanello
and
J.
Neugebauer
,
J. Chem. Phys.
135
,
234103
(
2011
).
88.
C.
Schober
,
K.
Reuter
, and
H.
Oberhofer
,
J. Chem. Phys.
144
,
054103
(
2016
).
89.
F.
Gajdos
,
S.
Valner
,
F.
Hoffmann
,
J.
Spencer
,
M.
Breuer
,
A.
Kubas
,
M.
Dupuis
, and
J.
Blumberger
,
J. Chem. Theory Comput.
10
,
4653
4660
(
2014
).
90.
E.
Fermi
,
Nuclear Physics: A Course Given by Enrico Fermi at the University of Chicago
(
University of Chicago Press
,
1950
).
91.
H.
Shinaoka
,
E.
Gull
, and
P.
Werner
,
Comput. Phys. Commun.
215
,
128
(
2017
).
92.
E.
Gull
,
A. J.
Millis
,
A. I.
Lichtenstein
,
A. N.
Rubtsov
,
M.
Troyer
, and
P.
Werner
,
Rev. Mod. Phys.
83
,
349
(
2011
).
93.
D.
Jacob
,
J. Phys.: Condens. Matter
27
,
245606
(
2015
).
94.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
3868
(
1996
).
95.
A.
Tkatchenko
and
M.
Scheffler
,
Phys. Rev. Lett.
102
,
073005
(
2009
).
96.
V.
Blum
,
R.
Gehrke
,
F.
Hanke
,
P.
Havu
,
V.
Havu
,
X.
Ren
,
K.
Reuter
, and
M.
Scheffler
,
Comput. Phys. Commun.
180
,
2175
2196
(
2009
).
97.
I. Y.
Zhang
,
X.
Ren
,
P.
Rinke
,
V.
Blum
, and
M.
Scheffler
,
New J. Phys.
15
,
123033
(
2013
).
98.
B. P.
Pritchard
,
D.
Altarawy
,
B.
Didier
,
T. D.
Gibson
, and
T. L.
Windus
,
J. Chem. Inf. Model.
59
,
4814
4820
(
2019
).
99.
G. S.
Michelitsch
and
K.
Reuter
,
J. Chem. Phys.
150
,
074104
(
2019
).
100.
M.
Leetmaa
,
M. P.
Ljungberg
,
A.
Lyubartsev
,
A.
Nilsson
, and
L. G. M.
Pettersson
,
J. Electron Spectrosc. Relat. Phenom.
177
,
135
157
(
2010
).
101.
P.
Virtanen
,
R.
Gommers
,
T. E.
Oliphant
,
M.
Haberland
,
T.
Reddy
,
D.
Cournapeau
,
E.
Burovski
,
P.
Peterson
,
W.
Weckesser
,
J.
Bright
,
S. J.
van der Walt
,
M.
Brett
,
J.
Wilson
,
K. J.
Millman
,
N.
Mayorov
,
A. R. J.
Nelson
,
E.
Jones
,
R.
Kern
,
E.
Larson
,
C. J.
Carey
,
İ.
Polat
,
Y.
Feng
,
E. W.
Moore
,
J.
VanderPlas
,
D.
Laxalde
,
J.
Perktold
,
R.
Cimrman
,
I.
Henriksen
,
E. A.
Quintero
,
C. R.
Harris
,
A. M.
Archibald
,
A. H.
Ribeiro
,
F.
Pedregosa
,
P.
van Mulbregt
,
A.
Vijaykumar
,
A. P.
Bardelli
,
A.
Rothberg
,
A.
Hilboll
,
A.
Kloeckner
,
A.
Scopatz
,
A.
Lee
,
A.
Rokem
,
C. N.
Woods
,
C.
Fulton
,
C.
Masson
,
C.
Häggström
,
C.
Fitzgerald
,
D. A.
Nicholson
,
D. R.
Hagen
,
D. V.
Pasechnik
,
E.
Olivetti
,
E.
Martin
,
E.
Wieser
,
F.
Silva
,
F.
Lenders
,
F.
Wilhelm
,
G.
Young
,
G. A.
Price
,
G.-L.
Ingold
,
G. E.
Allen
,
G. R.
Lee
,
H.
Audren
,
I.
Probst
,
J. P.
Dietrich
,
J.
Silterra
,
J. T.
Webber
,
J.
Slavič
,
J.
Nothman
,
J.
Buchner
,
J.
Kulick
,
J. L.
Schönberger
,
J. V.
de Miranda Cardoso
,
J.
Reimer
,
J.
Harrington
,
J. L. C.
Rodríguez
,
J.
Nunez-Iglesias
,
J.
Kuczynski
,
K.
Tritz
,
M.
Thoma
,
M.
Newville
,
M.
Kümmerer
,
M.
Bolingbroke
,
M.
Tartre
,
M.
Pak
,
N. J.
Smith
,
N.
Nowaczyk
,
N.
Shebanov
,
O.
Pavlyk
,
P. A.
Brodtkorb
,
P.
Lee
,
R. T.
McGibbon
,
R.
Feldbauer
,
S.
Lewis
,
S.
Tygier
,
S.
Sievert
,
S.
Vigna
,
S.
Peterson
,
S.
More
,
T.
Pudlik
,
T.
Oshima
,
T. J.
Pingel
,
T. P.
Robitaille
,
T.
Spura
,
T. R.
Jones
,
T.
Cera
,
T.
Leslie
,
T.
Zito
,
T.
Krauss
,
U.
Upadhyay
,
Y. O.
Halchenko
, and
Y.
Vázquez-Baeza
,
Nat. Methods
17
,
261
(
2020
).
102.
O.
Karis
,
A.
Nilsson
,
M.
Weinelt
,
T.
Wiell
,
C.
Puglia
,
N.
Wassdahl
,
N.
Mårtensson
,
M.
Samant
, and
J.
Stöhr
,
Phys. Rev. Lett.
76
,
1380
1383
(
1996
).
103.
C.
Keller
,
M.
Stichler
,
G.
Comelli
,
F.
Esch
,
S.
Lizzit
,
W.
Wurth
, and
D.
Menzel
,
Phys. Rev. Lett.
80
,
1774
1777
(
1998
).
104.
L.
Gross
,
N.
Moll
,
F.
Mohn
,
A.
Curioni
,
G.
Meyer
,
F.
Hanke
, and
M.
Persson
,
Phys. Rev. Lett.
107
,
086101
(
2011
).
105.
A. G.
Borisov
,
J. P.
Gauyacq
,
E. V.
Chulkov
,
V. M.
Silkin
, and
P. M.
Echenique
,
Phys. Rev. B
65
,
235434
(
2002
).
106.
J.
Repp
,
G.
Meyer
,
S. M.
Stojković
,
A.
Gourdon
, and
C.
Joachim
,
Phys. Rev. Lett.
94
,
026803
(
2005
).
107.
W.
Dou
and
J. E.
Subotnik
,
J. Chem. Phys.
145
,
054102
(
2016
).
108.
W.
Dou
and
J. E.
Subotnik
,
J. Chem. Phys.
146
,
092304
(
2017
).
109.
W.
Dou
and
J. E.
Subotnik
,
J. Chem. Phys.
148
,
230901
(
2018
).
110.
A.
Cao
,
V. J.
Bukas
,
V.
Shadravan
,
Z.
Wang
,
H.
Li
,
J.
Kibsgaard
,
I.
Chorkendorff
, and
J. K.
Nørskov
,
Nat. Commun.
13
,
2382
(
2022
).
111.
S.
Vijay
,
G.
Kastlunger
,
K.
Chan
, and
J. K.
Nørskov
,
J. Chem. Phys.
156
,
231102
(
2022
).
112.
S.
Bhattacharjee
,
U. V.
Waghmare
, and
S.-C.
Lee
,
Sci. Rep.
6
,
35916
(
2016
).
113.
P.
Chen
,
D.
Fan
,
A.
Selloni
,
E. A.
Carter
,
C. B.
Arnold
,
Y.
Zhang
,
A. S.
Gross
,
J. R.
Chelikowsky
, and
N.
Yao
,
Nat. Commun.
14
,
1460
(
2023
).
114.
G.
Mandi
and
K.
Palotas
,
Appl. Surf. Sci.
304
,
65
72
(
2014
).
115.
H.-N.
Xia
,
E.
Minamitani
,
R.
Žitko
,
Z.-Y.
Liu
,
X.
Liao
,
M.
Cai
,
Z.-H.
Ling
,
W.-H.
Zhang
,
S.
Klyatskaya
,
M.
Ruben
, and
Y.-S.
Fu
,
Nat. Commun.
13
,
6388
(
2022
).
Published open access through an agreement with Technische Informationsbibliothek

Supplementary Material