Noncollinear magnetic ordering in the Shastry-Sutherland Kondo lattice model: Insulating regime and the role of Dzyaloshinskii-Moriya interaction
NNoncollinear magnetic ordering in the Shastry-Sutherland Kondo lattice model:Insulating regime and the role of Dzyaloshinskii-Moriya interaction
Munir Shahzad and Pinaki Sengupta
School of Physical and Mathematical Sciences, Nanyang Technological University, 21 Nanyang Link, Singapore 637371 (Dated: December 5, 2017)We investigate the necessary conditions for the emergence of complex, noncoplanar magneticconfigurations in a Kondo lattice model with classical local moments on the geometrically frus-trated Shastry-Sutherland lattice and their evolution in an external magnetic field. We demonstratethat topologically nontrivial spin textures, including a new canted flux state, with nonzero scalarchirality arise dynamically from realistic short-range interactions. Our results establish that a fi-nite Dzyaloshinskii-Moriya (DM) interaction is necessary for the emergence of these novel magneticstates when the system is at half filling, for which the ground state is insulating. We identify theminimal set of DM vectors that are necessary for the stabilization of chiral magnetic phases. Thenoncoplanarity of such structures can be tuned continually by applying an external magnetic field.This is the first part in a series of two papers; in the following paper the effects of frustration,thermal fluctuations, and magnetic field on the emergence of novel noncollinear states at metallicfilling of itinerant electrons are discussed. Our results are crucial in understanding the magneticand electronic properties of the rare-earth tetraboride family of frustrated magnets with separatespin and charge degrees of freedom.
I. INTRODUCTION
The study of strongly interacting quantum many-bodysystems with independent spin and charge degrees of free-dom on frustrated lattices has attracted heightened in-terest in the recent past. The interplay between geo-metric frustration and strong interaction between itin-erant electrons and localized moments in these systemsresults in novel quantum phases and phenomena thatare not observed in their nonfrustrated counterparts .These competing interactions, together with crystal elec-tric fields and coupling to the itinerant electrons, oftenstabilize noncoplanar ordering of these moments .When an electron moves through such background spintextures, it picks up a Berry phase which underlies severalnovel transport phenomena such as the topological (or ge-ometric) Hall effect and unconventional magnetoresistiveproperties . The interest in these systems is drivenby the desire both to understand the underlying mech-anism driving the novel phenomena and to control theiremergence by external tuning fields in order to harnesstheir unique functionalities for practical applications.In this paper, we study the Kondo lattice model(KLM) on the geometrically frustrated Shastry-Sutherland lattice (SSL) with classical spins where thestandard (antiferromagnetic) Heisenberg interaction be-tween the local moments is supplemented by an ad-ditional Dzyaloshinskii-Moriya (DM) interaction. TheSSL is a paradigmatic geometry to study the effects ofcompeting interactions in the presence of frustration .The Shastry-Sutherland Kondo lattice model (SS-KLM)has previously been studied with S = 1 / , where quantum fluctuations of the local mo-ments play a crucial role in determining the character ofthe ground state. In the present study, we revisit thismodel, but with the local moments treated as classicalspins. This is not simply of academic interest. There ex- ists a complete family of rare-earth tetraborides ( R B , R =Tm, Er, Ho, Dy), quasi-two-dimensional metallicmagnets in which the magnetic-moment-carrying rare-earth ions are arranged in a SSL in the layers. Due tostrong spin-orbit coupling, the rare-earth ions in thesecompounds carry large magnetic moments and, conse-quently, can be treated as classical spins. They act as ef-fective local fields that interact strongly with the electronspin . In this study our goal is to construct a mini-mal model where topologically nontrivial chiral magneticphases can be realized from physically relevant interac-tions and investigate their evolution in an external mag-netic field. In particular, we explore the role of differentcomponents of the DM interaction in stabilizing differentaspects of the local-moment configurations. What are theminimal DM vectors required to stabilize a tunable non-coplanar spin configuration? How does an applied fieldaffect the noncoplanarity of the spin configuration? Doesthe nature of the chiral spin state change in the presenceof an external field? These are some of the questions weaddress in this work. Our results reveal that multiplenoncoplanar spin arrangements (characterized by differ-ent values of the scalar spin chirality) with long-rangemagnetic order are stabilized over an extended range ofparameters. Not surprisingly, we find that DM interac-tions play a crucial role in stabilizing chiral spin configu-rations. Furthermore, we are able to tune the noncopla-narity (equivalently, the topological character) of the spintextures, changing and suppressing the net chirality, byapplying an external magnetic field. This is in contrastto most previous studies in which the noncoplanar tex-tures of the local moments were imposed by extraneousfactors (e.g., crystal electric field in pyrochlores) and assuch cannot be changed easily. In an accompanying pa-per , we follow this up by studying the role of thermalfluctuations, frustration, and magnetic field in stabilizingthe noncollinear magnetic states and phase transitions a r X i v : . [ c ond - m a t . s t r- e l ] D ec at one-quarter and three-quarter filling of itinerant elec-trons, for which the ground state is metallic. Moreover,we found out that unlike the insulating case (discussed inthis paper), the noncollinear textures emerge for metal-lic densities without a DM interaction between the localmoments. II. MODEL
The Hamiltonian describing the SS-KLM with addi-tional DM interactions is given byˆ H = ˆ H e + ˆ H c , (1)where ˆ H e represents the electronic Hamiltonian,ˆ H e = − (cid:88) (cid:104) i,j (cid:105) ,σ t ij ( c † iσ c jσ + H.c.) − J K (cid:88) i S i · s i . (2)The first term is the kinetic energy of the itinerantelectrons; (cid:104) i, j (cid:105) represents the Shastry-Sutherland bonds(viz., first neighbors along the principal axes and thesecond neighbors along the alternate diagonals), and t ij are the transfer integrals for these bonds. The sec-ond term is the on-site Kondo-like interaction betweenthe spin of the itinerant electrons s i and localized mo-ments S i . The conduction electron spin is defined as s i = c † i,α σ αβ c i,β , where σ αβ is the vector element of theusual Pauli matrices. As mentioned in the Introduction,we treat the localized spins as classical vectors with unitlength ( | S i | = 1). In this limit, the sign of J K (ferro-magnetic or antiferromagnetic) is irrelevant since eigen-states that correspond to different signs are related bya global gauge transformation. The states of the local-ized spins are specified by the angular components as S i = (sin θ i cos φ i , sin θ i sin φ i , cos θ i ). The second partof the Hamiltonian (1) represents the classical localizedspin part: ˆ H c = ˆ H ex + ˆ H DM + ˆ H H , (3)where ˆ H ex expresses the classical Heisenberg interactionbetween the localized spins, ˆ H ex = (cid:80) (cid:104) i,j (cid:105) J ij S i · S j . ˆ H DM describes the DM interaction, ˆ H DM = (cid:80) (cid:104) i,j (cid:105) D ij . ( S i × S j ) , where D ij is the DM vector which is determined bythe crystal symmetry of the lattice. The precise values(directions and magnitude) of DM vectors will depend onthe details of the crystal symmetry of each compound.In this study, we choose a generic set of DM vectors andidentify the minimal interactions that are necessary forstabilizing noncoplanar spin textures. The explicit formof the DM vectors on the different bonds is given in thecaption of Fig. 1. The last term in the Hamiltonian (3) isthe Zeeman term for the localized spins due to an exter-nal (longitudinal) magnetic field, ˆ H H = − h z (cid:80) i S zi . A Zeeman term for the itinerant electrons is not includedexplicitly since the instantaneous spin orientation of theelectrons is determined completely by the local momentsin the large J K limit that we consider in this study. Here-after, the parameters with primes represent the interac-tions on diagonal bonds, while the unprimed ones referto the axial bonds. (cid:74) (cid:78)(cid:74)(cid:78) (cid:78) (cid:78)(cid:78)(cid:74) (cid:74) (cid:74) (cid:78) (cid:74) y x (cid:74) z (cid:74) D ⊥ out-of-plane (cid:78) D ⊥ into-plane D (cid:107) D (cid:48) FIG. 1. (Color online) DM interaction defined on the unitcell of SSL where the direction of the arrow from site i to site j indicates the direction of cross product S i × S j . The redarrows represent the parallel components of D , while (cid:74) and (cid:78) represent the out-of-plane and into-plane components of D . Blue arrows indicate the components of D (cid:48) on the diagonalbonds. The directions of x , y , and z axes are also indicated. III. METHOD AND OBSERVABLES
To investigate the above model, we use an unbiasedMonte Carlo (MC) method that has been used previouslyin the study of similar models . A brief review ofthis method is presented here closely following Refs. 34and 35. The dynamics of large localized moments of therare-earth ions is slow compared to itinerant electrons,and accordingly, we can decouple their dynamics fromthat of the itinerant electrons. While studying the lat-ter, we treat the local moments as static classical fieldsat each site. The electronic part of the Hamiltonian isbilinear in fermionic operators. Using the single-electronbasis, ˆ H e can be represented as a 2 N × N matrix for afixed configuration of classical localized spins, where N is the number of sites.In order to explore the thermodynamic properties, wewrite the partition function for the whole system by tak-ing two traces, Z = Tr S Tr I exp {− β [ ˆ H e ( { x r } ) − µ ˆ n e ] }× exp[ − β ( ˆ H c )] , (4)where Tr S and Tr I represent the traces over the classicallocalized spins denoted by { x r } and the itinerant-electrondegrees of freedom, respectively. The trace over itinerant-electron degrees of freedom can easily be calculated bythe numerical diagonalization of Hamiltonian matrix ˆ H e for a fixed configuration of localized spins { x r } ,Tr I exp {− β [ ˆ H e ( { x r } ) − µ ˆ n e ] }≡ (cid:89) ν (1 + exp {− β [ ε ν ( { x r } ) − µ ] } ) , (5)where µ is the chemical potential, β = 1 /k B T is the in-verse temperature, and ˆ n e = N (cid:80) iσ c † iσ c iσ is the numberdensity operator for conduction electrons. The partitionfunction for the whole system then takes the form Z = Tr S exp[ − S eff ( { x r } ) − β ( ˆ H c )] . (6)The corresponding effective action is S eff ( { x r } ) = (cid:80) ν F ( y ), where F ( y ) = − ln[1 + exp {− β ( y − µ ) } ]. Thegrand-canonical trace over localized spin degrees of free-dom is evaluated by sampling the spin configurationspace using a MC method. The probability distributionfor a particular configuration of localized spins { x r } canbe written as P ( { x r } ) ∝ exp[ − S eff ( { x r } ) − β ( ˆ H c )] . (7)The thermodynamic quantities that depend on local-ized spins are calculated by the thermal averages of spinconfigurations, while the quantities that are associatedwith itinerant electrons are calculated from the eigenval-ues and eigenfunctions of ˆ H e ( { x r } ). We start the sim-ulations with a random configuration of localized spins { x r } and calculate Boltzmann action S eff ( { x r } ) for thisconfiguration. The spin configuration is updated via theMetropolis algorithm based on the change in the effec-tive actions of the configurations resulting from randomupdates, ∆ S eff = S eff ( { x (cid:48) r } ) − S eff ( { x r } ). To identifymagnetic orderings we calculate the magnetization perunit site as well as the spin structure factor, which is theFourier transform of the spin-spin correlation function, S ( q ) = 1 N (cid:88) i,j (cid:104) S i · S j (cid:105) exp[ i q · r ij ] , (8)where r ij is the position vector from the i th to j th siteand (cid:104)·(cid:105) represents the thermal average over the grand-canonical ensemble. In order to describe the evolutionof the ground state under the influence of an externalmagnetic field, we calculate the magnetization per site,defined as m = (cid:118)(cid:117)(cid:117)(cid:116)(cid:42)(cid:18) (cid:80) i S i N (cid:19) (cid:43) . (9)To elucidate the difference between topological trivialand nontrivial states we evaluate the local scalar spin chirality. On a triangle the chirality is defined as χ (cid:52) = S i · ( S j × S k ) . (10)We use the total static spin chirality χ = N u (cid:80) (cid:52) χ (cid:52) (where the sum is over all the triangles formed on theplaquettes with diagonal bonds and N u is the numberof unit cells of SSL) as a quantitative measure of chi-ral order. This quantity is zero for collinear or copla-nar magnetic states such as ferromagnetic (FM), antifer-romagnetic (AFM), and pure flux states, whereas it isnonzero for noncoplanar configurations, e.g., all-out andthree-in–one-out states observed in pyrochlores. Finally,as an additional characterization of the chiral nature ofthe spin configurations, we measure the circulation ofthe in-plane components around each square plaquetteas f m = (cid:80) (cid:3) S i · r ij , where S i is the spin at site i and r ij are the vectors connecting sites i and j around the squareplaquette in a counterclockwise direction. A nonzero cir-culation identifies a flux configuration of the local mo-ments. IV. RESULTS
Simulations of the Hamiltonian (1) are performed inlattices of dimension L × L , with L = 8 −
16, over awide range of parameters. For smaller lattices, we usethe exact-diagonalization Monte Carlo (ED-MC) methodwhere the full Hamiltonian is diagonalized to calculatethe effective action for each MC step. For the larger lat-tices, we use the traveling cluster approximation (TCA)method : a 6 × T = 0 . T = 0 .
08. We repeat this process with a step of tempera-ture ∆ T = 0 .
02, finally calculating the thermal averagesof the observables at temperature T = 0 .
02. For the lat-tice sizes studied, the thermal gap to the lowest excitationis greater than the finite-size gap at T = 0 .
02. In otherwords, T = 0 .
02 is sufficiently small that ground-stateestimates of the measured observables can be reliably ob-tained. Measurements are done for 50000 MC steps after60 000 steps discarded in total for thermalization.With its multidimensional parameter space, the Hamil-tonian (1) is expected to support a rich array of ground-state phases over different ranges of the parameters. Inthe present work, we restrict our attention to the mag-netic behavior at an electronic filling factor (cid:104) n e (cid:105) = 1 / t = 1 .
0; this sets the energy scale for theproblem. The diagonal hopping is set to t (cid:48) = 1 .
2, and thevalues of the exchange interactions along the axial anddiagonal bonds are set at J = 0 . J (cid:48) = 0 .
12. Thischoice is motivated by the experimental observation ofapproximately equal bond lengths in the rare-earth tetra-boride family of compounds. In most materials of rele-vance to the present model, strong DM interactions exist.While the exact nature of DM interaction depends on thecrystal symmetries, we have chosen a generic form of DMinteraction for our study. Indeed, investigating the role ofDM interaction in stabilizing noncoplanar spin configu-rations is a central goal of the present study. Finally, fol-lowing the experimental observation in other frustratedmetallic magnets such as pyrochlores, the strength of theKondo coupling is chosen to be the strongest energy scalein the problem, J K = 8 . A. Effect of DM interaction
In the first part of the study, a systematic variationof the different components of the DM vectors is per-formed to identify the minimal set of vectors necessaryfor noncoplanar configurations of the local moments. Westudy the effects of the DM vectors normal to the planeof the lattice D ⊥ and two in-plane components, D (cid:107) and D (cid:48) (Fig .1). It is found that while D ⊥ is essential forthe emergence of noncollinearity, D (cid:48) and D (cid:107) stabilize thenoncoplanarity of the ground state.
1. Role of D ⊥ We start our discussion by analyzing the nature of themagnetic ground state at zero field. The evolution of theground state as D ⊥ is varied is shown in Fig. 2(a), wherewe have plotted the magnitude of the peaks in the spinstructure factor at q = ( π, π ) and ( π, D ⊥ = 0,the ground state has predominantly longitudinal AFMorder as the magnitude of the peak at ( π, π ) is almost 1and the value of the magnetization is almost zero. Withincreasing the value of D ⊥ , the ground state remains inthe same phase up to a critical value D c ⊥ ≈ .
07, beyondwhich there is a discontinuous transition to a state char-acterized by a large magnitude of the peak at ( π, π,
0) and (0 , π ). Nominally, suchfeatures in the structure factor point towards a cantedAFM state. However, the strong Kondo-like interactionprecludes such ordering in the double-exchange model, asdescribed in Ref. 40. The true nature of the ground stateis illustrated by a snapshot of the real-space (periodic)equilibrium spin configuration obtained from the simula-tions and shown schematically in Fig. 2(b). The in-planecomponents of the local moments are arranged in a near- ⊥ S ( π , π ) / N L = 8 L = 12 L = 16 S ( π , ) / N L = 8 L = 12 L = 16 (a) Ù ← FIG. 2. (Color online) (a) Evolution of the magnitude ofthe peaks in the static structure factor at q = ( π, π ) and( π,
0) as a function of D ⊥ for 8 ×
8, 12 ×
12, and 16 × × D ⊥ = 0 .
11, illustrating aflux state with a vanishingly small chirality. The results areobtained at T = 0 .
02 while keeping t = 1 . t (cid:48) = 1 . J = 0 . J (cid:48) = 0 . J K = 8 . D (cid:107) = 0 . D (cid:48) = 0 .
0, and h z = 0 . ideal flux pattern along vanishingly small chirality; thatis, the magnetic ground state is a noncollinear flux state.The net magnetization and the chirality remain vanish-ingly small across the range of D ⊥ studied, further con-firming the coplanar character of the flux state. Suchcomplex spin textures are essential ingredients for theobservation of unusual transport and electronic phenom-ena as noncollinear and noncoplanar magnetic orderingof localized spins behave like emergent electromagneticfields. χ L = 8L = 12L = 16 ′ m / m s (a)(b) FIG. 3. (Color online) (a) Chirality per unit cell and (b)magnetization per unit site as a function of D (cid:48) for 8 × ×
12, and 16 ×
16 lattice sizes. The MC results are obtainedat T = 0 . t = 1 . t (cid:48) = 1 . J = 0 . J (cid:48) = 0 . J K = 8 . D (cid:107) = 0 . D ⊥ = 0 .
11, and h z = 0 .
2. Role of D (cid:48) After finding the minimum value of D ⊥ that causesthe phase transition to a noncollinear state, we discussthe effect of D (cid:48) on stabilizing the noncoplanarity of theground state. The results for chirality and magneti- zation per unit site as a function of D (cid:48) are shown inFigs. 3(a) and 3(b). For weak D (cid:48) ( | D (cid:48) | (cid:46) . D (cid:48) , the data are converged with system size and con-vincingly indicate a nonzero chirality for the magneticground state. A strong value of D (cid:48) results in an enlargedout-of-plane component of the localized spins making theground state more canted. With increasing D (cid:48) the non-coplanarity of the magnetic ordered state increases mono-tonically [see Fig. 3(a)]. The same effect is observed inthe behavior of the magnetization per unit site: m/m s decreases with increasing system sizes at weak D (cid:48) but isfinite (and nonzero) for | D (cid:48) | (cid:38) .
03 and increases mono-tonically with D (cid:48) . The enlarged out-of-plane componentof the spins contributes to an increase in zero-field mag-netization [shown in Fig. 3(b)]. The static spin structurefactor S ( q ) at D (cid:48) = − .
10 (not shown here) exhibits asubdominant peak at q = (0 , D (cid:48) ; we call thisa canted flux state. It is worth mentioning that D (cid:48) can-not induce a phase transition to a noncollinear flux stateon its own; one always needs a nonzero value of D ⊥ forthat. In contrast to pyrochlores in which the tetrahedralordering of the local moments is fixed by the crystal-field effects, the canted flux state in our study arises dy-namically from competing interactions in the presence ofgeometric frustration. This enables us to control thesecomplex magnetic orderings continually via an externalmagnetic field.
3. Role of D (cid:107) ( h z = 0) The other in-plane component of the DM vector D (cid:107) simply reinforces noncoplanar configurations driven by D ⊥ (increased χ ) and increases the uniform magnetiza-tion by enhancing the canting of the local moments out-of- xy plane. B. Effects of an external magnetic field
One of the most intriguing features of the canonical(purely magnetic) Shastry-Sutherland model is the ap-pearance of magnetization plateaus in an applied mag-netic field. DM interactions are expected to stronglymodify the plateau structure. Our results show thatmagnetization plateaus are completely suppressed in therange of simultaneous DM and Kondo-like interactionsstudied in this work. For D ⊥ ≥ .
08 (which is neces-sary to stabilize the noncoplanar spin textures that weare interested in) and in the absence of D (cid:107) , the magne-tization increases monotonically all the way to satura- tion. The ground state is in a canted flux state at zeromagnetic field. The peak at ( π,
0) in the static struc-ture factor, S ( π, /N , decreases continuously to zero ata saturation field strength of h zs ≈ .
5. With increasingfield, the canting of the local moments (which is nonzerobut small at zero field due to the DM interaction) in-creases continuously until the local moments are fullypolarized. Interestingly, the approach to saturation isdistinct from that of a standard two-dimensional Heisen-berg AFM in the absence of an increased slope in themagnetization curve, further underscoring the difference m / m s D || = 0.0 D || = -0.10 D || = -0.15 h z χ (a) (b) FIG. 4. (Color online) (a) Magnetization per unit site and(b) chirality per unit cell as a function of external magneticfield for different values of D (cid:107) . The results are obtained forthe 12 ×
12 lattice at T = 0 . t = 1 . t (cid:48) = 1 . J = 0 . J (cid:48) = 0 . J K = 8 . D ⊥ = 0 .
11, and D (cid:48) = − . of the canted flux state from a conventional canted AFMstate. For nonzero D (cid:107) , there is a sharp increase in magne-tization upon the application of a small longitudinal field( h z (cid:46) .
1) to a state with a finite m/m s that depends onthe strength of D (cid:107) . Whether or not this happens via adiscontinuous transition at h z = 0 is not clear from ourresults. Upon further increasing the longitudinal field,the magnetization increases monotonically up to satura-tion, with the saturation field increasing with increasing D (cid:107) . The nature of the magnetic ground state remains acanted flux state throughout the field range. It remainsto be seen if magnetization plateaus can be stabilized forany range of DM interaction and coupling between thelocal moments and the itinerant electrons.Finally, we discuss the topological nature of the mag-netic ground state as it is tuned by an external field.Figure 4(b) shows the evolution of the net static spin chirality with increasing magnetic field for a representa-tive set of DM vectors where the zero-field ground stateis a canted flux state. With the increase of D (cid:107) the cant-ing of the localized spins increases even at zero field. Theground-state spin texture remains noncoplanar in natureover the entire range of applied field strength up to sat-uration. The chirality increases monotonically up to anintermediate value of the applied field and then decreasescontinuously to zero at saturation. The change in chiral-ity is simply driven by the increasing canting of the spinsalong the direction of the applied field (and a subsequentdecrease in the magnitude of the in-plane component ofthe local moments). Once again, a snapshot of local spinconfigurations elucidates the true nature of the magneticground state in an applied field [Figs. 5(a) and 5(b)]. Thein-plane components of the local moments are arrangedin a flux pattern on alternating plaquettes, whereas thereis a net canting of the longitudinal component parallel tothe applied field. The transition to saturation is markedby the complete breaking of the flux pattern driven bythe polarization of all the spins in the direction of mag-netic field. The qualitative features of Figs. 5(a) and 5(b)can be quantified partially in terms of the circulation ofthe in-plane components around each square plaquette, f m = (cid:80) (cid:3) S i · r ij . A nonzero circulation identifies a fluxconfiguration of the local moments. Figures 5(c) and 5(d)present plots of the local circulation for the same sets ofparameters as in Figs. 5(a) and 5(b). For h z = 0 . h c ≈ . V. SUMMARY
To summarize, we have studied the Kondo latticemodel with additional DM interaction on the Shastry-Sutherland lattice. Our results show that complex, non-coplanar spin configurations can be generated dynami-cally from purely short range interactions and couplingto itinerant electrons. We conclude that DM interactionsare necessary for the emergence of chiral spin configura-tions when the electronic spectrum is gapped; that is,the system is in an insulating state. We have carefullyidentified the minimal DM vectors necessary for the sta-bilization of noncoplanar configurations of the local mo-ments. Furthermore, such noncoplanar structures canbe tuned continually by applying an external magneticfield. These results provide insight into the origin and na- −3−2−10123 −3−2−10123 (d)(c)
FIG. 5. (Color online) The real-space configurations of the localized spins for the 8 × h z = 0 . h z = 1 .
6, when local spins become almost polarized to the direction of magnetic field. Thecolor bars beside the top plots indicate the out-of-plane component of the spin vector S zi . (c) and (d) Snapshots depicting thecirculation of flux, clockwise or counterclockwise, on each plaquette for the same values of magnetic field. The MC simulationsare performed at T = 0 . t = 1 . t (cid:48) = 1 . J = 0 . J (cid:48) = 0 . J K = 8 . D (cid:107) = 0 . D ⊥ = 0 .
11, and D (cid:48) = − . ture of topologically nontrivial magnetic phases in metal-lic magnets. They will also be crucial in understandingthe magnetic and electronic properties of the rare-earthtetraboride family of frustrated magnets. ACKNOWLEDGMENTS
It is a pleasure to thank G. Alvarez and Y. Kato foruseful discussions on the development of the numericalmethod. This work is partially supported by Grant No.MOE2014-T2-2-112 from the Ministry of Education, Sin-gapore. I. Martin and C. D. Batista, Phys. Rev. Lett. , 156402(2008). K. Barros, J. W. F. Venderbos, G.-W. Chern, and C. D.Batista, Phys. Rev. B , 245119 (2014). H. Ishizuka and Y. Motome, Phys. Rev. B , 081105(2013). G.-W. Chern, SPIN , 1540006 (2015). Y. Kato, I. Martin, and C. D. Batista, Phys. Rev. Lett. , 266405 (2010). Y. Machida, S. Nakatsuji, Y. Maeno, T. Tayama, T. Sakak-ibara, and S. Onoda, Phys. Rev. Lett. , 057203 (2007). S. Nakatsuji, Y. Machida, Y. Maeno, T. Tayama, T. Sakak-ibara, J. van Duijn, L. Balicas, J. N. Millican, R. T.Macaluso, and J. Y. Chan, Phys. Rev. Lett. , 087204(2006). M. Udagawa and R. Moessner, Phys. Rev. Lett. ,036602 (2013). G.-W. Chern, Phys. Rev. Lett. , 226403 (2010). Y. Akagi and Y. Motome, J. Phys. Conf. Ser. , 012059(2011). Y. Akagi, M. Udagawa, and Y. Motome, Phys. Rev. Lett. , 096401 (2012). J. W. F. Venderbos, M. Daghofer, J. van den Brink, andS. Kumar, Phys. Rev. Lett. , 166405 (2012). R. Karplus and J. M. Luttinger, Phys. Rev. , 1154(1954). J. Ye, Y. B. Kim, A. J. Millis, B. I. Shraiman, P. Majum-dar, and Z. Teˇsanovi´c, Phys. Rev. Lett. , 3737 (1999). D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. ,1959 (2010). Y. Taguchi, Y. Oohara, H. Yoshizawa, N. Nagaosa, andY. Tokura, Science , 2573 (2001). B. S. Shastry and B. Sutherland, Phys. B+C (Amsterdam,Neth.) , 1069 (1981). P. Coleman and A. H. Nevidomskyy, J. Low Temp. Phys. , 182 (2010). B. H. Bernhard, B. Coqblin, and C. Lacroix, Phys. Rev.B , 214427 (2011). J. H. Pixley, R. Yu, and Q. Si, Phys. Rev. Lett. ,176402 (2014). L. Su and P. Sengupta, Phys. Rev. B , 165431 (2015). K. Siemensmeyer, E. Wulf, H.-J. Mikeska, K. Flachbart,S. Gab´ani, S. Mat’aˇs, P. Priputen, A. Efdokimova, andN. Shitsevalova, Phys. Rev. Lett. , 177201 (2008). K. Wierschem, S. S. Sunku, T. Kong, T. Ito, P. C. Can-field, C. Panagopoulos, and P. Sengupta, Phys. Rev. B , 214433 (2015). S. S. Sunku, T. Kong, T. Ito, P. C. Canfield, B. S. Shas-try, P. Sengupta, and C. Panagopoulos, Phys. Rev. B ,174408 (2016). T. Suzuki, Y. Tomita, N. Kawashima, and P. Sengupta,Phys. Rev. B , 214404 (2010). S. Mat’a, K. Siemensmeyer, E. Wheeler, E. Wulf, R. Beyer,T. Hermannsdrfer, O. Ignatchik, M. Uhlarz, K. Flachbart,S. Gabni, P. Priputen, A. Efdokimova, and N. Shitseval-ova, J. Phys. Conf. Ser. , 032041 (2010). T. Suzuki, Y. Tomita, and N. Kawashima, Phys. Rev. B , 180405 (2009). F. Iga, A. Shigekawa, Y. Hasegawa, S. Michimura, T. Tak-abatake, S. Yoshii, T. Yamamoto, M. Hagiwara, andK. Kindo, J. Magn. Magn. Mater. , e443 (2007). S. Michimura, A. Shigekawa, F. Iga, T. Takabatake, andK. Ohoyama, J. Phys. Soc. Jpn. , 024707 (2009). M. Shahzad and P. Sengupta, Phys. Rev. B , 224402(2017). H. Ishizuka and Y. Motome, Phys. Rev. Lett. , 257205(2012). H. Ishizuka and Y. Motome, Phys. Rev. B , 081105(2013). H. Ishizuka and Y. Motome, Phys. Rev. B , 085110(2015). Y. Motome and N. Furukawa, J. Phys. Soc. Jpn. , 3853(1999). N. Furukawa and Y. Motome, J. Phys. Soc. Jpn. , 1482(2004). S. Kumar and P. Majumdar, Eur. Phys. J. B , 571(2006). A. Mukherjee, N. D. Patel, C. Bishop, and E. Dagotto,Phys. Rev. E , 063303 (2015). S. Kumar and P. Majumdar, Phys. Rev. Lett. , 016602(2006). S. Kumar and P. Majumdar, Phys. Rev. Lett. , 136601(2005). M. Yamanaka, W. Koshibae, and S. Maekawa, Phys. Rev.Lett.81