Application of polynomial-expansion Monte Carlo method to a spin-ice Kondo lattice model
We present the results of Monte Carlo simulation for a Kondo lattice model in which itinerant electrons interact with Ising spins with spin-ice type easy-axis anisotropy on a pyrochlore lattice. We demonstrate the efficiency of the truncated polynomial expansion algorithm, which enables a large scale simulation, in comparison with a conventional algorithm using the exact diagonalization. Computing the sublattice magnetization, we show the convergence of the data with increasing the number of polynomials and truncation distance.
Department of Applied Physics, University of Tokyo, Japan
Max-Planck-Institut für Physik komplexer Systeme, Dresden, Germany
Interplay between localized spins and itinerant electrons has been one of the major topics in the field of strongly correlated electrons. It triggers various interesting phenomena, for instance, a variety of magnetic orderings induced by electron-mediated spin interactions, such as Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction [RKKYInteraction] and the double-exchange interaction [DEMechanism]. The impact of the spin-charge interplay is not limited to the magnetic properties, but also brings about peculiar electronic states and transport phenomena, such as non Fermi liquid behavior in the quantum critical region in rare-earth systems [Stewart2001] and the colossal magneto-resistance in perovskite manganese oxides [Dagotto2002].
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 interplay between itinerant electrons and localized spins has attracted much interest, and extensive number of studies on these compounds have been reported. For example, interesting features were reported in Molybdenum compounds, such as an unconventional anomalous Hall effect [Taguchi2001] and emergence of a peculiar diffusive metallic phase [Iguchi2009]. A characteristic resistivity minimum was also observed in an Iridium compound [Nakatsuji2006]. For these phenomena, the importance of local spin correlations inherent to the strong frustration has been suggested, but comprehensive understanding 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 [Motome2010-1, Motome2010-2, Motome2010-3]. In their studies, however, the accessible system size was limited to small sizes because the MC simulation was performed by a conventional algorithm using the exact diagonalization (ED) [Yunoki1998], in which the computational amount increases with . Larger size calculations are highly desired to further discuss the peculiar magnetic and transport phenomena in the frustrated spin-charge coupled systems. For this purpose, here we apply another faster algorithm, the polynomial expansion method (PEM) [Motome1999, Furukawa2004]. We test the efficiency of the algorithm for a variant of Kondo 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
Here, () denotes an annihilation (creation) operator of electrons at site with spin ; and represent the localized spin and itinerant electron spin, respectively. The model is defined on a pyrochlore lattice, three-dimensional frustrated lattice consisting of a corner-sharing network of tetrahedra. The sum is taken over the nearest-neighbor (n.n.) sites on the pyrochlore lattice. We set as the energy unit.
We assume that the localized spins are Ising spins with , whose anisotropy axes depend on the four-sublattice sites on each tetrahedron: the axes are set along the directions, namely, the directions connecting the centers of neighboring tetrahedra. The situation is the same as in the spin ice [Harris1997, Ramirez1999]. Although there is no bare interaction between the localized spins in the model (2), the spins communicate with each other through an effective interaction mediated 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 two spins 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-out configurations become energetically stable. We note that the sign of the coupling is irrelevant since the Hamiltonian is unchanged for and .
We apply MC calculations with PEM to the model (2). In the PEM, the density of states (DOS) for electrons under a spin configuration is obtained using the Chebyshev polynomial expansion [Motome1999], which then is used to calculate the action in the weight for the given spin configuration. This algorithm reduces the computational cost compared to a conventional MC method based on ED [Yunoki1998]. We also implement the truncation method which further reduces the computational cost [Furukawa2004]. We carry out the truncation by a real space distance, not by a magnitude of the matrix element in the original scheme; namely, we introduce a truncation distance defined by the Manhattan distance from a flipped spin in calculating the Chebyshev moments. The total computational amount for each MC step is largely reduced from for the conventional algorithm to [Furukawa2004]. This enables us to access to larger system sizes.
The calculations are done with varying the number of polynomials and truncation distance . In the following, we take and restrict ourselves to two typical electron densities; a low electron density (the chemical potential is fixed at ), for which the n.n. RKKY interaction is FM, and an intermediate electron density (), for which the n.n. RKKY interaction is AFM. All the following calculations were done for with periodic boundary conditions, typically with 2700 MC steps after 500 steps of thermalization. The results are divided into three bins of 900 steps to estimate the statistical error with 100 steps intervals in between each bin. MC simulation costs about 15 hours with 8 cpu parallelization of Intel Xeon E5502 processors [note].
Here, we show the MC results for the square of the sublattice magnetization, , with different and . We calculate by the diagonal component of the spin structure factor , as . Here, and are the indices of tetrahedra and are the indices of sites in a tetrahedron.
Figure LABEL:fig:all shows the PEM results of at an intermediate density for three typical temperatures(). Figure LABEL:fig:all(a) shows dependence at and figure LABEL:fig:all(b) is for dependence at . At this electron density, the localized spins align in the all-in/all-out configuration (alternative arrangement of all-in and all-out tetrahedra) at low since the n.n. RKKY interaction is dominantly AFM. The horizontal solid lines denote the results obtained by the ED MC method. The ED result at shows , which implies the system to be in the AFM ordered phase. At , , which indicates a disordered paramagnetic state. The intermediate is presumed to be around the critical point. In figure LABEL:fig:all(a), the PEM results at and converge to the ED results when . The data at , however, show slower convergence: It appears to require for reliable calculation. The convergence as to shows a similar tendency, as shown in figure LABEL:fig:all(b); is sufficient at and , while a larger appears to be necessary at .