Spin liquids in graphene

# Spin liquids in graphene

Minh-Tien Tran and Ki-Seok Kim 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
Department of Physics, Pohang University of Science and Technology, Pohang, Gyeongbuk 790-784, Korea
###### Abstract

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 Z 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 Z spin liquid turns out to be singlet pairing, but time reversal symmetry is preserved, taking in one valley and 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 Z spin liquid are shown to be continuous while the transition from the Z spin liquid to the antiferromagnetic Mott insulator turns out to be the first order. We emphasize that both algebraic spin liquid and Z 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 Z spin liquid.

## I Introduction

Interplay between the topological band structure and interaction drives one direction of research in modern condensed matter physics,Review_TI (); Review_Graphene () 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.1D_Book () 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.Soliton_TextBook ()

A recent study based on the quantum Monte Carlo simulation QMC_Graphene () has argued that essentially the same phenomenon as electron fractionalization in the Tomonaga-Luttinger 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, PSG (); Ran (); Sondhi () thus identified with a Z spin liquid. We point out other scenarios U1SR (); Wen () 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, PSG (); Ran (); Sondhi () 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, Kim_SR () is at the heart of the methodology, where exchange 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 theory Lee_Wen_Nagaosa () for the t-J model.

We find two kinds of spin liquids, identified with an algebraic spin liquid and a Z spin liquid, respectively, between the semi-metal and antiferromagnetic phases. We argue that the algebraic spin liquid Hermele1 (); Hermele2 () 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.Wen_Symmetry (); Tanaka_SO5 () 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 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.SC1_Graphene (); SC2_Graphene () However, we argue that time reversal symmetry is preserved, taking the singlet pairing to one valley while the pairing to another. This assignment turns out to be essential in order to have a fully gapped spectrum because the singlet pairing order parameter in one valley vanishes in the other valley, allowing the gapless Dirac spectrum. We propose the quantized thermal valley Hall effect QTHE (); QVHE () for the fingerprint of this gapped Z 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 Z 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 dependence between antiferromagnetic and valence bond correlations. The hallmark of the Z 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

### ii.1 Formulation

We start from the Hubbard model on the honeycomb lattice

 H=−t∑⟨ij⟩σc†iσcjσ+H.c.+U∑ini↑ni↓, (1)

where () is the annihilation (creation) operator for an electron at site with spin . is the hopping integral, and is the on-site Coulomb interaction, where represents the density of electrons with spin .

Introducing the Nambu-spinor representation

 ψi=⎛⎝ci↑c†i↓⎞⎠,

and performing the Hubbard-Stratonovich transformation for the pairing, density (singlet) and magnetic (triplet) interaction channels, we obtain an effective Lagrangian

 L=∑iψ†i(∂τ1−μσz)ψi−t∑⟨ij⟩ψ†iσzψj+H.c. −i∑i[ΦRi(ψ†iσxψi)+ΦIi(ψ†iσyψi)+φi(ψ†iσzψi)] +32Uκc∑i[(ΦRi)2+(ΦIi)2+(φi)2] +12Uκs∑im2i−∑imi(ψ†iψi). (2)

Here, and are associated with pairing-fluctuation and density-excitation potentials, introduced to decouple the charge channel. is an effective magnetic field, which decouples the spin channel. and 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 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 Z spin liquid phase transitions are shown not to depend on such phenomenological parameters, where both the critical points are determined with the combination between and . Only the Z 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 representation Kim_SR () means to write down an electron field as a composite field in terms of a charge-neutral spinon field and a spinless holon field

 ψi=Z†iFi, (3)

where is a fermion operator in the Nambu representation, and is an SU(2) matrix

 Zi=⎛⎜⎝zi↑−z†i↓zi↓z†i↑⎞⎟⎠. (4)

Here, is a boson operator, satisfying the unimodular (rotor) constraint, .

The key point of the slave-rotor representation SR () 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 () 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.SR () 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, (5) L0 = t∑⟨ij⟩Tr[XijY†ij+YijX†ij]+12Uκs∑im2i, (6) LF = ∑iF†i(∂τ1−iΩi⋅σ)Fi (7) −t∑⟨ij⟩(F†iXijFj+H.c.)−∑imi(F†iFi), LZ = 34Uκc∑iTr[Ωi⋅σ−iZi∂τZ†i+iμZiσzZ†i]2 (8) −t∑⟨ij⟩Tr[ZiσzZ†jY†ij+H.c.].

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 and , and shifting as

 Ωi⋅σ+iZi∂τZ†i−iμZiσzZ†i,

