Phase transitions in a frustrated biquadratic Heisenberg model with coupled orbital degrees of freedom for iron-based superconductors
W. Z. Zhuo,1M. H. Qin,1,*S. Dong,2X. G. Li,3and J.-M. Liu1,4,†
1Institute for Advanced Materials and Guangdong Provincial Key Laboratory of Quantum Engineering and Quantum Materials, South China Normal University, Guangzhou 510006, China
2Department of Physics, Southeast University, Nanjing 211189, China
3Hefei National Laboratory for Physical Sciences at Microscale, Department of Physics, University of Science and Technology of China, Hefei 230026, China
4Laboratory of Solid State Microstructures, Nanjing University, Nanjing 210093, China (Received 8 December 2015; revised manuscript received 26 February 2016; published 18 March 2016) In this paper, we study a biquadratic Heisenberg model with coupled orbital degrees of freedom by using a Monte Carlo simulation to investigate the phase transitions in iron-based superconductors. The antiferroquadrupolar state, which may be related to the magnetism of FeSe [R. Yu and Q. Si,Phys. Rev. Lett.115, 116401(2015)], is stabilized by the anisotropic biquadratic interaction induced by a ferro-orbital-ordered state. It is revealed that the orbital and nematic transitions occur at the same temperature for all the cases, supporting the mechanism of the orbital-driven nematicity as revealed in most recent experiments [S. H. Baek, D. V. Efremov, J.
M. Ok, J. S. Kim, J. van den Brink, and B. B¨uchner,Nat. Mater.14,210(2015)]. In addition, it is suggested that the orbital interaction may lead to the separation of the structural and magnetic phase transitions, as observed in many families of iron pnictides.
DOI:10.1103/PhysRevB.93.094424
I. INTRODUCTION
In the past few years, fascinating structural and magnetic phase transitions in iron-based superconductors have drawn extensive attention both experimentally and theoretically [1–5]. Experimentally, the collinear (π, 0) antiferromagnetic (AFM) order is developed at low temperatures (T) in most iron pnictides and usually accompanied by a tetragonal- to-orthorhombic distortion [6–8]. The structural transition temperature TS is either equal to or higher than the AFM transition temperatureTAFM: TSTAFM. To understand this phenomenon, two different mechanisms emphasizing the es- sential role of spin [9,10] and orbital [11,12] fluctuations have been proposed, respectively. Specifically, several theoretical calculations on models for pnictides suggest that the nematic order can be developed before the stabilization of the AFM order and drives the structural phase transition [13–16]. The experimental observations of the spin excitation spectrum can be well explained based on this mechanism [17,18]. However, the magnetism-based origin is strongly challenged by the recent experiments in FeSe [19,20].
As one of the most famous iron chalcogenide superconduc- tors, FeSe shows a very highTCsuperconductivity in its single- layer limit [21–23]. Interestingly, lowingTin bulk FeSe leads to a tetragonal-to-orthorhombic structural transition, while no long-range AFM transition has been observed, suggesting that the structural transition of FeSe may not have the magnetic origin [24]. Most recently, a theoretical study considering only the localized spins revealed a nematic quantum paramagnetic (PM) phase caused by quantum fluctuations and strongly frustrated exchange interactions, which may contribute to the nematicity in FeSe [25]. In addition, a frustrated bilinear-
biquadratic Heisenberg model describing the magnetism of FeSe was studied, and the antiferroquadrupolar (AFQ) and Ising-nematic orders were identified at low T [26]. The structural phase transition in FeSe was suggested to correspond to the Ising-nematic transition in the model. Furthermore, it was proposed that the Goldstone modes of the AFQ order may contribute to the low-energy dipolar magnetic fluctuations observed in the nuclear magnetic resonance measurements [24,27]. In addition, an unusual magnetic frustration was pro- posed in a theoretical work and was suggested to suppress mag- netic order and trigger ferro-orbital order in the nematic phase, consistent with the pressure dependence ofTCin FeSe [28].
On the other hand, recent experiments clearly demonstrated that the orbital transition may play a crucial role in the onset of nematic transition in bulk FeSe, strongly supporting the mechanism of orbital-driven nematicity [27]. In addition, it was observed that the shear-modulus softening above TS in FeSe and underdoped BaFe2As2 is nearly identical, suggesting a common origin for the structural transition in these materials [24]. As a matter of fact, the important role of orbital fluctuations has been addressed in several models for pnictides [29–32]. For example, the nature of phase transitions in pnictides has been investigated based on a spin-orbital model with coupled spin and orbital degrees of freedom using Monte Carlo (MC) simulations [33]. However, the biquadratic interaction, which is believed to be very important in iron superconductors [34], is neglected in this model, while a coupling to orbital degrees of freedom is not considered in the localized spin models describing the magnetism of FeSe discussed above. Thus, more model calculations studying the role of orbital fluctuations in iron pnictides and chalcogenides are still urgently needed to elucidate the physical origins for the nematicity and some other experimental observations.
In this paper, we study a classical biquadratic Heisenberg model with coupled orbital degrees of freedom. Several experimental observations in iron-based superconductors can
be qualitatively explained in our MC simulations. In detail, both the AFM and AFQ orders are predicted in the phase diagram of the model. The orbital phase transition temperature TO is always the same as or higher than the AFM (AFQ) transition temperature TAFM(TAFQ), and a nematic ordering can also be developed atT∗≈TO, supporting the mechanism of orbital-driven nematicity. The remainder of this paper is organized as follows: in Sec.IIthe model and the simulation method will be described. Section III is attributed to the simulation results and discussion, and the conclusion is presented in Sec.IV.
II. MODEL AND METHOD
We study a classical biquadratic Heisenberg spin (S=1) model with coupled spin and orbital degrees of freedom on a two-dimensional square lattice. It was reported that a ferro-orbital order in which the dxz/dyz orbital dominates can be developed in iron pnictides and FeSe as Tdecreases [18,27,35]. For simplicity, the orbital parameterOis set to be Ising variables with values 1 (occupation ofdxzorbital) or−1 (occupation ofdyzorbital). The model Hamiltonian is stated as
H =
i
J1+xi ·JOSSi· Si+x
+
i
J1+yi ·JOSSi· Si+y+J2
[ij]
Si· Sj
−
i
K1+xi ·KOS
(Si· Si+x)2
−
i
K1+yi ·KOS
(Si· Si+y)2−JO
ij
Oi·Oj,
(1) with
xi =
1 Oi=Oi+x =1
0 Oi= −1 or Oi+x= −1 yi =
1 Oi=Oi+y = −1 0 Oi=1 or Oi+y =1.
(2)
The first two terms are the nearest neighbor (NN) exchange interactions along the x andy directions. Following earlier work, the spin and orbital degrees of freedom are considered to be coupled by a Kugel-Khomskii-like mechanism [33,36], with the coupling magnitudeJOS. The third term denotes the next NN exchange interaction with couplingJ2. The fourth and fifth terms are the NN biquadratic interactions along thexand y directions, respectively, which are expected in any model calculations for iron superconductors [34]. In fact, it has been proved that the biquadratic interaction is also affected by the Fe-Fe bond length [37], which can be modulated by orbital orderings. Thus, the biquadratic interactions should be also dependent on orbital orders (with couplingKOS). At last, the orbital interaction with couplingJOis considered, which has been involved in several earlier models [12].
Generally, the introduce of the coupled orbital degrees of freedom and biquadratic interaction in the Hamiltonian can be interpreted to change the linear-response exchange constant,
which can be defined by [13]
Jij ≡ −(Si· Sj)−1∂2H
∂θ2 , (3)
whereθ is the angle between spins at sites i andj. For the (π, 0) AFM ground state normally accompanied by a dxz
orbital polarization, as experimentally reported in several iron pnictides [35], we haveJijx =J1+xi(JOS+2KOS)+2K1. Furthermore, it is suggested that the magnitude of Jijx may be increased when the sitesiandj are occupied by thedxz
orbital [12]. Thus, we choose KOS>0 and JOS>0 in the whole paper and setKOS=1 as the energy unit, for simplicity.
In addition, we fixJ2=0.7J1 andJOS=0.5J1, which are comparable with experimental reports [38]. Furthermore, we take the other coupling parameters as variables and study the phase diagram of the model using the standard Metropolis algorithm and temperature exchange method [39,40]. Unless stated elsewhere, the simulation is performed on a 36×36 lattice with periodic boundary conditions. To study the finiteT phase transitions to long-range magnetic orders, we introduce a spin-space anisotropy:
Si· Sj =SizSjz+λ
SixSjx+SiySyj
. (4)
To explore the phases in the system, the dipolar and quadrupolar magnetic structure factors are calculated by [26]
d(q)= 1 N2
ij
Si· Sjexp[iq·(ri− ri)], (5)
and
q(q) = 1 N2
ij
Qi· Qjexp[iq·(ri− rj)], (6) respectively. Here,Nis the number of sites, and. . .is the ensemble average. For classical Heisenberg spins,
Qi· Qj =2(Si· Sj)2−2
3Si2· S2j. (7) Then, the AFM and AFQ order parameters (mandmq) and the corresponding nematic order parameters (NAFMandNAFQ) are calculated by [13,26]
m2=d(π,0)+d(0,π)
mq2=q(π,0)+q(0,π), (8) NAFM=
1 2N
N
i
(Si· Si+x− Si· Si+y)
(9) NAFQ=
1 2N
N
i
(Qi· Qi+x− Qi· Qi+y)
.
The transition temperatures (TO, T∗, TAFQ, andTAFM) can be estimated from the positions of the peaks in the correspond- ing susceptibilityχO, χN, χAFQ, andχAFM, respectively.
III. SIMULATION RESULTS AND DISCUSSION A. The AFM and AFQ orders
First, we study the spin orders developed at low T for various J1 and K1. Figure 1 shows the simulated phase
FIG. 1. Calculated phase diagram in the (J1,K1) parameter space at T =0.01 for λ=0.95 and JO=0. The red line shows the boundary atT =0.
diagram in the (J1, K1) plane at T =0.01 for λ=0.95 and JO=0. The phase diagram is occupied separately by the AFM order and AFQ order. It is indicated that the strong competition between the exchange interaction and the biquadratic interaction can stabilize the AFQ order. Here, note that in the bilinear-biquadratic spin model, an isotropic NN biquadratic interaction alone cannot stabilize the AFQ order [26]. To understand the physics underlying our simulated results on the AFQ order, a qualitative discussion in the energy landscape is helpful.
Figures 2(a) and 2(b) show the distribution of dipolar and quadrupolar magnetic structure factors at T =0.01 for J1=0.1, K1= −0.5,JO=0 andλ=0.6, as an example. A dominant (π, 0) AFQ correlation is clearly shown, with the three spin components of the AFQ order plotted in Figs.2(c) and2(d), respectively. For an AFQ order, the NN spins along
(a) (b)
(c) (d)
FIG. 2. Distribution of the dipolar (a) and quadrupolar (b) magnetic structure factors at T =0.01 for J1=0.1, K1= −0.5, JO=0, and λ=0.6. Thezcomponent (c) and spin configuration on thexyplane (d) of the AFQ order.
thex direction are perpendicular to each other, while those along they direction are antiparallel with each other. In this case, the dyz ferro-orbital-ordered state is developed at TO (not shown here), resulting in the anisotropic NN biquadratic interaction (K1along thex direction andK1+1 along they direction). Thus, the NN spins along thex direction tend to be perpendicular to each other to save the negative biquadratic interactionK1 at the cost of the AFM exchange interaction J1, while those along theydirection are antiparallel with each other to satisfy the AFM exchange and positive biquadratic interactions, leading to the AFQ state. With increasing
|K1|(K1<0), the AFQ state can be further stabilized, resulting in the enlargement of theJ1 region with the AFQ order. On the other hand, the collinear AFM order can be stabilized by a positive biquadratic interaction (K1>0), and a positiveK1 never stabilizes the AFQ state, as confirmed in our simulations.
B. The AFQ phase transitions
Given the AFQ order in the phase diagram, we subsequently draw our attention onto the expected nematic order and the finiteTphase transitions. The order parameters for these phase transitions include orbital parameterOz= Oj/N, NAFQ, andmq. Figure3(a)shows these parameters as a function ofT forJ1=0.1, K1= −0.5,JO=0 andλ=0.6, as an example.
At highT, the system is in the PM state, and all these order parameters are rather small. WhenTfalls down to the critical points, these parameters increase from the baseline, fingering the development of orbital order, nematic order, and AFQ order, respectively. These transition points (TO,T∗, andTAFQ) are estimated from the positions of peaks in the calculated susceptibility on these parameters, as shown in Fig.3(b). It is clearly shown that the nematic order is developed at the same temperature as the orbital order does, givingT∗∼TO, consistent with experimental observations [27]. Furthermore, the AFQ order is stabilized at a slightly lowerT, givingTAFQ<
TO. Similar behaviors are confirmed for different lattice sizes, and the conclusion regarding the nematic state is reliable.
The effects of spin anisotropy parameter λ, orbital in- teraction JO (i.e., ferro-orbital coupling), and biquadratic interactionK1 on the phase transitions are also investigated, respectively, and the simulated phase diagrams are plotted in Fig.4. Figure4(a)shows the phase diagram on the (λ,T) space forJ1=0.1, K1= −0.5, andJO=0. A reducedλ(increased spin anisotropy) downshifts all the three transitions points.
FIG. 3. The calculated order parameters OZ, NAFQ, and mq (a) and their susceptibilities (b) as a function ofTforJ1=0.1, K1=
−0.5,JO=0 andλ=0.6.
(a) (b) (c)
FIG. 4. The simulated phase diagram (a) in the (λ, T) plane forJ1=0.1, K1= −0.5, andJO=0, (b) in the (JO,T) plane forJ1= 0.1, K1= −0.5, andλ=0.4, and (c) in the (K1,T) plane forJ1=0.1,J o=0, andλ=0.6.
This is understandable since the spin anisotropy favors the AFM order rather than the AFQ order. At the same time, the coupling between the spin and orbital degrees of freedom, in turn, destabilizes the orbital/nematic order. The nematic order in advantage of the AFQ order can be destroyed eventually by the spin anisotropy, leading to narrowed|TAFQ−T∗|with increasing (1−λ), i.e., the nematic order region is suppressed.
On the other hand, the orbital order can be further stabilized by a ferro-orbital coupling JO, as shown in Fig. 4(b) for the phase diagram on the (JO,T) parameter space forJ1= 0.1, K1= −0.5, andλ=0.4. Interestingly, the nematic order region is substantially enlarged asJOincreases, andT∗ ∼TO applies for all these cases, again demonstrating the important role of the orbital order in driving the nematicity in FeSe.
Furthermore, Fig.4(c)shows the phase diagram on the (K1, T) parameter space for J1=0.1, JO=0, and λ=0.6. The three transition temperatures decrease as the magnitude of K1increases. Note that the energy gain from the biquadratic interaction due to the phase transition from the PM state to the AFQ state is extensively decreased with the increasing|K1|, leading to the destabilization of the AFQ state. In turn, the nematic state and the orbital order are also destabilized, as shown in our simulations.
C. The AFM phase transitions
For integrity, we also study the collinear AFM phase transition for the case of dominantJ1. The order parameters and their susceptibilities at T =0.01 for J1 =1, K1=0.5, JO=0, and λ=0.95 are given in Figs. 5(a) and 5(b), respectively, clearly showing that the three phase transitions occur at the same temperature, i.e.,TAFM∼T∗ ∼TO. In fact, this phenomenon was reported in some iron pnictides, such as undoped 122 compounds CaFe2As2 and SrFe2As2 [41,42].
Most recently, the ferro-orbital order has been proved to generate the exchange anisotropy and stabilize the collinear AFM order in iron-based superconductors by spin-wave analysis of the spin-fermion model [43,44], well consistent with our simulations [45].
Similarly, the phase transitions can be also modulated by spin anisotropy, orbital interaction, and biquadratic interaction, and the corresponding results are given in Fig.6. As discussed earlier, a collinear AFM order is favored by the spin anisotropy, and TAFM will increase with decreasing λ, as confirmed in
Fig.6(a), which gives the phase diagram on the (λ,T) space for J1=1.0, K1 =0.5, andJO=0.2. Furthermore, the nematic order region is also suppressed with increasing (1−λ), and all the three transitions occur at the same temperature forλ <
0.8. Figure6(b)shows the phase diagram on the (JO,T) space forJ1 =1.0, K1=0.5, andλ=0.95. All the three transition temperatures increase with increasing JO, and a separation between TAFQ and T∗ can be observed for JO>0.1. The phase diagram on the (K1,T) space for J1 =1.0, λ=0.95, andJO=0.2 is shown in Fig.6(c). All the three transitions shift toward high T, as K1 increases due to the increasing positive biquadratic interaction.
The nature of these transitions is also investigated by the analysis of the Binder ratios. Figure7shows the Binder ratios for orbital gO, nematic gN, and spin gS(gO= Oz4/Oz22, as an example) vs T for J1=1.0, K1=0.5, and λ=0.95.
The clear divergences in the three curves forJO=0 indicate that these transitions are of the first order. Interestingly, the divergences ofgOandgN vanish whenJOis increased above 0.1, as shown in Fig.7(b), suggesting a trend to the second order orbital and nematic transitions. In addition, there is still a weak divergence in the spin Binder ratiogS, which is indicative of a weak/pseudo-first-order transition [46]. The nature of the AFM transition should be further checked by simulations on larger lattice sizes [47]. However, it is observed that the nematic transition always occurs at the same temperature as that of the ferro-orbital ordering, suggesting that the same mechanism for nematicity may also work in iron pnictides, in some extent.
FIG. 5. The calculated order parametersOZ,NAFM, andm(a) and their susceptibilities (b) as a function ofT forJ1=1.0, K1=0.5, JO=0, andλ=0.95.
(a) (b) (c)
FIG. 6. The simulated phase diagram (a) in the (λ, T) plane for J1=1.0, K1=0.5, and JO=0.2 and (b) in the (JO, T) plane for J1=1.0, K1=0.5, andλ=0.95, and (c) in the (K1,T) plane forJ1=1.0, λ=0.95, andJO=0.2.
D. Relevance to iron superconductors
Most recently, it was experimentally reported that spin- lattice relaxation rates are not affected at the nematic transition point, strongly suggesting the orbital origin for the nematicity in FeSe [27]. This viewpoint is supported in our simulations when the orbital order is suggested to be the primary driving mechanism for the finite temperature transitions. Furthermore, the AFQ order has been observed in the bilinear-biquadratic model, and its Goldstone modes are suggested to contribute to dipolar magnetic fluctuations (in the absence of any AFM order) revealed in experiments [26]. However, a rather strong biquadratic interaction between the next NNs, which has not been confirmed in first-principles calculations [28], is proved to be necessary to stabilize the AFQ order. In this paper, it is suggested that the AFQ order can be stabilized by an anisotropic NN biquadratic interaction induced by the development of orbital order. Furthermore, the AFQ state can be replace by the AFM state by tuning the exchange and biquadratic interactions, as shown in the phase diagram in Fig.1, which may contribute to the experimental observations of the appearance of magnetic order under external pressure in FeSe [48,49]. In addition, the nematic order accompanied by the orbital order may be developed in advance of the AFQ order when the orbital-orbital interaction is increased, further demonstrating the important role of the orbital fluctuations in the phase transitions in FeSe.
On the other hand, in iron pnictides CaFe2As2 and SrFe2As2, both the structural and magnetic phase transitions
(a) (b)
FIG. 7. CalculatedgO, gN, andgS as a function ofT forJ1= 1.0, K1=0.5, andλ=0.95 at (a)JO=0 and (b)JO=0.25.
occurs at a same temperature, which is reproduced in our simulations for smallJO<0.1. As a matter of fact, the phase transition in iron pnicitides have been studied based on a similar spin-orbital model, and several behaviors are captured by the simulated phase diagram. For example, even for very weak spin anisotropy (λ∼0.95), both the orbital and AFM transition points coincide with each other, consistent with our simulations forJO=0. Interestingly, our paper suggests that a ferro-orbital interaction may separate the two transitions and modulate the nature of the transitions. This behavior may contribute to some of the experimental observations in many families of iron pnictides. For example, bothTS andTAFMin BaFe2As2 are increased, and TS becomes well above TAFM when a compressive stress is applied [50,51]. To some extent, this phenomenon may be related to the increase of the orbital- orbital interaction, which deserves to be checked further.
IV. CONCLUSION
In conclusion, we have studied phase transitions in a biquadratic Heisenberg model with coupled orbital degrees of freedom. It is suggested that the development of the ferro-orbital-ordered state induces the anisotropic biquadratic interaction between NNs, and, in turn, stabilizes the AFQ state, which may be related to the magnetism of FeSe. For all the cases, the orbital and nematic transitions occur at a same temperature, supporting the mechanism of the orbital- driven nematicity as revealed in experiments. Furthermore, the orbital-orbital interaction may contribute to the separation of the structural and magnetic phase transitions as observed in many families of iron pnictides.
ACKNOWLEDGMENTS
The paper is supported by the Natural Science Foundation of China (51332007, 51322206, 11274094), and the National Key Projects for Basic Research of China (2015CB921202 and 2015CB654602), and the International Science & Tech- nology Cooperation Platform Program of Guangzhou (Grant No. 2014J4500016), and the Science and Technology Plan- ning Project of Guangdong Province, China (Grant No.
2015B090927006).
[1] D. C. Johnston,Adv. Phys.59,803(2010).
[2] G. R. Stewart,Rev. Mod. Phys.83,1589(2011).
[3] P. C. Dai, J. P. Hu, and E. Dagotto,Nat. Phys.8,709(2012).
[4] E. Dagotto,Rev. Mod. Phys.85,849(2013).
[5] P. C. Dai,Rev. Mod. Phys.87,855(2015).
[6] C. de la Cruz, Q. Huang, J. W. Lynn, J. Y. Li, W. Ratcliff, J. L.
Zarestky, H. A. Mook, G. F. Chen, J. L. Luo, N. L. Wang, and P. C. Dai,Nature (London)453,899(2008).
[7] S. Nandi, M. G. Kim, A. Kreyssig, R. M. Fernandes, D. K. Pratt, A. Thaler, N. Ni, S. L. Bud’ko, P. C. Canfield, J. Schmalian, R.
J. McQueeney, and A. I. Goldman,Phys. Rev. Lett.104,057006 (2010).
[8] S. Kasahara, H. J. Shi, K. Hashimoto, S. Tonegawa, Y.
Mizukami, T. Shibauchi, K. Sugimoto, T. Fukuda, T. Terashima, A. H. Nevidomskyy, and Y. Matsuda,Nature (London)486,382 (2012).
[9] C. Fang, H. Yao, W. F. Tsai, J. P. Hu, and S. A. Kivelson,Phys.
Rev. B77,224509(2008).
[10] R. M. Fernandes, L. H. VanBebber, S. Bhattacharya, P. Chandra, V. Keppens, D. Mandrus, M. A. McGuire, B. C. Sales, A. S.
Sefat, and J. Schmalian,Phys. Rev. Lett.105,157003(2010).
[11] C. C. Lee, W. G. Yin, and W. Ku,Phys. Rev. Lett.103,267001 (2009).
[12] W. C. Lv, J. S. Wu, and P. Phillips,Phys. Rev. B80,224506 (2009).
[13] A. L. Wysocki, K. D. Belashchenko, and V. P. Antropov,Nat.
Phys.7,485(2011).
[14] J. P. Hu, C. Setty, and S. Kivelson,Phys. Rev. B85,100507(R) (2012).
[15] S. H. Liang, A. Moreo, and E. Dagotto,Phys. Rev. Lett.111, 047004(2013).
[16] M. H. Qin, S. Dong, H. B. Zhao, Y. Wang, J. M. Liu, and Z. F.
Ren,New J. Phys.16,053027(2014).
[17] J. Zhao, D. T. Adroja, D. X. Yao, R. Bewley, S. Li, X. F. Wang, G. Wu, X. H. Chen, J. P. Hu, and P. C. Dai,Nat. Phys.5,555 (2009).
[18] M. Yi, D. Lu, J.-H. Chu, J. G. Analytis, A. P. Sorini, A. F.
Kemper, B. Moritz, S.-K. Mo, R. G. Moore, M. Hashimoto, W.-S. Lee, Z. Hussain, T. P. Devereaux, I. R. Fisher, and Z.-X.
Shen,Proc. Natl. Acad. Sci. USA108,6878(2011).
[19] T. M. McQueen, A. J. Williams, P. W. Stephens, J. Tao, Y. Zhu, V. Ksenofontov, F. Casper, C. Felser, and R. J. Cava,Phys. Rev.
Lett.103,057002(2009).
[20] S. Medvedev, T. M. McQueen, I. A. Troyan, T. Palasyuk, M. I.
Eremets, R. J. Cava, S. Naghavi, F. Casper, V. Ksenofontov, G.
Wortmann, and C. Felser,Nat. Mater8,630(2009).
[21] Q. Y. Wang, Z. Li, W. H. Zhang, Z. C. Zhang, J. S. Zhang, W.
Li, H. Ding, Y. B. Ou, P. Deng, K. Chang, J. Wen, C. L. Song, K. He, J. F. Jia, S. H. Ji, Y. Y. Wang, L. L. Wang, X. Chen, X.
C. Ma, and Q. K. Xue,Chin. Phys. Lett.29,037402(2012).
[22] J. F. Ge, Z. L. Liu, C. Liu, C. L. Gao, D. Qian, Q. K. Xue, Y.
Liu, and J. F. Jia,Nat. Mater.14,285(2015).
[23] D. H. Lee,Chin. Phys. B24,117405(2015).
[24] A. E. B¨ohmer, T. Arai, F. Hardy, T. Hattori, T. Iye, T. Wolf, H.
v. L¨ohneysen, K. Ishida, and C. Meingast,Phys. Rev. Lett.114, 027001(2015).
[25] F. Wang, S. A. Kivelson, and D. H. Li,Nat. Phys.11,959(2015).
[26] R. Yu and Q. Si,Phys. Rev. Lett.115,116401(2015).
[27] S. H. Baek, D. V. Efremov, J. M. Ok, J. S. Kim, J. van den Brink, and B. B¨uchner,Nat. Mater.14,210(2015).
[28] J. K. Glasbrenner, I. I. Mazin, H. O. Jeschke, P. J. Hirschfeld, R. M. Fernandes, and R. Valenti,Nat. Phys.11,953(2015).
[29] C.-C. Chen, B. Moritz, J. van den Brink, T. P. Devereaux, and R. R. P. Singh,Phys. Rev B80,180418(R) (2009).
[30] C.-C. Chen, J. Maciejko, A. P. Sorini, B. Moritz, R. R. P. Singh, and T. P. Devereaux,Phys. Rev B82,100504(R) (2010).
[31] H. Kontani, T. Saito, and S. Onari,Phys. Rev. B84, 024528 (2011).
[32] W. C. Lee, W. Lv, J. M. Tranquada, and P. W. Phillips,Phys.
Rev. B86,094516(2012).
[33] R. Applegate, R. R. P. Singh, C. C. Chen, and T. P. Devereaux, Phys. Rev. B85,054411(2012).
[34] J. K. Glasbrenner, J. P. Velev, and I. I. Mazin,Phys. Rev. B89, 064509(2014).
[35] T. Shimojima, K. Ishizaka, Y. Ishida, N. Katayama, K. Ohgushi, T. Kiss, M. Okawa, T. Togashi, X. Y. Wang, C. T. Chen, S.
Watanabe, R. Kadota, T. Oguchi, A. Chainani, and S. Shin, Phys. Rev. Lett.104,057002(2010).
[36] I. K. Kliment and D. I. Khomskii,Sov. Phys. Usp. 25, 231 (1982).
[37] A. L. Wysocki, K. D. Belashchenko, L. Ke, M. van Schilfgaarde, and V. P. Antropov,J. Phys.: Conf. Ser.449,012024(2013).
[38] R. A. Ewings, T. G. Perring, J. Gillett, S. D. Das, S. E. Sebastian, A. E. Taylor, T. Guidi, and A. T. Boothroyd,Phys. Rev. B83, 214519(2011).
[39] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H.
Teller, and E. Teller,J. Chem. Phys.21,1087(1953).
[40] K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn. 65, 1604 (1996).
[41] N. Ni, S. Nandi, A. Kreyssig, A. I. Goldman, E. D. Mun, S. L.
Bud’ko, and P. C. Canfield,Phys. Rev. B78,014523(2008).
[42] J. Q. Yan, A. Kreyssig, S. Nandi, N. Ni, S. L. Bud’ko, A.
Kracher, R. J. McQueeney, R. W. McCallum, T. A. Lograsso, A. I. Goldman, and P. C. Canfield,Phys. Rev. B 78, 024516 (2008).
[43] W.-G. Yin, C.-C. Lee, and W. Ku,Phys. Rev. Lett.105,107004 (2010).
[44] W. C. Lv, F. Kruger, and P. Phillips,Phys. Rev. B82,045125 (2010).
[45] Y. T. Tam, D. X. Yao, and W. Ku,Phys. Rev. Lett.115,117001 (2015).
[46] S. B. Jin, A. Sen, and A. W. Sandvik,Phys. Rev. Lett. 108, 045702(2012).
[47] A. Kalz and A. Honecker,Phys. Rev. B86,134410(2012).
[48] T. Imai, K. Ahilan, F. L. Ning, T. M. McQueen, and R. J. Cava, Phys. Rev. Lett.102,177005(2009).
[49] M. Bendele, A. Amato, K. Conder, M. Elender, H. Keller, H.
H. Klauss, H. Luetkens, E. Pomjakushina, A. Raselli, and R.
Khasanov,Phys. Rev. Lett.104,087003(2010).
[50] C. Dhital, Z. Yamani, W. Tian, J. Zeretsky, A. S. Sefat, Z. Q.
Wang, R. J. Birgeneau, and S. D. Wilson,Phys. Rev. Lett.108, 087001(2012).
[51] M. H. Qin, S. Dong, J.-M. Liu, and Z. F. Ren,New J. Phys.17, 013011(2015).
