Dynamic susceptibilities of the single impurity Anderson model within an enhanced non-crossing approximation

0.53 MB·English

by  Sebastian Schmitt

Dynamic susceptibilities of the single impurity Anderson model within an enhanced
non-crossing approximation
Sebastian Schmitt,1,2 Torben Jabben,1 and Norbert Grewe1
1Institut fu¨r Festko¨rperphysik, Technische Universita¨t Darmstadt, Hochschulstr. 6, D-64289 Darmstadt, Germany
2Lehrstuhl fu¨r Theoretische Physik II, Technische Universita¨t Dortmund,
Otto-Hahn Str. 4, D-44221 Dortmund, Germany
(Dated: January 19, 2010)
ThesingleimpurityAndersonmodel(SIAM)isstudiedwithinanenhancednon-crossingapprox-
imation (ENCA). This method is extended to the calculation of susceptibilities and thoroughly
0
tested,alsoinordertoprepareapplicationsasabuildingblockforthecalculationofsusceptibilities
1
and phase transitions in correlated lattice systems. A wide range of model parameters, such as
0
impurityoccupancy,temperature,localCoulomb repulsionandhybridizationstrength,arestudied.
2
Results for the spin and charge susceptibilities are presented. By comparing the static quantities

n to exact Bethe ansatz results, it is shown that the description of the magnetic excitations of the
a impuritywithintheENCAisexcellent,eveninsituationswithlargevalencefluctuationsorvanish-
J ing Coulomb repulsion. The description of the charge susceptibility is quite accurate in situations

