Grassmann tensor renormalization group for one-flavor lattice Gross-Neveu model with finite chemical potential
aa r X i v : . [ h e p - l a t ] J a n Preprint number: KANAZAWA-14-10
Grassmann tensor renormalization group forone-flavor lattice Gross-Neveu model withfinite chemical potential
Shinji Takeda ∗ and Yusuke Yoshimura ∗ Institute for Theoretical Physics, Kanazawa University, Kanazawa 920-1192, Japan ∗ E-mail: [email protected]; [email protected] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
We apply the Grassmann tensor renormalization group (GTRG) to the one-flavor latticeGross-Neveu model in the presence of chemical potential. We compute the fermionnumber density and its susceptibility and confirm the validity of GTRG for the finitedensity system. We introduce a method analogous to the reweighting method for MonteCarlo method and test it for some parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Subject Index typeset using PTP
TEX.cls . Introduction
The tensor renormalization group (TRG) is one of the numerical renormalization methods.It was originally introduced by Levin and Nave in the triangular lattice Ising model [1], andthen has been applied to bosonic models on square lattice: X-Y model[2], O (3) model[3] and φ theory[4]. Recently, Xie et al. developed the second TRG[5] to improve this method byusing global optimization instead of local one. An extension to higher dimensional system,named higher order TRG, was also introduced by Xie et al. and it was examined in the 3DIsing model[6]. Furthermore, a generalization to fermion system called Grassmann tensorrenormalization group (GTRG) was proposed by Gu et al. [7, 8]. Then it has been appliedto the two-dimensional QED[9], which is a gauge-fermion system, and a study including θ -term [10] was also given by Shimizu and Kuramashi.An advantage of the TRG is that this method can be applied to any systems suffering fromthe sign problem with the Monte Carlo method. The sign problem occurs in for example,finite fermion density systems, θ -term included systems, lattice chiral gauge theories and soon. The purpose of this paper is to apply the GTRG to a simple finite fermion density systemon the lattice, namely the Gross-Neveu model [11] containing four-fermion interaction in thepresence of chemical potential with the Wilson fermion lattice formulation. This modelis known to share important properties, asymptotically free and spontaneous symmetrybreaking, with QCD and considered to be its toy model. Usually, large N -expansion is usedto analyze the model. In this paper, nevertheless, we restrict to one flavor just for a simplicityalthough generalization to many flavors is straightforward. This work provides a benchmarkfor future study of complicated and higher-dimensional model, say QCD with finite densityeventually.This paper is organized as follows. In Sec. 2, after defining the lattice Gross-Neveu modelwith chemical potential, we review the derivation of the tensor network representation forthe model and explain the GTRG procedure as well as how to implement the anti-periodicboundary condition in this representation. Numerical results are presented in Sec. 3. In Sec.4, we propose a method analog to the reweighting method in the Monte Carlo method. Sec.5 is devoted to summary and outlook.
2. TRG for the Lattice Gross-Neveu Model
The Lagrangian density for the Gross-Neveu model [11] in two-dimensional Euclideancontinuum space-time is given by L GN = ¯ ψ ( ∂ + m ) ψ − g N h(cid:0) ¯ ψψ (cid:1) + (cid:0) ¯ ψiγ ψ (cid:1) i , (1)where ψ = ( ψ , ψ ) T is 2-component spinor field with N different flavors and m , g denotethe mass and the coupling constant respectively.The lattice version of the Lagrangian density is defined by L GNLat = X n ′ ¯ ψ n D n,n ′ ψ n ′ + g N h(cid:0) ¯ ψ n ψ n (cid:1) + (cid:0) ¯ ψ n iγ ψ n (cid:1) i , (2)2here n = ( n , n ) is a lattice site. The Wilson-Dirac operator[12, 13, 14] D n,n ′ includingthe chemical potential[15] is explicitly given by D n,n ′ = ( m + 2) δ n,n ′ − X ν, ± e ∓ µδ ν, (1 ± γ ν ) δ n,n ′ ± ˆ ν . (3)In the following, we consider 2-dimensional lattice box N × N . Let us first express the partition function in terms of the tensor network representation. Ageneral procedure to derive the tensor network representation is 1): expanding the integrand(Boltzmann weight) and then accompanying discrete variables describing the ordering ofexpansion become new degree of freedom (index of tensor), 2) integrating out the originaldegree of freedom ( ψ in this case) and then the elements of tensor are determined. Althoughthe derivation was already given in the previous work[9], we re-derive it in a slightly differentway to make this paper self-contained and hope that this is useful for readers. In this andnext subsections, we temporarily consider a system where periodic boundary condition isimposed for all directions, while the anti-periodic boundary condition will be discussed inSec. 2.4.In the following, we restrict to N = 1 and choose the representation of the gamma matrices γ = σ = ! , γ = σ = − ! , γ = iγ γ = σ = − ii ! . (4)The hopping terms for 2-direction is diagonal in spinor space,12 e − µ ¯ ψ n (1 + γ ) ψ n − ˆ2 = e − µ ¯ ψ n, ψ n − ˆ2 , , (5)12 e µ ¯ ψ n (1 − γ ) ψ n +ˆ2 = e µ ¯ ψ n, ψ n +ˆ2 , , (6)while the hopping terms for 1-direction¯ ψ n (1 ± γ ) ψ n ∓ ˆ1 = ¯ ψ n ± ± ! ψ n ∓ ˆ1 , are not diagonalized. One, however, can make them diagonal by introducing another basis: χ n, = 1 √ ψ n, + ψ n, ) , χ n, = 1 √ ψ n, − ψ n, ) , (7)¯ χ n, = 1 √ (cid:0) ¯ ψ n, + ¯ ψ n, (cid:1) , ¯ χ n, = 1 √ (cid:0) ¯ ψ n, − ¯ ψ n, (cid:1) , (8)which yields 12 ¯ ψ n (1 + γ ) ψ n − ˆ1 = ¯ χ n, χ n − ˆ1 , , (9)12 ¯ ψ n (1 − γ ) ψ n +ˆ1 = ¯ χ n, χ n +ˆ1 , . (10)(11) As discussed in Ref.[13, 14], the coupling constants for each four-fermion interaction ( ¯ ψψ ) and( ¯ ψiγ ψ ) should be treated independently for Wilson fermion formulation but as we will consideronly one-flavor theory where any kind of four-fermion interaction terms give a unique form, we donot distinguish between them in this paper. { χ n,i , ¯ χ n ′ ,j } = { χ n,i , χ n ′ ,j } = { ¯ χ n,i , ¯ χ n ′ ,j } = 0 , ∀ n, n ′ , i, j. (12)After the change of variable for the 1-direction hopping term, the Lagrangian density iswritten by L GNLat = ( m + 2) ¯ ψ n ψ n − g ¯ ψ n, ψ n, ¯ ψ n, ψ n, − ¯ χ n, χ n − ˆ1 , − ¯ χ n, χ n +ˆ1 , − e − µ ¯ ψ n, ψ n − ˆ2 , − e µ ¯ ψ n, ψ n +ˆ2 , , (13)and then corresponding the partition function with periodic boundary condition is given by Z P = Z DψD ¯ ψ exp − X n L GNLat ! = Z DψD ¯ ψ Y n e − ( m +2) ¯ ψ n, ψ n, e − ( m +2) ¯ ψ n, ψ n, e g ¯ ψ n, ψ n, ¯ ψ n, ψ n, · e ¯ χ n +ˆ1 , χ n, e ¯ χ n, χ n +ˆ1 , e e − µ ¯ ψ n +ˆ2 , ψ n, e e µ ¯ ψ n, ψ n +ˆ2 , , (14)where we have shown the spinor components explicitly. When expanding each exponentialfactor in the integrand, terms of order 0 and 1 give non-vanishing contribution due toGrassmann nature, Z P = X { s,t,x =0 , } Z DψD ¯ ψ Y n (cid:2) − ( m + 2) ¯ ψ n, ψ n, (cid:3) s n, (cid:2) − ( m + 2) ¯ ψ n, ψ n, (cid:3) s n, (cid:0) g ¯ ψ n, ψ n, ¯ ψ n, ψ n, (cid:1) s n, · (cid:16) ¯ χ n +ˆ1 , χ n, (cid:17) x n, (cid:16) ¯ χ n, χ n +ˆ1 , (cid:17) x n, (cid:16) e − µ ¯ ψ n +ˆ2 , ψ n, (cid:17) t n, (cid:16) e µ ¯ ψ n, ψ n +ˆ2 , (cid:17) t n, , (15)where the accompanying discrete variables { s, t, x } which describe the ordering of expansioncould be a candidate of new degree of freedom .The next step is to integrate out the original degree of freedom ψ, ¯ ψ and obtain the tensornetwork representation. When performing the integration, it is better to keep the fermionpair structure, which has a common exponent for each fermion field in the pair, to controlsign factors originating from the anticommutation relations. For that purpose, we introduce Actually the discrete variables { s } will not be a new degree of freedom since they will be integratedout as we will see later. { ξ, ¯ ξ, η, ¯ η } , for instance, when integrating ¯ ψ n , ψ n on a site n , (cid:16) ¯ χ n +ˆ1 , χ n, (cid:17) x n, = ( χ n, dη n, ) x n, (cid:16) ¯ χ n +ˆ1 , η n, (cid:17) x n, , (16) (cid:16) ¯ χ n, χ n +ˆ1 , (cid:17) x n, = ( ¯ χ n, d ¯ η n, ) x n, (cid:16) ¯ η n, χ n +ˆ1 , (cid:17) x n, , (17) (cid:16) e − µ ¯ ψ n +ˆ2 , ψ n, (cid:17) t n, = (cid:16) e − µ ψ n, dξ n, (cid:17) t n, (cid:16) e − µ ¯ ψ n +ˆ2 , ξ n, (cid:17) t n, , (18) (cid:16) e µ ¯ ψ n, ψ n +ˆ2 , (cid:17) t n, = (cid:16) e µ ¯ ψ n, d ¯ ξ n, (cid:17) t n, (cid:16) ¯ ξ n, e µ ψ n +ˆ2 , (cid:17) t n, , (19) (cid:16) ¯ χ n, χ n − ˆ1 , (cid:17) x n − ˆ1 , = ( ¯ χ n, d ¯ η n, ) x n − ˆ1 , (cid:16) ¯ η n, χ n − ˆ1 , (cid:17) x n − ˆ1 , , (20) (cid:16) ¯ χ n − ˆ1 , χ n, (cid:17) x n − ˆ1 , = ( χ n, dη n, ) x n − ˆ1 , (cid:16) ¯ χ n − ˆ1 , η n, (cid:17) x n − ˆ1 , , (21) (cid:16) e − µ ¯ ψ n, ψ n − ˆ2 , (cid:17) t n − ˆ2 , = (cid:16) e − µ ¯ ψ n, d ¯ ξ n, (cid:17) t n − ˆ2 , (cid:16) ¯ ξ n, e − µ ψ n − ˆ2 , (cid:17) t n − ˆ2 , , (22) (cid:16) e µ ¯ ψ n − ˆ2 , ψ n, (cid:17) t n − ˆ2 , = (cid:16) e µ ψ n, dξ n, (cid:17) t n − ˆ2 , (cid:16) e µ ¯ ψ n − ˆ2 , ξ n, (cid:17) t n − ˆ2 , . (23)Point here is that we can separate the original degree of freedoms at different site in a differentfermion pair. By collecting all contributions related with ψ n and ¯ ψ n in the partition function(without arising any sign factors) and then the integration for this part is given by Z dψ n, d ¯ ψ n, dψ n, d ¯ ψ n, · X s n, ,s n, ,s n, (cid:0) − ( m + 2) ¯ ψ n, ψ n, (cid:1) s n, (cid:0) − ( m + 2) ¯ ψ n, ψ n, (cid:1) s n, (cid:0) g ¯ ψ n, ψ n, ¯ ψ n, ψ n, (cid:1) s n, · ( χ n, dη n, ) x n, ( ¯ χ n, d ¯ η n, ) x n, (cid:16) e − µ ψ n, dξ n, (cid:17) t n, (cid:16) e µ ¯ ψ n, d ¯ ξ n, (cid:17) t n, · ( ¯ χ n, d ¯ η n, ) x n − ˆ1 , ( χ n, dη n, ) x n − ˆ1 , (cid:16) e − µ ¯ ψ n, d ¯ ξ n, (cid:17) t n − ˆ2 , (cid:16) e µ ψ n, dξ n, (cid:17) t n − ˆ2 , . (24)Note that there is no original fermion field at other sites n ′ = n . This integration can bedone manually and the result depends on the configuration of exponents in an abbreviatedform x n = ( x n, , x n, ) , t n = ( t n, , t n, ) , x n − ˆ1 = ( x n − ˆ1 , , x n − ˆ1 , ) , t n − ˆ2 = ( t n − ˆ2 , , t n − ˆ2 , ) . (25)We rewrite eq.(24) and define the bosonic part T x n t n x n − ˆ1 t n − ˆ2 as follows,(24) = T x n t n x n − ˆ1 t n − ˆ2 d ¯ η x n, n, dη x n, n, d ¯ ξ t n, n, dξ t n, n, dη x n − ˆ1 , n, d ¯ η x n − ˆ1 , n, dξ t n − ˆ2 , n, d ¯ ξ t n − ˆ2 , n, . (26)The explicit form of T x n t n x n − ˆ1 t n − ˆ2 is given in appendix A. By repeating this operation for allother sites, the partition function is finally written in the tensor network representation by Z P = X { x,t } Z Y n T x n t n x n − ˆ1 t n − ˆ2 , (27) In this integration, of course one has to break the pair structure and sign factors appear but thisis still tolerable. T x n t n x n − ˆ1 t n − ˆ2 = T x n t n x n − ˆ1 t n − ˆ2 d ¯ η x n, n, dη x n, n, d ¯ ξ t n, n, dξ t n, n, dη x n − ˆ1 , n, d ¯ η x n − ˆ1 , n, dξ t n − ˆ2 , n, d ¯ ξ t n − ˆ2 , n, · (cid:16) ¯ η n +ˆ1 , η n, (cid:17) x n, (cid:16) ¯ η n, η n +ˆ1 , (cid:17) x n, (cid:16) ¯ ξ n +ˆ2 , ξ n, (cid:17) t n, (cid:16) ¯ ξ n, ξ n +ˆ2 , (cid:17) t n, . (28) In this subsection, we explain the GTRG for the tensor network representation (28). In thesame way as the usual TRG, we decompose the bosonic part on a site n by the SVD andtruncate at D cut , T x n t n x n − ˆ1 t n − ˆ2 ≃ D cut X x n ∗− ˆ1 ∗ ,b =1 U x n t n ,x n ∗− ˆ1 ∗ ,b σ x n ∗− ˆ1 ∗ ,b U ∗ x n − ˆ1 t n − ˆ2 ,x n ∗− ˆ1 ∗ ,b = D cut X x n ∗− ∗ ,b =1 S x n t n x n ∗− ˆ1 ∗ S x n − ˆ1 t n − ˆ2 x n ∗− ˆ1 ∗ , (29)where U , are unitary matrix, σ is singular value and n ∗ is the coarse-grained lattice sitewith unit vectors ˆ1 ∗ = ˆ1 + ˆ2 , ˆ2 ∗ = ˆ1 − ˆ2. See Figure 1 for a graphical representation of thisdecomposition. For Grassmann part, we separate the Grassmann variables into two partswith new Grassmann variables ¯ η n ∗ , η n ∗ − ˆ1 ∗ on the coarse-grained lattice to control the signfactors, d ¯ η x n, n, dη x n, n, d ¯ ξ t n, n, dξ t n, n, dη x n − ˆ1 , n, d ¯ η x n − ˆ1 , n, dξ t n − ˆ2 , n, d ¯ ξ t n − ˆ2 , n, · (cid:16) ¯ η n +ˆ1 , η n, (cid:17) x n, (cid:16) ¯ η n, η n +ˆ1 , (cid:17) x n, (cid:16) ¯ ξ n +ˆ2 , ξ n, (cid:17) t n, (cid:16) ¯ ξ n, ξ n +ˆ2 , (cid:17) t n, = (cid:0) D x n t n d ¯ η x n ∗− ˆ1 ∗ ,f n ∗ (cid:1) (cid:16) D x n − ˆ1 t n − ˆ2 dη x n ∗− ˆ1 ∗ ,f n ∗ − ˆ1 ∗ (cid:17) (cid:0) ¯ η n ∗ η n ∗ − ˆ1 ∗ (cid:1) x n ∗− ˆ1 ∗ ,f , (30)where D and D are defined by D x n t n = d ¯ η x n, n, dη x n, n, d ¯ ξ t n, n, dξ t n, n, · (cid:16) ¯ η n +ˆ1 , η n, (cid:17) x n, (cid:16) ¯ η n, η n +ˆ1 , (cid:17) x n, (cid:16) ¯ ξ n +ˆ2 , ξ n, (cid:17) t n, (cid:16) ¯ ξ n, ξ n +ˆ2 , (cid:17) t n, , (31) D x n − ˆ1 t n − ˆ2 = dη x n − ˆ1 , n, d ¯ η x n − ˆ1 , n, dξ t n − ˆ2 , n, d ¯ ξ t n − ˆ2 , n, , (32)and new exponent x n ∗ − ˆ1 ,f with one-component is introduced with constraints, x n ∗ − ˆ1 ∗ ,f = X i ( x n,i + t n,i ) mod 2 = X i (cid:16) x n − ˆ1 ,i + t n − ˆ2 ,i (cid:17) mod 2 . (33)Then the tensor T x n t n x n − ˆ1 t n − ˆ2 is decomposed and approximated as T x n t n x n − ˆ1 t n − ˆ2 ≃ D cut X x n ∗− ˆ1 ∗ ,b =1 1 X x n ∗− ˆ1 ∗ ,f =0 Z S x n t n x n ∗− ˆ1 ∗ S x n − ˆ1 t n − ˆ2 x n ∗− ˆ1 ∗ (cid:0) ¯ η n ∗ η n ∗ − ˆ1 ∗ (cid:1) x n ∗− ˆ1 ∗ ,f · δ P i ( x n,i + t n,i ) mod 2 ,x n ∗− ˆ1 ∗ ,f δ P i ( x n − ˆ1 ,i + t n − ˆ2 ,i ) mod 2 ,x n ∗− ˆ1 ∗ ,f , (34)where x n ∗ = ( x n ∗ ,b , x n ∗ ,f ) and S x n t n x n ∗− ˆ1 ∗ = S x n t n x n ∗− ˆ1 ∗ ,b D x n t n d ¯ η x n ∗− ˆ1 ∗ ,f n ∗ , (35) S x n − ˆ1 t n − ˆ2 x n ∗− ˆ1 ∗ = S x n − ˆ1 t n − ˆ2 x n ∗− ˆ1 ∗ ,b D x n − ˆ1 t n − ˆ2 dη x n ∗− ˆ1 ∗ ,f n ∗ − ˆ1 ∗ . (36)6 ig. 1 The decomposition of tensor. The horizontal (vertical) axis corresponds to 1-direction (2-direction).For another decomposition rotated 90 degree (See Figure 1), by introducing new variables¯ ξ n ∗ and ξ n ∗ − ˆ2 ∗ , it is similarly given by T x n +ˆ2 t n +ˆ2 x n − ˆ1+ˆ2 t n ≃ D cut X t n ∗− ˆ2 ∗ ,b =1 1 X t n ∗− ˆ2 ∗ ,f =0 Z S t n x n +ˆ2 t n ∗− ˆ2 ∗ S t n +ˆ2 x n − ˆ1+ˆ2 t n ∗− ˆ2 ∗ (cid:0) ¯ ξ n ∗ ξ n ∗ − ˆ2 ∗ (cid:1) t n ∗− ˆ2 ∗ ,f · δ P i ( t n,i + x n +ˆ2 ,i ) mod 2 ,t n ∗− ˆ2 ∗ ,f δ P i ( t n +ˆ2 ,i + x n − ˆ1+ˆ2 ,i ) mod 2 ,t n ∗− ˆ2 ∗ ,f , (37)where S and S are defined by S t n x n +ˆ2 t n ∗− ˆ2 ∗ = S t n x n +ˆ2 t n ∗− ˆ2 ∗ ,b D t n x n +ˆ2 d ¯ ξ t n ∗− ˆ2 ∗ ,f n ∗ , (38) S t n +ˆ2 x n − ˆ1+ˆ2 t n ∗− ˆ2 ∗ = S t n +ˆ2 x n − ˆ1+ˆ2 t n ∗− ˆ2 ∗ ,b D t n +ˆ2 x n − ˆ1+ˆ2 dξ t n ∗− ˆ2 ∗ ,f n ∗ − ˆ2 ∗ , (39)with D t n x n +ˆ2 = dξ t n, n +ˆ2 , d ¯ ξ t n, n +ˆ2 , d ¯ η x n +ˆ2 , n +ˆ2 , dη x n +ˆ2 , n +ˆ2 , · (cid:16) ¯ η n +ˆ1+ˆ2 , η n +ˆ2 , (cid:17) x n +ˆ2 , (cid:16) ¯ η n +ˆ2 , η n +ˆ1+ˆ2 , (cid:17) x n +ˆ2 , , (40) D t n +ˆ2 x n − ˆ1+ˆ2 = d ¯ ξ t n +ˆ2 , n +ˆ2 , dξ t n +ˆ2 , n +ˆ2 , dη x n − ˆ1+ˆ2 , n +ˆ2 , d ¯ η x n − ˆ1+ˆ2 , n +ˆ2 , · (cid:16) ¯ ξ n +2 · ˆ2 , ξ n +ˆ2 , (cid:17) t n +ˆ2 , (cid:16) ¯ ξ n +ˆ2 , ξ n +2 · ˆ2 , (cid:17) t n +ˆ2 , . (41)The bosonic part S and S are determined by the SVD as follows M t n x n +ˆ2 ,t n +ˆ2 x n − ˆ1+ˆ2 = ( − t n, + t n, T x n +ˆ2 t n +ˆ2 x n − ˆ1+ˆ2 t n = X t n ∗− ˆ2 ∗ ,b S t n x n +ˆ2 t n ∗− ˆ2 ∗ ,b S t n +ˆ2 x n − ˆ1+ˆ2 t n ∗− ˆ2 ∗ ,b . (42)7 ig. 2 The contraction of original indices. The broken lines are the original lattice andthe bold lines are the coarse-grained lattice. The solid lines indicate contracted indices.A coarse-grained tensor is obtained by T x n ∗ t n ∗ x n ∗− ˆ1 ∗ t n ∗− ˆ2 ∗ = Z X { x n ,t n } S x n t n x n ∗− ˆ1 ∗ S t n x n +ˆ2 t n ∗− ˆ2 ∗ S x n +ˆ2 t n +ˆ1 x n ∗ S t n +ˆ1 x n t n ∗ · (cid:0) ¯ η n ∗ +ˆ1 ∗ η n ∗ (cid:1) x n ∗ ,f (cid:0) ¯ ξ n ∗ +ˆ2 ∗ ξ n ∗ (cid:1) t n ∗ ,f · δ P i ( x n,i + t n,i ) mod 2 ,x n ∗− ˆ1 ∗ ,f δ P i ( t n,i + x n +ˆ2 ,i ) mod 2 ,t n ∗− ˆ2 ∗ ,f · δ P i ( x n +ˆ2 ,i + t n +ˆ1 ,i ) mod 2 ,x n ∗ ,f δ P i ( t n +ˆ1 ,i + x n,i ) mod 2 ,t n ∗ ,f = T x n ∗ t n ∗ x n ∗− ˆ1 ∗ t n ∗− ˆ2 ∗ dη x n ∗ ,f dξ t n ∗ ,f d ¯ η x n ∗− ˆ1 ∗ ,f d ¯ ξ t n ∗− ˆ2 ∗ ,f · (cid:0) ¯ η n ∗ +ˆ1 ∗ η n ∗ (cid:1) x n ∗ ,f (cid:0) ¯ ξ n ∗ +ˆ2 ∗ ξ n ∗ (cid:1) t n ∗ ,f . (43)Note that constraint δ ,x n ∗ ,f + t n ∗ ,f + x n ∗− ˆ1 ∗ ,f + t n ∗− ˆ2 ∗ ,f mod 2 is imposed for the coarse-grainedtensor. Figure 2 shows the contraction for the original indices in this renormalization step.Repeat this renormalization step until the number of lattice point reaches 2 ×
2, namelyfour reduced tensors. From these tensors, the parition function is computed by full indexcontractions.Computational costs of a standard SVD routine are proportional to the third power ofthe matrix size, thus the cost of the decomposition of a tensor is of order D . On theother hand, the cost of the contraction is of order D . Therefore the total cost of GTRGis proportional to D . From here, we consider a system where the anti-periodic (periodic) boundary conditionis imposed for the 2-direction (1-direction). This is taken into account by modifying the8artition function Z = X { x,t,t ′ } Z Y n T x n t n x n − ˆ1 t ′ n B t ′ n t n − ˆ2 , (44)where the full tensor T is the same as before and the new matrix B is given by B t ′ n t n − ˆ2 = ( ( − t ′ n, + t ′ n, δ t ′ n ,t n − ˆ2 if n = 0 ,δ t ′ n ,t n − ˆ2 else . (45)As a result, a coarse-grained tensor contracted on a site with n = N − A x n ∗ t n ∗ x n ∗− ˆ1 ∗ t n ∗− ˆ2 ∗ = Z X { x n ,t n } ( − t n, + t n, + t n +ˆ1 , + t n +ˆ1 , S x n t n x n ∗− ˆ1 ∗ S t n x n +ˆ2 t n ∗− ˆ2 ∗ S x n +ˆ2 t n +ˆ1 x n ∗ S t n +ˆ1 x n t n ∗ · (cid:0) ¯ η n ∗ +ˆ1 ∗ η n ∗ (cid:1) x n ∗ ,f (cid:0) ¯ ξ n ∗ +ˆ2 ∗ ξ n ∗ (cid:1) t n ∗ ,f · δ P i ( x n,i + t n,i ) mod 2 ,x n ∗− ˆ1 ∗ ,f δ P i ( t n,i + x n +ˆ2 ,i ) mod 2 ,t n ∗− ˆ2 ∗ ,f · δ P i ( x n +ˆ2 ,i + t n +ˆ1 ,i ) mod 2 ,x n ∗ ,f δ P i ( t n +ˆ1 ,i + x n,i ) mod 2 ,t n ∗ ,f = Z X { x n ,t n } ( − t n, + t n, + t n +ˆ1 , + t n +ˆ1 , S x n t n x n ∗− ˆ1 ∗ S t n x n +ˆ2 t n ∗− ˆ2 ∗ S x n +ˆ2 t n +ˆ1 x n ∗ S t n +ˆ1 x n t n ∗ · (cid:0) ¯ η n ∗ +ˆ1 ∗ η n ∗ (cid:1) x n ∗ ,f (cid:0) ¯ ξ n ∗ +ˆ2 ∗ ξ n ∗ (cid:1) t n ∗ ,f · δ P i ( x n,i + t n,i ) mod 2 ,x n ∗− ˆ1 ∗ ,f δ P i ( t n,i + x n +ˆ2 ,i ) mod 2 ,t n ∗− ˆ2 ∗ ,f · δ P i ( t n,i + t n +ˆ1 ,i ) mod 2 , ( x n ∗ ,f + t n ∗− ˆ2 ∗ ,f ) mod 2 δ P i ( t n +ˆ1 ,i + x n,i ) mod 2 ,t n ∗ ,f = ( − x n ∗ ,f + t n ∗− ˆ2 ∗ ,f T x n ∗ t n ∗ x n ∗− ˆ1 ∗ t n ∗− ˆ2 ∗ dη x n ∗ ,f dξ t n ∗ ,f d ¯ η x n ∗− ˆ1 ∗ ,f d ¯ ξ t n ∗− ˆ2 ∗ ,f · (cid:0) ¯ η n ∗ +ˆ1 ∗ η n ∗ (cid:1) x n ∗ ,f (cid:0) ¯ ξ n ∗ +ˆ2 ∗ ξ n ∗ (cid:1) t n ∗ ,f . (46)Therefore the once renormalized partition function is obtained by Z (1) = X { x,t } Z Y n A x n t n x n − ˆ1 t n − ˆ2 = X { x,x ′ ,t,t ′ } Z Y n T x ′ n t n x n − ˆ1 t ′ n B x ′ n x n B t ′ n t n − ˆ2 (47)where the site indices are replaced by n ∗ → n for a readability and another boundary matricesare given by B x ′ n x n = ( ( − x ′ n,f δ x ′ n ,x n if n = n ,δ x ′ n ,x n else , (48) B t ′ n t n − ˆ2 = ( ( − t ′ n,f δ t ′ n ,t n − ˆ2 if n = n ,δ t ′ n ,t n − ˆ2 else . (49)9imilarly, a twice coarse-grained tensor contracted on n = n is obtained by A x n ∗ t n ∗ x n ∗− ˆ1 ∗ t n ∗− ˆ2 ∗ = Z X { x n ,t n } ( − x n +ˆ2 ,f + t n,f S x n t n x n ∗− ˆ1 ∗ S t n x n +ˆ2 t n ∗− ˆ2 ∗ S x n +ˆ2 t n +ˆ1 x n ∗ S t n +ˆ1 x n t n ∗ · (cid:0) ¯ η n ∗ +ˆ1 ∗ η n ∗ (cid:1) x n ∗ ,f (cid:0) ¯ ξ n ∗ +ˆ2 ∗ ξ n ∗ (cid:1) t n ∗ ,f · δ ( x n,f + t n,f ) mod 2 ,x n ∗− ˆ1 ∗ ,f δ ( t n,f + x n +ˆ2 ,f ) mod 2 ,t n ∗− ˆ2 ∗ ,f · δ ( x n +ˆ2 ,f + t n +ˆ1 ,f ) mod 2 ,x n ∗ ,f δ ( t n +ˆ1 ,f + x n,f ) mod 2 ,t n ∗ ,f = Z X { x n ,t n } ( − x n +ˆ2 ,f + t n,f S x n t n x n ∗− ˆ1 ∗ S t ′ n x n +ˆ2 t n ∗− ˆ2 ∗ S x n +ˆ2 t n +ˆ1 x n ∗ S t n +ˆ1 x n t n ∗ · (cid:0) ¯ η n ∗ +ˆ1 ∗ η n ∗ (cid:1) x n ∗ ,f (cid:0) ¯ ξ n ∗ +ˆ2 ∗ ξ n ∗ (cid:1) t n ∗ ,f · δ ( x n,f + t n,f ) mod 2 ,x n ∗− ˆ1 ∗ ,f δ ( t n,f + x n +ˆ2 ,f ) mod 2 ,t n ∗− ˆ2 ∗ ,f · δ ( x n +ˆ2 ,f + t n +ˆ1 ,f ) mod 2 ,x n ∗ ,f δ ( x n +ˆ2 ,f + t n,f ) mod 2 , ( x n ∗ ,f + t n ∗ ,f + x n ∗− ˆ1 ∗ ,f ) mod 2 = ( − t n ∗− ˆ2 ∗ ,f T x n ∗ t n ∗ x n ∗− ˆ1 ∗ t n ∗− ˆ2 ∗ dη x n ∗ ,f dξ t n ∗ ,f d ¯ η x n ∗− ˆ1 ∗ ,f d ¯ ξ t n ∗− ˆ2 ∗ ,f · (cid:0) ¯ η n ∗ +ˆ1 ∗ η n ∗ (cid:1) x n ∗ ,f (cid:0) ¯ ξ n ∗ +ˆ2 ∗ ξ n ∗ (cid:1) t n ∗ ,f . (50)and the twice renormalized partition function is obtained by Z (2) = X { x,t, } Z Y n A x n t n x n − ˆ1 t n − ˆ2 = X { x,t,t ′ } Z Y n T x n t n x n − ˆ1 t ′ n B t ′ n t n − ˆ2 (51)where B t ′ n t n − ˆ2 = ( ( − t ′ n,f δ t ′ n ,t n − ˆ2 if n = 0 ,δ t ′ n ,t n − ˆ2 else . (52)Therefore, in this formulation, the boundary condition returns to the original one every 2renormalization steps.
3. Numerical Results
First, we compare the numerical results of ln Z with the exact value ln Z exact in the freemassless case. Figure 3 shows the relative deviation δ ( D cut ) = ln Z ( D cut ) − ln Z exact ln Z exact , (53)as a function of D cut . The convergence behavior is roughly observed although it is not sosmooth. The convergence rate at µ = 1 is slower than that of µ = 2. For µ = 2, lattice volumedependence is not seen while for µ = 1 larger volume is strongly affected by truncation error.To see the convergence issue in more detail, we investigate the spectrum of bosonic tensorin Figure 4. Clear hierarchy is observed for µ = 2 while nearly degenerated structure isseen for µ = 1 especially after several iterations. Figure 5 shows the relative deviation as afunction of µ with fixed D cut = 64. The deviation rapidly increases around µ ≈ . | δ | D cut m=0,g=0 µ =1,N =N =8 µ =1,N =N =32 µ =2,N =N =8 µ =2,N =N =32 Fig. 3
The relative deviation δ as a function of D cut for free massless case. σ i / σ im=0,g=0 µ =1,3 iters. µ =1,7 iters. µ =2,3 iters. µ =2,7 iters. Fig. 4
Spectrum of the bosonic tensor for free massless case.Next, we compute the fermion number density defined as n = 1 N N ∂ ln Z∂µ . (54)11 | δ | µ m=0,g=0,D cut =64 N =N =4N =N =8N =N =16N =N =32 Fig. 5
The relative deviation δ as a function of µ with fixed D cut = 64 for free masslesscase. Around µ = 0 . N = N = 4, the TRG resultbecomes exact, thus the relative deviation is exactly zero up to machine precision thus thisshows a validity of our calculation.Figure 6 plots the fermion number density as a function of µ for some non-trivial sets ofparameters. Since the model is in two dimensional system with one-flavor, the saturationdensity for fermion number is one. We observe that the fermion density saturates to thisvalue for larger chemical potential.Finally, we perform the finite size scaling analysis for the quark number susceptibilitydefined as χ = 1 N N ∂ ln Z∂µ . (55)The susceptibility as a function of µ is shown in Figure 7 for various spatial volumes withtwo sets of parameter ( m, g ) = (0 ,
0) and (0 , . µ = 1 and the peak height shows no volume dependence, therefore we concludethat this transition is cross-over. For lower µ . .
6, the TRG results develop some peaksfor both couplings. In order to check whether these peaks are fake or not, we compare withthe exact results at g = 0 shown as curves for each volume N = 32 , ,
96 where for largervolume the peaks disappear in the lower µ region. From the comparison, we find that theTRG results at g = 0, shown as dots, tend to deviate from these curves for larger volume.Thus we conclude that these peaks at g = 0 of TRG results especially with larger volume arefake. For g = 0 .
7, since we cannot directly compare with the exact results, we are contentwith being comparing two results obtained by different resolution of the chemical potentialin the numerical derivative. And then the difference is barely seen thus we expect that thepeak around µ = 0 . g = 0 . µ ≈ n µ N =N =32,D cut =64 m=0,g=0m=0,g=0.7m=0.5,g=0.7saturation density Fig. 6
Fermion number density n as a function of µ with fixed N = N = 32 and D cut =64. For larger chemical potential, the number density for all cases of ( m, g ) we investigatedsaturates to unity as expected.relative deviation in Figure 5 becomes large. It has been known that the approximation forTRG gets worse near a critical point, while we observe that such a behavior occurs even forcross-over case. Needless to say, on another parameters ( N , N , m, g ), real phase transitioncan occur and the strength of transition may change, thus the source of the loss of accuracywe observed here could be a remnant of the real phase transition.
4. Reweighting Method
In the TRG calculation, one usually computes the partition function at several parameterpoints (mesh). Then numerical derivative of partition function with respect to the parameteris made by using a few points and one needs a fine mesh to reduce a discretization error.To reduce the computational time, we propose a method to obtain an approximated coarse-grained tensor at one parameter by using another set of singular values at different parameter.Using an analogy from Monte Carlo method, we refer to this method as the reweightingmethod.Let the bosonic part T x n t n x n − ˆ1 t n − ˆ2 at the original parameter and its SVD is given by T x n t n x n − ˆ1 t n − ˆ2 = X x n ∗− ˆ1 ∗ ,b ,x ′ n ∗− ˆ1 ∗ ,b U x n t n ,x n ∗− ˆ1 ∗ ,b σ x n ∗− ˆ1 ∗ ,b δ x n ∗− ˆ1 ∗ ,b ,x ′ n ∗− ˆ1 ∗ ,b U ∗ x n − ˆ1 t n − ˆ2 ,x ′ x ∗− ˆ1 ∗ ,b . (56)13 χ µ m=0,N =32,D cut =64 exact,N =32exact,N =64exact,N =96g=0,N =32g=0,N =64g=0,N =96g=0.7,N =32g=0.7,N =64g=0.7,N =96 Fig. 7
Finite size scaling of the fermion number susceptibility for fixed ( m, N , D cut ) =(0 , , N -dependence is not observed even in the presence of interaction g = 0 . T ′ x n t n x n − ˆ1 t n − ˆ2 can be written as T ′ x n t n x n − ˆ1 t n − ˆ2 = X x,t,x ′ ,t ′ ,x n ∗− ˆ1 ∗ ,b ,x ′ n ∗− ˆ1 ∗ ,b (cid:16) U x n t n ,x n ∗− ˆ1 ∗ ,b U ∗ xt,x n ∗− ˆ1 ∗ ,b (cid:17) T ′ xtx ′ t ′ (cid:16) U x ′ t ′ ,x ′ n ∗− ˆ1 ∗ ,b U ∗ x n − ˆ1 t n − ˆ2 ,x ′ n ∗− ˆ1 ∗ ,b (cid:17) = X x n ∗− ˆ1 ∗ ,b ,x ′ n ∗− ˆ1 ∗ ,b U x n t n ,x n ∗− ˆ1 ∗ ,b Σ x n ∗− ˆ1 ∗ ,b ,x ′ n ∗− ˆ1 ∗ ,b U ∗ x n − ˆ1 t n − ˆ2 ,x ′ n ∗− ˆ1 ∗ ,b , (57)where the new matrix Σ is given byΣ x n ∗− ˆ1 ∗ ,b ,x ′ n ∗− ˆ1 ∗ ,b = X x,t,x ′ ,t ′ U ∗ xt,x n ∗− ˆ1 ∗ ,b T ′ xtx ′ t ′ U x ′ t ′ ,x ′ n ∗− ˆ1 ∗ ,b . (58)By truncating the indices x n ∗ − ˆ1 ∗ ,b , x ′ n ∗ − ˆ1 ∗ ,b at D cut , the decomposition of T ′ can be formally defined by T ′ x n t n x n − ˆ1 t n − ˆ2 ≃ D cut X x n ∗− ˆ1 ∗ ,b ,x ′ n ∗− ˆ1 ∗ ,b =1 U x n t n ,x n ∗− ˆ1 ∗ ,b Σ x n ∗− ˆ1 ∗ ,b ,x ′ n ∗− ˆ1 ∗ ,b U ∗ x n − ˆ1 t n − ˆ2 ,x ′ n ∗− ˆ1 ∗ ,b . (59)Similarly, for another decomposition M t n x n +ˆ2 ,t n +ˆ2 x n − ˆ1+ˆ2 = ( − t n, + t n, T x n +ˆ2 t n +ˆ2 x n − ˆ1+ˆ2 t n = X t n ∗− ˆ2 ∗ ,b ,t ′ n ∗− ˆ2 ∗ ,b U t n x n +ˆ2 ,t n ∗− ˆ2 ∗ ,b σ t n ∗− ˆ2 ∗ ,b δ t n ∗− ˆ2 ∗ ,b ,t ′ n ∗− ˆ2 ∗ ,b U ∗ t n +ˆ2 x n − ˆ1+ˆ2 ,t ′ n ∗− ˆ2 ∗ ,b , (60) This decomposition is not optimal since this is not SVD of T ′ . is defined byΣ t n ∗− ˆ2 ∗ ,b ,t ′ n ∗− ˆ2 ∗ ,b = X x,t,x ′ ,t ′ ( − t ′ + t ′ U ∗ t ′ x,t n ∗− ˆ2 ∗ ,b T ′ xtx ′ t ′ U tx ′ ,t ′ n ∗− ˆ2 ∗ ,b . (61)By using the singular vectors U , , , at the original parameter , an intermediate tensor isdefined by˜ T x ′ n ∗ t ′ n ∗ x n ∗− ˆ1 ∗ t n ∗− ˆ2 ∗ = Z X { x n ,t n } U x n t n x n ∗− ˆ1 ∗ U t n x n +ˆ2 t n ∗− ˆ2 ∗ U x n +ˆ2 t n +ˆ1 x ′ n ∗ U t n +ˆ1 x n t ′ n ∗ · (cid:0) ¯ η n ∗ +ˆ1 ∗ η n ∗ (cid:1) x ′ n ∗ ,f (cid:0) ¯ ξ n ∗ +ˆ2 ∗ ξ n ∗ (cid:1) t ′ n ∗ ,f · δ P i ( x n,i + t n,i ) mod 2 ,x n ∗− ˆ1 ∗ ,b δ P i ( t n,i + x n +ˆ2 ,i ) mod 2 ,t n ∗− ˆ2 ∗ ,b · δ P i ( x n +ˆ2 ,i + t n +ˆ1 ,i ) mod 2 ,x ′ n ∗ ,b δ P i ( t n +ˆ1 ,i + x n,i ) mod 2 ,t ′ n ∗ ,b = ˜ T x ′ n ∗ t ′ n ∗ x n ∗− ˆ1 ∗ t n ∗− ˆ2 ∗ dη x ′ n ∗ ,f dξ t ′ n ∗ ,f d ¯ η x n ∗− ˆ1 ∗ ,f d ¯ ξ t n ∗− ˆ2 ∗ ,f · (cid:0) ¯ η n ∗ +ˆ1 ∗ η n ∗ (cid:1) x ′ n ∗ ,f (cid:0) ¯ ξ n ∗ +ˆ2 ∗ ξ n ∗ (cid:1) t ′ n ∗ ,f , (62)where U x n t n x n ∗− ˆ1 ∗ = U x n t n ,x n ∗− ˆ1 ∗ ,b D x n t n d ¯ η x n ∗− ˆ1 ∗ ,f n ∗ , (63) U t n x n +ˆ2 t n ∗− ˆ2 ∗ = U t n x n +ˆ2 ,t n ∗− ˆ2 ∗ ,b D t n x n +ˆ2 d ¯ ξ t n ∗− ˆ2 ∗ ,f n ∗ , (64) U x n +ˆ2 t n +ˆ1 x ′ n ∗ = U x n +ˆ2 t n +ˆ1 ,x ′ n ∗ ,b D x n +ˆ2 t n +ˆ1 dη x ′ n ∗ ,f n ∗ , (65) U t n +ˆ1 x n t ′ n ∗ = U t n +ˆ1 x n ,t ′ n ∗ ,b D t n +ˆ1 x n dξ t ′ n ∗ ,f n ∗ . (66)From this intermediate tensor, we can obtain not only the original coarse-grained tensor T x n ∗ t n ∗ x n ∗− ˆ1 ∗ t n ∗− ˆ2 ∗ = X x ′ n ∗ ,t ′ n ∗ σ x n ∗ ,b δ x n ∗ ,x ′ n ∗ σ t n ∗ ,b δ t n ∗ ,t ′ n ∗ ˜ T x ′ n ∗ t ′ n ∗ x n ∗− ˆ1 ∗ t n ∗− ˆ2 ∗ , (67)but also a coarse-grained tensor at different parameter T ′ x n ∗ t n ∗ x n ∗− ˆ1 ∗ t n ∗− ˆ2 ∗ = X x ′ n ∗ ,t ′ n ∗ Σ x n ∗ ,b ,x ′ n ∗ ,b δ x n ∗ ,f ,x ′ n ∗ ,f Σ t n ∗ ,b ,t ′ n ∗ ,b δ t n ∗ ,f ,t ′ n ∗ ,f ˜ T x ′ n ∗ t ′ n ∗ x n ∗− ˆ1 ∗ t n ∗− ˆ2 ∗ , (68)which is not optimal but is approximately fine if the difference of the original and targetparameter is small. If the intermediate tensor and the singular vectors at the original param-eter have been stored, a coarse-grained tensor for another parameters is obtained by only D order computational cost. The relative deviation δ = ln Z RW − ln Z exact ln Z exact , (69)between ln Z RW computed by using the reweighting method and the exact one is shown inFigure 8. The relative deviation increases as the distance from original parameters and thelattice size. The deviation reweighting from nearly transition point ( µ = 1) quickly increasecompared with that of off-transition ( µ = 0 , µ = 2). The corresponding singular values are not included. δ µ m=0,g=0 N =N =8N =N =16N =N =32-0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.5 1 1.5 2 δ µ m=0,g=0 N =N =8N =N =16N =N =32-0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.5 1 1.5 2 δ µ m=0,g=0 N =N =8N =N =16N =N =32 Fig. 8
The relative deviation between reweighting method and exact value in eq.(69) asa function of µ for free massless case with fixed D cut = 64. From top to bottom, the originalvalue of µ is given by µ = 0, 1 and 2 respectively.16 . Summury and Outlook We have applied the GTRG to the one-flavor lattice Gross-Neveu model with chemicalpotential in the Wilson fermion formulation. At some non-trivial parameter set at finitedensity, we found a transition-like behavior and the finite size scaling study shows that thistransition is a cross-over but not a real phase transition. Furthermore, we observed thataround the “transition” point the approximation of TRG gets worse, although this is not acritical point.We introduced the reweighting method for TRG and demonstrated for some parameters.As a result, the errors increase as the distance from original parameters and the latticesize. Furthermore we observed that the reweighting from around “transition” point quicklydeteriorates compared with reweighting from off-transition region.This is the first application of the GTRG to finite density sistem. We hope that theformulation given in this work is extended another finite density systems.
ACKNOWLEDGMENTS
We would like to thank Y. Shimizu and D. Satou for helpful advice. S.T. is grateful to Y.Kuramashi for useful conversation. This work is partially supported by the Grants-in-Aid forScientific Research from the Ministry of Education, Culture, Sports, Science and Technology(Nos. 26800130).
A. Details of bosonic tensor
In this appendix, we show explicit elements of bosonic tensor T x n t n x n − ˆ1 t n − ˆ2 in eq.(26). T = ( m + 2) + 2 g , T = − ( m +2) √ e − µ , T = ( m +2) √ e µ , T = − ,T = − ( m + 2) , T = − ( m +2) √ e − µ , T = − √ e − µ , T = − ,T = ( m + 2) , T = − ( m +2) √ e µ , T = √ e µ , T = − ,T = 1 , T = √ e − µ , T = √ e µ , T = − ,T = − ( m +2) √ e − µ , T = − ( m + 2) e − µ , T = − e − µ , T = − √ e − µ ,T = − e − µ , T = ( m +2) √ e − µ , T = √ e − µ , T = e − µ ,T = 1 , T = − √ e − µ , T = √ e − µ , T = e − µ ,T = − ( m +2) √ e µ , T = ( m + 2) e µ , T = − e µ , T = √ e µ ,T = − ( m +2) √ e µ , T = √ e µ , T = 1 , T = e µ ,T = − √ e µ , T = − e µ , T = − √ e µ , T = e µ ,T = − , T = − √ e − µ , T = − √ e µ , T = 1 ,T = , T = √ e − µ , T = , T = − √ e µ ,T = − , others = 0 . Note that the four-fermion coupling g enters only in the first element. References [1] M. Levin and C. P. Nave, Phys. Rev. Lett. 99, 120601 (2007)[2] J. F. Yu, Z. Y. Xie, Y. Meurice, Y. Liu, A. Denbleyker, H. Zou, M. P. Qin, J. Chen, and T. Xiang,Phys. Rev. E 89, 013308 (2014)[3] J. Unmuth-Yockey, Y. Meurice, J. Osborn, H. Zou, arXiv:1411.4213 [hep-lat]
4] Y. Shimizu, Mod. Phys. Lett. A 27, 1250035 (2012)[5] Z. Y. Xie, H. C. Jiang, Q. N. Chen, Z. Y. Weng, T. Xiang, Phys. Rev. Lett. 103 160601 (2009)[6] Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang, T. Xiang, Phys. Rev. B 86, 045139 (2012)[7] Z.-C Gu, F. Verstraete, X.-G Wen, arXiv:1004.2563 [cond-mat.str-el][8] Z.-C. Gu, Phys. Rev. B 88, 115139 (2013)[9] Y. Shimizu and Y. Kuramashi, Phys. Rev. D. 90, 014508 (2014)[10] Y. Shimizu and Y. Kuramashi, Phys. Rev. D. 90, 074503 (2014)[11] D. J. Gross, A. Neveu, Phys. Rev. D 10, 3235 (1974)[12] T. Eguchi and R. Nakayama, Phys. Lett. B 126, 89 (1983).[13] S. Aoki and K. Higashijima, Prog. Theor. Phys. 76, 521 (1986).[14] T. Izubuchi, J. Noaki and A. Ukawa, Phys. Rev. D. 58, 114507 (1998)[15] J. B. Kogut, H. Matsuoka, M. Stone, H. W. Wyld, S. Shenker, J. Shigemitsu, D. K. Sinclair, Nucl.Phys. B225 [FS], 93 (1983) ; P. Hasenfratz F. Karsch, Phys. Lett. 125B, 308 (1983)4] Y. Shimizu, Mod. Phys. Lett. A 27, 1250035 (2012)[5] Z. Y. Xie, H. C. Jiang, Q. N. Chen, Z. Y. Weng, T. Xiang, Phys. Rev. Lett. 103 160601 (2009)[6] Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang, T. Xiang, Phys. Rev. B 86, 045139 (2012)[7] Z.-C Gu, F. Verstraete, X.-G Wen, arXiv:1004.2563 [cond-mat.str-el][8] Z.-C. Gu, Phys. Rev. B 88, 115139 (2013)[9] Y. Shimizu and Y. Kuramashi, Phys. Rev. D. 90, 014508 (2014)[10] Y. Shimizu and Y. Kuramashi, Phys. Rev. D. 90, 074503 (2014)[11] D. J. Gross, A. Neveu, Phys. Rev. D 10, 3235 (1974)[12] T. Eguchi and R. Nakayama, Phys. Lett. B 126, 89 (1983).[13] S. Aoki and K. Higashijima, Prog. Theor. Phys. 76, 521 (1986).[14] T. Izubuchi, J. Noaki and A. Ukawa, Phys. Rev. D. 58, 114507 (1998)[15] J. B. Kogut, H. Matsuoka, M. Stone, H. W. Wyld, S. Shenker, J. Shigemitsu, D. K. Sinclair, Nucl.Phys. B225 [FS], 93 (1983) ; P. Hasenfratz F. Karsch, Phys. Lett. 125B, 308 (1983)