Nonequilibrium Green's functions and atom-surface dynamics: Simple views from a simple model system

We employ Non-equilibrium Green's functions (NEGF) to describe the real-time dynamics of an adsorbate-surface model system exposed to ultrafast laser pulses. For a finite number of electronic orbitals, the system is solved exactly and within different levels of approximation. Specifically i) the full exact quantum mechanical solution for electron and nuclear degrees of freedom is used to benchmark ii) the Ehrenfest approximation (EA) for the nuclei, with the electron dynamics still treated exactly. Then, using the EA, electronic correlations are treated with NEGF within iii) 2nd Born and with iv) a recently introduced hybrid scheme, which mixes 2nd Born self-energies with non-perturbative, local exchange-correlation potentials of Density Functional Theory (DFT). Finally, the effect of a semi-infinite substrate is considered: we observe that a macroscopic number of de-excitation channels can hinder desorption. While very preliminary in character and based on a simple and rather specific model system, our results clearly illustrate the large potential of NEGF to investigate atomic desorption, and more generally, the non equilibrium dynamics of material surfaces subject to ultrafast laser fields.


Introduction
Since the advent of femtosecond (fs) spectroscopy, the experimental characterisation of ultrafast chemical reactions has grown spectacularly [1]. At the same time, with the development of ever Figure 1. The model system for adsorbate-dynamics investigated in this work, for a finite (5-site) and infinite 1D substrate (a bulk substrate is also shown). The adsorbate is coupled to the substrate via an hopping term V (x) which depends on the adsorbate-substrate distance x (a schematic effective adsorbate potential is also shown). The interaction U among electrons is present only in the lower valence level of the adsorbate (see Eq. (1) and afterwards for details). more sophisticated techniques to probe and control atomic and molecular processes far away from equilibrium, there has also been substantial progress in the theoretical description of the ultrafast dynamics of atoms and molecules in the gas phase [2,1,3,4,5]. At present increasing attention is directed towards the study of surface-adsorbate complexes, since these systems are key to understanding important chemical processes, [6,7,8,9], and also constitute the next natural paradigm for fundamental light-matter processes.
However, for surface-adsorbate systems the level of theoretical understanding is at the moment not on the level comparable with free atoms and molecules [10,11,12,13]. Even though there have been recent steps in this direction [14,15,16,17,18,19] using sophisticated methods such as Time-Dependent Density Functional Theory (TDDFT) [20,21,22] and Non-Equilibrium Green's Functions (NEGF) [23,24,25,13], the approximate treatment of particle interactions used in practice are not sufficient to accurately describe even the initial stages of desorption. Thus, at present, there is motivation to consider model approaches to ultrafast surface-adsorbate dynamics, that capture the highly non-perturbative situation of correlated electron-nuclear motion and that can be used to gain insight into the outcome of laser induced excitations [26,27].
Here, we perform an exploratory application of NEGF to adsorbate dynamics, an important physical mechanism at surfaces out of equilibrium. To this end, we apply NEGF to a simplified version of a recently introduced model for adsorbate dynamics induced by ultrashort laser pulses [28]. If solved exactly, the model can deal in a simplified way but on equal footing with electronelectron and electron-nuclear interactions, plasmon screening, real-space nuclear dynamics and core-hole relaxation. Thus it permits to study in a coherent way many competing surface response mechanisms and time scales. Here, for simplicity, core levels in the adsorbate and the plasmon response from the substrate are excluded from the outset.
To solve the model exactly, we must limit ourselves to a finite number of substrate sites, or equivalently a limited number of dissipation channels. As an alternative and a way to overcome this obvious limitation, we also consider a semi-infinite substrate within the framework NEGF, where however approximate treatment of interactions must be introduced. For this approach the finite-size version of our model provides an exact benchmark. To disentangle the effects of our various approximations, we perform them in sequence, starting from the numerically exact solution. The first step is to treat the nuclear dynamics in the Ehrenfest approximation while retaining the exact treatment of the electrons. In the second step the electron dynamics is treated within perturbation theory at the level of the Second-Born approximation (2BA). When both approximations are being employed we are free to extend the substrate to become semi-infinite by standard embedding procedures [29].
It is fair to say that a perturbative treatment of the NEGF (in our case the 2BA) may fail to describe ultrafast real-time dynamics of electrons, especially when dealing with strongly electron correlations. For such systems an accurate and computationally viable first-principles description like NEGF or TDDFT is currently lacking. They both rely on key ingredients [the exchange-correlation (XC) potential for TDDFT and the self-energy Σ for the NEGF] that in general are only approximately known and may be not sufficient in the strongly-correlated fast regimes. Recently, we suggested a step towards the solution to this problem using the strengths of both methods, by mixing a non-perturbative adiabatic-TDDFT XC potential with perturbative-based self-energy of NEGF [30]. Here, we test this hybrid NEGF/TDDFT method at the 2BA level and compare it to the pure NEGF 2BA in our simple model for adsorbate dynamics.
The rest of this paper is structured as follows: the model for adsorbate dynamics is introduced in Sect. 2; the different methods of solutions are described in Sect. 3; the practical details are mentioned in Sect. 4; the results and their discussion are presented in Sect. 5. Some conclusive remarks and outlooks are given in Sect. 6.

