Effective Ruderman-Kittel-Kasuya-Yosida-like interaction in diluted double-exchange model: self-learning Monte Carlo approach
EEffective Ruderman-Kittel-Kasuya-Yosida-like interaction in diluted double-exchangemodel: self-learning Monte Carlo approach
Hidehiko Kohshiro
1, 2 and Yuki Nagai
2, 3 Institute for Solid State Physics, The University of Tokyo, Kashiwa, Chiba, Japan CCSE, Japan Atomic Energy Agency, 178-4-4, Wakashiba, Kashiwa, Chiba, 277-0871, Japan Mathematical Science Team, RIKEN Center for Advanced Intelligence Project (AIP),1-4-1 Nihonbashi, Chuo-ku, Tokyo 103-0027, Japan (Dated: May 15, 2020)We study the site-diluted double exchange (DE) model and its effective Ruderman-Kittel-Kasuya-Yosida-like interactions, where localized spins are randomly distributed, with the use of the Self-learning Monte Carlo (SLMC) method. The SLMC method is an accelerating technique for Markovchain Monte Carlo simulation using trainable effective models. We apply the SLMC method to thesite-diluted DE model to explore the utility of the SLMC method for random systems. We check theacceptance ratios and investigate the properties of the effective models in the strong coupling regime.The effective two-body spin-spin interaction in the site-diluted DE model can describe the originalDE model with a high acceptance ratio, which depends on temperatures and spin concentration.These results support a possibility that the SLMC method could obtain independent configurationsin systems with a critical slowing down near a critical temperature or in random systems where afreezing problem occurs in lower temperatures.
I. INTRODUCTION
Interactions between itinerant electrons and localizedspins have attracted much attention in magnetic mate-rials such as Manganite . The double-exchange (DE)model, a fundamental model of the itinerant magnetism,has also attracted much attention from theoretical andexperimental points of view, since its itinerant natureyields various rich phases .The electron-spin interaction in weak coupling regimein the DE model is approximated by the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction between lo-calized spins, whose coupling constant oscillates and de-cays with the distance between localized spins . Inthe diluted magnetic alloys ( e.g. , AuFe and CuMn) , thespin glass (SG) order, where spins freeze in the spatiallyrandom manner, is explained by existence of the RKKYinteraction between randomly localized spins. The SGstate shows novel transport phenomena due to the ran-dom spin structure . The nature of SG, however, isstill under debate even if it is treated as a simplifiedclassical spin model, so-called Edwards-Anderson model,whose exchange coupling constants are given by a certainprobability distribution .Even treating localized spins as classical ones, the stud-ies beyond the weak coupling regime are not trivial sinceitinerant electrons play a more important role in the DEmodel due to its quantum nature. For example, the ex-istence of the SG order in the site-diluted DE modelhas not been theoretically confirmed yet. A standardmethod to simulate classical spins is the Markov chainMonte Carlo (MCMC) method. In the DE model, how-ever, computational complexity is huge due to a fermiondeterminant for itinerant electrons. For example, thisnaive calculation takes O ( N ) per sweep with the systemsize N since a full diagonalization is needed to estimate the fermion determinant for given spin configurations foreach Monte Carlo step. Many alternative algorithms tosimulate electrons coupled to classical fields have beendeveloped .In MCMC methods for physics, obtaining independentconfigurations sometimes becomes hard near a criticaltemperature or at low temperatures. A critical slowingdown of the autocorrelation time near the critical tem-perature is a famous effect when the update algorithmfor generating configurations is local. Cluster global up-date algorithms such as the Wolff and the Swendsen-Wang algorithms are usually adopted for reducing theautocorrelation time if one can consider a model wherethese algorithms are available. The freezing problem iscaused by the multi-valley structure of the energy land-scape manifested at low temperatures, and the states aretrapped in the local minima. The freezing problem ap-pears in random systems such as the SG systems andmakes difficult its MCMC simulations. In principle, onecan reach the true minima in finite-temperature systemsafter many Monte Carlo (MC) steps.Recently, an efficient MCMC algorithm combined witha machine-learning technique, the self-learning MonteCarlo (SLMC) method, was proposed . The SLMCmethod constructs effective models to generate the Boltz-mann weight by means of machine learning from gatheredoriginal Markov chains. Once an effective model pre-pared, we can perform MCMC simulations of the originalmodels using the effective models trained. It was shownthat SLMC reduces the number of the MC steps of theoriginal model even near the critical slowing down .The SLMC method accelerates a simulation of a modelwhere cluster algorithms are unavailable by using an ef-fective model where cluster algorithms are available suchas the Ising model . As shown in Fig. 1, with the use ofthe SLMC, we can bypass heavy numerical computationscoming from an estimation of the original model. This a r X i v : . [ c ond - m a t . d i s - nn ] M a y A B training by MCMC of the original model ℋ eff ( S ) 𝑊 ( 𝐴 ) 𝑊 ( 𝐵 ) 𝑊 eff ( 𝐴 ) 𝑊 eff ( 𝐵 )min { 𝑊 ( 𝐵 ) 𝑊 eff ( 𝐴 ) 𝑊 ( 𝐴 ) 𝑊 eff ( 𝐵 ) } proposal of B by MCMC of the effective modelproposal of B by MCMC of the original modelMetropolis-Hastings cheaper SLMC expensive! = exp ( − 𝛽ℋ eff ( 𝐴 ) ) = exp ( − 𝛽ℋ eff ( 𝐵 ) ) FIG. 1. A brief schematic summary of the SLMC procedure.The computational costs of the SLMC method are cheaperthan one of the naive Monte Carlo calculation. property is also useful for solving the freezing problem.It is not trivial whether the SLMC method is efficientfor random systems. If one uses a very complicated effec-tive model with deep neural networks, the original modelmight be reproduced by the effective model. However, interms of the computational complexity of actual MCMCsimulations, a simple effective model is needed to by-pass heavy numerical computations. In the week cou-pling regime, there is a simple effective model known asthe RKKY interaction, which was used also for the ran-domly distributed spin systems .In this paper, we show that there is a simple effectivemodel similar to the RKKY interaction even in strongcoupling regime. To prove it, We employed the SLMCmethod with the cumulative update for fermions coupledto the classical field . We first checked success or fail-ure through the acceptance ratio. Next, we investigatedthe properties of the effective models trained. The mainissues are the dependencies of temperature and distancebetween localized spins of the effective models trained.The temperature dependence was not investigated in thesite-diluted DE model and also in the regular DE modelpreviously.The letter is organized as follows. First, we introducethe site-diluted DE model. Next, we explain how to per-form the SLMC simulation and show obtained results forthe site-diluted DE model. Finally, we summarize andgive some discussions. II. MODEL AND METHOD
In this paper, we focus on the site-diluted DE modelon three-dimensional cubic lattice to examine the effi-ciency of the SLMC method. The site-diluted DE model( i.e. , site-diluted ferromagnetic Kondo model) is definedas follows:ˆ H = − t (cid:88) σ, (cid:104) i,j (cid:105) (cid:16) ˆ c † iσ ˆ c jσ + h . c . (cid:17) − J (cid:88) i p i S i · ˆ σ i − µ (cid:88) σ,i ˆ c † iσ ˆ c iσ , where S i and ˆ σ i are the classical Heisenberg spins andthe Pauli matrices respectively. t > J > p i is a binary randomnumber (0 or 1) which represents whether the site i isoccupied ( p i = 1) or not ( p i = 0) with a localized spin.When p i is always 1 for a arbitrary i , The model corre-sponds to the regular DE model. It is known the regularDE model holds the ferromagnetic order for enough large J/t . For simplicity, we utilize a semiclassical approx-imation where localized spins are classical ones and en-ergy of the electrons changes instantaneously with changeof the localized spins i.e., the localized spins are “ time-independent ”.The partition function of the model Z is given by Z = (cid:88) S W ( S ) , where W ( S ) ≡ (cid:81) n (cid:0) β ( µ − E n ( S )) (cid:1) is the Boltzmannweight with a spin configuration S . Here, { E n ( S ) } is theeigen spectrum.We show a summary of the procedure of the SLMCsimulation in Fig. 1. The SLMC method is based on twoprocedures: “learn” and “earn”. “Learn” is the trainingprocedure of the effective models. We generate W ( S ) ofthe original models as the learning data by performingMCMC simulation of the original model. We optimizethe effective model by using the generated learning dataso that it generates W eff ( S ) closer to original W ( S ). Weadopt a site-diluted two-body classical spin interactionup to 6th neighbor, as in the previous study for the reg-ular DE model . The effective Hamiltonian is expressedas H eff = E − (cid:88) (cid:104) i,j (cid:105) n J eff n p i p j S i · S j , where, { E , J eff n } are trainable parameters which we op-timize. The difference from the previous study is thatwe only consider the interactions between randomly ar-ranged spins. In the appendix , we explain details foroptimizing the effective model.“Earn” is the actual simulation of the original modelusing the effective model optimized. The acceptance ra-tio from current state S to a new state generated by theMCMC of the effective model S (cid:48) is given by p ( S → S (cid:48) ) = min (cid:26) , W ( S (cid:48) ) W eff ( S ) W ( S ) W eff ( S (cid:48) ) (cid:27) , which is known as the Metropolis-Hastings test. We notethat a proposed state S (cid:48) is cumulatively updated andthe length of the effective Markov chain is not limited.We can update states as long as we need to reduce auto-correlations since the acceptance ratio is converged as afunctional of the length of the effective Markov chain.In the effective MCMC, we use the heatbath and theoverrelaxation procedure, which is known as an effectiveupdate algorithm in classical spin systems. After 200 ef-fective MC steps, where one step includes one heatbathseeps and four overrelaxations, the Metropolis-Hastings T a cc e p t a n c e r a t e N = 64 N = 64 N = 64 N = 64 N spin = 64 N spin = 48 N spin = 32 N spin = 16 FIG. 2. The acceptance ratios versus the temperature withdifferent spin concentrations for localized spins. The accep-tance ratios keep at least 40 % until T = 0 . test determined by Eq. (II) is done with the proposed spinconfiguration. We use an annealing-like training processto obtain effective models with different temperatures.We first gather training data set and train the effectivemodel H eff at a higher temperature (at T = 0 . t as unity and J/t = 16 and µ = 8, which possesthe ferromagnetic order in the regular case and keep nearthe quarter filling. We used the 3d cubic lattice N = 4 and we averaged over 5 different realizations with thefixed spin concentration and investigated with decreasingspin concentrations. III. RESULTSA. Acceptance ratios
We investigate the averaged acceptance ratio of theSLMC. The acceptance ratio shows how many proposedspin configurations are accepted during simulations. Ifthe learning procedure is failed, the effective model can-not fit the original Boltzmann weight W ( S ) with a givenspin configuration S . Such as bad proposal is almost re-jected and the acceptance ratio becomes very low, whichmeans a failure of a simulation. The averaged acceptanceratio is estimated as p ∼ e −√ MSE with the use the mean-squared error (MSE) between the original and effectiveenergy given as MSE = 1 N S (cid:88) S | log W ( S ) − log W eff ( S ) | , where N S is the number of spin configurations. A highacceptance ratio is important for reducing an autocor-relation time . We define the autocorrelation time τ A for the observable A in the original MC simulation. Ifthe number of effective MC steps is larger than τ A , theproposed spin configuration is uncorrelated. Since theproposed configuration is accepted with the probability p , an uncorrelated configuration is obtained after 1 /p times Metropolis-Hastings tests. Therefore, the accep-tance ratio represents the quality of the effective modeland criteria whether the simulations work well or not.The averaged acceptance ratios are shown in Fig. 2in systems with different classical spin concentrations.In the regular model where classical spins are locatedon all lattice points ( N spin = 64), there is a ferromag-netic transition at T c ∼ . . We find that there isno drastic change of the acceptance ratio around a crit-ical temperature. This is naturally explained by thatthe effective interactions do not change between differ-ent phases. The temperature dependence of the ac-ceptance ratio is explained by the temperature depen-dence of the MSE as follows. The temperature depen-dence of the effective action log W eff ( S ) is described bylog W eff ( S , β ) = − β H eff ( β ). Therefore, the averaged ac-ceptance ratio is estimated as p ∼ exp (cid:16) − (cid:112) β (cid:112) MSE H (cid:17) , where MSE H ≡ N S (cid:80) S (cid:12)(cid:12)(cid:12) log W ( S ) β − H eff ( S , β ) (cid:12)(cid:12)(cid:12) . Sincethe temperature dependence of the acceptance ratioshown in Fig. 2 is well explained by the function e − a √ /T ,the mean-squared error of the original and effectiveHamiltonian does not depend on the temperature much.We note that the effective Hamiltonian itself depends onthe temperature as shown in the latter part of this paper.The above analysis can be applied to cases of differentspin concentrations. The acceptance rates are kept 40%below to T = 0 .
1. This means the SLMC method workswell with the site-diluted DE model. We find that if asystem has less localized spins the acceptance ratio be-comes high. This suggests that there is an effective two-body classical spin-spin interaction in the strong couplingregime where
J/t = 16 in this paper.
B. Spin concentration and temperaturedependence of the effective models
We investigate the spin concentration and tempera-ture dependence of the effective model with the use ofthe SLMC. In the weak coupling regime, it is well knownthat the RKKY interaction does not depend on temper-ature since the RKKY interaction is derived from thesecond-order perturbation method. It is also known thatthe RKKY interaction does not depend on the spin con-figuration since this is a two-body interaction betweenspins. In the regular DE model reported in the previous
TABLE I. Effective long-range interactions T = 0 . N spin E J J J J J J
64 -55.253(1) 0.0691(1) -0.00739(9) 0.00234(5) -0.00158(6) 0.00047(4) -0.00033(3)56 -48.02(1) 0.0732(2) -0.00747(9) 0.0022(2) -0.0020(2) 0.00047(9) -0.00034(5)48 -40.85(2) 0.0781(2) -0.0072(2) 0.0017(2) -0.0023(1) 0.00051(9) -0.00027(8)40 -33.77(3) 0.0837(7) -0.0069(2) 0.00120(6) -0.0015(2) 0.00024(5) -0.00003(10)32 -26.86(3) 0.091(1) -0.0068(2) 0.0011(1) -0.0017(2) 0.00028(6) -0.00006(10)24 -20.02(2) 0.095(1) -0.0064(4) 0.0006(2) -0.0021(2) 0.0002(1) 0.00002(12)16 -13.36(2) 0.101(1) -0.0070(5) -0.00006(42) -0.0017(2) -0.00006(5) 0.00009(11)TABLE II. Effective long-range interactions T = 0 . N spin E J J J J J J
64 -52.30(1) 0.0744(7) -0.0087(2) 0.0041(1) -0.0004(4) 0.0007(2) -0.0004(2)56 -45.33(1) 0.0794(7) -0.0099(4) 0.0043(5) -0.0015(2) 0.0007(2) -0.0004(2)48 -38.38(2) 0.0841(9) -0.0087(5) 0.0035(5) -0.0023(5) 0.0001(2) -0.0005(2)40 -31.58(6) 0.090(2) -0.0099(9) 0.004(1) -0.0017(8) 0.0006(1) -0.0006(5)32 -24.97(5) 0.101(3) -0.010(1) 0.0034(7) -0.0015(5) 0.0009(2) -0.0012(2)24 -18.46(2) 0.111(2) -0.0103(8) 0.0028(2) -0.0025(7) 0.0012(2) -0.0007(2)16 -12.25(3) 0.116(3) -0.013(2) 0.002(1) -0.0015(9) 0.0001(4) 0.0003(4) J e ff n N spin / N = 16/64 T = 0.1 T = 0.4 N spin / N = 32/64 N spin / N = 48/64 N spin / N = 64/64 distance FIG. 3. The coupling constants J eff n of the trained effectivemodels versus distances between localized spins at T = 0 . T = 0 . paper , it was found that the effective coupling con-stants oscillate and decay as a function of the distancebetween localized spins, which is similar to the RKKYinteractions.In Fig. 3, we show the coupling constants J eff n of thetrained effective models against distances between spinsat T = 0 . T = 0 . T = 0 .
4) and the Table II ( T =0 . J eff n . Here, n = 1 means the nearest neighbor coupling and n = 2 the next-nearest neighbor coupling, and so on.As shown in Fig. 3, we find the RKKY features, os-cillating and decaying behavior as a functional of thedistance between spins, for all spin concentrations atboth a higher temperature T = 0 . T = 0 .
1. However, the amplitude of the cou-plings depends on the spin concentrations and the tem-peratures. The spin concentration dependence suggeststhat this effective ”RKKY” interaction in the strong cou-pling regime can not be derived from second-order per-turbation theory. The high acceptance ratio is shownin Fig. 2 indicates that this effective two-body interac-tion well describes the original model and there are somerenormalizations from many-body effects.We show the temperature dependence of the effectivecoupling constants J eff n in the trained effective modelsfor each distance n with different spin concentrations, asshown in Fig. 4. The characteristics are different for thenearest neighbor n = 1 and other cases. For the nearestneighbor n = 1, temperature dependencies are relativelysmaller with small error bars. On the other hand, forother cases, temperature dependencies are nonmonotonicand relatively larger. The number of further interactionsis fewer than nearer ones, therefore, they are relativelyirrelevant and can easily fluctuate. Since the temperaturedependence of the effective coupling constants is not sostrong, the ”annealing” procedure that we use to obtainthe effective model in lower temperature works well inthe site-diluted DE model.As seen in the acceptance ratios calculated, the tem-perature dependence of the effective model is small. Onthe other hand, we observed temperature-dependent be- J e ff n n = 1 J e ff n n = 1 J e ff n n = 1 J e ff n n = 1 N spin = 64 N spin = 48 N spin = 32 N spin = 16 n = 2 n = 2 n = 2 n = 2 n = 3 n = 3 n = 3 n = 3 n = 4 n = 4 n = 4 n = 4 n = 5 n = 5 n = 5 n = 5 n = 6 n = 6 n = 6 n = 6 T N = 64 N = 64 N = 64 N = 64 FIG. 4. The coupling constants J eff n of the trained effective models versus temperature. The error bars are defined as standarderror over different localized spin realizations havior in the effective coupling constant trained. Theeffective model H eff is optimized minimizing MSE H andsatisfying H eff ( S , β ) (cid:39) log W ( S ,β ) β . This means an effec-tive model should mimicks not energy but local “free en-ergy” at finite temperature. The “entropy”, therefore,may cause temperature-dependent behavior of effectivemodels.There is another possible explanation for the existenceof temperature dependence. In terms of the quantumfield theory, one has to integrate out the electron degreesof freedom to obtain effective classical interactions. Twoclassical spins interact with each other through electronGreen’s functions in the effective Lagrangian. At zerotemperature, such a Green’s function has been obtainedas a result of the exact summation of the Born series .At finite temperature, the electron Green’s function de-pends on temperatures. The SLMC imitates this effectivemodel in the Hamiltonian formalism with the use of thetemperature-dependent coupling constant.We show that the SLMC works well in random sys-tems and there is a simple effective model. There aretwo ways to study the SG transition in the site-dilutedDE model. One is the MC simulation with effective two-body interactions obtained by the SLMC. We find that,even in the strong coupling regime, there is an effectivetwo-body interaction to describe the original DE model.Therefore, without doing the Metropolis-Hastings test,one might have good accuracy in the classical MC sim-ulation with this effective interaction. The other is the SLMC with the use of the two-body classical spin inter-actions. To grasp the tail of the SG transition, one hasto calculate systems in various sizes with fixed spin con-centration, which takes much time. By accelerating theMCMC simulations with the use of the SLMC, one canstudy whether the SG transition occurs in the DE modelor not, which is still an open question. We note thata transition temperature of the SG transition would bemuch lower than one of the ferromagnetic transition dueto spin frustrations. For the challenge to the problems ofthe SG transition in the itinerant electron systems, largersystem sizes and much lower-temperature simulations areneeded. The temperature exchange method is usuallyused for the classical localized SG models and will makeour simulation better, but combinations of the tempera-ture exchange method and the SLMC method are futurework. IV. SUMMARY
We performed the self-learning Monte Carlo methodsto the semiclassical site-diluted double-exchange model.With decreasing localized spin concentrations, we ob-served modest acceptance rates and the effective modelskept RKKY behavior, which oscillates and decays withdistance. As well as the regular DE model, the SLMCmethod works well with the diluted DE model. We foundthat effective RKKY-like interaction depends on the spinconcentrations and temperature, while the RKKY inter-action derived in the weak coupling regime does not de-pend on both. We showed that the SLMC works well inrandom systems and there is a simple effective model ina strong coupling regime.
ACKNOWLEDGMENT
The calculations were partially performed by the su-percomputing systems SGI ICE X at the Japan AtomicEnergy Agency. This work was partially supported byJSPS–KAKENHI Grant Numbers 18K03552 to Y.N.,and the “Topological Materials Science” (No. 18H04228)JSPS–KAKENHI on Innovative Areas to Y.N..
Appendix: details of optimization procedure ofeffective models
In this appendix, we will explain datails how to opti-mize the effective models. We assume the effective modelas H eff ( S , β ) = (cid:80) (cid:104) i,j (cid:105) n J eff n ( β ) S i · S j ( S ) and consider upto N J eff th neighbor interactions at inverse temperature β . We prepare M states S , . . . , S M as training data and C n ( S ) ≡ (cid:80) (cid:104) i,j (cid:105) n S i · S j ( S ). A matrix C and a vector J are defined as follows: C ≡ C ( S ) · · · C N J eff ( S )... ... . . . ...1 C ( S M ) · · · C N J eff ( S M ) , J ≡ E J eff1 ... J eff N J eff . We also define a vector H ≡ t ( H ( S ) , · · · , H ( S M )) isaligined energies calculated via states as traning data.For the acutual simulation, you should restore C usedthe following procedure. J which give a minimum of (cid:107) H − C J (cid:107) is a solutionof the normal equation t CC J = t C H and it becomes J = ( t CC ) − t C H . The number of learning parame-ters to train is a few (only N J eff parameters), a directcalculation of the N J eff + 1 dimensitonal inverse matrix( t CC ) − is not expensive. Of cource, if the number oflearning parameters is large, you sholud use an iterativemethod. C. Zener, Phys. Rev. , 403 (1951). P. W. Anderson and H. Hasegawa, Phys. Rev. , 675(1955). K. Kubo and N. Ohata, Journal of the Physical Society ofJapan , 21 (1972), https://doi.org/10.1143/JPSJ.33.21. T. A. Kaplan and S. D. Mahanti,
Physics of manganites (Springer US, 2002). M. B. Salamon and M. Jaime, Rev. Mod. Phys. , 583(2001). K. Nagai, T. Momoi, and K. Kubo, Journalof the Physical Society of Japan , 1837 (2000),https://doi.org/10.1143/JPSJ.69.1837. S. Kumar and P. Majumdar, Phys. Rev. Lett. , 246602(2003). D. Pekker, S. Mukhopadhyay, N. Trivedi, and P. M. Gold-bart, Phys. Rev. B , 075118 (2005). S. Kumar and J. van den Brink, Phys. Rev. Lett. ,216405 (2010). S. Minami and H. Kawamura, Journal of thePhysical Society of Japan , 044702 (2015),https://doi.org/10.7566/JPSJ.84.044702. Y. Akagi and Y. Motome, Journal of thePhysical Society of Japan , 083711 (2010),https://doi.org/10.1143/JPSJ.79.083711. M. A. Ruderman and C. Kittel, Phys. Rev. , 99 (1954). T. Kasuya, Progress of Theoretical Physics ,45 (1956), https://academic.oup.com/ptp/article-pdf/16/1/45/5266722/16-1-45.pdf. K. Yosida, Phys. Rev. , 893 (1957). J. A. Mydosh,
Spin glasses: an experimental introduction (Taylor & Francis, 1993). G. Tatara and H. Kawamura, Journal of thePhysical Society of Japan , 2613 (2002),https://doi.org/10.1143/JPSJ.71.2613. Y. Niimi, M. Kimata, Y. Omori, B. Gu, T. Ziman,S. Maekawa, A. Fert, and Y. Otani, Phys. Rev. Lett. ,196602 (2015). K. Binder and A. P. Young, Rev. Mod. Phys. , 801(1986). J. A. Mydosh, Reports on Progress in Physics , 052501(2015). H. Kawamura and T. Taniguchi (Elsevier, 2015) pp. 1 –137. J. Alonso, L. Fernndez, F. Guinea, V. Laliena, andV. Martn-Mayor, Nuclear Physics B , 587 (2001). N. Furukawa, Y. Motome, and H. Nakata, ComputerPhysics Communications , 410 (2001), conference onComputational Physics 2000: ”New Challenges for theNew Millenium. N. Furukawa and Y. Motome, Journal of thePhysical Society of Japan , 1482 (2004),https://doi.org/10.1143/JPSJ.73.1482. G. Alvarez, C. en, N. Furukawa, Y. Motome, andE. Dagotto, Computer Physics Communications , 32(2005). G. Alvarez, P. K. V. V. Nukala, and E. D’Azevedo, Journalof Statistical Mechanics: Theory and Experiment ,P08007 (2007). A. Weiße, Phys. Rev. Lett. , 150604 (2009). K. Barros and Y. Kato, Phys. Rev. B , 235101 (2013). U. Wolff, Phys. Rev. Lett. , 361 (1989). R. H. Swendsen and J.-S. Wang, Phys. Rev. Lett. , 86(1987). J. Liu, Y. Qi, Z. Y. Meng, and L. Fu, Phys. Rev. B ,041101 (2017). J. Liu, H. Shen, Y. Qi, Z. Y. Meng, and L. Fu, Phys. Rev.B , 241104 (2017). R. W. Walstedt and L. R. Walker, Phys. Rev. Lett. ,1624 (1981). K.-C. Zhang, Y.-F. Li, G.-B. Liu, and Y. Zhu, PhysicsLetters A , 1898 (2012). Y. Motome and N. Furukawa, Journal of thePhysical Society of Japan , 3186 (2001),https://doi.org/10.1143/JPSJ.70.3186. H. Shen, J. Liu, and L. Fu, Phys. Rev. B , 205140(2018). Y. Nagai, M. Okumura, and A. Tanaka, Phys. Rev. B , 115111 (2020). T. M. Rusin and W. Zawadzki, Phys. Rev. B , 205201(2020). K. Hukushima and K. Nemoto, Journal of thePhysical Society of Japan65