Sum-Rule Conserving Spectral Functions from the Numerical Renormalization Group

Sum-Rule Conserving Spectral Functions from the Numerical Renormalization Group

Andreas Weichselbaum and Jan von Delft

Physics Department, Arnold Sommerfeld Center for Theoretical Physics, and Center for NanoScience, Ludwig-Maximilians-Universita¨t Mu¨nchen, D-80333 Mu¨nchen, Germany

(Received 1 August 2006; published 16 August 2007)

We show how spectral functions for quantum impurity models can be calculated very accurately using a complete set of discarded numerical renormalization group eigenstates, recently introduced by Anders and Schiller. The only approximation is to judiciously exploit energy scale separation. Our derivation avoids both the overcounting ambiguities and the single-shell approximation for the equilibrium density matrix prevalent in current methods, ensuring that relevant sum rules hold rigorously and spectral features at energies below the temperature can be described accurately.

DOI:10.1103/PhysRevLett.99.076402 PACS numbers: 71.27.+a, 73.21.La, 75.20.Hr

Quantum impurity models describe a quantum system with a small number of discrete states, the ‘‘impurity,’’

coupled to a continuous bath of fermionic or bosonic excitations. Such models are relevant for describing trans- port through quantum dots, for the treatment of correlated lattice models using dynamical mean field theory, or for the modeling of the decoherence of qubits.

The impurity’s dynamics in thermal equilibrium can be characterized by spectral functions of the typeA^BC! R_dt

2e^i!thBt^ Ci^ T. Their Lehmann representation reads A^BC! X

a;b

hbjC^jaie^Ea

Z hajB^jbi!E_ba; (1) withZP

ae^EaandE_ba E_bE_a, which implies the sum rule R

d!A^BC! hB^ C^i_T. In this Letter, we de- scribe a strategy for numerically calculatingA^BC!that, in contrast to previous methods, rigorously satisfies this sum rule and accurately describes both high and low frequencies, including!&T, which we test by checking our results against exact Fermi-liquid relations.

Our work builds on Wilson’s numerical renormalization group (NRG) method [1]. Wilson discretized the environ- mental spectrum on a logarithmic grid of energies ⁿ (with >1, 1nN ! 1), with exponentially high resolution of low-energy excitations, and mapped the im- purity model onto a ‘‘Wilson tight-binding chain,’’ with hopping matrix elements that decrease exponentially as ⁿ⁼² with site index n. Because of this separation of energy scales, the Hamiltonian can be diagonalized itera- tively: adding one site at a time, a new ‘‘shell’’ of eigen- states is constructed from the new site’s states and theM_K lowest-lying eigenstates of the previous shell (the so-called

‘‘kept’’ states), while ‘‘discarding’’ the rest.

Subsequent authors [2–10] have shown that spectral functions such as A^BC! can be calculated via the Lehmann sum, using NRG states (kept and discarded) of those shellsnfor which!ⁿ⁼². Though plausible on heuristic grounds, this strategy entails double-counting

ambiguities [5] about how to combine data from successive shells. Patching schemes [9] for addressing such ambigu- ities involve arbitrariness. As a result, the relevant sum rule is not satisfied rigorously, with typical errors of a few percent. Also, the thermal density matrix (DM) ^ e^{H^}=Z has until now been represented rather crudely using only the single N_Tth shell for which T’ ^1=2NT1[8], with a chain of length NN_T, resulting in inaccurate spectral information for!&T. In this Letter we avoid these problems by using in the Lehmann sum an approximate but complete set of eigenstates, introduced recently by Anders and Schiller (AS) [11].