Model Hamiltonian for adsorbate dynamics
The model we use is a simplified version of the one introduced in Ref. [28]. It consists of a rigid atom chain (the substrate) connected to a mobile adsorbate (Fig. 1). The Hamiltonian iŝ whereĤ a describes the adsorbate,Ĥ s the substrate andĤ as the adsorbate-substrate interaction. To bring the system out of equilibrium a laser fieldΛ(t) is applied. The adsorbate part H a is wherep is the momentum of a mobile atom with mass M . The operator a † v,σ creates an electron with spin σ at the valence orbital v of the adsorbate with energy v , andn v,σ = a † v,σ a v,σ . In the adsorbate level v 1 , valence electrons mutually interact with interaction strength U . The substrate is modelled as a non-interacting tight-binding chain with one orbital per site: where c † R,σ creates an electron at site R with spin σ. In the following, the hopping amplitude V between neighbouring sites R and R is taken as the energy unit, i.e. V = 1. The adsorbatesubstrate interaction Hamiltonian H as is given by where we denote the surface site of the substrate by S and x is the adsorbate-surface distance. The first term gives a repulsive interaction between the atomic cores, while the second term is attractive and due to the exchange of electrons between the adsorbate and the surface. We consider κ and λ as phenomenological parameters that can be adjusted to give reasonable values for the binding energy E b , vibrational frequency ω ph and effective hopping amplitude V = g e −λ(x−1) . The shape of the fieldΛ(t) in Eq.(1) depends on the experiment considered, and is switched on at t = 0 so thatΛ = 0 for t ≤ 0. In the dipole approximation, and for a field coupling only the adsorbate valence levels, is an arbitrary function of time with a vanishing integral.

Exact and approximate treatments 3.1. Finite substrate: Exact wavefunction solution
An exact solution of the Schrödinger equation for the model presented above is viable only for a short substrate, since the size S of the Hilbert space H grows exponentially with the number of orbitals. In more detail, we note that H = H el ⊗ H n , with H el the electronic and H n the atomic (nuclear) Hilbert spaces, respectively. Due the explicit form of H in Eq. (1), the size of S n grows linearly, while, for a spin-compensated system with L electronic orbitals, we find an exponential scaling of S. With and a space-grid subspace for the nucleus of size S n 500, substrate lengths L 6 can be comfortably afforded.
Specifically, the ground state |g of the model is then determined with exact diagonalization (ED) , by expanding the state vector in the basis {|n iσ , x k }, where n iσ is the occupation on site i of electrons with spin σ and x k is the k:th mesh-point on a uniform grid in the interval [0, x max ]. Assuming that the system is in the state |g at t = 0, we then propagate the Schrödinger equation for t > 0 using the short iterated Lanczos algorithm.

