Decision making by Berry's flux
DDecision making by Berry’s flux
Benjamin O. Sung and Michael J. Lawler
2, 1 Dept. of Physics, Cornell University, Ithaca, NY 14853 Dept. of Physics, Applied Physics and Astronomy, Binghamton University, Vestal, NY 13850 (Dated: November 5, 2018)Order by disorder is a decision making process for frustrated systems but often leads to simpleanswers. We study order by disorder in the kagome Kondo model known for its complexity seekingrich decision making capabilities. At half filling and large Kondo coupling to hopping ratio J K /t ,the full manifold of 120 o kagome ground states are degenerate at second order in t/J K . We showthis degeneracy lifts at sixth order when a fermion can hop around a hexagon and feel the Berryflux induced by a given spin texture. Using Monte-Carlo we then seek the ground state of this sixthorder Hamiltonian and find in a 4x4 unit cell system that a co-planar 2x4 unit cell order is selectedover the √ × √ q = 0 and cuboc1 states, a result that survives even in the thermodynamic limit.This state is selected for its SU(2) flux properties induced by the spin texture and is analogous tothe integer quantum Hall effect. Given the existence of numerous quantum Hall plateaus in othersystems, the existence of this 2x4 unit cell state suggests that complex decision making is possibleon the manifold of 120 o states and achievable in different kagome-Kondo models. I. INTRODUCTION
Order by disorder is a decision making process: dif-ferent ground states will be selected from the same de-generate manifold depending on the type of fluctuationsit experiences. The ground state with the highest en-tropy, for example, is selected when the thermal order bydisorder effect is active. The ground state with the leastzero point energy, on the other hand, is selected whenthe quantum order by disorder effect is active. These ex-amples are pedagogically discussed in a recent review .In effect, the system computes the ground state accord-ing to the active order by disorder mechanism. In thisway, the order by disorder effect fits into the subject ofcomplexity science .Yet, in most cases, order by disorder, either thermal,quantum or another mechanism, selects a simple groundstate. Naturally, this is expected in the simplest casessuch as the J1-J2 model on the square lattice . Butfor kagome antiferromagnets the coplanar √ × √ and quantumfluctuations . This state has a tripled unit cell but isamong the simplest in the massively degenerate kagomeground state manifold. Even quantum order by disorderin pyrochlore Heisenberg antiferromagnets selects a col-inear state (though which colinear state is not clear atpresent ). Again, this is simpler than a general state inthe massively degenerate manifold would suggest. So or-der by disorder seems to act as a de-complexifying mech-anism.In this light, the discovery of order by disorder in clas-sical Kondo lattice models on highly frustrated lattices is interesting. Complex orders, some of which have mul-tiple wave vectors, are non-coplanar and incommensu-rate, arise in these models on the square lattice , cubiclattice , triangluar lattice , between the triangular andkagome lattices and kagome lattice . Further, theorder-by-disorder effect on the kagome lattice model does not select either the √ ×√ q = 0 but possibly the cuboc1 non-coplanar 120 state . So it seems possi-ble that order by disorder due to the fermion hopping inthese models may give rise to complex selection amonga highly degenerate ground state manifold and that it isnot de-complexifying.Given the potential for complex orders, order-by-disorder in classical Kondo models could also be inter-esting should it produce an integer quantum Hall effect.This is possible and indeed has provided much ofthe motivation for the study of these models . In 2Dan effect is particularly expected at finite temperatureshould the complex order have a non-zero scalar spinchirality . So, if the finite spin chirality cuboc1state were the selected state in the kagome case order-by-disorder would provide a mechanism for the stabilityof a state with an integer quantum Hall effect.In this paper, we revisit the order by disorder problemin the kagome Kondo lattice model with classical spinsat half filling. This problem is characterized by a smallparameter t/J K , a gap to electronic excitations and anSU(2) flux variable U felt by the electrons as they hoparound in a background classical spin texture. By carry-ing out a perturbative expansion in t/J K to sixth order,we show that the selection of a 120 o state is due to theflux felt by an electron as it hops around a hexagon. Us-ing Monte-carlo we then show that the state selected in a4 x x √ × √ q = 0 state:it has a spin origami sheet that is fully folded inone direction and perfectly flat in the other. We haveverified that it beats the √ × √ q = 0 and cuboc1state in the thermodynamic limit. Further, its SU(2)flux properties are also special: yielding energetic bene-fits both for hopping around hexagons and on bow-ties(pairs of triangles). Finally, we have computed the elec-tronic band structure and verified that the absence of aninteger quantum Hall effect as expected due to the van- a r X i v : . [ c ond - m a t . s t r- e l ] D ec ishing scalar spin charility. We conclude with an outlookon how these results may generalize to enable selectionof other complex ordering patterns within the kagome120 o states and thereby achieve complex decision mak-ing among this manifold of states. II. MODEL
We begin with the Kondo Hamiltonian given by H = H hop + H kondo = − t (cid:88) (cid:104) ij (cid:105) c † iσ c jσ − J k (cid:88) i S i · s i (1)Here, the kondo term in the above Hamiltonian gives thecoupling between the on-site classical spin vector S i withthe local electron spin s i = c † i τ c i with coupling constant J k . The first term above describes the electron hoppingalong adjacent sites on the given lattice with amplitude t . It is well known that in the double exchange model, i.e.the limit J k → ∞ , the ferromagnetic state dominatesover all other spin states at all but half filling. There,instead of ferromagnetism, a gap in the fermion spectraopens up and antiferromagnetism is found . The leadingantiferromagnetic term turns out to be just a nearestneighbor Heisenberg model.A quick exploration of the landscape of magnetism iseasily achieved using a variational calculation that takesinto account just a few possible states: the ferromagnetic, q = 0, and q = √ × √ ×
24 kagome lattice with periodic boundaryconditions, we obtain phase diagram of Fig. 1.In fig. 1, we have only portrayed the simplest of statesof many well-known states on the Kagome lattice witha focus on the general competition between ferromag-netism and antiferromagnetism. The figure suggests thatfor J k t (cid:29)
1, all states tend to converge to the same energyat half filling in the double exchange model. However,one can slightly weaken this condition, and only considerthe approximation J k t (cid:29) J k t , all 120 o degree states are degenerate .In addition to q = 0 and √ × √ o spin configurations, bothcoplanar and non-coplanar, that one may place on an ar-bitrarily large, finite kagome lattice, and to second orderin perturbation theory are all degenerate in energy. Inorder to extract a true, unique ground state out of thedegenerate 120 o spin state manifold, we must proceed tohigher orders in perturbation theory. III. VARIATIONAL RESULTS
In this paper, we will show that the degeneracy is liftedexactly at sixth order and exhibit a newly found state
FIG. 1. Simplified variational phase diagram of kagomeKondo model with classical spins. Green: Ferromagnetic,Blue: q = 0, Red: q = √ × √
3. Over 1728 sites wherethere are 24 triangles along the x-axis and 24 triangles alongthe y-axis. The horizontal axis represents the fraction of elec-trons occupying a site. The vertical axis is the ratio J k t . with a 4 × o states. Here, we canbegin to understand this result with numerical evidencethat the degeneracy is lifted at sixth order.Our numerical argument follows by presenting theenergy data for three well-known 120 o states. Theenergy data for these spin configurations were obtainedby exact diagonalizing the hamiltonian given in equa-tion 1 and summing the lower half of the eigenvaluespectrum, corresponding to half filling. On an 24 × J k are as in the table below: J k q = 0 q = √ × √ J k , the energy data for thethree 120 o states differ slightly, due to non-trivial sub-leading terms in higher-order perturbation theory. How-ever, as we increase J k , terms of order n in perturbationtheory are suppressed by the factor J n − k , and hence theenergies begin to converge. By J k /t = 5, the cuboc1state is winning but only in the 5th significant figure. FIG. 2. Graph of J k vs Energy. The energy values from J k = 1 to 100 are plotted. The blue dots correspond to thedifference between the states q = 0 and q = √ × √ q = √ × √ J k as J k goesfrom 1 to 100, and the y-axis are the energy values. Thelinearity of the two plots shows that the degeneracy breaksfirst at sixth order in perturbation theory. Energy calculationswere run on a 24x24 kagome lattice. We now extend this numerical evidence for degener-acy splitting among the 120 o states at sixth order. Byacquiring energy data as in the above for the spectrum J k = 1 to 100, and taking the differences in energies, weobtain figure 2 and clear evidence that the degeneracy islifted at sixth order.Having shown numerically that degeneracy breaksat sixth order, one may calculate the contributionsfrom each order in perturbation theory via numericalfitting. However, due to noise, it is nearly impossi-ble to calculate the numerical coefficients to higherorders above 4. Further, even though we expect thatthe coefficients at third and fifth order vanish, thefifth order term as calculated numerically suggestsa probable non-trivial coefficient. By carrying out anasymptotic fitting of the energy data we obtain the table:Order of J k Numerical Coefficients1 -2591.732 45.49793 2751.64 -233.15This numerical data would not be sufficient to determinethe exact formula for any contributions from higher or-ders in perturbation theory. But it does provide us witha check on our analytic calculation below.
IV. FEYNMAN DIAGRAM APPROACHA. Roadmap of Feynman Diagram Calculation
Here, we provide a roadmap for our feynman diagramcalculation. We proceed through the canonical method,writing down the path integral Z using our free Hamil-tonian eqn 1 without the hopping term. As usual, com-pleting the square yields the propagator for the free the-ory. Adding the interaction (hopping term) and tak-ing functional derivatives yields the full propagator forour theory. By taylor expanding in the perturbation H hop , we obtain the interaction U ij , which may be readilycomputed via unitary diagonalization as a 2 × B. Derivation of Feynman Rule for Propagator
We begin our diagrammatic procedure by deriving thefeynman rules for our theory using the canonical ap-proach. Proceeding via path integral quantization withgrassman variables, we take the kondo part of our hamil-tonian H kondo = − J k (cid:88) (cid:104) ij (cid:105) S i · s i and obtain Z = (cid:90) DcD ¯ c exp ( i (cid:90) dt [ (cid:88) i i ¯ c iσ ∂ t c iσ − H ])= (cid:90) DcD ¯ c exp ( i (cid:90) dt [ (cid:88) i i ¯ c iσ ∂ t c iσ + J k (cid:88) i S i · ¯ c iσ τ σσ (cid:48) c iσ (cid:48) + (cid:88) i (¯ η iσ c iσ + ¯ c iσ η iσ )) (2)To calculate the propagator for our free theory, weFourier transform to momentum space and complete thesquare, obtaining the green’s function given below G ( k, ω ) = 1 iω + J k τ − i(cid:15) (3)Calculating the propagator for the interacting theory fol-lows as follows: we consider the hopping term H hop = − t (cid:88) (cid:104) ij (cid:105) tc † iσ c jσ Taking functional derivatives of the path integral withrespect to the grassman valued sources ¯ η iσ and η iσ allowsus to make the replacement H i → (cid:88) (cid:104) ij (cid:105) δδη iσ U σσ (cid:48) ij δδ ¯ η jσ (cid:48) (4)where U σσ (cid:48) ij describes the hopping between nearest neigh-bor sites. This gives us the full green’s function G σ σ j i = δ σ σ j i ω − J k τ + i(cid:15)sgn ( ω ) (5)where the upper indices denote spin indices, and thelower indices denote points in position space. C. Derivation of Feynman Rule for InteractionVertex
Now that we have calculated the propagator for ourtheory, we simply need to determine the feynman rulefor our interaction vertex. To do this, note that we mustcalculate the amplitude for hopping between two sitesgiven by some unitary 2 × U ij , corresponding toup and down spin states. We begin with equation 1 inthe matrix representation H = (cid:18) − J k S z − J k S x − iJ k S y − J k S x + iJ k S y J k S z (cid:19) (6)Diagonalizing as usual, we obtain the unitary matrices U = (cid:32) − S x − iS y √ S z S x + iS y √ − S z S z +1 √ S z − S z +1 √ − S z (cid:33) , U † = (cid:32) − S x + iS y √ S z S z +1 √ S z S x − iS y √ − S z − S z +1 √ − S z (cid:33) (7)Expressing the hopping H hop in terms of unitary matri-ces, we obtain H = − t (cid:88)
In full, the above calculation gives us the followingfeynman rules: G σσ (cid:48) ij i, σ j, σ (cid:48) U σσ (cid:48) ij i, σ j, σ (cid:48) To calculate an n th order process, the Taylor expansiontells us to connect n G’s and n U’s in an alternatingorder. Intuitively, one may think of the propagator asthe virtual electron hopping between nearest neighbors,and the interaction U as imposing the constraint that thepath traversed by the virtual electron must be connected.Before moving to sixth order, we give a sample calcula-tion to first and second order. We begin with first order,to check that our path integral does indeed give zero foronly one hopping. i, σ j, σ (cid:48) The expression to calculate the energy is given by E = − iT ln ( Z ) − iT ln ( (cid:104) exp ( − i (cid:90) dtH I ) (cid:105) )= E (0) G − iT ln ( (cid:104) exp ( − iT (cid:90) dtH I (cid:105) (10)Hence, proceeding with the calculation for first order pro-cesses, we obtain E (cid:48) = (cid:90) dτ − t (cid:88) ij (cid:104) c iσ ( t ) U σσ (cid:48) ij c jσ (cid:48) ( t ) (cid:105) = tT r (cid:90) t t dτ (cid:88) ij G ij ( t − t ) U ij = tT (cid:88) ij (cid:90) dω π G ij ( ω ) U ij (11)Now, noting that the green’s function gives the deltafunction G ij ∝ δ ij and that U ij vanishes identically for i = j , we see that the correction vanishes.We now consider the integral to second order. In thiscase, we have the diagram i, σ i k, σ k j, σ j l, σ l Here, we have two vertical lines corresponding to twovertices on the kagome lattice. Although not indicated,we choose the convention with the direction of propa-gation to the right, hence, we may only connect, usinga propagator, a vertex at the left with a vertex at theright. For instance, we may connect ( k, σ k ) with ( j, σ j ),but not ( k, σ k ) with ( i, σ i ). Note that there are only twodistinct ways to connect using propagators as indicatedbelow: i, σ i k, σ k j, σ j l, σ l i, σ i k, σ k j, σ j l, σ l As shown, the second diagram is disconnected, and henceit does not contribute to the energy correction at secondorder. Evaluating the first diagram as usual, one canverify that we obtain the correction E (cid:48) = T r ( (cid:90) dω π (cid:88) ijkl G li U ij G jk U kl )= − J k (12)where the spin indices have been suppressed where thesecond line has been evaluated for the case of 120 o states.This answer agrees with our calculation using the unitarymatrices and the usual quantum mechanical perturbativeenergy formula.We may carry this procedure out to sixth order, withwhich we will obtain the expression E (cid:48) = T r ( (cid:90) dω π (cid:88) GU GU GU GU GU GU ) (13)where spatial and and spin indices have been suppressed.We remark that there is no need to connect all the possi-ble lines in a single six order feynman diagram since the energy sums over all possible spatial indices. Hence, withthese feynman rules, we may derive perturbative energycorrections to any order with a single feynman diagram.
E. Numerical Calculations
Clearly, it is necessary to numerically evaluate the an-alytical expression obtained from feynman diagrammatictechniques for energy corrections of order n >
2. To dothis, we separate equation 13 into a linear combinationof Pauli matrices as in the below.We define the variables µ = 12 ( 1 ω − J k + i(cid:15) + 1 ω + J k − i(cid:15) ) ν = 12 ( 1 ω − J k + i(cid:15) − ω + J k − i(cid:15) ) (14)and α = 12 ( 1 (cid:112) S x (cid:112) S y ( S x · S y − iS x S y + iS x S y + S x + S y + 1)+ 1 (cid:112) − S x (cid:112) − S y ( S x · S y − iS x S y + iS x S y − S x − S y + 1)) δ = 12 ( 1 (cid:112) S x (cid:112) S y ( S x · S y − iS x S y + iS x S y + S x + S y + 1) − (cid:112) − S x (cid:112) − S y ( S x · S y − iS x S y + iS x S y − S x − S y + 1)) β = 12 ( 1 (cid:112) S x (cid:112) − S y ( − S x · S y + iS x S y − iS y S x + S x − S y + 1)+ 1 (cid:112) − S x (cid:112) S y ( − S x · S y + iS x S y − iS y S x − S x + S y + 1)) γ = 12 ( 1 (cid:112) S x (cid:112) − S y ( − S x · S y + iS x S y − iS y S x + S x − S y + 1) − (cid:112) − S x (cid:112) S y ( − S x · S y + iS x S y − iS y S x − S x + S y + 1)) (15)and rewrite the green’s function and unitary matrices as1 ω − J k τ σ σ = µτ + µτ U σσ (cid:48) ij = ατ + βτ + δτ + γτ (16)Consequently, we may rewrite equation 13 by substitut-ing in these linear combinations, and evaluate it numeri-cally in mathematica. We find to third order and fifth or-der that the energy correction vanishes as expected fromthe numerical results.Intuitively, we may think of the numerical result thatdegeneracy is first lifted at sixth order as follows. Sec-ond order in perturbation theory selects out a degeneratemanifold of 120 o states. Since third order and fifth ordervanishes, the fourth order term is the only possible non-trivial contribution. However, if we think of the electronas hopping around paths with the condition that it be-gins and ends on the same vertex, it is easy to see thatthere are no trivial loops that could be achieved via fourhoppings. This is only viable at sixth order, in which theelectron may hop across bowtie loops and hexagon loops.In particular, we verify numerically using the above pro-cedure that fourth order terms are identical for the well-known 120 o states on the kagome lattice. We attributethe non-trivial sixth order contribution as resulting fromthe Berry’s phase.We further proceeded via a simulation on a 4 × q = 0 state at sixth order as expected.Having derived this expression, this motivates us to carryout a Monte-Carlo simulation to find a global minimumof this theory. In the next section, we detail the resultsthat we obtained using a non-linear local optimizationmethod. V. NUMERICAL NONLINEAR LOCALOPTIMIZATION
Using the expression derived above, we may numeri-cally evaluate the sixth order contributions from any spinconfiguration on the kagome lattice. To do this, we cre-ate an ensemble of 1000 random spin configurations onthe 4 × × × o condition, i.e.,neighboring on-site spin vectors must have an inner prod-uct of − using the NMinimize method in Mathematica.We then calculate the sixth order contribution of eachof these 120 o states. Upon doing this, we encountered anew state, with energy lower than any other well-knownstate, which we will detail below.By calculating the sixth order contributions from non-trivial loops on each of the states on the 4 × FIG. 3. Spin pattern plots and their associated hexagon fluxsthat contribute at sixth order. The colors correspond to thevalues of the fluxes through the loops, i.e. ”hotter” colors aregreater and ”colder” colors are smaller. Top sub-figure is the q = 0 state, middle sug-figure is the cuboc1 state and thebottom sub-figure is the ”snake” state. state” dominates over both the q = 0 and cuboc q = 0, q = √ ×√ cuboc ×
24 kagomelattice. The reader should note that this verification iscompletely independent of our perturbation theory. Thiswas calculated using only the Hamiltonian for the sys-tem and inputting the relevant data for the classical spin
FIG. 4. Spin origami plots that map a 120 o state to afolded sheet of paper. Top state is the q = 0 pattern that mapsto a flat sheet of paper. Middle is the √ ×√ × q = 0 and q = √ × √ vectors for the “snake” state.We now make one further remark regarding the spinplots displayed in figure 3. For coplanar 120 o states onthe kagome lattice, it can be shown that only hexagonalfluxes contribute to the breaking of degeneracies sincethe bowtie loops all contribute the same energy. As de-tailed in the table below, the contributions of hexag-onal fluxes is completely consistent with our observa-tions that at sixth order, we have the ordering of states: FIG. 5. The newly found state is called the ”snake” state dueto the above behavior of alternating colors. It is a coplanarstate with rgb colors corresponding to the usual spins on the q = 0 and q = √ × √ J k = 1 to 100 between the ”snake” and q = √ ×√ o states are broken, and our ”snake” state isthe winner. Energy calculations were run on a 24x24 kagomelattice. snake < cuboc1 < q = 0 < q = √ × √ q = √ × √ q = 0 0.158203 -0.251953 0.158203Snake 0.158203 -0.251953 -0.193359 VI. ELECTRONIC BAND STRUCTURE OFSNAKE STATE
Finally, we have computed the electronic band struc-ture of the snake state and its associated Chern number.The band structure for the filled bands is shown in 7. Wehave further computed the Chern number following Ref.28. We find both for the bottom group of 16 bands andthe top group of 8 bands (together making the 24 filledbands out of 48 bands present in the snake state) that the Γ X K Γ Y K98.599.099.5100.0100.5101.0101.5102.0 E ne r g y [ e V ] FIG. 7. Electronic band structure for the electrons hoppingin the background of magnetic ordering of the snake state.Here the Brillouin zone is rectangular with sides at the X andY points and corners at the K point.
Chern number is trivial. So order by disorder, at leastas determined within our 4x4 unit cell calculation, is notselecting a state that could support an integer quantumHall effect.
VII. CONCLUSION
In this paper, we explored the problem of state se-lection at half filling of the kagome Kondo model withclassical spins. We were motivated to derive an analyt-ical expression for higher orders in perturbation theoryin order to better understand the degenerate 120 o statemanifold and achieved this to sixth order.Proceeding via a numerical non-linear local optimiza-tion algorithm, we find, out of an ensemble of 1000 ran-dom 120 o states, a state that beats all the other well-known state as readily verified numerically via fig 6. Inparticular, combining the 120 o state minimization along with the sixth order expression in the algorithm, we findthat many of the runs readily converge to this state. Fur-ther work in this direction would include working with alarger ensemble with a more powerful machine in an at-tempt to find an even better ground state. In particular,we would like to fully understand the contribution of thefluxes to the sixth order correction. While we are drawnto the conclusion that the hexagon fluxes are responsi-ble for the relative correction among the coplanar 120 o states, we are not completely sure as to how the bowtieand hexagon fluxes contribute to the non-coplanar 120 o states.It is remarkable that the snake state we find has a unitcell with 24 spins (2x4 unit cells) but was found in acalculation with 48 spins (4x4 unit cells). This suggests,another state with an even larger unit cell may ultimatelywin the order by disorder competition. However, even ifthis is not the case, order by disorder due to fermionhopping and associated Berry flux has selected a 120 state that to our knowledge has never been consideredbefore. Further it has more spins in its unit cell eitherthe √ × √ o manifold of states.The answer is likely yes: one could increase the spin rep-resentation of the fermion degree of freedom from spin1 / S as in the study of Ref. 16 who alsoconsider the kagome lattice. This would enable the orderby disorder effect to occur at other fillings than 1/2 withpotentially different Berry flux desires. Spin represen-tation could therefore introduce a hierarchy of order-by-disorder mechanisms each possibly selecting a differentstate. So, order by disorder via fermion hopping withBerry flux could provide a rich set of decision makingcapabilities on the kagome 120 o manifold and other suchmanifolds common in highly frustrated magnetism. J. Villain, R. Bidaux, J.-P. Carton, and R. Conte, J. Phys., , 1263 (1980), ISSN 0302-0738. J. T. Chalker, in
Introd. to Frustrated Magn. , edited byC. Lacroix, P. Mendels, and F. Mila (Springer, 2011) pp.3–22. M. Mitchell,
Complexity: A Guided Tour (Oxford Univer-sity Press, 2009). P. Chandra and B. Doucot, Phys. Rev. B, , 9335 (1988),ISSN 01631829. C. L. Henley, Phys. Rev. Lett., , 2056 (1989), ISSN00319007. J. T. Chalker, P. C. W. Holdsworth, and E. F. Shender,Phys. Rev. Lett., , 855 (1992), ISSN 0031-9007. A. Chubukov, J. Appl. Phys., , 5639 (1993). C. Henley and E. Chan, J. Magn. Magn. Mater., ,1693 (1995), ISSN 03048853. C. L. Henley, Phys. Rev. Lett., , 1 (2006), ISSN00319007, arXiv:0509005 [cond-mat]. U. Hizi and C. L. Henley, Phys. Rev. B, , 1 (2009), ISSN1098-0121. S. Ghosh, P. O. Brien, C. L. Henley, and M. J. Lawler,arXiv:1407.5354 (2014), ISSN 1550235X, arXiv:1407.5354. R. Ozawa, S. Hayami, K. Barros, G.-W. Chern, Y. Mo-tome, and C. D. Batista, , 1 (2015), arXiv:1510.06830. K. Pradhan and P. Majumdar, EPL (Europhysics Lett., , 37007 (2009), ISSN 0295-5075. H. Ishizuka and Y. Motome, Phys. Rev. Lett., , 237207(2012), ISSN 0031-9007. Y. Akagi and Y. Motome, J. Korean Phys. Soc., , 405(2013). M. Udagawa, H. Ishizuka, and Y. Motome, , 6(2013), arXiv:1310.4580. K. Barros, J. W. F. Venderbos, G.-W. Chern, and C. D.Batista, Phys. Rev. B, , 245119 (2014), ISSN 1098-0121,arXiv:1407.5369. L. Messio, B. Bernu, and C. Lhuillier, Phys. Rev. Lett., , 1 (2012), ISSN 00319007, arXiv:1110.5440v1. J. Ye, Y. B. Kim, a. J. Millis, B. I. Shraiman, P. Majumdar,and Z. Teˇsanovi´c, Phys. Rev. Lett., , 3737 (1999), ISSN0031-9007, arXiv:9905007 [cond-mat]. K. Ohgushi, S. Murakami, and N. Nagaosa, Phys. Rev. B, , 4 (1999), ISSN 0163-1829, arXiv:9912206 [cond-mat]. N. Nagaosa, J. Phys. Soc. Japan, , 042001 (2006), ISSN0031-9015. I. Martin and C. D. Batista, Phys. Rev. Lett., , 156402(2008), ISSN 0031-9007. L. N. Bulaevskii, C. D. Batista, M. V. Mostovoy, and D. I.Khomskii, Phys. Rev. B - Condens. Matter Mater. Phys., , 024402 (2008), ISSN 10980121, arXiv:0709.0575. S.-S. Gong, W. Zhu, L. Balents, and D. N. Sheng, Phys.Rev. B, , 1 (2015), ISSN 1098-0121, arXiv:1412.1571. E. F. Shender, V. B. Cherepanov, P. C. W. Holdsworth,and A. J. Berlinsky, Phys. Rev. Lett., , 3812 (1993),ISSN 00319007, arXiv:9303031 [cond-mat]. P. Chandra, P. Coleman, and I. Ritchey, J. Phys. I, , 591(1993), ISSN 1155-4304. We chose the system size of a 4 × × T. Fukui, Y. Hatsugai, and H. Suzuki, J. Phys. Soc. Japan,74