Skip to main content

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).