Non-equilibrium Transport in the Anderson model of a biased Quantum Dot: Scattering Bethe Ansatz Phenomenology
aa r X i v : . [ c ond - m a t . s t r- e l ] M a y Non-equilibrium Transport in the Anderson model of a biased Quantum Dot:Scattering Bethe Ansatz Phenomenology
Sung-Po Chao and Guillaume Palacios Center for Materials Theory, Department of Physics and Astronomy, Rutgers University, Piscataway, NJ 08854 Instituut voor Theoretische Fysica, Universiteit van Amsterdam,Valckenierstraat 65, 1018 XE Amsterdam, The Netherlands
We derive the transport properties of a quantum dot subject to a source-drain bias voltage atzero temperature and magnetic field. Using the Scattering Bethe Anstaz, a generalization of thetraditional Thermodynamic Bethe Ansatz to open systems out of equilibrium, we derive resultsfor the quantum dot occupation in and out of equilibrium and, by introducing phenomenological spin- and charge-fluctuation distribution functions in the computation of the current, obtain thedifferential conductance for large U Γ . The Hamiltonian to describe the quantum dot system is theAnderson impurity Hamiltonian and the current and dot occupation as a function of voltage areobtained numerically. We also vary the gate voltage and study the transition from the mixed valenceto the Kondo regime in the presence of a non-equilibrium current. We conclude with the difficultywe encounter in this model and possible way to solve them without resorting to a phenomenological method. PACS numbers: 72.63.Kv, 72.15.Qm, 72.10.Fk
I. INTRODUCTION
The past few years have witnessed a spectacularprogress in the fabrication and exploration of nano-structures giving experimentalists unprecedented controlover the microscopic parameters governing the physicsof these systems. Nano-structures, beyond their practi-cal applications, display an array of emergent phenom-ena stemming from their reduced dimensionality whichenhances quantum fluctuations and strong correlations.Often, experiments are carried out under non-equilibriumconditions, with currents passing through the structures.The measurements are performed over a wide range of pa-rameters, such as temperature and applied bias, allowingexperimental exploration of the interplay between non-equilibrium dynamics and strong correlation physics .A canonical example is the non-equilibrium Kondo effectobserved in a quantum dot attached to two leads held atdifferent chemical potentials µ i . The voltage difference V = µ − µ induces a non-equilibrium current I ( V )through the dot, interfering with and eventually destroy-ing the Kondo effect as the voltage is increased.The standard theoretical description of the transporttrough a quantum dot is the two-lead Anderson impuritymodel under a bias voltage. The 1- or 2-lead Hamilto-nian at zero bias is exactly solvable via Bethe Ansatz .Using this exact solution as well as NRG calculations forexample, the thermodynamics of the model have beenstudied in great detail. But the non-equilibrium situa-tion, namely when the two leads experience each a dif-ferent chemical potential, is much more difficult. This isdue to the subtle interplay between the non-equilibriumaspect of the problem on one hand and the presence ofstrong interactions on the other hand. Technically speak-ing it is very non trivial task to find a basis of states thatdiagonalize simultaneously the voltage term and the in-teraction term in the Hamiltonian. Nevertheless, a lot of efforts have been put forward to study this model but sofar, only approximate ways of dealing with the voltageand/or the interactions have been developed .In this paper we develop a phenomenological ap-proach to the problem, based on the Scattering BetheAnsatz (SBA), recently developed by P. Mehta andN. Andrei (MA) , a non-perturbative implementationof the Keldysh formalism to construct the current-carrying, open-system scattering eigenstates for the two-lead nonequilibrium Anderson impurity mode. The basicidea of the SBA is to construct scattering eigenstates ofthe full Hamiltonian defined directly on the infinite lineand match the incoming states by two Fermi seas describ-ing the initial state of the leads. The non-equilibriumsteady state transport properties of the system are thenexpressed as expectation values of the current or dotoccupation operators in these eigenstates. This pro-gram has been implemented for the Interacting Reso-nance Level Model (IRLM), a spinless interacting model,described in Ref. 7 where the zero temperature results forcurrent and dot occupation h ˆ n d i for all bias voltages werepresented. Another exact solution of this model at theso-called self-dual point by E. Boulat, H. Saleur andP. Schmitteckert in Refs. 25,26 uses conformal field the-ory techniques and compares successfully with t-DMRGresults.The main motivation of the present paper is to testthe very interesting ideas behind the SBA framework ona physically more relevant model such as the Andersonimpurity model and to focus on the phenomenology thatcan be extracted from it. Carrying out the program forthe non-equilibrium Anderson model we find difficultiesin the direct application of the SBA approach due to thefact that the ground state in the Bethe basis consists ofbound pairs of quasi-particles, leading to problems in thecomputation of the scattering phase shifts for the quasi-particles with complex momenta. This problem is notpresent in the IRLM when the Bethe momenta are belowthe impurity level and no bound states can be formed.We circumvent this difficulty by means of the follow-ing argument: The transport property computed in theIRLM is related to the single particle phase shift acrossthe impurity in the Bethe basis. Based on the same ideawe develop a phenomenological approach to describe thetransport property in the Anderson impurity model. Weidentify two types of possible phase shifts across impu-rity, which we refer to as ”spin-fluctuation” and ”charge-fluctuation” types to label two phenomenological phaseshifts akin to the fundamental excitations described inthe traditional Bethe Ansatz in this model. The phe-nomenological Ansatz is checked against exact results onthe dot occupation in equilibrium and the Friedel sumrule , in the linear response regime. Subsequently,we discuss our results for the out of equilibrium current,conductance and dot occupation. The scaling relationsfor the conductance, predicted from the Fermi liquid pic-ture of the problem at strong and weak coupling, are alsodiscussed.The paper is organized as follows. We start with aformal construction of scattering eigenstates in the two-lead Anderson impurity model. Then we discuss howwe impose boundary conditions, which serve as initialcondition in the time dependent picture, on the elec-trons within the leads. Next we shall discuss our re-sults for the dot occupation in equilibrium and the con-ductance in the linear response regime. Based on thechecks in equilibrium we then extend our computationto the out of equilibrium regime. The difficulty we en-counter for complex momenta and the way we handleit will also be addressed there. Comparison with an-other attempt of exact solution for this model by R.M. Konik et al. with the idea of dressed excitationsabove Fermi energy in the Bethe Ansatz picture, firstconsidered for the exact conductance of point contactdevice in the FQHE regime , will be discussed. Wewill also comment on the validity and implication of ournumerical results, among them the charge susceptibil-ity, in the out-of-equilibrium regime. Qualitative agree-ment between our theory and experimental result is thenpresented. The limit of U → ∞ is also summarized inthe last section based on the same phenomenological ap-proach. Finally, we summarize our results and concludewith some issues on the SBA approach to this model, andstate how they could be overcome. II. THE SCATTERING BETHE ANSATZAPPROACHA. Scattering state construction
In this section we apply the SBA approach to con-struct the scattering states of the full Hamiltonian. The (unfolded) 2-lead Anderson impurity Hamiltonian reads,ˆ H = X i =1 , Z d x ψ † iσ ( x )( − i∂ x ) ψ iσ ( x ) + ǫ d d † σ d σ + t i ( ψ † iσ (0) d σ + d † σ ψ iσ (0)) + U d †↑ d ↑ d †↓ d ↓ (1)where summation over the spin indices σ is implied. Thefields ψ iσ ( x ) describe chiral, right-moving electrons fromlead i , U is the on-site Coulomb repulsion between elec-trons on the dot, t i is the coupling between the dot andthe lead i , and ǫ d is the gate voltage. We have set theFermi velocity v F = 1.The model’s equilibrium properties have been studiedin great detail via the traditional Thermodynamic BetheAnsatz (TBA) . The SBA exploits in a new way theintegrability of the Anderson Model to construct current-carrying scattering eigenstates on the open line. Thereare two main requirements: One is the construction ofscattering eigenstates with the number of electrons ineach lead conserved prior to scattering off the impurity.Another is the asymptotic boundary condition: that thewave function of the incoming electrons, i.e. in the region( x ≪ . All information about the external bias ap-plied to the system is encoded in the boundary conditionby appropriately choosing the chemical potential of theincoming Fermi seas. As in all Bethe-Ansatz construc-tions, the full multi-particle wavefunction is constructedfrom single particle eigenstates (now on the infinite openline) and the appropriate two-particle S-matrices. Wefirst rewrite Eq. (1) in the even-odd basis asˆ H = ˆ H e + ˆ H o ˆ H e = X σ Z d x ψ † eσ ( x )( − i∂ x ) ψ eσ ( x ) + ǫ d d † σ d σ + t ( ψ † eσ (0) d σ + d † σ ψ eσ (0)) + U d †↑ d ↑ d †↓ d ↓ ˆ H o = X σ Z d x ψ † oσ ( x )( − i∂ x ) ψ oσ ( x )With ψ eσ ( x ) = t ψ σ ( x ) + t ψ σ ( x ) p t + t ψ oσ ( x ) = t ψ σ ( x ) − t ψ σ ( x ) p t + t and t = p t + t . In what follows we consider thecase t = t = t √ for simplicity. The single par-ticle solution for even and odd basis is: | e, pσ i = R d x ( e ipx g p ( x ) ψ † eσ ( x ) + e p δ ( x ) d † σ ) | i and | o, pσ i = R d x e ipx h p ( x ) ψ † oσ ( x ) | i , with | i the vacuum state and g p ( x ), h p ( x ), e p independent of spin and given by g p ( x ) = θ ( − x ) + e iδ p θ ( x ) + s ep θ ( x ) θ ( − x ) ,h p ( x ) = θ ( − x ) + θ ( x ) + s op θ ( x ) θ ( − x ) , (2) e p = t (1 + e iδ p + s ep / p − ǫ d ) . Here δ p ≡ − ( Γ ǫ d − p ) is the usual single particle scat-tering phase shift of the electrons off the impurity ob-tained when setting s ep = 0. Γ ≡ t is the width of theresonance level. We adopted a symmetric regularizationscheme θ ( ± x ) δ ( x ) = δ ( x ) and imposed | p | ≤ D , D be-ing the bandwidth cut-off . The s ( x ) = θ ( x ) θ ( − x ) termis a local constant ( ∂ x s ( x ) = 0) in this scheme and it isincluded in the odd channel function to allow the sametwo particle S-matrices, Eq.(4), in all channels . The θ ( x ) θ ( − x ) term in the even channel wave function is in-troduced in order to modify the original (when s ep = 0)single particle phase shift across the impurity. The choiceof s op and s ep will be addressed later. In the lead ba-sis, | i, pσ i , the single-particle scattering eigenstates withthe incoming particle incident from lead i , can be re-stored by taking a proper linear combination of even-oddstates. For example, | , pσ i = √ ( | e, pσ i + | o, pσ i ) = R d x e ipx α † ,pσ ( x ) | i is written as | , pσ i = Z d x e ipx n [ θ ( − x ) + 12 ( e iδ p + 1) θ ( x )] ψ † σ ( x )+ 12 ( e iδ p − θ ( x ) ψ † σ ( x ) + e p d † σ δ ( x ) + s † pσ ( x ) o | i (3)with | , pσ i = √ ( | e, pσ i − | o, pσ i ) = R d x e ipx α † ,pσ ( x ) | i and s † ipσ ( x ) related to the θ ( x ) θ ( − x ) terms by s † pσ ( x ) = (cid:18) s ep + s op √ ψ † σ ( x ) + s ep − s op √ ψ † σ ( x ) (cid:19) ×× θ ( x ) θ ( − x )and s † pσ ( x ) = (cid:18) s ep − s op √ ψ † σ ( x ) + s ep + s op √ ψ † σ ( x ) (cid:19) ×× θ ( x ) θ ( − x ) . These states have a single incoming particle ( x < i , that is reflected back into lead i with ampli-tude, R p = ( e iδ p + 1) / T p = ( e iδ p − /
2. Similar singleparticle states are discussed in Ref. 7.The multi-particle Bethe-Ansatz wave-function is con-structed by means of the two-particle S-matrix, S ( p, k ),describing the scattering of two electrons with momenta p and k . By choosing s op = − s ep will be discussed in section B and does not affect theresult here) in the single particle states we can constructthe same two-particles S-matrix for all combinations ineven-odd basis (see Appendix. B). The two-particles so-lution for both particles coming from lead 1 in spin singlet state takes the following form | k, ↑ ; 1 p, ↓i = Z d x d x A{ e i ( kx + px ) Z kp ( x − x ) α † k, ↑ ( x ) α † p, ↓ ( x ) }| i = n Z d x d x A { g ( x , x ) ψ † e ↑ ( x ) ψ † e ↓ ( x )+ h ( x , x ) ψ † o ↑ ( x ) ψ † o ↓ ( x ) + j ( x , x )( ψ † e ↑ ( x ) ψ † o ↓ ( x ) − ψ † e ↓ ( x ) ψ † o ↑ ( x )) } + Z d xA ( e ( x )( ψ † e ↑ ( x ) d †↓ − ψ † e ↓ ( x ) d †↑ )+ o ( x )( ψ † o ↑ ( x ) d †↓ − ψ † o ↓ ( x ) d †↑ )) + Am d †↑ d †↓ o | i Here A is the antisymmetrizer. A is an overall normal-ization factor and g ( x , x ) = Z kp ( x − x ) g k ( x ) g p ( x )+ Z kp ( x − x ) g k ( x ) g p ( x ) j ( x , x ) = Z kp ( x − x ) g k ( x ) h p ( x )+ Z kp ( x − x ) h k ( x ) g p ( x ) h ( x , x ) = Z kp ( x − x ) h k ( x ) h p ( x )+ Z kp ( x − x ) h k ( x ) h p ( x ) e ( x ) = Z kp ( − x ) g p ( x ) e k + Z kp ( x ) g k ( x ) e p o ( x ) = Z kp ( − x ) h p ( x ) e k + Z eokp ( x ) h k ( x ) e p m = ˜ Z kp (0) e k e p with Z kp ( x ) = e − iφ kp θ ( − x ) + e iφ kp θ ( x ) and ˜ Z kp (0) ≡ k + p − ǫ d k + p − U − ǫ d Z kp (0). Here tan( φ kp ) = − Ut ( k − p )( p + k − U − ǫ d ) .The derivation and more general form of two particlescase is written in Appendix. B. To include spin tripletcase we denote Z k i ,k j ( x i − x j ) ≡ Z k i ,k j ( x i − x j ) a ′ i a ′ j a i a j = I a ′ i a ′ j a i a j θ ( x j − x i ) + S a ′ i a ′ j a i a j ( k i , k j ) θ ( x i − x j ) where a i is thespin index before the scattering and a ′ i the spin indexafter the scattering. I a ′ i a ′ j a i a j is the identity matrix. TheS-matrices must satisfy the Yang-Baxter equations S a ′ a ′ a a ( k , k ) S a ′ a ′ a a ( k , k ) S a ′ a ′ a a ( k , k )= S a ′ a ′ a a ( k , k ) S a ′ a ′ a a ( k , k ) S a ′ a ′ a a ( k , k )for such a construction to be consistent. The two-particles S-matrix for this two-lead Anderson model isgiven by S a ′ i a ′ j a i a j ( k, p ) = ( B ( k ) − B ( p )) I a ′ i a ′ j a i a j + i U Γ P a ′ i a ′ j a i a j B ( k ) − B ( p ) + i U Γ (4)with B ( k ) = k ( k − ǫ d − U ), P = ( I · I + ~σ · ~σ ) thespin exchange operator with a i and a j representing theincoming spin indices. Since the S-matrix is the samefor all even-odd combinations the S-matrix does not de-pend on the lead index i , and the number of electrons ina lead, N i , can change only at the impurity site. Thiscircumstance allows us to construct the fully-interactingeigenstates of our Hamiltonian characterized by the in-coming quantum numbers, N and N the numbers ofincident electrons from lead 1 and 2 respectively. Thesequantum numbers are subsequently determined by thechemical potentials µ and µ .To complete the construction of the SBA current-carrying, scattering eigenstate, | Ψ , µ i i , we must stillchoose the ”Bethe-Ansatz momenta” { p l } N + N l =1 of thesingle particles states to ensure that the incoming parti-cles look like two Fermi seas in the region x <
0. Thisrequirement translates into a set of ”free-field” SBA equa-tions for the Bethe-Ansatz momenta-density of the par-ticles from the two leads . The argument is as follows:Away from the impurity | i, pσ i reduces to ψ † iσ ( x ) with theinter-particle S-matrix Eq. (4) present. Thus the scatter-ing eigenstates describing non-interacting electrons arein the Bethe basis rather than in the Fock basis of planewaves. The existence of many basis for the free elec-tron is due to their linear spectrum which leads to de-generacy of the energy eigenvalues. The wave function e ip x + ip x [ θ ( x − x ) + S θ ( x − x )] A is an eigenstateof the free Hamiltonian for any choice of S with, in par-ticular, S = defining the Fock basis and S given inEq. (4) defining the Bethe basis. The Bethe basis is thecorrect ”zero order” choice of a basis in the degenerateenergy space required in order to turn on the interac-tions. We proceed to describe the leads (two free Fermiseas) in this basis.We consider the system at zero temperature and zeromagnetic field in this paper. To describe the two Fermiseas on the leads translates to a set of Bethe Ansatzequations whose solution in this case consists of com-plex conjugate pairs: p ± ( λ ) = x ( λ ) ± iy ( λ ) in the λ -parametrization with x ( λ ) = ˜ ǫ d − s λ + ˜ ǫ d + p ( λ + ˜ ǫ d ) + U Γ y ( λ ) = − s − ( λ + ˜ ǫ d ) + p ( λ + ˜ ǫ d ) + U Γ . with ˜ ǫ d = ǫ d + U/
2. Each member of a pair can be ei-ther in lead 1 or in lead 2, since the S-matrix is unity inthe lead space. There are, therefore, two possible con-figurations for these bounded pairs. One possible wayof forming bounded pairs is described by four types ofcomplex solutions whose densities we denote σ ij ( λ ) with { ij } = { , , , } indicating the incoming electronsfrom lead i and lead j . The other possibility, which isperhaps more intuitive in comparing with the free elec-tron in the Fock basis, is to include only { ij } = { , } .These two types of states give the same results whenevaluating the expectation value of the dot occupation inequilibrium. However when we turn on the bias voltage,the results obtained from a 4-bound states descriptionshow some charge fluctuations even way below the impu-rity level which is not expected from the non-interacting ( U →
0) theory (shown in Appendix A). Thus we shalldisregard the 4-bound states solution on physical groundsand focus on the 2-bound states description in the follow-ing discussion.To describe in the Bethe basis the two leads as twoFermi seas filled up to µ and µ , respectively, these den-sities must satisfy the SBA equations,2 σ i ( λ ) = − π d x ( λ )d λ θ ( λ − B i ) − X j =1 , Z ∞ B j d λ ′ K ( λ − λ ′ ) σ j ( λ ′ ) (6)with K ( λ ) = π U Γ(2 U Γ) + λ .The SBA equations are derived from imposing bound-ary condition in the free leads (incoming state) regionand the value of momenta is connected with spin rapid-ity λ by using the quantum inverse scattering method.The Bethe Ansatz equations solved with periodic bound-ary conditions at the free lead region with total numberof particles N ( N = N + N as sum of particle num-ber from lead 1 and 2) and the total spin projection S ( S = S + S = N/ − M with M = M + M as numberof down spin particles from lead 1 and 2) are given by e ik lj L = M Y α =1 B ( k lj ) − λ α + iU Γ B ( k lj ) − λ α − iU Γ (7) Y l =1 , N l Y j =1 B ( k lj ) − λ α − iU Γ B ( k lj ) − λ α + iU Γ = M Y β = α λ α − λ β + 2 iU Γ λ α − λ β − iU Γwith total energy E = E + E and E l = P j k lj indicat-ing the energy of the electrons within the lead l at zerotemperature.The spectrum of Eq.(7) for one lead case has been an-alyzed by N. Kawakami and A. Okiji where they foundthat the ground state at zero temperature is composedof real λ i and complex k lj in the thermodynamic limitfor U >
0. The same situation also occurs in the speciallimit where U → ∞ where P. Schlottmann has donealso in the one lead case. The proof for two leads groundstate is similar to the one lead case and is shown explic-itly for the finite temperature calculation for the infiniteU case in Ref. 36.As has been mentioned above in the zero temperaturezero magnetic field ground state all λ i are real (and dis-tinct) and k lj form bound state for j = 1 , .., M withbound state momenta given by the poles or zeros in theS-matrix defined in Eq. (4) B ( k l ± ( λ j )) = λ j ± iU Γ = B ( x ( λ j ) ± iy ( λ j )) + γ ± ( λ j )(8)where γ ± = O (exp( − L )) and x ( λ ) and y ( λ ) are shownin the Eq. (5).Note that the bound state can be formed from fourpossible configurations for B < λ α < ∞ which we de-note bound state from lead i and lead j quasi momentadenoted as λ ijα . The bound state between B < λ α < B can only be formed by quasi momenta both coming fromlead 1. As already mentioned the four bound state dis-tribution does not give physically sensible results for thecharge susceptibility as shown in Appendix A and there-fore we will limit our discussion to two types of boundstate distribution here. Below we surpass the index oflead in λ and put back the index dependence in the endfor simplification. Inserting Eq. (8) into Eq. (7) we get e ik + α L = M Y β =1 λ α − λ β + 2 iU Γ λ α − λ β + γ + α (9) e ik − α L = M Y β =1 λ α − λ β + γ − α λ α − λ β − iU Γ (10) M Y β =1 λ β − λ α + γ + β λ β − λ α + γ − β = 1 (11)Thus for L → ∞ from multiplication of Eq. (9) andEq. (10) we have e ix ( λ α ) L = Y β λ α − λ β + 2 iU Γ λ α − λ β − iU Γ (12)Taking the logarithm of Eq. (12) we have:2 πJ α = − x ( λ α ) L − X β (cid:18) θ (cid:18) λ α − λ β U Γ (cid:19) + π (cid:19) (13)with θ n ( x ) ≡ tan − (2 x/n ) and { J α } a set of integer num-bers. We can extend the definition of J α to include inte-gers or half integers and rewrite Eq. (13) as πL J α = − x ( λ α ) − L X β θ (cid:18) λ α − λ β U Γ (cid:19) (14)Now let us put back the dependence in lead indices.Starting from Eq. (14) it can be shown that there isone-to-one correspondence between the λ α ’s and the J α ’sand that all λ α ’s have to be different. Thus the setof rapidities { λ ijα } , characterizing an eigenstate of theHamiltonian, is uniquely determined by one specific setof { J α } . For instance, the ground state of the Hamilto-nian H in the presence of a bias voltage is simply ob-tained by packing two ”Fermi seas” of non-consecutiveintegers (Pauli principle in lead space) up to certain”Fermi points” (see Fig. 1 ) corresponding to the B and B in the continuum limit. For notational simpli-fication we relabel { ij } = { , } as { l } = { , } . Nowdefining P ij σ ( λ ijα ) = L d J α d λ α ≡ P l σ ( l ) ( λ α ) and using ∂ x θ n ( x ) = /n x/n ) we can write Eq. (14) in the contin-uum limit (by taking L → ∞ and differentiate Eq. (14)with respect to λ ). Doing so we shall distinguish twodifferent domains:For B < λ < ∞ the particles are fully packed andstates are labeled by a different lead index l . In this FIG. 1: Sketch of the configuration of Bethe momenta corre-sponding to the ground state of H with an additional biasvoltage i.e. two Fermi seas at different chemical potential. domain, the SBA equations in the continuum limit takesthe form X l =1 σ ( l ) ( λ ) = − π d x ( λ )d λ − Z ∞ B d λ ′ K ( λ − λ ′ ) σ (2) ( λ ′ ) − Z ∞ B d λ ′ K ( λ − λ ′ ) σ (1) ( λ ′ ) . (15)For B < λ < B we can see from Fig. 1 that the lead 2states are unoccupied. We shall introduce a distributionof holes for the lead 2 that we will denote ˜ σ (2) ( λ ). Thecontinuum SBA equations in this regime are given by σ (1) ( λ )+˜ σ (2) ( λ ) = − π d x ( λ )d λ − Z ∞ B d λ ′ K ( λ − λ ′ ) σ (2) ( λ ′ ) − Z ∞ B d λ ′ K ( λ − λ ′ ) σ (1) ( λ ′ ) (16)Since ˜ σ (2) ( λ ) obeys the same equation as σ (2) ( λ ) as maybe seen from subtracting Eq. (16) and Eq. (15) we cancombine Eq. (15) and Eq.(16) together to get2 σ ( λ ) = − π d x ( λ )d λ − Z ∞ B d λ ′ K ( λ − λ ′ ) σ ( λ ′ ) − Z B B d λ ′ K ( λ − λ ′ ) σ ( λ ′ ) (17)with B < λ < ∞ for lead 1 and B < λ < ∞ for lead2 Bethe momenta density distributions. Each density isdefined on a domain extending from B i to the cutoff D - to be sent to infinity. The B i play the role of chemicalpotentials for the Bethe-Ansatz momenta and are deter-mined from the physical chemical potentials of the twoleads, µ i , by minimizing the charge free energy , F = X i ( E i − µ i N i ) = 2 X i Z ∞ B i d λ ( x ( λ ) − µ i ) σ i ( λ )(18)with σ the lead 1 particle density and σ the lead 2particle density. Note that σ and σ obeys the sameintegral equation Eq. (6) with different boundary ( σ ( λ )with λ ⊂ ( B , ∞ ) and σ ( λ ) with λ ⊂ ( B , ∞ )). Solvingthe SBA equations subject to the minimization of thecharge free energy fully determines the current-carryingeigenstate, | Ψ , µ i i and allows for calculation of physicalquantities by evaluating expectation value of the corre-sponding operators. In the following we shall discuss ourresults from equilibrium cases to non-equilibrium ones,starting with the expression for various expectation valueof physical quantities. B. Expectation value of current and dot occupation
For µ = µ all B i are equal to some equilibriumboundary B fixed by the choice of µ i . The dot occupationis given by the expectation value P σ h Ψ , µ i | d † σ d σ | Ψ , µ i i .Taking the limit L → ∞ ( L being the size of thelead) one can express n d as an integral over the den-sity of λ and the corresponding matrix element ν ( λ ) ≃ h p + ( λ ) p − ( λ ) | P σ d † σ d σ | p + ( λ ) p − ( λ ) ih p + ( λ ) p − ( λ ) | p + ( λ ) p − ( λ ) i taken to order L . Herethe state | p + ( λ ) p − ( λ ) i denotes a pair (or bound states) ofquasi-particles with complex momenta given by Eq. (5).The reason why n d is governed solely by one-bound statematrix element instead of a complicated many-particle object is because Bethe wave-functions are orthogonal toeach other for different pairs of Bethe momenta under thecondition that the size of the leads L taken to infinity.Here we address the different choice of s ep (with s op = − ν ( λ ). We shall first discussthe ” natural ” choice s ep = 0 (i.e. absence of θ ( − x ) θ ( x )terms) and show it reproduces the exact result for thedot occupation in equilibrium. While in checking thesteady state condition, i.e. d h n d i /dt = 0, for out-of-equilibrium situation, the choice of s ep = 0 fails. Toremedy this issue we propose s ep = 0 (i.e introducingcounter-intuitive θ ( − x ) θ ( x ) terms) schemes to circum-vent this difficulty. We check this proposed phenomeno-logical scheme in equilibrium against the exact dot occu-pation obtained in s ep = 0 case in the second part of thediscussion as a benchmark for our approach. First let usdiscuss the result for s ep = 0:(1) s ep = 0: We choose s ep = 0 as in the case of the 1-lead Anderson impurity model. Denote ν ( λ ) = ν SBA ( λ )in this choice. The dot occupation expectation value inequilibrium is given by n d = h Ψ , µ = µ | P σ ˆ d † σ ˆ d σ | Ψ , µ = µ ih Ψ , µ = µ | Ψ , µ = µ i = 2 Z ∞ B d λ σ ( λ ) ν SBA ( λ ) (19)where the factor 2 in front of the integral accounts forthe spin degeneracy. The matrix element of the operator d † σ d σ in the SBA state is given by ν SBA ( λ ) = 2Γ˜ x ( λ ) + ˜ y ( λ ) + 16 y ( λ )Γ [˜ x ( λ ) + ˜ y − ( λ )][˜ x ( λ ) + ˜ y ( λ )] (cid:18) ˜ x ( λ )2˜ x ( λ ) − U (cid:19) . where we introduced, for simplified notations, the func-tions ˜ x ( λ ) = x ( λ ) − ǫ d and ˜ y ± ( λ ) = y ( λ ) ± Γ.Eq. (19) can be proved to be exact by comparing itwith the traditional Bethe Ansatz (TBA) result. In thelatter, n d is computed as the integral of the impuritydensity. This observation that the SBA and TBA resultsfor n d agree in equilibrium shows the connection betweenthe dot occupation and the dressed phase shift acrossthe impurity. The dressed phase shift mentioned here isequivalent to the impurity density as can be seen in theEq.(C1) in Appendix C. The proof of the equivalencebetween TBA and SBA in equilibrium is also given inAppendix C.To describe the out-of-equilibrium state we first checkif the steady state condition d h ˆ n d i dt = 0 (or equivalently, d h ˆ N + ˆ N i dt = 0) is satisfied in this basis. As mentioned ear- lier these scattering states are formed by bounded quasi-particles with complex momenta and therefore the singleparticle phase across the impurity is not well defined inthe sense that | e iδ p ± | 6 = 1. This problem begins to sur-face as we set out to evaluate transport expectation valueand renders d h ˆ n d i dt = Z B B dλσ b ( λ )∆( λ ) = 0 (20)with ∆( λ ) = y ( λ )Γ [˜ x ( λ ) + ˜ y − ( λ )][˜ x ( λ ) + ˜ y ( λ )] . Thus it appears that using this basis the steady statecondition is not observed. This problem does not appearwhen the momenta are real as in the IRLM case .(2) s ep = 0: To remedy this problem we redefine thesingle particle phase shifts across the impurity, in anal-ogy to the results for the IRLM , through the choiceof nonzero s ep in Eq.(3). With a suitable choice of s ep we may restore a well defined single particle phase | e i ˜ δ p ± | = 1 with ˜ δ p ± denoting this new phase. The waywe judge whether we make the correct choice for the newphases ˜ δ p ± is to compare the dot occupation n d in equi-librium before and after the redefined phase. The ex-plicit form of s ep and phase ˜ δ p ± will be motivated belowbut first we shall show that a single redefined phase isnot sufficient to satisfy the constraint of dot occupationcomparison.Again the choice of new phases is constrained by therequirement that we shall obtain the same result for h P σ d † σ d σ i as given by ν SBA ( λ ) in equilibrium. Basedon this constraint it can be shown explicitly that a singlewell defined phase (in the sense of | e i ˜ δ p ± | = 1) is not suffi-cient to reproduce the equilibrium ν SBA ( λ ) as following:The new dot amplitude ˜ e p + and ˜ e p − have to satisfy | ˜ e p + | + | ˜ e p − | = 4Γ˜ x ( λ ) + ˜ y ( λ )) , | ˜ e p + | | ˜ e p − | = 4Γ [˜ x ( λ ) + ˜ y ( λ )][˜ x ( λ ) + ˜ y − ( λ )] . As both | ˜ e p + | and | ˜ e p − | are positive we see that a sin-gle redefined phase cannot satisfy the above constraintssimultaneously. Therefore we have to choose at least twosets of redefined phases ˜ δ ip ± (with i = s, h denoting spin-fluctuation or charge-fluctuation to be addressed later)and, along with them, some distribution functions f i toset the weight for these phases.To motivate the idea of searching the correct phaseshifts we shall come back to the derivation of dot occupa-tion in traditional Bethe Ansatz (TBA) picture. In TBAthe total energy of the system is described by energy ofthe leads electrons and energy shifts from the impurity, E = X j p j = X j (cid:18) πn j L + 1 L δ j (cid:19) (21)Based on Feynman-Hellman theorem, which is applicablein equilibrium (closed) system, we have h ˆ n d i = ∂E∂ǫ d = 1 L X j ∂δ j ∂ǫ d = 1 L X j ∂ ( δ p + j + δ p − j ) ∂ǫ d (22)The result for Eq. (22) agrees with those obtained fromEq. (C2) and can be viewed as a third approach to ob-tain the expectation value of the dot occupation. Thekey observation here is that this quantity is related tothe bare phase shift δ p + + δ p − and therefore the redefinedphases must be proportional to this quantity. Amongthem there are two likely candidates with redefined phaseshift given by δ p + + δ p − , describing the tunneling of abounded pair, and δ p + + δ p − , describing the tunneling of a single quasi-particle. In a sense this is the echo for the ele-mentary excitations above the Fermi surface in the Bethebasis characterized by N. Kawakami and A. Okiji ascharge-fluctuation excitation, which describes boundedpair quasi-particles excitation, and spin-fluctuation exci-tation, which describes one quasi-particle excitation. An-other similar picture is the spin-fluctuation and charge-fluctuation two fluids picture proposed by D. Lee et al albeit in a different context. We identify the phase de-fined by ˜ δ p − = ˜ δ p + = δ p + + δ p − ≡ ˜ δ sp (with s ep ± ≡ s sep ± = ( i ( p ± − ǫ d ) − Γ)( e i ( δp + + δp − ) − δ p − = ˜ δ p + = δ p + + δ p − ≡ ˜ δ hp (with s ep ± ≡ s hep ± = ( i ( p ± − ǫ d ) − Γ)( e i ( δ p + + δ p − ) − I with h ˆ I i definedby h ˆ I i = −√ iet ~ h X σ (( ψ † σ (0 ± ) − ψ † σ (0 ± )) d σ − h.c. ) i (23)in the state | Ψ , µ i i . Notice that ψ † iσ (0 ± ) ≡ lim ǫ → ( ψ † iσ ( − ǫ ) + ψ † iσ (+ ǫ )) / δ sp and˜ δ hp we have the expression for current as I ( µ , µ ) = h Ψ , µ , µ | ˆ I | Ψ , µ , µ i = 2 e ~ Z B B d λ σ b ( λ )( f s ( λ ) J s ( λ ) + f h ( λ ) J h ( λ )) (24)The corresponding spin-fluctuation and charge-fluctuation matrix element of the current operator basedon the spirit of Landauer transport, denoted as J s ( λ ) and J h ( λ ) with J α ( λ ) = | T p ( λ ) | = | e i ˜ δαp − | ( α = { s, h } ) de-pending on redefined phase shift ˜ δ αp only, are given by J s ( λ ) = 1 + sgn(˜ x ( λ ))(˜ x ( λ ) + y ( λ ) − Γ ) p (˜ x ( λ ) + y ( λ ) − Γ ) + 4Γ ˜ x ( λ ) (25) J h ( λ ) = 2Γ ˜ x ( λ )(˜ x ( λ ) + Γ ) − y ( λ )(Γ − ˜ x ( λ )) + y ( λ ) . (26)Here sgn( x ) = x | x | is the sign function. It is introducedin order to pick up the correct branch when taking thesquare root in denominator of Eq. (25). This way we en-sure that J s ( λ ) has the proper limit when U is sent toinfinity (cf Section III). Other than the motivations men-tioned above for identifying spin and charge fluctuationphase shifts the functional forms of J s ( λ ) and J h ( λ ) as afunction of bare energy x ( λ ) can also be used to identifythese two type of phase shifts (See Fig. 12 in Section IIIfor infinite U Anderson model, the finite U is similar).Next we shall choose the appropriate weight for eachtype of phase shift. So far we have not yet been ableto deduce the form of these weight functions f s ( λ ) and f h ( λ ) and we introduce them phenomenologically . Letus define phenomenological spin-fluctuation and charge-fluctuation weight functions as f s ( ε ( λ )) = D s ( ε ( λ )) D s ( ε ( λ )) + D h ( ε ( λ )) (27)and f h ( ε ( λ )) = D h ( ε ( λ )) D s ( ε ( λ )) + D h ( ε ( λ )) . (28)Here D s ( ε ( λ )) is the spin-fluctuation density of state, D h ( ε ( λ )) is the charge-fluctuation density of state as de-fined in Ref. 40, and ε ( λ ) is the corresponding dressedenergy i.e. the energy required to produce these spin-and charge-fluctuation excitations above the Fermi level.Here dressed energy refers to the sum of the bare en-ergy of adding/removing one bound state, as in chargefluctuation, or single quasi particle, as in spin fluctuation,and the energy shift from other quasi particles due to thischange. The equation that solves a single quasi-particle’sdressed energy ε ( λ ) reads ε ( λ ) = ( x ( λ ) − µ ) − Z ∞ B d λ ′ K ( λ − λ ′ ) ε ( λ ′ ) . (29)We wish to compare at this point our approach to theone taken by Konik et al . The authors’ Landauerapproach is based on an ensemble of renormalized exci-tations, the holons and spinons, and the conductance isexpressed in terms of their phase shift crossing the impu-rity. However, the leads are built of bare electrons andthus one faces the difficult problem of how to constructa bare electron out of renormalized excitations in orderto be able to impose the voltage boundary condition.The basic approximation adopted, electron ≈ antiholon+ spinon , is valid only when the electron is close to theFermi surface (see N. Andrei ), and therefore the ap-proach is trustworthy only for very small voltages. Nev-ertheless, the dressed excitations framework seems to give at least qualitatively good results when another energyscale (such as the temperature or an external field) isturned on . In contrast we construct the eigenstates ofthe Hamiltonian directly in terms of the bare electronfield and can therefore impose the asymptotic boundarycondition that the wave function tend to a product oftwo free Fermi seas composed of bare electrons. Whilewe do not have a mathematically rigorous derivation ofthe weight functions we introduced, the validity of thescattering formalism is not restricted to any energy win-dow other than energy cutoff. C. Results for equilibrium and linear response
In the numerical computation, for the practical pur-pose, we assumed Kondo limit ( U = − ǫ d , U Γ ≫
1) formof the spin-fluctuation and charge-fluctuation distribu-tions, i.e. D s ( ε ( λ )) ≃ π T k ε ( λ ) + T k (30)and D h ( ε ( λ )) ≃ √ U Γ Γ ( ε ( λ ) + ǫ d ) + Γ (31)with T k being the Kondo scale derived in Ref. 40 as T k = √ U Γ π e π ǫd ( ǫd + U )+Γ22 U Γ . (32)As we use the Kondo limit in our expression for thespin-fluctuation and charge-fluctuation distributions, weexpect our phenomenological approach works better forlarge U/ Γ. We also take ε ( λ ) ≃ x ( B ) − x ( λ ) for numer-ical convenience with B denoting the Bethe momentaboundary given by µ = µ = 0. The dot occupation h P σ d † σ d σ i evaluated by these new phases is given by h X σ d † σ d σ i = 2 Z ∞ B d λ σ b ( λ )( ν s ( λ ) f s ( λ )+ ν h ( λ ) f h ( λ ))+ Z ∞ B d λ σ b ( λ )( ν s ( λ ) f s ( λ ) + ν h ( λ ) f h ( λ )) ! (33)with ν s ( λ ) and ν h ( λ ) given as ν s ( λ ) = 1Γ " − (˜ x ( λ ) + y ( λ ) − Γ ) p (˜ x ( λ ) + y ( λ ) − Γ ) + 4Γ ˜ x ( λ ) × " y ( λ ) 1Γ − (˜ x ( λ ) + y ( λ ) − Γ ) p (˜ x ( λ ) + y ( λ ) − Γ ) + 4Γ ˜ x ( λ ) ! (cid:18) ˜ x ( λ )2˜ x ( λ ) − U (cid:19) (34) ν h ( λ ) = (cid:20) x ( λ )(˜ x ( λ ) + Γ ) − y ( λ )(Γ − ˜ x ( λ )) + y ( λ ) (cid:21) × " y ( λ )Γ˜ x ( λ )(˜ x ( λ ) + Γ ) − y ( λ )(Γ − ˜ x ( λ )) + y ( λ ) (cid:18) ˜ x ( λ )2˜ x ( λ ) − U (cid:19) (35)respectively. We may check whether this choice of phe-nomenological distribution functions satisfy the conditionin equilibrium that h X σ d † σ d σ i = 4 Z ∞ B d λ σ b ( λ ) ν SBA ( λ )= 4 (cid:18)Z ∞ B d λ σ b ( λ )( ν s ( λ ) f s ( λ ) + ν h ( λ ) f h ( λ )) (cid:19) . (36)We can see from the Top of Fig. 2 that the compari-son between the phenomenological and the exact resultfor the dot occupation in equilibrium is good deep intothe Kondo regime ( ǫ d ≃ − U ) and far away from it( ǫ d ≫
0) but is worse when we are in mixed valenceregion ( ǫ d ≃ D s ( ε ) and D h ( ε ), may go awayif we took more realistic form of D s ( ε ( λ )) and D h ( ε ( λ ))also in mixed valence regime as suggested in Fig. 2. How-ever the numerical procedure is much more complicatedthere. We confine ourself to this simpler limit in ourphenomenological approach.Another check on our result in equilibrium is to find thelinear response conductance through our formulation andcompare with the exact linear result given by the Friedelsum rule . The Friedel sum rule, which relates theequilibrium dot occupation to the phase shift experiencedby electrons crossing the dot, is related to zero voltageconductance by dIdV | V =0 = 2 sin ( π h ˆ n d i / by noting that at low-voltage eV = µ − µ ≃ πL ( N − N ) = 4 π R B B σ b ( λ ) dλ . By taking B ≃ B = B in theexpression for the current across the impurity Eq. (24)we get the zero bias conductance expressed as dIdV (cid:12)(cid:12)(cid:12) V =0 = e h (cid:2) f s ( B ) J s ( B ) + f h ( B ) J h ( B ) (cid:3) (37)Here B = B ( µ, ǫ d , Γ , U ) is determined by µ = µ = 0.The comparison between Friedel sum rule (FSR) re-sult and the conductance given by Eq. (37) (denoted as (pSBA)) is shown at the Bottom of Fig. 2. It displaysthe consequence of the equilibrium Kondo effect in thequantum dot set up: due to the formation of the Kondopeak attached to the Fermi level the Coulomb blockadeis lifted and a unitary conductance is reached for a rangeof gate voltages ǫ d around − U/
2. Again we see that thecomparison is good for large U/ Γ but poorer in mixedvalence regime for smaller U/ Γ, which is consistent withthe observation we made when evaluating h ˆ n d i as shownin top figure of Fig. 2. Having checked our results inequilibrium we shall go on to compute the current andthe dot occupation in the out-of-equilibrium regime. D. Results Out-Of-Equilibrium
Now let us begin to investigate the current and dot oc-cupation change as we turn on the voltage. We start withthe discussion on current vs voltage for various regime.The current vs voltage is plotted in the inset of figure ofFig. 3 for different values of U and at the symmetric point ǫ d = − U/
2. Note that we use an asymmetric bias voltagewhen solving numerically the integral equations originat-ing from Eq. (6) with constraint of minimizing the chargefree energy Eq. (18): Namely we fix µ ≃ − − − ) and lower µ . Therefore, a direct confronta-tion between the results obtained from real-time simula-tions of the Anderson model out-of-equilibrium isdifficult but the main features of our calculation matchthe predicted results: a linear behavior of the I - V charac-teristics at low-voltage, the slope being obtained from theFSR (2 in units of e /h at the symmetric point), and anon-monotonic behavior at higher voltage, the so-callednon-linear regime. In particular, our calculations showclearly that the current will decrease as U/ Γ is increasedwhich is in agreement with other numerical approaches(e.g. cf Fig. 2 of Ref. 20 for a comparison).The plots of the differential conductance vs sourcedrain voltage for different dot levels, ǫ d , tunnelingstrengths Γ and interaction strengths U are shown inFig. 3 and Fig. 5. Two major features emerge from these0 -4 -2 0 2 40.00.20.40.60.81.0 < n d > d =0.5(TBA) =0.5(pSBA) =0.25(TBA) =0.25(pSBA) =0.1(TBA) =0.1(pSBA) -4 -2 0 2 40.00.51.01.52.0 d I/ d V ( e / h ) d =0.5(FSR) =0.5(pSBA) =0.25(FSR) =0.25(pSBA) =0.1(FSR) =0.1(pSBA) FIG. 2: Top: h ˆ n d i as a function of ǫ d from the exact result(dotted line) and from Eq. (36) (solid line). Bottom: Thedifferential conductance in the linear-response regime, as afunction of ǫ d from the phenomenological Scattering BetheAnsatz (pSBA) and exact linear response conductance fromFriedel sum rule (FSR) for Γ = 0 .
5, 0 .
25, 0 .
1, and U = 8. plots: 1) A narrow peak around zero bias reaching max-imal value of 2 e /h (the unitary limit) for values of thegate voltage close to the symmetric point ( ǫ d ≃ − U/ ǫ d is raised from the Kondo regime, ǫ d = − U/
2, to themixed valence regime, ǫ d = 0, with the Kondo effect dis-appearing. As a function of the bias the various curvesdescribing the Kondo peak for different values of the pa-rameters can be collapsed onto a single universal functiond I/ d V = d I/ d V ( V /T ∗ k ) as shown in Fig. 4. Here T ∗ k is defined as T ∗ k = c √ U Γ π e ǫd ( ǫd + U )+Γ22 U Γ (38)with c = 0 . T ∗ k was ex-tracted from the numerics by requiring that the functiond I/ d V ( V /T ∗ k ) decreases to half its maximal value when V ≃ T ∗ k . The expression for T ∗ k as given by Eq. (38) dif-fers from the thermodynamic T k as defined in Eq. (32).The difference of prefactor in the exponential is certainlyrelated to the unusual choice of regularization scheme inthe SBA . The other possible implication for this dif-ferent formulation for the Kondo scale is also addressedlater when we discuss the experiment done by L. Kouwen-hoven et al . U=8, =1 U=12, =1 U=16, =1 U=24, =1 d I/ d V ( e / h ) V/ I V/ U/ =8 U/ =12 U/ =16 U/ =24 Slope=2
FIG. 3: dI/dV vs V /
Γ for Γ = 1, ǫ d = − U/
2, and various U . Inset: Steady state current vs voltage curves for Γ = 1, ǫ d = − U/
2, and various U . Dashed line is a line with constantconductance e h plotted for comparison. The small voltage behavior for differential conductancein symmetric case, i.e. ǫ d ≃ − U , is expected to be dIdV (cid:12)(cid:12)(cid:12) V ≪ T ∗ k ≃ e h − α V (cid:18) VT ∗ k (cid:19) ! and allows us to identify the constant α V from thequadratic deviation from 2 e /h . The quadratic fit of theuniversal curve around V ≃
0, as shown in Fig. 4, gives α V ≃
1. It is also expected for T ∗ k ≪ V ≪ U that thetail of the peak decays logarithmically as dIdV ∼ e h ( VT ∗ k ) . The latter behavior is observed (see inset of Fig. 4 ) in theregime U Γ ≫ < VT ∗ k < with the logarithmic1 -4 -4 -4 U=15 U=17 U=18 d I/ d V ( e / h ) V/ d I/ d V ( e / h ) V/T *K U=15 U=17 U=18 Fit: Y=2-2X -1 U=15 U=17 U=18 d I/ d V ( e / h ) V/T *K d I/ d V ( e / h ) *K )) U=15 U=17 U=18 Y=0.01+0.055x Y=0.025+0.055x Y=0.095+0.055x FIG. 4: Top: Zoomed in picture of the differential conduc-tance vs voltage nearby zero voltage. Inset shows the univer-sality in conductance vs voltage scaled by T ∗ k when VT ∗ k ≤ VT ∗ k < . T ∗ k nearby the Kondo peak structure. Inset showsthe logarithmic behavior when VT ∗ k ≫
1. Γ = 0 . function given by dIdV = e h " f (cid:18) U Γ (cid:19) + c ln ( VT ∗ k ) with the parameter c = 0 . f ( U Γ ) is simplya constant (in V ) shift. As suggested from the bottomplot of Fig. 4 (see also Fig. 14 for the infinite U case) thecharge fluctuation side peak does not fall into the samescaling relation but the strong correlations shift the cen-ter of the side peak closer to V = 0 (see Fig. 3 and Fig. 5).In other words the position of the resonance in the dI/dV curve naively expected around V = | ǫ d | is renormalized by the presence of interactions. In the inset of Fig. 5 weshow the logarithm of the voltage obtained at half widthhalf maximum (HWHM) of the zero voltage peak and -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 -4.5-4.0-3.5-3.0-2.5-2.0-1.5-1.0 e =-4 e =-2.5 e =-1 e =-0.5 e =-0 e =0.5 d I/ d V ( e / h ) V/4 a r b un i t Log(T *K )+6.215 Log(HWHM) d /4 FIG. 5: dI/dV vs V /
4Γ for U = 8, Γ = 0 .
25 and various ǫ d from Kondo ( ǫ d = −
4) to mixed valence regime ( ǫ d ≃ T ∗ k ) − ln( c ) and ln( V HWHM ) as afunction of impurity level ǫ d . Here V HWHM is the voltagedifference estimated at half value of differential conductanceat zero voltage. The constant shift − ln( c ) is chosen to givethe best fit in the data away from ǫ d = − U . compare it withln T ∗ k = ǫ d ( ǫ d + U ) + Γ U Γ + ln c √ U Γ π ! (after subtracting the constant ln c ). What is impor-tant and universal is that both quantities (ln V HWHM and ln T ∗ k ) exhibit a quadratic behavior in the gate volt-age ǫ d . Similar results had been found experimentallyby L. Kouwenhoven et al when they compare the fullwidth half maximum of dI/dV (from which they obtaina Kondo scale T k at finite voltage) with the tempera-ture dependence of the linear response differential con-ductance (from which another Kondo scale T k is ex-tracted). It is suggested from our numerical results thatboth ln T k (in analogy with our T k ) and ln T k (whichis our T ∗ k ) follows similar quadratic behavior in ǫ d butdiffer in their curvatures by a factor of π . In Ref. 5 thecurvatures of the quadratic behavior differ by a factor ofaround 2 (see Fig.3B in Ref. 5) which is attributed todephasing of spin fluctuations at finite voltage.Notice that in all the numerical data shown for cur-rent vs voltage we have chosen U Γ ≥ phenomenological distribution functions intro-duced to control the relative weight for spin- and charge-fluctuation contributions work is much better in the large U Γ regime (cf. Fig. 2).Next let us study the change in the dot occupation asa function of the voltage. The extension of the com-putation of the dot occupation out of equilibrium isstraightforward. Suppose we find the correct distribu-tion functions f s ( λ ) and f h ( λ ) then we have ν SBA ( λ ) =2 U=8U=12U=16U=24 < nd > V/ - d < nd > / d V V/ FIG. 6: h ˆ n d i vs V /
Γ for different U with ǫ d = − U and Γ = 1case. Inset: The corresponding nonequilibrium charge sus-ceptibility. A small peak shows up nearby V = 0 for all thesecurves. - d < nd > / d V V/4 d =-4 d =-3.5 d =-2 d =-0.5 d =0 d =0.5 d =-4 d =-3.5 d =-2 d =-0.5 d =0 d =0.5 < nd > V/4
FIG. 7: − d h ˆ n d i dV vs V /
4Γ for Γ = 0 . U = 8, and various ǫ d from Kondo to mixed valence regime. We see that the smallpeak nearby V = 0 only appears when ǫ d → − U . Inset: Thecorresponding h ˆ n d i vs V / ν s ( λ ) f s ( λ )+ ν h ( λ ) f h ( λ ). Under this assumption ν SBA ( λ )retains its form in and out-of-equilibrium and the generalexpression for h ˆ n d i is n d ( µ , µ ) = h Ψ , µ , µ | ˆ n d | Ψ , µ , µ i (39)= 2 Z ∞ B d λ σ b ( λ ) ν SBA ( λ ) + Z ∞ B d λ σ b ( λ ) ν SBA ( λ ) ! As the form for ν SBA ( λ ) is proved to be exact in equilib-rium, we shall regard Eq. (39) as an exact result for h ˆ n d i even out of equilibrium and valid in all different range of U , ǫ d , Γ under the assumption that the integrand doesnot change its form for in and out of equilibrium, whichis the case for general results of SBA. In the numeri- -4 -3 -2 -1 00.00.10.20.30.40.50.60.7 -4 -3 -2 -1 00.10.20.30.40.50.60.70.80.91.0 - < nd > / d d /4 V=0 V=1 V=1.5 V=2 V=3 V=4 < nd > d /4 V=0 V=1 V=1.5 V=2 V=3 V=4
FIG. 8: − ∆ h ˆ n d i ∆ ǫ d for various fixed voltages as a function of ǫ d for Γ = 0 . U = 8. Inset shows h ˆ n d i vs ǫ d for various fixedvoltage. cal results shown hereafter we shall use this expression,Eq.(39), for matrix element of dot occupation rather thanEq. (36). We adopt the same voltage drive scheme by fix-ing µ and lowering µ .By using this result we do not need to confine ourselffor large U Γ . The case for different U Γ with ǫ d = − U and for U = 8 , Γ = 0 .
25 with different ǫ d are shown inFig. 6 and Fig. 7. The main features of these plots are arelatively slow decrease of the dot occupation at low volt-age followed by an abrupt drop of h n d i . The decrease of h n d i takes place within a range of voltage of the orderof Γ. Then as we increase the voltage further anotherplateau develops. Note that, as expected, the bigger U isthe higher the voltage needed to drive the system out ofthe h n d i = 1 plateau. In a sense the charge fluctuationsare strongly frozen at large U and it costs more energyto excite them. The voltage where the abrupt drop in h n d i occurs corresponds to the energy scale at which the”charge fluctuation peak” was observed in the conduc-tance plots. This can be seen by comparing the positionof the broader peak in Fig. 5 with that of the abrupt dotoccupation drop in Fig. 7.Similar to the differential conductance we may definethe nonequilibrium charge susceptibility as χ c ( V ) | ǫ d = − ∂ h ˆ n d i ∂V that we obtain by taking a numerical derivative of the dotoccupation data with respect to the voltage. In the caseof U = − ǫ d / V ≃ U Γ . We identify this peak asa small remnant of the charge fluctuations in the Kondoregime. This statement is confirmed by noticing that thispeak goes away as U Γ increases, vanishing when U → ∞ U Andersonmodel is discussed. The second peak is located at thesame voltage as the charge fluctuation peak observed inthe conductance plots and is therefore associated to theresponse of the renormalized impurity level to the chargesusceptibility. This can be seen when comparing Fig. 5and Fig. 7.Another interesting quantity, the usual charge suscep-tibility , defined by χ c ( ǫ d ) | V = − ∂ h ˆ n d i ∂ǫ d , can also be quali-tatively described. In Fig. 8 we plot − ∆ h ˆ n d i ∆ ǫ d as a functionof ǫ d as we only have a few points in fixed ǫ d for finitevoltage. Notice that χ c ( ǫ d ) | V tends to be an universalcurve in large voltage, indicating charge on the dot re-mains at some constant value in the steady state withlarge voltage. This constant value at large voltage, aspointed out by C. J. Bolech, is around 0 .
65 for ǫ d = − U case. In preparing this article we noticed that a simi-lar computation, adopting the same asymmetric voltagedrive protocol as we have here, is carried out by R. V.Roermund et al for the dot occupation out of equilib-rium by using equation of motion method. We do get asimilar value for the dot occupation at large voltage. Thisvalue is different from the dot occupation value n d ≃ . U is turned off asshown in Fig. 13. This difference might have to do withthe 0 . inhigh temperature (temperature is high compared withthe Kondo scale but still small compared with phononmodes or electronic level) and zero magnetic field as thelinear response conductance given by n d = 0 .
65 by usingFriedel sum rule is around 0 .
73. In a sense the voltageseems to play a similar role to the temperature on theway it influences the dot occupation. Further connectionbetween these two behaviors could be clarified by com-puting the decoherence factor as in Ref. 21. This deco-herence factor is related to the dot correlation functionout of equilibrium which can be computed in three-leadsetup by using our approach. E. Comparison with other theoretical andexperimental results
In most of the other theoretical approaches the symmetric voltage drive ( µ = − µ ) is usually as-sumed to preserve particle-hole symmetry in symmetriccase ( ǫ d = − U ). It is thus difficult for us to make anydefinite comparison with other theoretical results. Thequalitative feature, as shown by the black curves in Fig. 9done by D. Matsumoto by using perturbation expan-sion in U at strong coupling fixed point, is similar to ourresults in the sense that the height of the charge fluctua-tion side peak and width are almost the same. The majordifferences are in the shape of Kondo peak and the posi-tion of the charge fluctuation side peak. A clear signatureof renormalized dot level ǫ d as hinted in renormalizationcomputation is clearly seen in our result. The shapeof Kondo resonance nearby zero voltage deviates from its quadratic behavior expected from Fermi liquid picture atsmaller voltage in our case as is expected for asymmetricvoltage drive . FIG. 9: Comparison of our theory with perturbation expan-sion in U done by D. Matsumoto on dI/dV (y-axis in unit of2 e /h ) vs V /U (x-axis). Our data (Blue, purple, and brownlines correspond to Γ U = 0 . , . , .
063 respectively. ∆shown in inset is Γ in our notation. EQ in the inset is con-ductance computed by equilibrium density of state which isnot relevant to our discussion here.) is shown as the mainfigure and Fig.8 in Ref 22 is shown in the inset. In Ref 22the voltage is driven symmetrically, i.e. µ = − µ , renderingthe factor of two difference in the voltage (i.e. VU = 0 . eVU = 1 in the inset. e = 1 in ourconvention.) in comparing our result with that in Ref 22. - - - FIG. 10: Comparison of theory with experiment of dI/dV (y-axis in unit of e /h ) vs V (x-axis in unit of mV ). Insetis the original data graph published in Ref. 6. The red dotsare given by our theory for U Γ = 8 with voltage rescaled tofit with original data in unit of mV . The value of differentialconductance (experiment data in black line) is rescaled from(0 . , .
3) to (0 ,
2) in unit of e h . We can also compare our results with experiments. Asshown in the inset of Fig. 10 is the dIdV vs V measuredin Co ion transistor by J. Park et al. . We rescaled the4differential conductance and superimposed our numeri-cal results on the data graph. The measurement wasdone by using an asymmetric drive of the voltage (bykeeping µ = 0 and changing µ to be larger or smallerthan zero) and thus there is an asymmetry in the differ-ential conductance as a function of voltage as illustratedin the data curve. In our numerics we only compute thescenario for µ = 0 and lowering µ (only for V >
V < V = 0 axis which illus-trates the case of µ = 0 and lowering µ . To comparewith the correct voltage setup on the V <
V > . To describe these type of transistors we shallstart with the Anderson-Holstein Hamiltonian. We arecurrently exploring the possibility of solving this modelby the Bethe Ansatz approach. III. INFINITE U ANDERSON MODEL
In the limit of U Γ → ∞ the finite U two-lead Ander-son impurity Hamiltonian becomes the two-lead infinite U Anderson model. The latter model is closely related,via the Schrieffer-Wolff transformation , to the notori-ous Kondo model, a model of spin coupled to a Fermiliquid bath. The reason for that is simple: since U → ∞ the charge fluctuations are essentially frozen out and onlythe spin fluctuations dominate the low-energy physics.The Hamiltonian is given byˆ H = X i =1 , Z d x ψ † iσ ( x )( − i∂ x ) ψ iσ ( x ) + ǫ d d † σ d σ + t i ( ψ † iσ (0) b † d σ + d † σ bψ iσ (0)) (40)Here the bosonic operator b is introduced to conserve b † b + P σ d † σ d σ = 1 and by applying the slave boson tech-nique we project out the phase space of double occupancyoccurring in finite U case. The corresponding Bethe mo-menta distribution function for the infinite U Andersonmodel is given by2 σ (Λ) = 1 π − Z B −∞ d Λ ′ K (Λ − Λ ′ ) σ (Λ ′ ) − Z B −∞ d Λ ′ K (Λ − Λ ′ ) σ (Λ ′ ) (41)with K (Λ) = π +(Λ − Λ ′ ) .Eq. (41) can be derived directly following the proce-dures in the finite U Anderson model. It can also bederived from the finite U result, Eq. (6), by taking the large U limit ( U ≫ ǫ d , U ≫ Γ): x ( λ ) U → − vuut λU + + q ( λU + ) + Γ U → − s λU + + | λU + | → −
12 (1 + 2 λU + . . . ) → − λU = Λ Uy ( λ ) U → s − ( λU + ) + (( λU + ) + Γ U ) / → vuut ( λU + )( − ( Γ U ) ( λU + ) ) / )2 (43) →
14 ( Γ U ) ! / + O ( U − ) ≃ Γ U with Λ ≡ − λU . Similar procedures as in Appendix C givethe matrix element ν SBA ∞ (Λ) for the dot occupation inthe infinite U Anderson model in equilibrium to be ν SBA ∞ (Λ) = 2Γ(Λ − ǫ d ) + (2Γ) . (44)In going to the out-of-equilibrium regime ( µ = µ ) wefollow the same phenomenological method as for the fi-nite U case. The result for the spin-fluctuation andcharge-fluctuation contributions to the dot occupationare given by ν s ∞ (Λ) = 1Γ − ǫ d − Λ p ( ǫ d − Λ) + 4Γ ! ν h ∞ (Λ) = 2Γ(Λ − ǫ d ) + (2Γ) . (45)We shall again check the consistency with the exact resultfor the dot occupation in equilibrium, namely h X σ d † σ d σ i = 4 Z BD dΛ σ b (Λ) ν SBA ∞ (Λ)= 4 Z BD dΛ σ b (Λ)( ν s ∞ (Λ) f ∞ s (Λ) + ν h ∞ (Λ) f ∞ h (Λ)) . Here D is related to the bandwidth and B is determinedby the equilibrium Fermi energy µ = µ = 0. f ∞ s (Λ)and f ∞ h (Λ) are expressed as f ∞ s (Λ) = T ∞ k /π (Λ − B ) + ( T ∞ k ) f ∞ h (Λ) = 2Γ(Λ − B − ǫ d ) + (2Γ) . Here the Kondo scale T ∞ k used in f s (Λ) takes the form T ∞ k = p | D | Γ π e − π | ǫd | Γ . U case are shown in Fig.11. Again wesee a nice match between our phenomenological approachand the exact result for | ǫ d Γ | 6 = 0 and some mismatch inthe mixed valence region | ǫ d Γ | ≃
0. This is consistent withthe results for finite U . -20 -10 0 10 200.00.20.40.60.81.0 TBA pSBA < n d > d / -20 -10 0 10 200.00.51.01.52.0 d / FSR pSBA d I/ d V ( e / h ) FIG. 11: Top: h ˆ n d i vs ǫ d Γ for exact TBA result and pSBA.Bottom: Linear response conductance dI/dV | V → vs ǫ d Γ forexact result (FSR) and pSBA in the infinite U Andersonmodel. D Γ = − U the com-parison nearby mixed valence region ( ǫ d ≃
0) is poorer.
The corresponding spin and charge fluctuation matrixelement for current, J s ∞ (Λ) and J h ∞ (Λ), are given by J s ∞ (Λ) = 1 − ǫ d − Λ p ( ǫ d − Λ) + 4Γ J h ∞ (Λ) = 2Γ (Λ − ǫ d ) + (2Γ) (46)The current expectation value is given by h ˆ I i = 2 e ~ Z B B d Λ σ (Λ)( J s ∞ (Λ) f ∞ s (Λ) + J h ∞ (Λ) f ∞ h (Λ)) -10 -8 -6 -4 -2 0 2 40.00.51.01.52.0 a r b un i t Js( ) Jh( )
FIG. 12: J s (Λ) and J h (Λ) vs Bethe momenta Λ (scaled by Γ)in infinite U Anderson model. ǫ d Γ = − x ( λ ). where B and B are related to µ and µ by minimizingcharge free energy FF = 2 Z B D dΛ σ (Λ)(Λ − µ ) + Z B D dΛ σ (Λ)(Λ − µ ) ! . Before we proceed to discuss the numerical results forcurrent vs voltage in this infinite U model let us look atthe structure of J s ∞ (Λ) and J h ∞ (Λ) as a function of Λ asshown in Fig. 12. Λ here represents the bare energy ofthe quasi-particle and plays the same role as x ( λ ) in thefinite U Anderson model. J s ∞ (Λ) alone would reproducethe main feature in the Friedel sum rule for ǫ d ≪
0. Inthis region the linear response conductance comes mainlyfrom the spin fluctuations. The upper plot of Fig. 12 fixes ǫ d and shows J s ∞ (Λ) vs Λ. We may also fix Λ = 0 (in thesense of choosing the equilibrium Fermi surface energy atΛ = 0) and plot J s ∞ ( ǫ d ) vs ǫ d . In this way we can see that J s ∞ ( ǫ d ) vs ǫ d reproduces the overall structure of the linearresponse conductance from the Kondo region ( ǫ d ≤
0) tothe mixed valence regime ( ǫ d ≃ δ p + + δ p − , contributing to J s ∞ (Λ), as thephase shift related to spin-fluctuation. J h ∞ (Λ) gives a Lorentz shape in bare energy scale Λ.This structure is akin to the charge fluctuation side peakwith peak position at energy scale around ǫ d as seen fromlower plot of Fig. 12. Thus we identify the phase shift δ p + + δ p − , contributing to J h ∞ (Λ), as the phase shift re-lated to charge-fluctuation. These structures also applyto the case of the finite U Anderson model.Now let us discuss the out of equilibrium numericalresults. The voltage is again driven asymmetrically byfixing µ ≃ µ . The exact dot occupationvs voltage for different ǫ d for infinite U and U = 0, ǫ d Γ = − d / =-6 (U=0) d / =-6 d / =-5 d / =-4 < n d > V/ FIG. 13: h ˆ n d i vs V Γ in infinite U Anderson model (for Red,Blue, and Purple dots corresponding to ǫ d Γ = − , − , −
4. TheBlack dots are U = 0 and ǫ d Γ = − D Γ = −
100 in this graph. the dot occupation decreases slowly at low voltage anddevelops an abrupt drop at a voltage scale correspondingto impurity level ǫ d . Also notice the apparent differencebetween the U = 0 plot (black dots) and the U → ∞ case(red dots) and for the same value of ǫ d Γ . For U → ∞ , thedot occupation at large voltage is around 0 .
65 for ǫ d Γ ≪ U casewhen U Γ is large (cf. Section II D). In contrast the non-interacting case ( U = 0) shows that h n d i → . phenomenological current vs voltage and the cor-responding differential conductance vs voltage are plot-ted in the top figure of Fig. 14. Again we see the zerobias anomaly and a broad charge fluctuation side peakin the differential conductance vs voltage. The scalingrelation of differential conductance vs voltage expectedin small voltage region can also be extracted by rescalingthe voltage by T ∞∗ k as shown in bottom figure of Fig. 14.Here T ∞∗ k is given by T ∞∗ k = p | D | Γ π e − π | ǫd | . Notice this T ∞∗ k differs from T ∞ k with a factor of twowithin the exponent. This factor of two difference rep-resents the difference in the curvature of the parabola asfunction of ǫ d (the logarithm of half width at half maxi-mum of the Kondo peak vs ǫ d shows parabolic curve asin inset of Fig. 6 for finite U case). This factor of two ra-tio bears even closer resemblance to the results shown inRef. 5. Note that in bottom figure of Fig. 14 the positionsof the side peak are different and show no universality inthat region. It shows universality for VT ∗ k ≤ d / =-6 d / =-5 d / =-4 d / =-3 d I/ d V ( e / h ) V/ d / =-6 d / =-5 d / =-4 d / =-3 I V/ V Tk * FIG. 14: Top: dIdV vs V Γ in infinite U Anderson model. Insetshows the I − V curves for these parameters. D Γ = −
100 in thisgraph. Bottom: dIdV vs VT ∗ k shows the scaling relation nearbyzero voltage for ǫ d Γ = − , − , − IV. CONCLUDING REMARKS
In this article we have explicitly computed the non-equilibrium transport properties in the Anderson modelfor all voltages using the Scattering Bethe Ansatz. In thecase of equilibrium we have also shown the equivalence oftraditional Bethe Ansatz and Scattering Bethe Ansatz byevaluating dot occupation in equilibrium. For the expres-sion of current we have introduced phenomenological dis-tribution functions to set the weight for spin-fluctuationand charge-fluctuation contributions to the current. Theresult shows correct scaling relation in Kondo regime aswell as satisfying the Friedel sum rule for linear responsefor large U Γ .Other interesting quantities, such as the nonequilib-rium charge susceptibility or the usual charge suscep-tibility, are computed numerically via exact expressionfor dot occupation as a function of voltage and impuritylevel. We believe this is the first report of an exact com-putation of the dot occupation out-of-equilibrium and itmay have interesting application in quantum computingas we understand more the dephasing mechanism. We7have also compared our results with perturbation calcu-lation and experimental measurement of nonlinear differ-ential conductance of a quantum dot.The major difficulty we encounter by using SBA comesfrom the single particle phase shift for complex momentawhich leads to a breakdown of steady state conditionwhen out of equilibrium. One possible issue resultingin this is the local discontinuity at odd channel s op ,the choice we made to enable us to construct a scat-tering state with fixed particles from lead 1 and lead2. It can be proved that without this choice we can-not write down fixed number of particles incoming fromeach lead in this Anderson impurity model and simi-larly for IRLM. The other issue in the study for Ander-son model is whether we shall include all possible boundstates in the ground state construction. From the math-ematical structure we shall choose 4 type of bound statesbut the results from charge susceptibility seems to sug-gest 2 type of bound states is the correct choice. To checkwhether this is in general correct we plan to come backto study the whole spectrum, which include bound statewhen Bethe energy higher than impurity level, of IRLMas this model bares structure similarity to the Andersonmodel described in this article. Following the SBA onIRLM there are lots of numerical approach and differentexact methods developed for this model and detailedcomparison for different approaches is desired for betterunderstanding its physics and scaling relation. By learn-ing how to deal with complex momenta in this model wemay also find the rule which may lead us to the exact expression for current in this Anderson impurity model. Acknowledgment
We are grateful to Kshitij Wagh, Andres Jerez, CarlosBolech, Pankaj Mehta, Avi Schiller, Kristian Haule, andPiers Coleman for many useful discussions and most par-ticularly to Chuck-Hou Yee for his important help withthe numerics and to Natan Andrei for numerous discus-sions and fruitful ideas. S. P. would also like to thankDaniel Ralph and Joshua Park for permission to use theirdata and discussion. G. P. acknowledges support fromthe Stichting voor Fundamenteel Onderzoek der Materie(FOM) in the Netherlands. This research was supportedin part by NSF grant DMR-0605941 and DoEd GAANNfellowship.
Appendix A: Discussion of 2 strings vs 4 strings
As we have discussed in the main text the boundedpair, formed by p ± ( λ ) = x ( λ ) ∓ iy ( λ ), can be formed byquasi-momenta from lead 1 or lead 2. We have shown theresults for two type of strings (bound states). Namely thestrings are formed by { ij } = { , } with i , j denotingincoming lead indices. In this section we discuss the caseof 4 type of strings and show thier corresponding numer- ical results in out of equilibrium regime (In equilibriumthe 2 strings and 4 strings give the same result for dotoccupation).The density distribution for the Bethe momenta (ra-pidities) is denoted by σ ij ( λ ) with { ij } = { , , , } indicating the incoming electrons from lead i and lead j .The σ ij ( λ ) is given by4 σ ij ( λ ) = − π d x ( λ )d λ − X i,j =1 , Z ∞ B ij d λ ′ K ( λ − λ ′ ) σ ij ( λ ′ )(A1)The factor of 4 indicates 4 type of possible configurationsand the constraint of exclusions in rapidities λ in solv-ing the quantum inverse scattering problem. The idea isthat in equilibrium four type of distributions are equallypossible for each bound state bare energy 2 x ( λ ). The B ij play the role of chemical potentials for the Bethe-Ansatzmomenta and are determined from the physical chemicalpotentials of the two leads, µ i , by minimizing the chargefree energy, F = X i ( E i − µ i N i ) = X i Z ∞ B ij d λ ( x ( λ ) − µ i ) σ ( i ) ( λ ) dλ with σ (1) ≡ σ + σ + σ the lead 1 particle densityand σ (2) ≡ σ + σ + σ the lead 2 particle density.In the case of µ > µ we have B < B = B < B for this finite U Anderson model but the equation for σ ij ( λ ) is the same for different combination of i and j .The reason is we put a quasi-hole state, rather than aquasi-particle, in the integral equation Eq.(A1) similarto the treatment of Wiener-Hopf approach. For example,for B < λ < B there could be three type of quasi-particle state { ij } = { , , } and we put { ij } = { } state as quasi-hole state. This hole state still count oneweight of the probability of 4 distributions and thereforethe factor of 4 on the left hand side of Eq.(A1) retainseven out of equilibrium. Similar idea is also applied intwo type of bound state (strings) solution.Other than their differences in the density distributionthe computations for the current and dot occupation ex-pectation value are quite similar to the two strings case.We show their numerical results in the following.The differential conductance vs voltage as shown inFig.15, obtained by taking numerical derivative on cur-rent vs voltage data, essentially gives the same picture asin two strings case, namely a sharp Kondo peak nearby V = 0 and a broad side peak corresponding to chargefluctuations. In the case of h n d i vs V , however, there isan additional feature occurring at an energy scale higherthan the energy scale of the charge fluctuation side peak(corresponding to the voltage position of 2nd peak shownin the inset) as shown in Fig. 16. This is especially ap-parent if we looked at the nonequilibrium charge suscep-tibility as shown in inset of Fig.16.As we do not expect there should be any further chargefluctuations, we rule out, by physical argument, the pos-sibility of 4 strings configuration.8 d =-4 d =-3 d =-2 d =-1 d = 0 d = 1 d I/ d V ( e / h ) V/4
FIG. 15: dIdV vs V for U = 8, Γ = 0 .
25 and various ǫ d from ǫ d = − U to ǫ d = 1. The inset is the enlarged region nearbyzero voltage. U=4,U=8,U=16, < n d > V/4
U=4 U=8 U=16
FIG. 16: h n d i vs V for different U , Γ = 0 .
25 and ǫ d = − U .The inset is − ∂ h n d i ∂V | ǫ d vs V voltage. A third peak shows upin U = 4 case. Appendix B: Two particles solution and choice of s op For the two particles solution we follow similar con-struction in P B Wiegmann and A M Tsvelick’s work and the Scattering Bethe Ansatz approach developed byP. Mehta and N. Andrei . Since Eq.(1) is rotational in-variant the spin quantum number is conserved. We showthe solution with both particles with spin singlet incom-ing from lead 1 as an example in the following. Spinquantum number in z direction S z is a good quantumnumber and we can write the two particle solution of S z = 0 state as: | Ψ i = n Z dx dx { Ag ( x , x ) ψ † e ↑ ( x ) ψ † e ↓ ( x )+ Ch ( x , x ) ψ † o ↑ ( x ) ψ † o ↓ ( x )+ Bj ( x , x )( ψ † e ↑ ( x ) ψ † o ↓ ( x ) − ψ † e ↓ ( x ) ψ † o ↑ ( x )) } + Z dx ( Ae ( x )( ψ † e ↑ ( x ) d †↓ − ψ † e ↓ ( x ) d †↑ )+ Bo ( x )( ψ † o ↑ ( x ) d †↓ − ψ † o ↓ ( x ) d †↑ )) + Amd †↑ d †↓ o | i Here
A, B, C are arbitrary constants to be determinedlater. To satisfy ˆ H | Ψ i = E | Ψ i = ( k + p ) | Ψ i we have:0 = [ − i ( ∂ x + ∂ x ) − E ] g ( x , x )+ t [ δ ( x ) e ( x ) + δ ( x ) e ( x )] (B1)0 = [ − i ( ∂ x + ∂ x ) − E ] h ( x , x ) (B2)0 = [ − i ( ∂ x + ∂ x ) − E ] j ( x , x ) + tδ ( x ) o ( x )(B3)0 = ( − i∂ x − E + ǫ d ) e ( x ) + tg (0 , x ) + tδ ( x ) m (B4)0 = ( − i∂ x − E + ǫ d ) o ( x ) + tj (0 , x ) (B5)0 = ( U + 2 ǫ d ) m + 2 te (0) − Em (B6)For U = 0 the model becomes non-interacting and thetwo particles solution becomes direct product of two oneparticle solutions. | Ψ i = | ψ k ↑ i ⊗ | ψ p ↓ i + | ψ p ↑ i ⊗ | ψ k ↓ i = Z dx dx { ( g k ( x ) ψ † e ↑ ( x )+ h k ( x ) ψ † o ↑ ( x )+ e k d †↑ δ ( x ))( g p ( x ) ψ † e ↓ ( x ) + h p ( x ) ψ † o ↓ ( x ) + e p d †↓ δ ( x ))+ ( g p ( x ) ψ † e ↑ ( x ) + h p ( x ) ψ † o ↑ ( x ) + e p d †↑ δ ( x ))( g k ( x ) ψ † e ↓ ( x ) + h k ( x ) ψ † o ↓ ( x ) + e k d †↓ δ ( x )) }| i Therefore at U = 0 we have: g ( x , x ) = g k ( x ) g p ( x ) + g k ( x ) g p ( x ) h ( x , x ) = h k ( x ) h p ( x ) + h k ( x ) h p ( x ) j ( x , x ) = g k ( x ) h p ( x ) + h k ( x ) g p ( x ) e ( x ) = e k g p ( x ) + e p g k ( x ) o ( x ) = e k h p ( x ) + e p h k ( x ) m = 2 e p e k Now for U = 0 we shall derive the solution of this form g ( x , x ) = Z kp ( x − x ) g k ( x ) g p ( x )+ Z kp ( x − x ) g k ( x ) g p ( x ) (B7)Plug Eq.(B7) into Eq.(B1) we get e ( x ) = Z kp ( − x ) g p ( x ) e k + Z kp ( x ) g k ( x ) e p (B8)9Plugging above two results into Eq.(B4) into Eq.(B6) weget for m = 2 ˜ Z kp (0) e k e p we have:( − i∂ x Z kp ( − x )) g p ( x ) e k + ( − i∂ x Z kp ( x )) g k ( x ) e p − tZ kp ( − x ) e p δ ( x ) e k − tZ kp ( x ) e k δ ( x ) e p + 2 t ˜ Z kp (0) e k e p = 0 (B9)2 ˜ Z kp (0) e k e p = 2 t ( Z kp (0) g p (0) e k + Z kp (0) g k (0) e p ) p + k − U − ǫ d (B10)Now take Z kp ( x ) = e − iφ kp θ ( − x ) + e iφ kp θ ( x )we get tan( φ kp ) = − Ut ( k − p )( p + k − U − ǫ d ) and˜ Z kp (0) = k + p − ǫ d k + p − U − ǫ d Z kp (0). Define Γ ≡ t and B ( k ) ≡ k ( k − ǫ d − U ) as in Ref. 34 we can rewritetan( φ kp ) = − U Γ( B ( k ) − B ( p )) .From Eq.(B2) we can write h ( x , x ) as: h ( x , x ) = Z ookp ( x − x ) h k ( x ) h p ( x )+ Z ookp ( x − x ) h k ( x ) h p ( x ) (B11)with arbitrary Z ookp ( x − x ). Now write j ( x , x ) as: j ( x , x ) = Z eokp ( x − x ) g k ( x ) h p ( x )+ Z eokp ( x − x ) h k ( x ) g p ( x ) (B12)again with Z eokp ( x − x ) undetermined. Plug Eq.(B12)into Eq.(B3) we get o ( x ) is written as: o ( x ) = Z eokp ( − x ) h p ( x ) e k + Z eokp ( x ) h k ( x ) e p (B13)Now if we choose Z eokp ( x − x ) = Z kp ( x − x ) and plugEq.(B12) and Eq.(B13) into Eq.(B5) we get:( − k + ǫ d ) Z kp ( − x ) h p ( x ) e k + ( − p + ǫ d ) Z kp ( x ) h k ( x ) e p + t ( Z kp ( − x ) h p ( x ) g k (0) + Z kp ( x ) h k ( x ) g p (0))+( − i )( ∂ x Z kp ( − x )) h p ( x ) e k + ( − i )( ∂ x Z kp ( x )) h k ( x ) e p = − φ kp )( h p (0) e k − h k (0) e p ) = 0 (B14)To satisfy Eq.(B14) we can set h p (0) = 0 for arbitrary p .This can be done by choosing s op = − Z ookp ( x − x ) is arbitrary we can choose Z ookp ( x − x ) = Z kp ( x − x ). Also from Eq.(B10) we have˜ Z kp (0) = p + k − ǫ d p + k − U − ǫ d Z kp (0) (B15)Since the Hamiltonian in Eq.(1) has rotational invariancethe general form of scattering matrix for particles withmomentum k, p and spins a , a is given by: S a ′ a ′ a a ( k, p ) = b ( k, p ) + c ( k, p ) ˆ P a ′ a ′ a a (B16)where ˆ P a ′ a ′ a a = (1 a ′ a · a ′ a + ~σ a ′ a · ~σ a ′ a ) is the permutationoperator in spins. For antiparallel spins (singlet state) as shown above ˆ P a ′ a ′ a a = − b ( k, p ) − c ( k, p ) = Z kp ( x > Z kp ( x < B ( k ) − B ( p ) − i U Γ B ( k ) − B ( p ) + i U Γ (B17)For the triplet state ( ˆ P a ′ a ′ a a = 1) the interaction termwith the impurity is absent and the particles passingthrough each other without changing their phase b ( k, p ) + c ( k, p ) = 1 (B18)Thus from Eq.(B17) and Eq.(B18) we get the two particleS-matrix as:ˆ S ( k, p ) a ′ i a ′ j a i a j = ( B ( k ) − B ( p )) I a ′ i a ′ j a i a j + i U Γ P a ′ i a ′ j a i a j B ( k ) − B ( p ) + i U Γ (B19)Thus the integrability of two lead with Anderson typedot system is the similar to the integrability of one leadAnderson model.The choice of identical two particles S-matrices (bychoosing s op = −
4) enables us to construct the scat-tering state labeled by lead indices by choosing appro-priate
A, B, C in this even-odd basis. For example, ifboth particles are coming from lead 1, we shall choose(
A, B, C ) = A ( t t , − t t t , t t ) such that the amplitude ofincoming state from lead 2 is zero ( A being an overallrenormalization constant). We can therefore label theeigenstate by the incoming state from lead i and/or lead j . Without this s op term we cannot write back from even-odd basis to lead indices basis in this two leads Andersonmodel and similarly in IRLM in Ref. 7. Appendix C: Equivalence of TBA and SBA inequilibrium
Eq.(19) can be proved to be exact by comparingwith the traditional Bethe Ansatz where h P σ d † σ d σ i =2 R ∞ B d λ σ imp ( λ ) with impurity density σ imp ( λ ) given by σ imp ( λ ) = δ p + + δ p − π − Z ∞ B d λ ′ K ( λ − λ ′ ) σ imp ( λ ′ ) (C1)The driving term (first term) of Eq.(C1) is expressed bybare phase shift δ p + + δ p − and thus we can view σ imp ( λ ) asthe dressed phase shift across the impurity. By compar-ing Eq.(C1) and Eq.(6) in equilibrium ( σ i ( λ ) = σ b ( λ ) de-scribing bulk quasi-particle density when B = B = B .)we get Z ∞ B d λ σ imp ( λ ) (cid:18) − π d x ( λ )d λ (cid:19) =2 Z ∞ B d λ σ b ( λ ) (cid:18) δ p + + δ p − π (cid:19) (C2)0by noting that the integration kernel K ( λ − λ ′ ) is sym-metric in λ and λ ′ . Since the equality is true for arbitrary B we can also rewrite Eq.(C2) as Z ∞ B d λ σ imp ( λ ) = 2 Z ∞ B d λ σ b ( λ ) δ p + + δ p − − d x ( λ )d λ ! ≡ Z ∞ B d λ σ b ( λ ) ν T BA ( λ )and the resulting ν T BA ( λ ) is given by ν T BA ( λ ) = − ˜ x ( λ ) y ′ ( λ ) x ′ ( λ ) − ˜ y − ( λ )˜ x ( λ ) + ˜ y ( λ )+ ˜ x ( λ ) y ′ ( λ ) x ′ ( λ ) + ˜ y + ( λ )˜ x ( λ ) + ˜ y ( λ ) (C3)Now let us show the computation for ν SBA ( λ ). First wewrite one particle state of Eq.(1) in even channel (with s ek = 0 for the moment) as | k, σ i = Z e ikx α † ek,σ ( x ) dx | i (C4)= Z e ikx { (¯ θ + A k θ ) ψ † eσ + B k d † σ δ ( x ) } dx | i Solving ˆ H | k, σ i = k | k, σ i we get − i ( − A k ) + B k t = 0 ǫ d B k + t A k kB k Thus we get A k = k − ǫ d − i t k − ǫ d + i t and B k = tk − ǫ d + i t . We mayalso define g k ( x ) = e ipx (¯ θ + A k θ ) and e k = B k to haveeasier comparison with Wiegmann and Tsvelick’s work .The two particles state is obtained by constructing prod-uct of two α † ep,σ ( x ) particles state with appropriate twoparticles S-matrix expressed in Z k + k − ( x − x ).In principle we shall use | Ψ , N , N i as the many bodystate to compute expectation value. However the simpli-fication here, similar to the case of IRLM in Ref.7, is thatdifferent λ (corresponding to different p ( λ )) are orthogo-nal to each other in L → ∞ limit. Thus the many bodyexpectation value can be obtained via two body compu-tation and the rest just get canceled by normalizationfactor. To put it more explicitly let us denote p i as thereal part of the complex pair p ± i . Different | p i i is orthog-onal to each other under the condition of size of the leadstaken to infinity, or h p i | p j ih p i | p i i → L → ∞ for i = j . Thusthe evaluation of matrix element for operator ˆ o is givenby h p , p , . . . | ˆ o | p , p , . . . ih p , p , . . . | p , p , . . . i = X p i h p i | ˆ o | p i ih p i | p i i . Based on this result we demonstrate the explicit compu-tation for dot occupation by two particles wavefunctionsin the following.Denote | Ψ i as the two particles solution. We may writespin singlet state as | Ψ i = Z dx dx A n e i ( kx + px ) Z kp ( x − x ) α † ek, ↑ ( x ) α † ep, ↓ ( x ) o | i = Z dx dx n Z kp ( x − x ) { g k ( x ) g p ( x ) ψ †↑ ( x ) ψ † e ↓ ( x ) + g k ( x ) e p ψ †↑ ( x ) d †↓ δ ( x )+ e k g p ( x ) d †↑ δ ( x ) ψ †↓ ( x ) + e k e p d †↑ d †↓ δ ( x ) δ ( x ) } − Z kp ( x − x ) { g k ( x ) g p ( x ) ψ † e ↓ ( x ) ψ † e ↑ ( x )+ g k ( x ) e p ψ † e ↓ ( x ) d †↑ δ ( x ) + e k g p ( x ) d †↓ δ ( x ) ψ † e ↑ ( x ) + e k e p d †↓ d †↑ δ ( x ) δ ( x ) } o | i = n Z dx dx [ Z kp ( x − x ) g k ( x ) g p ( x ) + Z kp ( x − x ) g k ( x ) g p ( x )] ψ † e ↑ ( x ) ψ † e ↓ ( x )+ Z dx [ Z kp ( x ) g k ( x ) e p + Z kp ( − x ) g p ( x ) e k ]( ψ † e ↑ ( x ) d †↓ − ψ † e ↓ ( x ) d †↑ ) + 2 e k e p ˜ Z kp (0) d †↓ d †↑ o | i With A denoting anti-symmetrization and ˜ Z kp (0) = k + p − ǫ d k + p − U − ǫ d Z kp (0). Solving ˆ H | k, σ ; p, − σ i = ( k + p ) | k, σ ; p, − σ i we obtain Z kp ( x − x ) = θ ( x − x ) + ( k − p )( k + p − ǫ d − U ) − iU t ( k − p )( k + p − ǫ d − U ) + iU t θ ( x − x )1For the case of bound state the two particle S-matrix isgiven by Z k + k − ( x − x ) = θ ( x − x ) ≡ θ x . The nor- malization factor and matrix element of dot occupationgiven by the even channel two particles wavefunction are h Ψ | Ψ i = Z dy dy Z dx dx ( θ y g k + ( y ) g k − ( y ) + θ y g k + ( y ) g k − ( y )) ∗ × ( θ x g k + ( x ) g k − ( x ) + θ x g k + ( x ) g k − ( x )) δ ( x − y ) δ ( x − y )+2 Z dy Z dx [ θ ( y ) g k + ( y ) e k − + θ ( − y ) g k − ( y ) e k + ] ∗ [ θ ( x ) g k + ( x ) e k − + θ ( − x ) g k − ( x ) e k + ] δ ( x − y )+4( e k + e k − ˜ Z k + k − (0)) ∗ ( e k + e k − ˜ Z k + k − (0)) X σ h Ψ | ˆ d † σ ˆ d σ | Ψ i = 2 Z dy Z dx [ θ ( y ) g k + ( y ) e k − + θ ( − y ) g k − ( y ) e k + ] ∗ [ θ ( x ) g k + ( x ) e k − + θ ( − x ) g k − ( x ) e k + ] δ ( x − y )+8( e k + e k − ˜ Z k + k − (0)) ∗ ( e k + e k − ˜ Z k + k − (0))= 2 (cid:26)Z dx [ θ ( x ) | g k + ( x ) e k − | + θ ( − x ) | g k − ( x ) e k + | ] + 4 | e k + e k − ˜ Z k + k − (0) | (cid:27) Note that the even channel bound state can be written assum over bound state of { , , , } (4 strings type)or { , } (2 strings type) with the same real part ofenergy k = x ( λ ). This can be viewed as the consistency counting from Fock basis to Bethe basis as electrons inlead 1 and lead 2 has 4 fold degeneracies in its initialstate (2 different spins in each lead). Also note that Z dx dx θ x | g k + ( x ) g k − ( x ) | = Z dx dx | e i ( k + x + k − x ) (¯ θ + θ A k + )(¯ θ + θ A k − ) | θ = Z dx dx e − ξ k ( x − x ) | ¯ θ ¯ θ θ + θ ¯ θ θ A k + + θ θ θ A k + A k − | = (cid:18) L ξ k − − e − ξ k L (2 ξ k ) (cid:19) (cid:0) | A k + A k − | (cid:1) + (cid:18) − e − ξ k L ξ k (cid:19) | A k + | Z dx θ ( x ) | g k + ( x ) e k − | = Z dx θ ( x ) | e i ( k + iξ k ) x ( θ ( − x ) + A k + θ ( x )) e k − | = Z L dx e − ξ k x | A k + e k − | = 12 ξ k (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) k − ǫ d + iξ k − i Γ k − ǫ d + iξ k + i Γ tk − ǫ d − iξ k + i Γ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 12 ξ k (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) tk − ǫ d + iξ k + i Γ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Z dx θ ( − x ) | g k − ( x ) e k + | = Z dx θ ( − x ) | e i ( k − iξ k ) x ( θ ( − x ) + A k − θ ( x )) e k + | = Z − L dx e ξ k x | A k − e k + | = 12 ξ k (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) tk − ǫ d + iξ k + i Γ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) with ˜ Z k + k − (0) = k − ǫ d )2( k − ǫ d ) − U Z k + k − (0) and Z k + k − (0) = based on our regularization scheme. By expressing k = x ( λ ) and ξ k = y ( λ ) and taking L → ∞ thus preserving L terms only we get2 h Ψ | P σ ˆ d † σ ˆ d σ | Ψ ih Ψ | Ψ i = 1 L ν
SBA ( λ ) (C5)= 1 L ( x ( λ ) + ˜ y ( λ ) + 16 y ( λ )Γ (˜ x ( λ ) + ˜ y − ( λ ))(˜ x ( λ ) + ˜ y ( λ )) (cid:18) ˜ x ( λ )2˜ x ( λ ) − U (cid:19) ) . By expressing ν T BA ( λ ) and ν SBA ( λ ) in λ explicitly wesee that ν T BA ( λ ) = ν SBA ( λ ). Since h P σ d † σ d σ i =2 R ∞ B dλσ imp ( λ ) in TBA we have proved that the expecta- tion value evaluated by the state we constructed is exactand the equivalence of SBA and TBA in equilibrium inthis two-lead Anderson model. D. Goldhaber-Gordon, H. Strickman, D. Mahalu, D.Abusch-Magder, U. Meirav and M. A. Kastner, Nature , 156 (1998) M. Grobis, I. G. Rau, R. M. Potok, H.Shtrikman, and D. Goldhaber-Gordon, Phys. Rev. Lett. , 246601 (2008) J. Schmid, J. Weis, K. Eberl, and K. Von Klitzing, Physica
B258 , 182 (1998) S. M. Cronenwett, T. H. Oosterkamp and L. P. Kouwen-hoven, Science , 540 (1998) S. M. Cronenwett, H. J. Lynch, D. Goldhaber-Gordon,L. P. Kouwenhoven, C. M. Marcus, K. Hirose, N. S.Wingreen, and V. Umansky, Phys. Rev. Lett. , 226805(2002) W. G. van der Wiel, S. De Franceschi, T. Fujisawa, J. M.Elzerman, S. Tarucha, L. P. Kouwenhoven, Science J. Park, A. N. Pasupathy, J. I. Goldsmith, C. Chang,Y. Yaish, J. R. Petta, M. Rinkoski, J. P. Sethna, H. D.Abruna, P. L. McEuen, and D. C. Ralph, Nature , 722(2002) P. Mehta and N. Andrei, Phys. Rev. Lett. ,216802 (2006) ibid. , 086804 (2008). See also cond-mat/0702612 for more detailed discussion of techniques ofScattering Bethe Ansatz and cond-mat/0703426. T. K. Ng and P. A. Lee, Phys. Rev. Lett. , 1768 (1988) Glazman, L. I. and Raikh, M. E., JETP Letters, , 452(1988) Y. Meir and N. Wingreen, Phys. Rev.
B 49 , 11040 (1994) M. Pustilnik and L. I. Glazman, Phys. Rev. Lett. ,216601 (2001) M H. Hettler and H. Schoeller, Phys. Rev. Lett. , 4907(1995) A. Oguri, Phys. Rev.
B 64 , 153305 (2001) A. Oguri, J. Phys. Soc. Jpn. , 110 (2005) Z. Ratiani and A. Mitra, Phys. Rev.
B 79 , 245111 (2009) A. Rosch, J. Kroha, and P. W¨olfle, Phys. Rev. Lett. ,156802 (2001) K. S. Thygesen1 and A. Rubio, Phys. Rev.
B 77 , 115333(2008) L. G. G. V. Dias da Silva1, F. Heidrich-Meisner, A. E.Feiguin, C. A. Busser, G. B. Martins, E. V. Anda, and E.Dagotto, Phys. Rev.
B 78 , 195317 (2008) J. Eckel, F. Heidrich-Meisner, S. G. Jakobs, M. Thorwart,M. Pletyukhov and R. Egger, cond-mat/1001.3773 (2010) F. Heidrich-Meisner, A. E. Feiguin, and E. Dagotto, Phys. Rev.
B 79 , 235336 (2009) R. V. Roermund, S-Y Shiau and M. Lavagna, cond-mat/1001.3873 (2010) D. Matsumoto, J. Phys. Soc. Jpn. , 1449 (2000) C. D. Spataru, M. S. Hybertsen, S. G. Louie, and A. J.Millis , Phys. Rev.
B 79 , 155110 (2009) A. Schiller and N. Andrei, cond-mat/0710.0249 (2007). E. Boulat and H. Saleur, Phys. Rev.
B 77 , 033409 (2008) E. Boulat, H. Saleur and P. Schmitteckert, Phys. Rev. Lett. , 140601 (2008) J. S. Langer and V. Ambegaokar, Phys. Rev. , 1090(1961) D. C. Langreth, Phys. Rev. , 516 (1966) R. M. Konik, H. Saleur, and A. Ludwig, Phys. Rev. Lett. , 236801 (2001) R. M. Konik, H. Saleur, and A. Ludwig, Phys. Rev.
B 66 ,125304 (2002) P. Fendley, A.W.W. Ludwig, H. Saleur, Phys. Rev.
B 52 ,8934 (1995) P. Fendley, A.W.W. Ludwig, H. Saleur, Phys. Rev. Lett. , 3005 (1995) P. B. Wiegman and A. M. Tsvelik, J. Phys.
C. 16 , 2281(1983); P. B. Wiegman and A. M. Tsvelik, Adv. in Phys. , 453 (1983) N. Kawakami and A. Okiji, J. Phys. Soc. Jap. 51, 1145(1982); N. Kawakami and A. Okiji, Solid State Commun.43, 365 (1982) P. Schlottmann, Z. Phys.
B 52 , 127 (1983) The thermodynamic Bethe Ansatz proof of ground stateconfiguration for two leads case was done by C. J. Bolechand then by S. P. Chao. The idea is to write down finitetemperature free energy for two leads system at differentchemical potentials and find the lowest energy state as tem-perature is taken to zero. The ground state is shown tobe formed by complex solutions originated from poles (ze-ros) of two particles S-matrices. For details see S. P. Chao,Ph.D. thesis (2010). A renormalizable Hamiltonian such as the Anderson modelrequires regularization and a cut-off scheme to define it.The results are universal once the cut-off is removed. Inintermediate stages as the cut off is finite it is importantto adopt a scheme that does not break integrability. Thescheme adopted here satisfies this requirement. N.Andrei,K. Furuya and J. H. Lowenstein, Rev. Mod. Phys. ,331 (1983), section VI. In our regularization scheme the locally discontinuous function s ( x ) ≡ θ ( x ) θ ( − x ) satisfies ∂ x s ( x ) = 0. Without the regularization factor in the odd sector themodel is still integrable. However, it is not possible to iden-tify the fixed number of incoming particles by this choice,A. Nishino and N. Hatano in J. Phys. Soc. Jpn. , 063002(2007). It’s also possible to construct scattering eigenstate not inthe Bethe Ansatz form. See J. T. Shen and S. Fan, Phys.Rev. Lett. , 153003 (2007) and A. Nishino, T. Imamuraand N. Hatano, Phys. Rev. Lett. , 146803 (2009) forIRLM, and the two particles state for Anderson model inT. Imamura, A. Nishino and N. Hatano, Phys. Rev. B 80 ,245323 (2009). It is, however, difficult to find a consistentway to write down N particles eigenstate through this ap-proach. N. Kawakami and A. Okiji, Phys. Rev.
B 42 , 2383 (1990) D. K.K. Lee and P. A. Lee, Physica
B 259-261 , 481(1999) This dressed energy is the sum of the dressed energy ofspinon and that of antiholon nearby equilibrium Fermi sur-face. N. Andrei, Phys. Lett.
A 87 , 299 (1982). A. O. Gogolin, R.M. Konik, A.W.W. Ludwig, and H.Saleur, Ann. Phys. (Leipzig) , 678 (2007). The equality can be proved analytically in linear response.In numerics we also see very good agreement, especially forlarge voltage, with voltage computed by µ − µ given byfree energy and voltage computed by difference in particlenumber. F. D. M. Haldane, Phys. Rev. Lett. , 416 (1978) A.C. Hewson, A. Oguri, and D. Meyer, Eur. Phys. J.
B40 , 177 (2004) J. Paaske and K. Flensberg, Phys. Rev. Lett. , 176801(2005) A. C. Hewson,
The Kondo Problem to Heavy Fermions ,Cambridge Studies in Magnetism (1993). Here we adopted the Kondo scale as in the article by N.S. Wingreen and Y. Meir, Phys. Rev.
B 49 , 11040 (1994).We put a factor of √
10 to increase this scale so the Kondopeak can be observed more easily in numerics. E. Lebanon and A. Schiller, Phys. Rev.