Continuous-Time Quantum Monte Carlo Method for the Coqblin-Schrieffer Model
aa r X i v : . [ c ond - m a t . s t r- e l ] A ug Typeset with jpsj2.cls < ver.1.2 > Full Paper
Continuous-Time Quantum Monte Carlo Method for the Coqblin-Schrieffer Model
Junya
Otsuki ∗ , Hiroaki Kusunose , Philipp Werner and Yoshio Kuramoto
Department of Physics, Tohoku University, Sendai 980-8578, Japan Department of Physics, Faculty of Science, Ehime University, Matsuyama 790-8577, Japan Columbia University, 538 West, 120th Street, New York, New York 10027, USA (Received November 30, 2018)
An impurity solver based on a continuous-time quantum Monte Carlo method is developed forthe Coqblin-Schrieffer model. The Monte Carlo simulation does not encounter a sign problemfor antiferromagnetic interactions, and accurately reproduces the Kondo effect. Our algorithmcan deal with an arbitrary number N of local degrees of freedom, becomes more efficient forlarger values of N , and is hence suitable for models with orbital degeneracy. The dynamicalsusceptibility and the impurity t -matrix are derived with the aid of the Pad´e approximationfor various values of N , and good agreement is found with other methods and available exactresults. We point out that the Korringa-Shiba relation needs correction for a finite value of theexchange interaction. KEYWORDS: continuous-time quantum Monte Carlo method, Coqblin-Schrieffer model, Kondo model, t -matrix, dynamical susceptibility, Pad´e approximation, Korringa-Shiba relation
1. Introduction
Strong correlations among localized and conductionelectrons lead to the Kondo effect in impurity systems,and heavy fermions in a periodic lattice. Two contrast-ing approaches may be used to deal with lattice mod-els theoretically: one is to solve the model on a finitecluster by a method such as exact diagonalization, theother involves the solution of an effective impurity sys-tem within the framework of dynamical mean field the-ory (DMFT). The former approach is more suitable forlow-dimensional systems, while the latter becomes exactin infinite-dimensional systems. For actual heavy fermionmaterials in three dimensions, the DMFT is the simplestapproach which can take local correlations into account.But even within this approximate framework, perturba-tive calculations fail in the most interesting parameterrange, where fierce competition arises between local andinter-site correlations. Therefore, numerical approachesare required to reliably solve the effective impurity prob-lem.In this paper, we present a new impurity solver basedon the continuous-time quantum Monte Carlo method(CT-QMC) for fermion systems, which has been origi-nally proposed by Rubtsov et al. The CT-QMC eval-uates the infinite sum of multiple integrals in a per-turbation expansion by means of a Monte Carlo proce-dure. The CT-QMC does not require a Trotter decom-position in contrast with auxiliary field quantum MonteCarlo methods such as Hirsch-Fye. While the CT-QMCwas first formulated as a perturbation expansion withrespect to the Coulomb interaction, an alternative ex-pansion with respect to hybridization around the atomiclimit has been developed. For the Anderson model, thelatter approach does not encounter a sign problem andempirically, it is found that even more complicated mul-tiorbital and cluster models can be simulated withoutencountering negative weight configurations. Therefore ∗ E-mail address: [email protected] it is a powerful impurity solver for DMFT and its exten-sions, and has been applied to several models.
3, 5–8
The CT-QMC solvers
2, 3 are not suitable for verystrong correlations. In the “segment” picture of ref. 3, thesimulation of the Anderson model with strong Coulombrepulsion and a deep local level involves short excur-sions to high-energy states (short overlapping segmentsor anti-segments). Without special precautions, the ac-ceptance rate of the Monte Carlo moves will be very low.However, the empty and doubly occupied states play anessential role and give rise to the effective exchange in-teraction. Since these high-energy states are easily dealtwith by ordinary perturbation theory, it is reasonableto employ an effective model which eliminates these vir-tual processes from the start. It is well known that theCoqblin-Schrieffer (CS) model corresponds to such a lo-calized limit of the Anderson model. The purpose of thispaper is to formulate the CT-QMC in a form which isapplicable to the CS model.Our scheme, to be presented below, has the followingfeatures:(i) As long as interactions are antiferromagnetic, thescheme is free from the minus sign problem. Thisis consistent with the fact that the CS model withthe antiferromagnetic exchange is derived from theAnderson model.(ii) Acceptance probabilities for the random walk arehigher compared to the case of the Anderson model.This is because the formulation excludes the chargedegree of freedom. The remaining magnetic ex-change processes are sampled efficiently.(iii) An arbitrary number N of local states can be dealtwith in a simple manner. Larger values of N makethe computation faster in general, since the numberof operators in each component decreases. Hence thescheme can easily handle orbital degeneracy.Since the algorithm does not encounter the minus signproblem, highly accurate dynamics can be derived. The Full Paper J. Otsuki et al. algorithm is therefore especially useful for large- N sys-tems, where other methods such as exact diagonalizationare not practical due to the rapidly increasing number ofbasis states.This paper is organized as follows. In the next section,the CT-QMC formalism is presented for the CS model.Section 3 discusses the Monte Carlo sampling procedure.The algorithm is also applied to the Kondo model whichis discussed in §
4. Numerical results for static quantitiesare shown in § N . Single-particleand two-particle spectra are evaluated with the aid of thePad´e approximation in §
6. The summary of the paper isgiven in §
2. Formulation
In order to present our new algorithm, it is conve-nient to summarize the general formulation of the CT-QMC.
2, 3
The partition function for the Hamiltonian H = H + H is given by Z = Tr ( T τ e − βH exp " − Z β d τ H ( τ ) , (1)where H ( τ ) = e βH H e − βH is an operator in the inter-action picture. Expanding the exponent we obtain an-other expression for the partition function: Z = ∞ X k =0 ( − k k ! Z β d τ · · · Z β d τ k × Tr { T τ e − βH H ( τ ) · · · H ( τ k ) } . (2)In the CT-QMC, this sum of multiple integrals is sam-pled stochastically.We consider the Coqblin-Schrieffer model with N com-ponents H = H c + H f = X k α ǫ k c † k α c k α + X α ( E α + J αα ) X αα ,H = X αα ′ J αα ′ X αα ′ ( − c α c † α ′ ) , (3)where ǫ k is the energy of conduction electrons with re-spect to the chemical potential, and X αα ′ is the X -operator on the local states | α i defined by X αα ′ = | α ih α ′ | . c α = N − / P k c k α , with N being number ofsites, is the annihilation operator for the conduction elec-tron at the impurity site. The order of the conductionoperators in H specifies the definition of the equal-timeGreen function used in the perturbation expansion. Inorder to distinguish the impurity contribution to the par-tition function, we factorize Z as Z = Z f Z c , where Z c isthe partition function without the impurity. The impu-rity contribution Z f is given by Z f = ∞ X k =0 ( − k Z β d τ · · · Z βτ k − d τ k X α α ′ · · · X α k α ′ k × J α α ′ · · · J α k α ′ k × s Y α h T τ c † α ( τ ′ ) c α ( τ ′′ ) · · · c † α ( τ ′ k α ) c α ( τ ′′ k α ) i c × Tr f { T τ e − βH f X α α ′ ( τ ) · · · X α k α ′ k ( τ k ) } , (4)where conduction electron operators are grouped bycomponent index, and averaged separately by h· · · i c = Z − Tr c { e − βH c · · · } . A resultant sign in the permutationis represented by s . k α is the number of operators c † α c α for each component α , and P α k α = k . Using Wick’stheorem for the conduction electrons, we obtain Z f = Z D[ k ] W k ,W k =( − k J α α ′ · · · J α k α ′ k · s Y α det D ( k α ) α × Tr f { T τ e − βH f X α α ′ ( τ ) · · · X α k α ′ k ( τ k ) } , (5)where R D[ k ] denotes the sum of k and { α i } and themultiple integrals over { τ i } . The k α by k α matrix D ( k α ) α is defined by ( D ( k α ) α ) ij = g α ( τ ′′ i − τ ′ j ) with g α ( τ ) = −h T τ c α ( τ ) c † α i c .In the CT-QMC, statistical averages are evaluated bymeans of, for example, the Metropolis algorithm with theweight of Monte Carlo configurations given by eq. (5).The implementation of the random walk will be dis-cussed in §
3. During the simulations, it is enough to storeonly the inverse matrix M ( k α ) α = ( D ( k α ) α ) − . The matrix M ( k α ) α is utilized to evaluate determinant ratios and isefficiently updated using fast-update formulae. t -matrix The conduction electron Green function G α ( τ, τ ′ ; { τ i } )for each configuration { τ i } is represented as follows: G α ( τ, τ ′ ; { τ i } ) = g α ( τ − τ ′ ) − X ij g α ( τ − τ j )( M ( k α ) α ) ji g α ( τ i − τ ′ ) . (6)This equation can be derived by using the fast-updateformula for M ( k α ) α . An average over the Monte Carloensemble gives the physical Green function: G ( τ, τ ′ ) = h G ( τ, τ ′ ; { τ i } ) i MC . After the Fourier transform, we ob-tain G α (i ǫ n ) = g α (i ǫ n ) + g α (i ǫ n ) t α (i ǫ n ) g α (i ǫ n ) , (7) t α (i ǫ n ) = − T *X ij ( M ( k α ) α ) ji e i ǫ n ( τ j − τ i ) + MC , (8)where ǫ n = (2 n + 1) πT is the fermion Matsubara fre-quency. Since numerical summations over i and j aretime consuming, it is more convenient to measure inimaginary time as follows: t α ( τ ) = − T *X ij ( M ( k α ) α ) ji δ ( τ, τ j − τ i ) + MC . (9)In numerical calculations, the δ -function is replaced by arectangular function with a finite width, and τ is sampledin τ > . Phys. Soc. Jpn. Full Paper J. Otsuki et al. 3
We should note that the t -matrix in frequency spaceincludes a constant term t (1) α in the first Born approxi-mation t (1) α = J αα h X αα i . (10)Although M ( k α ) α contains information of t (1) α at τ j = τ i ,it is convenient to add it after the Fourier transformationrather than measure it in τ -space. The constant is rele-vant for the calculation of the Kondo model as will bediscussed in §
4. If the t -matrix is sampled in frequencyspace, the constant term is automatically included in t α (i ǫ n ). We consider two-particle correlation functions χ αα ′ ( τ )defined by χ αα ′ ( τ ) = h T τ ˜ X H αα ( τ ) ˜ X α ′ α ′ i , where the super-script H indicates the Heisenberg operator, and the tildeindicates deviation from the expectation value: ˜ X αα = X αα − h X αα i . In the CT-QMC, h T τ X H αα ( τ ′ ) X H α ′ α ′ ( τ ′′ ) i isevaluated by * Tr f { T τ e − βH f X α α ′ ( τ ) · · · X α k α ′ k ( τ k ) X αα ( τ ′ ) X α ′ α ′ ( τ ′′ ) } Tr f { T τ e − βH f X α α ′ ( τ ) · · · X α k α ′ k ( τ k ) } + MC . (11)In actual computations, it is not necessary to evaluatethe matrix products for the trace. Instead, it is suffi-cient to judge whether X αα ( τ ′ ) and X α ′ α ′ ( τ ′′ ) are onsegments of α and α ′ components respectively. In otherwords, we test whether X αα ( τ ′ ) and X α ′ α ′ ( τ ′′ ) are per-mitted by the X -operators in front and behind of them.By averaging over τ ′ with τ = τ ′ − τ ′′ , χ αα ′ ( τ ) is ob-tained with high accuracy. The equal-time correlationcan be represented in terms of the mean occupation by χ αα ′ (0) = h X αα i ( δ αα ′ − h X α ′ α ′ i ).We also consider response functions for M = P α m α X αα χ ( τ ) = h T τ M H ( τ ) M i . (12)By choosing proper m α with P α m α = 0, we can dealwith magnetic, quadrupole and other moments within N states. Provided the local levels are degenerate andthe system has SU( N ) symmetry, the susceptibility isgiven by the Curie constant C N = N − P α m α . Thedynamical susceptibility χ ( τ ) is then given in terms of χ αα ′ ( τ ) by χ ( τ ) C N = X α χ αα ( τ ) − N − X α = α ′ χ αα ′ ( τ ) . (13)The sums over α and α ′ improve the accuracy of sam-pling, although χ αα and χ αα ′ are identical for all α andfor any combination of α = α ′ , respectively. The equal-time correlation is χ (0) /C N = 1. The static susceptibil-ity is evaluated by integrating χ ( τ ). Although the staticsusceptibility can also be obtained by measuring the ex-pectation value of M in the presence of a small externalfield, the evaluation using χ ( τ ) gives results with higherprecision. Thermodynamic quantities can be evaluated from thesingle-particle Green function. The expression for the in-ternal energy is obtained from the equation of motion for G ( τ, τ ′ ) in the limit τ ′ → τ + 0. A contribution E imp of the impurity to the internal energy is given in termsof the impurity t -matrix t α (i ǫ n ) for forward scattering asfollows: E imp = h H i − h H c i c = X α " E α h X αα i + T X n i ǫ n N X k g k α (i ǫ n ) ! t α (i ǫ n )e i ǫ n δ , (14)where δ is a positive infinitesimal quantity.The specific heat C is evaluated from the difference of E imp at different temperatures. Assume that the inter-nal energies at temperatures T and T ( T < T ) areobtained as E and E , respectively. The specific heat at( T + T ) / C = E − E T − T , (15)and its statistical error ∆ C is estimated by∆ C = p (∆ E ) + (∆ E ) T − T , (16)where ∆ E and ∆ E are standard deviations of E and E , respectively. In this derivation, Gaussian distribu-tions have been assumed. We note that the specific heatcomputed in eq. (15) includes, in addition to the statis-tical errors, an error due to the finite differences, whichis proportional to powers of T − T . Hence, in orderto obtain a reasonable numerical accuracy, the statisti-cal errors of the internal energy need to be smaller than E − E when the temperature difference T − T is de-creased.
3. Monte Carlo Procedure
In the CT-QMC method, we evaluate the statisticalaverage of physical quantities by samplings with respectto the weight W k in eq. (5). The random walk in con-figuration space must satisfy the ergodicity and detailedbalance conditions.Updates which change the order of J are required forergodicity, and updates which shift one of the operatorsincrease sampling efficiency. In this section, we introducetwo different algorithms for the random walk. The first,which changes the perturbation order by ± N = 2 model, the perturbation order must be changed by ±
2, because only even powers of J contribute to the par-tition sum. In the general case, if only some coupling con-stants are finite, one needs to manipulate several opera-tors (up to N ) in one update. The algorithm is presentedin the latter of this section. Either algorithm should bechosen depends on the model. J. Phys. Soc. Jpn.
Full Paper J. Otsuki et al. b t k t i a i a k a k a i - a t ......... ...... Fig. 1. A diagram representing the configuration of { τ i } and { α i } of order J k . t n t n + a n a n + t n t n + a n a n + a t(cid:223)(cid:13) Fig. 2. Illustration of an insertion of a segment.
A certain configuration of order J k is represented interms of { τ i } = ( τ , · · · , τ k ) and { α i } = ( α , · · · , α k ).These variables define the sequence of X -operators X α k α k − ( τ k ) · · · X α i α i − ( τ i ) · · · X α α k ( τ ) . (17)The corresponding c -operators are given by( − k +1 c † α k ( τ ) c α k ( τ k ) · · · c † α i ( τ i +1 ) c α i ( τ i ) × · · · c † α ( τ ) c α ( τ ) . (18)We represent the above configuration by a diagram asshown in Fig. 1. The most efficient update for { τ i } and { α i } is the addition or removal of a single element. Weconsider the process of adding τ and α , which are ran-domly chosen in the range [0 , β ) and from the N compo-nents, respectively. If τ satisfies τ n +1 > τ > τ n , { τ i } and { α i } change into ( τ , · · · , τ n , τ, τ n +1 , · · · , τ k ) and( α , · · · , α n , α, α n +1 , · · · , α k ), respectively. Then one ofthe X -operators is altered as X α n +1 α n ( τ n +1 ) → X α n +1 α ( τ n +1 ) X αα n ( τ ) , (19)which corresponds to the change illustrated in Fig. 2.Namely, the segment α is inserted between α n and α n +1 with shortening of the segment α n . In the correspondingremoval process, we erase one randomly chosen segment.According to the detailed balance condition, the ratioof the transition probabilities should be p ( k → k + 1) p ( k + 1 → k ) = W k +1 W k N βk + 1 . (20)The factors N and β are due to the random choices of α and τ , respectively, and k + 1 due to that in the removalprocess. Since W k has the dimension of J k , W k +1 β/W k isa dimensionless quantity. For k = 0, the ratio W k +1 /W k is given by W k +1 W k = J α n +1 α J αα n J α n +1 α n exp[ − l ( E α − E α n )] × det D (+) α det D α det ˜ D α n det D α n , (21) where l = τ n +1 − τ is the length of the new segment. D (+) α is the matrix with c † α ( τ n +1 ) c α ( τ ) added to D α , and ˜ D α n is the matrix with one of the operators shifted in timeaccording to c † α n ( τ n +1 ) → c † α n ( τ ). The ratio of determi-nants can be evaluated using fast-update formulae. If α = α n in Fig. 2, the change is just an addition of adiagonal element X αα ( τ ), so that eq. (21) is reduced to W k +1 W k = − J αα det D (+) α det D α . (22)Here D (+) α is a matrix in which c † α ( τ ) c α ( τ +0) is added tothe original one. The equal-time Green function in D (+) α should be g α (+0) to keep the probability positive. Thisis verified by taking the CS limit in the correspondingformula of the Anderson model (see Appendix A). In thecase of k = 0, we should average over the local statessince all states contribute to the trace. Then W /W isgiven by W W = − J αα ρ α g α (+0) , (23)where ρ α is defined by ρ α =exp( − βE α ) / P α ′ exp( − βE α ′ ).The ratios of the weights in eqs. (21)–(23) changetheir signs depending on the signs of the coupling con-stants. We have confirmed by numerical calculations thatthe probability remains positive in the case of antiferro-couplings, i.e., J αα ′ >
0. This is consistent with the factthat the CS model with antiferro-couplings is derivedfrom the Anderson model, where the minus sign problemdoes not appear. We also note that we could consider an-other process in which α is inserted between α n − and α n . It corresponds to a diagram where α and α n are ex-changed in Fig. 2. However, the ergodicity condition canbe satisfied without this process. If some coupling constants are zero, the ergodicity mayrequire updates that change the diagrams by several per-turbation orders. At most N operators should be insertedand removed in one update. To this end, we define irre-ducible sets of operators, which consist of κ operatorswith a set of κ distinct components ( α ′ , · · · , α ′ κ ). Thenthe irreducible sets of X -operators are given by X α ′ κ α ′ κ − ( τ ′ κ ) · · · X α ′ j α ′ j − ( τ ′ j ) · · · X α ′ α ′ κ ( τ ′ ) . (24)We consider a process of insertion and removal of theirreducible operator, which is illustrated in Fig. 3.To determine an irreducible operator, we begin by ran-domly choosing κ from 1 to N and τ ′ from the range[0 , β ). If τ ′ satisfies τ n +1 > τ ′ > τ n , α ′ κ is fixed as α ′ κ = α n . The remaining components α ′ j (1 ≤ j < κ )are determined so as to be different from each other.The imaginary times τ ′ j (1 < j ≤ κ ) are randomly cho-sen from the restricted range ( τ ′ j − , τ n +1 ). In the case ofa removal process, we choose one, say α n , from all thesegments. If two identical components appear betweenthe selected segment and the next segment with compo-nent α n , the operator set is not irreducible and hence theupdate is rejected at this stage. If this is not the case, weproceed to attempt a removal of the operators. . Phys. Soc. Jpn. Full Paper J. Otsuki et al. 5 t n t n + a n a n + t n t n +1 a n a n + a ' t ' (cid:223)(cid:13)a ' k t ' k ...... ( =a n ) Fig. 3. Illustration of an addition of an irreducible operator.
From the detailed balance condition, we obtain theupdate probability of the above process for k = 0 as p ( k → k + κ ) p ( k + κ → k ) = W k + κ W k N !( N − κ )! 1 k + κ κ Y j =1 l (max) j , (25)where l (max) j is given by l (max)1 = β and l (max) j = τ n +1 − τ ′ j − for j >
1. The factor N ! / ( N − κ )! is the inverse of theprobability that the corresponding irreducible operator ischosen. The ratio W k + κ /W k is given by W k + κ W k = − κ Y j =1 J α ′ j α ′ j − exp[ − l j ( E α ′ j − E α n )] det D (+) α ′ j det D α ′ j , (26)where l j is a length of the segment α ′ j . D (+) α ′ j is a matrix towhich c † α ′ j ( τ ′ j +1 ) c α ′ j ( τ ′ j ) has been added, with τ ′ κ +1 = τ ′ .In the case of k = 0, the above procedure is applicablewith τ n = 0 and τ n +1 = β . One then has to choose α ′ κ at random, and take the trace over the local states. Thetransition probability for k = 0 becomes p (0 → κ ) p ( κ →
0) = W κ W N · N !( N − κ )! κ Y j =1 l (max) j . (27)The additional factor N originates in the random choiceof α ′ κ . We note that the factor 1 / ( k + κ ) in eq. (26) isnot necessary for k = 0 because no choice is made in theremoval process. The ratio W κ /W is given by W κ W = − Z − f κ Y j =1 n J α ′ j α ′ j − exp( − l j E α ′ j ) det D (+) α ′ j o , (28)where we define Z f by Z f = P α exp( − βE α ).
4. The Kondo Model
When we deal with S = 1 / H = X k σ ǫ k c † k σ c k σ + J S · σ c , (29) where σ c is the Pauli matrix for conduction electrons.The Hamiltonian is rewritten in terms of the CS inter-action as follows: H = X k σ ǫ k c † k σ c k σ + v X σ c † σ c σ + J X σσ ′ f † σ f σ ′ c † σ ′ c σ , (30)where v = − J/
2. Therefore the algorithm for the CSmodel is applicable by including the potential scatteringterm into H .We introduce the Green function including the poten-tial scattering, ˜ g ( z ), which is a scalar in the spin indices.In terms of the bare Green function g ( z ), ˜ g ( z ) is ex-pressed as ˜ g = g − vg . (31)In the simulation of the CS model, g ( z ) is replaced by˜ g ( z ). Within the CT-QMC, the impurity t -matrix t J ( z )is computed with respect to ˜ g ( z ). To obtain the t -matrix t ( z ) of the Kondo model, eq. (29), the contribution ofthe potential scattering should be subtracted from ˜ g ( z ).The full Green function G ( z ) can be expressed as G = ˜ g + ˜ gt J ˜ g = g + gtg. (32)Solving eq. (32) with respect to t ( z ), we obtain t = v − vg + t J (1 − vg ) , (33)where the first term is the t -matrix due to the poten-tial scattering. Concerning the two-particle correlationfunction, no modification is required, because it does notdepend on the selection of the basis set in the perturba-tion expansion.
5. Imaginary-Time Data and Static Quantities
We apply our algorithm to a model with SU( N ) sym-metry. In this case, both types of updates introduced in § − π Im g α ( ǫ + i δ ) = ρ θ ( D − | ǫ | ) , (34)where θ ( ǫ ) is a step function and ρ = 1 / D . In theMatsubara formalism, the Green function is representedas g (i ǫ n ) = − ρ arctan( D/ǫ n ) and N − P k g k (i ǫ n ) = − ( ǫ n + D ) − . We choose the unit D = 1 in this pa-per. The standard deviations in the MC ensembles areevaluated from 20 averages. In MC simulations, we haveobserved minus sign probabilities at temperatures lowerthan the Kondo temperature. However, they appear onlyat the rate of about one in 10 updates, and may be dueto rounding errors.In the following, we show numerical results for physi-cal quantities without analytic continuation, such as thestatic susceptibility and the specific heat. These resultsare, within error bars, exact. Spectral function in realfrequencies will be shown in the next section. J. Phys. Soc. Jpn.
Full Paper J. Otsuki et al. P ( k ) k (a) J = 0.20.40.60.81.0-0.06-0.05-0.04-0.03-0.02-0.01 0 0 5 10 15 20 25 30 35 40 45 50 t ( i e n ) – t ( ) n (b) J = 0.2ReImexact Fig. 4. (a) Probability distribution P ( k ) for terms of order J k ,and (b) t -matrix t (i ǫ n ) for N = 1. The temperature is chosen as T = 0 . N = 1 : potential scattering The Hamiltonian, eq. (3), is trivially solvable when N = 1, which corresponds to potential scattering. Theexact t -matrix for the potential scattering is given by thefirst term on the right-hand side of eq. (33) with v = J ,which takes into account an infinite sequence of scat-tering events. We apply our algorithm to the potentialscattering to compare with the exact solution, prior toapplications to N >
1. Figure 4(a) shows the probabil-ity of appearance of perturbation terms of order J k inthe Monte Carlo simulations. Although the exact resultsinclude infinite series of scattering, only terms of finiteorder of J are significant in Monte Carlo simulations.The center of the distribution increases as J increases.Results for the t -matrix in the Monte Carlo ensemble av-erage are shown in Fig. 4(b) together with the analyticalresult. The term t (1) given by eq. (10) is subtracted inthis figure. The exact results are completely reproduced,which demonstrates the convergence of the perturbationseries, and shows that sampling finite perturbation or-ders is sufficient. N = 8 : large- N case The present algorithm can handle arbitrary N . We firsttake N = 8 and give exemplary results of static quanti-ties as well as raw data in the imaginary-time domain.Figure 5 shows the dynamical susceptibility χ ( τ ) for sev-eral values of T . We have plotted with lines because the c ( t ) / C N t T = 0.0010.0050.010 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 t / b Fig. 5. Dynamical susceptibility χ ( τ ) in the imaginary-time do-main for N = 8 and J = 0 . c ( T ) / C N TJ = 0.0700.0750.0800.0900.100 CT-QMCNCA Fig. 6. Temperature dependence of the static susceptibility χ for N = 8 and several values of J . Dashed lines are results computedin the NCA. intervals between data points are fine enough, and theassociated errors are negligible. It turns out that the re-duction of the correlation, which starts at χ ( τ = 0) = 1,becomes more rapid for lower temperatures. This im-plies Kondo screening at low temperatures. The insetshows χ ( τ ) against τ /β . At T = 0 . τ = β/
2, while correlations remainfor T = 0 . T = 0 . χ ( ω = 0) can be evaluatedby integrating χ ( τ ). Figure 6 shows the temperature de-pendence of the static susceptibility for several values of J for N = 8. For comparison, we plot results computedin the non-crossing approximation (NCA), which incor-porates terms up to the next-leading order in the 1 /N expansion by integral equations.
11, 12
The CT-QMC re-sults are in excellent agreement with the NCA at hightemperatures. On the other hand, deviations are visibleat lower temperatures. This is due to the inaccuracies ofthe NCA at low temperatures and low frequencies. Wehave confirmed in fact that the incorrect upturn appearsin the NCA results for smaller N , where the accuracy ofthe 1 /N expansion diminishes. On the other hand, the . Phys. Soc. Jpn. Full Paper J. Otsuki et al. 7 – t ( t ) t T = 0.0010.0050.010 0 0.01 0.016 0 0.2 0.4 0.6 0.8 1 t / b Fig. 7. The impurity t -matrix t α ( τ ) in the imaginary-time do-main for N = 8 and J = 0 . CT-QMC produces proper values of the static suscepti-bility at all temperatures.We next show results for the impurity t -matrix. Fig-ure 7 shows the t -matrix t ( τ ) in the imaginary-time do-main for the same parameters as in Fig. 5. Statisticalerrors are so small that we have plotted with lines asin the case of χ ( τ ). As temperature decreases, − t ( τ ) in-creases. Correspondingly, the discontinuity at the bound-ary t ( − − t (+0) increases as temperature decreases.Since the discontinuity is given by − π − R Im t ( ω )d ω , theincrease indicates enhanced scattering due to the impu-rity. This corresponds to the formation of the Kondo sin-glet.The impurity t -matrix with imaginary frequency de-scribes thermodynamic quantities as shown in eq. (14).The temperature dependence of the internal energy E ( T )is shown in Fig. 8(a). The statistical errors are invisiblein this figure. The data give the specific heat C ( T ) viaeq. (15). Figure 8(b) shows the temperature dependenceof C ( T ) plotted together with results from the NCA. Statistical errors become larger at low temperatures dueto the small mesh of temperatures. The CT-QMC agreeswith the NCA within error bars, and reproduces a peakdue to the Kondo effect. N So far we have presented imaginary-time data and re-sultant static quantities for N = 8. Next we examinedifferent values of N . Coupling constants are chosen soas to fix N J to yield almost equal values of the Kondotemperature. We define the Kondo temperature T K bymeans of the static susceptibility χ by T K = C N /χ atlow enough temperature. Table I shows χ at T = 0 . T K for N = 2, 4, 6 and 8. For all N , T K is less than 1/30 of D , and therefore can be consid-ered small compared to D . Figure 9 shows probabilities P ( k α ) of appearance of k α for different N at T = 0 . D α is equal to k α . It turns out thatthe peak of P ( k α ) shifts to higher values of k α for smaller N . This is because the power of J is divided into N com-ponents, resulting in the decrease of k α . Consequently, alarger value of N reduces the computational burden, and -0.065-0.06-0.055-0.05-0.045-0.04-0.035-0.03-0.025-0.02-0.015-0.01 0.001 0.01 0.1 1 E ( T ) T (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.001 0.01 0.1 1 C ( T ) T (b) CT-QMCNCA Fig. 8. Temperature dependence of the internal energy E ( T ) andthe specific heat C ( T ) for N = 8 and J = 0 . χ and the Kondo temperatures T K estimated by T K = C N /χ at T = 0 . N J χ ( T = 0 . /C N T K makes it possible to reach temperatures much lower than T K .For static quantities, exact solutions have been ob-tained based on the Bethe ansatz. Temperature depen-dences of all physical quantities are given through a sin-gle energy scale. Hence we can compare our results withthe exact solutions, using the T K determined above. Wenote that the characteristic energy T in ref. 15 relatesto our Kondo temperature T K by T K = 2 πT /N . Fig-ure 10(a) shows the static susceptibility as a function of T /T K for the parameters listed in Table. I. Dashed linesare the Bethe ansatz solutions for N = 2 and 8. Compar-ison between our results and the exact solution revealsthat T K determined by T K = C N /χ systematically de-viates from the exact results. This is ascribed to the fi-nite cutoff of the conduction band in our model, which isshown to enhance the static susceptibility. Consequently, T K determined by χ tends to be smaller than the Bethe J. Phys. Soc. Jpn.
Full Paper J. Otsuki et al. P ( k a ) k a N = 8 N = 6 N = 4 N = 2 Fig. 9. Probability distributions P ( k α ) at T = 0 . T K are listed in Table I. c ( T ) T K / C N T / T K (a) N = 8 N = 2 CT-QMCBA 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.01 0.1 1 10 100 C ( T ) T / T K (b) N = 864 N = 2 CT-QMCBA Fig. 10. Temperature dependence of (a) the static susceptibilityand (b) the specific heat for N = 2, 4, 6 and 8. Parametersand T K are listed in Table I. Dashed lines are the Bethe ansatzsolution. ansatz value. The effect of the finite band width willbe discussed in detail later. Temperature dependencesof the specific heat are shown in Fig. 10(b). Large peaksaround T /T K ∼ . T /T K ∼
10 originate from the cutoffof the conduction band. We recognize systematic devia-tions between our results and the Bethe ansatz solutions,corresponding to those found in the static susceptibility.However, the overall behavior agrees well if we scale tem- -0.7-0.6-0.5-0.4-0.3-0.2-0.1 0 0 5 10 15 20 25 30 35 40 45 50 t ( i e n ) – t ( ) n Kondo model (Re)(Im)CS model (Re)(Im)
Fig. 11. t -matrix t (i ǫ n ) of the Kondo model and the CS modelfor J = 0 . T = 0 . perature independently of the static susceptibility. We can deal with the Kondo model in the way pre-sented in §
4. In applying our algorithm for the CS modelto the Kondo model, we use unperturbed Green func-tions without particle-hole symmetry, although the orig-inal model is particle-hole symmetric. Therefore it is astrict check of our algorithm whether the particle-holesymmetry is recovered after solving the CS interactions.Figure 11 shows the t -matrix t (i ǫ n ) for J = 0 . T = 0 . t α (i ǫ n ) = 0indicates symmetry with respect to particle-hole excita-tions. It is verified that the result for the Kondo modelkeeps the particle-hole symmetry, while the CS modeldoes not due to the potential scattering.
6. Dynamical Quantities with Analytic Contin-uation
We proceed to spectral functions in real frequencies.To perform analytic continuations, we employ two kindsof conventional methods: the Pad´e approximation andthe maximum entropy method (MEM). The Pad´e ap-proximation fits data at the Matsubara frequencies inthe upper-half plane with use of a rational function, andextrapolates onto the real axis. Since the Fourier trans-form drives away noises in imaginary-time data to highfrequencies, reasonable accuracy is expected at low fre-quencies. However the extrapolation is sensitive to sta-tistical and numerical errors, especially at high energieswhich are far from the imaginary axis. If the data includesignificant errors, the Pad´e approximation does not workin general. In this case, the statistical errors can be takeninto account by the MEM. The MEM derives spectralfunctions from data in the imaginary-time representa-tion without Fourier transformation. In this subsection,we shall show results by the Pad´e approximation, sinceour data have been obtained with sufficient accuracy. . Phys. Soc. Jpn.
Full Paper J. Otsuki et al. 9 – I m t ( w ) w (a) N = 2 Pade app.NRG ( T = 0) 0 0.1 0.2 0.3 0.4 0.5 0.6-0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 – I m t ( w ) w (b) N = 4 Pade app. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7-0.3 -0.2 -0.1 0 0.1 0.2 0.3 – I m t ( w ) w (c) N = 6 Pade app.NCA 0 0.1 0.2 0.3 0.4 0.5 0.6-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 – I m t ( w ) w (d) N = 8 Pade app.NCA Fig. 12. The t -matrix − Im t ( ω +i δ ) computed in the Pad´e approx-imation at T = 0 .
001 and comparison with the NRG ( N = 2)and the NCA ( N = 6 and 8). Parameters are listed in Table I.Marks at ω = 0 indicate exact values obtained from the Friedelsum rule, eq. (35). t -matrix Figure 12 shows − Im t ( ω ) computed from the Pad´e ap-proximation at T = 0 .
001 for the parameters listed inTable I. We plot N components separately to check theaccuracy of the analytic continuations. We find an excel-lent agreement between all components around the Fermi – I m t ( w ) w Pade app.NRG ( T = 0) Fig. 13. The t -matrix − Im t ( ω +i δ ) of the Kondo model computedin the Pad´e approximation. Parameters are J = 0 . T =0 . T = 0. level, while at higher energies the spectra show some dif-ferences. The extent of the difference can be regarded asa measure of the error in the approximation. For com-parison, we plot results computed using the numericalrenormalization group (NRG)
18, 19 at T = 0 for N = 2,and results from the NCA for N = 6 and 8. It turns outthat for all N the overall shapes obtained by the differentmethods are in good agreement. Results for the Kondomodel are shown in Fig. 13. In computing the spectrum,we have set the real part of t (i ǫ n ) to 0 neglecting tinystatistical errors. We have confirmed that the spectralfunction is symmetic with respect to ω = 0 and agreeswith the result from NRG.At T = 0 the Friedel sum rule relates the t -matrix atthe Fermi level with the occupation number of the localstate. For the constant density of states represented ineq. (34), − Im t α (0 + i δ ) is given by − Im t α (0 + i δ ) = 1 πρ sin ( πn α ) , (35)where n α = h X αα i is the mean occupation of the state α . The exact values of − Im t (0 + i δ ) have been marked inFigs. 12 and 13. We remark that our results computed inthe Pad´e approximation agree with the Friedel sum rulefor all N . We conclude that analytic continuation by thePad´e approximation works for the CT-QMC data. In thecase where there is an energy gap or several excitationsin the spectrum, however, the Pad´e approximation maynot provide this level of accuracy. We next present spectra of two-particle response func-tions. Figure 14 shows Im χ ( ω + i δ ) /ω computed inthe Pad´e approximation for the same parameters as inFig. 12. Since a comparison between N components,which is made in Fig. 12, is not available in this case,we have confirmed reproducibility of the spectra betweendifferent Monte Carlo ensembles. Spectra for N = 2 arealmost Lorentzian for all plotted temperatures, while for N ≥ t -matrix in Fig. 12,where the Kondo resonance shifts to higher energy with Full Paper J. Otsuki et al. I m c ( w ) / w C N w (a) N = 2 T = 0.0010.0050.010 0 100 200 300 400 500 600 700 800 900 0 0.02 0.04 0.06 0.08 0.1 0.12 I m c ( w ) / w C N w (b) N = 4 T = 0.0010.0050.010 0 200 400 600 800 1000 1200 0 0.02 0.04 0.06 0.08 0.1 0.12 I m c ( w ) / w C N w (c) N = 6 T = 0.0010.0050.010 CT-QMCNCA 0 200 400 600 800 1000 1200 0 0.02 0.04 0.06 0.08 0.1 0.12 I m c ( w ) / w C N w (d) N = 8 T = 0.0010.005 0.010 CT-QMCNCA Fig. 14. Im χ ( ω ) /ω for several temperatures. Parameters arelisted in Table I. Marks at ω = 0 indicate values deduced from theKorringa-Shiba relation, eq. (36), with a = 1 and χ ( T = 0 . N = 6 and 8, results from the NCA are plotted by a dashedline. increasing N . We also plot, for reference, results com-puted in the NCA for N = 6 and 8. In Figs. 14(c) and(d), we can see excellent agreement between two resultsat T = 0 .
01, which is of the order of T K . At lower tem-peratures, on the other hand, both results do not agreearound ω = 0 due to the unphysical peak of the NCA. a N J r N = 28 Fig. 15. The parameter a in eq. (36) as a function of NJρ for N = 2 and 8. In the Fermi liquid regime, the Korringa-Shiba relation(KSR) connects the imaginary part of the dynamical sus-ceptibility with the static ones.
20, 21
It has been provenfor the orbitally degenerate Anderson model in the wideband limit. In models with finite band width, however,the static susceptibility is enhanced and does not satisfythe KSR. This is due to the fact that the cutoff of the con-duction band enhances the static susceptibility over theuniversal value, while it does not affect the low-energyimaginary part. Appendix B demonstrates the deviationfor the non-interacting Anderson model. To account fordeviations from the wide band limit, we introduce a pa-rameter a and modify the original KSR aslim ω → a χ Im χ ( ω + i δ ) ω = πN C N . (36)Here a = 1 corresponds to the KSR, but a is greater thanunity in the case of finite band width. We have marked,in Fig. 14, values of Im χ ( ω + i δ ) /ω | ω =0 deduced from χ at T = 0 .
001 and a = 1 in eq. (36). Comparing withspectra at T = 0 . N . For N = 2, χ with a = 1 leads to about 50%deviation from the actual value of Im χ ( ω + i δ ) /ω | ω =0 .Since the width of the Kondo resonance is not negligibleas compared with D especially for N = 2, as shown inFig. 12, the effect of the finite cutoff cannot be neglected.To see the validity of the KSR in the wide band limit,or equivalently, in the weak coupling limit, we plot theparameter a against N Jρ in Fig. 15. In this figure, errorbars have been estimated from the statistical errors of χ .We note that the errors actually come from inaccuraciesof Im χ ( ω + i δ ) as well due to the analytic continuation.For N = 8, the parameter a deviates from unity only byabout 5%. For N = 2, on the other hand, a varies lin-early against N Jρ . The deviation from unity is almostproportional to N Jρ in the range of our Monte Carlosimulations. Consequently, in order to satisty the KSR, N Jρ should be much smaller than unity for N = 2.
7. Summary
We have extended the CT-QMC method for fermionsto the Coqblin-Schrieffer model. An arbitrary number N of local states can be treated by our algorithm and for . Phys. Soc. Jpn. Full Paper J. Otsuki et al. 11 antiferromagnetic interactions, the scheme is free fromthe minus sign problem. Therefore it is a powerful toolto study heavy fermion systems within the framework ofDMFT.We have computed spectral functions using the Pad´eapproximation for the analytic continuation. The impu-rity t -matrix shows an excellent correspondence betweenthe N components at low frequencies, which demon-strates the good quality of the imaginary-time data ob-tained using the new scheme. The accuracy of the mag-netic spectra has been examined by comparison with theNCA for large N . Both results are in good agreement attemperatures of the order of T K . At lower temperatures,where the NCA cannot accurately reproduce the low-energy excitations, the present scheme still works prop-erly.The spectra have been tested using available Fermiliquid relations. We have found deviations from theKorringa-Shiba relation, whereas the t -matrix satisfiesthe Friedel sum rule for all N . The deviations are dueto the finite cutoff of the conduction band, which affectsthe static susceptibility but not the low-energy spectrum.These results reveal a weak point of the CT-QMC: it isdifficult to approach the universal regime with regard tothe static susceptibility.The present algorithm can easily be applied to DMFTsimulations of the Coqblin-Schrieffer lattice model.Namely, the effective medium is optimized with use ofthe impurity t -matrix. It is possible to address the for-mation of heavy quasi-particles through the conductionelectron Green function as well as the t -matrix of thelattice systems. Acknowledgment
One of the authors (J. O.) is supported by ResearchFellowships of the Japan Society for the Promotion ofScience for Young Scientists.
Appendix A: Derivation of update probabilitiesin the CS model from the Ander-son model
The CS model is derived from the Anderson modelby taking the localized limit with strong correlations. Namely, charge fluctuations are suppressed by infiniteCoulomb repulsion and the deep local level. We take thelimits ǫ f → −∞ and V → ∞ with J = − V /ǫ f fixed,where ǫ f and V denote the energy of the local level andstrength of hybridization, respectively. Similarly, updateprobabilities in the CS model can be derived from thecorresponding expressions in the Anderson model. In thisappendix, we demonstrate the way to take the limit inthe Monte Carlo formalism.We consider the Anderson model of spin-less fermionsfor simplicity. In this case, the localized limit leads to apotential scattering, or equivalently, the CS model with N = 1. The update probability for cutting a segment(addition of an anti-segment) is given by p ( k → k + 1) p ( k + 1 → k ) = V e lǫ f (cid:18) − det D (+) det D (cid:19) βl max k + 1 , (A ·
1) where l denotes a length of a segment which will be re-moved. D (+) is the matrix obtained by adding the opera-tors c † ( τ ) c ( τ + l ) to the end of the matrix D . In the limitof ǫ f → −∞ , the probability of such an update with l finite becomes 0 due to the factor e lǫ f . Hence l shouldbe reduced as the inverse of | ǫ f | in the update. For thispurpose, we introduce a cutoff length l defined by l = λ − ǫ f , ( λ ≫ . (A · l ≥ l , the update is negligible due to the factor e − λ .Hence we can restrict the length to l < l . The restrictionfor l replaces the factor l max with l in eq. (A · l = xl with 0 < x <
1, the update probability is givenin terms of J by p ( k → k + 1) p ( k + 1 → k ) = Jλ e − λx (cid:18) − det D (+) det D (cid:19) βk + 1 . (A · x in the limit of ǫ f → −∞ , we integrate out x as follows: Z d xλ e − λx ≃ . (A · λ → ∞ , which isrealized in the limit ǫ f → −∞ . As a result, we obtain theupdate probability that the localized electron is removedfor an infinitesimal time as follows: p ( k → k + 1) p ( k + 1 → k ) = J (cid:18) − det D (+) det D (cid:19) βk + 1 . (A · J αα X αα ( − c α c † α ) is added. It is obvious fromthis derivation that D (+) is the matrix with c † ( τ ) c ( τ +0) added to the original one, and that the equal-timeGreen function should be g (+0). In a similar manner,all formulae of transition probabilities in the CS modelcan be derived from the corresponding processes in theAnderson model. Appendix B: The Korringa-Shiba relation in amodel with finite band width
The Korringa-Shiba relation connects Im χ ( ω +i δ ) /ω | ω =0 to χ (0) by a universal value. The equationhas been proven in the wide band limit. Hence it may notbe satisfied in numerical calculations for systems with fi-nite band width, as in the present study. In order to clar-ify the deviation from the universal value, we considerthe non-interacting Anderson model with N -fold degen-eracy. Assuming a constant density of states, eq. (34),the Matsubara Green function is given by G − f (i ǫ n ) = i ǫ n − ǫ f + (2i∆ /π ) arctan( D/ǫ n ) ≃ i ǫ n a − − ǫ f + i∆sgn( ǫ n ) , (B · πV ρ and a − = 1 − πD . (B · /D ≪
1. The wide band limit corre-sponds to a = 1, and a finite band produces a correctionproportional to ∆ /D . Full Paper J. Otsuki et al.
We evaluate the dynamical susceptibility using theabove Green function including the effect of the finiteband width. The dynamical susceptibility is given by χ (i ν n ) = − N C N T X n ′ G f (i ǫ n ′ ) G f (i ǫ n ′ + i ν n ) , (B · ν n = 2 nπT is the boson Matsubara frequency.Evaluating the imaginary part at T = 0, we obtain thewell-known relation Im χ ( ω + i δ ) /ω | ω =0 = πN C N ρ f (0).It turns out that the quantity a does not affect the low-energy spectrum. On the other hand, the real part isinfluenced by the cutoff of the band. The sum over theMatsubara frequency is replaced by an integral along theimaginary-frequency axis, and evaluates to χ (0) = aN C N ρ f (0) . (B · a = 1, the Korringa-Shiba relation is sat-isfied. On the other hand, if a >
1, the static suscepti-bility is enhanced, and the Korringa-Shiba relation doesnot hold. In the Anderson model with U = 0 or in theCoqblin-Schrieffer model, ∆ is replaced by the width ofthe Kondo resonance, and therefore is of the order of theKondo temperature.
1) A. Georges, G. Kotliar, W. Krauth and M. J. Rozenberg: Rev.Mod. Phys. (1996) 13.2) A.N. Rubtsov, V.V. Savkin and A.I. Lichtenstein: Phys. Rev.B (2005) 035122. 3) P. Werner, A. Comanac, L.de’ Medici, M. Troyer and A.J.Millis: Phys. Rev. Lett. (2006) 076405; P. Werner and A.J.Millis: Phys. Rev. B (2006) 155107.4) The absence of a sign problem can be demonstrated in analogyto the proof for the Hirsch-Fye method in Yoo et al., J. Phys.A: Math. Gen. (2005) 10307. We thank R. Kaul for bringingthis to our attention.5) E. Gull, P. Werner, A.J. Millis and M. Troyer: cond-mat/0609438.6) P. Werner and A.J. Millis: Phys. Rev. B (2007) 085108.7) K. Haule: Phys. Rev. B (2007) 155113.8) P. Werner and A.J. Millis: cond-mat/0701730.9) B. Coqblin and J.R. Schrieffer: Phys. Rev. (1969) 847.10) A.L. Fetter and J.D. Walecka: Quantum Theory of Many-Particle Systems (McGraw-Hill, New York, 1971).11) Y. Kuramoto: Z. Phys. B (1983) 37; H. Kojima, Y. Ku-ramoto and M. Tachiki: Z. Phys. B (1984) 293.12) N.E. Bickers: Rev. Mod. Phys. (1987) 845.13) Y. Kuramoto and H. Kojima: Z. Phys. B (1984) 95.14) J. Otsuki, H. Kusunose and Y. Kuramoto: J. Phys. Soc. Jpn. (2006) Suppl. 256.15) V.T. Rajan: Phys. Rev. Lett. (1983) 308.16) H.J. Vidberg and J.W. Serene: J. Low Temp. Phys. (1977)179.17) M. Jarrell and J.E. Gubernatis: Phys. Rep. (1996) 133.18) K.G. Wilson: Rev. Mod. Phys. (1975) 773.19) O. Sakai, Y. Shimizu and T. Kasuya: J. Phys. Soc. Jpn. (1989) 3666.20) A.C. Hewson: The Kondo Problem to Heavy Fermions (Cam-bridge University Press, 1993).21) H. Shiba: Prog. Theor. Phys.54