Application of polynomial-expansion Monte Carlo method to a spin-ice Kondo lattice model
aa r X i v : . [ c ond - m a t . s t r- e l ] J u l Application of polynomial-expansion Monte Carlomethod to a spin-ice Kondo lattice model
Hiroaki Ishizuka , Masafumi Udagawa , , and Yukitoshi Motome Department of Applied Physics, University of Tokyo, Japan Max-Planck-Institut f¨ur Physik komplexer Systeme, Dresden, GermanyE-mail: [email protected]
Abstract.
We present the results of Monte Carlo simulation for a Kondo lattice model inwhich itinerant electrons interact with Ising spins with spin-ice type easy-axis anisotropy ona pyrochlore lattice. We demonstrate the efficiency of the truncated polynomial expansionalgorithm, which enables a large scale simulation, in comparison with a conventional algorithmusing the exact diagonalization. Computing the sublattice magnetization, we show theconvergence of the data with increasing the number of polynomials and truncation distance.
1. Introduction
Interplay between localized spins and itinerant electrons has been one of the major topics in thefield of strongly correlated electrons. It triggers various interesting phenomena, for instance, avariety of magnetic orderings induced by electron-mediated spin interactions, such as Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction [1] and the double-exchange interaction [2]. Theimpact of the spin-charge interplay is not limited to the magnetic properties, but also bringsabout peculiar electronic states and transport phenomena, such as non Fermi liquid behaviorin the quantum critical region in rare-earth systems [3] and the colossal magneto-resistance inperovskite manganese oxides [4].Recent studies on metallic pyrochlore oxides have opened yet another aspect in the field,namely, the geometrical frustration. The effect of geometrical frustration on the interplaybetween itinerant electrons and localized spins has attracted much interest, and extensivenumber of studies on these compounds have been reported. For example, interesting featureswere reported in Molybdenum compounds, such as an unconventional anomalous Hall effect [5]and emergence of a peculiar diffusive metallic phase [6]. A characteristic resistivity minimumwas also observed in an Iridium compound [7]. For these phenomena, the importance of localspin correlations inherent to the strong frustration has been suggested, but comprehensiveunderstanding is not reached yet.One of the authors and his collaborator recently reported unbiased Monte Carlo (MC)calculations of a Kondo lattice model on a pyrochlore lattice [8, 9, 10]. In their studies,however, the accessible system size N was limited to small sizes because the MC simulationwas performed by a conventional algorithm using the exact diagonalization (ED) [11], in whichthe computational amount increases with O ( N ). Larger size calculations are highly desired tofurther discuss the peculiar magnetic and transport phenomena in the frustrated spin-chargecoupled systems. For this purpose, here we apply another faster algorithm, the polynomialxpansion method (PEM) [12, 13]. We test the efficiency of the algorithm for a variant ofKondo lattice models on a pyrochlore lattice.
2. Model and method
We here consider a Kondo lattice model on a pyrochlore lattice, whose Hamiltonian is given by H = − t X h i,j i ,σ ( c † i,σ c j,σ + H . c . ) − J X i S i · σ i . (1)Here, c i,σ ( c † i,σ ) denotes an annihilation (creation) operator of electrons at site i with spin σ (= ↑ , ↓ ); S i and σ i represent the localized spin and itinerant electron spin, respectively. Themodel is defined on a pyrochlore lattice, three-dimensional frustrated lattice consisting of acorner-sharing network of tetrahedra. The sum h i, j i is taken over the nearest-neighbor (n.n.)sites on the pyrochlore lattice. We set t = 1 as the energy unit.We assume that the localized spins S i are Ising spins with | S i | = 1, whose anisotropyaxes depend on the four-sublattice sites on each tetrahedron: the axes are set along the h i directions, namely, the directions connecting the centers of neighboring tetrahedra. The situationis the same as in the spin ice [14, 15]. Although there is no bare interaction between the localizedspins in the model (1), the spins communicate with each other through an effective interactionmediated by the kinetic motion of electrons (RKKY interaction). Similar to the spin ice case,when the n.n. interaction is dominantly ferromagnetic (FM), a local spin configuration with twospins pointing in and the other two pointing out (two-in two-out) is favored in each tetrahedron.On the other hand, when the interaction is dominantly antiferromagnetic (AFM), all-in or all-outconfigurations become energetically stable. We note that the sign of the coupling J is irrelevantsince the Hamiltonian is unchanged for J → − J and S i → − S i .We apply MC calculations with PEM to the model (1). In the PEM, the density of states(DOS) for electrons under a spin configuration is obtained using the Chebyshev polynomialexpansion [12], which then is used to calculate the action in the weight for the given spinconfiguration. This algorithm reduces the computational cost compared to a conventional MCmethod based on ED [11]. We also implement the truncation method which further reducesthe computational cost [13]. We carry out the truncation by a real space distance, not by amagnitude of the matrix element in the original scheme; namely, we introduce a truncationdistance defined by the Manhattan distance from a flipped spin in calculating the Chebyshevmoments. The total computational amount for each MC step is largely reduced from O ( N ) forthe conventional algorithm to O ( N ) [13]. This enables us to access to larger system sizes.The calculations are done with varying the number of polynomials 15 ≤ m ≤
50 andtruncation distance 3 ≤ d ≤
7. In the following, we take J = 2 and restrict ourselves totwo typical electron densities; a low electron density n e = P i,σ h c † i,σ c i,σ i /N ∼ .
03 (the chemicalpotential is fixed at µ = − . n e ∼ .
35 ( µ = − . N = 4 × with periodic boundary conditions, typicallywith 2700 MC steps after 500 steps of thermalization. The results are divided into three binsof 900 steps to estimate the statistical error with 100 steps intervals in between each bin. MCsimulation costs about 15 hours with 8 cpu parallelization of Intel Xeon E5502 processors [16].
3. Results
Here, we show the MC results for the square of the sublattice magnetization, m s , withdifferent m and d . We calculate m s by the diagonal component of the spin structure factor S αβ ( q ) = P n,l h S αn S βl i exp(i q · r nl ) /N , as m s = 4 S αα ( q = ) /N . Here, n and l are the indicesof tetrahedra and α, β = 1 , , , d (b) T =0.07 T =0.08 T =0.09 (a) T =0.07 T =0.08 T =0.09 m s Figure 1.
MC resultsby the PEM for m s at anintermediate density n e ∼ .
35. We take d = 6 in (a) and m = 40 in (b). The resultsand errors of the ED methodare shown by horizontal solidlines and shades. (a) m m s d T =0.012 T =0.020 T =0.030 0.00.20.40.60.81.0 3 4 5 6 7 (b) Figure 2.
MC results bythe PEM for m s at a lowdensity n e ∼ .
03. We take d = 6 in (a) and m = 40 in(b). The results and errorsof the ED method are shownby horizontal solid lines andshades.Figure 1 shows the PEM results of m s at an intermediate density n e ∼ .
35 for threetypical temperatures( T ). Figure 1(a) shows m dependence at d = 6 and figure 1(b) is for d dependence at m = 40. At this electron density, the localized spins align in the all-in/all-outconfiguration (alternative arrangement of all-in and all-out tetrahedra) at low T since the n.n.RKKY interaction is dominantly AFM. The horizontal solid lines denote the results obtained bythe ED MC method. The ED result at T = 0 .
07 shows m s > .
8, which implies the system to bein the AFM ordered phase. At T = 0 . m s ∼ .
1, which indicates a disordered paramagneticstate. The intermediate T = 0 .
08 is presumed to be around the critical point. In figure 1(a),the PEM results at T = 0 .
07 and T = 0 .
09 converge to the ED results when m &
30. The dataat T = 0 .
08, however, show slower convergence: It appears to require m & ∼
40 for reliablecalculation. The convergence as to d shows a similar tendency, as shown in figure 1(b); d & T = 0 .
07 and T = 0 .
09, while a larger d & T = 0 . T range becomes lower; the RKKY energy scale gets smaller for lowerdensity. Another reason is that the Fermi energy comes close to the band edge (bottom); thePEM is an expansion technique of DOS, and hence, larger m is necessary for good convergencewhen DOS changes rapidly as in the band edge. Furthermore, when the electron density is small,a fluctuation of the density in MC measurements will harm the precision; in our calculations, n e has typically an error of ∆ n e ≃ .
01, which might affect the results when n e ∼ ∆ n e .Such situation is illustrated in figure 2 for an extremely low density n e ∼ .
03. In this region,the lowest T state is characterized by a FM order of two-in two-out tetrahedra with aligningthe net moment of each tetrahedra along a h i direction. Note that the T range in figure 2is much smaller than in figure 1. The data show much poor convergence to the ED results:Nevertheless, m &
40 and d & m ∼ −
40 and d ∼ − m ≤ J . In the previous study, J was takento be infinity, which greatly simplifies the band structure, while we take a much smaller value J = 2 in the present study. The resultant complicated band structure requires a larger numberof m to reproduce its fine details. Another possibility is the effect of geometrical frustrationin the present model. In general, the frustration leads to degeneracy among different spinconfigurations in a small energy window. For such situation, a small systematic error in thecalculation of the MC weight could lead to considerable errors in observables. Furthermore,the frustration suppresses all the energy scales, which might also require much larger efforts toobtain converged data.
4. Summary
To summarize, we examined the efficiency of the Monte Carlo simulation using the polynomialexpansion method for a frustrated Kondo lattice model with spin-ice type Ising spins. Bycomparison with the results by the exact diagonalization method, our result indicated that thepolynomial-expansion Monte Carlo calculation with m ∼ −
40 and d ∼ − m and d are necessaryfor convergence near the critical temperature and for very low electron density. Our resultsdemonstrate that the polynomial expansion method is a practical method for investigating thephysics of spin-charge coupled systems even when J is much smaller than the bandwidth. MonteCarlo study of the phase diagram by systematic analysis up to larger system sizes is in progress.The authors thank T. Misawa for fruitful discussions. H.I. is supported by Grant-in-Aid forJSPS Fellows. This work is supported by KAKENHI (Grants No. 19052008, 21340090, 21740242,and 22540372), the Global COE Program “the Physical Sciences Frontier”, and HPCI StrategicProgram, from MEXT, Japan. References [1] Ruderman M A and Kittel C 1954
Phys. Rev.
99; Kasuya T 1956
Prog. Theor. Phys.
45; Yosida K1957
Phys. Rev.
Phys. Rev. Phys. Rev.
Rev. Mod. Phys. Nanoscale Phase Separation and Colossal Magnetoresistance (Berlin:Springer-Verlag) and references there in[5] Taguchi Y, Ohhara Y, Yoshizawa H, Nagaosa N, and Tokura Y 2001
Science
Phys.Rev. Lett.
Phys. Rev. Lett.
96 087204[8] Motome Y and Furukawa N 2010
Phys. Rev. Lett.
J. Phys.: Conf. Ser.
Phys. Rev. B R060407[11] Yunoki S, Hu J, Malvezzi A L, Moreo A, Furukawa N, and Dagotto E 1998
Phys. Rev. Lett. J. Phys. Soc. Jpn. J. Phys. Soc. Jpn. Phys. Rev. Lett. Nature N .[17] Motome Y and Furukawa N 2003 Phys. Rev. B68