Wilson’s truncation scheme.—The Wilson chain’s ze- roth site represents the bare impurity Hamiltonianh^₀with a set ofd₀ impurity statesj₀i. It is coupled to a fermionic chain, whose nth site (1nN) represents a set ofd statesj_ni, responsible for providing energy resolution to the spectrum at scaleⁿ⁼². For a spinful fermionic band, for example, _n2 f0;";#;"#g, hence d4. (Bosonic chains can be treated similarly [10].) The Hamiltonian H^ H^_N for the full chain is constructed iteratively by adding one site at a time, usingH^_nH^n1h^_n (acting in adⁿd₀-dimensional Fock spaceFnspanned by the basis statesfj_ni j₀ig), whereh^_nlinks sitesnandn 1 with hopping strength ⁿ⁼². Since the number of eigenstates of H^_n grows exponentially with n, Wilson proposed the following iterative truncation scheme to nu- merically diagonalize the Hamiltonian: Let n₀ be the last iteration for which a complete setfjsi^K_n0gof kept eigenstates of H^_n0 can be calculated without trunction. For n > n₀, construct the orthonormal eigenstatesfjsi^X_ngofH^_n(thenth

‘‘shell’’), with eigenvalues Eⁿ_s, as linear combinations of the kept eigenstates jsi^K_n1 ofH^n1 and the statesj_niof siten,

js⁰i^X_n X^K

_ns

j_ni jsi^K_n1A_KXⁿ_ss⁰; (2) with coefficients arranged into a matrix A_KXⁿ whose ele- PRL99,076402 (2007) P H Y S I C A L R E V I E W L E T T E R S week ending

17 AUGUST 2007

0031-9007=07=99(7)=076402(4) 076402-1 © 2007 The American Physical Society

ments are labeled by ss⁰. The superscript XK or D indicates that the new shell has been partitioned into

‘‘kept’’ states (say the M_K lowest-lying eigenstates of H^_n) to be retained for the next iteration and ‘‘discarded’’

states (the remaining ones). Since h^_n acts as a weak per- turbation (of relative size ¹⁼²) on H^n1, the d-fold degeneracy of the states j_ni jsi^X_n1 is lifted, resulting in a characteristic energy spacing ⁿ⁼² for shell n.

Iterating until the spectrum of low-lying eigenvalues has reached a fixed point (fornN, say), one generates a set of eigenstatesfjsi^X_ng with the structure of matrix product states [12] (Fig. 1). The states generated for the last Nth shell will all be regarded as discarded [11].

Anders-Schiller basis.—Recently, AS have shown [11]

that the discarded states can be used to build a complete basis for the whole Wilson chain: the statesfjsi^X_ngdescrib- ing the nth shell are supplemented by a set of d^Nn degenerate ‘‘environmental’’ states fje_ni j_Ni j_n1igspanning the rest of the chain to construct the set of statesfjsei^X_n je_ni jsi^X_ng. These reside in the complete Fock spaceFN of the full chain, spanningFN ifnn₀. Ignoring the degeneracy-lifting effect of the rest of the chain, these states become approximate eigenstates of the HamiltonianH^_Nof thefullchain (‘‘NRG approximation’’), H^_Njsei^X_n ’Eⁿ_sjsei^X_n; (3) with eigenenergiesindependentof the (d^Nn)-fold degen- erate environmental index e_n. (This will facilitate tracing out the environment below.) By construction, we have

Dmhsejs⁰e⁰i^D_{n mnene}⁰_nss⁰ and

Kmhsejs⁰e⁰i^D_n

0; mn _ene⁰_nA_KK^m1. . .A_KDⁿ _ss⁰; m < n: (4) The discarded states of shellnare orthogonal to the dis- carded states of any other shell, and to the kept states of that or any later shell. Combining the discarded states from all shells thus yields a complete set of NRG eigenstates of H^_N, the ‘‘Anders-Schiller basis,’’ that span the full Fock spaceFN (P

nhenceforth stands forPN n>n₀):

1^d0dNX

se

jsei^K_n0^K_n0hsej X

n

X

se

jsei^D_n^D_nhsej: (5) Local operators.—Let us now consider a ‘‘local’’ opera- tor B^ acting nontrivially only on sites up to n₀. Two

particularly useful representations are B^ X

ss⁰e

jsei^K_n0Bⁿ_KK⁰_ss^0K_n0hs⁰ej X

n

X

KK

XX⁰

