Failure of mean-field approach in out-of-equilibrium Anderson model
Bertalan Horváth, Bence Lazarovits, Olivier Sauret, Gergely Zaránd
aa r X i v : . [ c ond - m a t . s t r- e l ] D ec Failure of mean-field approach in out-of-equilibrium Anderson model
B. Horv´ath , B. Lazarovits , , O. Sauret , G. Zar´and Theoretical Physics Department, Institute of Physics,Budapest University of Technology and Economics, Budafoki ´ut 8, H-1521 Hungary Research Institute for Solid State Physics and Optics of the Hungarian Academy of Sciences,Konkoly-Thege M. ´ut 29-33., H-1121 Budapest, Hungary
To explore the limitations of the mean field approximation, frequently used in ab initio molecularelectronics calculations, we study an out-of-equilibrium Anderson impurity model in a scatteringformalism. We find regions in the parameter space where both magnetic and non-magnetic solutionsare stable. We also observe a hysteresis in the non-equilibrium magnetization and current as afunction of the applied bias voltage. The mean field method also predicts incorrectly local momentformation for large biases and a spin polarized current, and unphysical kinks appear in variousphysical quantities. The mean field approximation thus fails in every region where it predicts localmoment formation.
PACS numbers: 73.63.Kv, 75.20.Hr, 71.23.An, 73.23.-b
The Anderson impurity model (AIM) has been thesubject of great theoretical and experimental interest inthe past decades (for a review see Ref. 2). There is anumber of experimental systems including quantum dots,or single atoms and molecules contacted by leads, whichprovide experimental realizations of various versions ofthe AIM under out-of-equilibrium conditions. These sys-tems are not just prototypes of out of equilibrium sys-tems but a theoretical understanding of them would becrucial for future molecular electronics and mesoscopicapplications.Anderson constructed his famous model in Ref. 1 to de-scribe local moment formation and solved it within themean-field (MF) approximation. Within this approxi-mation, he found a phase transition to a state wheremagnetic moments are formed. Further work revealedthat, in reality, quantum fluctuations of this local mo-ment lead to the formation of a Kondo-singlet betweenthe impurity and conduction electrons at low tempera-ture, where the impurity spin is thus completely screened.The spontaneous symmetry breaking found by Andersonis thus an artifact of the mean field approximation. Nev-ertheless, the MF treatment indicates clearly the regionsof strong correlations, and it can also serve as a startingpoint for accurate approximations as in the local momentapproach (LMA) or interpolative perturbation theory (IPT). The latter approach can easily be generalized tonon-equilibrium situations using Keldysh formalism .In lack of more accurate methods, the mean field ap-proximation is also used in molecular electronics calcula-tions, where LDA or eventually Hartree-Fock equationsare solved in a scattering state or Keldysh approach todescribe moment formation . However, it is not clear atall, how reliable these approximations are. The purposeof this paper is to shed some light on the weaknesses ofthe non-equilibrium mean field approach on the simplestpossible test case, the out of equilibrium Anderson model,and to show where usual ab initio calculations should fail.Our conclusion is that the mean field approach fails qual-itatively and quantitatively essentially everywhere where it predicts local moment formation. Our study, which isbased on the scattering state formalism, is complemen-tary to the recent work of Komnik and Gogolin , whoused a Green’s function formalism to study the mean fieldequations of the non-equilibrium Anderson model. As weshall see, in the strongly correlated regions several arti-facts emerge such as non-equilibrium driven spontaneoussymmetry breaking as well as multiple stable solutionswhich lead to the appearance of hysteresis. These in-stabilities are probably also parts of the reasons, whynon-equilibrium IPT suffers from all kinds of instabili-ties. These instabilities were avoided in previous worksby applying spin-independent approximations or usingan interpolative self-energy or both .The non-equilibrium AIM Hamiltonian consists of fourparts. The first part describes a single impurity level withenergy ε d and an on-site Coulomb interaction ( U ) H d = X σ = ↑ , ↓ ε d d † σ d σ + U n ↑ n ↓ , (1)where d † σ and d σ are the creation and annihilation op-erators of the impurity electrons corresponding to spinstate σ and n σ = d † σ d σ . The second and third terms de-scribe the left (L) and right (R) leads which we model bytight-binding chains, H α = X k,σ (cid:0) − t cos k + µ α (cid:1) c † kασ c kασ . (2)Here α ∈ ( L, R ), c † kασ and c kασ are the creation andannihilation operators of a conduction electron of wavenumber k ∈ { , π } in lead α , ˜ t is the hopping along theleads, and µ α is the chemical potential of the left or theright lead. The chemical potentials of the two leads aredifferent due to a finite bias voltage leading to a non-equilibrium situation. The fourth part, H t , describes thetunneling between the leads and the impurity H t = V X σ (cid:2) d † σ ( V − c − σ + V + c σ ) + h.c. ] . (3)Here V ∓ are the hybridization matrix elements betweenthe impurity and the left and the right leads, and c l = ± σ denote the conduction electron annihilation operators onthe sites next to the impurity; sites along the left andright chains are labelled by l = {−∞ , . . . , − } and l = { , . . . , ∞} , respectively. For the sake of simplicity, herewe study a symmetrical situation, V − = V + = V , but ourconclusions are rather independent of this assumption.To study the Hamiltonian above we used a mean-fieldapproximation and replaced the impurity term as H d → H MFd = X σ ( ε d + U h n − σ i ) d † σ d σ . (4)The non-equilibrium expectation value of the occupationnumbers, h n σ i , in Eq. (4) can be obtained by solvingself-consistent equations discussed later. The expression ε d + U h n − σ i can be regarded as an effective energy level ofspin- σ electron. The hybridization between impurity andconduction electrons gives rise to a finite lifetime for im-purity states, reflected in the broadening of the effectiveimpurity levels with a finite width, Γ = 2 πV ρ = V / ˜ t where ρ is the density of states (DOS) of the conductionelectrons at the Fermi-level of the half-filled leads.To evaluate the non-equilibrium expectation values, h n σ i , we shall use a scattering formalism. The annihila-tion operator of the left-coming scattering state of energy ε and spin σ can be expressed as c σ ( ε ) = X l< c lσ e ikR l + α σ ( ε ) X l< c lσ e − ikR l ++ β σ ( ε ) X l> c lσ e ik ′ R l + γ σ ( ε ) d σ . (5)Here c lσ is the annihilation operator of the l th site and α σ ( ε ), β σ ( ε ) and γ σ ( ε ) are the reflection, transmissionand dot coefficients, respectively, which we obtain bysolving the corresponding Schr¨odinger equation. Sincethe energy is conserved in course of the scattering pro-cess, the wave numbers k and k ′ are connected by thedispersion relation ε = − t cos k + µ L = − t cos k ′ + µ R .Waves coming from the right hand side can be con-structed in a similar way. The occupation numbers arethen calculated from the dot coefficients, γ σ ( ε ) and γ ′ σ ( ε )corresponding to the states coming from the left and theright, respectively, h n σ i = ρ Z ∞−∞ d ε (cid:0) | γ σ ( ε ) | f L ( ε ) + | γ ′ σ ( ε ) | f R ( ε ) (cid:1) , (6)where the DOS of leads is approximated by its value atFermi-level, and f α ( ε ) ≡ f ( ε − µ α ) is the Fermi function.In the large bandwidth approximation (˜ t → ∞ , Γ = V / ˜ t finite) the occupation numbers become h n σ i = Γ2 π ∞ Z −∞ f L ( ε ) + f R ( ε )( ε − ε d − U h n − σ i ) + Γ d ε , (7) Γ / U µ / U increasing bias ( µ c2 )decreasing bias ( µ c1 )stability eq. µ / U m MC P µ c1 µ c2 FIG. 1: Magnetic (M), coexistence (C) and paramagnetic (P)regions as a function of µ/U and Γ /U values for ε d /U = − / /U = 0 .
05 (along the dotted linein the main figure) which simplifies further for zero temperature as h n σ i = X α ∈ ( L,R ) π cot − (cid:18) ε d + U h n − σ i − µ α Γ (cid:19) . (8)These self-consistent equations can also be obtained us-ing the Keldysh formalism . Eqs. (7) or (8) are solvediteratively by assuming that the applied bias is symmet-rical µ L = µ/ µ R = − µ/
2. In these calculations thebias voltage is increased or decreased gradually and thelocal stability of the solutions is always checked.First, let us discuss the bias-dependence of the occupa-tion numbers at T = 0. In the inset of Fig. 1 the magneti-zation m = h n ↑ i − h n ↓ i is plotted both for increasing anddecreasing bias voltage at a fixed ratio of Γ /U = 0 . ε d /U = − . µ > .
5) the stable solution isparamagnetic. Between the two limiting cases a regionappears, µ c < µ < µ c , where both the magnetic andthe non-magnetic solutions are stable. We shall refer tothis region as a coexistence region. The existence ofmultiple stable solutions for the occupation numbers inthis region is reflected in the hysteresis of the magneti-zation too. The sharp decay of the magnetization shownin the inset of Fig. 1 and the existence of a hysteresisbetween the critical fields indicate clearly a first ordertransition, predicted incorrectly by the MF solution.The parameter space can thus be divided into magnetic( M ), paramagnetic ( P ) and coexistence ( C ) regions. Inthe paramagnetic regions a single stable solution existsonly ( h n ↑ i = h n ↓ i ), while in the magnetic one two sta-ble magnetic (corresponding to magnetizations ± m ) andan unstable paramagnetic solution can be found. In thecoexistence region two magnetic and one paramagnetic µ / U Γ /U00.30.60.91.2 µ / U Γ /U increasing biasdecreasing bias MC P M PM P M P a) b)c) d)
FIG. 2: ”Phase diagrams” as a function of µ/U and Γ /U for T = 0 and different ε d values, a ) ε d /U = − . b ) − . c )0 and d ) 0 .
25. Notation of the regions and the lines are thesame as Fig. 1 stable solutions and two unstable magnetic solutions ex-ist. Therefore, these regions can be distinguished by thenumber of the solutions of Eq. (8). In the symmetriccase when h n ↑ i + h n ↓ i = 1, the easiest way to construct a”phase diagram” is to sweep possible values of h n ↑ i−h n ↓ i pairs, substitute them into Eq. (8) and count the numberof solutions for different parameter sets.In the paramagnetic region one can exploit the factthat Eq. (8) has a non-magnetic solution for every possi-ble parameter set. Therefore one can substitute h n ↑ i = h n ↓ i into Eq. (8) and search only for non-magnetic states.The region of stability for this solution can then be de-termined through a linear stability analysis .Fig. 1 shows a typical a magnetic ”phase diagram”as a function of µ/U and Γ /U for the symmetric non-equilibrium AIM at T = 0. At the “upper critical line” µ c the magnetic solution becomes unstable. In the mag-netic case, the effective levels corresponding to differentspins are not equally occupied and lie at different ener-gies. The magnetic solution becomes unstable when thevalue of the bias voltage reaches approximately the effec-tive energy of one of the two differently occupied levels.The critical line µ c in Fig 1 marks, on the otherhand, the border of stable paramagnetic solutions. Inthe special case, ǫ d = − U/ h n ↑ , ↓ i = 0 . µ +16Γ − U Γ /π = 0. This analytical result is nicely re-produced by our numerical stability analysis (see Fig. 1).The values of µ c and µ c are functions of Γ /U and ε d /U ; for increasing Γ /U , the coexistence and magneticregions disappear and only the non-magnetic solutionsurvives. In Fig. 2 the ”phase diagrams” can be seenas a function of µ/U and Γ /U at T = 0. Depending onthe value of ε d we can distinguish four different regions: empty regime ( ε d > mixed valence regime ( ε d ≈ local moment regime ( − U ≤ ε d ≤
0) and a doubly oc-cupied regime ( ε d < − U ) which behaves similarly to the empty regime by particle-hole symmetry. µ / U Γ / U µ / U Γ / U increasing biasdecreasing bias MC P MC PM P PM a) b)c) d) C FIG. 3: ”Phase diagrams” as a function of µ/U and Γ /U for ε d /U = − .
5. The temperature a ), T /U = 0 . b ), 0 . c )0 .
151 and d ) 0 . In the local moment regime , shown in Fig. 1 andFig. 2 a , for ε d ≈ − U/
2, there is approximately one elec-tron on the impurity forming a local spin moment. Inequilibrium, this finite magnetic moment is predicted onthe dot below a critical value of Γ /U , and this momentis destroyed by a large enough bias voltage. The coex-istence region only appears in this local moment regime and vanishes above ε d ≈ − U/ b ).In the empty regime ( ε d ≥ c and 2 d , the equilibrium magnetization completely dis-appears. However, surprisingly, the MF solution pre-dicts the appearance of local moments for small Γ’s and2 ε d ≤ µ ≤ ε d + U biases. This intriguing local mo-ment formation has a simple physical meaning: For largevalues of U/ Γ and 2 ε d ≤ µ ≤ ε d + U , the bias volt-ages are large enough to inject an electron to the emptylevel, however, they are not large enough to overcomethe Coulomb energy of injecting a second electron to thelocal level. Therefore, electrons pass through the dot oneby one, and a fluctuating magnetic moment appears onthe dot. Note that in this regime, the magnetization isinduced exclusively by the finite bias voltage.The overall effect of the temperature in the appliedMF-approximation is to destroy the magnetic moment onthe impurity and drive the system to be paramagnetic.Fig. 3 shows the temperature dependence of the ”phasediagram” in the symmetric case. The coexistence regiongradually vanishes for increasing temperatures, while themagnetic and paramagnetic regions get larger. Increasingthe temperatures further the magnetic region disappearstoo.The previous results are summarized in Fig. 4 for afixed ratio Γ /U = 0 .
05 that is already sufficiently smallto find the coexistence and magnetic regions in the localmoment regime. The whole ”phase diagram” is symmet-ric to ε d /U = − . ε d /U ≈ − .
5. For large asym-metries (double or zero equilibrium occupation) a non- -1.0 -0.5 0 ε d / U µ / U decreasing biasincreasing bias CM MP PP
FIG. 4: Magnetic (M), coexistence (C) and paramagnetic (P)regions as a function of µ/U and ε d /U values, for Γ /U = 0 . equilibrium magnetic solution appears, while in equilib-rium only the paramagnetic solution exists. This mag-netic region has also been observed although not dis-cussed in detail in Ref. 10.The mean-field solution also leads to the appear-ance of non-physical features in the transport proper-ties of the impurity. We calculated the transport prop-erties within the MF approximation using the Landauer-B¨uttiker formalism , but similar results can be obtainedapplying the Keldysh formalism . The current can beevaluated as I σ ≡ eh Z ∞−∞ d ε ( f L ( ε ) − f R ( ε )) | t σ ( ε ) | . (9)Here the transmission coefficient t σ = β σ p v k ′ /v k is nor-malized to the flux, with v k and v k ′ the velocities of inci-dent and transmitted electrons with wave numbers k and k ′ . Applying the large bandwidth approximation again,the current can be written for finite temperatures as I σ = e π ~ ∞ Z −∞ d ε Γ ( f L ( ε ) − f R ( ε ))( ε − ε d − U h n − σ i ) + Γ , (10)with h n σ i the non-equilibrium occupation numbers ob-tained from Eq. (7). This integral can be trivially evalu-ated at T = 0 temperature.Fig. 5 shows the current as a function of bias for twodifferent level positions in the strongly correlated regime.For ε d /U = − . ε d /U = 0 .
25, on the other hand, no hysteresis appearsbut the current shows a two-step behavior as a functionof bias voltage. In this empty regime , the MF equations thus account qualitatively correctly for the charging ofthe local level, but the kinks appearing in the I ( µ ) curveare due to the incorrectly predicted symmetry breakingand are thus again artifacts of the MF solution.In the inset of Fig. 5 we have plotted the polarizationof the current, P I ≡ ( I ↑ − I ↓ ) / ( I ↑ + I ↓ ), as a functionof µ . In the symmetric case the current is not polarized -2.0 -1.0 0 1.0 2.0 µ / U -0.2-0.100.10.2 I ( U e / h ) decreasing biasincreasing bias -1.0 0 1.0 µ/ U P I FIG. 5: The current as a function of the bias voltage for ε d /U = − . ε d /U = 0 .
25 (dot-ted line) for Γ /U = 0 .
02. Hysteresis can be observed in thecurrent in the local moment regime . Inset: the polarization ofthe current for ε d /U = − .
35 and Γ /U = 0 . due to the electron-hole symmetry, but for small electron-hole asymmetries, a hysteresis appears in the polarizationtoo. In general, the polarization is finite whenever a localmoment appears on the impurity and there is no electron-hole symmetry.To conclude, we have studied the Anderson model outof equilibrium in the framework of the scattering for-malism combined with a mean-field approximation. Thismethod, frequently used in molecular transport calcu-lations, incorrectly predicts a magnetic phase transitionas well as a bias-induced magnetic moment formation,accompanied by hysteresis in various physical quantitiesand the coexistence of multiple solutions. The MF ap-proach thus fails whenever correlations become impor-tant. These artifacts of the mean field approach shouldalert physicists who study transport through stronglycorrelated and magnetic molecules, and urge one to usemore sophisticated methods that avoid spontaneous sym-metry breaking and account for dynamical effects. Acknowledgment.
We are grateful L´aszl´o Szunyogh forvaluable discussions and remarks. This research has beensupported by Hungarian Grants No. OTKA NF061726,T046303, and F68726. P.W. Anderson, Phys. Rev. , 41 (1961). A. C. Hewson, in The Kondo Problem to Heavy Fermions(Cambridge University Press, Cambridge, England, 1993). D. E. Logan, M. P. Eastwood, and M. A. Tusch, J. Phys.: Condens. Matter , 2673 (1998). A. Mart´ın-Rodero, M. Baldo, F. Flores, and R. Pucci, SolidState Commun. , 911 (1982). A. A. Aligia and L. A. Salguero, Phys. Rev. B , 075307 (2004). A. L. Yeyati, A. Mart´ın-Rodero, and F. Flores, Phys. Rev.Lett. , 2991 (1993). A. A. Aligia, Phys. Rev. B , 155125 (2006). J. Rammer and H. Smith, Rev. Mod. Phys. , 323 (1986). D. Waldron, L. Liu and H. Guo, Nanotechnology ,424026 (2007) and references therein; D. Waldron, P.Haney, B. Larade, A. MacDonald and H. Guo, Phys. Rev.Lett. , 166804 (2006); K. Palot´as, B. Lazarovits, L.Szunyogh and P. Weinberger, Phys. Rev. B , 134421(2004) A. Komnik and A. O. Gogolin, Phys. Rev. B , 153102(2004). B. Horvatic, D. Sokcevic and V. Zlatic, Phys. Rev. B ,675 (1987). H. Kajueter and G. Kotliar, Phys. Rev. Lett. , 131(1996). B. Horv´ath, Diplome thesis (Budapest University of Tech-nology and Economics, Budapest, Hungary) (2006). M. B¨uttiker, Y. Imry, R. Landauer, and S. Pinhas, Phys.Rev. B31