Effective Equilibrium Description of Nonequilibrium Quantum Transport I: Fundamentals and Methodology
aa r X i v : . [ c ond - m a t . s t r- e l ] J u l Effective Equilibrium Description of Nonequilibrium Quantum Transport I:Fundamentals and Methodology
Prasenjit Dutt , Jens Koch , J. E. Han , Karyn Le Hur Departments of Physics and Applied Physics, Yale University, New Haven, CT 06520, USA and Department of Physics, State University of New York at Buffalo, Buffalo, NY 14260, USA (Dated: August 23, 2018)The theoretical description of strongly correlated quantum systems out of equilibrium presents several chal-lenges and a number of open questions persist. In this paper we focus on nonlinear electronic transport through aquantum dot maintained at finite bias using a concept introduced by Hershfield [Phys. Rev. Lett. , 2134 (1993)]whereby one can express such nonequilibrium quantum impurity models in terms of the system’s Lippmann-Schwinger operators. These scattering operators allow one to reformulate the nonequilibrium problem as aneffective equilibrium problem associated with a modified Hamiltonian, thus facilitating the implementation ofequilibrium many-body techniques. We provide an alternative derivation of the effective Hamiltonian of Hersh-field using the concept of an “open system”. Furthermore, we demonstrate the equivalence between observablescomputed using the Schwinger-Keldysh framework and the effective equilibrium approach. For the study oftransport, the non-equilibrium spectral function of the dot is identified as the quantity of principal interest andwe derive general expressions for the current (the Meir-Wingreen formula) and the charge occupation of the dot.We introduce a finite temperature formalism which is used as a tool for computing real time Green’s functions.In a companion paper we elucidate a generic scheme for perturbative calculations of interacting models, withparticular reference to the Anderson model. PACS numbers: 72.10.Bg, 73.63.Kv, 72.10.Fk
I. INTRODUCTION
Quantum systems exhibiting an interplay of interactionsand out-of-equilibrium effects are of significant interest andconstitute an active area of research.
In contrast to thesituation of equilibrium physics, however, there is currentlyno unifying theoretical framework for describing the dynam-ics of a generic quantum system out of equilibrium. Thereal-time Schwinger-Keldysh formalism has been suc-cessful in the treatment of specific systems, but occasion-ally gives rise to pathological perturbative expansions whichsuffer from infrared divergences. This can be viewed asan artifact of the infinite contour used in the integration,and one generally requires additional relaxation mechanismsto regulate the theory. While quantum impurity modelsin equilibrium have been extensively studied, probing andunderstanding properties of nonequilibrium steady states isa far more subtle task, many aspects of which are yet tobe explored. Nevertheless, there have been significant ad-vances in our understanding via several distinct approaches,which include the scattering Bethe Ansatz, field the-ory techniques, time-dependent density matrix renormal-ization group (RG), time-dependent numerical RG, perturbative RG, Hamiltonian flow equations, func-tional RG, strong-coupling expansions, diagram-matic Monte Carlo, and imaginary-time nonequilibriumquantum Monte Carlo. Mesoscopic quantum objects such as quantum dots (“arti-ficial atoms”), characterized by a set of discrete energy lev-els, in contact with reservoir leads, are described by quan-tum impurity models. In this context, the Anderson model,which describes a single discrete level coupled via tunnelingto a Fermi-liquid sea, is of particular interest. For a levelwith sufficiently low energy and with strong on-site interac- tion, the Anderson model mimics the situation of a Coulomb-blockaded quantum dot. In this regime, charge degrees offreedom on the level are frozen out and the Anderson modelbecomes intimately related to the Kondo model. The latterdescribes a magnetic impurity (manifested by the spin of thehighest occupied level for an odd number of electrons on thequantum dot) entangled with the spins of a Fermi sea. (Inthe mesoscopic setting, the Fermi sea is replaced by the reser-voir leads.) The Kondo entanglement then produces a promi-nent Abrikosov-Suhl resonance in the density of states of thequantum dot at the Fermi level, which results in perfect trans-parency when the quantum dot is symmetrically coupled toits leads (the source and the drain).
Nanoscale systemscan be routinely driven out of equilibrium by applying a biasvoltage between the two reservoir leads. Among the variousissues which arise out of equilibrium, the precise fate of theAbrikosov-Suhl resonance when applying a finite bias volt-age, remains a delicate issue. In general, it seems essential toelaborate theoretical and numerical methods which will allowaccess to the full current-voltage characteristics of interactingmesoscopic systems.The present paper constitutes Part I of a two-part series.Its main purpose is to reformulate electronic nonequilibriumtransport in quantum impurity models in terms of an effectiveequilibrium steady-state density matrix. This concept was ini-tiated by Hershfield, who proposed that the appropriatedensity matrix in Boltzmann form could be constructed ex-plicitly by invoking the Lippmann-Schwinger operators of thesystem. Recently, there have been several numerical efforts atimplementing this scheme . However, the theoreticalfoundations of this approach were yet to be firmly established.In this work, we concretely establish the validity of this effec-tive equilibrium approach by making use of the open systemlimit, which guarantees that the system relaxes to a steadystate. This well-controlled mathematical contrivance mimicsthe role of relaxation mechanisms, and does not require theexplicit inclusion of the bath degrees of freedom at the levelof the Hamiltonian. Furthermore, we demonstrate the equiv-alence between observables computed using the Schwinger-Keldysh framework and the effective equilibrium approach.It is important to emphasize the fact that the effective equi-librium description encompasses the case when one includesinteractions on the dot. We then propose an imaginary-timeformulation within this framework, which allows us to buildgeneral formulae and establish the mathematical machineryfor evaluating quantities relevant to transport, and to sys-tematically implement perturbative and non-perturbative tech-niques familiar from finite temperature equilibrium field the-ory. In Part II, we employ the formalism in the context ofthe nonequilibrium Anderson model and perform a systematicperturbative expansion in the repulsive interaction on the dot.Other aspects of the method, including numerical approaches,have been recently investigated in Refs. 21,23,24,40,41,53,54.In addition, we note that the non-interacting resonant levelmodel may also be a pertinent starting point in tackling sit-uations with strong interactions, as for example the large- N Anderson model or the Toulouse limit of the Kondomodel. Our paper is organized as follows. In Section II, we in-troduce the model and define the “open-system limit”.
This concept provides a transparent framework that guaran-tees the existence and uniqueness of a steady-state densitymatrix and we will use it recurrently in our discussion. InSection III, we present the effective equilibrium approach andjustify the validity of the method. We provide an alterna-tive derivation of the effective equilibrium density matrix pro-posed by Hershfield and establish its equivalence with theformal expression for the density matrix proposed by Doyonand Andrei . Furthermore, we present useful recursion rela-tions for the expectation values of observables which underpinthe equivalence between this approach and the time dependentprescription. In Section IV, we demonstrate the equivalencebetween observables computed within this framework and theSchwinger-Keldysh formalism which enables us to arrive atthe Meir-Wingreen formula for the current . We then de-rive a compact expression for the charge occupation on thedot, whereby we identify the spectral function of the dot as thecentral quantity of interest. It is instructive to note that this ef-fective equilibrium formalism avoids the (often complex) cou-pled Dyson’s equations for the various Green’s functions onthe Keldysh contour, and it is in principle possible to com-pute the Green’s function of interest directly. In Section V wedevelop the methodology appropriate for the effective equi-librium approach in terms of a finite temperature imaginary-time formalism, which can be used to analytically computereal time Green’s functions. II. MODEL
As a prototype for our discussion of Hershfield’s approachwe consider a system consisting of two Fermi-liquid leads,
PSfrag replacements ∼ k B T Φ µ − µ Γ t − t +1 H L H D H L ( α = −
1) ( α = +1) FIG. 1: Schematic setup (upper panel) and energy diagram (lowerpanel) of a generic quantum impurity model out of equilibrium. Thesystem is at a temperature T , and Φ = µ − µ − ≡ eV denotes thevoltage bias between the source and drain leads. The (bare) energyof the dot level is given by ǫ d . The energy broadening of this level,given by the width Γ , is due to the tunnel coupling between dot andleads. coupled by tunnel junctions to a central system with a numberof discrete levels, the “dot”, see Fig. 1. The Hamiltonian ofsuch a generic system can be written in the form H = H L + H D + H T . The first term H L = X αkσ ǫ αk c † αkσ c αkσ (1)describes the left and right leads ( α = ± ), where c αkσ ( c † αkσ ) annihilates (creates) an electron (strictly speaking aFermi-liquid quasiparticle) in state k with spin projection σ in lead α . The corresponding energy dispersion is denotedby ǫ αk . These leads couple to the dot, whose Hamiltonian isgiven by H D = X σ ǫ d d † σ d σ + H int (2)with d σ ( d † σ ) annihilating (creating) an electron with spin σ in the discrete level with energy ǫ d . Any electron interactionsare lumped into the contribution H int , and we will assume thatthese interactions are localized on the dot (as it is the case forthe Anderson model). Finally, the tunneling of electrons be-tween leads and dot is captured by the tunneling Hamiltonian H T = 1 √ Ω X αkσ t αk (cid:16) c † αkσ d σ + h.c. (cid:17) , (3)where Ω is the lead volume (assumed identical for both leads)and t αk specifies the tunneling matrix element for electrontransfer between state k in lead α and the discrete dot state. Inthe presence of a bias voltage, realized as a chemical potentialdifference Φ = µ − µ − between left and right lead, thetunneling induces an electric current.Typically, the steady state is reached after a short time de-termined by relevant relaxation rates in the leads. For the pur-pose of calculations aiming at steady-state quantities, it is con-venient to avoid the consideration of microscopic relaxationmechanisms and instead invoke the so-called “open-systemlimit”. In short, this approach proceeds as follows: Ini-tially, up to some time t = t < in the early past, the tunnel-ing term is absent and the Hamiltonian of the system is givenby H = H L + H D , such that leads and dot are decoupled.For t < t the system is hence described by the separabledensity matrix ρ = exp " − β H − Φ2 X α αN α ! . (4)Here, β = ( k B T ) − denotes the inverse temperature and N α = P kσ c † αkσ c αkσ is the number operator of electrons inlead α . Between times t < t < , the tunneling is then‘switched on’ adiabatically, i.e. , H = H + H T e ηt θ ( t − t ) ,where the parameter η → + defines a slow switch-on rate.At time t = 0 the tunneling has reached its full strength, thesystem is in its steady state and observables can be evaluated.As demonstrated in Refs. 51 and 18, the existence anduniqueness of steady state is tied to the validity of the inequal-ities v F /L ≪ | t | − ≪ η , where v F denotes the Fermi veloc-ity and L the linear system size. Intuitively, these inequalitiesensure that hot electrons hopping onto a given lead at time t will not be reflected back and return to the junction be-fore the measurement process, and further that the process ofswitching on H T remains adiabatic. We also note that the en-ergy scale of switch-on | t | − suffices to smear out the energylevel spacing v F /L . In this sense, the openness of the sys-tem provides the “dissipation” mechanism necessary for thesteady state, allowing the high-energy electrons to escape toinfinity and thus, effectively relax. On a more technical level,the inequalities result in the factorization of long-time corre-lation functions, which facilitates the proof of existence anduniqueness of steady state.The crucial ingredient in our discussion is that L/v F deter-mines the largest time scale in our problem. The exact proto-col by which we switch on the tunneling is irrelevant for theformation of steady state. Let us take as an example the non-interacting case, and for the sake of convenience assume thatwe have identical leads, i.e. , t αk = t . Here, instead of adi-abatically turning on the tunneling suppose we do a quench,one observes that the transients decay with a relaxation time ∼ Γ , where Γ = 2 πt ν denotes the linewidth of the dot. Here ν symbolizes the density of states, which for simplicityis assumed to be a constant. To make the argument concrete,assume that the tunneling Hamiltonian is absent for t < ,and the distribution functions in the leads are given by Fermifunctions f ( ǫ − α Φ2 ) ( α = +1 , − specify the source anddrain reservoirs respectively). Furthermore, suppose that thedot is initially unoccupied, i.e. n in d = 0 denotes the initialcharge on the dot. The occupation of the dot at a time t( > ) is given by n d ( t ) = n d ss (cid:0) e − t (cid:1) − e − Γ t π × Z dω (cid:20) f ( ǫ − Φ2 ) + f ( ǫ + Φ2 ) (cid:21) cos [( ω − ǫ d ) t ]Γ + ( ω − ǫ d ) . (5)Here n d ss = Γ π Z dω f ( ǫ − Φ2 ) + f ( ǫ + Φ2 )Γ + ( ω − ǫ d ) (6)denotes the steady state expectation value of the dot occupa-tion, after the transients have decayed. It is interesting to notethat if we turn the coupling of the QD to the left and rightleads adiabatically to zero at the same rate, (formally this im-plies taking Γ → in Eq. (6)) then the final QD occupation isgiven by n fin d = f ( ǫ d − Φ2 ) + f ( ǫ d + Φ2 ) . (7)The fact that n in d = n fin d , clearly illustrates the irreversibilityof the turning on process. It encapsulates precisely how thesystem is driven out of equilibrium and is the reason we arerequired to analytically continue time to the Keldysh contour.This demonstrates that the steady state formalism cannot becontinued from zero to finite tunneling perturbatively. There-fore, in our effective equilibrium approach, we include thecontribution of the tunneling terms non-perturbatively. Inter-actions on the other hand can be turned on adiabatically. Thisprocedure is reversible and does not exhibit similar anomalies. III. DERIVATION OF THE EFFECTIVE EQUILIBRIUMFORM OF THE STEADY-STATE DENSITY MATRIX
The description of a nonequilibrium problem in terms ofan effective equilibrium density matrix was first proposed byHershfield in Ref. 49. To facilitate such a description, thenonequilibrium steady-state density matrix ρ is rewritten inthe usual Boltzmann form, ρ = exp [ − β ( H − Y )] , (8)at the cost of introducing a correction operator, which we willcall (following Hershfield’s convention) the Y operator. For-mally, the definition of this correction operator as Y = 1 β ln ρ + H (9)is always possible. However, in this form it is neitherparticularly elucidating nor useful for calculating nonequi-librium transport properties. Hershfield put forward theidea that Y can be expressed explicitly and compactly interms of Lippmann-Schwinger operators ψ αkσ , which arefermionic operators that diagonalize the full Hamiltonian, H = X αk ǫ αk ψ † αkσ ψ αkσ . (10)Hershfield proposed that the Y operator has the general form Y = Φ2 X αkσ αψ † αkσ ψ αkσ . (11)Since this operator encodes the entire Φ dependence, the Y operator is also called bias operator. Note that Y vanishes atzero bias and the steady-state density matrix correctly simpli-fies to the equilibrium density matrix.In this Section we provide a detailed proof of the repre-sentation of Y given in Eq. (11). In contrast to Hershfield,we do not invoke the presence of additional relaxation mech-anisms to reach steady state. Instead, we make systematicuse of the time-dependent open system approach, whichcircumvents ill-defined expressions and makes the proof rig-orous. A. Proof of the explicit form of the Y operator The structure of the proof is as follows. The starting point isthe expansion of the steady-state density matrix into a powerseries in the tunneling Hamiltonian H T . This can be accom-plished either by using the Hershfield form of the density ma-trix, or by employing the time-dependent framework of theopen-system limit where the steady-state density matrix is ob-tained by adiabatically switching on the coupling in the farpast. In both cases, one obtains analytical expressions for thepower series in H T . Order for order comparison of these ex-pansions then leads to a system of nested differential equa-tions for the bias operators. Finally, the solution to this systemof differential equations is then shown to be identical with therepresentation of Y in terms of Lippmann-Schwinger opera-tors, see Eq. (11). We should emphasize the fact that the ex-pansion in powers of H T is a purely formal procedure, whichwe adopt for the sake of a systematic comparison, and theproof below is non-perturbative in H T .
1. Expansion of ρ using the effective equilibrium representation In the derivation of the explicit form of the Y operator, it isconvenient to collect terms in orders of the tunneling Hamil-tonian H T . We thus start by expanding Y = P ∞ n =0 Y n intoa series in powers of the tunneling, such that Y n ∝ ( H T ) n .In the following, the index n will always be used for powercounting of the tunneling Hamiltonian H T . From Eq. (8) wethus obtain ρ = exp " − β ( ( H − Y ) + ( H T − Y ) − ∞ X n =2 Y n ) = exp " − β ∞ X n =0 X n . (12)Here, we have regrouped the Hamiltonian with the Y operatororder by order in the auxiliary operator X n ∼ ( H T ) n , which is hence defined as X n ≡ H − Y , n = 0 H T − Y , n = 1 − Y n , n ≥ . (13)One can expand the exponential operator to collect terms inpowers of the tunneling H T such that ρ = ∞ X l =0 ( − β ) l l ! ∞ X i ,...,i l =0 X i . . . X i l = ∞ X n =0 ∞ X l =0 ( − β ) l l ! X i + ··· + i l = n X i . . . X i l ≡ ∞ X n =0 ρ n , (14)where, following our general notation, ρ n denotes the order ( H T ) n contribution to the steady-state density matrix.
2. Expansion of ρ using the open-system approach Let us now employ the open-system approach (as outlinedin Section II). In this case, the steady-state density matrix isobtained by switching on the tunneling in the early past t < . Up to this time t , the leads are decoupled from the dot.The adiabatic switch-on of tunneling is facilitated by using H T e ηt θ ( t − t ) for the tunneling Hamiltonian (note that in thisscenario, the total Hamiltonian is therefore time dependent).At time t = 0 , transients have decayed and time evolutionhas turned the original density matrix ρ = ¯ ρ ( t = t ) intothe steady-state density matrix ρ = ¯ ρ ( t = 0) . We note thatthe condition of adiabaticity is strictly true only in the limit | t | ≪ η , where | t | → ∞ , which is assumed in the open-system limit. We emphasize that the actual evaluation of thislimit is deferred until the very end of all calculations. Keeping / | t | and η small but nonzero in the interim is crucial formathematical clarity and for avoiding ill-defined expressions.This time, it will be convenient to work in the interactionpicture (with respect to the tunneling), where in general O I ( t ) = e iH ( t − t ) O e − iH ( t − t ) . (15)denotes the interaction picture of the operator O . Note thatthe time t (not t = 0 ) has been chosen as the referencetime where Schr¨odinger and interaction pictures agree, O = O I ( t ) . [In this we differ from the conventions adoptedby Hershfield, who chose different reference times for theHeisenberg and interaction representation.]In the interaction picture, the density matrix satisfies theevolution equation i ddt ¯ ρ I ( t ) = [ H T,I ( t ) , ¯ ρ I ( t )] . (16)This is formally solved by ¯ ρ I ( t ) = ∞ X n =0 ( − i ) n n ! T (cid:26) n Y i =1 (cid:20)Z tt dt i (cid:21) (17) × [ H T,I ( t ) , [ H T,I ( t ) , [ . . . [ H T,I ( t n ) , ρ ]] . . . ] (cid:27) ≡ ∞ X n =0 ¯ ρ n,I ( t ) , where the n = 0 term simply fixes the boundary condition ¯ ρ ( t = t ) = ρ and T denotes time-ordering of operators.It should be noted that Eq. (17) already has the form of apower series in the tunneling and thus has been used to definethe n -th order contribution ¯ ρ n,I ( t ) to the interaction-picturedensity matrix. For later purposes it is useful to note that thecontributions can alternatively be obtained from ddt ¯ ρ n,I ( t ) = − i [ H T,I ( t ) , ¯ ρ n − ,I ( t )] , (18)i.e., a system of nested differential equations associated withthe boundary conditions ¯ ρ ,I ( t ) = ρ and ¯ ρ n,I ( t ) = 0 for n ≥ .
3. Derivation of differential equations for Y n For the results to be consistent, we require that the steady-state density matrix ρ be identical to the steady-state density matrix obtained in the open-system limit. This allows one todetermine the correct form of the Y operator, order for orderin the tunneling.We thus impose the identity of the steady-state density ma-trices e iH ( t − t ) ρ n e − iH ( t − t ) = ¯ ρ n,I ( t ) , (19)where we have transformed ρ n into the interaction picture.Eq. (19) is expected to hold for all t > , since at thatpoint the interaction has been fully switched on, and the time-independent Hamiltonian (used in Hershfield’s effective equi-librium approach) and the time-dependent Hamiltonian (fromthe open-system approach) are identical.We now differentiate the left-hand side of Eq. (19) with re-spect to time and use Eq. (14) to obtain ddt ρ n,I ( t ) = (20) ∞ X l =1 ( − β ) l l ! X i + ··· + i l = n l X k =1 X i ,I · · · dX i k ,I ( t ) dt · · · X i l ,I , where X i k ,I = e iH ( t − t ) X i k e − iH ( t − t ) denotes the interac-tion picture representation of X i k . Similarly, we may differen-tiate the right-hand side of Eq. (19). The resulting commuta-tor is given in Eq.(18) and contains ρ n − ,I ( t ) , for which wesubstitute the corresponding expression from Eq. (14). Thisway, we obtain ddt ρ n,I ( t ) = i [ ρ n − ,I ( t ) , H T,I ( t )] = i ∞ X l =1 ( − β ) l l ! X i + i + ... + i l = n − l X k =1 X i ,I ( t ) . . . [ X i k ,I ( t ) , H T,I ( t )] . . . X i l ,I ( t )= i ∞ X l =1 ( − β ) l l ! X i + i + ... + i l = n l X k =1 X i ,I ( t ) . . . [ X i k − ,I ( t ) , H T,I ( t )] . . . X i l ,I ( t ) , (21)where for n < we define X n = 0 . In the last step of Eq.(21), the summation constraint is shifted from n − to n to fa-cilitate the comparison with Eq. (20). This comparison yieldsthe relation ddt X n,I = i [ X n − ,I ( t ) , H T,I ( t )] . (22)Finally, utilizing the relation (13) between X and Y operators,one finds that the Y n operators also satisfy ddt Y n,I = i [ Y n − ,I ( t ) , H T,I ( t )] . (23)To obtain the Y operator from these differential equations, itis crucial to specify the boundary conditions at the initial time t = t . Given the relation e iH ( t − t ) e O e − iH ( t − t ) = e O I ( t ) , valid for any operator O , the boundary condition ¯ ρ ( t = t ) = ρ implies lim t ց t Y I ( t ) = Φ2 X α αN α . (24)It now remains to prove that the following interaction-pictureexpression of YY I ( t ) = Φ2 X αkσ αψ † αkσ,I ( t ) ψ αkσ,I ( t ) , (25)in terms of the Lippmann-Schwinger operators ψ † αkσ,I ( t ) represents the solution to the above initial-value problem.Once we have recapitulated the crucial properties of theseLippmann-Schwinger operators in the following subsection,it will be simple to confirm that this ansatz indeed solves thedifferential equation (23) subject to the boundary condition(24).
4. Properties of the Lippmann-Schwinger operators
It has been shown for the generic model defined in Sec-tion II that the Lippmann-Schwinger operators are fermionic, n ψ † αkσ , ψ α ′ k ′ σ ′ o = δ αα ′ δ kk ′ δ σσ ′ , (26)and diagonalize the full Hamiltonian (including all interac-tions), H = X αk ǫ αk ψ † αkσ ψ αkσ . (27)One can expand ψ † αkσ in powers of the tunneling Hamiltonian H T , i.e. , ψ † αkσ = ∞ X n =0 ψ † αkσ,n (28)such that ψ † αkσ,n ∝ ( H T ) n . In the interaction representation,the operators ψ † αkσ satisfy ddt ψ † αkσ,I ( t ) = i [ H , ψ † αkσ,I ( t )]= i [ H I ( t ) ,ψ † αkσ,I ( t )] − i [ H T,I ( t ) , ψ † αkσ,I ( t )] , (29)subject to the boundary condition lim t ց t ψ αkσ,I ( t ) = c αkσ . (30)Equation (29) can be further simplified, and collecting ordersof ( H T ) n , cast into the set of nested differential equations ddt ψ † αkσ,n,I ( t ) = iǫ αk ψ † αkσ,n,I ( t )+ i [ ψ † αkσ,n − ,I ( t ) , H T,I ( t )] . (31)As discussed in Ref. 53, the formal solution can be ex-pressed compactly as ψ † αkσ = c † αkσ + 1 ǫ αk − L + iη L T c † αkσ (32)in terms of Liouvillian superoperators L and L T . The ac-tion of such superoperators on any operator O is defined as L A O = [ A, O ] . Further, η → + is a regularization similarin spirit to the one utilized in the more familiar Lippmann-Schwinger states in scattering theory. From Eq. (32) one canwrite the detailed form of the operator ψ † αkσ,n ψ † αkσ = ∞ X n =0 (cid:20)(cid:18) ǫ αk − L + iη L T (cid:19) n c † αkσ (cid:21) ≡ X n ψ † αkσ,n . (33)Note that when the tunneling is set to zero, one finds indeedthat ψ † αkσ → c † αkσ .
5. Conclusion of the proof
It is now simple to verify that the ansatz for Y n,I ( t ) solvesthe set of differential equations given by Eq. (25), and obeysthe appropriate boundary condition (24): We start by writingthe ( H T ) n contribution to the Y operator as Y n,I ( t ) = X αkσ α Φ2 n X p =0 ψ † αkσ,p,I ( t ) ψ αkσ,n − p,I ( t ) . (34)Differentiating this with respect to time yields ddt Y n,I ( t ) = X αkσ α Φ2 n X p =0 (cid:20) (cid:18) ddt ψ † αkσ,p,I ( t ) (cid:19) ψ αkσ,n − p,I ( t )+ ψ † αkσ,p,I ( t ) (cid:18) ddt ψ αkσ,n − p,I ( t ) (cid:19) (cid:21) , (35)so that in conjunction with Eq. (29) one obtains ddt Y n,I ( t ) = X αkσ α Φ2 n − X p =0 [ ψ † αkσ,p,I ( t ) ψ αkσ,n − − p,I ( t )]= i [ Y n − ,I ( t ) , H t,I ] . (36)This concludes the proof.In summary, we find that the steady state dynamics of thesystem can be described by an effective equilibrium densitymatrix of the form ρ = e − β ( H − Y ) , (37)where the operator Y = Φ2 X αkσ αψ † αkσ ψ αkσ (38)encodes the entire nonequilibrium boundary condition of thesystem. B. Recursion relations for expectation values of observables
In addition to the previous proof, we follow Hershfield and underpin the equivalence of the adiabatic approach andthe effective equilibrium approach by showing that they leadto identical recursion relations for expectation values of ob-servables. Again, the systematic use of the open-system limitmakes the proof sound.Let O denote a generic observable. To derive the first recur-sion relation, we will decompose the steady-state expectationvalue hOi into a series counting the powers of the tunnelingHamiltonian H T . In the first step, we thus substitute the ex-pansion (14) of ρ , hOi = tr[ ρ O ]tr[ ρ ] = tr [ P ∞ n =0 ρ n O ]tr [ P ∞ m =0 ρ m ] . (39)Pulling out a factor of / tr[ ρ ] , expanding the denominatoras a geometric series, and collecting terms order for order in H T , we can rewrite this as hOi = ∞ X n =0 n X l =0 ( − l ∞ X ′ j ,...,j l =1 l Y s =1 (cid:20) tr[ ρ j s ]tr[ ρ ] (cid:21) tr[ ρ n − P ls =1 j s O ]tr[ ρ ] ≡ ∞ X n =0 hOi n , (40)where P ′ denotes a restricted summation, subject to the con-dition P ls =1 j s ≤ n . Eq. (40) allows one to prove the impor-tant recursion relation hOi n = tr[ ρ n O ]tr[ ρ ] − n X k =1 tr[ ρ k ]tr[ ρ ] hOi n − k , (41)which relates the n -th order term to the expectation value of O with respect to ρ n , and all lower-order terms hOi m ( m =0 , . . . , n − ).Now, we turn to the time-dependent representation assum-ing an adiabatic switch-on of the tunneling. The steady-stateexpectation value of an observable may now be obtained via hOi = tr[¯ ρ I (0) O I (0)]tr[ ρ ] , (42)where we have switched to the interaction picture. To arrive atthe desired recursion relation, we again decompose this into apower series of the tunneling Hamiltonian, hOi = tr[¯ ρ I (0) O I (0)]tr[ ρ ] = ∞ X n =0 hOi n , (43)where hOi n = ( − i ) n n ! 1tr ρ tr (cid:26) T n Y i =1 (cid:20)Z t =0 t dt i (cid:21) (44) × [ H T,I ( t ) , [ H T,I ( t ) , [ . . . [ H T,I ( t n ) , ρ ]] . . . ] O I (0) (cid:27) . This expression for hOi n allows one to derive the second re-cursion relation, which reads hOi n = tr[ ρ n,I (0) O I (0)]tr[ ρ ] − n X k =1 tr[ ρ k,I (0)]tr[ ρ ] hOi n − k . (45)The details of this derivation are given in Appendix A. Theproof makes explicit use of the factorization of two-time cor-relation functions in the limit of large time separation. Theagreement between Eqs. (41) and (45) underpins the equiva-lence between the time-dependent adiabatic approach and theeffective equilibrium approach.
IV. CORRESPONDENCE WITH THESCHWINGER-KELDYSH APPROACH
Having obtained the effective equilibrium form of thesteady state density matrix, we now illustrate how one can compute transport observables such as the current and thespectral function within this description. We propose an imag-inary time formulation for treating such nonequilibrium sys-tems in a manner similar to finite temperature field theory.In imaginary time, we define the propagation of an operatorby O ( τ ) = e τ ( H − Y ) O e − τ ( H − Y ) = e τ ( L−L Y ) O . (46)The nonequilibrium thermal Green’s function is defined on < τ < β as G O O ( τ ) = −hT [ O ( τ ) O (0)] i = −hO ( τ ) O (0) i . (47)Fourier transforming in imaginary time this results in G O O ( iω n ) = *( O , e iω n + iω n − L + L Y O )+ , (48)where ω n = (2 n + 1) π/β ( n ∈ Z ) denotes the fermionicMatsubara frequencies.Switching to real time, the Heisenberg representation of anoperator O is given by O ( t ) = e iHt O e − iHt . The nonequilib-rium real time retarded Green’s function can then be expressedas G ret O O ( t ) = − iθ ( t ) h{O ( t ) , O (0) }i = − iθ ( t ) tr (cid:2) e − β ( H − Y ) {O ( t ) , O (0) } (cid:3) tr (cid:2) e − β ( H − Y ) (cid:3) . (49)By using the spectral representation and then Fourier trans-forming, we obtain from this G ret O O ( ω ) = (cid:28)(cid:26) O , ω − L + iη O (cid:27)(cid:29) . (50)We emphasize that, despite the effective equilibrium characterof the Hershfield approach, there does remain one importantdifference between the effective equilibrium description andany regular equilibrium many-body theory. This distinctionarises from the different propagators in imaginary versus realtime; namely, in Hershfield’s effective equilibrium formalismall Heisenberg operators in imaginary time evolve under themodified Hamiltonian H − Y , whereas Heisenberg operatorsin real time evolve under the Hamiltonian H . As a result,imaginary time and real time Green’s function, G O O ( iω n ) and G ret O O ( ω ) are not simply related by a direct analytic con-tinuation iω n → ω + iη .As an observable of prime interest in transport, let us nowconsider the current flowing through the dot. It is given by I = I + I − − e (cid:28) d ( N ( t ) − N − ( t )) dt (cid:29) = − e X α α (cid:28) dN α ( t ) dt (cid:29) = i e X α α h [ N α ( t ) , H ] i = i X αkσ α et αk √ Ω D(cid:16) c † αkσ d σ − d † σ c αkσ (cid:17)E = i X αkσ α et αk √ Ω (cid:16) G d σ c † αkσ ( τ = 0) − G c αkσ d † σ ( τ = 0) (cid:17) = Im "X αkσ α et αk √ Ω G c αkσ d † σ ( τ = 0) . (51)where − e denotes the charge of the electron. Here, we haveused the fact that the system is in a steady state and the cur-rent is unchanged under time translations. As we shall seebelow, it will prove to be most convenient to express the cur-rent in terms of the Fourier representation of the imaginarytime Green’s function, namely I = Im " X αkσω n α et αk √ Ω 1 β G c αkσ d † σ ( iω n ) . (52)We now show the equivalence of this approach with theSchwinger-Keldysh formalism and recover the familiar Meir-Wingreen formula for the steady-state current. For sim- plicity, let us assume that the coupling to the dot is inde-pendent of k , i.e. t αk = t α , and that the leads are identical ǫ αk = ǫ k . (This is not a strict requirement and the proof canbe easily generalized.) Our starting point is Eq. (52) and wenow write explicitly β X ω n G c αkσ d † σ ( iω n ) = 1 β X ω n *( c αkσ , e iω n + iω n − L + L Y d † σ )+ . (53)Using Eq. (32) for c αkσ and evaluating the effect of L T , weobtain β X ω n G c αkσ d † σ ( iω n ) = 1 β X ω n (cid:20) *( ψ αkσ , e iω n + iω n − L + L Y d † σ )+ − *( ǫ k + L − iη d σ , e iω n + iω n − L + L Y d † σ )+ (cid:21) . (54)With a small trick, one can show that the second term doesnot need to be evaluated when computing the current I . Tosee this, we exploit the fact that in steady state, the currentin the left and right junctions must be identical and hence,we can write the current as a weighted average of the form I = (cid:0) t I − + t − I (cid:1) / (cid:0) t + t − (cid:1) . This eliminates the pres-ence of the second term in Eq. (54) from the expression forthe current. Proceeding with the remaining term we obtain inseveral steps, β X ω n *( ψ αkσ , e iω n + iω n − L + L Y d † σ )+ = 1 β X ω n *( e iω n + iω n + L − L Y ψ αkσ , d † σ )+ = 1 β X ω n *( e iω n + iω n − ǫ k + α Φ / ψ αkσ , d † σ )+ = f (cid:18) ǫ k − α Φ2 (cid:19) h{ ψ αkσ , d † σ }i = t α √ Ω f (cid:18) ǫ k − α Φ2 (cid:19) (cid:28)(cid:26) ǫ k + L − iη d σ , d † σ (cid:27)(cid:29) = t α √ Ω f (cid:18) ǫ k − α Φ2 (cid:19) G ret d σ d † σ ( ǫ k ) , (55)where the transition from the first to the second line is facil-itated by the general relation (B2). Also, in going from thesum over k to the integral over ǫ k we make the replacement P k → R ν dǫ k . Here we have followed the standard proce-dure and assumed the spectrum has been linearized. Further ν denotes the density of states, which is assumed to be a con-stant. Finally, this results in the well-known Meir-Wingreenexpression for the current, I =2 e Γ Γ − Γ + Γ − Z dǫ k A d ( ǫ k ) × (cid:20) f (cid:18) ǫ k + Φ2 (cid:19) − f (cid:18) ǫ k − Φ2 (cid:19)(cid:21) , (56) where A d ( ǫ k ) = − π X σ Im h G ret d σ d † σ ( ǫ k ) i (57)is the nonequilibrium spectral function of the dot and is gen-eral a function of the bias voltage Φ . We shall however sup-press this explicit bias dependence of the spectral function. Inthe above expression Γ α = πt α ν denotes the partial broaden-ing of the level due to the coupling to the lead α .Let us now derive the above identities via a slightly dif-ferent route and recover the traditional form of the MeirWingreen formula. Proceeding as in Eq. (54) one observes β X ω n G d σ c † αkσ ( iω n ) = 1 β X ω n (cid:20) *( d σ , e iω n + iω n − L + L Y ψ † αkσ )+ − *( e iω n + iω n + L − L Y d σ , ǫ k − L + iη d † σ )+ (cid:21) . (58)It is simple to show, in a sense similar to Eq. (55) that the 1st part of the above expression β X ω n *( d σ , e iω n + iω n − L + L Y ψ † αkσ )+ = t α √ Ω f (cid:18) ǫ k − α Φ2 (cid:19) G adv d σ d † σ ( ǫ k ) . (59)Thus, using Eq. (51) together with Eqs. (55) and (59) we get I = i X αkσ α et αk √ Ω (cid:16) G d σ c † αkσ ( τ = 0) − G c αkσ d † σ ( τ = 0) (cid:17) = i e π X ασ α Γ α Z ∞−∞ dǫ k (cid:26) f (cid:18) ǫ k − α Φ2 (cid:19) h G adv d σ d † σ ( ǫ k ) − G ret d σ d † σ ( ǫ k ) i − β X ω n (cid:20) *( e iω n + iω n + L − L Y d σ , ǫ k − L + iη d † σ )+ − *( ǫ k + L − iη d σ , e iω n + iω n − L + L Y d † σ )+ (cid:21)(cid:27) . (60)Using an identity from Appendix B, while integrating over ǫ k in the last 2 terms, we find that I = i e π X ασ α Γ α (cid:26) Z ∞−∞ dǫ k f (cid:18) ǫ k − α Φ2 (cid:19) h G adv d σ d † σ ( ǫ k ) − G ret d σ d † σ ( ǫ k ) i + 2 iπ G d σ d † σ ( τ = 0) (cid:27) . (61)Now iπ G d σ d † σ ( τ = 0) = 2 iπ (cid:10) d † σ d σ (cid:11) = 2 πG < ( t = 0) = R ∞−∞ dǫ k G < ( ǫ k ) . This helps us recover the traditional form of theMeir Wingreen formula: I = i e π X ασ α Γ α Z ∞−∞ dǫ k (cid:26) f (cid:18) ǫ k − α Φ2 (cid:19) h G adv d σ d † σ ( ǫ k ) − G ret d σ d † σ ( ǫ k ) i + G < ( ǫ k ) (cid:27) . (62)We proceed in a similar fashion to derive an expression forthe charge occupation on the dot n d = X σ h d † σ d σ i = 1 β X σ,ω n G d σ d † σ ( iω n )= 12 Z dǫ k A d ( ǫ k ) (cid:20) f (cid:18) ǫ k − Φ2 (cid:19) + f (cid:18) ǫ k + Φ2 (cid:19)(cid:21) = Z dǫ k A d ( ǫ k ) f eff ( ǫ k , Φ) , (63)where f eff ( ǫ k , Φ) = 12 (cid:20) f (cid:18) ǫ k − Φ2 (cid:19) + f (cid:18) ǫ k + Φ2 (cid:19)(cid:21) . (64)Note that on taking the limit Φ → , Eq.(63) reduces to aform well-known from equilibrium many-body theory. In Ap-pendix B we provide a straightforward derivation of this resultusing the spectral representation within the effective equilib-rium formulation. This result can also be obtained within theKeldysh framework in a slightly indirect manner as we haveshown in Appendix C. This perhaps accounts for its omis-sion in the existing literature, where the charge occupation istypically obtained using the lesser Green’s function G < . We thus identify the non-equilibrium spectral function ( ω, Φ) as the central quantity of interest, which one can useto calculate transport observables. In the effective equilib-rium scheme this can be done directly, without resorting tothe Dyson equations which couple the various Green’s func-tions, and can be complicated. It is instructive to note thatfor the non-interacting case the spectral function is indepen-dent of the bias, but the introduction of interactions generallyimparts to the spectral function a complex bias dependence. V. GENERATING FUNCTIONAL FRAMEWORK
As usual, the evolution in imaginary time of an oper-ator in Heisenberg picture is governed by the propagator e − τ ( H − Y ) , which closely matches the form of the density ma-trix e − β ( H − Y ) . This similarity makes the imaginary time de-scription very convenient for the evaluation of Green’s func-tions when using functional techniques. We will now outlinethe general methodology for constructing Green’s functionsusing the effective equilibrium density matrix within a func-tional framework.We construct the generating functional as a coherent state0functional integral over Grassmann variables, correspondingto the Lippmann-Schwinger operators Z [ J, J ∗ ] = Z Y αkσ ( D ψ ∗ αkσ D ψ αkσ ) e − S [ J,J ∗ ] , (65)where S denotes the action S = Z β dτ X αkσ (cid:20) ψ ∗ αkσ ( ∂ τ + ǫ αkσ − α Φ / ψ αkσ − (cid:0) J ∗ αkσ ψ αkσ + ψ ∗ αkσ J αkσ (cid:1)(cid:21) , (66)and for notational convenience, we have suppressed the factthat the Grassmann variables are functions of imaginary time.Any Green’s function in terms of the Lippmann-Schwingeroperators is now given by an appropriate functional derivativeof the generating functional G ψ α ′ k ′ σ ′ ψ † αkσ ( τ ) = −hT [ ψ αkσ ( τ ) ψ † α ′ k ′ σ ′ (0)] i = 1 Z δ δJ ∗ αkσ ( τ ) δJ α ′ k ′ σ ′ (0) Z ( J, J ∗ ) (cid:12)(cid:12)(cid:12)(cid:12) J =0 . (67)Here Z = Z (0) denotes the generating functional with thesource terms absent.We observe that the generating functional and consequentlythe Green’s functions have a trivial form in terms of theLippmann-Schwinger operators. However for transport wehave to compute Green’s functions in terms of the originaldegrees of freedom of the system such as G d σ d † σ and G c αkσ d † σ .The goal now is to express the c αkσ and the d σ in terms ofthese Lippmann-Schwinger operators. This is by no means atrivial task, since the precise form of the Lippmann-Schwingeroperators is in general not known. However for a systemgoverned by single particle dynamics, the exact form of thescattering states is readily calculated and the actual degreesof freedom of the system and these states bear a linear rela-tionship to each other. In this limit, it is straightforward topass from one set of operators to the other, and it will consti-tute the starting point of our calculations. Later this will laythe foundation for computing Green’s functions in interactingtheories, perturbatively in the interaction.For systems governed by single particle dynamics theLippmann-Schwinger operators are given as a linear combi-nation of lead electron and dot electron operators ψ † αkσ = X α ′ k ′ σ ′ Λ α ′ k ′ σ ′ c † α ′ k ′ σ ′ + X σ ′ κ σ ′ d † σ ′ , (68)where Λ α ′ k ′ σ ′ and κ σ ′ are appropriate c-numbers.For interacting theories, the form of the Lippmann-Schwinger operator turns out to be quite non-trivial and canbe written schematically as ψ † ς = X ′ υ i ,ξ j ,υ ′ k ,ξ ′ l Λ υ ...υ ′ ...ξ ...ξ ′ ... × Y i (cid:16) c † i (cid:17) υ i Y j (cid:16) d † j (cid:17) ξ j Y k ( c k ) υ ′ k Y l ( d l ) ξ ′ l , (69) where the summation P ′ is over υ i , ξ j , υ ′ k , ξ ′ l ∈ { , } withthe constraint P i υ i + P j ξ j − P k υ ′ k − P l ξ ′ l = 1 , so that intotal there is one more creation operator than there are anni-hilation operators. Here the coefficients Λ υ ...υ ′ ...ξ ...ξ ′ ... arec-numbers and we have used an abbreviated notation i, j, k, l for the labels of the single-particle lead and dot states. Thus, ina generic interacting model a simple way to pass between theLippmann-Schwinger basis and the original basis is in generalnot easily accomplished. The construction of the Lippmann-Schwinger operator order by order in the interaction is how-ever well defined and from this we can compute the Green’sfunctions in powers of the interaction strength, as shown inthe companion paper (part II). For the purposes of numer-ics, it is possible to circumvent explicit construction of theLippmann-Schwinger operator by working at the level of thedensity matrix. The functional formalism described above provides a con-venient platform for the implementation of non-perturbativetechniques. We are currently using this scheme for a large- N treatment of the infinite- U Anderson model out of equi-librium, following the slave-boson language of Read andNewns . The saddle point solution of the functional in-tegral in this case, mimics the non-interacting theory and it ispossible to systematically build corrections in powers of /N . VI. CONCLUSIONS
In this paper, we have provided a rigorous alternativederivation of the effective equilibrium density matrix ap-proach proposed by Hershfield, in the context of nonequilib-rium quantum transport using the concept of an “open sys-tem”. Furthermore, we have illustrated the methodology ofthis formulation for a quantum dot, in which the reservoirleads play the role of good (infinite) thermal baths, such that aunique steady state exists. The latter can be cast into the formof an effective equilibrium density matrix, where the associ-ated modified Hamiltonian can be explicitly written in termsof the Lippmann-Schwinger scattering state operators.We have detailed the foundations of the theory and demon-strated a rigorous correspondence between the observablescomputed using the Schwinger-Keldysh formalism and the ef-fective equilibrium approach. We emphasize that an advan-tage of the effective equilibrium approach is that it is con-ducive to the implementation of numerical methods, suchas numerical RG and Quantum Monte Carlo. Further-more, we have built a generating functional framework usingan imaginary-time formulation that will constitute the basicmathematical machinery for our computation of the steady-state dynamics of the system. In particular, we will de-velop an infrared-convergent systematic perturbation theoryscheme to treat interactions within the effective equilibriumapproach. This scheme enables a direct computation of thenon-equilibrium spectral function, from which we can in-fer transport observables, without resorting to (often compli-cated) coupled Dyson’s equations which arise in the Keldyshcontext.Finally, we underline that the non-interacting resonant level1model may serve as an appropriate starting point in thecase with strong interactions, e.g. , in the large- N Andersonmodel or in the close vicinity of the Toulouse limit inthe Kondo model. The effective equilibrium method mayalso be generalized to more sophisticated systems, such asone-dimensional mesoscopic systems embodied by the Lut-tinger paradigm where it is possible, in principle, to com-bine bosonization and the scattering picture. Acknowledgments
We acknowledge N. Andrei and H. Baranger for stimulat-ing discussions. This work is supported by the Department ofEnergy under grant DE-FG02-08ER46541 (P.D. and K.L.H.),by the Yale Center for Quantum Information Physics (NSFDMR-0653377, J.K. and K.L.H.) and by the NSF under grantDMR-0907150 (J.E.H.).
Appendix A: Derivation of the recursive form of O n In this appendix we provide the detailed derivation ofthe recursion relation (45) when working within the time-dependent approach, switching on the tunneling adiabatically.Our starting point for the derivation is the expression for theorder- ( H nT ) contribution to the steady-state expectation valueof some operator O , hOi n = ( − i ) n tr[ ρ ] tr (cid:26) Z t dt n Z t n t dt n − · · · Z t t dt [ H T,I ( t n ) , [ H T,I ( t n − ) , [ . . . [ H T,I ( t ) , ρ ]] . . . ] O I ( t = 0) (cid:27) , (A1)see Eq. (44) in the main text. Carrying out the t integrationand using Eq. (16), one obtains hOi n = ( − i ) n − tr[ ρ ] tr (cid:20) Z t dt n Z t n t dt n − · · · Z t t dt (A2) × [ H T,I ( t n ) , [ . . . [ H T,I ( t ) , ( ρ ,I ( t ) − ρ ,I ( t ))] . . . ] O I (0) (cid:21) . We now retain the term involving ρ ,I ( t ) , and proceed bycarrying out the t integration for the term containing ρ ,I ( t ) .This integration yields 2 terms: one with ρ ,I ( t ) which weretain, and the other with ρ ,I ( t ) which we subject to furtherintegration. Continuing in this fashion until all the variables t , . . . , t n in the latter term have been integrated out, we arriveat the result hOi n = tr [ ρ n,I (0) O I (0)]tr[ ρ ] (A3) − n X p =1 ( − i ) n − p tr[ ρ ] Z t n +1 =0 t dt n Z t n t dt n − · · · Z t p +2 t dt p +1 × tr (cid:26) [ H T,I ( t n ) , [ . . . , [ H T,I ( t p +1 ) , ρ p,I ( t )] . . . ] O I (0) (cid:27) . (A4) Let us focus on the p -th term of the above sum. One observesthat when p = n , none of the integrals are present and theargument of the trace is just ρ p,I ( t ) O I (0) . Suppressing theintegrals and other c-numbers we isolate the product of oper-ators within the trace tr (cid:26) [ H T,I ( t n ) , [ . . . , [ H T,I ( t p +1 ) , ρ p,I ( t )]] . . . ] O I (0) (cid:27) . We can cycle ρ p,I ( t ) in the above expression out of thenested commutators, thereby inverting the nested commuta-tor structure to re-express this as tr (cid:26) ρ p,I ( t )[ . . . [ O I (0) , H T,I ( t n )] , ...H T,I ( t p +1 )] (cid:27) . (A5)This allows us to rewrite hOi n in the form hOi n = tr [ ρ n,I (0) O I (0)]tr[ ρ ] (A6) − n X p =1 ( − i ) n − p tr[ ρ ] Z t dt n Z t n t dt n − · · · Z t p +2 t dt p +1 × tr (cid:20) ρ p,I ( t )[ . . . [ O I (0) , H T,I ( t n )] , . . . , H T,I ( t p +1 )] (cid:21) . To proceed further, let hh· · · ii = tr (cid:0) e − β ( H − Y ) · · · (cid:1) tr (cid:0) e − β ( H − Y ) (cid:1) , (A7)denote an expectation value with respect to the density ma-trix ρ . It has been shown by Doyon and Andrei in Ref. 51that, in the open system limit, long-time correlation functionsfactorize, i.e. hhO ( t ) O ( t ) ii → hhO ( t ) ii hhO ( t ) ii (A8)as | t − t | → ∞ . We will now employ this factorizationin Eq. (A6). Renaming the integration variables t ′ ≡ t p +1 , t ′ ≡ t p +2 , . . . , t ′ n − p ≡ t n and inserting ρ ρ − = hOi n = tr [ ρ n,I (0) O I (0)]tr[ ρ ] (A9) − n X p =1 ( − i ) n − p tr[ ρ ] Z t dt ′ n − p Z t ′ n − p t dt ′ n − p − · · · Z t ′ t dt ′ tr (cid:26) ρ ρ − ρ p,I ( t ) | {z } A ( t ) [[ . . . [ O I (0) , H T,I ( t ′ n − p )] , . . . ] , H T,I ( t ′ )] | {z } B ( t ′ ,...,t ′ n − p ,t =0) (cid:27) . We now recall that in the open system limit we take v F /L ≪ / | t | ≪ η with t → −∞ . This implies that the e ηt term in H T,I essentially cuts off the lower limit of the time integrals at a time ∼ /η much earlier than t . Thus, Eq. (A10) can be castinto an integral over hh A ( t ) B ( t ′ , . . . , t ′ n − p , t = 0) ii , where | t ′ − t | , . . . , | t ′ n − p − t | , | − t | → ∞ . Using Eq. (A8) we get hh A ( t ) B ( t ′ , . . . , t ′ n − p , t = 0) ii = hh A ( t ) ii hh B ( t ′ , . . . , t ′ n − p , t = 0) ii , (A10)and hence hOi n = tr [ ρ n,I (0) O I (0)]tr[ ρ ] (A11) − n X p =1 tr[ ρ p,I ( t )]tr[ ρ ] tr (cid:20)(cid:0) − i ) n − p Z t dt ′ n − p Z t ′ n − p t dt ′ n − p − · · · Z t ′ t dt ′ [[ . . . [ O I (0) , H T,I ( t ′ n − p )] . . . , ] H T,I ( t ′ )] ρ (cid:21) . Inverting the nested commutator structure as previously, we get hOi n = tr [ ρ n,I (0) O I (0)]tr[ ρ ] (A12) − n X p =1 tr[ ρ p,I ( t )]tr[ ρ ] tr (cid:20)(cid:0) − i ) n − p Z t dt ′ n − p Z t ′ n − p t dt ′ n − p − · · · Z t ′ t dt ′ [ H T,I ( t ′ n − p ) , [ . . . , [ H T,I ( t ′ ) , ρ ]] . . . ] O I (0) (cid:21) . Finally, comparing the second term in the latter equation withEq. (A1) and noting that tr[ ρ p,I ( t )] = tr[ ρ p,I (0)] , we arriveat the final recursion relation hOi n = tr [ ρ n,I (0) O I (0)]tr[ ρ ] − n X p =1 tr[ ρ p,I (0)]tr[ ρ ] hOi n − p . (A13) Appendix B: Derivation of n d using the spectral representation In this appendix we provide additional details on the equiv-alence between the effective equilibrium approach and theSchwinger-Keldysh formalism and explicitly show the iden-tity of the electron occupancy on the quantum dot as obtainedin the Schwinger-Keldysh formalism and as obtained with theeffective equilibrium approach. We start by stating two usefulidentities involving Liouvillian superoperators. The first iden-tity regards the imaginary-time Heisenberg representation ofan operator O and is given by O ( τ ) = e τ ( H − Y ) O e − τ ( H − Y ) = e τ ( L−L Y ) O . (B1)The second identity is a small trick which transfers a Liouvil-lian superoperator from one side of an anticommutator to the other, (cid:28)(cid:26) O , z − L O (cid:27)(cid:29) = (cid:28)(cid:26) z + L O , O (cid:27)(cid:29) , (B2)where z is a c-number.We now set out to demonstrate the equivalence of theSchwinger-Keldysh expression for the dot occupancy, n d = 12 X σ Z ∞−∞ dǫ A d σ ( ǫ ) (cid:20) f (cid:18) ǫ + Φ2 (cid:19) + f (cid:18) ǫ − Φ2 (cid:19)(cid:21) = 12 ν (0)Ω X σk A d σ ( ǫ k ) (cid:20) f (cid:18) ǫ k + Φ2 (cid:19) + f (cid:18) ǫ k − Φ2 (cid:19)(cid:21) (B3)and the expression obtained in the effective equilibrium ap-proach, n d = X σ h d † σ d σ i = X σ G d σ d † σ ( τ = 0) . (B4)It is crucial to recall the definition of the Green’s func-tions in the effective equilibrium approach (see Section IV).The Fourier representation of the (imaginary) time orderedGreen’s function is given by G O O ( iω n ) = *( O , e iω n + iω n − L + L Y O )+ , (B5)3and the retarded/advanced real time Green’s function in fre-quency space is obtained by G ret/adv O O ( ω ) = (cid:28)(cid:26) O , ω ∓ L ± iη O (cid:27)(cid:29) . (B6)To prove the equivalence between Eqs. (B3) and (B4), westart with the Keldysh expression and evaluate the spectralfunction A d σ ( ǫ ) = − π Im G ret d σ d † σ ( ǫ ) = 12 πi (cid:20) G adv d σ d † σ ( ǫ ) − G ret d σ d † σ ( ǫ ) (cid:21) . (B7)With Eq. (B6) we evaluate the retarded and advanced Green’sfunctions involved, G ret d σ d † σ ( ǫ k ) = (cid:28)(cid:26) d σ , ǫ k − L + iη d † σ (cid:27)(cid:29) = √ Ω t h{ d σ , ψ † αkσ }i , (B8) where, for simplicity, we are considering symmetric tunnelcouplings and have set t = t − = t . Similarly, one finds forthe advanced Green’s function G adv d σ d † σ ( ǫ k ) = (cid:28)(cid:26) d † σ , ǫ k − L + iη d σ (cid:27)(cid:29) = √ Ω t h{ d † σ , ψ αkσ }i . (B9)Utilizing these relations in the evaluation of the occupancy,we find n d = (B10) πi tν (0) √ Ω X αkσ (cid:0) h{ d † σ , ψ αkσ }i − h.c. (cid:1) f (cid:18) ǫ k − α Φ2 (cid:19) . Writing the Fermi function as a Matsubara sum, one finds n d = 14 πi tν (0) √ Ω 1 β X αkσω n e iω n + iω n − ǫ k + α Φ2 (cid:16) h{ d † σ , ψ αkσ }i − h{ d σ , ψ † αkσ }i (cid:17) = 14 πitν (0) √ Ω 1 β X αkσω n (cid:18) *( d † σ , e iω n + iω n + L − L Y ψ αkσ )+ − *( e iω n + iω n − L + L Y ψ † αkσ , d σ )+ (cid:19) = 14 πitν (0) √ Ω 1 β X αkσω n (cid:18) *( d † σ , e iω n + iω n + L − L Y c αkσ )+ + *( d † σ , e iω n + iω n + L − L Y t √ Ω 1 ǫ k + L − iη d σ )+ − *( e iω n + iω n − L + L Y c † αkσ , d σ )+ − *( e iω n + iω n − L + L Y t √ Ω 1 ǫ k − L + iη d † σ , d σ )+ (cid:19) . (B11)The first and third terms can be combined, and recalling Eq. (52) they can be shown to cancel, πi tν (0) √ Ω 1 β X αkσω n (cid:16) G d σ c † αkσ ( iω n ) − G c αkσ d † σ ( iω n ) (cid:17) = − πet ν (0) X α αI α = 0 , (B12)since in the steady state we must satisfy I = I − . Therefore, this results in: n d = 14 πi X αω n σ β Z ∞−∞ dǫ k (cid:18) *( d † σ , e iω n + iω n + L − L Y ǫ k + L − iη d σ )+ − *( e iω n + iω n − L + L Y ǫ k − L + iη d † σ , d σ )+ (cid:19) = 14 πi X αω n σ β Z ∞−∞ dǫ k (cid:18) *( e iω n + iω n − L + L Y d † σ , ǫ k + L − iη d σ )+ − *( e iω n + iω n − L + iη d † σ , ǫ k + L − L Y d σ )+ (cid:19) . (B13)In the next step, we will make use of the identity Z ∞−∞ dǫ ǫ ± L ∓ iη O = ± iπ O , (B14)and we briefly outline its proof. It is convenient to use thespectral representation for this purpose and thus enumerate the set of eigenstates of H by {| n i} . With this, one obtains Z ∞−∞ dǫ h m | ǫ ± L ∓ iη O| n i = Z ∞−∞ dǫ ǫ ± ( ǫ m − ǫ n ) ∓ iη h m |O| n i dǫ. (B15)4After shifting the integration variable ǫ → ǫ ± ( ǫ m − ǫ n ) , wecan further simplify the above expression to yield Z ∞−∞ dǫ ǫ ∓ iη h m |O| n i (B16) = (cid:18) P.V. (cid:20)Z ∞−∞ dǫǫ (cid:21) ± iπ (cid:19) h m |O| n i . Noting that the principal value integral vanishes concludes theproof of Eq. (B14). With this identity in place, Eq. (B13) canbe brought into its final form n d = 14 X αω n σ β (cid:20) *( d † σ , e iω n + iω n + L − L Y d σ )+ + *( e iω n + iω n − L + L Y d † σ , d σ )+ (cid:21) = 1 β X ω n σ G d σ d † σ ( iω n ) = X σ G d σ d † σ ( τ = 0) . (B17)This demonstrates that the Schwinger-Keldysh approach andthe effective equilibrium formulation reproduce the same ex-pression for the occupancy on the quantum dot. We note thatthe general expressions (53) and (B7) are particularly usefulsince they can be directly applied through numerical proce-dures such as numerical RG. Appendix C: Derivation of n d using the Keldysh approach The basic equation which will constitute the starting pointof our proof is given by Eq.(15) in Ref. 14. The current in lead α = ± for spin projection σ is given by I ασ ( t ) = − α e ~ Z t − t dt Z dǫ π Im (cid:26) e − iǫ ( t − t ) Γ α ( t, t ) × (cid:2) G <σσ ( t, t ) + f α ( ǫ ) G rσσ ( t, t ) (cid:3) (cid:27) , (C1)where we have restricted ourselves to a single level on the dothaving a spin index σ . The occupation number( n d σ ) of the dotin the state σ obeys the following differential equation: dn d σ ( t ) dt = 1 − e X α αI ασ . (C2)In steady state ( i.e. , t=0) we get dn d σ ( t ) dt (cid:12)(cid:12)(cid:12)(cid:12) t =0 = 0 . (C3)This implies X α Z −∞ dt Z dǫ π Im (cid:26) e − iǫt Γ α (0 , t ) (cid:2) G <σσ (0 , t ) + f α ( ǫ ) G rσσ (0 , t ) (cid:3) (cid:27) = 0 . (C4) Since in steady state there exists time translational invari-ance, we have G <σσ (0 , t ) = G <σσ (0 − t ) and G rσσ (0 , t ) = G rσσ (0 − t ) . Further we can set η → and t → −∞ , in thesense defined by the open system limit, right from the begin-ning. This in turn implies that Γ α (0 , t ) = Γ α = πν (0) t =Γ / . Let us focus on the first term in the expression above A = X α Z −∞ dt Z dǫ π Im (cid:26) e − iǫt Γ α G <σσ ( − t ) (cid:27) = − i X α Γ α Z −∞ Z dǫ π (cid:26) e − iǫt G <σσ ( − t ) − e iǫt (cid:2) G <σσ ( − t ) (cid:3) ∗ (cid:27) Using the fact [ G <σσ ( − t )] ∗ = − G <σσ ( t ) A = − i X α Γ α Z dǫ π G <σσ ( ǫ )= − i G <σσ (0)= 12 Γ (cid:10) d † (0) d (0) (cid:11) = 12 Γ n d σ . (C5)Similarly for the second term we get B = X α Z −∞ dt Z dǫ π Im (cid:26) e − iǫt f α ( ǫ )Γ α G rσσ ( − t ) (cid:27) = Γ2 X α Z dǫ π f α ( ǫ ) Im (cid:26) Z −∞ dt e iǫt G rσσ ( − t ) (cid:27) = Γ2 X α Z dǫ π f α ( ǫ ) Im (cid:26) G rσσ ( ǫ ) (cid:27) = Γ Z dǫ X α f α ( ǫ )2 | {z } f eff ( ǫ, Φ) π Im (cid:26) G rσσ ( ǫ ) (cid:27)| {z } − A dσ ( ǫ ) = −
12 Γ Z dǫf eff ( ǫ, Φ) A d σ ( ǫ ) . (C6)Combining the two terms and putting in Eq. (C4) we get n d σ = Z dǫf eff ( ǫ, Φ) A d σ ( ǫ ) , (C7)which implies that n d = X σ n d σ = Z dǫf eff ( ǫ, Φ) X σ A d σ ( ǫ )= Z dǫf eff ( ǫ, Φ) A d ( ǫ ) . (C8)This is identical to the expression obtained via the spectralrepresentation in Appendix B.5 R. M. Potok, I. G. Rau, H. Shtrikman, Y. Oreg, and D. Goldhaber-Gordon, Nature , 167 (2007). T. Delattre, C. Feuillet-Palma, L. G. Herrmann, P. Morfin, J.-M.Berroir, G. F`eve, B. Plac¸ais, D. C. Glattli, M.-S. Choi, C. Mora,and T. Kontos Nature Physics , 208 (2009). R. Egger, Nature Physics , 175 (2009). H. Steinberg, G. Barak, A. Yacoby, L. N. Pfeiffer, K. W. West,B. I. Halperin, and K. Le Hur, Nature Physics , 116 (2008). K. Le Hur, B. I. Halperin, and A. Yacoby, Annals of Physics (NY) , 3037 (2008). C. Mora and K. Le Hur, Nature Physics , 697 (2010). V. V. Deshpande, M. Bockrath, L. I. Glazman, and A. Yacoby,Nature , 209 (2010). Z. Zhong, N. M. Gabor, J. E. Sharping, A. L. Gaeta, and P. L.McEuen, Nature Nanotechnology , 201 (2008). T. Dirks, Y.-F. Chen, N. O. Birge, and N. Mason, Appl. Phys. Lett. , 192103 (2009). D. B. Gutman, Y. Gefen, and A. D. Mirlin, Phys. Rev. B ,085436 (2010). J. Schwinger, J. Math. Phys. , 407 (1961). L. V. Keldysh, Soviet Physics JETP , 1018 (1965). J. Rammer and H. Smith, Rev. Mod. Phys. , 323 (1986). A.-P. Jauho, N. S. Wingreen, and Y. Meir, Phys. Rev. B , 5528(1994). O. Parcollet and C. Hooley, Phys. Rev. B , 085315 (2002). A. C. Hewson,
The Kondo Problem to Heavy Fermions (Cam-bridge University Press, UK, 1997). P. Mehta and N. Andrei, Phys. Rev. Lett. , 216802 (2006). P. Mehta and N. Andrei, cond-mat/0702612 (2007). P. Fendley, A. W. W. Ludwig, and H. Saleur, Phys. Rev. Lett. ,3005 (1995). E. Boulat, H. Saleur, and P. Schmitteckert, Phys. Rev. Lett. ,140601 (2008). B. Doyon, Phys. Rev. Lett. , 076806 (2007). F. Heidrich-Meisner, A. E. Feiguin, and E. Dagotto, Phys. Rev. B , 235336 (2009). F. B. Anders, Phys. Rev. Lett. , 066804 (2008). S. Schmitt and F. B. Anders, Phys. Rev. B , 165106 (2010). A. Rosch, J. Paaske, J. Kroha, and P. W¨olfle, Phys. Rev. Lett. ,155603 (2003). C.-H. Chung, K. Le Hur, M. Vojta, and P. W¨olfle, Phys. Rev. Lett. , 216803 (2009). H. Schoeller, Eur. Phys. J. Special Topics , 179 (2009). S. Kehrein, Phys. Rev. Lett. , 056602 (2005). S. G. Jakobs, V. Meden, and H. Schoeller, Phys. Rev. Lett. ,150603 (2007). C. Karrasch, M. Pletyukhov, L. Borda, and V. Meden, Phys. Rev.B , 125122 (2010). H. Schmidt and P. W¨olfle, Ann. Phys. (Berlin) , 60 (2010). C. Mora, P. Vitushinsky, X. Leyronas, A. A. Clerk, and K. Le Hur,Phys. Rev. B , 155322 (2009). P. Vitushinsky, A. A. Clerk, and K. Le Hur, Phys. Rev. Lett. , 036603 (2008). C. Mora, X. Leyronas, and N. Regnault, Phys. Rev. Lett. ,036604 (2008). Z. Ratiani and A. Mitra, Phys. Rev. B , 24511 (2009). M. Schiro and M. Fabrizio, Phys. Rev. B , 155302 (2009). P. Werner, T. Oka, and A. J. Millis, Phys. Rev. B , 035320(2009). P. Werner, T. Oka, M. Eckstein, and A. J. Millis, Phys. Rev. B ,035108 (2010). L. Muehlbacher, D. F. Urban, and A. Komnik,arXiv:1007:1793(2010). J. E. Han and R. J. Heary, Phys. Rev. Lett. , 236808 (2007). J. E. Han, arXiv:1001.4989 (2010). M. Pustilnik and L. Glazman, in
Nanophysics: Coherence andTransport , edited by H. Bouchiat et al. (Elsevier, 2005), pp. 427–478. Single Charge Tunneling (eds. H. Grabert and M. H. Devoret(Plenum Press, New York), 1992). D. Goldhaber-Gordon, H. Shtrikman, D. Mahalu, D. Abusch-Magder, U. Meirav, and M. Kastner, Nature , 156 (1998). S. M. Cronenwett, T. H. Oosterkamp, and L. P. Kouwenhoven,Science , 540 (1998). J. Schmid, J. Weis, K. Eberl, and K. von Klitzing, Physica B , 182 (1998). M. Pletyukhov, D. Schuricht, and H. Schoeller, Phys. Rev. Lett. , 106801 (2010). R. Leturcq, L. Schmid, K. Ennslin, Y. Meir, D. C. Driscoll, andA. C. Gossard, Phys. Rev. Lett. , 126603 (2005). S. Hershfield, Phys. Rev. Lett. , 2134 (1993). A. Schiller and S. Hershfield, Phys. Rev. B , 14978 (1998). B. Doyon and N. Andrei, Phys. Rev. B , 245326 (2006). P. Dutt, J. Koch, J. E. Han, and K. Le Hur, arXiv:xxxx.yyyy(2010). J. E. Han, Phys. Rev. B , 125122 (2007). A. Oguri, Phys. Rev. B , 035302 (2007). N. Read and D. M. Newns, J. Phys. C: Solid State Phys. , L1055(1983). G. Toulouse, Phys. Rev. B , 270 (1970). Y. Meir and N. S. Wingreen, Phys. Rev. Lett. , 2512 (1992). Y. Meir, N. S. Wingreen, and P. A. Lee, Phys. Rev. Lett. , 3048(1991). T. L. Schmidt, P. Werner, L. M¨uhlbacher, and A. Komnik, Phys.Rev. B , 235110 (2008). M. Gell-Mann and M. L. Goldberger, Phys. Rev. , 398 (1953). B. A. Lippmann and J. Schwinger, Phys. Rev. , 469 (1950). Note that throughout the text, density matrices are not assumedto be normalized. Normalization is always introduced explicitlywhen evaluating expectation values. Note that due to the identity (19) we may, from here on, drop allbars on ρρ