B^ⁿ_XX⁰: (6) The left equality, writtenB^ B^ⁿ_KK⁰in brief, represents the operator in the complete basis set fjsei^K_n0g, with matrix elements known exactly numerically (possibly up to fer- mionic minus signs depending on the environmental states, but these enter quadratically in correlation functions and hence cancel). The right-hand side (RHS) of Eq. (6) ex- presses B^ in the AS basis and is obtained as follows:

starting from B^ⁿ_KK⁰, one iteratively refines the ‘‘kept- kept’’ part ofB^ from, say, the (n1)th iteration in terms of the NRG eigenstatesfjsei^x_ngof the next shell, including both kept and discarded states (XK; D),

B^ ⁿ¹KK X

XX⁰

X

ss⁰e

jsei^X_nBⁿ_XX⁰_ss^0X_n⁰hs⁰ej X

XX⁰

B^ⁿ_XX⁰; (7) thereby defining the operatorsB^ⁿ_XX⁰, with matrix elements Bⁿ_XX⁰_ss⁰ A_XK^nyBⁿ¹_KK A_KXⁿ0_ss⁰. Splitting off allXX⁰ KKterms (DD,KD,DK) and iteratively refining eachKK term untilnN, we obtain the RHS of Eq. (6). It has two important features. First, the matrix elements of the time- dependent operator Bt ^ e^{iHt^} Be^ ^{iHt^} , evaluated within the NRG approximation,B^XXn ⁰t_ss⁰ ’ B^XXn ⁰_ss⁰e^itEnsEns0, contain differences of eigenenergies from the same shell only, i.e., calculated with the same level of accuracy.

Second, by excluding KK terms it rigorously avoids the double-counting ambiguities and heuristic patching rules plaguing previous approaches [2–10].

Thermal averages. —To calculate thermal averages h. . .i_T Tr^. . ., we write the full density matrix (FDM) ^ e^{H^}=Z using the NRG approximation Eq. (3),

^ ’X

n

X

se

jsei^D_n e^Ens Z

Dnhsej X

n

w_n^ⁿ_DD; (8) wherew_nd^NnZ^D_n=ZandZ^D_n P_D

s e^Ens. The RHS of Eq. (8) expresses^as sum over^ⁿ_DD, the density matrix for the discarded states of shell n, properly normalized as Tr^ⁿ_DD 1, and entering with relative weightw_n, with P

nw_n1. Similarly, for spectral functions we have h. . .i_T X

n

w_nh. . .i_n; A! X

n

w_nAn!; (9) where the averages h i_n and spectral functions An! are calculated with respect to^ⁿ_DDof shellnonly.

Previous strategies [4–11] for thermal averaging amount to using a ‘‘single-shell approximation’’w_nnNTfor the density matrix and terminating the chain at a length N N_Tset byT’^1=2NT1. As a result, spectral features on scales !T, which would require a longer chain, are described less accurately [see Figs. 2(a) and 2(b)]. Our FIG. 1. Diagram for the kept (or discarded) matrix product

state js⁰i^K_n (or js⁰i^D_n): the nth box represents the matrix block A_KXⁿ, its left, bottom, and right legs carry the labels of the states jsi^K_n1,j_ni, andjs⁰i^K_n (orjs⁰i^D_n), respectively.

PRL99,076402 (2007) P H Y S I C A L R E V I E W L E T T E R S week ending 17 AUGUST 2007

076402-2

novel approach avoids these problems by using the full density matrix (FDM), summed overallshells, letting the weighting functionw_nselect the shells relevant for a given temperature yielding a smooth T dependence [see Fig.2(c)]. Sincew_n has a peak width of five to ten shells depending on,dandM_Kand peaks atnvalues somewhat above N_T [arrow Fig. 2(b)], spectral information from energies well belowTis retained.

Let us now consider the spectral functionA^BC!, for local operatorsB^ andC^. Equations (4), (6), (8), and (9) can be used to evaluate hB^tC^i_n. Fourier transforming the result we find (sums overss⁰and_nimplied)

A^BCn ! Xⁿ

m>n₀

X

KK

XX⁰

C^m_X⁰_X^mn_{XX s}⁰_sB^m_XX⁰_ss⁰!E^m_s0s;

^mn_{DD ss}⁰ _ss⁰e^Ens Z_n ;

