Phase transitions in the classical exchange-anisotropic Kitaev-Heisenberg model
X. M. Zhang,1,*R. M. Liu,1,*Z. Jin,1T. T. Liu,1D. Y. Chen,1Z. Fan,1M. Zeng,1X. B. Lu,1 X. S. Gao,1M. H. Qin ,1,†and J.-M. Liu1,2
1Guangdong Provincial Key Laboratory of Quantum Engineering and Quantum Materials, and Institute for Advanced Materials, South China Academy of Advanced Optoelectronics, South China Normal University, Guangzhou 510006, China
2Laboratory of Solid State Microstructures, Nanjing University, Nanjing 210093, China
(Received 21 July 2020; revised 1 October 2020; accepted 10 October 2020; published 27 October 2020) The Kitaev model on the honeycomb lattice has been receiving substantial attention due to the discovery of quantum spin liquid state associated with this model. Consequently, its classical partners such as the Kitaev-Heisenberg (KH) model and associated phase transitions become concerned. Specifically, an intermediate Kosterlitz-Thouless (KT) phase engaged in the transition from the high-temperature (T) disordered state to the low-Tsixfold degenerate state is predicted in the isotropic KH model [Phys. Rev. Lett.109, 187201 (2012)], but so far no sufficient experimental proof has been reported. In this work, we consider an essential extension of this KH model on the honeycomb lattice by including the Kitaev exchange anisotropy that is non-negligible in realistic materials. The associated phase transitions are thus investigated using the Monte Carlo simulations. It is found that such an anisotropy will result in a degradation of the sixfold degeneracy of the ground state in the isotropic KH model down to the fourfold or twofold degenerate ground state, and the finite-Tphase transitions will also be modified remarkably. Interestingly, the intermediate KT phase can be suppressed by this Kitaev exchange anisotropy. This work thus provides a more realistic description of the physics ingredient with the KH model and presents a possible explanation on absence of the intermediate phase in real materials where the Kitaev exchange anisotropy can be more or less available.
DOI:10.1103/PhysRevE.102.042132
I. INTRODUCTION
During the past few years, the Kitaev model on the honey- comb lattice hosting a quantum spin liquid ground state has attracted a great deal of attention due to its important role in basic physical research and its candidate platforms including the iridatesA2IrO3 (A=Li, Na) and relatedα-RuCl3 [1–7].
Consequently, the corresponding classical Kitaev-Heisenberg (KH) model with the Heisenberg term available in real ma- terials is also concerned, which exhibits fascinating physics and emergent phenomena attributing to the honeycomb lattice symmetry and competition between the spin-anisotropic Ki- taev and isotropic Heisenberg exchange interactions [8,9].
For examples, the ground-state phase diagram in the exchange parameter space exhibits several distinct phases including: the ferromagnetic phase with spin configuration shown in Fig.1(a), the antiferromagnetic (AFM) stripy phase in Fig.1(b), the zigzag phase in Fig.1(c), and the AFM Néel phase in Fig.1(d) [10–13]. The appearance of these phases adds novel physics to such highly frustrated systems on one hand, and on the other hand they can be used to explain the emergent phenomena observed in real materials. One case is iridatesA2IrO3 whose magnetic properties should be de- scribed by this KH model with involved contributions of these phases [14–17].
*These authors contributed equally to this work.
Indeed, earlier calculations based on the isotropic ver- sion of this KH model, i.e., the Kitaev and Heisenberg exchanges are isotropic, suggest that the competition between the nearest-neighbor AFM Kitaev and ferromagnetic (FM) Heisenberg exchanges is responsible for stabilizing the zigzag state observed experimentally in Na2IrO3[10], while further neighboring and off-diagonal exchanges play an important role in determining the magnetisms of Li2IrO3 [4]. More- over, the long-range zigzag state is rather fragile and can be replaced by a spin-glass state due to tiny impurity doping [18,19], which has been reproduced by the Monte Carlo (MC) simulations on the modified KH model [20,21]. Both the dis- ordered Kitaev and Heisenberg exchanges caused by magnetic substitution are suggested to be responsible for the doping effects such as fragile zigzag state.
In addition, interesting phase transitions and critical be- haviors of the isotropic KH model have been theoretically uncovered, contributing a great deal to the development of statistical mechanics [22,23]. For example, the finite-size scal- ing analysis revealed that a critical Kosterlitz-Thouless (KT) phase with continuously variable exponents is sandwiched by a low-temperature (T) ordered phase and a high-Tdisordered phase [12,24], and this intermediate KT phase and its con- nection with experimental observations in realistic materials have been often discussed. The phase transition of this KH model is rather similar to that of the six-state clock model to some extent, attributing to the sixfold degeneracy of the ordered phases in the two models [25].
With continuous attention to the phase transitions and crit- ical behaviors for this isotropic KH model, strong interest
FIG. 1. Spin configurations in (a) the ferromagnetic order, (b) the stripy order, (c) the zigzag order, and (d) the Néel order. Solid and empty circles represent the up spins and down spins, respectively.
(X,Y,Z) denote the three bonds of the honeycomb lattice.
goes to the transition behaviors of some modified KH models in order to better fit realistic physics. Along this line, one such modification refers to the so-called anisotropic Kitaev exchange. In fact, it is noted that the isotropic Kitaev ex- change should be an ideal extreme for a realistic system, and commonly the exchange should be anisotropic, although the anisotropy may be strong or weak. For example, the exchange anisotropy may be induced inA2IrO3by tiny ion substitution and/or uniaxial strain. The honeycomb lattice distortion, as reported in Li-doped Na2IrO3, can also certainly result in the so-called Kitaev exchange anisotropy [26]. More importantly, so far no explicit evidence with the existence of this inter- mediate KT phase, to our best knowledge, has been reported experimentally, raising concern with the reason why this well- defined intermediate KT phase predicted in the isotropic KH model is absent in real materials. While the possible reason is connected with the impact of inevitable defects and inhomo- geneity, what is the underlying mechanism for such defects or inhomogeneity? Is the exchange anisotropy as induced by the inhomogeneity the possible reason? Indeed, an involvement of such Kitaev exchange anisotropy may change the ground state and the stability of the intermediate KT phase, and an extensive investigation of this anisotropic KH model becomes urgently needed.
In fact, the present work will show that the sixfold degener- acy of the ordered phase in the isotropic KH model would be degraded down to fourfold or even twofold degeneracy with increasing Kitaev exchange anisotropy. This change is ex- pectable and similar prediction was given in our earlier work on the exchange-anisotropy induced variation of the phase transitions in the anisotropic triangular Ising model [27]. In such sense, the phase transitions and critical behaviors of the exchange-anisotropic KH model should be checked, and in particular the stability of the intermediate KT phase should be discussed carefully, besides additional ingredients of physics to understand those experimentally observed but not yet well understood phenomena.
In this work, we start from such a classical KH model with Kitaev exchange anisotropy, i.e., the exchange-anisotropic KH (EAKH) model, and investigate the ground-state and relevant phase-transition behaviors in the parameter space.
A major outcome is the disappearance of the intermediate KT phase predicted by the isotropic KH model. The phase transition from the high-T paramagnetic state to the low-T fourfold degenerate state is straightforward at finite temper- ature, thus explaining the absence of such intermediate KT phase in experimental observation. One believes that this work may be a substantial step towards an understanding of more emergent physics with the classical KH model.
II. THE EAKH MODEL AND METHOD OF SIMULATIONS We first present the EAKH model as a generalized version of the original isotropic KH model. The ground states will be calculated in this section, while the phase diagram and phase transitions will be presented in Secs.IIIandIV, respectively, for two different schemes, i.e., the simple scheme and ex- tended scheme to be addressed below. These two schemes can be used to investigate various phase transitions and relevant critical properties.
A. The EAKH Model
The Hamiltonian of the proposed EAKH model can be written as
H= −2JK
⎛
⎝
i jX
SxiSxj+
i jY
SyiSyj+μ
i jZ
SziSzj
⎞
⎠ +JH
i j
SiSj, (1)
where the first term is the Kitaev interaction with subscripts (X,Y,Z) denoting the three bonds (depicted in Fig.1) of the honeycomb lattice, Si=(Six,Siy,Siz) represents the Heisen- berg spin with unit length on siteiand three components along the cubic axesx,y, andzof the IrO6octahedra. Factorμis the anisotropy factor (0< μ <2) measuring the Kitaev exchange anisotropy defined along the zaxis. The second term is the Heisenberg interaction between the nearest-neighbor spins [2]. Here it is noted that the Heisenberg exchange anisotropy may be considered too. Nevertheless, one can see that such an anisotropy hardly affects the degeneracy of these states and thus is not considered here for simplicity, although it could be also available in real materials.
It is known that the isotropic KH model exhibits the sixfold degenerate ground state, and accessing one of the six possible states at low T is attributed to the competition between the Kitaev exchange JK and Heisenberg exchangeJH, in which the spins are oriented along one of the cubic directions which is spontaneously chosen by the system. When the Kitaev exchange anisotropy is introduced, the sixfold degeneracy of the ground state is decreased to a fourfold or twofold one, depending on the value of μ, i.e., the spins in the ground state are in-plane or out-of-plane oriented. As a matter of fact, the ground state and corresponding phase boundaries can be determined by comparing the total energy terms of these phases as functions ofJK,JH, andμ. These calculations are
TABLE I. Internal energies as functions ofαandμin the possi- ble ground states.
Fourfold degenerate Twofold degenerate Néel AF state ExyN=2JK−3JH EzN =2μJK−3JH
Stripy state ExyS = −2|JK| −JH EzS = −2μ|JK| −JH
Zigzag state ExyZ = −2|JK| +JH EzZ = −2μ|JK| +JH
Ferromagnetic state ExyF = −2JK+3JH EzF = −2μJK+3JH
straightforward and the results are summarized in TableI. For the fourfold degenerate states, the energy terms per site are (ExyN, ExyS, ExyZ, and ExyF) for the Néel state, the stripy state, the zigzag state, and the ferromagnetic state, respectively, noting the superscripts (N,S,Z,F) standing for the four states and the subscriptxy for the in-plane spins. For the twofold degenerate states, these energy terms are (EzN,EzS,EzZ,EzF), re- spectively, noting the subscriptzstanding for the out-of-plane spins.
B. Monte Carlo simulation
Subsequently, we investigate the finite-Tphase transition and critical behaviors of the EAKH model based on the Monte Carlo (MC) simulations on various lattices with a total number of sites equal to N =2×L×L using the over-relaxation [28], heat-bath [29], and temperature-exchange methods [30]. Typically, the initial 5×105 MC steps (over-relaxation sweeps) are discarded for equilibrium consideration and an- other 5×105MC steps are retained for statistical averaging.
We take an exchange sampling after every 10 MC steps, and perform a heat-bath sweep per 10 MC steps.
Following the earlier work [24], the order parame- ters for describing these states can be defined, respec- tively, as MF =
i(SiA)/N for the ferromagnetic state [Fig.1(a)],MS =
i(SiA−SiB+SiC−SiD)/Nfor the stripy state [the four sublattices are shown in Fig. 1(b)], MZ =
i(SiA+SiB−SiC−SiD)/N for the zigzag state [Fig.1(c)], and MN =
i(SiA−SiB)/N for the Néel state [Fig. 1(d)].
Furthermore, the corresponding susceptibilities and Binder’s cumulants are also calculated to estimate the transition tem- peratures [31].
III. PHASE TRANSITIONS TO THE STRIPY STATE AND NÉEL STATE
Following the earlier work, our major attention first goes to a simple scheme for the exchanges, for the convenience of comparison and discussion. In this simple scheme, a ratio factor α(0< α <1) is introduced to measure the relative magnitudes of ferromagnetic exchange JK and AFM ex- change JH, satisfying relation JK =α and JH =1−α [3].
The EAKH model is thus discussed within this scheme, including the as-generated phase transitions in the (α, μ) pa- rameter space and the associated critical behaviors.
In the absence of exchange anisotropy, i.e., μ=1, the ground state favors the two-sublattice AFM Néel state for 0< α <1/3, and the AFM stripy state for 1/3< α <1.
Subsequently, the degeneracy of the ground state is mainly
FIG. 2. Ground-state phase diagram of the EAKH model for the simple scheme.
determined by parameter μ, as shown in Fig. 2 where the calculated ground-state phase diagram is presented. In details, for a fixed μ, the critical line α=1/(2+μ) differentiates the Néel state and the stripy state, noting that the competi- tion between the Kitaev exchange and Heisenberg exchange determines the ground state. Moreover, the degeneracy of the ground stripy (Néel state) stabilized forα >1/(2+μ)[α <
1/(2+μ)] is degraded down to the fourfold (twofold) de- generacy forμ <1 and the twofold (fourfold) degeneracy for μ >1. It is noted that the zigzag state and FM state cannot be stabilized within this simple scheme, while they can be the ground states within the extended scheme to be discussed in Sec.IV.
A. Phase transition to the stripy state
The stripy state occupies the region of α >1/(2+μ) in the phase diagram. Following the earlier work [24], the histogram method is used to investigate the discreteness of the order parameter, considering the fact that the projections preserve the cubic symmetry of the model and well reflect the orientation of the order parameter. The corresponding results are shown in Fig.3where the simulated projection patterns of the stripy order parameterMSon the (111) plane forα=0.55 are presented. It is seen that for μ=1, i.e., the isotropic cases, the sixfold peaked pattern is observed, as shown in Fig.3(a), demonstrating the sixfold degeneracy of the stripy phase with order parameter aligned along one of the cubic
FIG. 3. The low-T distribution of the projections of the stripy order parameter on the (111) plane forα=0.55 for (a)μ=1, and (b)μ=0.8.
FIG. 4. Binder cumulantBs as a function ofTfor different sys- tem sizeLforα=0.55 for (a)μ=1, and (b)μ=0.8. The log-log plots of the order parametermSas a function ofLlog(mS) vs log(L) at various temperatures for (c)μ=1, and (d)μ=0.8. The solid curves indicate the linear behavior that corresponds to a power-law dependencemS ∼L−η/2.
axes. For the anisotropic cases, e.g., by settingμ=0.8, the stripy order with the spins oriented along thezaxis cannot be stabilized anymore, resulting in the fourfold peaked pattern of the distribution function of stripy order, as clearly shown in Fig.3(b).
The critical temperatureTc for the stripy state, at which the long-range order is destroyed, can be estimated from the crossing of Binder’s fourth-order cumulantBsfor different lat- tice sizesL, noting that theBs(T) curves for differentLusually intersect at the critical point for continuous phase transitions.
The simulated results for μ=1 andμ=0.8 are plotted in Figs.4(a)and4(b), respectively. It is clearly shown that the critical point Tc significantly shifts toward the high-T side when the Kitaev exchange anisotropy is induced, attributed to the reduced degeneracy, notingTc=0.12±0.006 forμ=1 (isotropy) andTc=0.212±0.006 forμ=0.8 (anisotropy).
Subsequently, we investigate the critical behavior of the phase transition using the finite-size scaling analysis and pay particular attention to the cases forT >Tc, considering the fact that the long-range ordered state is developed at T ∼ Tc. Similarly, a two-dimensional complex order parameter, mS =6
j=1|Mj,S|exp(ιθj) withθj=πnj/3,nj=0, . . . ,5 is defined to characterize the six possible ordered states (the or- der parameters of the Néel, zigzag, and FM states are similarly defined by parametersmN, mZ, and mF, respectively) [32].
The evaluated data are plotted in Figs. 4(c) and4(d) using the log-log plots of mS as functions of L for various T at μ=1 and μ=0.8, respectively. One sees immediately the distinctly different behaviors right aboveTcbetween the two casesμ=1 andμ=0.8: the well-defined linear dependence in a relatively wideTrange aboveTcforμ=1, the isotropic model, but no such range exists forμ=0.8, the anisotropic model. We address this difference below.
It is well known that for a KT phase, the order parameter follows the power-law dependence of lattice sizeL, i.e.,mS∼ L−η/2, making the linear dependence of log(mS) on log(L).
This linear dependence is one of the intrinsic characters for the existence of a KT phase, often used to identify the existence of an intermediate KT phase in the present model. Therefore, for the μ=1 case, the linear log(mS)∼log(L) behavior in the T window fromT =Tc∼0.12 toT =0.15 is clearly iden- tified. Furthermore, for a KT phase, the critical exponent η should be aT-dependent continuous variable. Our fitting of the data betweenT =0.12 andT =0.15, i.e., in theTwindow T =T −Tc=0.03, does show the continuous variation of exponentη fromη1=0.112 (T =0.12) to η2 =0.21 (T = 0.15), constituting the two boundary values for this KT phase.
Our results thus confirm the existence of the intermediate KT phase for theμ=1 case on one hand. On the other hand, the two boundary values, η1=0.112 andη2=0.21, agree well with the six-state clock model whereη1=1/9 andη2=1/4 revealed theoretically [33]. Here, without loss of generality, the linear behavior is identified when the calculatedRsquare for the best linear fitting satisfiesR2>0.995.
Then we turn to the anisotropic case withμ=0.8. Indeed, at the critical point Tc=0.212, the linear log(mS)∼log(L) behavior remains true, as shown in Fig.4(d). The best-fitted critical exponent at Tc is estimated to be η=0.12, only slightly deviating from η1 =1/9. Nevertheless, beyond the critical point Tc, the power law is clearly broken and the linear log(mS)∼log(L) behavior is no longer available. Here it is clearly identified from Fig. 4(d) that the data at T = 0.218 begin to deviate from the linear dependence, noting Tc=0.212 and T =T −Tc=0.006. At T =0.224 and T =T −Tc=0.012, the log(mS)∼log(L) curve becomes seriously curved. To this end, one can conclude that the inter- mediate KT phase existing in the isotropic KH model (μ=1) is completely suppressed by the induced Kitaev exchange anisotropy (μ=0.8).
B. Phase transition to the Néel state
The suppression of the intermediate KT phase by the Ki- taev exchange anisotropy, as revealed above for the stripy state as the ground state, is also applicable to the cases of other ground states. Here we show the results on the Néel ground state. The fourfold degenerate Néel state occupies the region of μ >1 and α <1/(2+μ) in the phase diagram.
Figure 5(a) shows the calculated Binder’s cumulant BN as a function of T for variousL atα=0.1 and μ=1.2. The critical point for this continuous phase transition is estimated to beTc=0.304±0.006. This estimatedTcis also confirmed by the log-log plots of order parametermN as a function of L, as shown in Fig.5(b). Specifically, the linear behavior is observed around Tc where the critical line is characterized byη=0.13. However, the data at other temperatures (dashed curves) show remarkable deviation from the linear behavior, demonstrating again the absence of the intermediate phase in the EAKH model with the fourfold degenerate ground state.
In short, it is suggested that the critical behavior for the phase transitions from the high-T disordered phase to the low-T long-range order depends on the degeneracy of the ground state, which can be effectively modulated by the
FIG. 5. Forα=0.1 andμ=1.2, (a) Binder cumulantBN as a function ofTfor differentL, and (b) the log-log plots ofmN as a function ofLat various temperatures.
exchange anisotropy available in real materials [34,35]. Un- like the isotropic KH model or the six-state clock model, the reported intermediate KT phase does not exist in the EAKH model, most likely attributing to the degeneracy reduction of the ground state. As a matter of fact, the phase transition to the twofold degenerate stripy state or Néel state has also been investigated based on the same method, and the intermediate phase is excluded, while the corresponding results are not shown here for brevity.
IV. PHASE TRANSITIONS TO THE ZIGZAG STATE AND FM STATE
Beyond this simple scheme, one needs to investigate a more extended scheme where the magnitudes ofJK andJH
are defined by a more general relation. This extension follows the experimental observations and becomes of high concern.
In fact, a zigzag state was experimentally reported below the critical point in Na2IrO3and such a state seems not available in the simple scheme discussed above. However, it can be also potentially explained by this EAKH model by defining the exchange parameters in a more general way asJK = −Asinϕ andJH =Acosϕ[36]. In this work, we call this definition the extended scheme.
It should be mentioned that the simple scheme can be somehow covered by setting the value ofϕwithJK = −Asinϕ andJH =Acosϕin the 3π/2< ϕ <2πwindow, correspond- ing to the simple scheme withJK =αandJH=1−αin the 0< α <1 window. In this sense, the two schemes address the two slightly different aspects of the same EAKH model, allowing us to investigate all the possible phases, including the zigzag state and FM state here. For example, setting ϕ=11π/16 and A=1, the experimentally reported zigzag ground state can be reproduced. If the anisotropy factor is set asμ=0.8, the state degeneracy is degraded from the sixfold to the fourfold.
The log-log plots of the simulated order parametermZ as a function ofLat variousTare presented in Fig.6(a). A linear behavior is observed at the critical point Tc=0.292 where the critical line is characterized byη=0.125. Similarly, the power-law relationship between mZ and L is not available at other temperatures, demonstrating again the direct phase transition from the high-Tdisordered state to the low-Tzigzag state, without any intermediate KT phase. To our best knowl- edge, no sufficient experimental results have been reported
FIG. 6. The log-log plots of (a)mZ forϕ=11π/16, and (b)mF
forϕ=9π/8 as a function ofLat various temperatures forμ=0.8.
to confirm the intermediate phase in the honeycomb lattice iridates. It was suggested that this absence is partially ascribed to the existence of defects and inhomogeneities. Certainly, these defects and inhomogeneities would at least introduce the deviation of the local Kitaev exchange from the isotropic property, disabling the sixfold degeneracy of the ground state [37,38].
We can also discuss the phase transition to the FM state for integrity of this work, and present the simulated results in Fig. 6(b)where the log-log plots of mF as a function of L at various T for ϕ =9π/8 and μ=0.8 are plotted. In this case, the fourfold degenerate FM ground state is favored and the data show a linear behavior only in the very narrow T range aroundTc=0.438 with critical exponentη=0.13, further confirming the earlier conclusion that no intermediate KT phase is available.
Finally, we tend to compare the present work with the phase transition in related models such as the q-state clock model. As reported earlier, the critical properties of the isotropic KH model are in agreement with the six-state clock model, i.e., the intermediate KT phase has its critical exponent η variable continuously from 1/9 to 1/4 [33,39]. However, earlier work clearly revealed that there is only one transition from the high-T disordered phase to the low-T long-range phase, indicating the absence of any intermediate phase for q<5 [40,41]. To some extent, the finite-Tphase transition to the fourfold degenerate state in the present EAKH model may be related to the transition in the four-state clock model or the four-state Potts model [40,42]. It should be mentioned that the critical behavior of phase transition toward the fourfold degenerate state seems to be much more complicated than a simple extension [34], which deserves further investiga- tion in the future. Moreover, the phase diagrams of extended KH models including extra bond-dependent anisotropic terms have been discussed in earlier works [43–47], but the critical behaviors remain to be further clarified.
Given these discussions, we believe that our work shows the significance of the exchange anisotropy that does result in a degradation of the sixfold degeneracy of the ground state in the isotropic KH model to the fourfold degenerate state and suppresses the intermediate phase. This effect can be an important clue for future study.
V. CONCLUSION
In conclusion, we have studied the phase transitions in the classical KH model with additional Kitaev exchange
anisotropy. The sixfold degeneracy of the ground state in the isotropic model is degraded to a fourfold or a twofold one by the introduced anisotropy, resulting in the suppression of the critical intermediate phase. This property is available in all the finite-T phase transitions to the fourfold degenerate phases including the stripy phase, the Néel phase, the zigzag phase, and the ferromagnetic phase, which may serve as one of the reasons for the missing of intermediate phase in experiments.
ACKNOWLEDGMENTS
The work is supported by the Natural Science Founda- tion of China (Grants No. 51971096 and No. 51721001), the Science and Technology Program of Guangzhou (Grant No.
2019050001), the Natural Science Foundation of Guangdong Province (Grant No. 2019A1515011028), and the Science and Technology Planning Project of Guangzhou in China (Grant No. 201904010019).
[1] A. Kitaev,Ann. Phys.321, 2 (2006).
[2] G. Jackeli and G. Khaliullin, Phys. Rev. Lett. 102, 017205 (2009).
[3] J. Chaloupka, G. Jackeli, and G. Khaliullin,Phys. Rev. Lett.
105, 027204 (2010).
[4] J. G. Rau, E. K.-H. Lee, and H.-Y. Kee,Phys. Rev. Lett.112, 077204 (2014).
[5] A. Banerjee, P. Lampen-Kelley, J. Knolle, C. Balz, A. A. Aczel, B. Winn, Y. H. Liu, D. Pajerowski, J. Q. Yan, C. A. Bridges, A.
T. Savici, B. C. Chakoumakos, M. D. Lumsden, D. A. Tennant, R. Moessner, D. G. Mandrus, and S. E. Nagler, npj Quant.
Mater.3, 8 (2019).
[6] J. Wang, B. Normand, and Z.-X. Liu, Phys. Rev. Lett. 123, 197201 (2019).
[7] J. S. Wen, S. L. Yu, S. Y. Li, W. Q. Yu, and J. X. Li,npj Quant.
Mater.4, 12 (2019).
[8] L. Janssen, E. C. Andrade, and M. Vojta,Phys. Rev. Lett.117, 277202 (2016).
[9] G. W. Chern, Y. Sizyuk, C. Price, and N. B. Perkins,Phys. Rev.
B95, 144427 (2017).
[10] J. Chaloupka, G. Jackeli, and G. Khaliullin,Phys. Rev. Lett.
110, 097204 (2013).
[11] Y. Singh, S. Manni, J. Reuther, T. Berlijn, R. Thomale, W.
Ku, S. Trebst, and P. Gegenwart,Phys. Rev. Lett.108, 127203 (2012).
[12] C. C. Price and N. B. Perkins,Phys. Rev. Lett.109, 187201 (2012).
[13] J. Oitmaa,Phys. Rev. B92, 020405(R) (2015).
[14] X. Liu, T. Berlijn, W.-G. Yin, W. Ku, A. Tsvelik, Y.-J. Kim, H.
Gretarsson, Y. Singh, P. Gegenwart, and J. P. Hill,Phys. Rev. B 83, 220403(R) (2011).
[15] F. Ye, S. X. Chi, H. Cao, B. C. Chakoumakos, J. A. Fernandez- Baca, R. Custelcean, T. F. Qi, O. B. Korneta, and G. Cao,Phys.
Rev. B85, 180403(R) (2012).
[16] J. Reuther, R. Thomale, and S. Rachel, Phys. Rev. B 90, 100405(R) (2014).
[17] S. C. Williams, R. D. Johnson, F. Freund, S. Choi, A. Jesche, I.
Kimchi, S. Manni, A. Bombardi, P. Manuel, P. Gegenwart, and R. Coldea,Phys. Rev. B93, 195158 (2016).
[18] S. Manni, Y. Tokiwa, and P. Gegenwart, Phys. Rev. B 89, 241102(R) (2014).
[19] H. Lei, W.-G. Yin, Z. Zhong, and H. Hosono,Phys. Rev. B89, 020409(R) (2014).
[20] E. C. Andrade and M. Vojta,Phys. Rev. B90, 205112 (2014).
[21] W. P. Cai, Z. R. Yan, R. M. Liu, M. H. Qin, M. Zeng, X. B. Lu, X. S. Gao, and J.-M. Liu,J. Phys.: Condens. Matter29, 405806 (2017).
[22] S. S. Zhang, Z. T. Wang, G. B. Halász, and C. D. Batista,Phys.
Rev. Lett.123, 057201 (2019).
[23] P. S. Anisimov, F. Aust, G. Khaliullin, and M. Daghofer,Phys.
Rev. Lett.122, 177201 (2019).
[24] C. Price and N. B. Perkins,Phys. Rev. B88, 024410 (2013).
[25] M. S. S. Challa and D. P. Landau, Phys. Rev. B 33, 437 (1986).
[26] G. Cao, T. F. Qi, L. Li, J. Terzic, V. S. Cao, S. J. Yuan, M.
Tovar, G. Murthy, and R. K. Kaul,Phys. Rev. B88, 220414(R) (2013).
[27] R. M. Liu, W. Z. Zhuo, J. Chen, M. H. Qin, M. Zeng, X. B. Lu, X. S. Gao, and J.-M. Liu,Phys. Rev. E96, 012103 (2017).
[28] M. Creutz,Phys. Rev. D36, 515 (1987).
[29] J. A. Olive, A. P. Young, and D. Sherrington,Phys. Rev. B34, 6341 (1986).
[30] K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn. 65, 1604 (1996).
[31] K. Binder,Phys. Rev. Lett.47, 693 (1981).
[32] G. Ortiz, E. Cobanera, and Z. Nussinov,Nucl. Phys. B854, 780 (2012).
[33] J. V. José, L. P. Kadanoff, S. Kirkpatrick, and D. R. Nelson, Phys. Rev. B16, 1217 (1977).
[34] A. Kalz and A. Honecker,Phys. Rev. B86, 134410 (2012).
[35] R. M. Liu, W. Z. Zhuo, S. Dong, X. B. Lu, X. S. Gao, M. H.
Qin, and J.-M. Liu,Phys. Rev. E93, 032114 (2016).
[36] D. Gotfryd, J. Rusnaˇcko, K. Wohlfeld, G. Jackeli, J. Chaloupka, and A. M. Ole´s,Phys. Rev. B95, 024426 (2017).
[37] Y. Singh and P. Gegenwart,Phys. Rev. B82, 064412 (2010).
[38] M. Majumder, R. S. Manna, G. Simutis, J. C. Orain, T. Dey, F. Freund, A. Jesche, R. Khasanov, P. K. Biswas, E. Bykova, N. Dubrovinskaia, L. S. Dubrovinsky, R. Yadav, L. Hozoi, S.
Nishimoto, A. A. Tsirlin, and P. Gegenwart,Phys. Rev. Lett.
120, 237202 (2018).
[39] S. K. Baek and P. Minnhagen, Phys. Rev. E 82, 031102 (2010).
[40] H. H. Roomany and H. W. Wyld,Phys. Rev. B23, 1357 (1981).
[41] J. Tobochnik,Phys. Rev. B26, 6201 (1982).
[42] R. J. Baxter,J. Phys. C6, L445 (1973).
[43] K. Shinjo, S. Sota, and T. Tohyama,Phys. Rev. B91, 054401 (2015).
[44] J. Chaloupka and G. Khaliullin, Phys. Rev. B 92, 024413 (2015).
[45] J. Chaloupka and G. Khaliullin, Phys. Rev. B 94, 064435 (2016).
[46] J. Schmidt, D. D. Scherer, and A. M. Black-Schaffer,Phys. Rev.
B97, 014504 (2018).
[47] D. G. Joshi,Phys. Rev. B98, 060405(R) (2018).
