Composite spin and quadrupole wave in the ordered phase of Tb 2+x Ti 2−x O 7+y
aa r X i v : . [ c ond - m a t . s t r- e l ] J a n Composite spin and quadrupole wave in the ordered phase of Tb x Ti − x O y H. Kadowaki, H. Takatsu, T. Taniguchi, B. F˚ak, and J. Ollivier Department of Physics, Tokyo Metropolitan University, Hachioji-shi, Tokyo 192-0397, Japan Institute Laue Langevin, BP156, F-38042 Grenoble, France (Dated: January 26, 2016)The hidden ordered state of the frustrated pyrochlore oxide Tb x Ti − x O y is possibly oneof the two electric multipolar, or quadrupolar, states of the effective pseudospin-1/2 Hamiltonianderived from crystal-field ground state doublets of non-Kramers Tb ions. These long-range ordersare antiparallel or parallel alignments of transverse pseudospin components representing electricquadrupole moments, which cannot be observed as magnetic Bragg reflections by neutron scattering.However pseudospin waves of these states are composite waves of the magnetic-dipole and electric-quadrupole moments, and can be partly observed by inelastic magnetic neutron scattering. Wecalculate these spin-quadrupole waves using linear spin-wave theory and discuss previously observedlow-energy magnetic excitation spectra of a polycrystalline sample with x = 0 .
005 ( T c = 0 . I. INTRODUCTION
Geometrically frustrated magnets have been activelystudied in recent years . In particular, pyrochloremagnets showing spin ice behavior have interestingfeatures such as finite zero-point entropy and emergentmagnetic monopole excitations . A quantum spin-liquidstate is theoretically predicted for certain spin-ice likesystems , where transverse spin interactions transformthe classical spin ice into quantum spin liquid. Thisquantum spin ice (QSI), or U(1) quantum spin liquid, ischaracterized by an emergent U(1) gauge field fluctuat-ing down to T = 0 and by excitations of gapped bosonicspinons and gapless photons . By changing the in-teractions of the QSI in some ways the system under-goes a quantum phase transition to long range ordered(LRO) states of transverse spin or pseudospin , beinginterpreted as Higgs phases . Experimental investi-gations of the U(1) quantum spin liquid and neighboringLRO states have been challenged by several groups .However it is difficult to characterize the quantum spinliquid states, which preclude standard techniques of ob-serving magnetic Bragg reflections and magnons.Among magnetic pyrochlore oxides , R Ti O (R= Dy, Ho) are the well-known classical Ising spin-iceexamples . A similar system Tb Ti O (TTO) has at-tracted much attention, because magnetic moments re-main dynamic with short range correlations down to50 mK . Since TTO has been thought to be close tothe classical spin ice, the low-temperature dynamical be-havior of TTO could be attributed to QSI . Inspiredby this intriguing idea, many experimental studies ofTTO have been performed to date (and referencesin Refs. 9 and 21). However the interpretation of exper-imental data has been a conundrum , partly owing tostrong sample dependence . Among these studies,our investigation of polycrystalline Tb x Ti − x O y showed that a very small change of x induces a quan-tum phase transition between a spin-liquid state ( x < − . x c ) and a LRO state with a hidden order pa-rameter ( x c < x ). It is important to clarify the origin of this order parameter, which becomes dynamical in thespin-liquid state ( x < x c ).In this and companion work, we try to reformu-late the problem of TTO and to reinterpret its puzzlingexperimental data based on the theoretically predicted electronic superexchange interactions. A novel ingredientof these interactions is the Onoda-type coupling be-tween neighboring electric quadrupole moments of non-Kramers Tb ions. The theory proposes an effec-tive pseudospin-1/2 Hamiltonian described by the Paulimatrices representing both magnetic-dipole and electric-quadrupole moments. Depending on the parameters ofthe Hamiltonian there are two electric quadrupole order-ing phases, which are candidates for the hidden orderof TTO. These electric quadrupolar orders do not bringabout observable magnetic Bragg peaks. However, theseorders can be detected by their elementary excitations(inelastic magnetic scattering), and by proper interpre-tation using a linear spin-wave theory.In this paper, starting from the crystal-field (CF)ground state doublet of TTO, we account for its single-site electric quadrupole moments, their LRO, and pseu-dospin wave excitations in the electric quadrupole LRO.A standard linear spin-wave theory predicts that thepseudospin wave in the electric quadrupole LRO is, inreality, a composite wave of magnetic-dipole and electric-quadrupole moments. We discuss this possibility forTb x Ti − x O y using previously observed low-energymagnetic excitation spectra of a polycrystalline samplewith x = 0 .
005 ( T c = 0 . II. CRYSTAL FIELD AND ELECTRICMULTIPOLE MOMENT
The CF states and inelastic neutron excitation spectraof TTO have been investigated by many authors ;readers are referred to Ref. 30 for details. In a low en-ergy range, there are four CF states: ground doubletstates and first-excited doublet states at E ∼
16 K. Sincethe interesting temperature range is below 1 K, we ne-glect the first-excited doublet states and consider onlythe ground state doublet, for simplicity.Among studies of CF, we adopt the CF parameters ofRef. 27 (or Ref. 28). The CF ground state doublet ofTTO can be written by | ± i D = A | ± i ∓ B | ± i + C | ∓ i ± D | ∓ i , (1)where | m i stands for the | J = 6 , m i state within a JLS -multiplet . The coefficients of Eq. (1) are A = 0 . B = 0 . C = 0 . D = 0 . of the crystallographic four sites are x = √ (1 , , ¯2) , y = √ (¯1 , , , z = √ (1 , ,
1) (2)for sites at t n + d with d = (0 , , x = √ (1 , ¯1 , , y = √ (¯1 , ¯1 , , z = √ (1 , ¯1 , ¯1) (3)for sites at t n + d with d = (0 , , x = √ (¯1 , , , y = √ (1 , , , z = √ (¯1 , , ¯1) (4)for sites at t n + d with d = (1 , , x = √ (¯1 , ¯1 , ¯2) , y = √ (1 , ¯1 , , z = √ (¯1 , ¯1 ,
1) (5)for sites at t n + d with d = (1 , , t n is anFCC translation vector.In the CF ground state doublet of Eq. (1),the magnetic-dipole and electric-multipole momentoperators are represented by 2 × σ x , σ y , σ z and the unit matrix. The magneticdipole moment operators within | ± i D are J x = J y = 0 ,J z = (4 A + B − C − D ) σ z = 3 . σ z , (6)which implies that Tb magnetic dipole moments be-have as Ising-like spins.As pointed out in Ref. 25, for non-Kramers ions inthe pyrochlore structure including Tb in TTO the CFground doublet states have additionally electric multipolemoments. These electric multipole moment operators arerepresented by σ x , σ y , and the unit matrix. Using theexplicit form of Eq. (1), the electric quadrupole momentoperators within | ± i D are expressed by [3 J z − J ( J + 1)] = 3 A − B − C + D = 3 . , √ [ J x − J y ] = (cid:16) − √ B + 9 √ AC (cid:17) σ x = 3 . σ x , √ [ J x J y + J y J x ] = (cid:16) − √ B + 9 √ AC (cid:17) σ y = 3 . σ y , √ [ J z J x + J x J z ] = (cid:18) − √ BC − q AD (cid:19) σ x = − . σ x , √ [ J y J z + J z J y ] = (cid:18) √ BC + 9 q AD (cid:19) σ y = 8 . σ y . (7) Similarly we can show that the electric 16-pole and64-pole moment operators , expressed by the Racahoperators ˜ O p,q ( J ) with p = 4 and 6, respectively (orStevens’s operators), are proportional to σ x ± iσ y or theunit matrix within | ± i D . Therefore within the CFground state doublet, pseudospin operators σ x and σ y represent the electric multipole moments. A single-siteCF ground state expressed by | ψ i = ( | i D , | − i D ) χ, (8)where χ is the pseudospin wave-function, has the largestexpectation of the magnetic dipole moment |h ψ | σ z | ψ i| =1 (and h ψ | σ x | ψ i = h ψ | σ y | ψ i = 0) for χ = (cid:18) (cid:19) or χ = (cid:18) (cid:19) . The other states expressed by χ = (cid:18) cos θ e − iφ/ sin θ e iφ/ (cid:19) (9)in which θ is in the range 0 < θ < π have finite expecta-tion values of the electric quadrupole moment operators; h ψ | σ x | ψ i 6 = 0 and/or h ψ | σ y | ψ i 6 = 0. These states haveslightly deformed f -electron charge densities from thatof the magnetic states with θ = 0 or θ = π . More specifi-cally, the approximate f -electron charge density of thestate | ψ i is given by h ψ | ρ ( r ) | ψ i ≃ ( − e )[ R f ( r )] h ψ | ρ e (ˆ r ) | ψ i π . (10)The angular dependence ρ e (ˆ r ) of this equation is ρ e (ˆ r ) = n + X p =2 , , q [4 π (2 p + 1)] / α p Y p,q (ˆ r ) ∗ ˜ O p,q ( J ) , (11)where ( α , α , α ) = ( α, β, γ ) are the Stevens factors , n = 8 is the number of f -electrons, and Y p,q (ˆ r ) are thespherical harmonics. By evaluating h ψ | ρ e (ˆ r ) | ψ i usingseveral spinors of Eq. (9), one can show that the defor-mation of the f -electron charge density is mainly deter-mined by the electric quadrupole moments. The electric16-pole and 64-pole moments have non-negligible contri-butions to the deformation similarly to the analyses ofthe CF states . In these meanings, the CF ground(psedospin-1/2) states | ψ i can be referred to as compos-ite spin and quadrupole states. III. EFFECTIVE PSEUDOSPIN-1/2HAMILTONIAN
The generic form of the effective pseudospin-1/2Hamiltonian for non-Kramers CF ground state doubletsof 4 f magnetic ions in the pyrochlore structure was de-rived in Ref. 25 by calculating the nearest-neighbor (NN)superexchange interaction. This Hamiltonian consists oftwo parts. The first part is the NN magnetic interaction H m,NN = J nn X h r , r ′ i σ z r σ z r ′ , (12)which represents the NN classical spin-ice model for J nn >
0. The second part is the NN quadrupolar in-teraction H q = J nn X h r , r ′ i [2 δ ( σ + r σ − r ′ + σ − r σ + r ′ )+ 2 q ( e iφ r , r ′ σ + r σ + r ′ + H.c.)] , (13)where σ ± r = ( σ x r ± iσ y r ) / σ α r ( α = x, y, z definedusing the local axes Eqs. (2)-(5)) stand for the Pauli ma-trices of the pseudospin at a site r . The phases φ r , r ′ are φ r , r ′ = 0, − π/
3, and 2 π/ i, i ′ ) = (0 , , (1 , i, i ′ ) = (0 , , (2 , i, i ′ ) = (0 , , (1 , r , r ′ ) = ( t n + d i , t n ′ + d i ′ ).For the magnetic interaction of TTO we probably haveto include the classical dipolar interaction, i.e., H m = H m,NN + Dr × X h r , r ′ i (cid:26) z r · z r ′ | ∆ r | − z r · ∆ r ][ z r ′ · ∆ r ] | ∆ r | (cid:27) σ z r σ z r ′ , (14)where the summation runs over all pairs of sites, r nn isthe NN distance, and ∆ r = r − r ′ . The parameter D isdetermined by the magnetic moment of the CF groundstate doublet. We adopt D = 0 .
29 K, correspondingto the experimental value of the magnetic moment 4.6 µ B23 . As discussed in Refs. 33 and 34, when the mag-netic interaction of Eq. (14) represents the dipolar spinice ( J nn + D nn > H m can be approximated by the NNclassical spin-ice Hamiltonian H m ≃ ( J nn + D nn ) X h r , r ′ i σ z r σ z r ′ , (15)where D nn = D = 0 .
48 K.In our computations we used an effective pseudospin-1/2 Hamiltonian of the form H eff = H m + H q . (16)We note that this is not very different from the origi-nal Onoda-type interaction ( D nn = 0) and results ofRefs. 25 and 7 can be approximately used at least in theelectric quadrupolar phases, in which xy -components ofthe pseudospin ( σ x r , σ y r ) show LRO and semi-classical the-oretical treatments are applicable. IV. PSEUDOPIN WAVE
The studies of the effective Hamiltonian ofEq. (16) showed that there are two electric quadrupolarstates: the PAF state (planar antiferropseudospin) andthe PF state (planar ferropseudospin) depending on thetwo parameters ( δ, q ) (see Fig. 7 in Ref. 25 and Fig. 3in Ref. 7 for details). In these states, the xy -componentsof the pseudospin show LRO with the modulation vector k = 0. It should be noted that this wave vector k = 0 is selected by quantum and thermal fluctuations forPAF, i.e., by an order-by-disorder mechanism.In order to calculate elementary excitations in the PAFand PF states, we choose one of the pseudospin structures( h σ x t n + d i i , h σ y t n + d i i ) = ( (0 , h σ y i ) ( i = 0 , − (0 , h σ y i ) ( i = 1 ,
2) (PAF)(17)and( h σ x t n + d i i , h σ y t n + d i i ) = (0 , h σ y i ) ( i = 0 , , ,
3) (PF) . (18)We apply the simple linear spin-wave theory, MF-RPA (mean field, random phase approximation), in the sameway as described in § h σ y i of Eqs. (17) and (18) is calculated by the MF approx-imation. For the present purpose, we are interested inelementary excitations only at low temperatures, and h σ y i = 1 is a good approximation. To obtain dis-persion relations of pseudospin waves, MF-RPA utilizesthe generalized susceptibility χ ( k , E ) and neutron mag-netic scattering intensity S ( Q , E ) . Useful examples ofMF-RPA computations including straightforward tech-nical extensions for pyrochlore structures are describedRefs. 35 and 36. General computational treatments ofMF-RPA are discussed in Refs. 37 and 38. Followingthese references , the generalized susceptibilityis given by χ ( k , E ) = [1 − χ ( E ) J ( k )] − χ ( E ) . (19)where k is a wave vector in the first Brillouin zone, χ ( E )and J ( k ) denote the single-site generalized-susceptibilityof the MF Hamiltonian and the Fourier transform of theexchange and dipolar coupling constants. The neutronmagnetic scattering intensity S ( Q = G + k , E ) is givenby S ( Q , E ) ∝ − e − βE X ρ,σ ( δ ρ,σ − ˆ Q ρ ˆ Q σ ) × X i,i ′ U ( i ) ρ,z U ( i ′ ) σ,z Im n χ i,z ; i ′ ,z ( k , E ) e − i G · ( d i − d i ′ ) o , (20)where only the local z -component of the pseudospin σ z r contribute to the scattering. If one assumes that all thepseudospin components represent a magnetic dipole mo-ment vector with an isotropic g -factor, virtual neutronscattering intensity S v ( Q , E ) is given by S v ( Q , E ) ∝ − e − βE X ρ,σ ( δ ρ,σ − ˆ Q ρ ˆ Q σ ) × X i,α,i ′ ,α ′ U ( i ) ρ,α U ( i ′ ) σ,α ′ Im n χ i,α ; i ′ ,α ′ ( k , E ) e − i G · ( d i − d i ′ ) o , (21)where U ( i ) ρ,α is the rotation matrix from the local ( α )frame defined at the sites t n + d i to the global ( ρ ) frame. S v ( Q , E ) is useful when displaying dispersion relations Γ X W L XK Γ E ( m e V ) E ( m e V ) (a)(b) 0105 FIG. 1. (Color online) Magnetic S ( Q , E ) (a) and virtual S v ( Q , E ) (b) of the PAF ordering (Eq. (17)) using interac-tion parameters J nn = 1 K, q = 0 .
85, and δ = 0. Γ X W L XK Γ E ( m e V ) E ( m e V ) (a)(b) 0105 FIG. 2. (Color online) Magnetic S ( Q , E ) (a) and virtual S v ( Q , E ) (b) of the PAF ordering (Eq. (17)) using interac-tion parameters J nn = 1 K, q = 0 .
5, and δ = 0 . of all pseudospin waves, because the amplitude of theelectric quadrupole moment are excluded for S ( Q , E ).In Fig. 1(a) we show the inelastic magnetic scatter-ing intensity S ( Q , E ) (Eq. (20)) of the the PAF ordering(Eq. (17)) along several symmetry directions in the FCCBrillouin zone using the interaction parameters J nn = 1K, q = 0 .
85, and δ = 0 adopted in Ref. 23. One can seetwo flat excitation branches in Fig. 1(a). We also show Γ X W L XK Γ E ( m e V ) E ( m e V ) (a)(b) 0105 FIG. 3. (Color online) Magnetic S ( Q , E ) (a) and virtual S v ( Q , E ) (b) of the PF ordering (Eq. (18)) using interactionparameters J nn = 1 K, q = 0 .
8, and δ = − . Γ X W L XK Γ E ( m e V ) E ( m e V ) (a)(b) 0105 FIG. 4. (Color online) Magnetic S ( Q , E ) (a) and virtual S v ( Q , E ) (b) of the PF ordering (Eq. (18)) using interactionparameters J nn = 1 K, q = 0, and δ = − . the virtual S v ( Q , E ) (Eq. (21)) in Fig. 1(b). This figureclearly shows that there are four excitation branches con-sistent with the k = 0 structure posessing four sites in theunit cell. These four pseudospin-wave branches are com-posite spin ( σ z r ) and quadrupole ( σ x r ) waves. Figs. 1(a)and (b) show that the amplitude of the spin compo-nents is strong and weak in the two lower- E and thetwo higher- E branches, respectively. Fig. 2 shows themagnetic S ( Q , E ) and virtual S v ( Q , E ) using differentparameters J nn = 1 K, q = 0 .
5, and δ = 0 . E excitation branchesbecome more dispersive by the finite value of δ comparedto Fig. 1.In Fig. 3 we show the magnetic S ( Q , E ) and virtual S v ( Q , E ) using parameters J nn = 1 K, q = 0 .
8, and δ = − .
6, which are in the PF phase (Eq. (18)). Comparedto the PAF cases, the difference between the magneticand virtual S ( Q , E ) becomes less pronounced. Fig. 4shows the magnetic S ( Q , E ) and virtual S v ( Q , E ) usingparameters J nn = 1 K, q = 0, and δ = − . q = 0, the two lower- E branches are more flattened and merge into almost onebranch. V. MAGNETIC SPECTRA OFPOLYCRYSTALLINE TB x TI − x O y Finally, we would like to compare the previouslyobserved inelastic magnetic neutron scattering spec-tra peaked around E = 0 . x Ti − x O y with x = 0 .
005 ( T c = 0 . .The neutron scattering experiment was performed onthe time-of-flight spectrometer ILL-IN5 operated with λ = 10 ˚A. Fig. 5(b) shows Q -dependent powder spec-tra taken at T = 0 . S ( Q , E ).Fig. 5(a) shows an example of this powder averaged S ( | Q | , E ) choosing the parameters J nn = 1 K, q = 0 . δ = 0, which are in the PAF phase (Eq. (17)). Wethink that these figures show reasonably good agreementbetween the calculation and the observation. In spiteof using the over-simplified model Hamiltonian for TTOand the crude linear-spin-wave theory for the frustratedquantum system, essential features of experimental spec-tra can be reproduced by the approximate calculation.The slight Q -dependence and the non-resolution limitedpeak-width, ∆ E ≫ (∆ E ) resolution = 0 .
01 meV, have beenone of the puzzling observations of TTO. The present in-terpretation using the composite spin-quadrupole wavecan be an answer . VI. SUMMARY
In this study, we try to reformulate the problemof Tb x Ti − x O y and reinterpret its puzzling ex- perimental facts based on the theoretically predicted pseudospin-1/2 Hamiltonian including the electronic su-perexchange interaction between electric quadrupole mo-ments. In this scenario, the hidden order in some TTOsamples is an electric quadrupolar LRO. Although thisLRO does not give rise to strong magnetic Bragg scat-tering, it can be observed by inelastic magnetic neutronscattering as a composite spin-quadrupole wave. We em- E ( m e V ) Q ( Å -1 ) Q ( Å -1 )(b)(a) 0105 FIG. 5. (Color online) (a) Powder averaged magnetic S ( | Q | , E ) of the PAF ordering (Eq. (17)) using interaction pa-rameters J nn = 1 K, q = 0 .
8, and δ = 0. (b) Inelastic neutronscattering spectra of polycrystalline Tb x Ti − x O y with x = 0 .
005 at T = 0 . T c . ploy a MF-RPA linear spin-wave theory and compare itscomputation with previously observed low-energy mag-netic excitation spectra of a polycrystalline sample with x = 0 .
005 ( T c = 0 . . This may possibly implythat Tb x Ti − x O y samples with x < x c are in theU(1) quantum spin-liquid phase . ACKNOWLEDGMENTS
We thank S. Onoda and Y. Kato for useful discussions.This work was supported by JSPS KAKENHI grant num-bers 25400345 and 26400336. The neutron scatteringperformed using ILL-IN5 (France) was transferred fromJRR3-HER (proposal 11567) with the approval of ISSP,Univ. of Tokyo, and JAEA, Tokai, Japan. C. Lacroix, P. Mendels, and F. Mila, eds.,
Introduction toFrustrated Magnetism (Springer, Berlin, Heidelberg, 2011). J. S. Gardner, M. J. P. Gingras, and J. E. Greedan,Rev. Mod. Phys. , 53 (2010). S. T. Bramwell and M. J. P. Gingras, Science , 1495(2001). C. Castelnovo, R. Moessner, and S. L. Sondhi,Nature , 42 (2008). M. Hermele, M. P. A. Fisher, and L. Balents,Phys. Rev. B , 064404 (2004). L. Savary and L. Balents,Phys. Rev. Lett. , 037202 (2012). S. Lee, S. Onoda, and L. Balents,Phys. Rev. B , 104412 (2012). Y. Kato and S. Onoda,Phys. Rev. Lett. , 077202 (2015). M. J. P. Gingras and P. A. McClarty,Rep. Prog. Phys. , 056501 (2014). O. Benton, O. Sikora, and N. Shannon,Phys. Rev. B , 075154 (2012). L.-J. Chang, S. Onoda, Y. Su, Y.-J. Kao, K.-D.Tsuei, Y. Yasui, K. Kakurai, and M. R. Lees,Nature Communications , 992 (2012). K. A. Ross, L. Savary, B. D. Gaulin, and L. Balents,Phys. Rev. X , 021002 (2011). J. S. Gardner, S. R. Dunsiger, B. D. Gaulin, M. J. P. Gin-gras, J. E. Greedan, R. F. Kiefl, M. D. Lumsden, W. A.MacFarlane, N. P. Raju, J. E. Sonier, I. Swainson, andZ. Tun, Phys. Rev. Lett. , 1012 (1999). H. R. Molavian, M. J. P. Gingras, and B. Canals,Phys. Rev. Lett. , 157204 (2007). H. Takatsu, H. Kadowaki, T. J. Sato, J. W.Lynn, Y. Tabata, T. Yamazaki, and K. Matsuhira,J. Phys. Condens. Matter , 052201 (2012). T. Taniguchi, H. Kadowaki, H. Takatsu, B. F˚ak, J. Ol-livier, T. Yamazaki, T. J. Sato, H. Yoshizawa, Y. Shimura,T. Sakakibara, T. Hong, K. Goto, L. R. Yaraskavitch, andJ. B. Kycia, Phys. Rev. B , 060408 (2013). S. Petit, P. Bonville, J. Robert, C. Decorse, and I. Mire-beau, Phys. Rev. B , 174403 (2012). T. Fennell, M. Kenzelmann, B. Roessli, H. Mutka, J. Ol-livier, M. Ruminy, U. Stuhr, O. Zaharko, L. Bovo,A. Cervellino, M. K. Haas, and R. J. Cava,Phys. Rev. Lett. , 017203 (2014). K. Fritsch, K. A. Ross, Y. Qiu, J. R. D. Copley, T. Guidi,R. I. Bewley, H. A. Dabkowska, and B. D. Gaulin,Phys. Rev. B , 094410 (2013). K. Fritsch, E. Kermarrec, K. A. Ross, Y. Qiu, J. R. D.Copley, D. Pomaranski, J. B. Kycia, H. A. Dabkowska,and B. D. Gaulin, Phys. Rev. B , 014429 (2014). S. Petit, S. Guitteny, J. Robert, P. Bonville,C. Decorse, J. Ollivier, H. Mutka, and I. Mirebeau,EPJ Web of Conferences , 03012 (2015). Y. Chapuis, Ph.D. thesis, Universit´e Joseph Fourier(2009), http://tel.archives-ouvertes.fr/tel-00463643/en/. H. Takatsu, S. Kittaka, A. Kasahara, Y. Kono, T.Sakakibara, Y. Kato, S. Onoda, B. F˚ak, J. Ollivier, J.W. Lynn, T. Taniguchi, M. Wakita, and H. Kadowaki,arXiv:1506.04545. M. Wakita, T. Taniguchi, H. Edamoto, H. Takatsu, and H.Kadowaki, arXiv:1509.04583. S. Onoda and Y. Tanaka, Phys. Rev. B , 094411 (2011). M. J. P. Gingras, B. C. den Hertog, M. Faucher, J. S.Gardner, S. R. Dunsiger, L. J. Chang, B. D. Gaulin, N. P.Raju, and J. E. Greedan, Phys. Rev. B , 6496 (2000). I. Mirebeau, P. Bonville, and M. Hennion,Phys. Rev. B , 184436 (2007). A. Bertin, Y. Chapuis, P. Dalmas de R´eotier, andA. Yaouanc, J. Phys. Condens. Matter , 256003 (2012). J. Zhang, K. Fritsch, Z. Hao, B. V. Bagheri, M. J. P. Gin-gras, G. E. Granroth, P. Jiramongkolchai, R. J. Cava, andB. D. Gaulin, Phys. Rev. B , 134410 (2014). A. J. Princep, H. C. Walker, D. T. Adroja,D. Prabhakaran, and A. T. Boothroyd,Phys. Rev. B , 224430 (2015). J. Jensen and A. R. Mackintosh,
Rare Earth Magnetism (Clarendon Press, Oxford, 1991). H. Kusunose, J. Phys. Soc. Jpn. , 064710 (2008). B. C. den Hertog and M. J. P. Gingras,Phys. Rev. Lett. , 3430 (2000). S. V. Isakov, R. Moessner, and S. L. Sondhi,Phys. Rev. Lett. , 217201 (2005). Y.-J. Kao, M. Enjalran, A. Del Maestro, H. R. Molavian,and M. J. P. Gingras, Phys. Rev. B , 172407 (2003). S. Petit, P. Bonville, I. Mirebeau, H. Mutka, andJ. Robert, Phys. Rev. B , 054428 (2012). M. Rotter, Comput. Mater. Sci. , 400 (2006).38