^m<n_{KK ss}⁰ A_KK^m1. . .A_KDⁿⁿⁿ_DDA_DK^ny. . .A_KK^m1y_ss⁰: (10)

Similarly, the static quantityhB^ C^i_n equals the first line’s RHS without the function. The matrix elements ^mn_{XX ss}⁰ P

eX

mhsej^ⁿ_DDjs⁰ei^X_m are given by the second and third lines, together with^mn_KK ^m<n_DD 0. After performing a ‘‘forward run’’ to generate all relevant NRG eigenenergies and matrix elements,A^BC!can be calcu- lated in a single ‘‘backward run,’’ performing a sum with the structure PN

m>n₀C^redB ^m, starting from m N. Here ^m;red_XX PN

nmw_n^mn_XX (updated one site at a time during the backward run) is thefullreduced density matrix for shell m, obtained iteratively by tracing out all shells at smaller scalesⁿ⁼² (nm).

Equations (8)–(10) are the main results of our ‘‘FDM- NRG’’ approach. They rigorously generalize Hofstetter’s DM-NRG [8] (which leads to similar expressions, but using w_nnNT and without excluding KK matrix ele- ments), and provide a concise prescription, free from double-counting ambiguities, for how to combine NRG data from different shells when calculatingA^BC!. The relevant sum rule is satisfiedidentically, since by construc- tion R

d!A^BCn ! hB^ C^i_n holds for every n and arbi- trary temperature and NRG parametersandM_K.

Smoothing discrete data. —We obtain smooth curves for A^BC!by broadening the discretefunctions in Eq. (10) using a broadening kernel that smoothly interpolates from a log-Gaussian form (of width) [2,4] forj!j*!₀, to a regular Gaussian (of width!₀) forj!j< !₀, where!₀is a

‘‘smearing parameter’’ whose significance is explained below. To obtain high-quality data, we combine small choices of with an average over N_z slightly shifted discretizations [3] (see [13] for more details).

Application to Anderson model. —We illustrate our method for the standard single-impurity Anderson model (SIAM). Its local Hamiltonian h^₀ P

₀c^y₀c₀

Uc^y_0"c0"c^y_0#c0# describes a localized state with energy ₀, with a Coulomb penalty U for double occupancy. It is coupled to a Wilson chain P

n_nc^y_n1c_nH:c:, which generates a local level width . We calculated A^<! A^cy0c0!, A^>! A^c0cy0! and AA^>A^<. An ‘‘improved’’ version A^im thereof can be obtained by calculating the impurity self-energy

!; T [6,13] via FDM-NRG, which is less sensitive to smoothening details and yields more accurate results for the Kondo peak heightAT’00at zero temperature.

Sum rules. —As expected, we find FDM-NRG to be significantly more accurate at lower computational cost

−2 −1 0 1 2

0 1

ω/T A(ω)/A0

T=TK=3.71⋅10−4 εd=−U/3, N

z=8, α=0.3, Λ=2 (a) FDM

DM−NRG NRG FDM (smeared)

10−4 10−2 100 102 104

0 1 2

Λ=3 Λ=2

T

ω/TK

A(ω)/A0 T/5

α=0.69, N z=1 ω0/T=0.005 (b) MK=8192

MK=512

0 0.02 0.04 0.06 0.08 0.1

0.95

1 N

z=4 (not averaged) G(0)=1.0044⋅ A0 ω0=n.a.

Tfit=0.027 TK cfit=6.17 ±0.26 (4.19%) [c≡π4/16=6.0881]

T/TK

G(T)/G(0)

(c) FDM−NRG

1−c⋅(T/TK)2 1−cfit⋅(T/TK)2

0 0.1 0.2 0.3

−0.1 0

ω/TK , • (T/TK)⋅π/√3, x δAT(ω), • δA(T)

T/TK 0.069 0.047 0.032 ...

0.003 0.002 0.001 Aim

0(0)=1.0064⋅A 0 Nz=24, α=0.3, ω0/T=0.5

(d) δAT(ω)

δA(T)

−(3c/2π2)⋅x2 (Fermi Liquid)

−0.2 0 0.2 0