where is the pseudospin potential field, we recover Eq. (2) exactly with an introduction of an electron field . This procedure is well described in the previous study.Kim_SR () 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 which results from on-site pairing (virtual) fluctuations, captured by the off diagonal variable of the SU(2) matrix field . We note that the diagonal rotor field corresponds to the zero sound mode, giving rise to the Mott transition via gapping of their fluctuations. The additional boson rotor variable 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 theory Lee_Wen_Nagaosa () of the t-J model.

### ii.2 Mean-field ansatz

We perform the mean-field analysis, taking the following ansatz

 Xij=(wδv∗δvδ−w∗δ)⋅σz, (9) Yij=(˜wδ˜v∗δ˜vδ−˜w∗δ)⋅σz, (10)

where denotes the bond between the nearest neighbor sites. In the honeycomb lattice there are three nearest neighbor bonds. We choose , , , and , where and are symmetric factors for the hopping parameter () and the pairing order parameter (). The choice for and depends on the symmetry of the considered phase. For example, the -wave pairing symmetry is given by , the -wave symmetry , and the -wave symmetry . For the magnetic order parameter we choose an antiferromagnetic ansatz .

Particle-hole symmetry at half filling results in while pairing potentials of and vanish in the mean-field level. Then, we obtain a general expression for the free energy

 F=−1β∑k,iωlog[(iω)2−t2w2|γ(k)|2−(tv|ζ(k)|+m)2]−1β∑k,iωlog[(iω)2−t2w2|γ(k)|2−(tv|ζ(k)|−m)2] +2β∑k,iνlog[(−34κcU(iν)2+λ)2−4t2˜w2|γ(k)|2−4t2˜v2|ζ(k)|2]+4tN∑δ(w˜w|γδ|2+v˜v|ζδ|2)+N2κsUm2−Nλ, (11)

where is the energy dispersion for spinons and holons, and is associated with the pairing potential. is performed in the unit cell. is a Lagrange multiplier field, introduced to keep the slave-rotor constraint. is the total number of sites.

