Transport theory of coupled quantum dots based on auxiliary operator method

APS T ransp ort theory of coupled quan tum dots based on auxiliary op erator metho d Jung Hyun Oh and D. Ahn Institute of Quantum Inf ormation Pr o c essing and Systems, University of Se oul, Se oul 130-743, Kor e a Vladimir Bubanja Industrial R ese ar ch Ltd, PO Box 31-31 0, L ower Hutt, New Ze aland (Dated: June 1, 2018) W e f orm ulate the theory of electron transp ort through coup led-quantum dots by extending the auxiliary operator representation. By using the g enerating functional technique, w e derive the exact expressions for currents, dot-occup ation num b ers and spin co rrelations, and examine them based on the non-equ ilibrium Green’s function meth od under the non-crossing approximation (NCA). Our form ulation generalizes th e previous NCA approac hes by allowi ng full o ccupation numbers with a finite Coulom b repulsion. P A C S num b ers: 73.63.Kv,73.23.Hk,72.15.Qm I. INTRO DUCTION T r ansp ort pro p e rties of the double quantum-dot system have b een extensively s tudied both exper iment ally and theoretically .[1] This ar tificial molecule, analo gous to the t w o-impurity Anderson pro blem, pr ovides a go o d platform for examining the e x citing ph ysics of the correla ted electro n b ehavior, such a s the Ko ndo effect.[2 – 4] The double quantu m- dot structur e is also a fascinating sub ject fro m the p oint of view of p os sible a pplica tions in quantum co mputation, where it is suggested as a basic building blo ck, with qubits b eing represented by electron spins in ea ch quan tum dot.[5]. The r ich v ariety of correla ted electron phenomena o f a double quant um-dot system emerg es fro m the inclusion of the electron-electro n in teractions.[6 – 12] In this regar d, man y theoretica l approaches have been dev elop ed that concentrate on a limited range of relev ant para meter v alues , such as infinite C o ulomb repulsion, a finite in teraction but under equilibrium tr ansp ort co nditions, or symmetric dot o c c upations. Ho w ever, in order to describ e the control and measurement o f quantum bits in de ta il,[13, 14] it is necessar y to dev elop the theo r y that deals w ith correlated electron b ehavior in a wide range of in teraction para meters, level o ccupancies in dots, as well as to include the time-dep e ndent per turbations. The need for such a theor y comes from the fact that in exp erimental studies o f spin blo ck ade in la teral coupled quantum dots, independent tunnel barrier tuning with arbitra ry dot occupations hav e been achiev ed.[1 5, 16] Similarly , initialization and manipulatio n of quantum bits requires description of sudden c hanges of energy le vels in the dots due to the time v arying bias a nd gate voltages. In this pap er, we derive the expressions for the current , densities of states, do t occupancies , and spin co rrelations of the double- do t sys tem. O ur a pproach enables the treatment of ar bitrary Coulo mb in tera ctions, o ccupatio n num ber s, finite temperatur e , as w ell a s time v ar ying v oltages. T o do this, we ex tend the auxiliary o p erator representation and apply the non-eq uilibr ium Green’s functions metho d a sso ciated with the gener ating functional descr ibe d b y the non-cros s ing approximation.[17, 18] In o rder to establish the v alidit y of our appr o ach, we compar e o ur results, in a v ariety of situations, with the previous NCA and exact metho ds lik e the numerical renor ma lization group (NRG) scheme. W e find that o ur formulation r epro duces the previous NCA res ults, but deviates from the NR G metho d. This is not surprising since it is w ell known that the NCA fails in describing the lo w-energy F ermi-liquid regime. Since the vertex corrections cur e the low temper ature transpo rt pr o p erties,[19] this work may b e used for mor e inv olved further studies. The pa p er is organiz e d a s follows: In Sec. II we int ro duce the Hamiltonian of the double-dot system, and reformulate it in terms o f the auxiliar y particle o pe r ators in order to calculate the non-equilibrium Green’s function as so ciated with the generating functiona l. W e derive the ex pressions for ph ysical qua ntities via the relev a nt pro jection in the auxiliary particle o ccupation n um b er subspace. In Sec. I I I, the transpor t pro pe r ties of the double-dot system a re examined by using the n umer ical calculations, and are co mpared with the previous results. W e summar ize o ur main results in Sec. IV. Some mathematical deta ils a re deferr e d to the the appe ndices. 2 II. CALCULA TION METHOD A. Hamiltonian W e mo del the system, c onsisting of tw o quantum dots connected in series to the left and right electro des , by the Hamiltonian, H = H dots + H leads + H T . (1) T a king the full electr on-electron interaction into account,[20 ] the Hamiltonian of the coupled q uantum dots is g iven by , H dots = X ασ ( ǫ ασ n ασ + t H c † ασ c ¯ ασ + U α 2 n ασ n α ¯ σ + 1 4  U I − J 2  X σ ′ n ασ n ¯ ασ ′ ) − J ~ S L · ~ S R , (2) where w e assume that each dot ( α = L, R ) has energy lev els ǫ ασ lab eled with spin index ( σ = ↑ , ↓ ), and is coher ently coupled to the o ther ( ¯ α ) with the tunneling matrix elemen t t H . The dot num be r o p erators are given in terms of the cr eation (annihilation) operato rs c † ασ ( c ασ ) by n ασ = c † ασ c ασ , and the spin op er ators are g iven by ~ S α = 1 2 P σσ ′ ~ σ σσ ′ c † ασ c ασ ′ , where ~ σ ar e the Pauli matrices . U α and U I are the Coulom b in tera ction parameters for electr ons on the dot α and inter-dot, resp ectively , while J is the exchange coupling constant. The second and the third term in E q. (1) describ e the Hamiltonian o f the leads and tunneling betw een the dots and the adjacent leads; H leads = X kασ ǫ kασ a † kασ a kασ , H T = X kασ n T α kσ a † kασ c ασ + T α ∗ kσ c † ασ a kασ o , (3) where a † kασ ( a kασ ) c r eates (annihilates ) an electro n at the lea d α , and the constants T α kσ provide the coupling strength betw een the do t and the adjacent lead. W e as sume that the energy levels of ea ch dot ar e controlled independently by the nearby gate electro des[1] and the chemical potential µ α at the lead α is adjusted by applied voltage difference ∆ µ to b e µ L = µ 0 − ∆ µ/ 2 , µ R = µ 0 + ∆ µ/ 2 , (4) with resp e c t to the equilibrium chemical p otential µ 0 . In or der to take into account all the pos sible o ccupancies of the double-dot s y stem (ra nging from zero to four), we extend the idea of the auxiliary particle r epresentation [19, 2 1 ] a nd in tro duce the auxilia ry op er ators d † m ( d m ) as, c † ασ = X mm ′ ξ ασ mm ′ d † m d m ′ . (5) Here, d † m | v ac i is chosen as the m − th basis v ecto r diagonalizing the iso lated coupled-quantum dots (o r molecular states) as sp ecified in T able I. The a uxiliary opera tors satisfy the commutation relatio ns [ d m ′ , d † m ] ς = δ m,m ′ and [ d m ′ , d m ] ς = [ d † m ′ , d † m ] ς = 0 , where o dd (even)-n umbered sta tes are assumed to be fermionic (b osonic) and ς = + if bo th m and m ′ denote fermionic sta tes , otherwise ς = − . The overlap matrix ξ ασ mm ′ = h m | c † ασ | m ′ i has a non- zero v alue for the c ombination o f bo son-fermion or fermio n-b oson op erator s to ens ure orig ina l comm utation relation of [ c ασ , c † α ′ σ ′ ] = δ α,α ′ δ σ,σ ′ under the constr a int of Q = P m d † m d m = 1 . The pro o f is g iven in Appendix A. In terms of the auxiliary o per ators the Hamilto nia n o f the coupled dots is given by , H dots = 15 X m =0 ǫ m d † m d m + H int ( { d m ′ , d † m } ) . (6) The term H int represents the interaction b etw ee n auxilia ry particles (applicable for Q ≥ 2) and w e omit it hereafter since it do esn’t affect our final r esults. 3 T able I: Schemati c representation of the eige nstates for t w o coupled quantum dots with a coupling strength t H . Here, energy splitting for one- and three-particle states are giv en by ∆ 1 σ = p ( ǫ Rσ − ǫ Lσ ) 2 / 4 + t 2 H and ∆ 3 σ = p ( ǫ Lσ + U L − ǫ Rσ − U R ) 2 / 4 + t 2 H , respectively . Energies of tw o-particle singl et states are solutions of a cubic equation, ∆ ǫ 3 +( U L + U R − U I − J )(∆ ǫ 2 − [ E L − E R ] 2 ) − ([ E L − E R ] 2 + 16 t 2 H )∆ ǫ = 0 and t h eir eigenstates are det ermined by u k = (∆ ǫ k + E L − E R ) w k / √ 8 t H and v k = (∆ ǫ k + U L + U R − U I − J ) u k / √ 8 t H − w k under a normalization of u 2 k + v 2 k + w 2 k = 1. W e also abbreviate E α = ǫ α ↑ + ǫ α ↓ + U α . 0-particle state | m = 0 i ≡| e i ǫ 0 = 0 1-particle state | m = 1 , 2 i = (cos φ σ c † Lσ − sin φ σ c † Rσ ) | e i ǫ 1 , 2 = ( ǫ Lσ + ǫ Rσ ) / 2 − ∆ 1 σ | m = 3 , 4 i = (sin φ σ c † Lσ + cos φ σ c † Rσ ) | e i ǫ 3 , 4 = ( ǫ Lσ + ǫ Rσ ) / 2 + ∆ 1 σ with tan 2 φ σ = 2 t H / ( ǫ Rσ − ǫ Lσ ) 2-particle state | m = 5 i = c † L ↑ c † R ↑ | e i ǫ 5 = ǫ L ↑ + ǫ R ↑ + ( U I − J ) / 2 | m = 6 i = c † L ↓ c † R ↓ | e i ǫ 6 = ǫ L ↓ + ǫ R ↓ + ( U I − J ) / 2 | m = 7 i = 1 √ 2 ( c † L ↑ c † R ↓ + c † L ↓ c † R ↑ ) | e i ǫ 7 = ( ǫ 5 + ǫ 6 ) / 2 | m = k i = ( u k ˆ S g + v k c † L ↑ c † L ↓ + w k c † R ↑ c † R ↓ ) | e i ǫ k = ( E L + E R ) / 2 + ∆ ǫ k / 2, where ˆ S g = 1 √ 2 ( c † L ↑ c † R ↓ − c † L ↓ c † R ↑ ) , k = 8 , 9 , 10 3-particle state | m = 11 , 12 i = (cos θ σ c Lσ − sin θ σ c Rσ ) | f i ǫ 11 , 12 = ( ǫ 15 + ǫ L ¯ σ + ǫ R ¯ σ ) / 2 − ∆ 3 σ | m = 13 , 14 i = (sin θ σ c Lσ + cos θ σ c Rσ ) | f i ǫ 13 , 14 = ( ǫ 15 + ǫ L ¯ σ + ǫ R ¯ σ ) / 2 + ∆ 3 σ with tan 2 θ σ = 2 t H / ( ǫ Rσ + U R − ǫ Lσ − U L ) 4-particle state | m = 15 i = c † L ↑ c † L ↓ c † R ↑ c † R ↓ | e i ≡| f i ǫ 15 = E L + E R + 2 U I − J B. Generating functional F o r the ease in ev aluating the exp ecta tion v alue of an y o p er ator O , we introduce a Lag range multiplier λ asso cia ted with the aux iliary pa rticle n umber Q as, H → H + λQ − X α µ α n leads α , (7) with n leads α = P kσ a † kασ a kασ . Then, the system bec omes the gr and canonical ensemble with resp ect to the auxiliary particle n um ber Q , i.e. , Q is now uncons trained. With the gr and canonica l ensemble, we define a generating functiona l W = − ln Z a s a n extensio n of the Gibbs free e ne r gy . Here, the gener alized par tition function Z , in terms o f the coherent path integral r epresentation is given by ,[22, 23] Z = I D [ c ∗ ασ , c ασ , a ∗ kασ , a kασ ] e − S/i ¯ h , (8) with the ac tion r epresented on a clos ed time contour as, S = I dτ h X m d ∗ m ( τ )  i ¯ h∂ τ − ǫ m − λ  d m ( τ ) + X kασ a ∗ kασ ( τ )  i ¯ h∂ τ − ǫ kασ + µ α  a kασ ( τ ) − X kασ  T α kσ a ∗ kασ ( τ ) c ασ ( τ ) + T α ∗ kσ c ∗ ασ ( τ ) a kασ ( τ )  i . ( 9) W e note that the F ermi (Bose) particle o p e rators are no w replaced b y the cor resp onding Grassman (complex) v ar ia bles { a ∗ ασ , a ασ , d ∗ m , d m } . F r om the unconstra ined genera ting functional W , the ex pe c tation v a lue in the Q = 1 ensemble can b e calc ulated easily by noting that the o p erator Q commutes with the tota l Hamiltonian, and Q is thus a go o d quantum num b er. This fact enables us to expa nd the pa rtition function in pow ers of ζ = e − λβ . The partition function belong ing to the Q = 1 subs pace can b e obtained by differentiating with res pe c t to ζ , i.e. , Z Q =1 = lim λ →∞ ∂ ∂ ζ Z . Based on this relation, we ca n ev aluate the expectation v alue of O by taking a functiona l deriv ative of Z with resp ect to its conjugate v ariable η as , hOi C = − 1 Z Q =1 δ δ η Z Q =1 = lim λ →∞  hOi GC + ( ∂ /∂ ζ ) hOi GC e β λ h Q i GC  . (10) where hO i C denotes the av erage ov er the canonical ensem ble, i.e., over the subspa c e Q = 1 while hOi GC ≡ δ W /δ η is the average over the grand canonical ensemble. When the ca nonical ex p ecta tion v alue of the op er ator O ha s a zer o 4 exp ectation v alue in the Q = 0 subs pa ce, the ab ove relation is further simplified to, hOi C = lim λ →∞ hOi GC h Q i GC (11) As seen in the following section, since the e xp ectation v alues of interest hav e a zero exp ectation v alue in the Q = 0 subspace w e her eafter focus on the a verage o ver the gr and canonical ensemble based on Eq. (11). The partition function in the grand ca nonical ensem ble is calculated f ollowing the standard series expa nsion pro cedure.[23] Firstly , we integrate the action o ver Grassman v a riables { a ∗ kασ , a kασ } a nd o bta in S = I dτ X m d ∗ m ( τ )  i ¯ h∂ τ − ǫ m − λ  d m ( τ ) − I dτ dτ ′ | T α kσ | 2 c ∗ ασ ( τ ) g kασ ( τ , τ ′ ) c ασ ( τ ′ ) . (12) Then w e expr ess the v aria bles { c ∗ ασ , c ασ } in Eq . (5) in terms of their functional deriv atives as, c ∗ ασ ( τ ) = X mm ′ ξ ασ mm ′ ς m δ δ η m ( τ ) δ δ η ∗ m ′ ( τ ) c ασ ( τ ) = X mm ′ ξ ασ ∗ mm ′ ς m ′ δ δ η m ′ ( τ ) δ δ η ∗ m ( τ ) (13) where η m and η ∗ m are their corr esp onding conjugate v ar ia bles. In that way the partition function is further rewritten as, Z = Z 0 exp ( 1 i ¯ h X ασ I dτ dτ ′ c ∗ ασ ( τ ) g ασ ( τ , τ ′ ) c ασ ( τ ′ ) ) e i ¯ h P m H dτ dτ ′ η ∗ m g m ( τ ,τ ′ ) η m ( τ ′ ) (14) where Z 0 is the unp er turb ed par tition function and g ασ ( τ , τ ′ ) ≡ P k | T α kσ | 2 g kασ ( τ , τ ′ ). H ere, the low e r c a se g ’s are the unpe rturb ed Gr een’s functions o f the lead-electr o ns and auxiliar y pa rticles; for instance, re tarded, Keldysh, a nd adv anced compo nents of the Gr een’s functions o f the m -th auxiliary par ticle a re given b y g R m ( t, t ′ ) = 1 i ¯ h θ ( t − t ′ ) h [ d m ( t ) , d † m ( t ′ )] − ς m i 0 GC = 1 i ¯ h θ ( t − t ′ ) e ( ǫ m + λ )( t − t ′ ) /i ¯ h g K m ( t, t ′ ) = 1 i ¯ h h [ d m ( t ) , d † m ( t ′ )] ς m i 0 GC = 1 i ¯ h  tanh β ( ǫ m + λ ) 2  − ς m e ( ǫ m + λ )( t − t ′ ) /i ¯ h g A m ( t, t ′ ) = g R ∗ m ( t ′ , t ) (15) resp ectively , where the sup er script ’0’ denotes the average in the case of the dots being deco upled from the leads, and ς m is − 1(+1 ) if the particle m is a fermion (b oson). T o simplify the expressio ns, it is sometimes convenien t to use gr eater, g > = ( g K + g R − g A ) / 2, le s ser, g < = ( g K − g R + g A ) / 2, a nd cor related, g C = g R − g A , Green’s functions int erchangeably . Next, we expand the expo nential function in Eq. (14) in p ow er series a nd obtain the partition function by collecting all the connected diag rams; Z = Z 0 e −{W (1) + W (2) + W (3) + ···} . (16) Here, W ( n ) is the collectio n of the | T α kσ | n -order diagrams. F or instance, W (1) and W (2) lo ok lik e, W (1) = − i ¯ h X m ς m -1.2 -0.7 -0.2 0.3 0.8 -1.2 -0.7 -0.2 0.3 0.8 m , W (2) = − i ¯ h 2 X m ς m   -1.2 -0.7 -0.2 0.3 0.8 -1.2 -0.7 -0.2 0.3 0.8 m + 2 -1.2 -0.7 -0.2 0.3 0.8 -1.2 -0.7 -0.2 0.3 0.8 m + 2 -1.2 -0.7 -0.2 0.3 0.8 -1.2 -0.7 -0.2 0.3 0.8 m + -1.2 -0.7 -0.2 0.3 0.8 -1.2 -0.7 -0.2 0.3 0.8 m   , (17) where the solid (dotted) lines denote the unp erturb ed Green’s functions g m ( g ασ ). Her e , large dots indicate the times at which the tunneling even ts o cc ur and the o verlap matrix is assumed to be multiplied as , -1.2 -0.7 -0.2 0.3 0.8 -1.2 -0.7 -0.2 0.3 0.8 m m’ ασ × ξ ασ ∗ m ′ m , -1.2 -0.7 -0.2 0.3 0.8 -1.2 -0.7 -0.2 0.3 0.8 m m’ ασ × ξ ασ mm ′ . (18) 5 Finally , the calculatio n of the generating functional W = − ln Z is done in a straightforward fashion via the Luttinger- W a r d functiona l Φ and rep eated terms;[24] W = − ln Z (0) + W (1) + W (2) + · · · = Φ + X p ς p T r  ln G − 1 p + Σ p G p  , p = ασ and m. (19) Here, the Luttinger-W ard functional Φ is the sum o f all the closed skeleton diag rams with a non-interacting Green’s functions ( g p ) replaced by the full Green’s functions ( G p ). Up to now, the generating functional of the co upled-dot system has b een der ived without a ny approximations, and th us the as so ciated Green’s functions g ive the exa ct expressio ns for the physical quan tities as sho wn in the Appendices B and C. In the next section w e describ e the approximation of the Luttinger-W ard functional Φ, and present the expressions for the physical qua nt ities in the static case. C. Non-crossing approximation and pro jection to Q = 1 Hereafter w e employ the non-cros s ing approximation, that is, w e confine our attention to the first skeleton diag ram originated from W (1) , and appr oximate the Luttinger-W ard functional by , Φ = − i ¯ h X m ς m -1.2 -0.7 -0.2 0.3 0.8 -1.2 -0.7 -0.2 0.3 0.8 m = − i ¯ h X mm ′ ασ ς m | ξ ασ m,m ′ | 2 I dτ dτ ′ G m ( τ ′ , τ ) G ασ ( τ , τ ′ ) G m ′ ( τ , τ ′ ) (20) where the thick lines re pr esent the full Green’s functions of the par ticles instead o f unp erturb ed ones (thin lines) in Eq. (17). Then, since the generating function W is statio nary with resp ect to G p , namely δ W /δ G p = 0 , the self-energies ca n be obta ined fro m, Σ p ( τ , τ ′ ) = − ς p δ Φ δ G p ( τ ′ , τ ) . (21) Using the NCA functional Eq. (20), this gives, Σ ασ ( τ , τ ′ ) = − i ¯ h X mm ′ ς m | ξ ασ m,m ′ | 2 G m ( τ , τ ′ ) G m ′ ( τ ′ , τ ) , Σ m ( τ , τ ′ ) = i ¯ h X m ′ ασ h | ξ ασ m,m ′ | 2 G ασ ( τ , τ ′ ) − | ξ ασ m ′ ,m | 2 G ασ ( τ ′ , τ ) i G m ′ ( τ , τ ′ ) , (22) for the self-energies o f electrons in the leads and auxiliary par ticles. As functions o f real time arg umen ts, the a b ov e expressions can b e rewr itten as, Σ R ασ ( t, t ′ ; λ ) = − i ¯ h X mm ′ ς m | ξ ασ m,m ′ | 2  G R m ( t, t ′ ; λ ) G < m ′ ( t ′ , t ; λ ) + G < m ( t, t ′ ; λ ) G A m ′ ( t ′ , t ; λ )  , Σ K ασ ( t, t ′ ; λ ) = − i ¯ h X mm ′ ς m | ξ ασ m,m ′ | 2  G < m ( t, t ′ ; λ ) G C m ′ ( t ′ , t ; λ ) + G C m ( t, t ′ ; λ ) G < m ′ ( t ′ , t ; λ )  (23) for the electr ons in the leads, and Σ R m ( t, t ′ ; λ ) = i ¯ h X m ′ ασ  | ξ ασ m,m ′ | 2 G > ασ ( t, t ′ ; λ ) − | ξ ασ m ′ ,m | 2 G < ασ ( t ′ , t ; λ )  G R m ′ ( t, t ′ ; λ ) , Σ < m ( t, t ′ ; λ ) = i ¯ h X m ′ ασ  | ξ ασ m,m ′ | 2 G < ασ ( t, t ′ ; λ ) − | ξ ασ m ′ ,m | 2 G > ασ ( t ′ , t ; λ )  G < m ′ ( t, t ′ ; λ ) . (24) for the aux ilia ry par ticles. 6 Next w e pro ject the self-ener g y to the Q = 1 ensemble (details can b e found elsewher e[17, 1 8, 25]). F o r this we exploit t wo fa cts ab out the lesser a nd greater comp onents o f the Green’s functions for the auxiliary particles. Fir stly , it is impo rtant to note that self-energies of electrons in the leads dep end o n lesser comp onents of Gr een function G < m . This makes the pro jection to Q = 1 subspace easy b ecause G < m ( t, t ′ ) = ς m i ¯ h h d m ( t ′ ) † d m ( t ) i GC = O ( e − β λ ) (25) means a zero exp ectatio n v alue in the Q = 0 subspace, and th us one ca n use Eq. (11) in ev a lua ting observ a bles. Secondly , the lesser and gr eater co mpo nent s o f the Gr een’s functions for the auxiliary particles are given by , G <,> m ( t, t ′ ) = Z dt 1 dt 2 G R m ( t, t 1 )Σ <,> m ( t 1 , t 2 ) G A m ( t 2 , t ′ ) (26) without dependence on g <,> m due to the los s of memory .[17, 20] Before the pro jection to the Q = 1 subspace, we eliminate the λ -dependence in g R,A m . Since λ is r elated to g R,A m ( t, t ′ ) only through the factor e λ ( t − t ′ ) /i ¯ h in Eq. (15), the elimination o f λ such as g m ( t, t ′ ; λ ) → g m ( t, t ′ ) ≡ g m ( t, t ′ ; λ = 0) results in the mo dified forms of self-e ne r gies in Eqs . (2 3) and (24) like, Σ p ( t, t ′ ; λ ) → Σ p ( t, t ′ ; λ ) e − λ ( t − t ′ ) /i ¯ h . By us ing Eq. (11) and taking the pro jection of λ → ∞ , the pro jected self-ener gies o f Eq. (23) become, Σ R ασ ( t, t ′ ) = − i ¯ h X mm ′ ς m | ξ ασ m,m ′ | 2  G R m ( t, t ′ ) G < m ′ ( t ′ , t ) + G < m ( t, t ′ ) G A m ′ ( t ′ , t )  , Σ K ασ ( t, t ′ ) = − i ¯ h X mm ′ ς m | ξ ασ m,m ′ | 2  G < m ( t, t ′ ) G C m ′ ( t ′ , t ) + G C m ( t, t ′ ) G < m ′ ( t ′ , t )  (27) where w e use the abbrevia ted notation of Σ R,K ασ ( t, t ′ ) ≡ lim λ →∞ e λβ Σ R,K ασ ( t, t ′ ; λ ) e − λ ( t − t ′ ) /i ¯ h and set a e β λ h Q i GC term aside in Eq. (11) for a while. Her e, G < m ( t, t ′ ) is defined in Eq. (26) with its se lf-e nergy Σ < m ( t, t ′ ) ≡ lim λ →∞ e λβ Σ < m ( t, t ′ ; λ ) e − λ ( t − t ′ ) /i ¯ h . Using Eq. (2 4), the self-energy is giv en by , Σ < m ( t, t ′ ) = i ¯ h X m ′ ασ  | ξ ασ m,m ′ | 2 g < ασ ( t, t ′ ) − | ξ ασ m ′ ,m | 2 g > ασ ( t ′ , t )  G < m ′ ( t, t ′ ) . (28) Whereas, the Dys on eq uation of G R,A m ( t, t ′ ) is G R,A m ( t, t ′ ) = g R,A m ( t, t ′ ) + Z dt 1 dt 2 g R,A m ( t, t 1 )Σ R,A m ( t 1 , t 2 ) G R,A m ( t 2 , t ′ ) with its se lf-energy defined by Σ R,A m ( t, t ′ ) ≡ lim λ →∞ Σ R,A m ( t, t ′ ; λ ) e − λ ( t − t ′ ) /i ¯ h ; from Eq . (2 4), Σ R m ( t, t ′ ) = i ¯ h X m ′ ασ  | ξ ασ m,m ′ | 2 g > ασ ( t, t ′ ) − | ξ ασ m ′ ,m | 2 g < ασ ( t ′ , t )  G R m ′ ( t, t ′ ) . (29) During the pro jection, we employ the rela tion, G <,> ασ ( t, t ′ ; λ ) = g <,> ασ ( t, t ′ ) + Z dt 1 dt 2 G R ασ ( t, t 1 )Σ <,> ασ ( t 1 , t 2 ; λ ) G A ασ ( t 2 , t ′ ) (30) and neg lect the s econd term, due to its O ( e − λβ ) dependence as seen from Eq. (24). On the other hand, the expec ta tion v a lue of the op erato r Q is given by lim λ →∞ e λβ h Q i GC = lim λ →∞ e λβ i ¯ h X m ς m G < m ( t, t ; λ ) = i ¯ h X m ς m G < m ( t, t ) (31) where the seco nd step can b e derived in a similar wa y to that of App e ndix B. Throughout this work, w e keep lim λ →∞ e λβ h Q i GC to be unity via the normalization a nd cons equently Eqs. (2 7)-(28) are also the av er aged v alues in the canonical ensemble. Eqs. (27)-(28) are the main results of this work, which can be applied to a double-do t system at arbitrary temper - ature, Coulomb in teraction, source-dra in and ga te voltage co nfigurations, including the time-dep endent pr oblems. 7 D. Physical q uant ities in static cases Since in the s tatic ca se, the Green’s functions dep end only on the time interv al, it b ecomes convenien t to us e the F o urier transform, G ( t, t ′ ) = 1 2 π ¯ h Z ∞ −∞ dE e E ( t − t ′ ) /i ¯ h G ( E ) . (32) By using the cut-off, ρ ασ c ( E ) = 2 π P k | T α kσ | 2 δ ( E − ǫ kσ ), the unpertur be d Green’s function o f e le ctrons in the lead α , is then given in the energy space a s, g <,> ασ ( E ) = ± i ¯ hρ ασ c ( E ) f ( ± ( E − µ α )) , (33) where f ( E ) = 1 / (1 + e β E ) is the F ermi-Dirac distribution function. By defining the s pe c tral function A m such that G < m ( t, t ′ ) = − 2 π iς m A m ( t, t ′ ), the self-e nergies o f Eq s. (27)-(28) ar e r ewritten as, Σ R ασ ( E ) = X mm ′ Z ∞ −∞ dE ′  | ξ ασ m ′ ,m | 2 G R m ′ ( E + E ′ ) − | ξ ασ m,m ′ | 2 G A m ′ ( E ′ − E )  A m ( E ′ ) , Σ K ασ ( E ) = X mm ′ Z ∞ −∞ dE ′  | ξ ασ m ′ ,m | 2 G C m ′ ( E + E ′ ) − | ξ ασ m,m ′ | 2 G C m ′ ( E ′ − E )  A m ( E ′ ) Σ R m ( E ) = i 2 π X ασm ′ Z ∞ −∞ dE ′  | ξ ασ m,m ′ | 2 g > ασ ( E − E ′ ) − | ξ ασ m ′ ,m | 2 g < ασ ( E ′ − E )  G R m ′ ( E ′ ) . (34) Here, the sp ectra l function A m ( E ) is deter mined fro m A m ( E ) = i 2 π | G R m ( E ) | 2 X ασm ′ Z ∞ −∞ dE ′  | ξ ασ m ′ ,m | 2 g > ασ ( E ′ − E ) − | ξ ασ m,m ′ | 2 g < ασ ( E − E ′ )  A m ′ ( E ′ ) with the nor malization c o ndition of P m R dE A m ( E ) = 1 fro m E q. (31). On the other hand, the expectatio n v alues of the ph ys ical q ua ntit ies in the Q = 1 ensemble can b e obtained b y combining the re s ults of the Appendix B with Eq. (3 4). W e summarize the results, in the energy space; for the current in the lead α , I α = q 2 π ¯ h ℜ X σ Z ∞ −∞ dE  g > ασ ( E )Σ < ασ ( E ) − g < ασ ( E )Σ > ασ ( E )  , (35) for the density of states, D OS ασ ( E ) = − 1 π Im G R ασ ( E ) = − 1 π ImΣ R ασ ( E ) , (36) for the o ccupa tion num ber , h n ασ i C = X m  ∂ ǫ m ∂ ǫ ασ  Z ∞ −∞ dE A m ( E ) , (37) and for the spin-spin corre lations, S 2 = h ~ S L · ~ S R i C = − X m ∂ ǫ m ( U I → U I + J / 2) ∂ J Z ∞ −∞ dE A m ( E ) . (38) II I. RESUL TS AND D ISCUSSION In this s ection, we illustrate the num erical solutions of Eq. (34), and the resulting ph ysical quantit ies of E qs. (35) - (3 8 ), as well as the accuracy of the present theory . 8 F o r the cut-off function we choose a L o rentzian mo del ρ ασ c ( E ) ρ ασ c ( E ) = Γ ασ W 2 ( E − µ α ) 2 + W 2 ) (39) with W b eing the half width of the conduction ba nd. In solving the Dyson’s equations with self-energies given in Eq. (34), we use the adaptive mes h scheme where mor e mesh p oints are inserted into a hig h weigh ted r egion for every interaction. The itera tion is rep eated until the following sum rules co nv erg e within 0.01%; − 1 π Z ∞ −∞ Im G R m ( E ) dE = 1 , X m Z ∞ −∞ A m ( E ) dE = 1 . T o achiev e this numerical accuracy , w e use abo ut 1000 mesh p oints for each Green’s function of a n auxiliary particle . F o r simplicity we consider the sy mmetr ic case, U L = U R = U , ǫ Lσ = ǫ Rσ = ǫ d , U I = 0, with Γ ασ = Γ, and J = 0, and all the energies are measured in units of Γ (in a n exper iment Γ is t ypically of the or de r of µ eV to meV). W e present results for tw o kinds of sy s tems; one is a single quan tum do t (that is, we take t H → ∞ ) and the other is a double q ua nt um dot (finite t H ). Although the s ing le qua ntu m dot case has b een extensively studied, we revisit the problem to s how that o ur formulation indeed encompas ses the previous results. A. Single quantum dot W e first consider a single quant um dot and ex amine the correlated q ua ntu m transp ort through it. T o do this we write ǫ d → ǫ d + t H , and take t H → ∞ . Then, there a re four lo w-lying states relev a nt to transp ort: | 0 i = | e i , | 1 , 2 i = 1 √ 2 ( c † Lσ − c † Rσ ) | e i , and | 8 i = 1 2 ( c † L ↑ − c † R ↑ )( c † L ↓ − c † R ↓ ) | e i , while their energies a r e given by ǫ 0 = 0 , ǫ 1 , 2 = ǫ d , and ǫ 8 = 2 ǫ d + 1 4 ( U L + U R + U I + J ), r esp ectively . In the limit U → ∞ , the s tate | 8 i can be discar ded further , a nd the present formalis m rec overs the r e sults o f Ref. [17]. F or this case, res ults of the system describ ed by t y pical parameters are shown in Fig. 1. In Fig . 1-(a), we plot the e q uilibrium densities of states fo r sev eral temp era tur es, where the broad p ea ks are caus ed by the usual tra nsitions betw een le vels (in this case | 0 i and | 1 , 2 i ), while the sharp ones (lo c ated a t E = µ = 2 . 0) a re the Kondo p eaks. The later ones incr ease a s temp era ture is lowered,[17 ] with satura tio n well b elow the Kondo temper ature [19] T K = min ( U √ I 2 π , r W ¯ h Γ 2 ) e − π /I , (40) where I = ¯ h Γ  1 | ǫ d − µ | + 1 ǫ d + U − µ |  . F r om this re la tion, T K = 2 . 8 × 1 0 − 3 is estimated in the case of Fig. 1 while the calculated Kondo temp era tur e, equal to its half width at half maximum, is 3 . 2 × 1 0 − 2 . The ov er- estimation of the Kondo tempera ture is a known consequence the NCA, as w ell as the Kondo p eak height.[17] In the figure Fig. 1-(b), we show the v aria tion of the dot o ccupation, n B = h n m =1 i + h n m =2 i as a function of the dot ener gy level, a nd in Fig. 1-(c), the relation b etw een the ele c tronic occupatio n n B and the heigh t of the density of state at E = µ . Actually , n B and the densit y of state are r elated throug h the F riedel s um rule: D OS ( E = µ 0 ) = 1 π Γ sin 2 ( π n B / 2) . (41) F o r the v alidity o f our calculatio ns, w e also plot previous res ults of the NCA and numerical r enormaliza tion g r oup (NR G) metho d from Ref. [18]. In the case of the occupation, w e find that our results are in go o d agr eement with the previous results, with a minor deviation resulting fro m the use of a different cut-off function. On the other hand, in the compariso n of the F riedel sum rule, a larg e deviation of our res ults are found fro m those of the exa c t result and NR G. As seen in res ults of Ref. [18], the pr evious NCA calculation also show nea r ly the same deviatio n. This fact leads us to the ov er -estimated Kondo pe ak with the NCA in the wider ra nge of the occupa tion. 9 -5 -3 -1 1 3 E 0.0 0.2 0.4 DOS (a) 1.99 2.00 2.01 0.0 0.2 0.4 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 ε d −µ 0 0.0 0.2 0.4 0.6 0.8 1.0 n B (b) 0.0 0.2 0.4 0.6 0.8 1.0 n B 0.0 0.1 0.2 0.3 0.4 DOS( µ 0 ) (c) Figure 1: In (a), the equilibrium densities of states are plotted fo r a single q uantum dot with an infinite rep ulsive p otentia l. The inset sho ws the blo w-u p around Kondo p eaks for temperatures T = 1 . 0 × 10 − 4 (solid), 5 . 0 × 10 − 4 (dotted), 1 . 0 × 10 − 3 (dashed), and 5 . 0 × 10 − 3 (dot-dashed), resp ectively . W e use the parameters of ǫ d = − 0 . 5, µ 0 = 2 . 0, ∆ µ = 0, and W = 100 . 0 whic h gives T K = 2 . 8 × 10 − 3 in Eq. (40) . In ( b ), w e compare the electronic occupation n B as a function of energy level with the previous results. Crosses represen t results from NRG, b oxes from the NCA of R ef. [18] at T = 0 and circles from our approac h at T = 1 . 0 × 10 − 4 , resp ectively . Since our calculation co de is n ot a v ailable at T = 0, w e c hoose sufficiently lo w temp erature for the comparison. In (c), th e sum rule is examined together with the exact result (solid line) from Eq. ( 41). When the Coulomb interaction is finite, all of the four states | 0 i , | 1 i , | 2 i , and | 8 i take part in electron transp ort. F o r the sa me quantum dot of Fig. 1, but with a finite p otential U , w e sho w the DOS a s a function of E in Fig. (2)-(a). As the Co ulo mb potential decr eases, the Kondo p eak is lowered be c ause the electr o n correlation is unimpo rtant. This is a lso predicted by the F riedel sum rule of Eq. (41). Since the tw o-particle sta te | 8 i beco mes ene r getically fav orable with smaller p o tent ial U , the electron o ccupation increases up to t wo. In Fig. 2-(b), we plot the heigh t of the Kondo pea k a s a function of o ccupation n B . It is found that the present result (circles ) exhibits the same decreasing behavior with larg er o ccupa tion as the F riedel s um rule, how ever still shows the ov er -estimation o f the K ondo p eak. In Fig. 2-(c), the s pin c orrela tion is shown as a function of the Coulo mb p otential U and is compared with the o cc upation n B . Since the spin co rrelatio n originates fr om only a tw o-particle state, re sults in the figur es a re prop ortional to the o ccupancy of | 8 i auxilia ry particle. Thus, one can se e that as the o ccupa tion n B approaches tw o, the spin corre la tion bec omes S 2 = h 8 | ~ S L · ~ S R | 8 i = − 3 / 8. Fig. 3 shows the linear co nductance (a), and the ele ctronic o ccupation n B (b) a s a function o f the sing le - particle energy for a finite Coulom b p otential U = 10 . 0. F o r this calculation, w e apply a small voltage b etw een the left and right leads of ∆ µ = 0 . 0 1 and the conducta nce is calcula ted a s the cur rent at a lead div ided by ∆ µ . As the temp era ture is low ered, the co nductance increa ses and approa ches 2 e 2 /h , whic h is in accorda nce with the exp erimental results reflecting the Kondo effect.[3] On the other hand, the electro nic o ccupation n B shows w ea k temp erature dependence, as shown in Fig. 3-(b). The conductance maximum are appr oximately at n B = 0 . 5 a nd 1 . 5, which coincides with the condition of mo s t pro bable sequen tial tunneling: µ = ǫ 1 , 2 and µ + ǫ 1 , 2 = ǫ 8 . B. Coupled-quantum dots When the coupling strength t H betw een the dots is finite, the system now r e presents double quantum dots a nd a ll 16 many-b o dy molecular states tak e part in the transp ort. As the first exa mple, we consider the cas e of U → ∞ and 10 -6 -2 2 6 E 0.0 0.5 1.0 1.5 DOS U=8.0 4.0 2.0 1.0 0.5 0.0 (a) 1.0 1.2 1.4 1.6 1.8 2.0 n B 0.0 0.1 0.2 0.3 0.4 DOS( µ 0 ) (b) 0 2 4 6 8 U 1.0 1.5 2.0 0.0 -0.2 -0.4 n B S 2 (c) Figure 2: In (a), we show the equ ilibrium density of states for the single qu antum dot with the same p arameters as Fig. 1, but with fi n ite potential U . In (b), the density of states at the chemical p otential is plotted as a function of the occup ation n B , where the solid line is the F riedel sum rule of Eq. (41) and circles is th e presen t result, resp ectively . In (c), we p lot the v ariation of th e occupation n B and spin correlation S 2 with resp ect to the p otential U . Dotted lines are the guide for ey es. U I = J = 0. Then, low-lying states are ǫ 0 = 0 , ǫ 1 , 2 = ǫ d − t H , ǫ 3 , 4 = ǫ d + t H , and ǫ 5 , 6 , 7 , 8 = 2 ǫ d from T able I. Due to the larg e Coulomb p o tential U , o ne ca n see that the do uble o cc upa tion on ea ch qua nt um dot is pro hibited. And one ex p ects that sequential tunneling occ ur s dominantly for tw o c o nditions o f ǫ d = t H and ǫ d = − t H . The former corres p o nds to the tra nsition betw e en | 0 i and | 1 , 2 i , and the latter is that b etw een | 1 , 2 i and | 5 , 6 , 7 , 8 i . In Fig. 4, we examine the conductance as a function of the chemical p otential differ e nce in the vicinit y of the latter case. F or a given ǫ d = − 2 . 5, w e compare calcula ted conductance for t H = 2 . 0, 2 . 6, and 3 . 2. Among three cases, ov er all conductance for t H = 2 . 0 shows the largest v a lue. It is interesting b ecause the la rgest one will b e the case t H = 2 . 6 acco rding to the sequential tunneling conditio n of ǫ d = − t H . W e attribute this to the level renormaliza tion owing to the electron correla tion. On the o ther hand, sharp peaks a re found a r ound ∆ µ = 0, whose height increa ses as tempera tures are lo wered. The p eaks are found to r esult from the K ondo effect a s inferred from the density of states in Fig. 4-(b). Actually , the similar calc ula tion is alrea dy p er fo rmed in Ref. [12], wher e double p eaks o f the conductance ar ound ∆ µ = 0 differently from the pr e s ent result are observed. W e a ttribute the discrepancy betw een a single peak and a double p ea k predicted in each work to the difference in the formulation o f the problem. While in Ref. [12], the lo calized basis such as c † ασ | e i is used, we use the dia gonalizing basis shown in T a ble I. Strictly sp eaking , the present work trea ts the double-dot system as a single-qua ntum dot with m ulti-le vel molecula r states, which leads to the mo dified coupling strengths b etw een the dots and the leads, weigh ted b y ξ ασ mm ′ in Eq. (5). Therefor e, a lthough bo th approaches ado pt the NCA, the details of the F ey nma n diagrams are differen t, and we expect that the results of b o th approaches w ould con verge by including more c rossing diag rams. Finally , w e consider the coupled quantum do ts with finite Coulomb p otential. When the Coulomb p otential becomes comparable to the dot-do t in ter action t H , their comp e tition giv es rise to the rich electronic structure and all the ener g y levels may b e relev ant to the tra nsp ort. I n Fig. (5)-(a), w e show the v a r iation of energy levels as a function of the Coulomb p o tential in the ca se o f iso lated co upled do ts . As the Coulom b p otential decr eases fr o m infinity , it is found that more lev els fall in to the r ange of r elev ant energ y . In o ther w ords, this means that v arious transitions b etw een states beco me a v ailable and a r e r esp onsible for mo r e p eak s in the density o f s ta tes as shown in Fig. (5)-(b). Due 11 -10 -6 -2 2 6 0.0 0.2 0.4 0.6 0.8 1.0 conductance(2e 2 /h) (a) -10 -6 -2 2 6 ε d - µ 0 0.0 0.5 1.0 1.5 2.0 n B (b) Figure 3: W e plot th e linear conductance (in unit of 2 e 2 /h ) as a function of single-particle energy in ( a ) and corresponding electronic o ccu pation (b) , for a fin ite Coulom b p otentia l U = 10 . 0 at temp eratures T = 0 . 00 3(solid), 0 . 03(dotted), and 0 . 3(dashed), resp ectively . W e u se W = 50 . 0 and µ 0 = 4 . 0. The linear conductance is calculated with a finite p otential difference ∆ µ = 0 . 01 through I α / ∆ µ . -0.04 -0.02 0.00 0.02 0.04 ∆µ 0.0 0.2 0.4 0.6 0.8 1.0 conductance(2e 2 /h) t H =2.0 2.6 3.2 (a) -10 -6 -2 2 6 E 0.0 0.1 0.2 0.3 0.4 DOS T=0.001 (b) -0.2 0.0 0.2 0.0 0.2 0.4 Figure 4: In (a), t he cond u ctance (in unit of 2 e 2 /h ) of the coupled-q uantum dot with U = ∞ and U I = 0 are show n as a function of the chemica l potential difference for temperatures T = 0 . 001 (thic k), 0 . 0 05 (medium), and 0 . 01 (thin), respectively . In (b), w e plot the d ensit y of states at ∆ µ = 0 . 0 and T = 0 . 0 01. The v alue of th e in terdot in teraction parameter for each of the curves is t H = 2 . 0 (solid), 2 . 6 (dotted), and 3 . 2 (dashed) in both panels and W = 10 . 0, ǫ d = − 2 . 5, and µ 0 = 0. to the detailed change of energy levels, the conductance is als o found to b e lar gely modified. In Fig. (5)-(c), the conductance are shown for three different Co ulomb po ten tials. Compared to that of the infinite Co ulomb po ten tial case in Fig. (4)-(c), calculated results ar e lar gely s uppressed. This is b eca us e the transition energies determined from the comp etition of v ario us in ter actions are to o large for electrons to tunnel through dots, which is similar to the Coulomb block ade effect for la rge p otential U . Nev er theless, one can see shar p pe aks in the calculated conductance at ∆ µ = 0. These p ea k s result from the Ko ndo effect as in the case of the infinite Coulomb p otential, meaning that correla ted trans po rt still occ urs even in small Coulomb p otential. W e find that the height s of the conductance at ∆ µ = 0 a re m uch la rger than thos e of the master equatio n approa ch, ho wev er , smaller than those of NRG (not sho wn 12 0.0 10.0 20.0 30.0 U -15 -5 5 15 ε m (a) t 21 t 21’ t 23 -10 -5 0 5 10 E 0.0 0.1 0.2 0.3 0.4 DOS(E) (b) t 23 t 21 t 21’ t 23 -0.04 -0.02 0.00 0.02 0.04 ∆µ 0.02 0.04 0.06 0.08 conductance(2e 2 /h) (c) U=10.0 3.0 5.0 Figure 5: F or the coupled qu antum dots of Fig. 4, the v ariation of energy levels is plotted in (a) as a fun ction of Coulomb p otentia l U where one-, tw o-, three-particle states are represen ted by solid, dotted, and dashed lines, respectively . The arrows indicate possible transitions into one- particle ( t 21 , t 21 ′ ) and th ree-particle ( t 23 ) state from the ground (tw o-particle single) state. In ( b ), w e sho w the density of state for finite Coulomb p otentials U = 3 . 0 ( solid), 5 . 0 (dotted), and 10 . 0 (dashed), respectively , at temperature T = 0 . 003 and displa y the transitions corresp onding to each p eak. In (c), the conductance (in unit of 2 e 2 /h ) is sho wn as a function of the c h emical p otential difference. here).[8, 11] This means that our a ppr oach accounts for co rrelated behavior of electro ns pa r tially . IV. SUMMAR Y In summary , we form ulate the electro n transp ort through tw o laterally coupled quantum dots by ex tending the auxiliary op erator metho d to a mult i-level case, and deriv e the non-equilibrium Green’s function in a conse rving w ay . By using the g enerating functional technique, we present exact expressions for the current through the system, as well as the de ns ities of states, occupa ncies, a nd spin correlations of the dots. T o obtain the Luttinger-W ar d functional, we include the first-or de r diagr am (non-crossing approximation). F or the v alidity o f our results, we examine v ario us situations and compa r e c a lculated results with those of pre vious NCA and exact NR G appr oaches. W e find that our formulation encompasses the previous NCA results success fully , how ever gives the devia ted b ehavior from the NRG metho d. This means that the present metho d accounts for the correla ted behavior partially and the vertex co rrection is needed for more accur ate des cription of transp o r t. Nevertheless, since the present theor y cop e with all the ra nges of the Coulomb energy and occupancie s as well as time-dependent voltages, it can b e applied to reveal transpo rt prop erties of v a rious double-dot problems. 13 Ackno wle dgments Numerical calculations ar e per formed on the sup erco mputer, BlueF ern, at the University of Canterbury , New Zealand. [1] W. G. v an der Wiel, S. De F ranceschi, J. M. Elzerma n, T. F ujisaw a, S. T arucha and L. P . Kouw enhov en, Rev. Mo d. Ph ys. 75 , 1 ( 2003); R. Hanson, L. P . Kouw enhov en, J. R. Petta, S. T aruc ha, L . M. K. V andersyp en, Rev. Mod. Phys. 79 , 1217 (2007). [2] T. K. Ng and P . A. Lee, Ph ys. Rev. Lett. 61 , 1768 (1988); L. I. Glazma n and M. E. Raikh, JETP Lett. 47 452 (1988). [3] D. Goldhab er-Gorden, H. Shtrikman, D. Mahalu, D. A busch-Magder, U . Meira v , and M. A. Kastner, Nature 391 , 156 (1998). [4] J. Schmid, J. W eis, K. Eb erl, and K. v. Klitzing, Phys. R ev . Lett. 84 , 5824 (2000); H. Jeong, A. M. Chang, and M. R. Melloch, Science 293 , 2221 (2001); U. Wilhelm, J. Schmid, J. W eis, and K. v . Klitzing, Physica E 14 , 385 (2002); A. W. Holleitner, R. H. Blic k , A. K. H ¨ uttel, K. Eberl, and J. P . Kotthaus, Science 297 , 70 (2002). [5] D. Loss and D. P . DiVincenzo, Phys. R ev . A 57 , 120 (1998). [6] B. A. Jones, C. M. V arma, and J. W. Wilkins, Phys. R ev. Lett. 61 , 125 (1988). [7] A. Georges and Y. Meir, Phys. Rev. Lett. 82 , 3508 (1999). [8] W. Izumida and O. Sak ai, Ph ys. Rev. B 62 , 10260 (2000). [9] R. Aguado and D. C. Langreth, Phys. Rev. Lett. 85 , 1946 (2000). [10] J. Mra v lje, A. Ram ˘ sak, and T. Rejec, Phys. Rev. B 73 , 241305(R) (2006). [11] M. Lee, M. Choi, R. L´ opez, R . Aguado, J. Martinek, and R . ˘ Zitko, cond -mat.mes-hall/0911. 0959. [12] R. Aguado and D. C. Langreth, Phys. Rev. B 67 , 245307 (2003). [13] J. R. Pe tta, A. C. Johnson, J. M. T aylor, E. A . Laird, A. Y acob y , M. D. Luk in, C. M. Marcus, M. P . Hanson, A. C. Gossard, Science 309 , 2180 (2005). [14] F. H . L. Kopp en s, C. Buizert, K. J. Tielro oij, I . T. Vink , K. C. Now ack, T. Meunier, L. P . K ouw enhov en, and L. M. K. V andersyp en, Nature 442 , 766 (2006). [15] A. C. Johnson, J. R. Petta, C. M. Marcus, M. P . Hanson, and A. C. Gossard, Ph ys. Rev. B 72 , 165308 (2005). [16] K. Ono, D. G. Austing, Y. T okura, and S. T aruc har, S cience 297 , 1313 (2002). [17] N. S. Wingreen and Y. Meir, Phys. Rev. B 49 , 11040 (1994). [18] T. A. Costi, J. Kroha, and P . W¨ olfle, Phys. Rev. B 53 , 1850 (1996). [19] K. Haule, S. Kirchner, J. Kroha, and P . W¨ olfle, Phys. Rev. B 64 , 155111 (2001). [20] G. D. Mahan, in Many-p article physics , 2nd ed. (Plenum, New Y ork, 1990). [21] Z. Zout and P . W. Anderson, Ph ys. Rev. B 37 , 627 (1988). [22] J. W. Negele and H. Orland, in Quantum Many-p article Systems , (Add ison-W esley , 1988); Y. Ustu mi, H . Imamura, M. Hay ashi, and H. Ebisaw a, Ph ys. Rev. B 66 , 024513 (2002). [23] J. H. Oh, D. Ahn, and S. W. Hwang, Ph y s. Rev. B 72 , 165348 (2005). [24] J. M. Luttinger and J. C. W ard, Phys. Rev . 118 , 1417 (1960). [25] P . Coleman, Ph ys. Rev. B 29 , 3035 (1984). App endix A: commutation rel ations In this a ppe ndix we show that the auxiliary pa r ticle r epresentation, c † ασ = X mm ′ ξ ασ m,m ′ d † m d m ′ (A1) gives the correct comm utation relations [ c ασ , c † α ′ σ ′ ] + = δ αα ′ δ σσ ′ . Here, a quasi-pa rticle oper ator d m is fermionic (bo sonic) when the n um b er of pa rticles in a state it repr esents is o dd (even), a nd is assumed to sa tisfy the commutation relation [ d m , d † m ′ ] ± = δ mm ′ . It is imp ortant to note tha t the e x pansion co efficient ξ ασ mm ′ = h m | c † ασ | m ′ i is nonze r o only if the num ber of particles in | m i is larger than that in | m ′ i by one. This means tha t Eq. (A1) is the combin ation of fermion a nd b oso n op er ators. Now, w e ca lculate the co mm utation r elation, [ c ασ , c † α ′ σ ′ ] + = X m 1 ,m 2 ,m 3 ,m 4 ξ ασ ∗ m 2 m 1 ξ α ′ σ ′ m 3 m 4 [ d † m 2 d m 1 , d † m 3 d m 4 ] + . (A2) 14 T o calculate the rig ht -handed side, it is convenien t to sepa rate the sums in to the b osonic ( B ) and fermionic ( F ) terms, X m 1 ,m 2 ,m 3 ,m 4 = X m 2 ,m 3 X m 1 ∈ F + X m 1 ∈ B ! X m 4 ∈ F + X m 4 ∈ B ! =   X m 1 ∈ F ,m 2 ∈ B + X m 1 ∈ B ,m 2 ∈ F     X m 4 ∈ F ,m 3 ∈ B + X m 4 ∈ B ,m 3 ∈ F   =   X m 1 ,m 3 ∈ F X m 2 ,m 4 ∈ B + X m 1 ,m 3 ∈ B X m 2 ,m 4 ∈ F   . (A3) Here, in the second line we us e the fact that c ασ and c † ασ are the pro ducts of fermionic and b oso nic op er ators, and in the third line w e dis play the co llection of non-zer o terms. By substituting Eq. (A3) into Eq. (A2), and using [ d m , d † m ′ ] ± = δ mm ′ , w e ar rive at, [ c ασ , c † α ′ σ ′ ] + = X mm ′ m ′′  ξ ασ ∗ m ′′ m ξ α ′ σ ′ m ′′ m ′ + ξ ασ ∗ m ′ m ′′ ξ α ′ σ ′ mm ′′  d † m d m ′ . ( A4) F ur thermore, since the expans io n co efficient s s atisfy the or thogonality relation, X m ′′  ξ ασ ∗ m ′′ m ξ α ′ σ ′ m ′′ m ′ + ξ ασ ∗ m ′ m ′′ ξ α ′ σ ′ mm ′′  = h m | [ c ασ , c † α ′ σ ′ ] + | m ′ i = δ mm ′ δ αα ′ δ σσ ′ , (A5) the comm utation r elation is simplified to [ c ασ , c † α ′ σ ′ ] + = δ αα ′ δ σσ ′ X m d † m d m = δ αα ′ δ σσ ′ Q. (A6) Thu s, in the subspace Q = 1 , the combination of quas i-particle op era tors leads to the corr e ct comm utation relatio n betw een c ασ and c † α ′ σ ′ . App endix B: Ex pressions for physical quantities In this section we derive the exact expressio ns fo r obser v ables in terms of the Gr een’s functions. As examples, we show the pro cedure for ev aluating the current, I ασ ( t ) = q d dt * X k a † kασ ( t ) a kασ ( t ) + GC = iq ¯ h X k D T α ∗ kσ c † ασ ( t ) a kασ ( t ) − T α kσ a † kασ ( t ) c ασ ( t ) E GC , (B1) Green’s functions, G R ασ ( t, t ′ ) = 1 i ¯ h θ ( t − t ′ ) h [ c ασ ( t ) , c † ασ ( t ′ )] + i GC , (B2) and the o ccupatio n num ber o f electrons in each quantu m dot, h n ασ ( t ) i =  c † ασ ( t ) c ασ ( t )  GC . (B3) T o do this, we attac h fictitious field e − ip ασ ( τ ) to T α kσ for curre nt, and add fictitious ener gy h ασ ( t ) to ǫ ασ for the av erage nu mber. Then, fro m the g e nerating functiona l W the ab ov e quantities ca n b e calcula ted as, I ασ ( t ) = − i q δ W δ ∆ p ασ ( t )     p ασ = h ασ =0 , G ασ ( τ , τ ′ ) = δ [ W − W (0) ] δ g ασ ( τ , τ ′ )     p ασ = h ασ =0 , h n σα ( t ) i = − i ¯ h δ W δ ∆ h m ( t )     p ασ = h ασ =0 , (B4) 15 where p ασ ( t ) = ± ∆ p ασ ( t ) / 2 and h ασ ( t ) = ± ∆ h ασ ( t ) / 2 are assumed o n the upper (+) and low er ( − ) branches of the Keldy s h contour.[22, 23] With the fictitious fields, the ev aluation of the g enerating functional is stra ightforw a rd bec ause bare Gr een’s function is simply changed as, g ασ ( τ , τ ′ ) → e ip ασ ( τ ) g ασ ( τ , τ ′ ) e − ip ασ ( τ ′ ) g − 1 m ( τ , τ ′ ) → ( i ¯ h∂ τ − ǫ m [ ǫ ασ + h ασ ( τ )] − λ ) δ ( τ − τ ′ ) . (B5) In order to ev aluate the functional de r iv atives, we expand the gener a ting functional in ser ies, W = P ∞ n =0 W ( n ) where W (0) = X p ς p T r ln[ g − 1 p /i ¯ h ] , W ( n ) = − X p ς p n I g p ( τ , τ ′ ) ˜ Σ ( n ) p ( τ ′ , τ ) dτ dτ ′ . (B6) Here, g p ( τ , τ ′ ) con tains the fictitious fields, and ˜ Σ ( n ) p ( τ ′ , τ ) represent all the proper and impro pe r n -th order self- energies. By p er forming the functiona l deriv atives we obta in, δ W δ ∆ h ασ ( t ) = − X m ς m ∂ ǫ m ∂ ǫ ασ I δ h ασ ( τ ) δ ∆ h ασ ( t ) h g m ( τ , τ ′ ) + g m ( τ , τ 1 ) ˜ Σ m ( τ 1 , τ 2 ) g m ( τ 2 , τ ′ ) i dτ dτ ′ , δ W δ ∆ p ασ ( t ) = i I δ p ασ ( τ ) δ ∆ p ασ ( t ) h g ασ ( τ , τ ′ ) ˜ Σ ασ ( τ ′ , τ ) − ˜ Σ ασ ( τ , τ ′ ) g ασ ( τ ′ , τ ) i dτ dτ ′ . (B7) By expressing ˜ Σ = P n ˜ Σ ( n ) = Σ + Σ g Σ + Σ g Σ g Σ + . . . = g − 1 G Σ = Σ Gg − 1 with prop e r se lf-e nergy Σ, a nd performing the Keldysh r otation for the pro jection onto the real time, we finally obtain I ασ ( t ) = q ℜ Z ∞ −∞ dt ′  G K ασ ( t, t ′ )Σ A ασ ( t ′ , t ) + G R ασ ( t, t ′ )Σ K ασ ( t ′ , t )  , G ασ ( τ , τ ′ ) = Σ ασ ( τ , τ ′ ) + I Σ ασ ( τ , τ 1 ) G ασ ( τ 1 , τ 2 )Σ ασ ( τ 2 , τ ′ ) dτ 1 dτ 2 , h n ασ ( t ) i = i ¯ h X m ς m  ∂ ǫ m ∂ ǫ ασ  G < m ( t, t ) . (B8) F o r static cases, sinc e Green’s functions depend only on the difference b etw een the time a rguments, the expr ession for the cur r ent is further re duce d in energy representation of Eq. (32), I ασ = q 2 π ¯ h ℜ Z ∞ −∞ dE  G > ασ ( E )Σ < ασ ( E ) − G < ασ ( E )Σ > ασ ( E )  , = q 2 π ¯ h ℜ Z ∞ −∞ dE  g > ασ ( E )Σ < ασ ( E ) − g < ασ ( E )Σ > ασ ( E )  , (B9) where in the second line w e ma ke use o f Eq. (3 0). App endix C: Current conserv ation The current conserv a tion can be shown b y concen trating on o ne of n -th order diagrams in the generating function, which consist of n conduction or 2 n auxiliary particle Green’s functions. Each diagram can b e expressed either in terms of n conduction Green’s functions, or in terms of 2 n auxiliary particle Green’s functions . Since they represent the same diag ram, we can write X ασ 1 n I g ασ ( τ , τ ′ ) ˜ Σ ( n ) ασ ( τ ′ , τ ) dτ dτ ′ = − X m ς m 2 n I g m ( τ , τ ′ ) ˜ Σ ( n ) m ( τ ′ , τ ) dτ dτ ′ . (C1) By s umming all diag rams in the genera ting functional, we arrive at X ασ I G ασ ( τ , τ ′ )Σ ασ ( τ ′ , τ ) dτ ′ = − X m I ς m 2 G m ( τ , τ ′ )Σ m ( τ ′ , τ ) dτ ′ . (C2) 16 Additionally , by applying the Keldysh r otation ont o rea l time, the ab ov e relation beco mes X ασ Z dt ′  G K ασ ( t, t ′ )Σ A ασ ( t ′ , t ) + G R ασ ( t, t ′ )Σ K ασ ( t ′ , t )  = − X m ς m 2 Z dt ′  G K m ( t, t ′ )Σ A m ( t ′ , t ) + G R m ( t, t ′ )Σ K m ( t ′ , t )  . (C3) Using this rela tion, the sum of currents through bo th tunneling barr iers can be written as X ασ I ασ ( t ) = − q ℜ X m ς m 2 Z dt ′  G > m ( t, t ′ )Σ < m ( t ′ , t ) − G < m ( t, t ′ )Σ > m ( t ′ , t ) + G ++ m ( t, t ′ )Σ ++ m ( t ′ , t ) − Σ ++ m ( t, t ′ ) G ++ m ( t ′ , t )  . In static cas e, we obta in, by using the energy representation, X ασ I ασ = − q 2 π ¯ h ℜ X m ς m 2 Z dE  G > m ( E )Σ < m ( E ) − G < m ( E )Σ > m ( E )  = 0 , where Eq. (2 6) is used.