Charge dynamics in molecular junctions: Nonequilibrium Green's Function approach made fast
S. Latini, E. Perfetto, A.-M. Uimonen, R. van Leeuwen, G. Stefanucci
aa r X i v : . [ c ond - m a t . m e s - h a ll ] N ov Charge dynamics in molecular junctions: Nonequilibrium Green’s Function approachmade fast
S. Latini, E. Perfetto, A.-M. Uimonen, R. van Leeuwen,
2, 3 and G. Stefanucci
1, 4, 3 Dipartimento di Fisica, Universit`a di Roma Tor Vergata,Via della Ricerca Scientifica 1, 00133 Rome, Italy Department of Physics, Nanoscience Center, FIN 40014, University of Jyv¨askyl¨a, Jyv¨askyl¨a, Finland European Theoretical Spectroscopy Facility (ETSF) INFN, Laboratori Nazionali di Frascati, Via E. Fermi 40, 00044 Frascati, Italy
Real-time Green’s function simulations of molecular junctions (open quantum systems) are typi-cally performed by solving the Kadanoff-Baym equations (KBE). The KBE, however, impose a se-rious limitation on the maximum propagation time due to the large memory storage needed. In thiswork we propose a simplified Green’s function approach based on the Generalized Kadanoff-BaymAnsatz (GKBA) to overcome the KBE limitation on time, significantly speed up the calculations,and yet stay close to the KBE results. This is achieved through a twofold advance: first we showhow to make the GKBA work in open systems and then construct a suitable quasi-particle prop-agator that includes correlation effects in a diagrammatic fashion. We also provide evidence thatour GKBA scheme, although already in good agreement with the KBE approach, can be furtherimproved without increasing the computational cost.
PACS numbers: 05.60.Gg,05.10.-a,73.63.-b,72.10.Bg
I. INTRODUCTION
Charge transfer through nanoscale interfaces is anubiquitous dynamical process in molecular electronics,photovoltaics, electroluminiscence and transient spec-troscopy, to mention a few emerging fields of research.
The complexity of the molecules (or molecular aggregate)and of the contacts to a source/drain electrode, as wellas the simultaneous interplay of Coulomb repulsion andvibrational effects make these research fields an interdis-ciplinary topic where physics, chemistry and engineeringmeet. Reliable theoretical predictions require an accuratedescription of the nuclear degrees of freedom, a carefulselection of the electronic basis functions and a propertreatment of correlation effects.Among the ab initio methods, Density FunctionalTheory (DFT) and its Time Dependent extension (TDDFT) stand out for the advantageous scaling of thecomputational cost with increasing the system size andthe propagation time. However, as for any other method,a (TD)DFT implementation is based on some approxima-tion and, at present, the available approximations are in-adequate to capture correlation effects like the Coulombblockade or the polarization-induced renormalizationof the molecular levels.
These effects are particularlyimportant in a donor-acceptor complex, in a molecularjunction in the weak-coupling regime and more generallywhen the transition rate for an electron to move fromone atom to another is small. Many-body approachesbased on Nonequilibrium Green’s Functions (NEGF)offer a promising alternative as the relevant scatteringprocesses to describe the aforementioned effects can beincorporated either through a proper selection of Feyn-man diagrams or through a decoupling scheme for thehigher order Green’s functions. Real-time simulationswithin the NEGF are performed by solving the Kadanoff- Baym equations (KBE), which are a set of cou-pled nonlinear integro-differential equations for the one-particle Green’s function. Unfortunately, the price topay in solving the KBE is that the computational timescales cubically with the propagation time (given the self-energy) whereas in TDDFT the scaling is linear (giventhe exchange-correlation potential).In the mid eighties Lipavsky et al. proposed anapproximation to scale down the computational time(from cubic to quadratic) of the KBE. This approxima-tion is known as the Generalized Kadanoff-Baym Ansatz(GKBA) and has been successfully applied to strongly in-teracting nuclear matter, electron plasma, , carrierdynamics of semiconductors, , optical absorptionspectra, quasi-particle spectra, and more recently ex-cited Hubbard clusters. In all these cases the systemis either a bulk periodic system or a finite system. It iscurrently unknown how the GKBA performs for nanos-tructures chemically bonded to or adsorbed on a surface(open system). In fact, in open systems a number ofissues have to be addressed before a GKBA calculationcan be carried out. For instance the GKBA remains anapproximation even in a noninteracting (or mean-field)treatment whereas in closed systems it is exact. Further-more the performance of the GKBA strongly depends onthe quality of the quasi-particle propagator and, as weshall see, in open systems the available approximationsperform rather poorly.This work contains a thorough study of the GKBA inopen systems. In Section II we derive the fundamen-tal equations and present a few exact properties. Herethe discussion is mainly focussed on noninteracting andmean-field electrons. Important aspects of the GKBAlike the construction of a mean-field propagator as well asissues related to relaxation and local thermalization areanalyzed and addressed. This preliminary investigationis particularly relevant since, as previously mentioned,the GKBA is an approximation already at the mean-fieldlevel. In the correlated case the GKBA simulations us-ing a mean-field propagator are far off the KBE results.In Section III we propose a couple of correlated propa-gators to remedy this deficiency. Our propagators havethe merit of scaling quadratically with the propagationtime and hence the computational gain of the GKBA ismaintained. The different GKBA schemes are comparedwith the full KBE approach in Section IV. We considertwo systems, a molecular junction under applied bias anda donor-acceptor complex under illumination, and calcu-late local currents and densities. Both systems constitutea severe test for the GKBA as the inclusion of correla-tions changes dramatically the mean-field picture. Theimportant message emerging from this study is that oneof the proposed GKBA schemes is in fairly, sometimesextremely, good agreement with the KBE approach. Wealso provide numerical evidence that the GKBA schemecan be further improved at the same computational cost.In conclusion, time-dependent simulations of open sys-tems within the NEGF framework can be made muchfaster.
II. GKBA IN OPEN SYSTEMS
In this Section we briefly review the KBE for opensystems and discuss in detail the simplifications broughtabout by the GKBA. The most general Hamiltonianwhich describes a molecular junction in contact with M electronic reservoirs has the formˆ H = M X α =1 ˆ H α + ˆ H J + ˆ H T . (1)In Eq. (1) the Hamiltonian of the α reservoir readsˆ H α = X kασ ǫ kα ˆ d † kασ ˆ d kασ (2)with ˆ d kασ the annihilation operator for electrons of spin σ and energy ǫ kα . The Hamiltonian of the molecularjunction is expressed in terms of the operators ˆ d iσ forelectrons of spin σ in the i -th localized molecular orbitalˆ H J = X ijσ h ij ˆ d † iσ ˆ d jσ + 12 X ijmnσσ ′ v ijmn ˆ d † iσ ˆ d † jσ ′ ˆ d mσ ′ ˆ d nσ (3)where h ij are the one-electron matrix elements of the one-body part (kinetic plus potential energy) and v ijmn arethe two-electron Coulomb integrals. The last term in Eq.(1) is the tunneling Hamiltonian between the differentsubsystems and readsˆ H T = X kασ X i (cid:16) T kα,i ˆ d † kασ ˆ d iσ + H . c . (cid:17) (4) with T kα,i the tunneling amplitude between the i -th stateof the molecular junction and the k state of the α reser-voir.Initially, say at time t = 0, the system is in equilib-rium at inverse temperature β and chemical potential µ .We assume that this equilibrium state can be reachedstarting from the uncontacted ( T kα,i = 0) and noninter-acting ( v ijmn = 0) system in the remote past, t = −∞ ,and then propagating forward in time with the full in-teracting and contacted Hamiltonian until t = 0. Thisamounts to assume that initial-correlation and memoryeffects are washed out. In our experience this assump-tion is always verified. At time t = 0 the systemis driven out of equilibrium by external electromagneticfields, ǫ kα → ǫ kα + V α ( t ) and h ij → h ij ( t ). We are in-terested in monitoring the evolution of the electronic de-grees of freedom through the calculation of observablequantities like, e.g., the local occupation and current. A. Green’s function and KBE
The building block of any diagrammatic many-bodyapproach is the Green’s function defined according to G ij ( z, z ′ ) = 1 i hT n ˆ d iσ,H ( z ) ˆ d † jσ,H ( z ′ ) o i . (5)In this definition the symbol “ h . . . i ” denotes a gran-canonical average, and T is the contour ordering act-ing on operators in the Heisenberg picture. The Green’sfunction has arguments z and z ′ on the contour γ goingfrom −∞ to ∞ (forward branch) and back from ∞ to −∞ (backward branch). On this contour G ij satisfiesthe equations of motion (in matrix form) (cid:20) i ddz − h HF ( z ) (cid:21) G ( z, z ′ ) = δ ( z, z ′ ) + Z γ d ¯ z Σ( z, ¯ z ) G (¯ z, z ′ )(6)and its adjoint. Let us describe the various quantitiesin this equation. The Hartree-Fock (HF) single-particleHamiltonian is the sum of h and the HF potential h HF ,ij = h ij + X mn (2 v imnj ρ nm − v imjn ρ nm ) (7)where ρ nm ( z ) ≡ − iG nm ( z, z + ) (8)is the time-dependent single-particle density matrix. Thekernel Σ = Σ em +Σ c is the sum of the so called embeddingself-energy and the correlation self-energy. The formercan be calculated directly from the parameters of theHamiltonian and readsΣ em ,ij ( z, z ′ ) = X kα T i,kα g kα ( z, z ′ ) T kα,j (9)where g kα ( z, z ′ ) = 1 i (cid:2) θ ( z, z ′ ) ¯ f ( ǫ kα ) − θ ( z ′ , z ) f ( ǫ kα ) (cid:3) e − iφ kα ( z,z ′ ) (10) Σ = + c ,ij i n m jp q r s G nm G pq G sr v i r pn v m q s j i n m jq s r p v irpn v mqsj G pm G nq G sr FIG. 1. Diagrams for the 2B correlation self-energy. is the Green’s function of the disconnected α reser-voir. In Eq. (10) f ( ǫ ) = 1 / ( e β ( ǫ − µ ) + 1) is the Fermifunction, ¯ f ( ǫ ) = 1 − f ( ǫ ) and the phase φ kα ( z, z ′ ) = R zz ′ d ¯ z ( ǫ kα + V α (¯ z )). The expression of the correlation self-energy depends on the choice of diagrams that we decideto include. In this work we consider the second-Born(2B) approximation which has been shown to produceresults very close to those of the GW approximation, and to those of numerically exact techniques in modelsystems. The 2B self-energy is given by the sum of thelowest order bubble diagram plus the second-order ex-change diagram, see Fig. 1, Σ c ,ij ( z, z ′ ) = X nmpqrs v irpn v mqsj × [2 G nm ( z, z ′ ) G pq ( z, z ′ ) G sr ( z ′ , z ) − G nq ( z, z ′ ) G sr ( z ′ , z ) G pm ( z, z ′ )] (11)To solve Eq. (6) we convert it into a set of coupledequations, known as the KBE, for real time (as opposedto contour time) quantities. This is done by letting z tovary on the forward (backward) branch and z ′ to vary onthe backward (forward) branch of the contour γ . Usingthe Langreth rules to convert contour-time convolu-tions into real-time convolutions we find (in matrix form) (cid:20) i ddt − h HF ( t ) (cid:21) G < ( t, t ′ ) = I < ( t, t ′ ) (12) G > ( t, t ′ ) " − i ←− ddt ′ − h HF ( t ′ ) = I > ( t, t ′ ) (13)with collision integrals I < ( t, t ′ ) = Z ∞−∞ d ¯ t (cid:2) Σ < ( t, ¯ t ) G A (¯ t, t ′ ) + Σ R ( t, ¯ t ) G < (¯ t, t ′ ) (cid:3) , (14) I > ( t, t ′ ) = Z ∞−∞ d ¯ t (cid:2) G > ( t, ¯ t )Σ A (¯ t, t ′ ) + G R ( t, ¯ t )Σ > (¯ t, t ′ ) (cid:3) . (15)Here the superscripts “ >, <, R , A” refer to the lesser,greater, retarded and advanced Keldysh components.Equations (12,13) are solved by a time stepping tech-nique, starting from a value G ≶ ( t in , t in ) at some initialtime t in < t and t ′ until a maximum propagation time t max . The time t in is chosen remotely enough in the past in order to havefull relaxation at t = 0, time at which the external fieldsare switched on. As I ≶ ( t, t ′ ) in Eqs. (14,15) involves integrals between t in (the self-energy vanishes for timessmaller than t in since the system is initially uncontactedand noninteracting) and either t or t ′ , the numerical ef-fort in solving the KBE scales like t . B. GKBA
The GKBA allows us to reduce drastically the com-putational time. The basic idea consists in obtaining aclosed equation for the equal time G < from which to cal-culate the time-dependent averages of all one-body ob-servables like, e.g., density, current, dipole moment, etc.The GKBA is therefore an ansatz for the density matrix ρ ( t ) = − iG < ( t, t ), not for the spectral function which hasto be approximated separately, see below.The exact equation for ρ ( t ) follows from the differencebetween Eq. (12) and its adjoint, and reads ddt ρ ( t ) + i [ h HF ( t ) , ρ ( t )] = − (cid:0) I < ( t, t ) + H . c . (cid:1) . (16)This is not a closed equation for ρ as the collision integralcontains the off-diagonal (in time) G ≶ . To close Eq. (16)we make the GKBA G < ( t, t ′ ) = iG R ( t, t ′ ) G < ( t ′ , t ′ ) − iG < ( t, t ) G A ( t, t ′ )= − G R ( t, t ′ ) ρ ( t ′ ) + ρ ( t ) G A ( t, t ′ ) (17)and similarly G > ( t, t ′ ) = G R ( t, t ′ )¯ ρ ( t ′ ) − ¯ ρ ( t ) G A ( t, t ′ ) (18)where ¯ ρ ( t ) = 1 − ρ ( t ) = iG > ( t, t ). However, the GKBAalone is not enough to close Eq. (16) since the quasi-particle propagator G R (and hence G A = [ G R ] † ), orequivalently the spectral function, remains unspecified.The possibility of using the GKBA in open systemsstrongly relies on the choice of G R . This is an importantpoint which we thoroughly address in the next Section.For the time being we observe that the numerical effortin solving Eq. (16) scales like t provided that thecalculation of G R does not scale faster.
1. Exact properties
Among the properties of the GKBA we mention thefulfillment of the relation G R − G A = G > − G < for anychoice of G R , and the fact that Eqs. (17,18) become anidentity in the limit t → t ′ since G R ( t + , t ) = − i . Anothervaluable feature (in systems out of equilibrium) is thatthe GKBA preserves the continuity equation. There is,however, an even more important property from whichthe physical contents of the GKBA become evident. Inclosed systems (Σ em = 0) and for HF electrons (Σ c = 0)the collision integrals vanish and Eqs. (17,18) are thesolution of Eqs. (12,13) provided that G R is the HFpropagator G R ( t, t ′ ) = − iθ ( t − t ′ ) T e − i R tt ′ d ¯ t h HF (¯ t ) (19)where T is the time-ordering operator. Therefore, themore the quasi-particle picture is valid the more theGKBA is accurate. A more exhaustive discussion on therange of applicability of the GKBA in closed systems canbe found in Refs. 19 and 39.In open systems the GKBA is not the solution of theHF equations since Σ em = 0 and hence the collision in-tegral is nonvanishing. The reliability of the GKBA inopen systems needs to be investigated already at the HFlevel. In HF the collision integrals are evaluated withΣ = Σ em and G R being the solution of (cid:20) i ddt − h HF ( t ) (cid:21) G R ( t, t ′ ) = δ ( t, t ′ )+ Z d ¯ t Σ Rem ( t, ¯ t ) G R (¯ t, t ′ ) . (20)In HF-GKBA the collision integrals are evaluated withΣ = Σ em , G < (¯ t, t ′ ) = ρ (¯ t ) G A (¯ t, t ′ ) and G A = [ G R ] † somesuitable propagator. If we calculate G R from Eq. (20)then the numerical advantage of the GKBA is lost sincethe computational cost of solving this equation scales like t . Thus the questions are: can a “computationallycheap” propagator be constructed for open systems? Ifso, how accurate is the solution of the HF-GKBA equa-tion?To answer these questions we consider a WideBand Limit (WBL) embedding self-energy Σ Rem ( t, t ′ ) = − ( i/ δ ( t − t ′ ) where Γ is a positive-semidefinite self-adjoint matrix. In this case the solution of Eq. (20) is G R ( t, t ′ ) = − iθ ( t − t ′ ) T e − i R tt ′ d ¯ t ( h HF (¯ t ) − i Γ / (21)which has the same mathematical structure of Eq. (19).In particular it has the group property G R ( t + δ, t ′ ) = iG R ( t + δ, t ) G R ( t, t ′ ) (22)and hence the number of operations to calculate G R forall t < t max and t ′ < t scales like t . The HF collisionintegral reads I < ( t, t ) = Z ∞−∞ d ¯ t Σ < em ( t, ¯ t ) G A (¯ t, t ) − i G < ( t, t ) , (23)whereas the HF-GKBA collision integral reads I < ( t, t ) = Z ∞−∞ d ¯ t Σ < em ( t, ¯ t ) G A (¯ t, t ) − i ρ ( t ) G A ( t − , t ) . (24)If in Eq. (24) we use for G A = [ G R ] † the HF re-sult in Eq. (21) then the collision integrals are iden-tical since G A ( t − , t ) = i and iρ ( t ) = G < ( t, t ). Weconclude that the G < ( t, t ) that solves the HF and HF-GKBA equations is the same provided that we use thesame G R of Eq. (21). This observation contains usefulhints on how to approximate the quasi-particle propaga-tor of open systems without paying a too high computa-tional price. We emphasize that the locality in time ofthe retarded embedding self-energy and of the HF self-energy Σ HF ( z, z ′ ) = δ ( z, z ′ )[ h HF ( z ) − h ( z )] are distinct T λ T τ T T n μ TT n T n T τ T λ μ ε n ε n ε n ε n FIG. 2. Example of an open system as described in the maintext with N τ = 9 transverse channels and a chain of four sites. and should not be lumped together. The former is purelyimaginary and hence Σ < em = 0 whereas the latter is purelyreal and hence Σ < HF = 0. Alternatively we can say thatΣ HF is local on the contour whereas Σ em is not. This is acrucial difference: in closed systems the off-diagonal HF-GKBA G < ( t, t ′ ) is the same as the HF G < ( t, t ′ ) whereasin open systems it remains an approximation even for aWBL embedding self-energy. Only the diagonal HF andHF-GKBA G < ( t, t ) are identical in this case.
2. An approximate propagator for mean-field electrons
In most physical situations the removal and additionenergies relevant to describe the electron dynamics of themolecular junction after the application of a voltage dif-ference or a laser pulse are well inside the continuumspectrum of the reservoirs. It is therefore natural to studyhow well the GKBA equation performs when G R is cho-sen as in Eq. (21) withΓ = i (cid:2) Σ Rem ( µ ) − Σ Aem ( µ ) (cid:3) . (25)In Eq. (25) the quantity Σ Rem ( µ ) is the Fourier trans-form of the equilibrium embedding self-energy evaluatedat the chemical potential. This choice of Γ is expected toyield accurate results whenever Σ Rem ( ω ) depends weaklyon ω for frequencies around µ . Let us address this issuenumerically. We consider a class of systems consistingof two reservoirs, α = L, R , with N τ transverse channelsand a nanostructure with a chain geometry, see Fig. 2.We use a tight-binding representation and characterizethe Hamiltonian of the reservoirs by a transverse hop-ping T τ and a longitudinal hopping T λ between nearestneighboring sites, and an onsite energy ǫ = µ (half-filledreservoirs). The molecular chain has matrix elements h ij = T n between nearest neighboring sites i and j and h ii = ǫ n on the diagonal. The left reservoir is contactedthrough its middle terminal site to the leftmost site ofthe chain while the right reservoir is contacted throughits middle terminal site to the rightmost site of the chain.We denote by T the corresponding matrix elements of theHamiltonian. In Fig. 3 we compare GKBA versus full KBE resultsfor noninteracting and HF electrons. In all cases the
GKBAKBE1dWBL t T λ = - 5 T λ = - 2 one-site chain one-site chain one-site chain T λ = - 9 t tn g n t t I V L -V R = 1.6V L -V R = 2.4 FIG. 3. Top panels: density n = ρ of a one-site chainconnected to leads with N τ = 1 after the sudden switch-onof a bias V L = 2 for different T λ = − , − , −
2. Bottomleft panel: HF density of site 1 of a two-site chain connectedto leads with N τ = 1 after the sudden switch-on of a bias V L = − V R = 1. Bottom right panel: HF current at the rightinterface of a four-site chain connected to leads with N τ = 9after the sudden switch-on of a bias V L = − V R = 0 . , . Coulomb integrals v ijmn = δ in δ jm v ij . The top panelsrefer to a system with N τ = 1 and a single-site chaindriven out of equilibrium by a bias V L = 2 and V R = 0.The parameters (in arbitrary units) are µ = ǫ n = 0, v ij = 0, and T = p γ | T λ | / γ = 0 .
4. From Eq.(25) we find Γ = 2 γ . The simulations have been per-formed at zero temperature for three different values of T λ = − , − , −
2, and are compared with exact numeri-cal results obtained using the algorithm of Ref. 41. Asexpected the agreement deteriorates with decreasing thebandwidth W = 4 | T λ | of the reservoirs since Σ Rem ( ω ) ac-quires a strong dependence on ω for ω in the bias window.The dashed lines indicate the steady-state value of n forone-dimensional reservoirs and for WBL reservoirs. KBEcorrectly approaches the one-dimensional steady-state inall cases whereas GKBA approaches the WBL steady-state only in the limit | T λ | → ∞ . In the bottom left panelwe consider a two-site chain driven out of equilibriumby a bias V L = − V R = 1 and again connected to one-dimensional reservoirs. In this case, however, the sys-tem is interacting and treated in the HF approximation.The chemical potential is chosen in the middle of theHOMO-LUMO gap of the disconnected chain with twoelectrons. For T n = − ǫ n = 0, and Coulomb integrals v = v = 2, v = v = 1 one finds µ = 2. The restof the parameters are T λ = − . T = − . ≃ .
67 for the GKBA simula-tions. Even though the HOMO-LUMO gap ∆ HL = 2 is -5 0 50.60.81 KBEGKBAWBLA -6 -4 -2 000.20.40.60.8 tn t in = 0t in = -2t in = -4t in = -6 n t FIG. 4. Results for the density n ( t ) of the one-site chainwith T λ = − n ( t in ) is varied. In theright panel n ( t in ) = 1 /
2, a bias V L = 2 is switched on at t = 0 and t in is varied. For clarity the curves with t in = − n , n = 0 , , , n/ not much smaller than the bandwidth W = 4 | T λ | = 6 westill observe a satisfactory agreement for the density ofsite 1 (a similar agreement is found for site 2, not shown).The damping time as well as the amplitude and frequencyof the transient oscillations are well reproduced; further-more the GKBA steady-state value differs by less than1% from the corresponding KBE value. The accuracy ofthe HF-GKBA is not limited to the diagonal matrix el-ements of the density matrix. This is exemplified in thebottom right panel where we show the current flowing atthe right interface of the four-site chain of Fig. 2 with N τ = 9 transverse channels, bias V L = − V R = 0 . , . µ = 2 .
26 (chosen in the middle of theHOMO-LUMO gap of the disconnected chain with 4 elec-trons), T n = − T λ = T τ = − T = − . ǫ n = 0 andCoulomb integrals v ii = v = 1 . v ij = ( v/ / | i − j | for i = j . The GKBA and KBE currents are in excel-lent agreement except for a slight overestimation of theGKBA steady-state value at small bias.In conclusion the GKBA equation with G R from Eq.(21) and Γ from Eq. (25) is a good approximation tostudy the HF dynamics of open systems provided thatthe embedding self-energy of the reservoirs has a weakfrequency dependence around the chemical potential.
3. Relaxation and local thermalization
For the GKBA results of Fig. 3 we started the propa-gation at time t in < ρ ( t ) thermalize in the absenceof external fields until t = 0 when a bias is switched on.For t in sufficiently remote in the past the density matrixattains a steady value ρ eq before the system is biased.By definition ρ eq is the static solution of Eq. (16) with dρ/dt = 0; therefore if we start with ρ ( t in ) = ρ eq then thedensity matrix remains constant in the interval ( t in , n ( t ) = ρ ( t ) = 1 / -5 0 50.50.60.70.80.91 Γ = 0.0Γ = 0.4Γ = 0.8 -5 0 5 10 -10 -5 0 5 10 t ttn FIG. 5. Time dependent occupation of the one-site chain fordifferent Γ. Left panel: t in = − t in = − V L = 2 switched on at t = 0.Right panel: t in = −
12 and bias V L = 2 switched on at t = 0. all t < n ( t in = −
6) = 1 / ρ eq ) bystarting the propagation at t = 0 with ρ (0) = ρ eq . Thisinitial condition guarantees the local thermalization of allone-time observables. However, in a fully relaxed systemany two-time correlator depends on the time-differenceonly, and to achieve this relaxation a “memory buffer”is needed. Suppose that we start the propagation with G < ( t in , t in ) = iρ eq . Then the equal-time G < ( t, t ) remainsconstant but the G < ( t, t ′ ) depends on t and t ′ separately.It is only for t, t ′ large enough that G < ( t, t ′ ) depends on t − t ′ . This concept is explained in the right panel of Fig.4 where we display n ( t ) when a bias V L = 2 is switchedon at t = 0. In all cases ρ ( t in ) = ρ eq = 1 / t in is varied. The absence of relaxation for toosmall | t in | is evident from the strong dependence of thetransient behavior on t in . The curves n ( t >
0) becomeindependent of t in only for t in . − ρ ex . If westart the propagation at time t = 0 with initial condition ρ (0) = ρ ex then the transient behavior is affected by spu-rious relaxation processes. The proper way of performingGKBA simulations consists in driving the relaxed systemtoward ρ ex with some suitable external fields.
4. Damping
For bulk systems like an electron gas the inclusion ofdamping in the propagator worsens the agreement withthe KBE results. In fact, the use of a non-Hermitianquasi-particle Hamiltonian h HF − i Γ / G R is a distinc-tive feature of open systems. Here we address how sensi-tive the results are to different values of Γ. We consider again the noninteracting one-site chain of Fig. 3 with T λ = −
9, for which Eq. (25) yields Γ = 0 .
8. In all caseswe set the initial condition n ( t in ) = 1. In Fig. 5 (leftpanel) we show the relaxation dynamics, starting from t in = −
8, of the unperturbed system for three differentΓ; the curves are essentially on top of each other. Thismay suggest that the dependence on Γ is weak. However,if we switch on a bias in the left lead V L = 2 at time t = 0(middle panel) we appreciate a strong Γ-dependence. Wemay argue that for small Γ the relaxation time is longerand hence that the curves with Γ = 0 . , . . t in . This is notthe case as clearly illustrated in the right panel where t in = −
12. The curve with Γ = 0 . .
8. The ap-parent weak Γ-dependence in the left panel is simply dueto the alignment of the on-site energy to the chemicalpotential, µ = ǫ n = 0. In general a proper choice of thequasi-particle damping is crucial for a correct descrip-tion of the system evolution. In HF theory the dampingis only due to embedding effects and the Γ of Eq. (25) isthe most accurate. The inclusion of correlation effects in-troduces an extra damping. Is it possible to maintain thesimple form in Eq. (21) for the quasi-particle propaga-tor and still have good agreement with the KBE results?In the next Section we discuss two different correlatedquasi-particle Hamiltonians to insert in Eq. (21). III. CORRELATED APPROXIMATIONS TOTHE PROPAGATOR
In the interacting case the exact equation of motionfor G R reads (cid:20) i ddt − h HF ( t ) (cid:21) G R ( t, t ′ ) = δ ( t, t ′ ) + Z d ¯ t Σ R ( t, ¯ t ) G R (¯ t, t ′ )(26)with Σ R = Σ Rem + Σ Rc . If we approximateΣ Rem ( t, t ′ ) ≃ − ( i/ δ ( t − t ′ ) , (27)with Γ from Eq. (25), we find the approximate equation (cid:20) i ddt − h (0) qp ( t ) (cid:21) G R ( t, t ′ ) = δ ( t, t ′ ) + Z d ¯ t Σ Rc ( t, ¯ t ) G R (¯ t, t ′ )(28)where h (0) qp ( t ) ≡ h HF ( t ) − i Γ / G R has to incorporate correlation ef-fects to some extent. Below we propose two schemes toapproximate the convolution Σ R G R and reduce Eq. (28)to a quasi-particle equation of the form (cid:20) i ddt − h qp ( t ) (cid:21) G R ( t, t ′ ) = δ ( t, t ′ ) . (30)The solution of Eq. (30) is G R ( t, t ′ ) = − iθ ( t − t ′ ) T e − i R tt ′ d ¯ t h qp (¯ t ) (31)and satisfies the group property of Eq. (22). Therefore,if we are successful in this task the calculation of G R willscale like t . A. Static correlation approximation
In open systems the correlation self-energy decays tozero when the separation between its time arguments ap-proaches infinity. If G R (¯ t, t ′ ) ≃ G R ( t, t ′ ) for t − ¯ t smallerthan the decay time of Σ Rc we can approximately write Z d ¯ t Σ Rc ( t, ¯ t ) G R (¯ t, t ′ ) ≃ (cid:20)Z d ¯ t Σ Rc ( t, ¯ t ) (cid:21) G R ( t, t ′ ) . (32)To evaluate the integral in the square brackets we makean adiabatic approximation on top of the GKBA, i.e., wereplace G R with the equilibrium propagator of a systemdescribed by the Hamiltonian ˆ H ( t ). Let us consider, forsimplicity, an interaction v ijmn = δ in δ jm v ij . Then theLangreth rules provides us with the following expres-sion of the retarded 2B self-energy, see Eq. (11),Σ Rc ,ij ( t, t ′ ) = 2 X kl v ik v jl [ G R ij ( t, t ′ ) G An alternative way to introduce correlation effects inthe propagator is again based on the adiabatic approx-imation but uses the concept of quasi-particles. Let usrepresent operators in the one-particle Hilbert space witha hat, e.g., ˆ h qp or ˆΣ Rc , and denote by | i i the basis ketof the molecular junction so that h i | ˆΣ Rc | j i = ˆΣ Rc ,ij , etc..For an isolated molecule in equilibrium the quasi-particleequation reads [ˆ h HF + ˆΣ Rc ( ǫ )] | ϕ i = ǫ | ϕ i (37)where ˆΣ Rc ( ǫ ) is the Fourier transform of the equilibriumself-energy. To lowest order in Σ Rc this equation impliesthat the correction to the HF energies ǫ HF ,n is ǫ qp,n = ǫ HF ,n + h ϕ n | ˆΣ Rc ( ǫ HF ,n ) | ϕ n i (38)where | ϕ n i is the eigenket of ˆ h HF with eigenvalue ǫ HF ,n .Equation (38) suggests to construct a quasi-particleHamiltonian in the following manner. We evaluate againthe GKBA form of Eq. (33) with the propagator of Eq.(34) and then calculate˜Σ( t, ω ) = Z dt e iω ( t − t ′ ) ˜Σ( t, t − t ′ ) . (39)From this quantity we construct the one-particle oper-ator ˆ˜Σ( t, ω ) = P ij | i i ˜Σ ij ( t, ω ) h j | and subsequently thediagonal self-energy operator in the HF basisˆ˜Σ( t ) = X n | ϕ n ih ϕ n | ˆ˜Σ( t, ǫ HF ,n ) | ϕ n ih ϕ n | . (40) I It t FIG. 6. Time-dependent current at the right interface of thefour-site junction with V L − V R = 1 . Imposing now that ˆ h qp ( t ) = ˆ h (0) qp ( t ) + ˆ˜Σ( t ) we get a self-consistent equation for ˜Σ( t ). We refer to this proce-dure as the GKBA + Quasi-Particle (QP) scheme. Asthe Fourier transform of ˜Σ( t, t − t ′ ) has to be evalu-ated in N different energies the calculation of ˜Σ( t ) inthe GKBA+QP scheme scales like N . IV. RESULTS In this Section we study the nonequilibrium correlated dynamics of the chain junction of Fig. 2 and of a modelphotovoltaic junction. We calculate local occupations,currents, and spectral functions using different GKBAschemes, and benchmark the results against full KBEsimulations. A clear-cut scenario will emerge in whichGKBA+SC is the most reliable scheme while all otherschemes suffer from some deficiencies. A. Chain junction Nonequilibrium correlation effects change drasticallythe HF picture of quantum transport. The applied biascauses an enhancement of quasi-particle scatterings andconsequently a substantial broadening of the spectralpeaks. The 2B steady current is larger (smaller) thanthe HF steady current at bias smaller (larger) than theHF HOMO-LUMO gap, see bottom-right panel of Fig. 4.In Fig. 6 we compare the current at small (left panel) andlarge (right panel) bias using KBE and different GKBAschemes. Even though the correlation-induced enhance-ment (at small bias V L − V R = 1 . ∼ . V L − V R = 2 . ∼ . 11) of the steady cur-rent relative to the HF values is qualitatively capturedby all GKBA schemes, quantitative differences emerge.GKBA0 is rather close to KBE during the initial tran- L + I R : GKBA+QP dN/dt : GKBA+QPI L + I R : GKBA+SC dN/dt : GKBA+SC t FIG. 7. Numerical evidence of the fulfillment of the continuityequation in the GKBA+QP and GKBA+SC schemes. sient but considerably overestimates the steady state.GKBA+QP corrects too much this deficiency and thesteady current is appreciably underestimated. Further-more the transient behavior worsens: the first peak isabsent and the current saturates too fast. This is due toa general problem of the GKBA+QP scheme. The equi-librium ˜Σ is too large or, equivalently, equilibrium cor-relations are overestimated. GKBA+SC gives an overallimprovement. The transient current reproduces severalKBE features (oscillation frequency and relative hight ofthe peaks) and the steady current is very close to theKBE value.By construction the GKBA schemes guarantee the sat-isfaction of the continuity equation. The rate of changeof the total number of electrons in the nanostructure, dN/dt , is equal to the sum of the currents flowing throughthe left and right interface, I L + I R . In Fig. 7 weshow that this analytic property is numerically confirmedwith high accuracy in the GKBA+QP and GKBA+SCschemes. B. Photovoltaic junction We consider a more complicated open system with thefeatures of a photovoltaic molecular junction. Inspiredby a paper by Li et al. we model the junction as adonor–acceptor complex connected to a left and right lh T A T A T A T DA T L,h T L RT hl e i ω t FIG. 8. Schematic illustration of the photovoltaic junctiondescribed in the main text. electrodes (reservoirs), see Fig. 8. The donor is describedby HOMO ( h ) and LUMO ( l ) levels and the LUMO isconnected to a chain of four acceptor sites each describedby a single localized orbital. These orbitals are mixed bythe acceptor Hamiltonian and form two valence and twoconduction levels. The junction is connected to the leftelectrode throught the HOMO with tunneling amplitude T L,h and to the right electrode through the rightmost ac-ceptor site with tunneling amplitude T ,R . The explicitform of the Hamiltonian of the donor–acceptor complexis ˆ H J = ǫ h ˆ n h + ǫ l ˆ n l + ǫ A X a =1 ˆ n a + T DA X σ (cid:16) ˆ d † lσ ˆ d σ + ˆ d † σ ˆ d lσ (cid:17) + T A X a =1 X σ (cid:16) ˆ d † aσ ˆ d a +1 σ + ˆ d † a +1 σ ˆ d aσ (cid:17) + U DA (ˆ n h + ˆ n l − X a =1 ˆ n a − a , (41)where ˆ n x = P σ ˆ d † xσ ˆ d xσ is the occupation operator for x = h, l, a . The interaction between the excess chargesof the donor and acceptor chain implicitly fixes the con-dition of charge neutrality. For one-dimensional reser-voirs with longitudinal hopping integral T λ = − 9, tun-neling amplitudes T L,h = T ,R = − . 3, donor levels ǫ h = − . ǫ l = − . 92, acceptor levels ǫ A = − . T DA = − . 1, intra-acceptor hop-ping T A = − . U DA = 0 . 5, the chemicalpotential µ = 0 . 04 is in the middle of the HF gap betweenthe valence and conduction acceptor levels. The equilib-rium system has HOMO and LUMO occupations 2 and0 respectively and the two valence levels of the accep-tor chain completely filled. The photovoltaic junction isdriven out of equilibrium by irradiation with monochro-matic light. For simplicity we assume that the light cou-ples only to the donor dipole moment and henceˆ H light ( t ) = s ( t ) T hl X σ (cid:16) e iωt ˆ d † hσ ˆ d lσ + e − iωt ˆ d † lσ ˆ d hσ (cid:17) , (42)where s ( t ) is a switching function. We consider T hl =0 . ω = 2 = | ǫ h − ǫ l | and study a pulse, s ( t ) = 1 for0 < t < π/T hl and zero otherwise, as well as continuousradiation, s ( t ) = 1 for t > s ( t ) = 0 for t < H J in the form of Eq. (3). The one-particleHamiltonian h ij with i, j = h, l, a then reads h = ˜ ǫ h ǫ l T DA T DA ˜ ǫ A T A T A ˜ ǫ A T A 00 0 0 T A ˜ ǫ A T A T A ˜ ǫ A (43) 012 n l n h KBE0 25 50 75 100 12500.511.5 n n n n cc up a ti on s t (a) tn l FIG. 9. Time-dependent occupations in the HF approxima-tion using GKBA (solid) and KBE (circles). The junction isperturbed by a monochromatic pulse. with ˜ ǫ h = ǫ h + U DA , ˜ ǫ l = ǫ l + U DA and ˜ ǫ A a = ǫ A + U DA − U DA /a . For the Coulomb integrals wefind v ijmn = δ in δ jm v ij with v = U DA / / / 40 0 1 1 / / / 41 1 0 0 0 01 / / / / / / . (44)Let us start with the mean-field analysis of the lightpulse. The duration π/T hl has been chosen to get a pop-ulation inversion of the HOMO and LUMO levels. InFig. 9 we show the HF occupations of the donor (toppanel) and acceptor (bottom panel) levels in GKBA andKBE. The impressive agreement is due to the fact thatfor T λ = − U DA . Dur-ing the pulse the HOMO level is partially refilled by theleft reservoir and the total charge on the donor over-comes 2. This excess charge is instantaneously felt byA1 which starts expelling electrons at a rate larger thanthe tunneling rate from LUMO to A1. We also observethat the charge transfer between LUMO and A1 is noteffective. The inset shows the LUMO occupation on alonger time scale. Electrons remain trapped and slosharound along the junction. In fact, in HF no steady-state is reached. The occurrence of self-sustained chargeoscillations in mean-field treatments has been observedin similar contexts and is most likely an artifact ofthe approximation. As we shall see, the correlated KBEresults are very different. Therefore the collision integraland the correlated propagator of the GKBA approach0 -60 -40 -20 00.010.0150.02 GKBA+SCHFGKBA+QPGKBA0 n h n l n n n n tt FIG. 10. Thermalization of the occupations in the correlatedcase. The correlated KBE value is represented by a dottedhorizontal line. have to correct the HF theory in a substantial manner.With the inclusion of correlations a deficiency of theGKBA+QP scheme emerges already during the thermal-ization process. In Fig. 10 the donor and acceptor oc-cupations are propagated within different schemes in theabsence of external fields using the HF value of the un-contacted system and ˜Σ( t in ) = 0 as initial conditions.Both GKBA0 and GKBA+SC thermalize, similarly toHF, to values very close to the equilibrium values of thecorrelated (2B) KBE approach (dotted horizontal line).In fact, in KBE the HF and 2B equilibrium occupationsare essentially the same since the correlation self-energy,except for a slight renormalization of the quasi-particleenergies (image charge effect), does not affect the widthof the spectral peaks. For the GKBA to reproduce theKBE thermalized values the imaginary part of ˜Σ has tobe small, and this is not the case in GKBA+QP. HereIm[ ˜Σ ll ( t )] and Im[ ˜Σ hh ( t )] tend to increase thus broad-ening the HOMO and LUMO spectral peaks. Hence theHOMO looses charge whereas the LUMO acquires chargeand the donor polarizability increases. This makes thefirst bubble diagram of the 2B self-energy larger, andtherefore the HOMO and LUMO spectral peaks morebroadened. In a separate calculation (not shown) we sim-ulated the GKBA+QP thermalization and found that thethermalization process is extremely slow, t in . − ∼ . 7, well above the KBE result.We are now ready to show the correlated results in thecase of a light pulse. The KBE occupations are shown inFig. 11 and are considerably different from the HF oc-cupations of Fig. 9. The GKBA+SC scheme is in fairlygood agreement with KBE for all occupations. To illus-trate the crucial role played by our correlated propagatorwe also display the LUMO and A1 occupations in theGKBA0 scheme (dotted-dashed line). Even though theinitial transient is acceptable the GKBA0 occupations GKBA+SCGKBA0KBE acce p t o r o cc up a ti on s n n n n tn l n h FIG. 11. Time-dependent occupations after a pulse in KBEand GKBA+SC. For n l and n we also show the results ofthe GKBA0 scheme (dotted-dashed). become soon inaccurate. Therefore the evaluation of theGKBA collision integral with HF propagator performsrather poorly in open systems. The GKBA+SC schemehas the merit of working both in and out of equilibrium.Like the only goal of TDDFT is to reproduce the den-sity of an interacting system, so the only goal of theGKBA is to reproduce the density matrix of an inter-acting system. The TDDFT or GKBA spectral function A ( t, t ′ ) = i [ G R ( t, t ′ ) − G A ( t, t ′ )] can be very different fromthe true one. This is, however, not always the case. InFig. 12 we show the time evolution of the KBE andGKBA+SC total spectral function defined according to A ( T, ω ) = − Z dτ e iωτ Tr h G R ( T + τ , T − τ i , (45)where T = ( t + t ′ ) / τ = t − t ′ is the relative time. Remarkably the twospectral functions have several common features. Themost important one is the broadening of the spectralpeaks after the pulse and the long elapsing time to relaxback to the equilibrium state. Another common featureis the drift of the acceptor peaks toward higher energyand the merging of the two middle peaks of the accep-tor chain. In GKBA0 the spectral peaks are sharp at alltimes whereas in GKBA+QP they are broadened at alltimes (not shown).To end our discussion on the performance of GKBAin open systems we consider in Fig. 13 the occupationsfor the continuous radiation. Here GKBA+SC is not asaccurate as in the case of the light pulse. However, theagreement with KBE remains satisfactory. The HOMOand LUMO occupations are essentially indistinguishablefrom the KBE values (top panel). The occupations of theacceptor sites next to the right electrode (A3 and A4)1 FIG. 12. Time-dependent KBE (top panel) and GKBA+SC(bottom panel) spectral functions for the photovoltaic junc-tion subject to a light pulse. The light pulse is switched onat time t = 40 for the KBE and t = 75 for the GKBA+SC.The parameters are the same as in Fig. 11. are slightly underestimated in GKBA+SC but the over-all trend (transient oscillations and steady-state value)are correctly reproduced (middle panel). A more quan-titative agreement is observed for the acceptor sites nextto the donor (A1 and A2). For the A1 occupation wealso show the GKBA0 occupation (bottom panel) andwe note again that after a short time the result deviatesconsiderably from the KBE result.Is there any possibility of improving over theGKBA+SC scheme using a different ˜Σ, or the only wayis to go beyond the GKBA? In the bottom panel of Fig.13 we display the A1 occupations for a hybrid schemein which ˜Σ is calculated from GKBA+SC at negativetimes (thermalization) and from GKBA+QP at positivetimes. The improvement up to times t ∼ 200 is impres-sive and extend to all acceptor occupations (not shown).Instead, for times t > 200 the KBE results are closerto those of the GKBA+SC scheme. More generally for t . 200 we observed that the hybrid scheme performsbetter than GKBA+SC for ω ∼ ǫ h − ǫ l (large current HOMO (GKBA+SC)LUMO (GKBA+SC)HOMO (KBE)LUMO (KBE) n n n GKBA+SCKBEGKBA0GKBA Hybrid t n n , n , n n h , n l FIG. 13. Time-dependent occupations in the presence of con-tinuous radiation in KBE and GKBA+SC. For n we alsoshow the results of the GKBA0 scheme as well as of the hy-brid scheme (thermalization with GKBA+SC, positive-timepropagation with GKBA+QP). in the junction) and worse otherwise (small current inthe junction). The purpose of this investigation is toprovide numerical evidence of the existence of a ˜Σ foraccurate GKBA simulations, and hence the possibility ofimproving the GKBA+SC scheme without increasing thecomputational cost. V. SUMMARY AND OUTLOOK We demonstrated that time-dependent NEGF simula-tions of molecular junctions (and more generally openquantum systems) can be considerably speeded up. Dif-ferent GKBA-based schemes have been proposed andsubsequently benchmarked against full KBE calculations.The GKBA+SC scheme turned out to be the most accu-rate both in and out of equilibrium, while still offering asignificant computational gain (for the longest propaga-tion ( t max = 300) of the photovoltaic junction the CPUtime is ∼ 10 minutes in GKBA+SC and ∼ 20 hours inKBE). We also showed that the GKBA+SC scheme can,in principle, be further improved without rising the com-putational price.All calculations have been performed within the 2Bapproximation for the correlation self-energy but theGKBA+SC scheme is completely general and not lim-ited to this special case. Clearly, in large nanostruc-tures screening is important and the interaction shouldbe treated, at least, within the GW approximation. An-other urgent extension of the GKBA+SC scheme is theinclusion of the interaction between electrons and nuclei.2This can be done either at the level of the Ehrenfestapproximation or by adding diagrams with electron-phonon vertices to the correlation self-energy. The GKBA+SC scheme, its extensions and refine-ments can be implemented in ab initio molecular codes to perform first principle time-dependent simulations ofopen nanostructures. Foreseeable applications are, e.g.,in the field of molecular photovoltaics and molecular elec-tronics. Here there is much interest in developing effi-cient quantum simulation methods for an accurate de-scription of the electron-hole formation, recombinationand separation as well as of charge transfer and possiblyionic reorganization or isomerization. In molecular pho-tovoltaics ab-initio studies have focused on the opticalspectra using linear response TDDFT or the Bethe-Salpeter equation. Real-time simulations remain, how-ever, the most powerful tool to resolve the different com-peting processes up to the ps time scale. State-of-the-art simulations treat the contacts as finite-size clusterswhile taking into account the full atomistic structure ei-ther semi-empirically or fully ab initio. However, these studies suffer from spurious boundary effects likethe formation of artificial electric fields and reflection ofcharge after a few tens of fs. There are no such limita-tions in the GKBA for open systems as the electrodes aredescribed in a virtually exact way through the embeddingself-energy. Furthermore the effects of the Coulomb in-teraction can be systematically included through the dia-grammatic expansion of the correlation self-energy. Theencouraging results presented in this work should fosteradvances in the development of a NEGF approach to ul-trafast processes at the nanoscale. Acknowledgements E. Perfetto and G. Stefanucci acknowledge fundingby MIUR FIRB Grant No. RBFR12SW0J. R. vanLeeuwen and A.-M. Uimonen acknowledge support fromthe Academy of Finland as well as the CSC IT center forproviding resources for scientific computing. X.-Y. Zhu, J. Phys. Chem. B , 8778 (2004). A. Nitzan, Chemical Dynamics in Condensed Phases (Ox-ford University Press, 2006). R. M. Dreizler, E. K. U. Gross, Density Functional The-ory: An Approach to the Quantum Many-Body Problem (Springer-Verlag, Berlin, 1990). R. G. Parr and W. Yang, Density-Functional Theory ofAtoms and Molecules (Oxford University Press, New York,1989). Time-Dependent Density Functional Theory , edited by M.A. L. Marques, C. A. Ullrich, F. Nogueira, A. Rubio, K.Burke and E. K. U. Gross (Springer, Berlin, 2006). C. Ullrich, Time-Dependent Density-Functional Theory:Concepts and Applications (Oxford University Press, Ox-ford, 2012). See e.g. Single Charge Tunneling: Coulomb Blockade Phe-nomena in Nanostructures , eds. H. Grabert and M. H. De-voret, (1992) NATO ASI Series: B, Physics, Vol. 234, NewYork, Plenum Press. H. Haug and A.-P. Jauho, Quantum Kinetics in Transportand Optics of Semiconductors (Springer, 2007). J. C. Cuevas and E. Scheer, Molecular Electronics: AnIntroduction to Theory and Experiments (World Scientific,2010). J. B. Neaton, M. S. Hybertsen and S. G. Louie, Phys. Rev.Lett. , 216405 (2006). K. Kaasbjerg and K. Flensberg, Nano Letters , 3809(2008). K. Thygesen and A. Rubio, Phys. Rev. Lett. , 046802(2009). J. M. Garcia-Lastra, C. Rostgaard, A. Rubio and K. S.Thygesen, Phys. Rev. B , 245427 (2009). P. My¨oh¨anen, R. Tuovinen, T. Korhonen, G. Stefanucciand R. van Leeuwen, Phys. Rev. B , 075105 (2012). G. Stefanucci and R. van Leeuwen, Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction (Cambridge University Press, 2013). L. P. Kadanoff and G. Baym, Quantum Statistical Mechan-ics (Westview Press, 1994). N.-H. Kwong and M. Bonitz Phys. Rev. Lett. , 1768(2000). N. E. Dahlen and R. van Leeuwen, Phys. Rev. Lett. ,153004 (2007). P. Lipavsk´y, V. ˇSpiˇcka and B. Velick´y, Phys. Rev. B ,6933 (1986). H. S. K¨ohler, Phys. Rev. E , 3145 (1996). M. Bonitz, D. Kremp, D. C. Scott, R. Binder, W. D.Kraeft and H. S. K¨ohler, J. Phys.: Condens. Matter ,6057 (1996). R. Binder, H. S. K¨ohler, M. Bonitz and N. H. Kwong,Phys. Rev. B N. H. Kwong, M. Bonitz, R. Binder, and H. S. K¨ohler,Phys. Stat. Sol. B , 197 (1998). M. Bonitz, D. Semkat and H. Haug, Eur. Phys. J. B , 309(1999). H. Haug, Phys. Stat. Sol. B , 139 (1992). A. Marini, J. Phys: Conf. Proc. , 012003 (2013). G. Pal, Y. Pavlyukh, W. H¨ubner, H. C. Schneider, Eur.Phys. J. B , 327 (2011). G. Pal, Y. Pavlyukh, W. H¨ubner, H. C. Schneider, Eur.Phys. J. B , 483 (2009). D. Hochstuhl and M. Bonitz, J. Phys: Conf. Proc. ,012007 (2013). K Balzer, S Hermanns and M Bonitz, J. Phys: Conf. Proc. , 012006 (2013). P. My¨oh¨anen, A. Stan, G. Stefanucci and R. van Leeuwen,EPL , 67001 (2008). Initial correlations can survive in low-dimensional and in-teracting reservoirs, see E. Perfetto, G. Stefanucci, and M.Cini, Phys. Rev. Lett. , 156802 (2010). P. My¨oh¨anen, A. Stan, G. Stefanucci and R. van Leeuwen,Phys. Rev. B , 115107 (2009). A.-M. Uimonen, E. Khosravi, A. Stan, G. Stefanucci, S.Kurth, R. van Leeuwen, and E. K. U. Gross, Phys. Rev. B , 115103 (2011). The second Born self-energy in Eq. (11) with all four-indextwo-electron integrals has been calculated in Ref. 18. D. C. Langreth, in Linear and Nonlinear Electron Trans-port in Solids , edited by J. T. Devreese and E. van Doren(Plenum, New York, 1976), pp. 3–32. Alternatively one could add a vertical track (imaginaryMatsubara axis) to the contour and starts the real-timepropagation at t = 0, see O. V. Konstantinov and V. I.Perel’s, Sov. Phys. JETP , 142 (1961); P. Danielewicz,Ann. Physics , 239 (1984). The GKBA is an ansatz for the lesser/greater Green’s func-tion. At present there exists no extension of the GKBA tothe left/right Green’s functions for the inclusion of initialcorrelations with the addition of a vertical track (imagi-nary Matsubara axis). Therefore, relaxation can only beachieved by propagating the unperturbed systems from t in < t = 0. M. Bonitz and D. Kremp, Phys. Lett. A , 83 (1996). For the calculation of Σ em in these geometries see, e.g.,Ref. 33 and S. Sanvito, C. J. Lambert, J. H. Jefferson, A.M. Bratkovsky, Phys. Rev. B , 11936 (1999). S. Kurth. G. Stefanucci, C. O. Almbladh, A. Rubio and E.K. U. Gross, Phys. Rev. B , 035308 (2005). K. S. Thygesen, Phys. Rev. Lett. , 166804 (2009). G. Li, A. Nitzan and M. A. Ratner, Phys. Chem. Chem.Phys. , 14270 (2012). E. Khosravi, A.-M. Uimonen, A. Stan, G. Stefanucci, S.Kurth, R. van Leeuwen and E. K. U. Gross Phys. Rev. B , 075103 (2012). C. Verdozzi, G. Stefanucci and C.-O. Almbladh, Phys. Rev.Lett. , 023001 (2007). C. A. Rozzi et al., Nature Comm. , 1602 (2013). For progresses in this direction see, e.g., A. Marini, C.Hogan, M. Gr¨uning and D. Varsano, Comput. Phys. Com-mun. , 1392 (2009). F. de Angelis, S. Fantacci and A. Selloni, Nanotech. ,424002 (2008). X. Blase, C. Attaccalite and V. Olevano, Phys. Rev. B ,115103 (2011). See, e.g., L. G. C. Rego and V. S. Batista, J. Am. Chem.Soc. , 7989 (2003). W. R. Duncan , W. M. Stier and O. V. Prezhdo, J. Am.Chem. Soc. , 7941 (2005). W. R. Duncan and O. V. Prezhdo, Annu. Rev. Phys.Chem.58