Minimizing the free energy, we obtain fully self-consistent equations for order parameters

 ˜w∑δ|γδ|2=−2tw4Nβ∑k,iω[|γ(k)|2(iω)2−t2w2|γ(k)|2−(tv|ζ(k)|+m)2+|γ(k)|2(iω)2−t2w2|γ(k)|2−(tv|ζ(k)|−m)2], (12)
 ˜v∑δ|ζδ|2=−24Nβ∑k,iω[(tv|ζ(k)|+m)|ζ(k)|(iω)2−t2w2|γ(k)|2−(tv|ζ(k)|+m)2+(tv|ζ(k)|−m)|ζ(k)|(iω)2−t2w2|γ(k)|2−(tv|ζ(k)|−m)2], (13)
 w∑δ|γδ|2=4t˜wN∑k|γ(k)|2(−34κcU(iν)2+λ)2−4t2˜w2|γ(k)|2−4t2˜v2|ζ(k)|2, (14)
 v∑δ|ζδ|2=4t˜vN∑k|ζ(k)|2(−34κcU(iν)2+λ)2−4t2˜w2|γ(k)|2−4t2˜v2|ζ(k)|2, (15)
 (16)
 1=4Nβ∑k,iν(−34κcU(iν)2+λ)(−34κcU(iν)2+λ)2−4t2˜w2|γ(k)|2−4t2˜v2|ζ(k)|2. (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 self-consistent mean-field equations at zero temperature

 ˜w∑δ|γδ|2=w8N/2∑k|γ(k)|2D(k,m)+w8N/2∑k|γ(k)|2D(k,−m), (18)
 (19)
 w∑δ|γδ|2=√κbU3˜w2N/2∑k|γ(k)|2E(k)(1√λ−2tE(k)−1√λ+2tE(k)), (20)
 v∑δ|ζδ|2=√κbU3˜v2N/2∑k|ζ(k)|2E(k)(1√λ−2tE(k)−1√λ+2tE(k)), (21)
 m=κsU2N/2∑kv|ζ(k)|+mtD(k,m)−κsU2N/2∑kv|ζ(k)|−mtD(k,−m), (22)
 1=√κbU31N/2∑k(1√λ−2tE(k)+1√λ+2tE(k)), (23)

where

 E(k)=√˜w2|γ(k)|2+˜v2|ζ(k)|2, (24) D(k,m)=√w2|γ(k)|2+(v|ζ(k)|+mt)2 (25)

are holon and spinon energy spectra in the presence of pairing and antiferromagnetism, respectively.

Considering symmetry, it is natural to take into account spatially uniform hopping

 |γ(k)|2=3+2cos(ky)+4cos(12ky)cos(√32kx). (26)

On the other hand, the -wave pairing potential is not allowed due to repulsive interactions. Counting the lattice symmetry of the honeycomb structure, the next candidate will be or for nearest neighbor singlet pairing SC2_Graphene (). We introduce a general combination of - and -wave pairing for the pairing term

 |ζ(k)|2=|cos(θ)ζx2−y2(k)+isin(θ)ζxy(k)|2, (27)

where is a combination factor, and () is the () -wave symmetry function

 ζx2−y2(kx,ky) = e−ikx√3−eikx2√3cos(ky2) ζxy(kx,ky) = ieikx2√3sin(ky2).

For this pairing symmetry becomes . We also consider the -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 and parameters, he can make this pairing symmetry a solution. However, this solution does not give the lowest free energy, compared with the pairing solution, consistent with earlier studies.SC1_Graphene (); SC2_Graphene ()

One may criticize the ansatz for uniform hopping in this paper because such an assumption excludes possible dimerized phases a priori. Actually, the Heisenberg model

 H=J1∑⟨ij⟩Si⋅Sj+J2∑⟨⟨kl⟩⟩Sk⋅Sl

has shown several types of dimerized phases when the ratio of is beyond a certain critical value Sondhi (); J1J2 (), approximately given by . 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 limit Hubbard_J1J2 (), where each parameter is given by Sondhi ()

 J1=4t{tU−4(tU)3},     J2=4t(tU)3

up to the fourth order process. Then, the ratio can be expressed in terms of as follows

 J2J1=1(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 term effectively, while the ring exchange term is expected to be small.Sondhi () However, the role of the ring exchange term has been also studied carefully.Ring_Exchange1 (); Ring_Exchange2 ()

An antiferromagnetic phase has been reported in .Sondhi (); J1J2 () This corresponds to , consistent with the result of the quantum Monte Carlo simulation.QMC_Graphene () Increasing frustration, the antiferromagnetic order disappears, and a paramagnetic Mott insulating state results, identified with a certain type of Z spin liquids. Such a spin-gapped state turns out to evolve into a dimerized phase with either translational or rotational symmetry breaking near .Sondhi (); J1J2 () It was reported that the spin liquid state turns into a dimerized phase with three-fold degeneracy around , where it breaks the symmetry but preserving the translational symmetry.Sondhi () On the other hand, the plaquette order was claimed to appear near before the dimerized phase, breaking the translational symmetry only.J1J2 () An important point is that if we translate the critical value in terms of of the Hubbard model, corresponds to and , . Comparing these critical values with the critical value for the semi-metal to spin liquid transition in the quantum Monte Carlo simulation QMC_Graphene (), the Mott critical value given by 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 Heisenberg model. In other words, the model seems to be an effective low energy model only in the limit of while physics of such a model will be different from that of the Hubbard model in the small 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 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.1 From semi-metal to algebraic spin liquid

The semi-metal phase is described by condensation of holons with . Considering the symmetry factor , the condensation occurs when the effective chemical potential given by the Lagrange multiplier field touches the maximum point of the holon dispersion, i.e., . These collective charge excitations become gapped when , and a Mott insulating state appears, characterized by with .

Taking with , we can determine the quantum critical point from the following mean-field equations

 ˜wc∑δ|γδ|2=14N/2∑k|γ(k)|, (29)
 wc∑δ|γδ|2=√κcUc6t˜wc12N/2∑k|γ(k)|(1√3−γ(k)−1√3+γ(k)), (30)
 1=√κcUc6t˜wc1N/2∑k(1√3−γ(k)+1√3+γ(k)). (31)

Inserting from Eq. (29) into Eq. (31), one obtains the critical value for the interaction strength

 κcUct=321N/2∑k|γ(k)|∑δ|γδ|2[1N/2∑k(1√3−γ(k)+1√3+γ(k))]2 =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.Hermele1 ()

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- limit,Hermele1 () where 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.Deconfinement_ASL () 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.Hermele2 (); Wen_Symmetry () 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.Tanaka_SO5 () 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 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.

### iii.2 From algebraic spin liquid to Z2 spin liquid

Increasing more from the semi-metal to algebraic spin liquid critical point , we find another paramagnetic Mott insulating phase, characterized by and with the pairing symmetry. Recall Eq. (27) and Eq. (28) for pairing symmetries that we checked explicitly. The algebraic spin liquid () to gapped spin liquid ( and with ) critical point is found with an ansatz of but . The 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 pairing symmetry. Actually, we checked self-consistency for various values of the angle parameter in Eqs. (27) and (28), and found in Eq. (27), corresponding to . In Fig. 2 we plot the spinon spectrum given by Eq. (25) for the pairing symmetry. It shows that if only the pairing order parameter is taken into account, the energy spectrum opens a gap at one Brillouin zone edge (for instance, the point), but it still keeps the Dirac cone at the other inequivalent Brillouin zone edge (the point). On the other hand, if we consider only the pairing symmetry, we see that the point remains gapless while only the point becomes gapped. This demonstration motivates us to assign the pairing symmetry to one valley () and the to the other (), 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 pairing in one valley is cancelled by that from the 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 . Note that . This intermediate phase between the semi-metal 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 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.QTHE (); QVHE () 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 Z spin liquid phase. This may be verified Hall_Numerics () by either quantum Monte Carlo simulation Hall_tU () or exact diagonalization Hall_tJ ().

It should be noted that our time reversal symmetry preserving Z spin liquid state is beyond the classification scheme based on the projective symmetry group because their possible Z spin liquids in the projective symmetry group are constrained with complete time reversal symmetric pairing.PSG (); Ran () In other words, 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.SC1_Graphene (); SC2_Graphene ()

### iii.3 From Z2 spin liquid to antiferromagnetic Mott insulator

Our last subject is to investigate the quantum phase transition from the Z spin liquid to the antiferromagnetic Mott insulator. Here, we should take into account two order parameters such as the 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 Z spin liquid phase can be determined by while and 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 (, ) and third (, ) 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 in the case of the pairing symmetry for the second decomposition scheme (, ).

The other quantum critical point is the Z spin liquid critical point inside the antiferromagnetic phase. It can be found when but 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 for the second decomposition (, ).

Our analysis for the quantum phase transition from the Z spin liquid to the antiferromagnetic Mott insulator shows that the nature of this transition depends on our phenomenological parameters of and . We could find the antiferromagnetic quantum critical point inside the Z spin liquid state for particular values of and while we could not obtain the Z spin liquid quantum critical point inside the antiferromagnetic Mott 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 Z spin liquid and antiferromagnetic phases. Actually, we could find only one solution for the Z spin liquid to algebraic spin liquid transition, given by the previous subsection. The remaining possibility is the first order transition between the Z 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 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-Ginzburg-Wilson forbidden continuous transition between the Z spin liquid with the pairing symmetry and the antiferromagnetic Mott insulator. Classification of Landau-Ginzburg-Wilson forbidden continuous transitions in two spatial dimensions has been performed in Ref. Ryu_Dirac_Mass, . 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-Ginzburg-Wilson forbidden continuous transition between the Z 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.Maslov_1D () The resulting electron Green’s function shows two kinds 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.1D_Book () 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.Soliton_TextBook () 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 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.Senthil_DQCP () 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.Hermele2 (); Wen_Symmetry (); Tanaka_SO5 () It was pointed out that the corresponding effective field theory would be given by an O(5) Wess-Zumino-Witten theory,Tanaka_SO5 () 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,1D_Book () describing critical dynamics of spinons.

Because the stability of the algebraic spin liquid is not guaranteed beyond the large- 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 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 singlet pairing order, which breaks time reversal symmetry. We would like to emphasize that time reversal symmetric combinations based on the -wave pairing symmetry turn out to give higher energies than the pairing order. We suspect that this time reversal symmetry breaking may be related with the Berry phase effect of the momentum space.Ryu_Berry () One way to verify this statement is to check how the 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 pairing to another valley. As a result, the edge state from the pairing is cancelled by that from the 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 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 Z 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.Ryu_Dirac_Mass ()

We would like to point out that the SU(2) slave-rotor theory seems to overestimate quantum fluctuations. If one sets as the order of , the critical strength of the Mott transition is the order of 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, and . Qualitatively the same situation also happens in the U(1) slave-rotor theory Hermele1 () 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 NaIrO Na2IrO3 (), a spin liquid state with a topological band structure KLe_Hur (), and the chiral spin liquid state with the anyon nature of excitations Kou_CSL (). These interesting proposals will be verified based on ”exact” numerical calculations Assaad_KMH ().

###### Acknowledgements.
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

 2HU = U6∑i(ψ†iσxψi)2+U6∑i(ψ†iσyψi)2 (33) + U6∑i(∑σniσ−1)2+U6∑i(∑σniσ−1) − U2∑iσ(σc†iσciσ)2+U2∑iσniσ.

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

 2HU = U6∑i(ψ†iσxψi)2+U6∑i(ψ†iσyψi)2 (34) + U6∑i(∑σniσ−1)2+U6∑i(∑σniσ−1) + U4∑i[(∑σc†iσciσ)2−(∑σσc†iσciσ)2].

This decoupling recovers the standard mean field theory for antiferromagnetism, but the coefficient of the term is , not equal to the coefficient of the term .

The third possible decoupling is

 2HU = U4∑i(ψ†iσxψi)2+U4∑i(ψ†iσyψi)2 (35) + U4∑i[(∑σc†iσciσ)2−(∑σσc†iσciσ)2].

The third decoupling scheme seems natural, but we introduce phenomenological parameters and .

## Appendix B The algebraic spin liquid to Z2 spin liquid critical point

The mean field equations for the algebraic spin liquid to Z spin liquid critical point are given by

 ˜wv∑δ|γδ|2=14N/2∑k|γ(k)|, (36)
 1ϑv∑δ|ζδ|2=14wv1N/2∑k|ζ(k)|2|γ(k)|, (37)
 (38)
 ϑv∑δ|ζδ|2=√κcUv6t˜wv12˜wvN/2∑k|ζ(k)|2γ(k)(1√˜λv−γ(k)−1√˜λv+γ(k)), (39)
 (40)

where is redefined.

The strategy for the critical interaction strength is to solve Eq. (40) with from Eq. (36). The point is how to find from other equations. From Eqs. (38) and (40) we obtain

 wv=12N/2∑k|γ(k)|(1√˜λv−γ(k)−1√˜λv+γ(k))∑δ|γδ|21N/2∑k(1√˜λv−γ(k)+1√˜λv+γ(k)). (41)

Inserting Eq. (41) into Eq. (37), we get

 ϑv = (42)

From Eqs. (39) and (40) we obtain

 (43)

Equations (42) and (43) with from Eq. (36) give

 [∑δ|ζδ|2]21N/2∑k|ζ(k)|2|γ(k)| = [∑δ|γδ|2]21N/2∑k|γ(k)|1N/2∑k|ζ(k)|2γ(k)(1√˜λv−γ(k)−1√˜λv+γ(k))1N/2∑k|γ(k)|(1√˜λv−γ(k)−1√˜λv+γ(k)). (44)

This equation determines . Once we find , we can obtain the critical value from Eq. (40) together with Eq. (36), given by

 κcUvt=31N/2∑k|γ(k)|2∑δ|γδ|2[1N/2∑k(1√˜λv−γ(k)+1√˜λv+γ(k))]2. (45)

## Appendix C To analyze the quantum phase transition from the Z2 spin liquid to the antiferromagnetic Mott insulator

### c.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

 ˜wm∑δ|γδ|2=14N/2∑kwm|γ(k)|2√w2m|γ(k)|2+v2m|ζ(k)|2, (46)
 ˜vm∑