Phase Diagram of the Frustrated Square-Lattice Hubbard Model: Variational Cluster Approach
JJournal of the Physical Society of Japan
FULL PAPERS
Phase Diagram of the Frustrated Square-Lattice Hubbard Model:Variational Cluster Approach
Kazuma Misumi, Tatsuya Kaneko, and Yukinori Ohta
Department of Physics, Chiba University, Chiba 263-8522, Japan
The variational cluster approximation is used to study the frustrated Hubbard model at half fillingdefined on the two-dimensional square lattice with anisotropic next-nearest-neighbor hopping parameters.We calculate the ground-state phase diagrams of the model in a wide parameter space for a variety oflattice geometries, including square, crossed-square, and triangular lattices. We examine the Mott metal-insulator transition and show that, in the Mott insulating phase, magnetic phases with N´eel, collinear, andspiral orders appear in relevant parameter regions, and in an intermediate region between these phases,a nonmagnetic insulating phase caused by the quantum fluctuations in the geometrically frustrated spindegrees of freedom emerges.
1. Introduction
The effect of geometrical frustration in strongly cor-related electron systems has been one of the major is-sues of condensed matter physics. In particular, a spin-liquid state caused by the frustration has been inter-preted as an exotic state of matter, where the magneticlong-range order is destroyed, yielding a quantum para-magnetic (or nonmagnetic) state at zero temperature or even exotic mechanisms of high-temperature super-conductivity. The Hubbard, Heisenberg, and relatedmodels defined on two-dimensional square and triangularlattices with geometrical frustration have been studied inthis respect to find novel quantum disordered states bya variety of theoretical methods.In the square-lattice cases, the J - J Heisenbergmodel with the nearest-neighbor ( J ) and next-nearest-neighbor ( J ) exchange interactions have been stud-ied for more than two decades. At J = 0,where the frustration is absent, the model is knownto have the N´eel-type antiferromagnetic long-rangeorder. With increasing J , the frustration increases,but at J = J , the model again has the groundstate with the collinear antiferromagnetic long-range or-der. The strongest frustration occurs around J /J =0 .
5, where nonmagnetic states such as a valencebond state
4, 6, 8, 10, 11, 14–16, 18, 22, 29) and a spin-liquidstate
12, 24–28) have been suggested to appear, the regionof which has recently been studied further in detail.
31, 32)
The t - t - U Hubbard model with the nearest-neighbor( t ) and next-nearest-neighbor ( t ) hopping parametersand the on-site repulsive interaction U has also beenstudied, where it has been shown that the critical in-teraction strength U c of the metal-insulator transitionincreases monotonically with increasing t /t andthat the ground state has the N´eel order at a small t /t and a collinear order around t = t .
33, 35)
Then, the non-magnetic insulating state appears between these orderedstates.
34, 36)
In the triangular-lattice cases, the anisotropic J - J (cid:48) triangular Heisenberg model has been studied. In theisotropic case ( J = J (cid:48) ), the 120 ◦ spiral ordered phaseis known to be stable. In the anisotropic case, theN´eel order is realized when J (cid:48) /J is small and the spi-ral order is realized around J (cid:48) /J = 1, and betweenthese phases, a dimer ordered phase or a spin-liquidphase
47, 48) has been predicted to appear. The anisotropic t - t (cid:48) - U triangular Hubbard model has also been stud-ied, where it has been shown that U c increaseswith increasing t (cid:48) /t : at a small t (cid:48) /t a metal-insulatortransition occurs from the metallic phase to the N´eelordered phase, whereas at t (cid:48) /t (cid:39)
55, 56)
Recently, the magnetic orders inthe triangular-lattice Heisenberg model with J and J have also been studied, where a nonmagnetic insulatingphase is shown to appear between the spiral and collinearphases. In this paper, motivated by the above developments inthe field, we study the frustrated square-lattice Hubbardmodel at half filling with the isotropic nearest-neighborand anisotropic next-nearest-neighbor hopping parame-ters and clarify the metal-insulator transition, the ap-pearance of possible magnetic orderings, and the emer-gence of a nonmagnetic insulating phase. The searchis made in a wide parameter space including square,crossed-square, and triangular lattices, as well as in weakto strong electron correlation regimes. We use the varia-tional cluster approximation (VCA) based on self-energyfunctional theory (SFT), which enables us to takeinto account the quantum fluctuations of the system, sothat we can study the effect of geometrical frustration a r X i v : . [ c ond - m a t . s t r- e l ] F e b . Phys. Soc. Jpn. FULL PAPERS on the spin degrees of freedom and determine the crit-ical interaction strength for the spontaneous symmetrybreaking of the model. We examine the entire regime ofthe strength of electron correlations at zero temperature,of which little detail is known. In particular, we compareour results in the strong correlation regime with those ofthe Heisenberg model, for which many studies have beencarried out. We also compare our results with those inthe weak correlation limit via the generalized magneticsusceptibility calculation and with those of the classicalHeisenberg model calculation where the quantum spinfluctuations are absent.We will thereby show that magnetic phases with N´eel,collinear, and spiral orders appear in relevant regions ofthe parameter space of our model and that a nonmag-netic insulating phase, caused by the quantum fluctua-tions in the frustrated spin degrees of freedom, emergesin a wide parameter region between the ordered phasesobtained. The orders of the phase transitions will also bedetermined. We will summarize our results as a ground-state phase diagram in a full two-dimensional parameterspace. This phase diagram will make the characterizationof the nonmagnetic insulating phase more approachable,although this is beyond the scope of the present paper.
Fig. 1. (Color online) (a) Schematic representation of thesquare-lattice Hubbard model with the isotropic nearest-neighborhopping parameter t and anisotropic next-nearest-neighbor pa-rameters t and t (cid:48) . (b) Isotropic square lattice at t = t (cid:48) = 0. (c)Crossed square lattice at t = t (cid:48) = t . (d) Isotropic triangular lat-tice at t = t and t (cid:48) = 0. The arrows represent the directions ofthe electron spins. The sublattices are indicated by different colors.
2. Model and Method
We consider the frustrated Hubbard model defined onthe two-dimensional square lattice at half filling as illus- trated in Fig. 1. The Hamiltonian is given by H = − t (cid:88) (cid:104) i,j (cid:105) (cid:88) σ c † iσ c jσ − t (cid:88) (cid:104)(cid:104) i,j (cid:105)(cid:105) (cid:88) σ c † iσ c jσ − t (cid:48) (cid:88) (cid:104)(cid:104) i,j (cid:105)(cid:105) (cid:48) (cid:88) σ c † iσ c jσ + U (cid:88) i n i ↑ n i ↓ − µ (cid:88) i,σ n iσ , (1)where c † iσ is the creation operator of an electron withspin σ at site i and n iσ = c † iσ c iσ . (cid:104) i, j (cid:105) indicates thenearest-neighbor bonds with an isotropic hopping pa-rameter t , and (cid:104)(cid:104) i, j (cid:105)(cid:105) and (cid:104)(cid:104) i, j (cid:105)(cid:105) (cid:48) indicate the next-nearest-neighbor bonds with anisotropic hopping param-eters t and t (cid:48) , respectively [see Fig. 1(a)]. U is the on-site Coulomb repulsion between electrons and µ is thechemical potential maintaining the system at half filling.In the large- U limit, the model can be mapped onto thefrustrated spin-1/2 Heisenberg model H = J (cid:88) (cid:104) i,j (cid:105) S i · S j + J (cid:88) (cid:104)(cid:104) i,j (cid:105)(cid:105) S i · S j + J (cid:48) (cid:88) (cid:104)(cid:104) i,j (cid:105)(cid:105) (cid:48) S i · S j (2)in the second-order perturbation of the hopping parame-ters with S i = (cid:80) α,β c † iα σ αβ c iβ /2, where σ αβ is the vec-tor of Pauli matrices. The exchange coupling constantsare given by J = 4 t /U , J = 4 t /U , and J (cid:48) = 4 t (cid:48) /U for the lattice shown in Fig. 1(a). We will compare ourresults of the Hubbard model in the strong correlationregime with those of the frustrated Heisenberg model,for which related studies have been carried out.We treat a wide parameter space of 0 ≤ t /t ≤ ≤ t (cid:48) /t ≤
1, including three limiting cases: (i)at t = t (cid:48) = 0 [square lattice, see Fig. 1(b)], where theN´eel order is realized, (ii) at t = t (cid:48) = t [crossed squarelattice, see Fig. 1(c)], where the collinear order is realized,and (iii) at t = t and t (cid:48) = 0 [triangular lattice, seeFig. 1(d)], where the 120 ◦ spiral order is realized. Wewill calculate how the above three ordered phases changewhen the hopping parameters are varied in the ranges0 ≤ t ≤ t and 0 ≤ t (cid:48) ≤ t . Fig. 2. (Color online) Left: twelve-site square-lattice cluster usedas a reference system in our analysis. Right: equivalent triangular-lattice cluster, where the three sites 1, 2, and 3 form an equilateraltriangle. The anisotropic triangular lattice is defined as t (cid:48) = 0 and t (cid:54) = t . We employ the VCA, which is a quantum clustermethod based on SFT, where the grand potential
2. Phys. Soc. Jpn.
FULL PAPERS
Ω of the original system is given by a functional of theself-energy. By restricting the trial self-energy to that ofthe reference system Σ (cid:48) , we obtain the grand potentialin the thermodynamic limit asΩ[Σ (cid:48) ] = Ω (cid:48) + Trln( G − − Σ (cid:48) ) − − Trln G (cid:48) , (3)where Ω (cid:48) and G (cid:48) are the exact grand potential and Greenfunction of the reference system, respectively, and G is the noninteracting Green function. The short-rangeelectron correlations within the cluster of the referencesystem are taken into account exactly.The advantage of the VCA is that the spontaneoussymmetry breaking can be treated within the frame-work of the theory. Here, we introduce the Weiss fieldsfor magnetic orderings as variational parameters. TheHamiltonian of the reference system is then given by H (cid:48) = H + H N + H C + H S with H N = h (cid:48) N (cid:88) i e i Q N · r i S zi (4) H C = h (cid:48) C (cid:88) i e i Q C · r i S zi (5) H S = h (cid:48) S (cid:88) i e a i · S i , (6)where h (cid:48) N , h (cid:48) C , and h (cid:48) S are the strengths of the Weiss fieldsfor the N´eel, collinear, and spiral orders, respectively.The wave vectors are defined as Q N = ( π, π ) for the N´eelorder and Q C = ( π,
0) or (0 , π ) for the collinear order.For the spiral order, the unit vectors e a i are rotated by120 ◦ to each other, where a i (= 1 , ,
3) is the sublatticeindex of site i . The variational parameter is optimized onthe basis of the variational principle ∂ Ω /∂h (cid:48) = 0 for eachmagnetic order. The solution with h (cid:48) (cid:54) = 0 corresponds tothe ordered state.We use the twelve-site cluster shown in Fig. 2 as thereference system. This cluster is convenient because wecan treat the two-sublattice states (N´eel and collinearstates) with an equal number of up and down spins and,at the same time, the three-sublattice state (spiral state)with an equal number of three sublattice sites. Note thatlonger-period phases such as a spiral phase mentioned ina different system and incommensurate ordered phasesare difficult to treat in the present approach.
3. Results of Calculations
First, let us discuss the phase diagram of our modelin the strong correlation regime
U/t = 60. The resultis shown in Fig. 3, where the result for our Hubbardmodel in the ( t /t , t (cid:48) /t ) plane as well as the sameresult converted to the Heisenberg model parameters( J /J , J (cid:48) /J ) are shown. We find three ordered phases:the N´eel ordered phase around ( t /t , t (cid:48) /t ) = (0 , t /t , t (cid:48) /t ) = (1 , t /t , t (cid:48) /t ) = (1 , Fig. 3. (Color online) (a) Calculated ground-state phase dia-gram of our model at
U/t = 60 in the ( t /t , t (cid:48) /t ) plane and(b) converted phase diagram in the ( J /J , J (cid:48) /J ) plane. The un-colored region corresponds to the nonmagnetic insulating phase.The transition to the collinear phase is of the first order (or dis-continuous) and the transitions to the N´eel and spiral phases areof the second order (or continuous). The phases along dashed lines(i), (ii), and (iii) shown in (a) are circumstantiated in Fig. 4. and (0 , E = Ω+ µ (per site) and magnetic order parame-ters M (per site) defined as M N = (2 /L ) (cid:80) i e i Q N · r i (cid:104) S zi (cid:105) for the N´eel order, M C = (2 /L ) (cid:80) i e i Q C · r i (cid:104) S zi (cid:105) for thecollinear order, and M S = (2 /L ) (cid:80) i e a i · (cid:104) S i (cid:105) for thespiral order, where (cid:104) (cid:105) stands for the ground-state expec-tation value and L is the number of sites in the system.In the following, we will circumstantiate the obtainedphases, particularly along lines (i), (ii), and (iii) drawnin Fig. 3(a), whereby we will discuss some details of ourcalculated results in comparison with other studies.Along line (i): The results are shown in the left pan-els of Fig. 4, where we assume t = t (cid:48) . At t = 0, the
3. Phys. Soc. Jpn.
FULL PAPERSFig. 4. (Color online) Calculated ground-state energies (upper panels) and order parameters (lower panels) for the N´eel, collinear,spiral, and nonmagnetic insulating phases as a function of t /t or t (cid:48) /t . The left, middle, and right panels correspond to lines (i), (ii),and (iii) in Fig. 3(a), where we assume t = t (cid:48) , t (cid:48) = 0, and t = t , respectively. The inset in (c) and (e) displays the energy differencebetween the spiral and nonmagnetic insulating phases ∆ E , and other insets enlarge the region near the phase boundary. ground state is the N´eel order, and with increasing t ,the energy of the N´eel order gradually approaches theenergy of the nonmagnetic state. At t /t = 0 .
73, theenergy of the N´eel order continuously reaches the energyof the nonmagnetic state and the N´eel order disappears.The calculated order parameter indicates a continuousphase transition. At t /t = 1, on the other hand, theground state is the collinear order. The ground-state en-ergy of the collinear order increases with decreasing t ,and at t /t = 0 .
79, it crosses to the nonmagnetic state,resulting in a discontinuous phase transition, as indicatedby the calculated order parameter. The nonmagnetic in-sulating state thus appears at 0 . < t /t < . . < J /J < . J - J square-lattice Heisenbergmodel, which have estimated the transition point be-tween the N´eel and nonmagnetic phases to be at J /J =0 . − .
17, 22, 24, 31, 32) our result slightly overestimatesthe stability of the N´eel order. This overestimation maybe caused by the cluster geometry used in our calcu-lations; if we use the 2 × J /J = 0 . which is in good agreement with the previous studies.The transition point between the collinear and nonmag-netic phases, on the other hand, has been estimated tobe at J /J = 0 . − .
17, 22, 24, 32) which is in goodagreement with our result.Along line (ii): The results are shown in the mid- dle panels of Fig. 4, where we assume t (cid:48) = 0. Withincreasing t from t = 0, at which the ground stateis the N´eel order, the energy of the N´eel order grad-ually approaches the energy of the nonmagnetic state,and at t /t = 0 .
88, the N´eel order disappears con-tinuously. The calculated order parameter indicates thecontinuous phase transition. At t /t = 1, on the otherhand, the ground state is the spiral order, although theenergy difference between the spiral and nonmagneticstates is very small [see the inset of Fig. 4(c)] due tothe strong geometrical frustration of the triangular lat-tice. With decreasing t from t /t = 1, the ground-stateenergy of the spiral order increases gradually and ap-proaches the energy of the nonmagnetic state, and at t /t = 0 .
89, the spiral order disappears continuously,in agreement with the calculated order parameter. Thus,the nonmagnetic phase appears in a very narrow regionof 0 . < t /t < .
89. The corresponding Heisenbergmodel parameters at which the N´eel and spiral orders dis-appear are around J /J = 0 .
79. The previous studies forthe anisotropic triangular-lattice Heisenberg model
42, 44) have given values around J /J = 0 . − .
87 for thetransition point, which are in good agreement with ourresult.Along line (iii): The results are shown in the right pan-els of Fig. 4, where we assume t = t . At t (cid:48) = 0, theground state is the spiral order, although the energy dif-ference from the nonmagnetic state is very small [see theinset of Fig. 4(e)]. With increasing t (cid:48) , the energy of the
4. Phys. Soc. Jpn.
FULL PAPERSFig. 5. (Color online) Calculated ground-state phase diagrams of our model in the ( t /t , t (cid:48) /t ) plane at (a) U/t = 10, (b) U/t = 6,and (c) U/t = 2. The uncolored regions in (a) and (b) correspond to the nonmagnetic insulating phase. The phases along dashed lines(i), (ii), and (iii) shown in (a) are circumstantiated in Fig. 6. Fig. 6. (Color online) Calculated results for the order parameters M (pink, green, and red dots) of the N´eel, collinear, and spiralphases and the single-particle gap ∆ /t (black crosses) at U/t = 10 (upper panels), U/t = 6 (middle panels), and U/t = 2 (lowerpanels). The left, middle, and right panels correspond to the lines (i), (ii), and (iii) defined in Fig. 5(a), where we assume t = t (cid:48) , t (cid:48) = 0,and t = t , respectively. spiral order gradually approaches the energy of the non-magnetic state, and at t (cid:48) /t = 0 .
34, the spiral order dis-appears continuously, in agreement with the calculatedorder parameter. On the other hand, with decreasing t (cid:48) from t (cid:48) /t = 1, at which the collinear order is stable,the ground-state energy of the collinear order increasesand crosses to the nonmagnetic state at t (cid:48) /t = 0 .
5. Phys. Soc. Jpn.
FULL PAPERS the calculated order parameter. The nonmagnetic statetherefore appears at 0 . < t (cid:48) /t < .
59, which cor-responds to the region 0 . < J (cid:48) /J < .
34 if we usethe Heisenberg model parameters. To our knowledge, nocomparable calculations have been made for the frus-trated Heisenberg model in this parameter region.
Next, let us discuss the phase diagram of our model inthe intermediate to weak correlation regime. The resultsat
U/t = 10, 6, and 2 are shown in Fig. 5. The detailedresults for the calculated single-particle gap and orderparameters are also shown in Fig. 6 along lines (i), (ii),and (iii) defined above.At U/t = 10, we find that the results are qualita-tively similar to those at U/t = 60, except for the tran-sition between the N´eel and spiral orders: the nonmag-netic insulating phase appears between these orders at U/t = 60 but a direct first-order transition occurs at U/t = 10
55, 56) with a double minimum structure in thegrand potential. The nonmagnetic insulating phase ap-pears at 0 . < t /t < .
84 along line (i), which is ingood agreement with the previous studies,
34, 36) wherethe values t /t = 0 . − .
77 for the transition betweenthe N´eel and nonmagnetic phases and t /t = 0 . − . . < t (cid:48) /t < .
64 along line (iii).At
U/t = 6, we find that the spiral phase disap-pears and a paramagnetic metallic phase appears inthe triangular lattice geometry at ( t /t , t (cid:48) /t ) (cid:39) (1 , , . < t /t < .
86 along line (i), the region of whichbecomes wider with decreasing value of
U/t from 60 to10 and 6, which is again in good agreement with theprevious studies.
33, 34)
The region of the nonmagnetic in-sulating phase also becomes wider along line (iii), whichoccurs between the collinear and paramagnetic metallicphases. Along line (ii), the nonmagnetic insulating phasewith a small charge gap appears again, which is betweenthe N´eel and paramagnetic metallic phases.
55, 56) At U/t = 2, the paramagnetic metallic phase over-whelms the collinear and nonmagnetic insulating phases,retaining only the N´eel ordered phase around t /t = t (cid:48) /t = 0. Within the N´eel phase, the charge gap opensonly at 0 < t /t < .
16 and the metallic N´eel orderedphase appears at 0 . < t /t < .
41 along line (i) andat 0 . < t /t < .
47 along line (ii). The perfect Fermisurface nesting at t /t = t (cid:48) /t = 0 and its deforma-tion away from t /t = t (cid:48) /t = 0 are responsible forthese results. The generalized magnetic susceptibility χ ( q ) calculated in the noninteracting limit of our model[Eq. (1)] explains this result (see Appendix B). The tran-sition between the N´eel ordered metallic phase and theparamagnetic metallic phase is continuous along line (i)and discontinuous along line (ii).
4. Summary
We have used the VCA based on SFT to study thetwo-dimensional frustrated Hubbard model at half fillingwith the isotropic nearest-neighbor and anisotropic next-nearest-neighbor hopping parameters. We have particu-larly focused on the effect of geometrical frustration onthe spin degrees of freedom of the model in a wide param-eter space including square, crossed-square, and triangu-lar lattices in a wide range of the interaction strengthat zero temperature. We have thereby investigated themetal-insulator transition, the magnetic orders, and theemergence of the nonmagnetic insulating phase, althoughthe phases with incommensurate orders or with longer-period orders than the cluster size used have not beentaken into account owing to the limitation of the VCA.We have also calculated the ground-state phase diagramof the corresponding classical Heisenberg model as wellas the generalized magnetic susceptibility in the nonin-teracting limit.We have thus determined the ground-state phase dia-gram of the model and found that, in the strong corre-lation regime, magnetic phases with the N´eel, collinear,and spiral orders appear in the parameter space, anda nonmagnetic insulating phase, caused by the effect ofquantum fluctuations in the frustrated spin degrees offreedom, emerges in the wide parameter region betweenthese three ordered phases. We have also found that thephase transition from the N´eel and spiral orders to thenonmagnetic phase is continuous (or a second-order tran-sition), whereas the transition from the collinear orderto the nonmagnetic phase is discontinuous (or a first-order transition). We have compared our results withthe results of the corresponding Heisenberg model calcu-lations that have been made so far and found that theagreement is good whenever the comparison is possible.We have also found that, in the intermediate correla-tion regime, the paramagnetic metallic phase begins toappear in the triangular lattice geometry, which over-whelms the collinear and nonmagnetic insulating phasesin the weak correlation regime, retaining only the N´eelordered phase in the square lattice geometry. We hopethat our results for the phase diagram obtained in thewide parameter space will encourage future studies onthe characterization of the nonmagnetic insulating phaseas well as on its experimental relevance.We thank S. Miyakoshi for useful discussions. T. K.acknowledges support from a JSPS Research Fellowshipfor Young Scientists. This work was supported in partby a Grant-in-Aid for Scientific Research (No. 26400349)from JSPS of Japan.
Appendix A: Ground-State Phase Diagram ofthe Classical Heisenberg Model
Here, we present the ground-state phase diagram ofthe classical Heisenberg model, which is defined as in
6. Phys. Soc. Jpn.
FULL PAPERSFig. A · (Color online) Calculated ground-state phase diagramof the corresponding classical Heisenberg model. Eq. (2) but its quantum spins S i are replaced by theclassical vectors ˜ S , so that quantum fluctuations of thesystem are completely suppressed although the frustra-tive features in the spin degrees of freedom are present.The Hamiltonian is given by H = (cid:80) q J ( q ) ˜ S − q · ˜ S q inmomentum space, where J ( q ) = J (cos q x + cos q y ) + J cos( q x + q y )+ J (cid:48) cos( q x − q y ) . (A · andthe phase diagram is obtained as shown in Fig. A · q = ( π, π )], collinear order [ q =( π, q = ( q, q ) and ( q (cid:48) , − q (cid:48) ) with q = cos − ( − J / J ) and q (cid:48) = cos − ( − J / J (cid:48) )]. There-fore, comparing with the results given in the main text,we may conclude that the quantum fluctuations in thegeometrically frustrated spin degrees of freedom are es-sential in the emergence of the nonmagnetic insulatingphase discussed in the main text. Appendix B: Generalized Susceptibility in theNoninteracting Limit
Here, we present the generalized magnetic susceptibil-ity (or Lindhard function) at zero frequency, χ ( q ) = 1 L (cid:88) k f ( (cid:15) k ) − f ( (cid:15) k + q ) (cid:15) k + q − (cid:15) k , (B · (cid:15) k is the corresponding noninteract-ing band disperson and f ( (cid:15) ) is the Fermi function. Thecalculated results at temperature 0 . t are shown inFig. B ·
1, where we find that a diverging behavior appearsonly at q = ( π, π ) in Fig. B · t /t = t (cid:48) /t = 0 . Fig. B · (Color online) Calculated generalized magnetic sus-ceptibility defined in Eq. (B ·
1) at (a) t /t = t (cid:48) /t = 0 .
0, (b) t /t = 1 . t (cid:48) /t = 0 .
0, (c) t /t = 1 . t (cid:48) /t = 0 .
8, and (d) t /t = t (cid:48) /t = 1 .
0. The corresponding Fermi surface is shown ineach panel. interaction strength U . There are characteristic featuresof χ ( q ) but no other diverging behaviors are found, in-dicating the absence of other magnetic orderings in theweak correlation limit.
1) L. Balents, Nature (London) , 199 (2010).2) P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod. Phys. ,17 (2006).3) P. Chandra and B. Doucot, Phys. Rev. B , 9335(R) (1988).4) S. Sachdev and R. N. Bhatt, Phys. Rev. B , 9323 (1990).5) M. J. de Oliveira, Phys. Rev. B , 6181 (1991).6) A. V. Chubukov and T. Jolicoeur, Phys. Rev. B , 12050(R)(1991).7) J. Oitmaa and Z. Weihong, Phys. Rev. B , 3022 (1996).8) M. E. Zhitomirsky and K. Ueda, Phys. Rev. B , 9007 (1996).9) A. E. Trumper, L. O. Manuel, C. J. Gazza, and H. A. Ceccatto,Phys. Rev. Lett. , 2216 (1997).10) R. R. P. Singh, Z. Weihong, C. J. Hamer, and J. Oitmaa, Phys.Rev. B , 7278 (1999).11) L. Capriotti and S. Sorella, Phys. Rev. Lett. , 3173 (2000).12) L. Capriotti, F. Becca, A. Parola, and S. Sorella, Phys. Rev.Lett. , 097201 (2001).13) G. M. Zhang, H. Hu, and L. Yu, Phys. Rev. Lett. , 067201(2003).14) K. Takano, Y. Kito, Y. Ono, and K. Sano, Phys. Rev. Lett. , 197202 (2003).15) J. Sirker, Z. Weihong, O. P. Sushkov, and J. Oitmaa, Phys.Rev. B , 144422 (2006).16) M. Mambrini, A. L¨auchli, D. Poilblanc, and F. Mila, Phys.Rev. B , 144422 (2006).17) R. Darradi, O. Derzhko, R. Zinke, J. Schulenburg, S. E. Kr¨uger,7. Phys. Soc. Jpn. FULL PAPERS and J. Richter, Phys. Rev. B , 214415 (2008).18) L. Isaev, G. Ortiz, and J. Dukelsky, Phys. Rev. B , 024409(2009).19) K. S. D. Beach, Phys. Rev. B , 224431 (2009).20) J. Ritcher and J. Schulenburg, Eur. Phys. J. B , 117 (2010).21) J. Reuther and P. W¨olfle, Phys. Rev. B , 144410 (2010).22) J. F. Yu and Y. J. Kao, Phys. Rev. B , 094407 (2012).23) L. Wang, Z. C. Gu, F. Verstraete, and X. G. Wen,arXiv:1112.3331.24) H.-C. Jiang, H. Yao, and L. Balents, Phys. Rev. B , 024424(2012).25) F. Mezzacapo, Phys. Rev. B , 045115 (2012).26) T. Li, F. Becca, W. J. Hu, and S. Sorella, Phys. Rev. B ,075111 (2012).27) L. Wang, D. Poilblanc, Z. C. Gu, X. G. Wen, and F. Verstraete,Phys. Rev. Lett. , 037202 (2013).28) W. J. Hu, F. Becca, A. Parola, and S. Sorella, Phys. Rev. B , 060402 (2013).29) R. L. Doretto, Phys. Rev. B , 104415 (2014).30) Y. Qi and Z. C. Gu, Phys. Rev. B , 235122 (2014).31) S.-S. Gong, W. Zhu, D. N. Sheng, O. I. Motrunich, and M. P.A. Fisher, Phys. Rev. Lett. , 027201 (2014).32) S. Morita, R. Kaneko, and M. Imada, J. Phys. Soc. Jpn. ,024720 (2015).33) T. Mizusaki and M. Imada, Phys. Rev. B , 014421 (2006).34) A. H. Nevidomskyy, C. Scheiber, D. S´en´echal, and A.-M. S.Tremblay, Phys. Rev. B , 064427 (2008).35) Z.-Q. Yu and L. Yin, Phys. Rev. B , 195122 (2010).36) S. Yamaki, K. Seki, and Y. Ohta, Phys. Rev. B , 125112(2013).37) P. Sindzingre, P. Lecheminant, and C. Lhuillier, Phys. Rev. B , 3108 (1994).38) J. Merino, R. H. McKenzie, J. B. Marston, and C. H. Chung,J. Phys.: Condens. Matter , 2965 (1999).39) Z. Weihong, R. H. McKenzie, and R. R. P. Singh, Phys. Rev.B , 14367 (1999).40) A. E. Trumper, Phys. Rev. B , 2987 (1999).41) S. Yunoki and S. Sorella, Phys. Rev. B , 014408 (2006); D.Heidarian, S. Sorella, and F. Becca, Phys. Rev. B , 012404(2009).42) M. Q. Weng, D. N. Sheng, Z. Y. Weng, and R. J. Bursill, Phys.Rev. B , 012407 (2006); A. Weichselbaum and S. R. White,Phys. Rev. B , 245130 (2011).43) O. A. Starykh and L. Balents, Phys. Rev. Lett. , 077205(2007). 44) R. F. Bishop, P. H. Y. Li, D. J. J. Farnell, and C. E. Campbell,Phys. Rev. B , 174405 (2009).45) J. Reuther and R. Thomale, Phys. Rev. B , 024402 (2011).46) A. Weichselbaum and S. R. White, Phys. Rev. B , 245130(2011).47) P. Hauke, T. Roscilde, V. Murg, J. Cirac, and R. Schmied, NewJ. Phys. , 075017 (2011).48) P. Hauke, Phys. Rev. B , 014415 (2013).49) J. Merino, M. Holt, and B. J. Powell, Phys. Rev. B , 245112(2014).50) H. Morita, S. Watanabe, and M. Imada, J. Phys. Soc. Jpn. ,2109 (2002).51) P. Sahebsara and D. Se´ne´chal, Phys. Rev. Lett. , 257004(2006).52) T. Watanabe, H. Yokoyama, Y. Tanaka, and J. Inoue, Phys.Rev. B , 214505 (2008).53) L. F. Tocchio, A. Parola, C. Gros, and F. Becca, Phys. Rev.B , 064419 (2009).54) L. F. Tocchio, H. Feldner, F. Becca, R. Valenti, and C. Gros,Phys. Rev. B , 035143 (2013).55) A. Yamada, Phys. Rev. B , 195108 (2014).56) M. Laubach, R. Thomale, C. Platt, W. Hanke, and G. Li, Phys.Rev. B , 245125 (2015).57) R. V. Mishmash, J. R. Garrison, S. Bieri, and C. Xu, Phys.Rev. Lett. , 157203 (2013).58) R. Kaneko, S. Morita, and M. Imada, J. Phys. Soc. Jpn. ,093707 (2014).59) P. H. Y. Li, R. F. Bishop, and C. E. Campbell, Phys. Rev. B , 014426 (2015).60) Z. Zhu and S. R. White, Phys. Rev. B , 041105(R) (2015).61) W.-J. Hu, S.-S. Gong, W. Zhu, and D. N. Sheng, Phys. Rev.B , 140403(R) (2015).62) M. Potthoff, M. Aichhorn, and C. Dahnken, Phys. Rev. Lett. , 206402 (2003).63) M. Potthoff, Eur. Phys. J. B , 429; , 335 (2003).64) C. Dahnken, M. Aichhorn, W. Hanke, E. Arrigoni, and M.Potthoff, Phys. Rev. B , 245110 (2004).65) J. M. Luttinger and L. Tisza, Phys. Rev. , 954 (1946).66) J. E. Hirsch, Phys. Rev. B , 4403 (1985).67) N. Bulut, D. J. Scalapino, and S. R. White, Phys. Rev. B ,2742 (1993).68) A. Sherman and M. Schreiber, Phys. Rev. B76