Spin liquids in graphene

Spin liquids in graphene Minh-Tien Tran1,2 and Ki-Seok Kim1,3 1 Asia Pacific Center for Theoretical Physics, Pohang, Gyeongbuk 790-784, Republic of Korea Institute of Physics, Vietnamese Academy of Science and Technology, P.O.Box 429, 10000 Hanoi, Vietnam 3 Department of Physics, Pohang University of Science and Technology, Pohang, Gyeongbuk 790-784, Korea arXiv:1011.1700v2 [cond-mat.str-el] 24 Jan 2011 2 We reveal that local interactions in graphene allow novel spin liquids between the semi-metal and antiferromagnetic Mott insulating phases, identified with algebraic spin liquid and Z2 spin liquid, respectively. We argue that the algebraic spin liquid can be regarded as the two dimensional realization of one dimensional spin dynamics, where antiferromagnetic correlations show exactly the same power-law dependence as valence bond correlations. Nature of the Z2 spin liquid turns out to be d + id′ singlet pairing, but time reversal symmetry is preserved, taking d + id′ in one valley and d − id′ in the other valley. We propose the quantized thermal valley Hall effect as an essential feature of this gapped spin liquid state. Quantum phase transitions among the semi-metal, algebraic spin liquid, and Z2 spin liquid are shown to be continuous while the transition from the Z2 spin liquid to the antiferromagnetic Mott insulator turns out to be the first order. We emphasize that both algebraic spin liquid and d ± id′ Z2 spin liquid can be verified by the quantum Monte Carlo simulation, showing the enhanced symmetry in the algebraic spin liquid and the quantized thermal valley Hall effect in the Z2 spin liquid. PACS numbers: I. INTRODUCTION Interplay between the topological band structure and interaction drives one direction of research in modern condensed matter physics,1,2 where emergence of Dirac fermions is at the heart of the interplay. The original example is the system of one dimensional interacting electrons, where interactions become enhanced at low energies, combined with one dimensionality, and electron fractionalization results, giving rise to a new state of matter, dubbed as the Tomonaga-Luttinger liquid.3 An interesting aspect is that such fractionalized excitations as spinons and holons are identified with topological excitations, carrying fermion quantum numbers associated with the topological structure of the Dirac theory.4 A recent study based on the quantum Monte Carlo simulation5 has argued that essentially the same phenomenon as electron fractionalization in the TomonagaLuttinger liquid may happen in two dimensions when local interactions are introduced in the graphene structure. This study claimed existence of a paramagnetic Mott insulator with a spin gap between the semi-metal and antiferromagnetic Mott insulating phases. Immediately, the nature of the spin gapped Mott state has been suggested to be an s-wave spin-singlet pairing order between next nearest neighbor spins,6–8 thus identified with a Z2 spin liquid. We point out other scenarios9,10 for the nature of the spin liquid state. In the present study we revisit this problem, the nature of possible spin liquids in the Hubbard model on the graphene structure. An important point of our study is to solve the Hubbard model directly beyond recent studies,6–8 where an additional energy scale was introduced to describe the spin-singlet pairing order. The SU(2) slave-rotor representation, invented by one of the authors,11 is at the heart of the methodology, where ex- change correlations via virtual processes are naturally caught to allow spin singlet-pairing. One may regard the SU(2) slave-rotor theory of the Hubbard model as an analogue of the SU(2) slave-boson theory12 for the t-J model. We find two kinds of spin liquids, identified with an algebraic spin liquid and a Z2 spin liquid, respectively, between the semi-metal and antiferromagnetic phases. We argue that the algebraic spin liquid13,14 can be regarded as the two dimensional realization of one dimensional spin dynamics, where antiferromagnetic correlations show exactly the same power-law dependence as valence bond correlations.15,16 Increasing interactions, pairing correlations between nearest neighbor sites become enhanced. As a result, the algebraic spin liquid is shown to turn into a gapped spin liquid state, where the spin gap results from d + id′ singlet pairing, believed to originate from the interplay between the topological structure and interaction. Actually, this pairing symmetry solution has been argued to be stable, based on an effective model in the weak coupling approach.17,18 However, we argue that time reversal symmetry is preserved, taking the d + id′ singlet pairing to one valley while the d − id′ pairing to another. This assignment turns out to be essential in order to have a fully gapped spectrum because the d + id′ singlet pairing order parameter in one valley vanishes in the other valley, allowing the gapless Dirac spectrum. We propose the quantized thermal valley Hall effect19,20 for the fingerprint of this gapped Z2 spin liquid. We discuss the nature of quantum phase transitions beyond the mean-field analysis. We would like to emphasize that appearance of both algebraic spin liquid and d ± id′ Z2 spin liquid can be verified by the quantum Monte Carlo simulation in principle. The fingerprint of the algebraic spin liquid is an enhanced symmetry, giving rise to the same power-law 2 dependence between antiferromagnetic and valence bond correlations. The hallmark of the d ± id′ Z2 spin liquid is the quantized thermal valley Hall effect, as mentioned above. We hope that the present study motivates quantum Monte Carlo simulation researchers to calculate such correlation functions. The present paper is organized as follows. In Sec. II we present the SU(2) slave-rotor theory of the Hubbard model. General mean field equations are also obtained. The mean field analysis of possible phase transitions is presented in Sec. III. In Sec. IV summary and discussion are presented. II. SU(2) SLAVE-ROTOR THEORY OF THE HUBBARD MODEL A. Formulation i ciσ (c†iσ ) is the annihilation (creation) operator where for an electron at site i with spin σ. t is the hopping integral, and U is the on-site Coulomb interaction, where niσ = c†iσ ciσ represents the density of electrons with spin σ. Introducing the Nambu-spinor representation   ci↑ , ψi = c†i↓ and performing the Hubbard-Stratonovich transformation for the pairing, density (singlet) and magnetic (triplet) interaction channels, we obtain an effective Lagrangian X † X † ψi σz ψj + H.c. ψi (∂τ 1 − µσz )ψi − t L= i hiji X † † † I −i [ΦR i (ψi σx ψi ) + Φi (ψi σy ψi ) + ϕi (ψi σz ψi )] i 3 X R 2 + [(Φi ) + (ΦIi )2 + (ϕi )2 ] 2U κc i 1 X 2 X mi (ψi† ψi ). mi − + 2U κs i i   ψi = Zi† Fi , Here, ziσ is a boson operator, satisfying the unimodular † † (rotor) constraint, zi↑ zi↑ + zi↓ zi↓ = 1. The key point of the slave-rotor representation21 is to extract out collective charge dynamics explicitly from correlated electrons. Such charge fluctuations are identified with zero sound modes in the case of short range interactions while plasmon modes in the case of long range interactions. Actually, one can check that the dispersion of the rotor variable (zi↑ ) is exactly the same as that of such collective charge excitations. In the slave-rotor theory the Mott transition is described by gapping of rotor excitations.21 Until now, the Mott transition has not been achieved successfully, based on the diagrammatic approach starting from the Fermi liquid theory. In this respect the slave-rotor theory can be regarded as an effective field theory, not bad for the Mott transition. Resorting to the SU(2) slave rotor representation in Eq. (3), we rewrite the effective Lagrangian Eq. (2) as follows L = L0 + LF + LZ , X † Tr[Xij Yij† + Yij Xij ]+ L0 = t hiji (2) LF = R(I) Φi and ϕi are associated with pairingHere, fluctuation and density-excitation potentials, introduced to decouple the charge channel. mi is an effective magnetic field, which decouples the spin channel. κc and κs are introduced for decoupling between singlet and triplet interactions in respect that we do not know how they become renormalized at low energies. One may regard these (3) fi↑ is a fermion operator in the Nambu † fi↓ representation, and Zi is an SU(2) matrix ! † zi↑ −zi↓ Zi = . (4) † zi↓ zi↑ where Fi = We start from the Hubbard model on the honeycomb lattice X X † ni↑ ni↓ , (1) ciσ cjσ + H.c. + U H = −t hijiσ two decoupling coefficients as phenomenological parameters to overcome the mean-field approximation. Several examples for decoupling are shown in appendix A. We emphasize that both the semi-metal to algebraic spin liquid and the algebraic spin liquid to Z2 spin liquid phase transitions are shown not to depend on such phenomenological parameters, where both the critical points are determined with the combination between U/t and κc . Only the Z2 spin liquid to antiferromagnetic transition turns out to depend on such parameters. We should be careful in determining this phase transition, comparing various cases (appendix A) with each other. The SU(2) slave-rotor representation11 means to write down an electron field as a composite field in terms of a charge-neutral spinon field and a spinless holon field X i −t LZ = Fi† (∂τ 1 − iΩi ·σ)Fi (5) 1 X 2 mi , (6) 2U κs i X X † mi (Fi† Fi ), (7) (Fi Xij Fj + H.c.) − i hiji X 3 Tr[Ωi ·σ−iZi ∂τ Zi† + iµZi σz Zi† ]2 4U κc i X Tr[Zi σ z Zj† Yij† + H.c.]. (8) −t hiji 3 a2 a d1 d3 d2 a1 FIG. 1: Graphene. A distance between two next nearest neighbor sites is chosen as the length unit. a1 and a2 are primitive translation vectors. δ1 , δ2 , and δ3 are three nearest neighbor bonds. It is not difficult to see equivalence between the SU(2) slave-rotor effective Lagrangian and the decoupled Hubbard model [Eq. (2)]. Integrating over field variables of Xij and Yij , and shifting Ωi ·σ as Ωi ·σ+iZi ∂τ Zi† − iµZi σz Zi† , I where Ωi = (ΦR i , Φi , ϕi ) is the pseudospin potential field, we recover Eq. (2) exactly with an introduction of an electron field Zi† Fi → ψi . This procedure is well described in the previous study.11 An important feature in the SU(2) slave-rotor description is appearance of pairing correlations between nearest neighbor electrons, given by off diagonal hopping in Xij which results from on-site pairing (virtual) fluctuations, captured by the off diago- F =− nal variable zi↓ of the SU(2) matrix field Zi . We note that the diagonal rotor field zi↑ corresponds to the zero sound mode, giving rise to the Mott transition via gapping of their fluctuations. The additional boson rotor variable zi↓ allows us to catch super-exchange correlations in the Mott transition. But, the appearance of pairing correlations does not necessarily lead to superconductivity because their global coherence, described by condensation of SU(2) matrix holons, is not guaranteed. The similar situation happens in the SU(2) slave-boson theory12 of the t-J model. B. Mean-field ansatz We perform the mean-field analysis, taking the following ansatz   wδ vδ∗ Xij = · σz , (9) vδ −wδ∗   w eδ veδ∗ Yij = · σz , (10) veδ −w eδ∗ where δ denotes the bond between the nearest neighbor sites. In the honeycomb lattice there are three nearest neighbor bonds. We choose wδ = wγδ , vδ = vζδ , w eδ = wγ e δ , and veδ = veζδ , where γδ and ζδ are symmetric factors for the hopping parameter w (w) e and the pairing order parameter v (e v ). The choice for γδ and ζδ depends on the symmetry of the considered phase. For example, the s-wave pairing symmetry is given by ζδ = (1, 1, 1), the dx2 −y2 -wave symmetry ζδ = (− 21 , − 21 , 1), and the dxy wave symmetry ζδ = (− 21 , − 21 , 0). For the magnetic order parameter mi we choose an antiferromagnetic ansatz mi = (−1)i m. Particle-hole symmetry at half filling results in µ + I iϕi = 0 while pairing potentials of ΦR i and Φi vanish in the mean-field level. Then, we obtain a general expression for the free energy 1X 1X log[(iω)2 − t2 w2 |γ(k)|2 − (tv|ζ(k)| + m)2 ] − log[(iω)2 − t2 w2 |γ(k)|2 − (tv|ζ(k)| − m)2 ] β β k,iω k,iω X 2X 3 N (iν)2 + λ)2 − 4t2 w e2 |γ(k)|2 − 4t2 e v 2 |ζ(k)|2 ] + 4tN m2 − N λ, (11) + log[(− (ww|γ e δ |2 + ve v |ζδ |2 ) + β 4κc U 2κs U δ k,iν P where γ(k) = δ γδ exp(irδ · k) is thePenergy dispersion for spinons and holons, and ζ(k) = δ ζδ exp(irδ · k) is associated with the pairing potential. δ is performed in the unit cell. λ is a Lagrange multiplier field, introduced to keep the slave-rotor constraint. N is the total number of sites. Minimizing the free energy, we obtain fully self-consistent equations for order parameters   X 2tw X |γ(k)|2 |γ(k)|2 2 w e |γδ | = − , (12) + 4N β (iω)2 − t2 w2 |γ(k)|2 − (tv|ζ(k)| + m)2 (iω)2 − t2 w2 |γ(k)|2 − (tv|ζ(k)| − m)2 P δ ve X δ k,iω |ζδ |2 = −   (tv|ζ(k)| − m)|ζ(k)| 2 X (tv|ζ(k)| + m)|ζ(k)| , + 4N β (iω)2 − t2 w2 |γ(k)|2 − (tv|ζ(k)| + m)2 (iω)2 − t2 w2 |γ(k)|2 − (tv|ζ(k)| − m)2 k,iω (13) 4 w X |γδ |2 = X |ζδ |2 = δ v δ m=− |γ(k)|2 4tw eX , 3 2 2 N (− 4κc U (iν) + λ) − 4t2 w e2 |γ(k)|2 − 4t2 ve2 |ζ(k)|2 k (14) |ζ(k)|2 4te vX , 3 N (− 4κc U (iν)2 + λ)2 − 4t2 w e2 |γ(k)|2 − 4t2 ve2 |ζ(k)|2 (15) k 2κs U X tv|ζ(k)| + m tv|ζ(k)| − m 2κs U X + , 2 2 2 2 2 Nβ (iω) − t w |γ(k)| − (tv|ζ(k)| + m) Nβ (iω)2 − t2 w2 |γ(k)|2 − (tv|ζ(k)| − m)2 k,iω (16) k,iω 1= (− 4κ3c U (iν)2 + λ) 4 X . Nβ e2 |γ(k)|2 − 4t2 ve2 |ζ(k)|2 (− 4κ3c U (iν)2 + λ)2 − 4t2 w k,iν (17) In this study our objective is to reveal the phase structure of the Hubbard model on the honeycomb lattice. It is convenient to take the zero temperature limit. Performing the Matsubara frequency summation, we obtain selfconsistent mean-field equations at zero temperature w e ve X δ w |ζδ |2 = X δ v X δ X w X |γ(k)|2 w X |γ(k)|2 + , 8N/2 D(k,m) 8N/2 D(k, −m) (18) 1 X (v|ζ(k)| − m 1 X (v|ζ(k)| + m t )|ζ(k)| t )|ζ(k)| + , 8N/2 D(k,m) 8N/2 D(k, −m) (19) δ |γδ |2 = k k |γδ | = r |ζδ |2 = r 2 k m= k κb U w e X |γ(k)|2 3 2N/2 E(k) k κb U ve X |ζ(k)|2 3 2N/2 E(k) k 1 p −p λ − 2tE(k) λ + 2tE(k) 1 1 p −p λ − 2tE(k) λ + 2tE(k) ! ! , (20) , (21) κs U X v|ζ(k)| + m κs U X v|ζ(k)| − m t t − , 2N/2 D(k,m) 2N/2 D(k, −m) k 1= 1 r (22) k κb U 1 X 3 N/2 k 1 1 p +p λ − 2tE(k) λ + 2tE(k) ! , (23) count spatially uniform hopping where p e2 |γ(k)|2 + e v 2 |ζ(k)|2 , E(k) = w r m D(k, m) = w2 |γ(k)|2 + (v|ζ(k)| + )2 t (24) (25) are holon and spinon energy spectra in the presence of pairing and antiferromagnetism, respectively. Considering symmetry, it is natural to take into ac- √ 1 3 |γ(k)|2 = 3 + 2 cos(ky ) + 4 cos( ky ) cos( kx ). (26) 2 2 On the other hand, the s-wave pairing potential is not allowed due to repulsive interactions. Counting the lattice symmetry of the honeycomb structure, the next candidate will be dx2 −y2 or dxy for nearest neighbor singlet pairing18 . We introduce a general combination of dx2 −y2 and dxy -wave pairing for the pairing term ζ(k) |ζ(k)|2 = | cos(θ)ζx2 −y2 (k) + i sin(θ)ζxy (k)|2 , (27) 5 where θ is a combination factor, and ζx2 −y2 (ζxy ) is the dx2 −y2 (dxy ) -wave symmetry function ζx2 −y2 (kx , ky ) = e kx −i √ 3 ζxy (kx , ky ) = ie i 2k√x3 −e i 2k√x3 sin( cos( ky ) 2 ky ). 2 For θ = ±π/3 this pairing symmetry becomes d ± id′ . We also consider the d + d′ -wave pairing symmetry |ζ(k)|2 = | cos(θ)ζx2 −y2 (k) + sin(θ)|ζxy (k)|2 , (28) but this pairing order turns out to be not a solution of the mean-field equations. If one tunes κc and κs parameters, he can make this pairing symmetry a solution. However, this solution does not give the lowest free energy, compared with the d + id′ pairing solution, consistent with earlier studies.17,18 One may criticize the ansatz for uniform hopping in this paper because such an assumption excludes possible dimerized phases a priori. Actually, the J1 − J2 Heisenberg model X X Sk · Sl S i · S j + J2 H = J1 hiji hhklii has shown several types of dimerized phases when the ratio of J2 /J1 is beyond a certain critical value8,22 , approximately given by J2 /J1 ≈ 0.2 ∼ 0.3. Here, the first term represents the exchange interaction between nearest neighbor spins, and the second expresses that between next nearest neighbor ones. This model can be derived from the Hubbard model, resorting to the degenerate perturbation theory in the t/U → 0 limit23 , where each parameter is given by8 J1 = 4t  t 3 o nt , −4 U U J2 = 4t  t 3 U antiferromagnetic order disappears, and a paramagnetic Mott insulating state results, identified with a certain type of Z2 spin liquids. Such a spin-gapped state turns out to evolve into a dimerized phase with either translational or rotational symmetry breaking near J2 /J1 ≈ 0.2 ∼ 0.3.8,22 It was reported that the spin liquid state turns into a dimerized phase with three-fold degeneracy around J2 /J1 ≈ 0.3, where it breaks the C3 symmetry but preserving the translational symmetry.8 On the other hand, the plaquette order was claimed to appear near J2 /J1 ≈ 0.2 before the dimerized phase, breaking the translational symmetry only.22 An important point is that if we translate the critical J2 /J1 value in terms of U/t of the Hubbard model, J2 /J1 ≈ 0.3 corresponds to U/t ≈ 2.7 and J2 /J1 ≈ 0.2, U/t ≈ 3.0. Comparing these critical values with the critical value for the semi-metal to spin liquid transition in the quantum Monte Carlo simulation5 , the Mott critical value given by U/t ≈ 3.5 turns out to be larger than those for dimerized phases. This means that the semi-metal phase will appear before reaching such dimerized phases in the Hubbard model owing to charge fluctuations, not introduced in the J1 −J2 Heisenberg model. In other words, the J1 − J2 model seems to be an effective low energy model only in the limit of U/t → ∞ while physics of such a model will be different from that of the Hubbard model in the small U/t case. However, it should be pointed out that these critical values cannot be guaranteed. Thus, we cannot exclude the possibility of dimerization near the Mott criticality of the Hubbard model completely. In addition, introduction of the next nearest neighbor hopping t′ will favor dimerization. In this respect it will be the best interpretation that the spin liquid physics may appear at finite temperatures at least, actually observed from the quantum Monte Carlo simulation. III. up to the fourth order process. Then, the J2 /J1 ratio can be expressed in terms of U/t as follows J2 1 = . J1 (U/t)2 − 4 It was argued that higher order terms such as third neighbor and ring exchange terms may be ignored because third neighbor exchange terms are not frustrating, just renormalizing the J1 term effectively, while the ring exchange term is expected to be small.8 However, the role of the ring exchange term has been also studied carefully.24,25 An antiferromagnetic phase has been reported in J2 /J1 < (J2 /J1 )AF ≈ 0.08.8,22 This corresponds to (U/t)AF ≈ 4.3, consistent with the result of the quantum Monte Carlo simulation.5 Increasing frustration, the A. SADDLE-POINT ANALYSIS From semi-metal to algebraic spin liquid The semi-metal phase is described by condensation of holons hzσ i 6= 0 with v = ve = m = 0. Considering the symmetry factor γδ = (1, 1, 1), the condensation occurs when the effective chemical potential given by the Lagrange multiplier field λ touches the maximum point of the holon dispersion, i.e., λc = 2tw ec max |γ(k)| = 6tw ec . These collective charge excitations become gapped when λ > λc , and a Mott insulating state appears, characterized by hzσ i = 0 with v = e v = m = 0. Taking λ = λc with v = ve = m = 0, we can determine the quantum critical point from the following mean-field equations 6 wc X δ |γδ |2 = 1= r r w ec X δ |γδ |2 = k κc U c 1 X 6tw ec N/2 k k 3−γ(k) (29) k κc U c 1 X |γ(k)| 6tw ec 2N/2 1 1 p −p 3 − γ(k) 3 + γ(k) 1 1 p +p 3 − γ(k) 3 + γ(k) Inserting w ec from Eq. (29) into Eq. (31), one obtains the critical value for the interaction strength P 1 |γ(k)| N/2 κc U c 3 k = 2   t 2P P 1 √ 1 +√ 1 |γδ |2 N/2 δ 1 X |γ(k)|, 4N/2 ! . ! , (30) (31) and the gapped spin liquid state because such a simulation should be performed at finite temperatures. But, calculations for correlation functions need not be done at zero temperature. It is sufficient to show equivalent correlation behaviors in the quantum critical region at finite temperatures. 3+γ(k) = 0.312. (32) It is interesting to notice that the resulting paramagnetic Mott insulator has all kinds of lattice symmetries. In particular, spin dynamics is described by gapless spinons. An effective field theory for spinon dynamics was proposed to be an SU(2) gauge theory with Dirac fermions.13 It is not at all straightforward to understand dynamics of such gapless spinons due to complexity of the SU(2) gauge theory. It has been shown that an interacting stable fixed point arises in the large-Nf limit,13 where Nf is the number of fermion flavors. Such a conformal invariant fixed point was also shown to appear in the U(1) gauge theory with gapless Dirac fermions.26 An interesting property of the stable fixed point is that the symmetry of the original microscopic model, here the Hubbard model, is enhanced, associated with special transformation properties of Dirac spinors.14,15 As a result, spin-spin correlations at an antiferromagnetic wave vector have exactly the same power-law dependence as valence bond-valance bond correlations, which means that the scaling dimension of the staggered spin operator is the same as that of the valence bond operator.16 This situation is completely unusual because scaling dimensions of these two operators cannot be the same in the level of the microscopic model. It is clear that one direct way to verify the algebraic spin liquid state is to observe the symmetry enhancement at low energies. If the staggered-spin correlation function turns out to display the same power-law behavior as the valence-bond correlation function, this will be an undisputable evidence for the algebraic spin liquid phase between the semi-metal phase and gapped spin liquid state. In the recent quantum Monte Carlo simulation data there seems to be uncertainty between the semi-metal phase B. From algebraic spin liquid to Z2 spin liquid Increasing κctU more from the semi-metal to algebraic spin liquid critical point κctUc , we find another paramagnetic Mott insulating phase, characterized by v 6= 0 and e v 6= 0 with the d ± id′ pairing symmetry. Recall Eq. (27) and Eq. (28) for pairing symmetries that we checked explicitly. The algebraic spin liquid (hzσ i = v = ve = m = 0) to gapped spin liquid (v 6= 0 and e v 6= 0 with hzσ i = m = 0) critical point is found with an ansatz of v = ve = 0 but v/e v ≡ ϑ 6= 0. The FIG. 2: (Color online) Spinon spectra ±D(k) given by Eq. (25) with w = 1, v = 0.1 and m = 0. Left figure: If we consider only the d + id′ pairing symmetry, spinon excitations are gapped only at the K point, but they remain gapless at the other K ′ point. Right figure: If we consider d + id′ in one valley (K) and d − id′ in the other valley (K ′ ), spinon excitations become fully gapped. In this figure we assign d + id′ to the upper band and d − id′ to the lower band, respectively, in order to realize the d ± id′ pairing symmetry. 7 mean-field equations to determine this critical point are given in appendix B. The system of equations is solved numerically, explicitly shown in appendix B. We find that the free energy reaches the lowest value for the d ± id′ pairing symmetry. Actually, we checked self-consistency for various values of the angle parameter θ in Eqs. (27) and (28), and found θ = ±π/3 in Eq. (27), corresponding to d ± id′ . In Fig. 2 we plot the spinon spectrum D(k) given by Eq. (25) for the d ± id′ pairing symmetry. It shows that if only the d + id′ pairing order parameter is taken into account, the energy spectrum opens a gap at one Brillouin zone edge (for instance, the K point), but it still keeps the Dirac cone at the other inequivalent Brillouin zone edge (the K ′ point). On the other hand, if we consider only the d − id′ pairing symmetry, we see that the K point remains gapless while only the K ′ point becomes gapped. This demonstration motivates us to assign the d + id′ pairing symmetry to one valley (K) and the d − id′ to the other (K ′ ), making the spinon spectrum fully gapped. Of course, this fully gapped state is energetically more favorable than the gapless state. In addition, this proposal resolves the problem of time reversal symmetry breaking at the same time. The edge state from the d+id′ pairing in one valley is cancelled by that from the d − id′ pairing in the other valley, preserving time reversal symmetry. One may regard this cancellation of such edge states as anomaly cancellation due to fermion doubling in condensed matter physics. The critical value turns out to be κc Uv /t = 0.315. Note that Uv > Uc . This intermediate phase between the semimetal and gapped spin liquid is the algebraic spin liquid with an enhanced symmetry, as discussed in the previous subsection. We would like to emphasize that this region of Uc < U < Uv is not wide at zero temperature. But, the quantum critical region at finite temperatures will not be so narrow, and it will not be so difficult to verify the algebraic spin liquid, considering staggered-spin correlations and valence-bond correlations. This pairing state can be verified by the quantized thermal valley Hall effect.19,20 The spinon number is not conserved due to particle-particle pairing, thus the charge Hall conductivity is not useful. On the other hand, both spin and energy (thermal) Hall coefficients are important probes. But, the spin Hall conductivity vanishes due to the different assignment between two valleys. The thermal valley Hall effect should be observed in this state, regarded as the fingerprint of our Z2 spin liquid phase. This may be verified27 by either quantum Monte Carlo simulation28 or exact diagonalization29. It should be noted that our time reversal symmetry preserving Z2 spin liquid state is beyond the classification scheme based on the projective symmetry group because their possible Z2 spin liquids in the projective symmetry group are constrained with complete time reversal symmetric pairing.6,7 In other words, d ± id′ singlet pairing orders are excluded from the first although these pairing orders are not only found but also argued to be stable in recent studies.17,18 C. From Z2 spin liquid to antiferromagnetic Mott insulator Our last subject is to investigate the quantum phase transition from the Z2 spin liquid to the antiferromagnetic Mott insulator. Here, we should take into account two order parameters such as the d ± id′ pairing and antiferromagnetic ones. Generically, we expect four possibilities. The first candidate is coexistence between such two orders, where the two critical lines cross each other. As a result, we have two critical points inside each phase. The second possibility is the multi-critical point, where the two critical points meet at one point. The third situation will be the first order transition between them. The last corresponds to an intermediate state without any ordering, where the two critical points do not meet. First of all, we exclude the last possibility because this phase is nothing but the algebraic spin liquid and there is no reason for this reentrant behavior. We start to examine the possibility of coexistence. The antiferromagnetic critical point inside the Z2 spin liquid phase can be determined by m = 0 while v and ve are finite, thus determined self-consistently. The mean field equations for this quantum critical point are given in appendix C-1. The strategy of solving the system of equations is how to reduce the number of self-consistent equations. Detailed calculations are provided in appendix C-1. As a result, we obtain two self-consistent equations for two unknown variables. These equations can be solved numerically. For the first (κc = 1, κs = 1) and third (κc = 3/2, κs = 1/2) decomposition schemes in appendix A, we could show that there are no mean field solutions at the transition point. On the other hand, we find Um /t = 0.360 in the case of the d + id′ pairing symmetry for the second decomposition scheme (κc = 1, κs = 1/2). The other quantum critical point is the Z2 spin liquid critical point inside the antiferromagnetic phase. It can be found when v = ve = 0 but m is finite, determined self-consistently. The mean field equations for this quantum critical point are given in appendix C-2. Solving the mean field equations self consistently, we could not find any solution. On the other hand, if the direct phase transition from the antiferromagnetic Mott insulator to the semi-metal is concerned, we find the critical point occurs at Um /t = 0.330 for the second decomposition (κc = 1, κs = 1/2). Our analysis for the quantum phase transition from the Z2 spin liquid to the antiferromagnetic Mott insulator shows that the nature of this transition depends on our phenomenological parameters of κc and κs . We could find the antiferromagnetic quantum critical point inside the Z2 spin liquid state for particular values of κc and κs while we could not obtain the Z2 spin liquid quantum critical point inside the antiferromagnetic Mott 8 insulating phase. We could not find the multi-critical point solution, either. As mentioned before, it is difficult to expect the algebraic spin liquid solution between the Z2 spin liquid and antiferromagnetic phases. Actually, we could find only one solution for the Z2 spin liquid to algebraic spin liquid transition, given by the previous subsection. The remaining possibility is the first order transition between the Z2 spin liquid and antiferromagnetic Mott insulator. We believe that the first order transition is the generic case for the phase transition between these two phases. The formal procedure will be to integrate over spinons and holons and to obtain an effective Landau-Ginzburg-Wilson free energy functional for both d ± id′ spin singlet pairing and antiferromagnetic order parameters. Based on the effective field theory, we can perform the renormalization group analysis and find the nature of the phase transition. This study is beyond the scope of the present study. One may ask the possibility of the Landau-GinzburgWilson forbidden continuous transition between the Z2 spin liquid with the d ± id′ pairing symmetry and the antiferromagnetic Mott insulator. Classification of LandauGinzburg-Wilson forbidden continuous transitions in two spatial dimensions has been performed in Ref. 30. Investigating the classification table carefully, we can find that this transition does not belong to any cases. The main reason is that the singlet pairing order parameter cannot be symmetrically equivalent to the antiferromagnetic order parameter. The classification scheme reveals that the Néel order parameter can form a hyper-vector with a triplet pairing order parameter. In this respect we are allowed to exclude the possibility of the Landau-GinzburgWilson forbidden continuous transition between the Z2 spin liquid and the antiferromagnetic Mott insulator. IV. DISCUSSION AND SUMMARY In this paper we investigated the phase structure of the Hubbard model on the honeycomb lattice. Physics of one dimensional interacting electrons is our reference. As well known, even if we start from weak interactions, they become enhanced at low energies, destabilizing the Fermi liquid state. In one dimension such quantum corrections can be summed exactly, resorting to the Ward identity.31 The resulting electron Green’s function shows two kinds SM ASL GSL AFM X Uc Uv Um U/t FIG. 3: Schematic phase diagram. Abbreviations: SM is the semi-metal phase, ASL is the algebraic spin liquid, GSL is the gapped spin liquid, and AFM is antiferromagnetism. The SMASL and ASL-GSL quantum phase transitions belong to the second order while the GSL-AFM quantum phase transition is the first order. of branch cuts, corresponding to collective charge and spin excitations. In this diagrammatic approach it is difficult to see the nature of such fractionalized excitations. But, the bosonization approach is helpful at low energies, revealing that spinons and holons are identified with topological solitons such as domain walls.3 One can interpret this phenomenon in another respect that topological solitons acquire fermion quantum numbers via fermion zero modes, regarded as realization of quantum anomaly.4 We believe that the spin-charge separation in one dimensional interacting electrons results from not only just interaction effects but also hidden topological properties of Dirac fermions. Then, the next natural question is whether we can find this physics in higher dimensions. The graphene structure is an ideal system for realization of Dirac fermions. The first observation in this Dirac fermion system is that the vanishing density of states needs a finite value of the interaction strength U for an antiferromagnetic order to be achieved. Then, the question is whether we can find intermediate phases between the semi-metal and antiferromagnetic Mott insulator, allowing fractionalized excitations as one dimensional interacting electrons. Indeed, we could find two kinds of paramagnetic Mott insulating phases, which show fractionalized excitations. The algebraic spin liquid appears from the semi-metal state via the Higgs transition, gapping of charge fluctuations. Although it is not clear how the topological nature of Dirac fermions is introduced to result in such a spin liquid state, spinon excitations in the algebraic spin liquid can be identified with topological excitations corresponding to meron (half skyrmion) excitations.32 The underlying mechanism is that the symmetry of the original microscopic model is enhanced at low energies, allowing a topological term to assign a fermion quantum number to such a topological excitation. The algebraic spin liquid turns out to have an O(5) symmetry in the physical case, where antiferromagentic correlations exhibit the same power-law dependence for distance as valence-bond correlations.14–16 It was pointed out that the corresponding effective field theory would be given by an O(5) WessZumino-Witten theory,16 identifying spinons with such topological excitations. Comparing the algebraic spin liquid with the Tomonaga-Luttinger liquid, there is one to one correspondence between them except that charge excitations are critical in the Tomonaga-Luttinger liquid. Actually, spin dynamics in one dimension is governed by the O(4) Wess-Zumino-Witten theory,3 describing critical dynamics of spinons. Because the stability of the algebraic spin liquid is not guaranteed beyond the large-Nf limit, we proposed how the quantum Monte Carlo simulation can prove the existence of such a phase. As discussed before, the symmetry enhancement can be verified, calculating both antiferromagnetic and valence-bond correlations at finite temperatures. If such correlations turn out to have the same scaling behavior, we have the algebraic spin liquid phase 9 just beside the semi-metal state. When interactions are increased more, pairing correlations between nearest neighbor sites become enhanced in the singlet channel, destabilizing the algebraic spin liquid. As a result, spinon excitations are gapped due to their pairing orders. An interesting point is that the nature of this gapped spin liquid state is given by the d + id′ singlet pairing order, which breaks time reversal symmetry. We would like to emphasize that time reversal symmetric combinations based on the d-wave pairing symmetry turn out to give higher energies than the d + id′ pairing order. We suspect that this time reversal symmetry breaking may be related with the Berry phase effect of the momentum space.33 One way to verify this statement is to check how the d + id′ pairing symmetry is changed, increasing the chemical potential from the Dirac point, where the Berry phase effect becomes weaken. Unfortunately, the quantum Monte Carlo simulation claimed that there is no time reversal symmetry breaking in the gapped spin liquid state. This inconsistency was resolved, taking d − id′ pairing to another valley. As a result, the edge state from the d + id′ pairing is cancelled by that from the d − id′ one. In addition to this time reversal symmetry, our proposal for the pairing order parameter turns out to be essential in order to have a fully gapped spectrum of spinon excitations, thus energetically more favorable than the case of only the d + id′ pairing, where spin excitations remain gapless. We suggested an experimental signature, that is, the quantized thermal valley Hall effect as the fingerprint of this gapped spin liquid. Finally, we investigated the quantum phase transition from the Z2 spin liquid to the antiferromagnetic Mott insulator. We concluded that the first order transition will take place generically. We argue that this first order transition is involved with two symmetrically unrelated order parameters, displaying different discrete symmetry properties, here time reversal symmetry. We claim that the Landau-Ginzburg-Wilson forbidden continuous transition will not appear, based on the existing classification scheme in the two dimensional Dirac theory on the honeycomb lattice.30 We would like to point out that the SU(2) slave-rotor theory seems to overestimate quantum fluctuations. If one sets κc as the order of 1, the critical strength of the Mott transition is the order of 10−1 for the critical value, compared with that from the quantum Monte Carlo simulation. This overestimation originates from strong band renormalization for spinons and holons, given by effective hopping integrals, Xij and Yij . Qualitatively the same situation also happens in the U(1) slave-rotor theory13 while the SU(2) slave-rotor theory seems to overestimate quantum fluctuations more. We believe that this aspect should be investigated more sincerely. Recently, the role of the spin-orbit interaction in the Hubbard model on the honeycomb lattice has been studied both extensively and intensively, where one purpose is to reveal the interplay between the topological band structure given by the spin-orbit coupling and strong correlation effect. Novel exotic phases have been suggested in this Kane-Mele-Hubbard model, some of which are quantum spin Hall effect in a transition metal oxide such as Na2 IrO3 34 , a spin liquid state with a topological band structure35 , and the chiral spin liquid state with the anyon nature of excitations36 . These interesting proposals will be verified based on ”exact” numerical calculations37 . Acknowledgments We would like to thank T. Takimoto for useful discussions. We were supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MEST) (No. 2010-0074542). M.-T. was also supported by the National Foundation for Science and Technology Development (NAFOSTED) of Vietnam. Appendix A: Decoupling scheme We discuss several decoupling schemes. The first example is UX † U X † 2HU = (ψi σx ψi )2 + (ψi σy ψi )2 6 i 6 i U XX UXX ( ( niσ − 1)2 + niσ − 1) + 6 i σ 6 i σ UX † U X − (σciσ ciσ )2 + niσ . (A1) 2 iσ 2 iσ Formally, this magnetic decoupling does not correspond to the conventional Hartree-Fock analysis for antiferromagnetism because the interaction strength is twice larger the standard mean field value. The second possible decoupling is UX † U X † 2HU = (ψi σx ψi )2 + (ψi σy ψi )2 6 i 6 i U XX UXX ( ( niσ − 1)2 + niσ − 1) + 6 i σ 6 i σ X † UX X † + [( ciσ ciσ )2 − ( σciσ ciσ )2 ]. (A2) 4 i σ σ This decoupling recovers the standard mean field theory for antiferromagnetism, but the coefficient of the term PP 5 U , not equal to the coefficient of the ( niσ − 1)2 is 12 σ i P † term (ψi σx ψi )2 . i The third possible decoupling is U X † U X † (ψi σx ψi )2 + (ψi σy ψi )2 2HU = 4 i 4 i X † U X X † + [( ciσ ciσ )2 − ( σciσ ciσ )2 ]. (A3) 4 i σ σ 10 The third decoupling scheme seems natural, but we in- troduce phenomenological parameters κc and κs . Appendix B: The algebraic spin liquid to Z2 spin liquid critical point The mean field equations for the algebraic spin liquid to Z2 spin liquid critical point are given by w ev X δ 1 X |γ(k)|, 4N/2 |γδ |2 = (B1) k 1 1 X |ζ(k)|2 1 X |ζδ |2 = , ϑv 4wv N/2 |γ(k)| δ X wv δ ϑv X δ 2 |γδ | = 2 |ζδ | = r 1= r   1 κc U v 1 X 1 −q , |γ(k)| q 6tw ec 2N/2 ev − γ(k) ev + γ(k) k λ λ X |ζ(k)|2 1 κc U v 6tw ev 2w ev N/2 γ(k) k r (B2) k  κc U v 1 X 6tw ev N/2 k  1 1 q −q ev − γ(k) ev + γ(k) λ λ 1 1 q +q ev − γ(k) ev + γ(k) λ λ  ,  , (B3) (B4) (B5) ev is redefined. where λv = 2tw ev λ The strategy for the critical interaction strength is to solve Eq. (B5) with w ev from Eq. (B1). The point is how to e find λv from other equations. From Eqs. (B3) and (B5) we obtain   P 1 1 1 √ √ |γ(k)| − 2N/2 ev −γ(k) ev +γ(k) λ λ k  . wv = (B6) P P 1 √ 1 +√ 1 |γδ |2 N/2 δ ev −γ(k) λ k ev +γ(k) λ Inserting Eq. (B6) into Eq. (B2), we get 2 ϑv = P δ P |ζδ | δ 1 2 |γδ | N/2 1 N/2 2 |γ(k)|  1 P √ N/2 P |ζ(k)|2 k |γ(k)| ϑv X δ |ζδ |2 = P |ζ(k)|2 k k γ(k) 1 2w ev N/2 P k    P δ 1 N/2 2 2 |ζδ | P |ζ(k)|2 = k |γ(k)|  P 2 2 |γδ | P 1 |γ(k)| N/2 δ k 1 ev −γ(k) λ √ 1 ev −γ(k) λ √ Equations (B7) and (B8) with w ev from Eq. (B1) give 1 N/2 1 N/2 √ 1 ev −γ(k) λ 1 ev −γ(k) λ k From Eqs. (B4) and (B5) we obtain 1 N/2  P P |ζ(k)|2 k P k γ(k) √ 1 ev +γ(k) λ − 1 λv +γ(k) + √e 1 λv +γ(k) − √e  .  1 ev −γ(k) λ − √e 1 λv −γ(k) − √e √  |γ(k)| √e .  + √e 1 λv +γ(k)   (B7) (B8) 1 λv +γ(k) 1 λv +γ(k)  . (B9) 11 ev . Once we find λ ev , we can obtain the critical value from Eq. (B5) together with Eq. This equation determines λ (B1), given by κc U v = t 1 3 N/2 2 P δ |γδ |2  1 N/2 P k  P k |γ(k)| √ 1 ev −γ(k) λ + √ 1 ev +γ(k) λ 2 . (B10) Appendix C: To analyze the quantum phase transition from the Z2 spin liquid to the antiferromagnetic Mott insulator 1. To find the antiferromagnetic quantum critical point inside the Z2 spin liquid state The mean field equations for this quantum critical point are given by w em vem X wm δ vm X δ 2 |γδ | = 2 |ζδ | = |γδ |2 = X |ζδ |2 = δ δ r r X 1 X wm |γ(k)|2 p , 2 |γ(k)|2 + v 2 |ζ(k)|2 4N/2 wm m vm |ζ(k)|2 1 X p , 2 |γ(k)|2 + v 2 |ζ(k)|2 4N/2 wm m k em X |γ(k)|2 κc U m w 3 2N/2 E(k) k κb Um vem X |ζ(k)|2 3 2N/2 E(k) k    1 1 p p , − λm − 2tE(k) λm + 2tE(k)  1 1 p , −p λm − 2tE(k) λm + 2tE(k) 1 κs U m 1 X p , 2 2 2 |ζ(k)|2 t N/2 wm |γ(k)| + vm 1= (C1) k (C2) (C3) (C4) (C5) k 1= r κc U m 1 X 3 N/2 k   1 1 p . +p λm − 2tE(k) λm + 2tE(k) Introducing xm = vm /wm and x em = vem /w em , we rewrite the above equations as w em wm X δ |γδ | 2 = r vem X |γδ |2 = X |ζδ |2 = δ δ |γ(k)|2 1 X p , 4N/2 |γ(k)|2 + x2m |ζ(k)|2 (C6) (C7) k xm |ζ(k)|2 1 X p , 4N/2 |γ(k)|2 + x2m |ζ(k)|2 (C8) k κc U m 1 X |γ(k)|2 p 6tw em 2N/2 |γ(k)|2 + x e2m |ζ(k)|2 k   1 1 , −q × q p p em − |γ(k)|2 + x em + |γ(k)|2 + x λ λ e2m |ζ(k)|2 e2m |ζ(k)|2 (C9) 12 vm X δ |ζδ |2 = r em X |ζ(k)|2 κc U m x p 6tw em 2N/2 |γ(k)|2 + x e2m |ζ(k)|2 k   1 1 q , − × q p p em − |γ(k)|2 + x em + |γ(k)|2 + x e2m |ζ(k)|2 e2m |ζ(k)|2 λ λ κs U m 1 X 1 p , twm N/2 |γ(k)|2 + x2m |ζ(k)|2 1= (C10) (C11) k 1= r κc U m 1 X 6tw em N/2 k  1 1 q +q p p em − |γ(k)|2 + x em + |γ(k)|2 + x λ λ e2m |ζ(k)|2 e2m |ζ(k)|2 em is redefined. From Eqs. (C7) and (C8) we get where λm = 2tw em λ P δ x em P δ |ζδ |2 |γδ |2 1 N/2 = xm 1 N/2 Similarly, Equations (C9) and (C10) give P δ xm P δ |ζδ |2 |γδ |2 1 N/2 =x em 1 N/2 P √ |ζ(k)|2 |γ(k)|2 +e x2m |ζ(k)|2 P √ |γ(k)|2 |γ(k)|2 +e x2m |ζ(k)|2 P √ k k   P k √ |ζ(k)|2 |γ(k)|2 +x2m |ζ(k)|2 k √ |γ(k)|2 |γ(k)|2 +x2m |ζ(k)|2 P  , (C12) . (C13) q √ 1 em − |γ(k)|2 +e λ x2m |ζ(k)|2 − q √ 1 em + |γ(k)|2 +e λ x2m |ζ(k)|2 q √ 1 2 |ζ(k)|2 em − |γ(k)|2 +e λ xm − q √ 1 2 |ζ(k)|2 em + |γ(k)|2 +e λ xm q √ 1 em − |γ(k)|2 +e λ x2m |ζ(k)|2 − q  . (C14) From Eqs. (C9) and (C12) we obtain wm X δ |γδ |2 = 1 2 1 N/2 k |γ(k)|2 |γ(k)|2 +e x2m |ζ(k)|2 1 N/2 P k   q √ 1 em − |γ(k)|2 +e λ x2m |ζ(k)|2 + q em + λ √ em + λ √ 1 |γ(k)|2 +e x2m |ζ(k)|2 1 |γ(k)|2 +e x2m |ζ(k)|2   . Taking both sides of Eq. (C9) to the square power with w em from Eq. (C7), we obtain " #2  P |γ(k)|2 1 1 1 q √ −q √ √ N/2 em − |γ(k)|2 +e em + |γ(k)|2 +e |γ(k)|2 +e x2m |ζ(k)|2 X λ λ x2m |ζ(k)|2 x2m |ζ(k)|2 k κc U m 2 wm . |γδ | = P |γ(k)|2 1 6twm √ δ N/2 k |γ(k)| +xm |ζ(k)| k λm − (C17) |γ(k)| +xm |ζ(k)| Equations (C15), (C16), and (C17) give us   1 κc 1 X |γ(k)|2 1 p q q − p p 3κs N/2 |γ(k)|2 + x e2m |ζ(k)|2 em − |γ(k)|2 + x em + |γ(k)|2 + x k e2m |ζ(k)|2 e2m |ζ(k)|2 λ λ 2 P P |γ(k)| 1 1 1 √ √ N/2 2 +x2 |ζ(k)|2 N/2 2 +x2 |ζ(k)|2 |γ(k)| |γ(k)| m m k .  k = P 1 1 1 q q + √ √ N/2 2 2 2 2 2 2 e e k (C16) |γ(k)|2 +x2m |ζ(k)|2 Inserting Um /wm from Eq. (C11) into this equation, we obtain !#2 " |γ(k)|2 1 1 1 P √ q q − √ √ N/2 |γ(k)|2 +e x2m |ζ(k)|2 em − |γ(k)|2 +e em + |γ(k)|2 +e X λ x2m |ζ(k)|2 λ x2m |ζ(k)|2 k κc 2 . wm |γδ | = P P 2 |γ(k)| 1 1 1 6κs √ √ δ N/2 2 2 2 2 2 N/2 2 k (C15) |γ(k)| +e xm |ζ(k)| λm + |γ(k)| +e xm |ζ(k)| (C18) 13 Inserting P δ x em = xm P δ |γδ |2 1 N/2 |ζδ |2 1 N/2 P |ζ(k)|2 |γ(k)|2 +x2m |ζ(k)|2 k √ k √ P (C19) |γ(k)|2 |γ(k)|2 +x2m |ζ(k)|2 from Eq. (C13) into Eqs. (C14) and (C18), we obtain two self-consistent equations for two unknown variables, xm em . These equations can be solved numerically. We fix xm first, and solve these two equations for λ em . Then, we and λ e obtain two functions, λm of xm . When two lines of these functions intersect, we obtain the solution of such equations. em are determined, the critical value of Um is also found from Eq. (C11). Once xm and λ 2. To find the Z2 spin liquid quantum critical point inside the antiferromagnetic Mott insulator The mean field equations at this critical point are given by X 1 X |γ(k)|2 q w ea |γδ |2 = , 4N/2 |γ(k)|2 + ( ma )2 δ k wa δ ϑa X δ twa X 1 |ζ(k)|2 1 X q |ζδ |2 = , ϑa 4wa N/2 |γ(k)|2 + ( ma )2 δ X (C20) 2 |γδ | = 2 |ζδ | = r r k twa   κc U a 1 X 1 1 q q − , |γ(k)| 6tw ea 2N/2 ea − γ(k) ea + γ(k) k λ λ X |ζ(k)|2 1 κc U a 6tw ea 2w ea N/2 γ(k) k 1=  r 1 1 q −q ea − γ(k) ea + γ(k) λ λ κs U a 1 X 1 q , twa N/2 2 |γ(k)| + ( ma )2 k 1= (C21) κc U a 1 X 6tw ea N/2 k   , (C22) (C23) (C24) twa  1 1 q +q , ea − γ(k) ea + γ(k) λ λ (C25) ea is redefined. We solve these equations basically in the same way as the previous case. First, where λa = 2tw ea λ we reduce the system of six equations into two equations of two unknown variables, and solve the two equations numerically. From Eqs. (C22) and (C25) we get   P 1 1 1 √ √ |γ(k)| − N/2 ea −γ(k) ea +γ(k) λ λ 1 k   . wa = P (C26) P 2 |γδ |2 1 1 1 √ √ + δ N/2 ea −γ(k) λ k ea +γ(k) λ Equations (C23) and (C25) together with w ea from Eq. (C20) give   P |ζ(k)|2 1 1 1 √ √ − e N/2 γ(k) ea −γ(k) λ λa +γ(k) 1 k   P ϑa = 2 P 2w ea |ζδ | 1 1 1 √ √ + δ N/2 P 2 |γδ | = Pδ |ζδ |2 δ 1 N/2 2 1 N/2 P k √ P |ζ(k)|2 k γ(k)  ea −γ(k) λ k 1 ea −γ(k) λ √ |γ(k)|2 1 ma 2 N/2 ) |γ(k)|2 +( tw a P k  ea +γ(k) λ 1 λa +γ(k) − √e 1 ea −γ(k) λ √  1 λa +γ(k) + √e . (C27) 14 From Eqs. (C21) and (C26) we obtain P 1 N/2 2 2 |ζδ | ϑa = Pδ |γδ |2 δ 1 N/2 Equations (C27) and (C28) leads to  P  δ P δ |γδ |2 |ζδ |2 2 1 N/2 2 = 1 N/2 P k P k k |γ(k)|  √ 1 ea −γ(k) λ 2 1 √ |ζ(k)| ma 2 N/2 ) |γ(k)|2 +( tw a |γ(k)| P |ζ(k)|2 k P γ(k)   1 ea −γ(k) λ √ P k  − √ 1 ea −γ(k) λ 1 λa +γ(k) − √e √ 1 ea +γ(k) λ √ 1 ea −γ(k) λ − P √ √ 1 ea +γ(k) λ   1 N/2 1 N/2 √ 1 ea +γ(k) λ + P k P k  . |γ(k)|2 ma 2 ) |γ(k)|2 +( tw a √ (C28) . (C29) |ζ(k)|2 ma 2 ) |γ(k)|2 +( tw a √ From Eqs. (C24) and (C25) we obtain 6κs w ea 1=  κc wa 1 N/2 1 N/2 P k  k 1 ma 2 ) |γ(k)|2 +( tw a 1 ea −γ(k) λ √ 1 λa +γ(k) + √e 2 . (C30) Inserting wa and w ea from Eqs. (C20) and (C26) into this equation, one obtains 3κs 1= κc 1 N/2 1 N/2 P k P k |γ(k)|2 ma 2 ) |γ(k)|2 +( tw a  |γ(k)| √e 1 λa −γ(k) 1 λa +γ(k) − √e Equations (C29) and (C31) are the last two equations ea and ma /twa . We fix ma /twa first, and determining λ ea numerically. Then, we solve the two equations for λ 1 2 3 4 5 1 N/2 √ X.-L. Qi and S.-C. Zhang, arXiv:1008.2026 (unpublished); M. Z. Hasan and C. L. Kane, arXiv:1002.3895 (unpublished); P. G. Silvestrov and E. G. Mishchenko, arXiv:0912.4658 (unpublished); S. Murakami, Prog. Theor. Phys. Suppl. 176, 279 (2008); M. Koenig, H. Buhmann, L. W. Molenkamp, T. L. Hughes, C.-X. Liu, X.-L. Qi, and S.C. Zhang, J. Phys. Soc. Jpn. 77, 031007 (2008). E. R. Mucciolo and C. H. Lewenkopf, J. Phys.: Condens. Matter 22, 273201 (2010); S. Das Sarma, S. Adam, E. H. Hwang, and E. Rossi, arXiv:1003.4731 (unpublished); D. S. L. Abergel, V. Apalkov, J. Berashevich, K. Ziegler, and T. Chakraborty, Advances in Physics, 59 (4), 261 (2010); A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009). A. O. Gogolin, A. A. Nersesyan, and A. M. Tsvelik, Bosonization and Strongly Correlated Systems (Cambridge University Press, New York, 2004). R. Rajaraman, Solitons and Instantons(Elsevier Science, New York, 2003). Z. Y. Meng, T. C. Lang, S. Wessel, F. F. Assaad, and A. Muramatsu, Nature 464, 847 (2010).  1 N/2 P k  P k √ √ 1 ma 2 ) |γ(k)|2 +( tw a 1 ea −γ(k) λ 1 λa +γ(k) + √e . (C31) ea of ma /twa . When two lines of obtain two functions, λ e λa and ma /twa intersect, we find the solution of these equations. 6 7 8 9 10 11 12 13 14 15 16 17 18 F. Wang, Phys. Rev. B 82, 024419 (2010). Y.-M. Lu and Y. Ran, arXiv:1005.4229 (unpublished); Y.M. Lu and Y. Ran, arXiv:1007.3266 (unpublished). B. K. Clark, D. A. Abanin, and S. L. Sondhi, arXiv:1010.3011 (unpublished). G. Wang, M. O. Goerbig, B. Gremaud, and C. Miniatura, arXiv:1006.4456 (unpublished). A. Vaezi and X.-G. Wen, arXiv:1010.5744 (unpublished). Ki-Seok Kim, Phys. Rev. Lett. 97, 136402 (2006); Ki-Seok Kim, Phys. Rev. B 75, 245105 (2007); Ki-Seok Kim and Mun Dae Kim, Phys. Rev. B. 81, 075121 (2010). P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod. Phys. 78, 17 (2006). M. Hermele, Phys. Rev. B 76, 035125 (2007), and references therein. M. Hermele, T. Senthil, and M. P. A. Fisher, Phys. Rev. B 72, 104404 (2005). Y. Ran and X.-G. Wen, arXiv:cond-mat/0609620v3 (unpublished). A. Tanaka and X. Hu, Phys. Rev. Lett. 95, 036402 (2005). C. Honerkamp, Phys. Rev. Lett. 100, 146404 (2008). A. M. Black-Schaffer and S. Doniach, Phys. Rev. B 75, 15 19 20 21 22 23 24 25 26 27 134512 (2007). O. Vafek, A. Melikyan, and Z. Tesanovic, Phys. Rev. B 64, 224508 (2001). D. Xiao, W. Yao, and Q. Niu, Phys. Rev. Lett. 99, 236809 (2007); W. Yao, D. Xiao, and Q. Niu, Phys. Rev. B 77, 235406 (2008). S. Florens and A. Georges, Phys. Rev. B 70, 035114 (2004). H. Mosadeq, F. Shahbazi, and S. A. Jafari, arXiv:1007.0127 (unpublished). A. H. MacDonald, S. M. Girvin, and D. Yoshioka, Phys. Rev. B 41, 2565 (1990). R. Kumar, D. Kumar, and B. Kumar, Phys. Rev. B 80, 214428 (2009). A. Banerjee, K. Damle, and A. Paramekanti, arXiv:1012.4546 (unpublished). M. Hermele, T. Senthil, M. P. A. Fisher, P. A. Lee, N. Nagaosa, and X.-G. Wen, Phys. Rev. B 70, 214437 (2004). Unfortunately, we could not find any references on numerically ”exact” calculations for thermal Hall effect. Instead, it has been performed to evaluate electrical Hall effect based on exact simulations28,29 . 28 29 30 31 32 33 34 35 36 37 F. F. Assaad and M. Imada, Phys. Rev. Lett. 74, 3868 (1995). J. O. Haerter and B. S. Shastry, Phys. Rev. B 77, 045127 (2008). S. Ryu, C. Mudry, C.-Y. Hou, and C. Chamon, Phys. Rev. B 80, 205319 (2009). D. L. Maslov, arXiv:cond-mat/0506035 (unpublished). T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and M. P. A. Fisher, Science 303, 1490 (2004); T. Senthil, L. Balents, S. Sachdev, A. Vishwanath, and M. P.A. Fisher, Phys. Rev. B 70, 144407 (2004). S. Matsuura and S. Ryu, arXiv:1007.2200 (unpublished). A. Shitade, H. Katsura, J. Kunes, X.-L. Qi, S.-C. Zhang, and N. Nagaosa, Phys. Rev. Lett. 102, 256403 (2009). S. Rachel and K. L. Hur, Phys. Rev. B 82, 075106 (2010). J. He, S.-P. Kou, Y. Liang, and S. Feng, arXiv:1012.0620 (unpublished). M. Hohenadler, T. C. Lang, F. F. Assaad, arXiv:1011.5063 (unpublished).