Microscopic models of classical degrees of freedom coupled to non-interacting
fermions occur in many different contexts. Prominent examples from solid state
physics are descriptions of colossal magnetoresistance manganites and diluted
magnetic semiconductors, or auxiliary field methods for correlated electron
systems. Monte Carlo simulations are vital for an understanding of such
systems, but notorious for requiring the solution of the fermion problem with
each change in the classical field configuration. We present an efficient,
truncation-free O(N) method on the basis of Chebyshev expanded local Green
functions, which allows us to simulate systems of unprecedented size N.Comment: 4 pages, 3 figure

Alexander Weiße

G. H. Golub

J. M. Ziman

English

arXiv

Weiße, Alexander

Max Planck Society eDoc Server

Green-function-based Monte Carlo method for classical fields coupled to fermions

Weiße, A.

MPG.PuRe

arXiv:0904.3100v1  [cond-mat.stat-mech]  20 Apr 2009Green-Function-Based Monte Carlo Method for Classical Fields Coupled to FermionsAlexander WeißeInstitut fu¨r Physik, Ernst-Moritz-Arndt-Universita¨t Greifswald, 17487 Greifswald, Germany∗(Dated: January 7, 2009)Microscopic models of classical degrees of freedom coupled to non-interacting fermions occur in many dif-ferent contexts. Prominent examples from solid state physics are descriptions of colossal magnetoresistancemanganites and diluted magnetic semiconductors, or auxiliary field methods for correlated electron systems.Monte Carlo simulations are vital for an understanding of such systems, but notorious for requiring the solutionof the fermion problem with each change in the classical field configuration. We present an efficient, truncation-free O(N) method on the basis of Chebyshev expanded local Green functions, which allows us to simulatesystems of unprecedented size N .PACS numbers: 02.70.Ss, 71.15.-m, 75.47.LxThe numerical simulation of quantum lattice models is akey tool in solid state research and many other fields ofphysics. One class of problems, which is notoriously dif-ficult to study, are fermions coupled to classical degrees offreedom. Such microscopic models can arise if parts of acomplex system are approximated classically. A prominentexample is the double-exchange model, which describes theferromagnetism of mixed-valence manganites on the basis ofclassical t2g-spins whose orientation affects the kinetic energyof eg valence electrons. [1, 2, 3] Another example are Mn-doped (III,V) semiconductors, where itinerant holes triggera ferromagnetic ordering of the Mn spins. [4] A completelydifferent route that leads to a coupling of fermions and classi-cal degrees of freedom are the auxiliary field methods, whichtackle the problem of interacting fermions. Here, a Hubbard-Stratonovich transformation is used to decouple the two-bodyinteraction into non-interacting fermions in an auxiliary field,which is summed over with Monte Carlo methods. [5, 6]For all these systems the cause of the numerical difficulty isthe requirement for a solution of the non-interacting fermionproblem whenever the classical field is varied in a MonteCarlo simulation. We propose an efficient local-update algo-rithm, which obtains the change in fermionic energy directlyfrom a few local Green functions. These Green functions caneasily be calculated by Chebyshev expansion, superseding es-timates of the density of states and thus trace calculations. Weillustrate the efficiency of the approach with simulations of thedouble-exchange model.The models we consider in this work are of the general formH =∑ijc†iAij(~φ)cj , (1)where c(†)i are fermion creation (annihilation) operators atlattice site i, and ~φ is a classical field with one or morecomponents at each site. For example, the double-exchangemodel [7, 8] is given by the HamiltonianH = −∑〈ij〉tijc†icj , (2)where the summation is over nearest-neighbor sites and thehopping tij depends on the orientation {θi, φi} of classicallocal spins at each site,tij = cosθi−θj2 cosφi−φj2 + i cosθi+θj2 sinφi−φj2 . (3)This complex matrix element is one for ferromagneticallyaligned spins and vanishes for anti-ferromagnetic alignment.At low temperature the system favors ferromagnetism, sinceit can gain kinetic energy.The thermodynamics is described by the partition functionZ = Trc Trf exp[−β(H(~φ)− µN)] (4)and its derivatives. Here the traces Trc and Trf sum over theclassical and fermionic degrees of freedom, respectively. Thefermionic trace can be rewritten in terms of the single-particleeigenvalues ǫi of H ,Z = Trc exp[−Seff(~φ)] , (5)such that the grand potential of the fermions times β,Seff(~φ) = −∑ilog(1 + exp[−β(ǫi(~φ)− µ)]), (6)defines an effective Euclidean action for the classical degreesof freedom. The second trace over the classical field can thenbe calculated with a standard Monte Carlo sampling, wherethe weight of a configuration ~φ is given byP (~φ) = exp[−Seff(~φ)]/Z . (7)However, a non-trivial problem remains: To calculate P (~φ)we need to know the spectrum {ǫi} of the non-interactingfermion system, or at least its change under a proposed MonteCarlo update ~φ→ ~φ′.There are different solutions to this problem: We could usebrute force and calculate all eigenvalues ofAij(~φ) whenever ~φis modified. For local update schemes this is, of course, veryexpensive and imposes severe restrictions on the accessiblesystem sizes. As a way out, hybrid approaches [9, 10, 11, 12]have been suggested, where updates of the whole field config-uration are calculated with an approximate dynamics. Then,the solution of the full fermion problem in the acceptance step2is required less frequently. However, if the approximate dy-namics does not closely match the exact one, the acceptancerate drops markedly, in particular for increased system size.Therefore, these approaches crucially depend on the qualityof the approximate action.Staying with local updates of the classical field one can tryto optimize the calculation of the fermion spectrum. Motomeand Furukawa [13, 14, 15] suggested a Chebyshev expansionof the fermion density of states, which can be calculated withan effort proportional to the square of the system size N . Afurther modification [16, 17], involving several truncations inthe moment calculation, reduced the effort to order N .In another recent approach [18] the evolution of the eigen-values of Aij(~φ) under small local changes of ~φ is tracked us-ing special techniques for low-rank matrix updates [19]. Thenagain, the full solution of the fermion problem is required onlyoccasionally. In the best case this leads to N logN scaling.In the present work we combine ideas from the last twoapproaches and directly calculate the change of the fermiondensity of states using a few real space Green functions. Re-lying on Chebyshev expansion these can be calculated with aneffort proportional to the system size N . Without any trunca-tions we arrive at an order-N algorithm.Let us start from the Hamiltonian H with the Hermitianhopping matrix A and ask how the spectrum changes, when alocal modification ∆ is added. Given a(E) = A−E1 and itsinverse G(E), i.e. the Green function with G(E)a(E) = 1,the spectrum of A+∆ follows from(A+∆)|ψ〉 = E|ψ〉 ,(a(E) + ∆)|ψ〉 = 0 . (8)Right multiplication with G(E) yieldsG(E)(a(E) + ∆)|ψ〉 = (1+G(E)∆)|ψ〉 = 0 . (9)The spectrum of A + ∆ is given by those values of E wherethe determinantd(E) := det(1+G(E)∆) (10)vanishes. Recalling introductory lectures on Green functions(see e.g. Ref. 20) we note that this N -dimensional determi-nant reduces to one with a dimension equal to the rank of ∆.For local Monte Carlo updates this is a small number. If wechange, for instance, the on-site potential, ∆ has a single non-zero matrix element and we merely need the local Green func-tion of the corresponding site. For the double exchange modela single spin flip affects the hopping between the site and itsnearest neighbors. Independent of the system size N or thespace dimension, Eq. (10) then reduces to a 2 × 2 problem,i.e., we need only four Green functions connecting the sitewith its environment and both to themselves.Before we go into the details of calculating d(E), let usfurther analyze the meaning of this quantity. Going back toEq. (9) we haved(E) = det[G(E)(a(E) + ∆)]= det[G(E)] det[a(E) + ∆] .(11)FIG. 1: (color online). (a) Labeling of the nearest neighbor siteson a cubic lattice. (b) Comparison of Im log(d(E + i ε)) calculatedfrom exact eigenvalues ǫi, ǫ′i and with Chebyshev expansion of orderM = 512 (ε = 0.0625) for a 63 site sample.This can be expressed in terms of the eigenvalues {ǫi} of Aand {ǫ′i} of A+∆,d(E) =∏i1ǫi − E∏i(ǫ′i − E) . (12)An important trick of our new approach consists of going tothe complex plane, E → z := E + i ε, and taking a logarith-mic derivative. We now observe that d(E) determines exactlywhat we need for the Monte Carlo update: the change in thedensity of states going from A to A+∆, i.e., from ~φ to ~φ′,1πIm limε→0d log(d(z))dz=1πIm limε→0∑i1ǫi − z −1ǫ′i − z=∑iδ(ǫi − E)− δ(ǫ′i − E)= ρ(E)− ρ′(E) .(13)The change in the effective action then readsSeff(~φ′)− Seff(~φ) ==∫log(1 + e−β(E−µ))(ρ(E)− ρ′(E))dE=βπ∫11 + eβ(E−µ)Im limε→0log(d(E + i ε))dE ,(14)where in the last line partial integration led to an integral overthe Fermi function. The key role of d(z) has already beennoted earlier [21, 22], but the evaluation of Eq. (10) becomesfeasible only if we can restrict ourselves to a minimal set ofGreen functions and an efficient, direct method for their cal-culation.Let us now explain this main part of our new approach forthe double exchange model on a cubic lattice. Here, a localupdate consists of rotating a single spin at a site o. This modi-fies the matrix element tij between o and its nearest neighborsto the north, east, south, and so on. Labeling the sites accord-ing to Fig. 1(a), a naı¨ve evaluation of Eq. (10) requires all 7×73Green functions Gij(z) with i, j ∈ {o, n, e, s, w, t, b} [23].However, we can do much better observing thatd(z) = det(1+G(z)∆)= [1 +∑j∈n.n.∆joGoj(z)][1 +∑j∈n.n.∆ojGjo(z)]−Goo(z)[∑j,k∈n.n.∆jo∆okGkj(z)](15)can be expressed in terms of only 2× 2 Green functions,d(z) = [1 +Gov(z)][1 +Gvo(z)]−Goo(z)Gvv(z) , (16)which connect the original site o and the environment state|v〉 = ∆|o〉 =∑j∈n.n.∆jo|j〉 . (17)All four Green functions can be calculated easily withthe Chebyshev expansion approach outlined in a recent re-view [24]. In a nutshell, diagonal elements Gii(z) are ex-panded in terms of the Chebyshev polynomials of first andsecond kind, Tm and Um, respectively,Gii(E + i ε) =i[µ0 + 2∑M−1m=1 µmTm(E/s)]√s2 − E2+ 2M−1∑m=1µmUm−1(E/s)= iµ0 + 2∑M−1m=1 µm exp[− im arccos(E/s)]√s2 − E2 . (18)The expansion coefficients µm are Chebyshev moments mod-ified by appropriate kernel factors, which improve the conver-gence of the truncated series and damp Gibbs oscillations,µm = 〈i|Tm(H/s)|i〉sinh[λ(1−m/M)]sinhλ. (19)The scaling factor s ensures that the spectrum of the Hamilto-nian H/s falls within the domain of the Chebyshev polynomi-als [−1, 1]. For the double exchange model we can choose s tobe a little larger than the bare bandwidth in the ferromagneticcase, s > 6. The kernel parameter λ regulates the resolu-tion of the method versus the damping of Gibbs oscillations.A good value is λ = 4. The resolution at which the Greenfunction is approximated is given by ε = λs/M . The limitε → 0 thus corresponds to infinite expansion order M . In anumerical simulation, of course, M is always finite and weneed to extrapolate data for different M to obtain the limitingvalue of Seff(~φ′) − Seff(~φ). Since in Eq. (14) we integrate afunction of Gij(z) over the Fermi function, the maximal res-olution should be better than the thermal broadening of theFermi step, which is of the order of 1/β. Low temperatures,therefore, require higher expansion orders.The most time consuming step of the whole simulation isthe calculation of the moments 〈i|Tm(H/s)|i〉. Using the re-cursion relationTm(x) = 2xTm−1(x) − Tm−2(x) (20)0 0.05 0.1 0.15 0.2 0.25T00.20.40.60.81ML = 6L = 8L = 12L = 16L = 24L = 300.12 0.13 0.14T0.40.50.60.7U 4FIG. 2: (color online). Main panel: Magnetization versus temper-ature for the double-exchange model on 3D clusters with periodicboundary conditions (N = L3). Inset: Binder ratio U4 showing acrossing near T ≈ 0.134.it reduces to sparse matrix-vector multiplications, the cost ofwhich scales linearly with the system size N . With a furthertrick based onT2m+i = 2TmTm+i − Ti with i = 0, 1 (21)we obtain two moments per matrix-vector multiplication.Moreover, half of the moments vanish due to the relation|v〉 = ∆|o〉 and the special structure of the Chebyshev re-cursion (20).At this point we should also note the advantage over ap-proaches based on a full expansion of the density of states:For the latter the moments are given by traces, µm ∼Tr{Tm(H/s)}, instead of simple expectation values. Unlesstruncated or otherwise approximated [16], their calculation re-quires O(N2) operations.So far we have discussed only the diagonal Green functionsGii. In Ref. 24 we showed that symmetric Green functionsGij = Gji can be expanded in the same manner. However,in the double exchange model the spins induce local magneticfields which break this symmetry. We therefore derive the off-diagonal Green functions Gov and Gvo from the two diagonalfunctions Go+v,o+v and Go+i v,o+i v. With all required mo-ments at hand we can evaluate the sums in Eq. (18) with fastFourier methods, calculate d(E) with Eq. (16), and finally in-tegrate over the Fermi function to obtain the change of Seff.In Fig. 1(b) we show a typical example of Im log(d(E + i ε))and compare the expansion with the exact result from a fulldiagonalization.Having explained the technical details of the approach letus now illustrate its applicability with a few results for thedouble-exchange model (2) at half filling, µ = 0. In the mainpanel of Fig. 2 we show the magnetization as a function oftemperature, where the latter is measured in units of the max-imal hopping amplitude t ≡ 1. As expected, we observe aphase transition to a ferromagnetically ordered phase belowT ≈ 0.14. A closer inspection based on the Binder parame-4FIG. 3: (color online). Structure of the 3D spin field at two differenttemperatures visualized with stream lines.ter [25]U4 = 1− 〈m4〉3〈m2〉2 (22)yields the estimate Tc ≈ 0.134, see the inset of Fig. 2. Thisagrees with previous estimates [11, 15, 26] of the critical tem-perature, which range between 0.128 and 0.139.The data in Fig. 2 is based on expansions of orderM = 256and averages over 6000 to 20000 Monte Carlo steps in the crit-ical region, where one step corresponds to N spin flips. It isthe low resource consumption which allows for these far moreprecise calculations compared to previous studies [15, 26],which were based on Chebyshev expansions of order M ≈20. Note also that we can handle much larger systems withN = 303 sites, i.e. a complex matrix dimension of 27000.To visualize the spatial structure of the spin field we imag-ine it as a liquid flow and draw curves tangent to the velocityfield. These stream lines are regular and parallel in the or-dered phase, but quite irregular and swirled near and abovethe phase transition, see Fig. 3. Of course, this visualizationmethod fails for truly disordered spin fields.In summary, we have presented an efficient local updatescheme for Monte Carlo simulations of classical fields cou-pled to fermions. At the core of the approach is an expres-sion which relates the change of the fermion spectrum to afew local Green functions. These can be calculated easilywith Chebyshev expansion. Compared to similar expansionapproaches a full trace over the fermion system is spared,which directly leads to a fast and precise order-N algorithm.Possibly our method can be further accelerated using the ap-proximations inherent in the truncated polynomial expansionmethod [16, 17]. The calculations presented were performedon the TeraFLOPS cluster of the Institute for Physics at Greif-swald university.∗ New affiliation: Max-Planck-Institut fu¨r Mathematik, P.O.Box7280, 53072 Bonn, Germany[1] C. Zener, Phys. Rev. 82, 403 (1951).[2] P. W. Anderson and H. Hasegawa, Phys. Rev. 100, 675 (1955).[3] P.-G. de Gennes, Phys. Rev. 118, 141 (1960).[4] J. Schliemann, J. Ko¨nig, and A. H. MacDonald, Phys. Rev. B64, 165201 (2001).[5] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys. Rev.D 24, 2278 (1981).[6] J. E. Hirsch and R. M. Fye, Phys. Rev. Lett. 56, 2521 (1986).[7] E. M. Kogan and M. I. Auslender, Phys. Status Solidi B 147,613 (1988).[8] A. Weiße, J. Loos, and H. Fehske, Phys. Rev. B 64, 054406(2001).[9] R. T. Scalettar, D. J. Scalapino, and R. L. Sugar, Phys. Rev. B34, 7911 (1986).[10] S. Duane, A. Kennedy, B. J. Pendleton, and D. Roweth, Phys.Lett. B 195, 216 (1987).[11] J. L. Alonso, L. A. Ferna´ndez, F. Guinea, V. Laliena, andV. Martı´n-Mayor, Nucl. Phys. B 596, 587 (2001).[12] A. Weiße, H. Fehske, and D. Ihle, Physica B 359–361, 702(2005).[13] Y. Motome and N. Furukawa, J. Phys. Soc. Jpn. 68, 3853(1999).[14] Y. Motome and N. Furukawa, J. Phys. Soc. Jpn. 69, 3785(2000).[15] Y. Motome and N. Furukawa, J. Phys. Soc. Jpn. 70, 3186(2001), erratum.[16] N. Furukawa and Y. Motome, J. Phys. Soc. Jpn. 73, 1482(2004).[17] G. Alvarez, C. Sen, N. Furukawa, Y. Motome, and E. Dagotto,Comp. Phys. Comm. 168, 32 (2005).[18] G. Alvarez, P. K. V. V. Nukala, and E. D’Azevedo, J. Stat.Mech. p. P08007 (2007).[19] G. H. Golub and C. F. van Loan, Matrix Computations (JohnsHopkins University Press, Baltimore, 1996), 3rd ed.[20] J. M. Ziman, Principles of the Theory of Solids (CambridgeUniversity Press, Cambridge, 1972), 2nd ed.[21] F. Fucito, E. Marinari, G. Parisi, and C. Rebbi, Nucl. Phys. B180, 369 (1981).[22] D. J. Scalapino and R. L. Sugar, Phys. Rev. Lett. 46, 519 (1981).[23] In Ref. 22 this is described as a reduction to a L× L problem.[24] A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, Rev. Mod.Phys. 78, 275 (2006).[25] K. Binder, Z. Phys. B 43, 119 (1981).[26] Y. Motome and N. Furukawa, J. Phys. Soc. Jpn. 72, 2126(2003).

Crossref

Green-Function-Based Monte Carlo Method for Classical Fields Coupled to Fermions

Green-Function-Based Monte Carlo Method for Classical Fields Coupled to
  Fermions

Green-Function-Based Monte Carlo Method for Classical Fields Coupled to Fermions

Abstract

Similar works

Full text

Available Versions

Max Planck Society eDoc Server

MPG.PuRe

Crossref