Finite and semi-infinite substrate: NEGF
If we are not interested in the full wavefunction of the system but we choose to directly look at single particle quantities, NEGF offer a natural general and in-principle-exact way to go beyond size limitations, and to treat electrons and nuclei exactly and on equal footing (in practice approximations must be introduced). For our present system, this would imply a twocomponent formulation of coupled Kadanoff-Baym (KBE) equations, one component for each type of degree of freedom [19,25]. However, to just explore the scope of NEGF for adsorbate dynamics, we here adopt a much cruder strategy, namely we use NEGF for electrons and a classical treatment for the nuclei. The most straightforward approximation is then to use the Born-Oppenheimer separation of electronic and atomic degrees of freedom, and for each atomic configuration extract effective potentials generated by the electrons and felt by the atoms, known as Born-Oppenheimer surfaces. Given these potential surfaces, there exists a number of schemes to propagate an initial atomic wave packet; here we consider the Ehrenfest approximation (EA) [31], which amounts to replace the nuclear wavepacket with a delta function, and to let the atomic position evolve under the average force from the electrons, or in other words demoting the operatorsx andp to classical variables 1 . For the model in Eq. (1) the mixed quantum-classical equations in the EA are Once reformulated in the language of NEGF, these become where the single particle Green's function is defined via the time contour ordering T as for a convenient choice of basis {ϕ i }, here taken as the site basis of our model. In the coupled equations above, Eqs. (7) governs the nuclear dynamics, whilst Eq. (8) is the standard Kadanoff-Baym equation for G, where dependence of the onebody Hamiltonian matrix h ij on x is shown, and an implicit dependence of Σ and G on x is understood.
If the system is composed of an interacting region connected to a non-interacting reservoir (as for example, the adsorbate and substrate parts in our system), the latter can be removed from the explicit description and still taken into account exactly via an embedding procedure [29]. In this case, the self-energy for the interacting region becomes Σ ij (z 1 , z 2 ) = Σ em ij (z 1 , z 2 ) + Σ M B ij (z 1 , z 2 ). The first term is the embedding self-energy that gives an exact coupling to the reservoir, and the second term is the many-body self-energy that contain all interactions between electrons in the active region. After finding the initial G ij (t 0 , t 0 ) on the imaginary part of the Keldysh contour, the KBE are propagated self-consistently within the same approximation for Σ ij [23]. With NEGF and the embedding self-energies, a much larger number of orbitals can be explicitly taken into account (plus the levels of the macroscopic reservoir) than e.g. with exact diagonalization, and in any dimensionality. The downside is that interactions now have in general to be treated in an approximate way, e.g. perturbatively: Σ M B ij (z 1 , z 2 ) ≈ Σ P T ij (z 1 , z 2 ).

A NEGF/TDDFT hybrid scheme
In our numerical calculations for adsorbate dynamics, correlations are treated in the 2nd Born approximation (2BA), Σ P T ≡ Σ 2BA . However, starting from this approximation, we also consider a second approach, recently introduced in Ref. [30], which adds to the Kadanoff-Baym equation an exchange-correlation potential v xc,N P of DFT which is nonperturbative in character. At the same time, we subtract a perturbative exchange-correlation potential v xc,P T , to avoid double counting of the interactions. The approach is conserving in the Kadanoff-Baym sense, and has the same numerical complexity of standard KBE. In this hybrid NEGF/TDDFT scheme, the modified Kadanoff-Baym equation reads (repeated indices are implicitly summed over) Both exchange-correlation potentials are taken in adiabatic local density approximation (ALDA) [32], i.e. the time-dependent XC-potential is approximated as a potential coming from a reference equilibrium system with the same density v xc i (z 1 ) ≈ v xc ref.
(n i (z 1 )). The reference system is chosen according to the model to be solved, as discussed in the next section.

Additional aspects of the NEGF/TDDFT scheme
In the NEGF/TDDFT method, we need suitable XC potentials v xc,N P and v xc,P T for the lower valence level. The geometrical structure of the model in Fig. 1 does not allow for a direct application of standard 1D ALDA (which is obtained from a 1D homogeneous Hubbard model where each site has a bond on the left and one on the right). However, both in equilibrium and for time-varying perturbations acting only on the adsorbate, symmetry arguments permit to re-interpret the adsorbate as a two-orbital central site connected to the parity-even states (but not to odd ones) of two identical 1D semi-infinite substrates, with a reduced site-substrate(s) coupling V sym = V / √ 2. For the lower valence level v 1 of the central site in the symmetric model, we employ the ALDA potential v xc ≈ v xc ref. n ν 1 , U/V sym using the homogeneous 1D Hubbard model as the reference equilibrium system [33,34,35]. The corresponding potentials v xc,N P and v xc,P T are shown in Fig. 2: they mutually differ in many ways in their behavior, the key difference being that v xc,N P (v xc,P T ) has (has not) a discontinuity at half-filling.