9 wherethesinglyoccupiedionicconfigurationistheunperturbedgroundstate; however,itseemsto
1 overestimate charge fluctuations in the asymmetric model at too low temperatures. The dynamic
spin excitation spectrum is dominated by the Kondo-screening of the impurity spin through the

] conduction band, i.e. the formation of the local Kondo-singlet. A finite local Coulomb interaction
l
e U leads to a drastic reduction of the charge response as processes involving the doubly occupied
- impurity state are suppressed. In the asymmetric model, the charge susceptibility is enhanced for
r
t excitation energies smaller than theKondo scale TK dueto theinfluenceof valence fluctuations.
s
.
t PACSnumbers: 75.40.Gb,71.27.+a, 72.10.Fk
a
m
- I. INTRODUCTION The total Hamiltonian is then the sum of these three
d
terms
n
o The single impurity Anderson model (SIAM) de-
Hˆ =Hˆ +Hˆ +Vˆ . (4)
c scribes an impurity of localized f-states with local c f
[
Coulomb interaction embedded into a metallic host of
2 non-interactingc-bandelectrons.1 In its simplest version beEsvoelnvedthoeuxgahctltyhewitthheirnmtohdeynBaemthicesaonfsatthzemmeotdheold,c2a–n5
v it discards the possibility of a complex orbital structure
dynamic quantities can in general not be obtained ex-
4 of the impurity and models the local f-states through a
actly and one has to rely on approximations. The
1 two-fold degenerate s-orbital. The Hamiltonian for the
3 SIAM has been extensively studied with various meth-
impurity reads
3 ods, including the numerical renormalization group
. (NRG),6–9the(dynamic)density-matrixrenormalization
5 U
0 Hˆf = ǫf fˆσ†fˆσ+ 2 nˆfσnˆfσ¯ , (1) group ((D-)DMRG),10–13 quantum Monte Carlo (QMC)
9 σ (cid:18) (cid:19) methods14–16 and direct perturbation theory with re-
X
0 specttothehybridization.17–20Especiallywiththedevel-
v: with fˆσ† (fˆσ) and nˆfσ = fˆσ†fˆσ the usual creation (annihi- opment of the dynamical mean-field theory (DMFT),21
i lation) and number operatorsfor f-electronwith spin σ, wherethesolutionofaneffectiveSIAMrepresentsthees-
X respectively. Thelocalone-particleenergyisgivenbyǫf, sentialstep towards the solutionof the correlatedlattice
r and the local Coulomb interaction is the usual density- system,theinterestinaccurateandmanageableimpurity
a
densityinteractionproportionaltothematrixelementU. solvers has increased.
Thenon-interactingconductionelectronsaremodeledby In this work, we extend the well established enhanced
a single band of Bloch states with crystal momentum k non-crossing approximation (ENCA)22–25 to the calcu-
characterized by the dispersion relation ǫc, lation of the static and dynamic susceptibilities of the
k
impurity. Like many other approximations formulated
Hˆ = ǫc cˆ† cˆ . (2) within the direct perturbation theory with respect to
c k kσ kσ the hybridization,23,26–31 the ENCA is thermodynami-
k,σ
X cally conservingin the sense ofKadanoffandBaym.32,33
These two parts mix via a hybridization amplitude V , It extends the usual non-crossing approximation (NCA)
k
to finite values of the Coulomb repulsion U via the in-
1 corporationof the lowestordervertex corrections,which
Vˆ = V fˆ†cˆ +h.c. . (3)
√N k σ kσ are necessary to produce the correct Schrieffer-Wolff ex-
0
Xk (cid:16) (cid:17) change coupling and the order of magnitude of the low
2
energy Kondo scale of the problem. From the NCA may lead to an incorrect distribution of spectral weight.
it is well known that some pathological structure ap- The ENCA can be solved quite effectively on simple
pears at the Fermi level below a pathology scale34,35 desktopcomputers,andisnumericallynotverydemand-
T 10−1 10−2T . TheENCAremovesthecuspsin ing. Spectral functions can be calculated within some
path K
≈ −
spectral functions associated with this pathology22 and seconds to minutes while dynamic susceptibilities may
only a slight overestimation of the height of the many- take up to some hours. Additionally, it contains no free
body resonance at very low temperatures remains. As parameters and no fine-tuning is necessary. This makes
it will be shown in this work, the skeleton diagrams se- it especially interesting for involved lattice calculations.
lected within the ENCA suffer from an imbalance be-
tween charge and spin excitations and overestimate the
influence of charge fluctuations. Other than that, it has II. THEORY
no further limitations.
Despite the known limitations of the NCA it has been A. Direct perturbation theory
widely applied to more complex situations due to the
forthright possibility of extensions. For the SIAM out of In direct perturbation theory with respect to the hy-
equilibrium it is one of the few methods to incorporate bridization term Vˆ the “unperturbed” system is repre-
nonequilibriumdynamicsaswellasmany-bodyeffects.36 sented by the uncoupled (V = 0) interacting impurity.
k
Complexorbitalmultipletscanbeincludedinastraight- This is diagonalized by the ionic many-body states M
| i
forward manner and connections to experimental data
canbe made.37 However,incorporatingfinitevaluesofU Hˆf = EM XˆM,M , (5)
maychangethe many-bodyfeaturesnearthe Fermilevel M
X
considerably.23,38 Therfore,awelltestedextensionofthe where the operators Xˆ = M M are projectors on
ENCA in order to calculate susceptibilities at finite U M,M | ih |
the eigenstates M and are diagonal versions of the so-
withthe sameaccuracyas the spectrais desired. In par- | i
ticular, the calculation of lattice susceptibilities within called ionic transfer (or Hubbard) operators XˆM,M′ =
M M′ . For a simple s-shell the quantum numbers M
DMFT39–41 needs a reliable strategy for an effective im- | ih |
characterize the empty 0 , singly occupied with spin σ
purity, and the treatment of one- and two-particle exci- | i
σ anddoublyoccupied 2 impuritystateswiththecor-
tations of the same footing is of paramount importance. |resiponding unperturbed|eiigenvalues E = 0, E = ǫf
Calculations of lattice susceptibilities and phase transi- 0 σ
and E =2ǫf +U, respectively. Furthermore, the parti-
tions of the Hubbard model with the ENCA as the im- 2
tion function and dynamic Green function are expressed
purity solver are presented elsewhere.41,42
in terms of a contour-integrationin the complex plane,
Comparedtothe “numericallyexact”schemeslikethe
renormalization group methods (RG), exact diagonal- Z =Tre−βHˆ = dz e−βz Tr z Hˆ −1 (6)
ization (ED) or QMC, the direct perturbation theory 2πi −
IC
hsyassteimtsaatidcvaenrtraogrsess:te(mi)mTinhge afrpopmroxtihmeadtiiosncrsetairzeatfiroene ooff GA,B(iη)= 1 dz e−βz Tr z Hˆ (cid:2)−1Aˆ (cid:3) (7)
Z 2πi − ·
the conduction band (RG and ED) or imaginary time IC n(cid:2) (cid:3)
(QMC). The continuum of band states is kept through- z+iη Hˆ −1Bˆ ,
· −
outthecalculation,anddynamicGreenfunctionsarefor-
(cid:2) (cid:3) o
mulated with continuous energy variables. Thus, there with iη either a fermionic or bosonic Matsubara fre-
are no discretization-artefacts43 and there is no need for quencydepending onthe type ofthe operatorsAˆandBˆ.
artificial broadening parameters44,45 or z-averaging.46,47 The contour encirclesallsingularitiesofthe integrand,
(ii) The coupled integral equations for dynamic quan- which are sitCuated on the Imz = 0 and Im(z +iη) = 0
tities, which have to be solved numerically, are formu- axesinamathematicalleypositivesense. Performingthe
lated on the real frequency axis, which renders the non- trace over the c-electrons first, the reduced f-partition
trivial numerical analytic continuation of a finite set function Z , the f-electron one-body Green function
f
of Fourier coefficients48 or deconvolution49 unnecessary. F G and generalized susceptibility χ
The continuous-time quantum Monte Carlo approach Gσ ≡ fσ,fσ†can be expressed as M,M′ ≡
(CTQMC)16 avoids a systematic Trotter-error (i) but it XˆMM,XˆM′M′
is still plagued with the occurrence of a negative sign- dz
Z = e−βz P (z) , (8)
problem.50 f 2πi M
M IC
The discretization errors (i) are of special importance X
forself-consistentcalculationsliketheDMFT.There,the
1 dz
solutionofanimpuritymodelisusedtoconstructaguess Fσ(iωn)= e−βz (9)
Z 2πi
fortheGreenfunctionofthelatticewhichisthenusedto f IC h
yielda new effective “conductionband” for the impurity P0(z)Pσ(z+iωn)Λ0,σ(z,iωn)
model. Errorsinthetreatmentofthecontinuumofband
+P (z)P (z+iω )Λ (z+iω ,iω ) ,
states will be propagated by the iterative solution and σ¯ 2 n 2,σ¯ n n
i
3
χ (iν )= 1 dz e−βzχ (z,z+iν ) coupled integral equations,22,23,25,51
M,M′ n M,M′ n
−Z 2πi
f IC
= 1 dz e−βzΓM,M′(z,z+iνn) Σ0(z)= dx∆+(x)Pσ(z+x)Λ0,σ(z,x) (13)
−Zf IC2πi · Xσ Z
·ΠM′(z,z+iνn) . (10) Σσ(z)= dx∆−(x)P0(z−x)Λ0,σ(z−x,x)
Z
+ dx∆+(x)P (z+x)Λ (z+x,x)
2 2,σ
Here Π (z,z′) = P (z)P (z′) and the ionic propaga- Z
M M M
tors Σ2(z)= dx∆−(x)Pσ(z x)Λ2,σ(z,x) ,
−
σ Z
X
P (z)= M z Hˆ +Hˆ −1 M (11) where the vertex functions are given by
M c c
h |h − i | i
1
(cid:2) (cid:3) , (12) Λ (z,z′)=1+ dx∆+(x)P (z+x)P (z+z′+x)
≡ z E Σ (z) 0,σ σ¯ 2
M M
− − Z
(14)
which describe the dynamics of an ionic state M , are Λ (z,z′)=1+ dx∆−(x)P (z x)P (z z′ x) ,
| i 2,σ σ¯ − 0 − −
introduced. In equation (12) the ionic propagator is ex-
Z
pressed with the help of the ionic self-energy Σ (z).
M and ∆±(x) = ∆(x)f( x). The hybridization function
h...ic = Z1cTrc e−βHˆc... indicates that the trace is ∆(x) is constructed fro±m the c-band electrons
to be taken over the c-states only, and Z represents the
c
(cid:0) (cid:1)
partition function of the isolated c-system. In equations 1
∆(x)= V 2δ(x ǫc) , (15)
(9) and (10), ΛM,M′ and ΓM,M′ represent vertex func- N0 | k| − k
tions to be specified later. These equations are graphi- Xk
cally represented in Figure 1.
and f(x) = 1/(eβx+1) is the Fermi function with β =
After re-writing the Hamiltonian in terms of the ionic 1/T the inverse temperature (note ~=c=k =1).
B
transfer operators XˆM,M′, the resolvent operator is ex- ForthegeneralizedsusceptibilitiesχM,M′(z,z′),anad-
panded with respect to Vˆ, z Hˆ −1 = z Hˆ ditional set of integral equations has to be set up which
f
Hˆc −1 ∞n=0 Vˆ z−Hˆf −Hˆc(cid:2)−1−n. (cid:3) Conseq(cid:2)ue−ntly in−- iasresh1o6wnsucghrapfuhniccatilolynsi,nbFuitguforre e2a.chInquparnintucmiplen,utmhbereer
stheret(cid:3)irnegpPrtehseenitdhaetniot(cid:2)inty(1ˆ11)=for ct;hM(cid:3)e|Miionii|ccihpcr|ohMpa|gaatnodrs,ustihnigs MχM,o2n(zly,zt′h)earfoeucrofuupnlcetdi.onDs uχeMt,0o(zt,hze′)c,oχnMse,rσv(izn,gz′n),ataunrde
P
perturbationtheorycanbeformulatedwithtime-ordered of the ENCA, these equations are closely related to the
Goldstone diagrams, representing the dynamics of the ENCA expressions for the ionic self-energies; only some
ionic states M and their hybridization processes79. additionalbosoniclinesenteringandleavingthesitehave
| i
to be introduced.
Approximations are then introduced for the self-
Whereas the lowest order vertex corrections of equa-
energies Σ and vertex functions Λ and Γ .
M M,M′ M,M′ tion (14) together with the ionic propagators (11) and
In order to be able to describe non-perturbative many-
self-energies (12), are sufficient to furnish a conserving
body phenomena like the Kondo effect, certain classes approximation, the vertices Γ (z,z′) for the suscep-
M,M′
of diagrams have to be re-summed up to infinite order,
tibilities(10)havetobeiterateduptoinfiniteorder! The
resulting in a formulation in terms of skeleton diagrams,
transcriptionofthesegraphsintoformulasisstraightfor-
andcorrespondingcoupledintegralequationsfortherele-
ward but lengthy and will be omitted here.
vantdynamicfunctions. Derivingtheseequationsconsis-
For very large Coulomb interactions the spectral
tentlythroughfunctionalderivativesfromoneLuttinger-
weight in the doubly-occupied propagator is completely
WardfunctionalΦrenderstheseapproximationsthermo-
movedto energiesofthe orderofU. This leadsto a neg-
dynamically consistent.
ligibleinfluenceintheregionofaccessibleenergiesofthe
order of the bandwidth due to the large energy denom-
inators, and all terms involving P effectively drop out
2
of the equations. Therefore, for infinitely large Coulomb
repulsion U , the ENCA reduces to the usual NCA
→∞
and all equations presented above approach the ones al-
B. ENCA for generalized dynamic susceptibilities
ready known from the literature.26,52 In this sense, the
ENCA is still referred to as non-crossing even though
Within the ENCA, the ionic self-energies and propa- crossing diagrams are included. For approximations in-
gatorshave to be determined from the well knownset of volving “real” crossing diagrams see Refs. 23,29,30,53.
4
z
0 Λ0,σ(z,x) Λ2,σ(z,x)
Σ0(z)= σ x,σ , Σσ(z)= + , Σ2(z)=
(cid:1) (cid:2) (cid:3) (cid:4)
0
z
z z
0 σ¯
iν
z+iν n
n
Fσ(iωn)= σ + 2 , χ0,σ(iνn)= iν
n
(cid:5)0 iω ,σ (cid:6)σ¯ iω ,σ (cid:7)
n n
z
FIG.1: Diagrammaticequationsfortheionicself-energies (upper),f-electronGreenfunction(lowerleft)anddynamicsuscep-
tibility(lowerright,χ0,σ(iν)asanexample). Thedashed,wiggly anddouble-wiggly linesrepresenttheempty,singly occupied
(σ¯representstheoppositeofthespinσ)anddoublyoccupiedionicpropagators,respectively. Thetrianglesandboxesrepresent
theappropriate vertexfunctions. The closed full lines in theuppergraphs represent uncorrelated propagations in thec-band.
The functions χ (z,z′) fulfill the symmetry rela- bosonicMatsubarafrequency canbe analyticallycontin-
M,M′
tions ued to the real axis, iν ν iδ ν± (δ >0 infinitesi-
n
→ ± ≡
mal). The result of this procedure reads
χ (z∗,z′∗)=χ (z,z′)∗ (16)
M,M′ M,M′
∞ dω
and χ (ν+)= Y (ω,ω+ν) (22)
M,M′ M,M′
2πi
χ (z,z′)=χ (z′,z) , (17) Z−∞ "
M,M′ M,M′
Y (ω,ω ν)∗ ,
which are revealed by inspection of the perturbation ex- − M,M′ − #
pansion. The susceptibilities also obey the sum rules
where the symmetry relations (16) were used and the
χ (iν )= Xˆ δ (18)
M′,M n h MMi iνn,0 defect quantities
M′
X
χM′,M(iνn)=hXˆM′M′iδiνn,0 (19) YM,M′(x,y)= e−βx χM,M′(x+,y+) χM,M′(x−,y+)
Z −
M f
Xχ (iν )= Xˆ δ =δ , (20) h (23i)
M′,M n h MMi iνn,0 iνn,0
MX,M′ XM were defined. Due to the appearance of the exponential
factor, a direct numerical evaluation of the defect quan-
which arise from the completeness relation of the lo-
titiesgiventheχ isnotpossible,andadditionalsets
cal ionic states, 1 = M M . These sum rules M,M′
f M| ih | of integralequations have to be solvedfor the Y (cf.
transform into equivalent statements for the functions MM′
χM′,M(z,z′), P Ref. 54 and 55).
In the following, only spin-symmetric situations are
P (z) P (z′) considered, and the propagators for opposite spins are
χM′,M(z,z′)=− M z−z′M . (21) identified, i.e. Pσ =Pσ¯ P↑, χσ,σ =χσ¯,σ¯ χ↑,↑, χσ¯,σ =
XM′ − χσ,σ¯ χ↑,↓ .... ≡ ≡
≡
It can be analytically checked that the ENCA respects
these sum rules. The form of equation (21) explicitly
1. Magnetic susceptibility
revealstheconservingnatureofthisapproximation,since
the sum of sets of two-particle correlation functions is
determinedbythecorrespondingone-particlecorrelation Therelevantquantityforthemagneticsusceptibilityis
function, i.e. the ionic propagator. Insofar the equations the auto-correlation function of the z-component of the
(21) resemble generalized Ward-Identities. spinoperatorSˆz nˆ↑ nˆ↓ =Xˆ↑,↑ Xˆ↓,↓. Thistranslates
∼ − −
Inordertoobtainthephysicalsusceptibilityasafunc- into the linear combination
tion of a real frequency, the contour integration of equa-
χ (z,z′)=χ (z,z′) χ (z,z′) , (24)
tion (10) has to be performed, and then the external mag ↑,↑ ↑,↓
−
5
= δ + + + +
M,0
(cid:2)
(cid:1)
(cid:3)
(cid:4) (cid:5) (cid:6)
= δ + + + +
M,σ
(cid:8)
(cid:7)
(cid:9)
(cid:10) (cid:11) (cid:12)
+ + + +
(cid:13)
(cid:14) (cid:15) (cid:16)
= δ + + + +
M,2
(cid:18)
(cid:17)
(cid:19)
(cid:20) (cid:21) (cid:22)
FIG.2: Systemofintegral equationsforχM,0(z,z′),χM,σ(z,z′),andχM,2(z,z′) withintheENCA.Zig-zaglines inthevertex
part represent the ionic state M, which can either be empty, singly occupied or doubly occupied. The arrows entering and
leaving each diagram carry theexternal bosonic frequency iν.
which needs to be determined. Setting up the equations The derivation of the corresponding defect equation for
for this linear combination using the general equations Y alongthe lines ofthe definition (23) is straightfor-
mag
depicted in Figure 2 leads to the cancellationof all spin- ward, but will be omitted here for brevity.
symmetric terms.41 The function χ (z,z′) decouples
mag
from the χ (z,z′) and χ (z,z′), and only one inte-
0,M 2,M
gral equation results,
χ (z,z′)=Π (z,z′) (25) 2. Charge susceptibility
mag ↑
n
1 dxdy∆+(x)∆−(y) P (z′+x)P (z y)
− 2 0 − For the charge susceptibility the relevant quantity is
Z h the auto-correlationfunction of the charge operator nˆ =
+P0(z′−y)P2(z+x) χmag(z+x−y,z′+x−y) . nˆ↑ +nˆ↓ = Xˆ↑,↑ +Xˆ↓,↓ +2Xˆ2,2 leading to the dynamic
i o
6
function
χ (z,z′)=χ (z,z′)+χ (z,z′) (26)
charge ↑,↑ ↑,↓
c (z,z′)=Π (z,z′) 1+2 dx∆−(x) Λ (z′,x)
0,2 2 2,↑
+2 χ2,2(z,z′)+χ2,↑(z,z′)+χ↑,2(z,z′) (− Z h
=χ (z(cid:16),z′) χ (z,z′)+χ (z,z′) χ ((cid:17)z,z′) +Λ (z,x) 1 c (z x,z′ x) (30)
2,2 0,2 0,0 2,0 2,↑ 0,↑
− − − − −
+S0,↑,2(z,z′) . i
+2 dxdy∆−(x)∆−(y)P (z′ x)P (z y)
↑ ↑
− − ×
Inthesecondequalitysumruleslike(21)wereused. The Z
functionS0,↑,2(z,z′)isnotrelevantfordynamicsuscepti- c0,0(z x y,z′ x y) .
bilities,sinceitonlycontributesinthecaseofavanishing × − − − − )
external frequency iν =0.
n
The hole susceptibility, i.e. the auto-correlation func- All three linear combinations are now coupled which
tionoftheholeoperatorhˆ =2 nˆ =2Xˆ0,0+Xˆ↑,↑+Xˆ↓,↓, makes the numerical solution of the full system in-
−
isgivenbythesameexpression,onlywithadifferentcon- evitable.
tributionfromthesumrules. Bothofthesecontributions
The corresponding defect equations are again deter-
do notchange the dynamic susceptibility since they pro-
mined in a straight-forward way, but are omitted here
duceonlythenecessaryconstantshiftbetweenthe static
due to their length.
(iν =0)chargeandholesusceptibilityafterthecontour
n
integration,
χ (iν )=χ (iν )+4(1 nˆ )δ . (27)
hole n charge n −h i iνn,0 C. Numerical evaluation
The dynamicalquantities ofinterestarebest obtained
by setting up the integral equations for the linear com- The ionic propagators and defect propagators are
binations c = χ χ , c = χ χ and strongly peaked at the ionic threshold, and they even
0,0 0,0 2,0 0,↑ 0,↑ 2,↑
c =χ χ , which−read − developasingularityatzerotemperature.23,54,56 Thepo-
0,2 0,2 2,2
− sitionofthis ionicthresholdisnotknowexactlyapriori,
so we use self-generating adaptive frequency meshes to
c (z,z′)=Π (z,z′) 1+2 dx∆+(x) Λ (z′,x) represent all functions numerically.
0,0 0 0,↑
( Z h After obtaining the ionic propagators from the set
of integral equations Eq. (13), Eqs. (25) and (28) –
+Λ (z,x) 1 c (z+x,z′+x) (28)
0,↑ − 0,↑ (30) can be solved to yield the two-particle quantities
i χ (ω±,ω ν+)andc (ω±,ω ν+). Withtheseat
+2 dxdy∆+(x)∆+(y)P (z′+x)P (z+y) mag ± M,M′ ±
↑ ↑ × hand, the corresponding defect equations can be solved
Z and physical susceptibilities along the lines of Eq. (22)
c (z+x+y,z′+x+y) can be extracted.
0,2
× ) All integral equations are solved via Krylov subspace
methods,57 where — starting from a sophisticated guess
for the unknown functions — the equations are iter-
ated until convergence is reached. In order to ac-
c (z,z′)=Π (z,z′) (29) celerate convergence and stabilize the iteration proce-
0,↑ ↑
( dure Pulay-mixing58 between different iterations is im-
plemented (see, for example, Ref. 59 or 60).
dx∆−(x) Λ (z′ x,x)
0,↑ − Unfortunately,the systemforthe defectquantitiesbe-
Z h comes singular for vanishing external frequency ν 0.
+Λ0,↑(z−x,x)−1 c0,0(z−x,z′−x) This is seen at the sum rule (21), which translates→into
i Y (ω,ω +ν) = 2πiξ (ω)/(ν +iδ) implying
+ dx∆+(x) Λ2,↑(z′+x,x) thaMt′thMe′,YM become o−f theMorder of 1 for small ν,
Z PwhileallteMrm,Ms′intheintegralequationsstaνyattheorder
+Λ(cid:2) (z+x,x) 1 c (z+x,z′+x)
2,↑ − 0,2 one. This becomes of some importance when extracting
+ dxdy∆+(x)∆−(y) P (z′+(cid:3) x)P (z y) static susceptibilities from calculations of the dynamic
2 0
− susceptibilities, where a very small but finite frequency
Z h
is used. The resulting convergence problems are proba-
+P0(z′ y)P2(z+x) c0,↑(z+x y,z′+x y) bly connected to the ones already mentioned in Ref. 61.
− − − ) More details on the numerics can be found in Ref. 41.
i
7
and W the half bandwidth of the conduction electron
π∆Aρf(ω)
1.2 DOS ρc(ω). The hybridization matrix element was as-
sumed to be local, i.e. momentum independent V V.
k
→
1 U=1.7 ∆A=0.2 Asalreadymentioned,thepathologyoftheNCAmani-
0.8 24 ǫf=−1 festsitselfintheoverestimationoftheheightoftheASR
∞ and a violation of Fermi liquid properties for too low
c-band temperatures as well as in in situations with large va-
0.6
lencefluctuations.26,34,35,54 Thispathologicalbehavioris
0.4 strongestforthe caseofaspin-onlydegeneracy(N =2),
which is considered in this work.
0.2 -0.03 0 0.03 Figure 3 shows the one-particle f-electron DOS
ρ (ω) = 1ImF (ω iδ) and the imaginary part of the
0 f π σ −
-4 -2 0 2 4 6 8 totalself-energyImΣtfot(ω−iδ)=−Im[1/Fσ(ω−iδ)]cal-
ω culated within the ENCA and the NCA (U = ). As
∞
ImΣtot(ω−iδ) canbe seen,the ENCAdoes notoverestimatethe height
f
6 oftheKondopeakorviolatetheFermiliquidpropertyof
the total self-energy ImΣtot(0 iδ) ∆ down to tem-
f − ≥ A
5 peratures of half the Kondo temperature. In contrast,
the U = NCA curves in the graphs do violate these
4 ∞
limits for the same parameter values.
ThefactthattheENCAperformsbetterthantheNCA
3
incomparablesituationsisaconsequenceofabetterbal-
2 -0.02 0 0.02 ancebetweendifferent kinds ofperturbationalprocesses.
The importance of such a balance has been pointed out
1 repeatedly.23,62,63
EventhoughtheperformanceoftheENCAisconsider-
0
-4 -3 -2 -1 0 1 2 3 4 ablyenhancedovertheNCA,iteventuallyoverestimates
ω the height of the ASR for even lower temperatures and
still misses to produce the correct T = 0 Fermi liquid
FIG.3: Rescaled f-electronspectra(uppergraph)andtotal relations. Further improvements can be archivedvia the
self-energy (lower graph) as a function of energy for various incorporation of higher-order diagrams.23
values of the Coulomb repulsion U. The temperature was
fixed at half the Kondo temperature T ≈TK/2 of each case,
the ionic level was chosen ǫf = −1 and the Anderson width
was ∆A = 0.2. The conduction band DOS was that of a A. Benchmarking the ENCA
three-dimensional tight binding lattice with half bandwidth
W = 3 (“c-band” in the upper graph). The insets show the
1. Static susceptibilities in the symmetric case
low energy intervalaround the Fermilevel, where theenergy
ismeasuredinunitsofthecorrespondingKondotemperature
and the horizontal lines indicate theT =0 limiting values. In order to obtain a better understanding of possible
shortcomings of the ENCA it is worthwhile to consider
charge and spin excitations separately and benchmark
III. RESULTS
them against some exactly known results. The thermo-
dynamics of the SIAM can be obtained exactly from the
As it is well known, the SIAM exhibits the Kondo ef- Bethe ansatz method.3,4,64–66 At zero temperature and
fect for sufficiently large Coulomb repulsion U and the forthesymmetriccase(ǫf = U)withaflatconduction
single particle level below the Fermi energy, ǫf < 0. A band of infinite bandwidth (W−2 ) the static suscep-
very prominent manifestation of that effect is the emer- tibilities can even be obtained in→cl∞osed form,66,67
genceoftheAbrikosov-Suhlresonance(ASR)intheone-
pteamrtpicelreatduernessitoyftohfesotradteesro(fDtOheS)KaotndtoheteFmeprmeriatleuvreelTfKor. χexact(T =0)= (gµB)2 1+ 1 π2∆UAdx e x−1π62x
Within the ENCA an order of magnitude estimation for mag 4TL √π Z0 (cid:8)√x (cid:9)!
T is given by22 (32)
K
min(W,U) 2U∆
T = √J e−π/J , J = A (31) with
K 2π −ǫf(ǫf +U)
∆ πU π∆
whereJ istheeffectiveantiferromagneticSchrieffer-Wolff T =U A exp + A (33)
L
exchange coupling, the Anderson width ∆A = πρc(0)V2 r2U (cid:26)−8∆A 2U (cid:27)
8
1000 exactT=0 10
TK
(0) 100 (1/2)TK
χmag 10 1 T=0.05
∆A=0.2
1 ǫf=−U/2
0.1
(0) 1 macghnaertgiec
χcharge 0.1 0.01 exaecxtacmtacghnaertgiceTT==00
0.01
0 1 2 3 4 5 6 0 0.2 0.4 0.6 0.8 1
U U
1+U
FIG.4: Staticmagnetic(upper)andcharge(lower) suscepti- FIG.5: Magneticandchargesusceptibilityforthesymmetric
bilityforthesymmetriccaseǫf =−U/2asafunctionofU for impurity (ǫf =−U/2) as a function of the Coulomb interac-
finitetemperaturesTK and TK/2 in comparison to theexact tion U at a fixed temperature T =0.05 and for an Anderson
T = 0 results of equations (32) and (34). The calculations width ∆A = 0.2. Notice, the x-axis is scaled as U/(1+U),
weredoneforaconstantc-electronDOSwithhalfbandwidth so that U = ∞ corresponds to U/(1+U) = 1. The arrow
W =10 and an Anderson width ∆A=0.3. on theright of the figurecorresponds to the valueof a Curie
susceptibility χmag/(gµB)2 = 1/(4T) = 5 of a free spin at
temperature T =0.05.
and
χexact (T =0)= 1 2 ∞dx e−π2∆UAx2 . of the order of 10−4 and decreases for larger U. The
charge π U∆ 2 susceptibilitycalculatedatν =10−5 thendoesnolonger
r A Z−∞ 1+ 2∆UA +x representthe staticlimitanymore,andthedecreaseseen
(cid:16) (cid:17) (34) inFigure4isproduced,whichisthereforenotindicating
a shortcoming of the ENCA method. The deviation can
Apart from the small coupling correction ∆A in the be cured by choosing a smaller value for the external
∼ U
exponent, T exactly coincides with the Kondo temper- frequency.
L
ature of equation (31). The charge susceptibility (lower graph in Figure 4)
In the graphs of Figure 4 the exact zero temperature showsnosignificanttemperaturedependenceforT =T
K
magnetic (upper) and charge (lower) susceptibilities of andT /2andliesrightontopofthe exactT =0result.
K
equations(32) and(34)(solidredlines)arecomparedto The deviations for U 4 are explained in the same way
≥
the ENCA susceptibilities (blue dots) for ǫf = U/2 as as for the magnetic susceptibility described above.
−
afunctionofU. FortheENCAcurvestwocharacteristic Thesereasonsforthe deviationatlargerU canbefur-
temperatures T and T /2 are chosen. ther confirmed by calculating the static susceptibilities
K K
The characteristic exponential U dependence of the for a fixed finite temperature T = 0.05 as a function of
magnetic susceptibility is essentially the same as for the the Coulomb interaction, which are shown in Figure 5.
exact Bethe ansatz result, but the absolute height is The figure compares the static magnetic (red dots) and
somewhat different. At T = T the magnetic suscep- charge (blue dots) susceptibilities for the symmetric sit-
K
tibility is roughly70%ofits zerotemperature saturation uation with the exact T =0 results.
value. However, at T = T /2 the susceptibility has al- Thechargesusceptibilitydecreasesmonotonicallywith
K
most saturated and the agreement with the T =0-value increasing Coulomb repulsion as expected. It follows
is quite good. theexactzerotemperaturesusceptibilityveryaccurately,
The deviations for U > 4 are due to the method used which again confirms its weak temperature dependence
to extract the static susceptibilities: The static limit is in symmetric situations at low temperatures.
obtained by evaluating the dynamic susceptibility at a Themagneticsusceptibility(reddots)doesagreewith
smallbutfinite externalfrequency,inourcaseν =10−5. the exact T = 0 solution for low Coulomb repulsions
Since the ENCA represents a conserving approximation, but deviates for U/(1+U) & 1/3 (U & 1/2). This is
the results obtained with this method agree with the a finite temperature effect, since for values of U & 1/2
static susceptibility obtainedfroma derivative ofa ther- the Kondo temperature is smaller than the chosen finite
modynamicpotential,orsolvingseparateequationsasin temperature of T = 0.05. Consequentlya, for U & 1/2
Otsuki et al.68. This is valid as long as the minimal fre- we are not in the low temperature regime, and the sus-
quency is negligible compared to the lowest energy scale ceptibility is not well described by its T = 0 value. The
in the problem. magneticsusceptibilitydoesnotgrowexponentiallywith
However, the Kondo temperature for U = 4 is only U as for T = 0, but instead saturates for U at
→ ∞
9
µ2 χ2mag(0)
eff 10
0.25 U=0.01
ǫf/V2=−01 T∆A==0.00.52 Uff=rree2ee
1
−2 U=4
free
0.1
0.125
0.01
0.001
0 -4 -3 -2 -1 0 1 2 3 4
0.001 0.01 0.1 1 ǫf+U/2
T
FIG. 7: Static magnetic susceptibility for a fixed T = 0.05
FIG. 6: Temperature dependent squared effective magnetic and ∆A =0.2 as functions of the ionic level position ǫf rela-
momentwithintheENCA(coloreddots)forafixedCoulomb tivetothehalf-fillingvalue−U/2forvariousvaluesofU. The
repulsion U = 4V2 and three different f-level positions conductionbandwaschosentobeconstantwithahalfband-
ǫf =0,−V2,−2V2. ThesolidgreycurvesaretheexactBethe widthofW =10. Curveswithoutdots(“free”)arecalculated
ansatz resultsfor thesame parameter values(afterOkiji and without two-particle interactions, i.e. with the particle-hole
Kawakami69). The calculations were done for a constant c- propagator of equation (36).
electron DOS with half bandwidth W =10 and ∆A=0.167.
one-particle DOS and the NCA-type of approximations
the asymptotic value of the Curie susceptibility of a free
wouldbe expectedtoyieldresultsoflowerquality. How-
spin, χ =1/(4T)=5,whichisindicatedbyanarrow
mag ever, as it can be seen, the magnetic excitations are still
on the right border of the graph.
described very accurately.
Calculationsoftheeffectivesquaredlocalmomentsfor
a wide range of Coulomb interactions and hybridization
2. Static susceptibilities in the asymmetric case
strength(notshown)resembletheexactsolutionsknown
from the literature7,8,66 and the characteristics of the
As a first test for the ENCA in the asymmetric sit-
different asymptotic regimes are well reproduced by the
uation Figure 6 compares the square of the effective
ENCA.
screened local magnetic moment of the impurity,
In Figure 7 the magnetic susceptibility is systemati-
1 cally examined as a function of the ionic level position
µ2 = Tχ (ν =0) , (35)
eff (gµ )2 mag ǫf and for fixed T and U. Also shown are the suscep-
B
tibilities without explicit two-particle interactions (lines
to the exact Bethe ansatz result for three different ionic
withoutdotslabeledas“free”),i.e.thelocalparticle-hole
level positions as a function of temperature. The solid propagator Pf(0) calculated via
grey lines are the exact Bethe ansatz solution, which is
taken from Ref. 69, while the colored points are ENCA ∞
calculations for the same parameter values. The half Pf(ν)= dωf(ω)ρf(ω) (36)
bandwidth was taken to be W = 10, which should be Z−∞
largeenoughto be comparableto the Bethe ansatz solu-
F (ω+ν+iδ)+F (ω ν iδ) ,
↑ ↑
tion where W = . × − −
∞
The ENCA slightly overestimates the squared effec- (cid:2) (cid:3)
tivemomentbutallcharacteristicfeaturesareessentially where F↑ represents the one-particle f-Green function
the same as for the Bethe ansatz. Especially the shape (cf.equation(9)),ρf(ω)thecorrespondingspectrumand
and the relative height of the curves is in remarkable f(ω) the Fermi function.
agreement: All three ENCAcurvescanbe broughtto lie All curves are symmetric around ǫf + 1U = 0 which
2
righton top of the exactBethe ansatz results when they just reflects the particle-hole symmetry of the model.
are rescaled with one single factor. This indicates that The particle-hole propagators shows the expected max-
theENCAproducesaslightlymodifiedKondoscale,but ima approximately situated at the positions of the Hub-
otherwisedescribesthestaticmagneticpropertiesalmost bard peaks in the one-particle DOS ǫf + 1U 1U.
2 ≈±2
exactly. This is especially remarkable for the intermedi- For very small U = 0.01, the susceptibility calculated
ate valence situation with ǫf = 0, where the empty and with the ENCA is indistinguishable from the particle-
singly occupied ionic configurations are almost degener- hole propagator, as a consequence of the near lack of
ate. In such situations stronger pathologies occur in the two-particle correlations.
10
limν→02π(gχµ(νB,T)2)2Imχν(ν,T) χcharge(0)
3.5
10
U=1.7
3 2
4
1
8
2.5 ∞
Korringa
2 0.1
U=1.7
U=2
1.5
U=4
0.01
1/(8T)
1
0.5 0.001
0.001 0.01 0.1 1 10 1e-05 0.0001 0.001 0.01 0.1 1 10
T/TK T/W
FIG. 8: limν→0 2πχ(mgµagB()ν2,T)2Imχmaνg(ν,T) evaluated as a func- FIG. 9: Temperature dependent static charge susceptibil-
tion of temperature and various Coulomb repulsions within ity χcharge(0) for various Coulomb repulsions U (∆A = 0.2,
the ENCA. For comparison the U =∞-NCA curve is shown ǫf =−1). ThearrowattheleftbordershowstheexactT =0
as well. A 3d-SCDOS was used for the conduction electrons limit for the symmetric case U = 2. The high temperature
and ∆A =0.3. asymptotic form of a non-interacting impurity (χ = 1/(8T))
is shown as well (thick line). The colored dots, which are
not connected by a line, denote the asymptotic susceptibili-
For larger Coulomb interactions, the ENCA suscepti- tiescalculatedwithoutanytwo-particleinteractions,i.e.from
bilityshowsonlyonebroadmaximumaroundǫf+1U =0 Pf(ν =0), seeequation (36). Theconductionband wascho-
2 sen to be that of a 3d-SC lattice with a half bandwidth of
(half filling), which grows in height and width with in-
W =3.
creasingU. Theenhancementofthesusceptibilityisdue
to the increasing local magnetic moment with larger U.
The plateau which develops around zero, is due to the
merical pre-factors might change slightly, but the left-
stability of the local moment as long as the singly occu-
handsideis stillexpectedto be ofthe orderofunity due
pied ionic configuration is stable, i.e. the lower Hubbard
to universality of the SIAM at low energies.
peak being below and the upper above the Fermi level,
For the NCA, the quantity lim Imχ (ν)/ν is
andthetemperatureisnottoolowcomparedtoT . But ν→0 mag
K known to diverge at T = 0,56 which is reproduced by
as soon as one of the Hubbard peaks extends over the
the U = curve in the figure. However, the ENCA
Fermi level, i.e. ǫf + 12∆A <0 or ǫf +U − 12∆A >0, the (finite U v∞alues) performs considerably better than the
moment is destabilized. For both Hubbard peaks below
NCA. The curves still slightly increase for temperatures
(above) the Fermi level, the impurity is predominantly
T /T /2, but they eventually saturate at a finite value
doubly occupied (empty) and the magnetic susceptibil- K
and do not diverge. This represents a considerable im-
ity drops drastically. The curve then rapidly approaches
provementof the qualitative behavior of the ENCA over
theparticle-holepropagator,indicatingthatexplicittwo-
the NCA.
particle interactions are unimportant.
The temperature dependent static charge susceptibil-
The reproduction of the correct results for the effec-
ity is shown in the Figure 9 for various values of U as
tively non-interacting limit at U = 0.01 as well as the
a function of temperature. For high temperatures and
empty- or fully occupied regimes is quite remarkable.
U < 4 the static susceptibility behaves effectively non-
As it was already mentioned earlier, the NCA does
interacting with χ =1/(8T) as expected.
violateFermiliquidpropertiesforverylowtemperatures. charge
ForU =4the susceptibility stillhas the characteristic
Another indication, in addition to the imbalance of the
1/T dependence for high T, but with a prefactor more
imaginarypart ofthe totalself energy at the Fermi level
closely to 1/16. This canbe understoodsince in this sit-
(-ImΣ(0+iδ)<∆ ),stemsfromthezerofrequencylimit
A
uation the upper Hubbard peak (incorporating roughly
of the imaginary part of the dynamic susceptibility. For
half the spectral weight) is energetically just above the
the Fermi liquid at T = 0 it has to obey the so called
Korringa-Shiba relation70 upper band edge of the c-band and therfore the accessi-
blespectralweightfortwo-particleexcitationsisapprox-
(gµ )2 Imχ (ν) imately halved.
B mag
lim =1 . (37)
ν→02πχmag(ν)2 ν The rapid drop of the susceptibility for U =4 at tem-
peratures T & W is attributed to inaccuracies in the
Thefunctiononthelefthandsideofthisrelationisshown numerics for solving the integral equations. However, in
inFigure8asafunctionoftemperatureforvariousvalues this effectively non-interacting regime the susceptibility
of U. The explicit form (37) holds for a flat infinitely- can be calculated without explicit two-particle interac-
wide conductionband DOS. For a different DOS the nu- tions via the particle-hole propagator of equation (36).