1 T/TK

=0.01

ω

20 30 40 50

0 0.5

Wilson shell n

wn NTΛ=3 NTΛ=2

FIG. 2 (color online). FDM-NRG results for the spectral func- tion A_T! of the SIAM, with U0:12, 0:01, _d U=2(T_K2:18510⁴),1:7, andM_K1024, unless indicated otherwise. Inset of (a): FDM-NRG result for AT! with !in units of bandwidth. For (a),(b), an unconventionally small smearing parameter was used, !00:005T [except for thick gray (red) curve in (a)], with !₀0:5T), leading to spurious low-frequency oscillations. These illustrate the differ- ences (a) between NRG (dashed green curve), DM-NRG [solid thin (blue) curve), and FDM-NRG (black curve) results for the regime!&T, and (b) between different choices ofM_K and for FDM-NRG, which yield different shapes for the weightsw_n [shown in inset of (b)]: larger reduces the scaleT at which oscillations set in, but yields less accurate values for the Kondo peak height in the regimeT&!&TK. (c),(d) Comparison of high-quality FDM-NRG data (dots, solid curves) with exact Fermi-liquid results (black dashed lines) for (c) the conductance GTforTT_K, and (d) forA^im_T !forT; !T_K. In (c),c_fit was found from a data fit to c_fitT=T_K² forT < T_fit (arrow).

In (d) we plot A_T! A^im_T! A^im_T0=A^im₀ 0vs!=T_K (curves) and AT A^im_T0=A^im₀ 0 1 vs T=TK= p3 (dots), for a set of 12 temperatures between 0.001 and 0:069T_K(with curves and dots having sameTin the same color), to illustrate the leading!andTbehavior ofA^im_T!; the dashed black line represents the expected Fermi-liquid behavior in both cases,3c=2²x²vsx.

PRL99,076402 (2007) P H Y S I C A L R E V I E W L E T T E R S week ending 17 AUGUST 2007

076402-3

than NRG or DM-NRG [8,15]. The sum rules Z d!A^cy0c0! hc^y₀c₀i_T; Z

d!A! 1 (11) hold exactly to 10¹⁵ for our discrete data, and to10⁴ after smoothing (due to numerical integration inaccura- cies). Moreover, even forM_Kas small as 256, our results forAT’00andA^imT’00typically agree to within 2% and 0.2%, respectively, with the Friedel sum rule, which re- quires A^exactT0 sin²hc^y₀c₀i₀. The exact relation A^<! f!A! (f is the Fermi function), which follows from detailed balance, is likewise satisfied well (though not rigorously so): the left-hand side of Eq. (11) typically equalsRd!f!A!to better than10⁴.

Low-frequency data. —Because of the underlying loga- rithmic discretization, all NRG-based schemes for calcu- lating finite-temperature spectral functions inevitably produce spurious oscillations at very low frequencies j!j T. The scale_T at which these set in can be under- stood as follows: the Lehmann sum in Eq. (1) is dominated by contributions from initial statesjaiwith energyE_a’T, represented by NRG shells withnnearN_T. The character- istic energy scale of these states limits the accuracy obtain- able for energy differences E_ba to accessible final states jbi. Thus the scale_T is set by those shells which contrib- ute with largest weightw_nto the density matrix.

We analyze this in more detail in Figs.2(a)and2(b)by purposefully choosing the smearing parameter to be un- conventionally small,!₀ T. The resulting spurious os- cillations are usually smeared out using !₀ *_T [Fig. 2(a), thick gray (red) curve], resulting in quantita- tively accurate spectral functions only forj!j*!₀’_T. For conventional NRG approaches, the ‘‘single-shell’’ ap- proximationw_nnNT typically leads to_T ’T, as can be seen in Fig. 2(a) [dashed (green) line and thin solid (blue) line]. In contrast, FDM-NRG yields a significantly reduced value of _T ’T=5 [Fig. 2(a), black line, and Fig. 2(b)], since the weighting functions w_n [inset of Fig.2(b)] retain weight over several shells below N_T, so that lower-frequency information is included.