Practical details 4.1. On the different treatments
In the actual calculations, we consider both finite (L = 5) and semi-infinite 1D substrates (see Eq. 1 and Fig. 1). For L = 5, the model is solved exactly with a full-quantum mechanical solution for the electrons and the nucleus (Sect. 3.1) and, in the EA, with different degrees of approximations for the electron dynamics (Sects. 3.2, 3.3). In the EA for the nuclei plus the exact dynamics of the electrons (EA+Ex e ), the electronic ground state |g e and nuclear equilibrium position x eq are found at the stationary points of Eq. (5-6), after which we propagate Eq. (5-6) in time using the Lanczos technique for the electrons and a velocity Verlet-type algorithm for the nuclei. When the EA is used in conjunction with NEGF, the electronic correlations are taken into account at the 2BA level (EA+2BA), or via the hybrid NEGF/TDDFT scheme (EA+Hyb). The ground state Green's function G ij (t 0 , t 0 ) and equilibrium position x eq are now found solving self-consistently Eq.  is shown to the left, and to the right is the corresponding perturbative potential v xc,P T ref.
in the 2BA. We plot the potentials as a function of uniform density n and the ratio of interaction and hopping U t for fixed interaction U = 4.
Within the NEGF treatments, we also consider a semi-infinite 1D substrate via an embedding self energy [29], with the solution method otherwise unchanged.

On the external fields
We used two time-dependent fields to perturb the system: Λ(t) ≡ Λ a (t) = Ae −(t−t 0 ) 2 /2σ 2 and Λ(t) ≡ Λ b (t) = Λ a (t) sin (ωt + π/2), with amplitude either A = 2 or 4. The electronic parameters of the model are chosen as V = 1, 1 = 1, 2 = 3, thus giving a cationic adsorbate in the ground state, and we consider the interaction strengths U = 4, 6 and 8. For the nucleus, κ = 2.4, g = 1.8 and λ = 1.2, giving a binding energy E b 1.5, an harmonic frequency ω ph 0.2 and an effective hopping amplitude V 1.8. These values are typical for chemisorption. The Gaussian pulse Λ a is an artificial perturbation and serves as a tool to investigate slower dynamics, while the sinusoidally modulated pulse Λ b represents a realistic pulse with a carrier wavelengh of 400nm (for a hopping amplitude V 2 eV the unit of time is 1/V 1.3 fs, which is the period of 400 nm radiation). Both perturbations extend from t = 0 to t = 8, with time measured in units of the inverse hopping 1/V . In this time interval the evolution is dominated by the perturbation and all the methods give basically the same result. For later times the dynamics is dominated by the memory of the system, and the solutions differ depending on the approximation used.

Exact vs Ehrenfest dynamics for the nuclei
Let us first look at the evolution of the nuclear subsystem, as shown in Fig. 3. With the full quantum mechanical solution we observe bond stretching for the pulse Λ a , but for large times the nuclear position returns to a region close to its equilibrium value. For the pulse Λ b the nuclear coordinate continues to increase also for large times, signalling a finite probability of desorption. This is typical of the full quantum solution, where in general one part of the initial wave packet remains bound and performs oscillations in the potential well, while another part overcomes the potential barrier and thereafter propagates quasi-freely (strictly the potential is zero only for x → ∞). We note that with increasing interaction U and increasing amplitude A, the desorption probability also increases. The former is a result of the interaction induced quenching of electronic motion between the adsorbate and surface, resulting in a lower average bond kinetic energy and therefore a weaker bond (cf. the form of electron-nuclear coupling in Eq. 4 and the definition K as = c † v 1 c S + c † v 2 c S + h.c. for the bond kinetic energy). Since the EA is a classical treatment of nuclei the only possible outcomes of an experiment is to have either full or no desorption, and in the intermediate regime the approximation will inevitably fail to reproduce the exact dynamics. This is verified in Fig. 3 where the exact results are best reproduced in the cases where the quantum solution predicts either a very small or very large desorption probability. The agreement is in general better for the pulse Λ a , which varies more slowly in time and is closer to the adiabatic situations where the EA is known to perform well.

