# Full text of "Comparative study of many-body perturbation theory and time-dependent density functional theory in the out-of-equilibrium Anderson model"

## See other formats

o (N < Comparative study of many-body perturbation theory and time-dependent density functional theory in the out-of-equihbrium Anderson model A.-M. Uimonen,i^2 £_ Khosravi,3.2 A. Staii,!'^ G. Stefanucci,-*' 5^ 2 S. Kurth,^''^'^ R. van Leeuwen/'^ and E. K. U. Gross^-^ ^ Department of Physics, Nanoscience Center, FI 4OOI4, University of Jyvaskyla, Jyvaskyld, Finland ^European Theoretical Spectroscopy Facility (ETSF) "' Max- Planck Institut fiir Mikrostrukturphysik, Weinberg 2, D-06120 Halle, Germany ^^ Diparttmento di Fisica, Universitd di Roma Tor Vergata, Via della Ricerca Scientifica 1, 00133 Rome, Raly ^INFN, Laboratori Nazionali di Frascati, Via E. Fermi 4O, 00044 Frascati, Raly ^ Nano-Bio Spectroscopy Group, Departamento de Fisica de Materiales, Universidad del Pais Vasco UPV/EHU, Gentro Fisica de Materiales CSIG-UPV/EHU, Avenida de Tolosa 72, E-20018 San Sebastian, Spain ^IKERBASQUE, Basque Foundation for Science, E-48011 Bilbao, Spain (Dated: August 12, 2011) We study time-dependent electron transport through an Anderson model. The electronic inter- actions on the impurity site are included via the self-energy approximations at Hartree-Fock (HF), second Born (2B), GW, and T-Matrix level as well as within a time-dependent density functional (TDDFT) scheme based on the adiabatic Bethe-Ansatz local density approximation (ABALDA) for the exchange correlation potential. The Anderson model is driven out of equilibrium by apply- ing a bias to the leads and its nonequilibrium dynamics is determined by real-time propagation. The time-dependent currents and densities are compared to benchmark results obtained with the time-dependent density matrix renormalization group (tDMRG) method. Many-body perturbation theory beyond HF gives results in close agreement with tDMRG especially within the 2B approxi- mation. We find that the TDDFT approach with the ABALDA approximation produces accurate results for the densities on the impurity site but overestimates the currents. This problem is found to have its origin in an overestimation of the lead densities which indicates that the exchange cor- relation potential must attain nonzero values in the leads. PACS numbers: 72. 10.Bg,71.10.-w,31.15.xm,31. 15.ee a -a o > O X I. INTRODUCTION The process of electron transport through molecules and nanostructures is part of a r apid ly growing research area in condensed matter physic^^^. On a fundamental level one has to deal with time-dependent processes in an open system where different scattering mechanisms such as electron-electron or electron-phonon interactions are of great importance. These factors make the transport problem not only difficult, but also very rich in physical phenomena. Most of the recent studies in molecular elec- tronics have focused on the description of steady-state transport while neglecting short-time dynamics such as transients and fast switching processes. However, these processes will become increasingly important since fast switching rates play a pivotal role in the operation of future devices. For the description of electron transport several nu- merical approaches have been developed that can deal with fully time-dependent systems^"^^. Among these are the time-dependent density matrix renormaliza- tion group (tDMRG) approacB ^, ti me-dependent den- sity functional theory (TDDFTf^^', and self-consistent many-body perturbation theory (MBP T) based on the Kadanoff-Baym (KB) equationpEMII. Each of these methods has its own advantages and disadvantages. To the best of our knowledge, a comparative study of these three approaches on an identical time-dependent sys- tem has not been carried out. Such a study would be very valuable for gaining insight into these methods and into the direction in which each method needs to be im- proved. In the steady-state regime of quantum trans- port such comparisons of many-body and benchmark ap- proaches were made by Wang et al)^ within the GW- approximation and by Schmitt and AnderalSl at second Born (2B) and GW level. In both cases good agreement with benchmark results was found in certain parameter ranges. We want to extend these comparisons to the transient regime as well. Let us give a brief description of the approaches that we use in this work. The tDMRG method is a numer- ical algorithm based on truncation of the Hilbert space of low-dimensional systema^SESl, Jn this work we did not carry out such calculations ourselves but we use pub- lished tDMRG results^ as a benchmark for both the TDDFT and MBPT appro aches In the TDDFT approactP^l^II a system of interacting electrons is mapped, in an exact manner, onto a sys- tem of noninteracting electrons moving in an effective time-dependent external potential known as the Kohn- Sham (KS) potential. The KS potential is functionally dependent on the electron density such that it produces a KS wave function with a density identical to the time- dependent density of the interacting system. It is im- portant to note that TDDFT yields in principle the ex- act tim e-dependent current through a molecular junc- tiorPI2l!. The use of one-particle equations in TDDFT allows for large scale first-principle calculations on re- alistic systems. In practice, however, approximations are unavoidable and the accuracy of a TDDFT calcula- tion crucially depends on the quality of the approximate exchange-correlation (XC) potential used. Most applic a- tions of TDDFT to quantum transport processeJ^nSElEII use the adiabatic approximation which assumes that the XC-potential instantaneously follows the density profile. This is a reasonable assumption when the density changes are slow on a time-scale of typical lead-to-molecule tun- neling rates, and also when the switch-on times of the ap- plied biases are small enough. However, it has also been pointed out that n on-adiabatic effects can have substan- tial influence22E3 on calculated properties. In such cases there is also a need to introduce spatial non-locality in the density functional because the non-localities in space and time are strongly related by conservation laws^. This re- lation is virtually unexplored within a quantum transport context. Gaining further insight into this issue is one of the goals of this work. The MBPT approach based on the KB equationtPI^ has been successfully applied to tinie-dep endent quan- tum transport for model systema^^i^IlMZl^ Xhe method offers the possibility of including relevant physical pro- cesses by means of selection of Feynman diagrams for the self-energy. The electron-electron correlations are thus considered via the many-body self-energy term which is treated perturbatively to infinite order by summa- tion of infinite classes of diag rams. Furthermore by us- ing conserving approximation^^Slso] g^ch as the Hartree- Fock (HF), second Born (2B), GW, and T-Matrix ap- proximations we can guarantee that conservation laws are obeyed, which has shown to be very important in quantum transporlp2Bll. In this approach one has di- rect access to quantities like quasiparticle spectra, life- times, and screened interactions which provide insight into the effects of electron correlation. In particular, the non-locality in time of the 2B, GW, and T-Matrix ap- proximations allows for a description of memory effects and quasi-particle broadening. We use the partition-free scheme where the device is initially contacted to the leads and the whole system in thermal equilibrium^^. In this approach both the transient and steady-state currents have a direct physical meaning as these currents are in- duced by the physical switch-on of a bias. In the par- titioned approaches they are instead induced by switch- on of a device-lead coupling which does not correspond to the standard experimental situation. We finally like to point out that the MBPT approach can be used to derive new improved time-dependent density functionals with memory and conserving properties'*'^ . This has been done successfully within the linear response regime^^Ml, Since both TDDFT and MBPT require the use of ap- proximations it is important to have independent bench- mark results. For the Anderson impurity model such re- sults in the time domain have recently been obtained with tDMRG^^. Therefore we will use this system as a test case for our comparative study of MBPT and TDDFT. The paper is organized as follows: in Sec. In] we intro- duce the model used in our investigation. In Sections |IIA| andUHlwe describe the MBPT and the TDDFT meth- ods used. In Sec. |III| we present numerical results and the last section summarizes our conclusions. II. THE MODEL We study an Anderson impurity modeF^ described by the Hamiltonian Hit)=Hc+J2Hc.{t) + HT, (1) where Hq , Ha , and H^ respectively describe the impurity region, the leads a {= L,R), and the tunneling between the impurity region and the leads. The Hamiltonian for the impurity site reads where cJ.,C(j are fermionic creation and annihilation op- erators and CT, a' are the spin indices, Sq is the on-site energy of the interacting site and U is the interaction term or the charging energy. The Hamiltonian Ha(t), describing the leads is oo -EE(^"4aC.+ l.a+i:^.C.), (3) where Sa is the on-site energy in the leads, Wa is the bias on the lead a and Vq, is the hopping between neighboring lead sites. The tunneling Hamiltonian describes the cou- pling between the impurity site and the leads, and has the form ^T = - ^ ( Vlink Cqo-CIo-L + Mink cl^CiaR + H.C. j , (4) where Mink is the hopping from the leads to the impurity site and vice versa. A. Kadanoff-Baym equations The nonequilibrium properties of the system are stud- ied with the aid of nonequilibrium Green function theory and TDDFT described later in section [Tl B[ The nonequi- librium Green function is defined as the expectation value with respect to the initial state of the contour-ordered product of creation and annihilation operatora^ G tcT^ja' {z,z') = ~^{T[cHM^)c^HM'')])^ (5) where i,j are the site-indices, T denotes the time- ordering operator along the Keldysh contour^, and where the contour variables z and z' specify the posi- tion on the contour^. The subscript H refers to opera- tors in the Heisenberg picture with respect to the time- dependent Hamiltonian iJ(z}22142l. The Green function of the whole system satisfies the equation of motion [id,l ~ U{z)]G{z,z')^6{z,z')l dzS™[G](z,2)G(z,z'), (6) where we introduced the many-body self-energy I]^™[G] which accounts for all the exchange and correlation ef- fects^ and where we suppressed spatial indices. The self- energy is a functional of the Gree n fun ction which in practice is defined diagrammatically^MZl, In this work we solve the equation of motion of the Keldysh Green function fully self-consistentlj0SlEl] using the four approximations of the many-body self-energy I]''™[G] shown in Fig. fll The self-consistent HF approx- imation is time-local and includes the Hartree and the exchange potential. The self-consistent 2B approxima- tion consists of the two diagrams to second order in the interactioEp2l. It describes dynamical screening of the electron-electron interaction via a simple bubble diagram and includes a vertex contribution via the second order exchange diagram. The fully self-consistent GW approxi- matiorpS incorporates the dynamical screening effects via the infinite summation of bubble diagrams^'^l In this ap- proximation, the Coulomb interaction is replaced by the screened potential W . The last approximation w e use is the fully self-consistent T-matrix approximatiorpSES. It contains the 2B diagrams and an infinite summation of the ladder diagrams. The GW and T-matrix approx- imations are complementary since the GW approxima- tion accounts for dynamical screening in infinite systems with long-range Coulombic interactions whereas the T- matrix approximation is known to be important in de- scribing infinite systems with a short range hard-core interactioiPim, When we describe a system attached to noninteracting leads the equation of motion of the Green function for the whole system can be folded into an effective equati on of motion of the Green function for the central regiorP^^. In the case of the impurity model that we consider this gives [id,^H{z)]G{z,z')^5{z,z') (7) + //^- {[Se„,(z,z)-FS"-[G](z,z)]G(z,/)}, where the embedding self-energy Yjcra{z,z') accounts for the tunneling of electrons between leads and the impu- rity site. The many-body self-energy depends only on the GW: T-Matrix: FIG. 1: Diagrammatic representation of the conserving many- body approximations to the self-energy. Wiggly lines denote the many-body interaction. All Green function lines (directed solid lines) are are fully dressed. Green function of the central site as the many-body in- teraction is restricted to the central site only. This Green function has only one spatial index. For the time-dependent observables calculated on the real axis we denote the contour parameter z by the real time t. The time-dependent density for the impurity site is given by no(t) = -^G<(i,i+), (8) where t'^ approaches t from an infinitesimally later time t+ — t + 5. The current through the lead a — (L,R) can be expr essed in terms of the so-called Keldysh Green functions as^I^^ I„(t) = 2Re|^ di[G<(t,i)E*„,,„(i",t) dtG-(i,t)E<^_Jt,i)] dfGT(i,f)sL,Jr,i)l, (9) where we integrated on the Keldysh contour and where the superscripts A, R and < refer to the advanced, re- tarded, and lesser components of the Green function and the self-energy. Further, ] and [ are the mixed compo- nents having one time argu men t on the imaginary axis and another on the real axis^SEfil. The initial many-body correlations and embedding effects are taken into account by the last term in equation ^ which is an integral over the vertical track of the Keldysh contouil^. If we assume that in the i ^- oo limit the terms with components on the imaginary track vanish and that the Green function and the self-energy depend only on i — i' then we can Fourier transform (|9) and we obtain the Meir-Wingreen formula for the steady-state currentP^l I^^-i did , r„(a;) G<{u})-2iTTfaiu})A{uj) (10) where T^iu!) is the imaginary part of the embedding self-energy, fa{uj) is the Fermi function and A{uj) is the steady-state spectral function^. Hence, equation ^ is a generalization of the Meir-Wingreen formula^. We further define the nonequilibrium spectral function A{T,u;) -Im ^e-[G>-G<](T+^,r-^), (11) where t = f — t' is a relati ve time and T = {t + t')/2 is an average time-coordinate^E^I^. In equilibrium, this function is independent of T and has peaks below the Fermi level at the electron removal energies of the sys- tem, while above the Fermi level it has peaks at the elec- tron addition energies. If the time-dependent external field becomes constant after some switching time, then also the spectral function becomes independent of T af- ter some transient period and has peaks at the addition and removal energies of the biased systenPH. B. Time-dependent density functional theory Within TDDFT the complication brought forward by considering an open system can be resolved in a very similar manner as in MBPT (see Section II A I , with the aid of an embedding self-energy. The equation of motion for the k-th single-particle orbital is projected onto the Anderson impurity site and reads [idt - i?"^ (i)] yJk (t) - / dt E^g (i, t}^k (t) Jo + E^-k.9L(t,0)^fc,«(0), (12) where S^g(i, t) is the KS embedding self-energy and g^^ is the retarded lead Green function. This expression is, in principle, exact. If we now assume that the exchange- correlation potential is zero in the leads then E^g(i,i) can be replaced by S^^(t,i) of Eq.([6|. We will assume this in the following. Then for the Anderson impurity model the KS Hamiltonian H^^ (t) has the following from H^'^t) = VKsit) = eo{t) + -Unoit)+v^c[n]{t). (13) The approximation for the XC-potential in this work is based on the local density approximation (LDA) for the static, non-uniform one-dimensional Hubbard model derived from the Bethe ansatz (Bethe ansatz LDA, BALDA) which has been suggested in Ref. [59, and fur- ther been developed in Ref. [60] . The adiabatic versiorpH of this functional ( ABALDA) makes w^c ["] local in both space and time. The modified version of ABALDA for the transport setujp2l jg taking into account the differ- ent hopping between the impurity site and the leads and reads explicitly where "xc(") = --Un-2Viink TTn\ cos I — 1 — cos nn T Here, ^ is a parameter determined by the equation — sin(7r/0 TT da; Jo{x)Ji{x) x[l + exp{Ux/{2Viink))V (15) (16) and Ji=o,i(^) ^^^ Bessel functions. A particularly in- teresting property of the BALDA is its discontinuity at half-miing63 z;,,(l+)-w,,(l-) = C/-4yii„kC0s(|). For the parameters used in this work (see Fig. [2|, the dis- continuity is both positive and negative. However, even if the physical gap should be positive, the results ap- pear not to be affected by this change in signSS. New parametrizations that alleviate this issue are currently being developecP^. The adiabatic approximation implies that 5n{t') S{t-t')U,{n{t)), (17) where /xc = dvxc{n)/dn, meaning the XC-response kernel is local in time (and in space). This local and instanta- neous approximation becomes valid for Hubbard systems in the limit of slowly varying density both in space and in time. These conditions are not satisfied for the quantum transport system under consideration. Despite this fact, reasonable densities were obtained using the BALDA for finite Hubbard chains^^ and it is therefore worthwhile to try the approximation for quantum transport phenom- ena. In the present case this approximation for v^c is only used on the impurity interacting site, since no inter- actions in the leads are present. For future reference we make a connection with the many-body approach of the previous section. The fully 0.4 r 0.2- ._.. u= . u= =1.4 V|^^^ =0 6, V =0.2 =0.5 =0.3535 .... u= u= "■5. V„„ FIG. 2: The BALDA XC-potential as a function of the density for parameters used in the subsequent sections. vl^^^^ir (1 - n)v<Jn) - 0in - 1)«< (2 - n), (14) self-consistent Green function for the whole system (i.e. leads plus impurity) satisfies the following equation Gij{z, z) — Gjj (z, z ) + J2 fdzdz'Gfk'{z,z)[EkiM^,^') kl •' - 5{z,z')5kiVk.^c{z)]Gij{z\z'). (18) Since the exact density is given by both the KS and the exact Green function, i.e., ^^(z) = —iGkk{z,z^) = —iG^f.{z,z'^), it follows that y^ / AzG'^^{z,z)vk,-^c{z)Gki{z,z) ^ k •'c Y, I d^dz' Gf§{z, z)^ki.Uz. z')Gu{z', z), (19) kl ■'c where Sxc is the many-body self-energy with the Hartree potential subtracted. If the self-energy is exact then the corresponding XC-potential that solves this Sham- Schliiter equatioii^ yields the exact density of the sys- tem. We see that the integral kernel on the left hand side of this equation is nonlocal in space and time. Hence, the solution of this integral equation for tife^xc will in general have values on any site k. This has been confirmed by recent work of Schenk et al^^. It is important to note that this is true even if the many-body interactions are restricted to the impurity site only. We therefore make an approximation if we set the XC-potential to zero in the leads. We will discuss the validity of this approximation in the results section. A. Equilibrium results We start by considering a system with interaction U ~ 1 and coupling to leads Mink = 0.5. The Fermi energy of the system is ej? = (half- filling). In Fig. IS] we display the ground-state density n^ on the correlated site for all values of the on-site energy, Eq, for the den- sity functional BALDA and the many-body HF, 2B, GW, and T-Matrix approximations. For Eq = — C//2 the sys- tem is invariant under the particle-hole transformation dja — ?► {—ydlj^ and therefore the exact density on the impurity site equals Ug = 1. This remains valid in all the approximation schemes employed. If we increase the gate potential Eq away from the particle-hole symmetric point, the density on the impurity site decreases almost linearly in all approximations. In order to enhance the differ- ences between the approximations, in the bottom panel we plot ri(eo) — niin(eo) where riiin(£o) = o,£o + b and the constants a and b are chosen such that niin(— [//2) = 1 and niin{U/2) = 0.35. In the vicinity of the particle-hole symmetric point, the BALDA has a cusp that is respon- sible for correlation induced density fluctuations on the impurity site. This gives a time-dependent description of the Coulomb blockade^. The HF approximation can describe the Coulomb blockade provided we allow the spin symmetry to be broken. The many-body approx- imations that we use here do not seem to be able to describe the Coulomb blockade without spin-symmetry breakingiS although the onset of the Coulomb blockade is observed^. It can be concluded from the above obser- vations that BALDA yields the Coulomb blockade with- out spin symmetry breaking^. For Eq < — C//4 the XC- potential is close to zero and BALDA consequently differs III. TRANSPORT THROUGH A WEAKLY COUPLED CORRELATED SITE We perform many-body and density-functional trans- port calculations for the Anderson impurity model. The one-site model is fully specified by three parameters: the Hubbard interaction (or charging energy) U, the on-site energy Eq and the hopping Viink connecting the inter- acting impurity site to leads. The leads on-site ener- gies are El — £r = and the hopping in the left and right lead Vl ^ Vjj ^ V. All parameters are given in units of the lead hopping V. For times i < the contacted system is in equilibrium at zero temperature (/i = Ef) and Fermi energy ef- A constant bias Wa in lead a = (L,R) is suddenly switched on at i = after which the time-dependent observables are calculated. We only consider weak coupling to the leads, i.e., Mink <C V, since in this regime the role of correlation effects is en- hanced. The equilibrium Green function is obtained as the self-consistent solution of the Dyson equation^ for different approximate many-body self-energies. In the TDDFT calculations the initial state is obtained by a self- consistent static DFT calculatiorP. For the XC-potential we use the modified BALDA defined in Section fllBl it=«; 0.9 0.8 ,0.7 0.6 0.5 0.4 0.3 Oi HF — 2B B--E1 GW T-Matrix i BALDA H 1 1 1 1 1 1 1 1 1 1 1 1 1 1 i U/2 FIG. 3: Ground-state density Ug on the correlated site versus the on-site energy, eo, for [7 = 1, Mink = 0.5 and e^ = 0. In the bottom panel we subtracted niin(eo) = aeo -I- 6 in order to enhance the difference between the curves. The constants a and b are such that niin(— C//2) = 1 and niin{U/2) — 0.35. 0.5 I 1.5 W/U FIG. 4: Steady-state density n^ (left) and current / (right) for a symmetrically applied bias Wl ~ —Wr — W/2 and for three different values of the on-site energy eo- The rest of the parameters are U — 1, Viink = 0.5 and ef = 0. FIG. 5: Steady-state density n^ (left) and current / (right) for a asymmetrically applied bias Wl = W, Wr — and for three different values of the on-site energy eo- The rest of the parameters are f/ = 1, Viink = 0.5 and sf = 0. substantially from the correlated MBPT results and fol- lows more closely the HF curve. When £q attains positive values, the correlation potential is large and negative, fa- voring charge accumulation (see Fig. [2]). Consequently, the BALDA deviates from HF and follows the correlated MBPT results, in particular with the GW results for Eq around U/2. As a general feature, we find that corre- lations favor the presence of electrons on the interacting site, since the density in the BALDA and the many-body approaches is larger than the HF density for all values of the on-site energy. B. Nonequilibrium steady-state results We now shift our attention to the nonequilibrium case. In the left panel of Fig. |4] we display the steady-state density and current (within ABALDA, HF, and 2B) for a symmetrically applied bias Wl — —Wr — W/2 and for three different values of the on-site energy Eq = — U/2 , , U/2. To improve the clarity of the plot we do not display the results for GW and T-Matrix as they are, in this parameter range, in close agreement to those ob- tained within 2B. In the left panel of Fig. |4] we see that the 2B, HF, and ABALDA densities are generally in good agreement with each other. For the corresponding steady-state current, bench- mark results are available from tDMRG calculations (see Ref. [M| . In the right panel of Fig. [4] we plot the cur- rents as a function of the bias W/U. Because the current is proportional to the overlap of the energy bands of the leads, for higher biases, i.e., W/U > 1.5, the steady-state current decrease with increasing bias. We note that for small bias values all the approxima- tions yield values for the current which are on top of the numerically exact tDMRG results for all on-site en- ergies considered. However, for higher biases, only the current obtained within 2B follows closely the tDMRG values for all on-site energies. Therefore, in this range of parameters, we will use the 2B results for benchmark- ing the other approximations. For Eq = U/2, the HF and ABALDA results follow closely the tDMRG and 2B curves, and for the whole bias-range. For higher biases and smaller on-site energies, i.e. Eq = and Eq — —U/2, they considerably overestimate the exact results. How- ever, the conductances, i.e. the initial slopes of the I- V curves in Fig. [4j still remain in close agreement with the 2B approximation and the tDMRG approach. This agrees with the Friedel sum rule that relates the conduc- tance to the densitySS. For the results displayed in Fig. [5] we considered the same system parameters and plotted the density (left panel) and the current (right panel) for an asymmet- rically applied bias Wl — W, Wr — 0. The overlap between the lead energy bands starts to decrease for W/U > 1 and, consequently, the currents decrease with increasing bias. The steady-state densities behave sim- ilarly to the case of symmetric biases (see Fig. [4]): the ABALDA and the HF results are in agreement with 2B results except for the case of gate potential Eq = at high bias. For the steady-state current (left panel) we also see the same trend: the ABALDA results are close to the HF results and overestimate the 2B results. We observe the same trends as in the case of symmetric bias (see Fig. [4]) which indicates that the 2B approximation also here gives a description close to the exact result. W = 0.4 W=1.0 HF 2B Q--H GW T-Matrix ABALDA o — o tDMRG U = 0.6 U= 1.4 o FIG. 6: Transient currents for different values of the applied bias Wl = -Wr = W/2, U = 0.5 and Vnnk = 0.3535. In the upper panels, eo = —U/2 corresponds to the particle-hole symmetric point. In the lower panels Eq = U 12. C. Time-dependent results: Adiabatic effects 50 40 30 FIG. 7: (Color online) Time-dependent density no(t) for a system with Fermi energy ef = 0, and eo = 0.2, Viink = 0.2 and for different values of the charging energy t/ = 0.6 (left column), 1.4 (right column). The system is driven out of equilibrium by an external bias Wl ~ 0.4 and Wr — 0. The constant Go — e^ /{2-Kh) = 1/{2-k) is the quantum of conduction in atomic units. We novif study the performance of the different approx- imations in the description of transient phenomena. The results are compared to the numerically exact tDMRG data of Ref. HH obtained for a lead-impurity hopping parameter Vlink = 0.3535. This decrease in the hop- ping parameter amounts to a slight enhancement of cor- relations as compared to the steady-state results of the previous section. The tDMRG calculations^ were done for a particle-hole symmetric situation with £o = —U/2. In addition, we compare the many-body results with ABALDA for the on-site potential eo = U/2, which is away from the discontinuity of the Wxc- In the upper panels of Fig. |6] we display the transient currents as a function of time for the various many-body approaches and the ABALDA as compared to the bench- mark tDMRG data. Since the tDMRG calculations are performed on finite systems, one sees the influence of re- flections at the system boundaries after a sufficiently long propagation time. The many-body results beyond HF are all in good agreement with the tDMRG results, the most accurate one being the 2B approximation. Not only the values of the steady-state current but also the charac- teristic bump in the transient is well-reproduced. The ABALDA and the HF approximations perform very sim- ilarly; they overestimate the values of the steady-state current and for a bias value of W^ = 0.4 they under- estimate the height of the transient bump. Also the many-body approximations underestimate the height of the bump somewhat. However, the best agreement is again found for the 2B approximation. It is difficult to pinpoint the origin of the different behavior of the tran- sient bump in the ABALDA and HF when compared to results obtained within correlated approximations. It is worth emphasizing, however, that in time-local approx- imations such as HF the terms responsible for the ini- tial correlation in the current formula of Eq. Q, i.e., the terms with components on the vertical track of the Keldysh contour, are lacking. In general these terms lead to damping and hence time-local approaches such as HF tend to overshoot the bump in the transient currentP. In the upper left panel of Fig. [6] such overshoots for the HF and ABALDA are probably masked by the fact that the final steady state current goes to a value that is too large. We finally like to point out that in systems with more levels the transient structure has a more rich oscillatory time-dependence which can be used to analyze the level structure of the central molecule^. In these cases the dif- ferences between the HF and the correlated approaches become more visible. In the lower panels of Fig. [6] we display the transient currents the on-site energy on the impurity site being Eq = U/2. The transients show a more pronounced os- cillatory behavior because of the increased energy-gap between the impurity level and the Fermi level of the right lead. This determines the oscillation frequency in the transient current (see Ref. 184). The many-body ap- proaches agree well with each other whereas the HF ap- proximation underestimates the value of the steady-state current for lower biases. In this case the ABALDA re- sults agree closely with the correlated many-body results. Due to the increased on-site energy Uxc becomes negative favoring charge accumulation on the Anderson impurity site (see Fig. [5] and Section III A I . In order to increase the effects of correlation we now reduce the hopping between the interacting site and leads to Vjink = 0.2 and consider two different charging energies U = 0.6 and U — 1.4. We also set eo — 0.2 and the Fermi energy to ej;- = 0. The system is driven out of equilibrium by a sudden switch-on of a constant, asymmetric bias Wl = 0.4, Wr = 0. In the upper row of Fig. [7]we show the time-dependent density for the interacting site. For U — 0.6, all re- sults obtained within correlated approximations are in close agreement to each other since if the interaction ap- proaches zero, all MBPT approximations become homol- ogous. As discussed before, in this regime ABALDA and HF are close to the MBPT approximations. By increas- ing the interaction the correlated MBPT approximations and ABALDA start to detach from HF. For stronger in- teractions, i.e., U = 1.4 (right panel-column), the HF density deviates considerably from the ABALDA and the many-body results. In the middle panel of Fig. [7] we show the time- dependent current through the right interface (from the interacting site to the lead) . As expected from the discus- sion in Section [HIB] the ABALDA systematically overes- timates the current given by tDMRG and 2B. The devi- ation from the 2B increases with increasing the interac- tion. The GW approximation also shows a smaller but noticeable deviation from the 2B approximation. The agreement between the ABALDA and the many-body results deteriorates gradually with an even further in- crease of the charging energy. For the MBPT results, the differences in the currents when increasing the in- teraction can be explained with the aid of the spectral function. We display the steady-state spectral functions in the lower panel of the Fig. [7| Since the current is ap- proximately proportional to the integral of the spectral function over the bias window (see Eq. ([lO|), the highest current is given by the approximation which has most spectral weight inside the bias window. On the other hand, the ABALDA spectral function being very close to the HF spectral function does not explain the rather large overestimation of the ABALDA current. As in the case of the site densities, for small charging energies the spectral functions of all the approximations remain very close to each other. The spectral functions of correlated MBPT approxima- tions are broadened compared to the HF spectral func- tions. This is because many-body interactions lead to a fast decay of many-body states generated by adding and removing particles. More precisely, the states |^(i)) = ^l/(*)l*o) and |$(i)) = d//(i)|*o) in which we add or remove a particle at time t to the impurity in the presence of a bias have decreased survival probabilities |(*(i)|*(t'))l' and |($(i)|$(t'))|' for \t-t'\ -^ c3o when we include interactions. This process is often referred to as quasi-particle scattering. When the charging energy is increased quasi-particle scattering broadens the spectral functions and lowers the intensity of the spec tral pe ak in the case of correlated MBPT approximation&^*^2I5II_ The broadening of the HF spectral function is indepen- dent of U due to the absence of quasi-particle scattering and depends only on the embedding to the leads. The same holds true for the ABLADA spectral function which remains very close to the HF spectral function when in- creasing the interaction. It should be noted, however, that the ABALDA spectral function is the one of the KS system and should not be regarded as an approxi- mation to the true spectral function. The clear broaden- ing of the MBPT spectral functions as compared to HF demonstrates the importance of non-adiabatic effects in the transient regime. Therefore, memory must be taken into account for a proper description of ultrafast time- dependent processes. In the next section we show that memory is, however, not enough to improve the results of the steady-state current and we identify a second impor- tant direction in which to go to improve the ABALDA. D. Time-dependent lead densities and non-locality In order to gain some insight on how to cure the de- ficiencies of the ABALDA XC-potential, so as to yield an improved time-dependent current, we argue as fol- lows: In equilibrium, the density deep inside the leads is the same in all approximations and it is uniquely de- termined by the Fermi energy ep- Let us denote with Ug (5=ground state) the density at a site with index j^ deep inside, say, the right lead, such that rij = Ug for all 0.515 0.51 FIG. 8: Time-dependent density in the right lead within the 2B approximation for a system with Fermi energy ef = 0, and Eo ~ 0.2, Viink = 0.2 and U = 0.6). The system is driven out of equilibrium by an external bias Wl ~ 0.4 and Wr — 0. A density wave entering the lead can clearly be observed. j > jd- If we plot the current Id{t) to the right of jd^ no difference will be observed in the site-density until after a time td = jd/v, where v is the velocity of the density wave-front moving into the right lead. This is clearly illustrated in Fig. [8] where we show the time-dependent lead densities obtained from a 2B calculation at interac- tion strenght U = 0.6 for the first 20 sites in the right lead. In the lower side of the figure we clearly see a wave front moving into the right lead. Let us then consider an interval of the right lead that extends from jd to jd + Nd with Nd ^ 1. In equilibrium the number of electrons in this interval is simply UgNd- At the time t ^ td the current wave-front reaches the site jd, enters inside the interval (jd,Jd + Nd) and after a time Td = Nd/v it goes out through the site jd + Nd- For times i > t^ + Td an equal amount of electrons en- ters in and exits from the interval, and a local steady-state is reached. The number of electrons in the considered in- terval is then given by n.Nd rigNd td+Ti dtld{t) UgNd + IsTd, (20) with Is the value of the steady-state current. Taking into account that Td — Nd/v we conclude that the steady- state density deep inside the leads must be = ng+ IJv. (21) From Fig. |8] we see that for our 2B calculation the ve- locity V has the value v = 1.88. Given the value of the current of Is — 0.034 for this case {U — 0.6) we find that the density difference Ug is approximately 0.018 which is in good agreement with the value in the upper panel of Fig. |9] Also for the case of the C/ = 1.4 interac- tion strength we see from the lower panel of Fig. [9] that a' 0.025 0.02 0.015 -^ ■ ■♦ ' 1 ,„.■■;:■--;•;;-:;■■.... 1 ■ ♦ .,,1.,.,^ ♦■■■■"■ ■■■,.■ ""■■A-"' '"-■A A ▲ HF • •2B ■■.. ■ y' • '• ■ a GW - ABALDA 1 1 1 V - .-■♦■-. ■■-♦■■■"■" ...,.-♦ A. A-,.-. ■;*:■.... -■■"-- '"-A ■■■■" "■■■•■■■■' ▲ ■■■• " 1 1 1 1 0.02 ? 0.015 a" 0.01 0.005 12 3 4 5 Site number FIG. 9: Difference between steady-state and ground-state density in the right lead for a system with Fermi energy Ef = 0, and eo ~ 0.2, l^ink = 0.2 and for different values of the charging energy U = 0.6 (top panel) and U — lA (bot- tom panel). The system is driven out of equilibrium by an external bias Wl =0.4 and Wr — 0. the ratio of the density differences Ug n, for ABALDA and 2B is the same as the corresponding ratio for the currents in Fig. [7] We note that the value v is close to the Fermi velocity in the lead at half-filling as ob- tained from a semi-classical calculation. This is given by V = 2V . Equation (21) shows that if different approxi- mations yield different values of the steady-state current they must also yield different values of the steady-state density deep inside the leads. This can indeed be seen in Fig. [91 where we plot Us — Ug for the various approx- imations for the first five sites in the right lead. The ordering of the density differences is identical to that of the currents in Fig. [7] Therefore, the ABALDA overes- timates the difference between steady-state and ground state densities in the leads. However, ABALDA gives a quite good description of the density on the impurity site, comparable to those obtained within the many-body approximations. We thus conclude that the ABALDA xc-potential is quite accurate on the impurity site but that setting the potential to zero in the leads is a too crude approximation. As was discussed in relation to the Sham-Schliiter equation (see Eq. 19) the XC-potential will in general have values in the leads even when the in- teraction is localized on the impurity site only. Hence, in order to obtain accurate values for the current within a TDDFT approach one needs an XC-potential that has a nonzero value in the leads. We wish to observe that this nonlocality is different in nature from the non-local de- pendence of the XC-potential on the density. The latter is already implied by the conclusions of the previous sec- 10 tion since non-locality in time and space are intimately related by conservation laws. IV. CONCLUSIONS AND OUTLOOK We study electron transport through an interacting Anderson impurity model within TDDFT and MBPT frameworks. Results obtained in the ground-state, tran- sient and steady-state regimes are compared with numer- ically exact tDMRG values. In the ground state, we find that for large values of the on-site energy, the density obtained using the ABALDA XC-functional is close to the densities obtained within correlated MBPT approximations. However, for smaller values of the on-site energy, the difference between the ABALDA and the correlated MBPT densities is signif- icant, ABALDA being closer to HF in this parameter range. In all the cases where benchmark tDMRG results are available we find that the MBPT approximations be- yond HF which we considered give densities and cur- rents close to the benchmark ones for the entire parame- ter range considered. This is true for both the transient and steady-state regimes. We find that in particular the 2B approximation performs very well. The transients obtained within the 2B approximation are the closest to the tDMRG ones, while the HF and ABALDA transients deviate significantly. This indicates that it is important to include memory or retardation effects to properly de- scribe quasi-particle scattering in non-equilbrium trans- port. Regarding the TDDFT approach we find that the ABALDA performs very well and yields accurate den- sities on the interacting site but in many cases overes- timates the steady-state currents. This problem can be linked to an overestimation of the lead-densities within the ABALDA. The results strongly suggest that it is nec- essary to go beyond the local approximation and that one especially needs to take into account XC-potentials that are nonlocal and that are non-zero within the leads. Improved functionals should therefore be nonlocal func- tionals in space. As has been clearly pointed out by Vignal^^, this implies that the functionals also need to be nonlocal in time in order to satisfy basic conserva- tion laws. The construction of such functionals is a clear challenge for the future. One way to proceed would be to make connections to many-body theory with conserving approximationa^. A. Acknowledgements Part of the calculations were performed at the CSC - IT Center for Science Ltd administered by the Min- istry of Education, Science and Culture, Finland. Stefan Kurth acknowledges funding by the "Grupos Consolida- dos UPV/EHU del Gobierno Vasco" (IT-319-07) and the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement No. 211956. Adrian Stan acknowledges funding by the Academy of Finland under grant no. 140327/2010. We acknowledge the support from the European Theoretical Spectroscopy Facility. ^ S. Datta, Electronics Transport in Mesoscopic Systems (Campridge University Press, New York, 1995). ^ M. Di Ventra, Electrical Transport in Nanoscale Systems (Cambridge University Press, 2008). ^ S. Gurviz and Y. Prager, Phys. Rev. B. 53, 15932 (1996). * G. Stefanucci and C.-O. Almbladh, Europhys. Lett. 67, 14 (2004). ^ V. Moldoveanu, V. Gudmundsson, and A. Manolescu, Phys. Rev. B 76, 165308 (2007). ® P. Pokes, F. Corsetti, and R. W. Godby, Phys. Rev. Lett. 101, 046402 (2008). ^ D. Ryndyk, R. Gutierrez, B. Song, and G. Cuniberti, Green function techniques in the treatment of quantum transport at the molecular scale, in Energy Transfer Dy- namics in Biomaterial Systems, edited by I. Burghardt, V. May, D. Micha, and E. Bittner, , Springer Series in Chemical Physics Vol. 93, p. 213, Springer- Verlag, 2009. ® P. Myohanen, A. Stan, G. Stefanucci, and R. van Leeuwen, Phys. Rev. B 80, 115107 (2009). ^ S. Andergassen, V. Meden, H. Schoeller, J. Splettstoesser, and M. Wegewijs, Nanotechnology 21, 272001 (2010). ^° X. Zheng et al, J. Chem. Phys. 133, 114101 (2010). " X. Zheng, F. Wang, C. Yam, Y. Mo, and G. Chen, Phys. Rev. B 75, 195127 (2007). ^^ M. von Friesen, C. Verdozzi, and C. Almbladh, Phys. Rev. Lett. 103, 176404. ^^ M. von Friesen, C. Verdozzi, and C. Almbladh, Phys. Rev. B 82, 155108 (2010). ^^ A. Brandschadel, G. Schneider, and P. Schmitteckert, An- nalen der Physik 522, 657 (2010). ^^ S. Kurth, G. Stefanucci, C.-O. Almbladh, A. Rubio, and E. K. U. Gross, Phys. Rev. B 72, 035308 (2005). ^® E. Boulat, H. Saleur, and P. Schmitteckert, Phys. Rev. Lett. 101, 140601 (2008). ^"^ P. Myohanen, A. Stan, G. Stefanucci, and R. van Leeuwen, Europhys. Lett. 84, 67001 (2008). ^* X. Wang, C. Spataru, M. Hybertsen, and A. MiUis, Phys. Rev. B 77, 045119 (2008). ^^ S. Schmitt and F. Anders, Phys. Rev. B 81, 165106 (2010). 2° U. SchoUwock, Rev. Mod. Phys. 77, 259 (2005). ^^ C. Karrasch et al, Eur. Phys. Lett. 90, 30003 (2010). ^^ C. Karrasch, M. Pletyukhov, L. Borda, and V. Meden, Phys. Rev. B 81, 125122 (2010). ^^ W. Metzner, M. Salmhofer, C. Honerkamp, V. Meden, and K. Schonhammer, arXiv:1105.5289v l [cond-mat.str- el] (2011). ■ ^"' F. Heidrich-Meisner, A. E. Feiguin, and E. Dagotto, Phys. Rev. B 79, 235336 (2009). 11 ^^ R. Dreizler and E. Gross, Density Functional Theory: An Approach to the Quantum Many-Body Problem (Springer- Verlag, New York, 1990). ^^ M. Marques et ai, Time- Dependent Density Functional Theory (Springer Heidelberg, 2006). ^■^ M. Di Ventra and T. Todorov, J.Phys: Condens.Matter 16, 8025 (2004). ^* F. Evers, F. Weigend, and M. Koentopp, Phys. Rev. B 69, 235411 (2004). ^^ R. Baer, S. I. T. Seideman, and D. Neuhauser, J. Chem. Phys. 120, 3387 (2004). ^° G. Stefanucci and C.-O. Almbladh, Phys. Rev. B 69, 195318 (2004). ^^ K. Burke, R. Car, and R. Gebauer, Phys. Rev. Lett. 94, 146803 (2005). ^^ N. Sai, N. Bushong, R. Hatcher, and M. Di Ventra, Phys. Rev. B 75, 115410 (2007). ^^ N. Sai, M. Zwolak, G. Vignale, and M. Di Ventra, Phys. Rev. Lett. 94, 186810 (2005). ^'^ G. Vignale and M. Di Ventra, Phys. Rev. B 79, 014201 (2009). ^^ G. Vignale and W. Kohn, Phys. Rev. Lett. 77, 2037 (1996). ^^ L.P. Kadanoffand G. Baym, Quantum Statistical Mechan- ics (Benjamin, New York, 1962). ^■^ P. Danielewicz, Ann. Phys. (New York) 152, 239 (1984). 3* G. Baym and L. Kadanoff, Phys. Rev. 124, 252 (1961). ^^ G. Baym, Phys. Rev. 127, 1391 ((1962)). *° K. Thygesen and A. Rubio, Phys. Rev. B 77, 115333 (2008). *^ M. Strange, C. Rostgaard, H. Hakkinen, and K. Thygesen, Phys. Rev. B 83, 115108 (2011). *2 M. Cini, Phys. Rev. B 22, 5887 (1980). *^ U. von Barth, N. Dahlen, R. van Leeuwen, and G. Ste- fanucci, Phys. Rev. B 72, 235109 (2005). '*'' M. Hellgren and U. von Barth, Phys. Rev. B 78, 115107 (2008). •^^ M. Hellgren and U. von Barth, J. Chem. Phys. 132, 044101 (2010). '"' P. W. Anderson, Phys. Rev. 124, 41 (1961). ^'^ R. van Leeuwen, N. E. Dahlen, G. Stefanucci, C.-O. Alm- bladh, and U. von Barth, Introduction to the keldysh for- malism, in Time-Dependent Density Functional Theory, edited by M. Marques et ai, , Lecture Notes in Physics Vol. 706, pp. 33-59, Springer- Verlag, 2006. ** N. E. Dahlen, R. van Leeuwen, and A. Stan, Journal of Physics, Conf.Ser. 35, 340 (2006). "^^ N. E. Dahlen and R. van Leeuwen, Phys. Rev. Lett. 98, 153004 (2007). ^° A. Stan, N. E. Dahlen, and R. van Leeuwen, J. Chem. Phys. 130, 224101 (2009). ^^ A. Stan, N. E. Dahlen, and R. van Leeuwen, J. Chem. Phys. 130, 114105 (2009). ^^ N. E. Dahlen and R. van Leeuwen, J. Chem. Phys. 122, 164102 (2007). ^^ L. Hedin, Phys. Rev. 139, 796 (1965). ^^ A. L. Fetter and J. D. Walecka, Quantum Theory of Many - Particle Systems (Dover Publications, 1971). ^^ Y. Meir and N.S. Wingreen, Phys. Rev. Lett. 68, 2512 (1992). ^^ N. Dahlen, A. Stan, and R. van Leeuwen, Journal of Physics, Conf.Ser. 35, 324 (2006). ^^ A.-M. Uimonen, E. Khosravi, G. Stefanucci, S. Kurth, R. van Leeuwen, and E.K.U. Gross, Journal of Physics, Conf.Ser. 220, 012018 (2010). ^* P. Myohanen, A. Stan, G. Stefanucci, and R. van Leeuwen, Journal of Physics, Conf.Ser. 220, 012017 (2010). ^^ K. Schonhammer, O. Gunnarsson, and R. M. Noack, Phys. Rev. B 52, 2504 (1995). ^° N. A. Lima, M. F. Silva, L. N. Oliveira, and K. Capelle, Phys. Rev. Lett. 90, 146402 (2003). ^1 C. Verdozzi, Phys. Rev. Lett. 101, 166401 (2008). ®^ S. Kurth, G. Stefanucci, E. Khosravi, C. Verdozzi, and E.K.U. Gross, Phys. Rev. Lett. 104, 236801 (2010). ®^ N. A. Lima, L. N. Oliveira, and K. Capelle, Europhys. Lett. 60, 601 (2002). ^'^ V. Franca, D. Vieira, and K. Capelle, |arXiv:1102.5018 (2011). ^^ R. van Leeuwen, Phys. Rev. Lett. 76, 3610 (1996). ^^ S. Schenk, P. Schwab, M. Dzierzawa, and U. Eckern, Phys. Rev. B 83, 115128 (2011). ®'^ H. Mera, K. Kaasbjerg, Y. M. Niquet, and G. Stefanucci, Phys. Rev. B 81, 035110 (2010).