Fermi-liquid relations. —To illustrate the accuracy of our low-frequency results, we calculated A^imT ! for

!; TT_K for the symmetric SIAM, and made quantita- tive comparisons to the exact Fermi-liquid relations [14],

A_T! ’A₀

1c 2

T T_K

₂ 3c

2² !

T_{K 2}

;

GT Z¹

1

d!A!; T

@f

@!

’A₀

1c T

T_{K 2}

: HereA₀1=,c⁴=16, and the Kondo temperature T_K is defined via the static magnetic susceptibility [4]

0jT0 1=4T_K. Figures 2(c) and 2(d) show the FDM- NRG data [gray (colored) dots and lines] to be in remark- ably good quantitative agreement with these relations (black dashed curves). The results for the ‘‘conductance’’

GT, being a frequency integrated quantity obtained by summing over discrete data directly without the need for broadening, are more accurate than for A^im_T !, and re- produce the prefactor c with an accuracy consistently within 5% (until now, accuracies of the order of 10%–

30% had been customary). The smoothness of the data in Fig.2(c), obtained using temperatures not confined to the logarithmic grid ⁿ⁼² [gray vertical lines in Fig. 2(b)], together with the remarkable stability with respect to dif- ferentzshifts illustrate the accuracy of our approach.

Conclusions. —Our FDM-NRG method offers a trans- parent framework for the calculation of spectral functions of quantum impurity models, with much improved accu- racy at reduced complicational cost. Its results satisfy frequency sum rules rigorously and give excellent agree- ment with other consistency checks such as the Friedel sum rule, detailed balance, or Fermi-liquid relations, including the regime!&T.

We thank F. Anders, R. Bulla, T. Costi, T. Hecht, W.

Hofstetter A. Rosch, and G. Za´rand for discussions, and the KITP in Santa Barbara for its hospitality. The work was supported by DFG (No. SFB 631 and No. De-730/3-1,3-2), and in part by the NSF (No. PHY99-07949).

Note added. —Just before completion of this work we learned that Peters, Pruschke, and Anders had followed up on the same idea [15].

[1] K. G. Wilson, Rev. Mod. Phys. 47, 773 (1975); H. R.

Krishnamurti, J. W. Wilkins, and K. G. Wilson, Phys.

Rev. B21, 1003 (1980);211044 (1980).

[2] O. Sakai, Y. Shimizu, and T. Kasuya, J. Phys. Soc. Jpn.58, 3666 (1989).

[3] M. Yoshida, M. A. Whitaker, and L. N. Oliveira, Phys.

Rev. B41, 9403 (1990).

[4] T. A. Costi, A. C. Hewson, and V. Zlatic, J. Phys. C6, 2519 (1994).

[5] T. A. Costi, Phys. Rev. B55, 3003 (1997).

[6] R. Bulla, A. C. Hewson, and T. Pruschke, J. Phys. C10, 8365 (1998).

[7] R. Bulla, T. Costi, and T. Pruschke, arXiv:cond-mat/

0701105.

[8] W. Hofstetter, Phys. Rev. Lett.85, 1508 (2000).

[9] R. Bulla, T. A. Costi, and D. Vollhardt, Phys. Rev. B64, 045103 (2001).

[10] R. Bulla, N. H. Tong, and M. Vojta, Phys. Rev. Lett.91, 170601 (2003).

[11] F. B. Anders and A. Schiller, Phys. Rev. Lett.95, 196801 (2005); Phys. Rev. B74, 245113 (2006).

[12] F. Verstraeteet al., arXiv:cond-mat/0504305.

[13] See EPAPS Document No. E-PRLTAO-99-025733 for appendices which give more details. For more information on EPAPS, see http://www.aip.org/pubservs/epaps.html.

[14] A. Hewson, The Kondo Problem to Heavy Fermions (Cambridge University Press, New York, 1993).

[15] R. Peters, T. Pruschke, and F. B. Anders, Phys. Rev. B74, 245114 (2006). We recommend this paper for a more detailed comparison of the new and previous approaches.

PRL99,076402 (2007) P H Y S I C A L R E V I E W L E T T E R S week ending 17 AUGUST 2007

076402-4