Exact vs NEGF electronic correlations within the Ehrenfest dynamics.
We now turn our attention to the electron dynamics, as depicted in Fig. 4. When comparing the exact solution with the EA+Ex e dynamics we see clear differences, that are most pronounced in the case of partial desorption. In particular the EA shows oscillations in the electron density of greater amplitude, that persist even for long time. This can be explained in part by the inability of the EA to transfer energy from the electronic to the nuclear subsystem, due to neglect of electron-nuclear correlations, that in the exact simulation act as a dissipation channel for the electrons. In part it is also due to the fact that the effective adsorbate-surface hopping amplitude is larger in the EA, making is more probable for the electrons to tunnel between the substrate and the adsorbate. With quantum mechanical nuclei, the hopping amplitude (proportional to e −λx , see Eq. 4) is decreased by the presence of a split-off part of the nuclear wave packet having propagated to large distances.
Comparing the exact and the EA+2B electron evolution we observe that the oscillations in density seen in the EA are reduced compared to the EA+Ex e result, and that this difference increases with stronger interactions. The damping of electron oscillations is a generic feature of approximate solutions of the KBE in finite systems, and was first discussed in [38]. Briefly, the reason for damping is that self-consistent approximations (like the 2B) include infinitely many electron-hole excitations (through the infinite set of diagrams contained in the self-energy), which cannot occur in finite systems. Thus, even though the agreement between the exact and the EA+2B solution seems better than EA+Ex e , it should be kept in mind that this result is due to unphysical effects.
With the EA+Hyb method we observe a slightly better agreement in phase with the EA+Ex e treatment, but for the model considered here the improvement of local correlations brought by the exchange-correlation potential have a rather small overall effect. This points to the importance of non-local correlations for the desorption process, and to capture more closely the exact dynamics approximations beyond 2B (or even inclusion of non-perturbative effects) are necessary. As a first step in this direction we are currently working on including correlations in the T -matrix approximation.

The effect of a semi-infinite substrate
Finally, we contact the system to a semi-infinite lead and compare the resulting dynamics within EA+2B and EA+Hyb to the corresponding finite system simulations (see Fig. 5). Interestingly, the addition of the particle reservoir seems to have a rather small impact on the short-time dynamics, both for the electronic and nuclear subsystem. At larger times the desorption probability decreases in the presence of a lead, since there is now an additional energy dissipation channel for the electrons that reduce the amount of energy transferred to the nuclei. However, for the nuclear part there is not much difference between the EA+2B and EA+Hyb methods. The electron dynamics of the finite and extended systems show a qualitative agreement, the major difference being that the amplitude of the density oscillations decay faster in the latter case. For a contacted system the damping of electron oscillations is an inherent physical effect, in contrast to the artificial behavior discussed in the previous section. For the faster pulse Λ b we also observe a slight dephasing of the EA+2B and EA+Hyb solutions, favoring the EA+Hyb method in both the finite and contacted case.
The difference between the EA+2B (or EA+Hyb) and the exact solution is typically greater than the difference between the finite and contacted solution within a given approximation scheme. This might imply that to model adsorbate-surface systems, and in particular nonperturbative processes such as desorption, it is more important to properly account for the electron-electron and electron-nuclear correlations close to the surface, than for the infinite extent of the system. This conclusion is strongly dependent on the parameters chosen, but for a strong electron-nuclear coupling (so that the system is in the surface molecule limit) we expect it to hold.  Figure 5. Time-dependent densities of the lower valence level (top) and interatomic distance (bottom), for a finite chain in the exact quantum treatment (red) and the approximate schemes EA+2B (blue) and EA+Hyb (orange), as well as for a contacted system within EA+2B (light blue) and EA+Hyb (light orange). The first column shows results for the external field Λ a (t) of amplitude A = 4.0 with the interaction strength U = 4, and the second column gives the results for the external field Λ b (t) of amplitude A = 2.0 with the interaction strength U = 6.

Conclusions
We present a model to treat real-time surface-adsorbate system dynamics induced by ultrafast laser pulses, and solve it exactly numerically in its finite size realisation. The exact results are compared to a treatment within the framework of NEGF, on the level of the second Born approximation and also to a recently introduced hybrid NEGF/TDDFT approach. By the standard embedding procedure this allows us to also solve the semi-infinite version of our model. In the finite size case we treat the nuclear coordinate exactly using first quantisation, and compare it to the semi-classical case of having exact electronic dynamics and nuclear motion within the Ehrenfest approximation. For the process of desorption, interesting due to its non-perturbative character, we find good qualitative agreement between the exact and Ehrenfest nuclear dynamics when the exact result predicts either a very small or a very large desorption probability. In the intermediate case, characterized by a splitting of the nuclear wave packet into a bound and a quasi-free part of comparable magnitude, the Ehrenfest treatment fails to reproduce the exact results. For the case of chemisorption considered here, we further find that the inclusion of a semi-infinite substrate have a limited effect but predicts a slightly lower desorption probability. Our results also indicate the importance of non-local correlation effects beyond second Born.