Restoration of chiral symmetry in cold and dense Nambu--Jona-Lasinio model with tensor renormalization group
Shinichiro Akiyama, Yoshinobu Kuramashi, Takumi Yamashita, Yusuke Yoshimura
PPrepared for submission to JHEP
UTHEP-752, UTCCS-P-134
Restoration of chiral symmetry in cold and denseNambu–Jona-Lasinio model with tensorrenormalization group
Shinichiro Akiyama, a Yoshinobu Kuramashi, b Takumi Yamashita, c Yusuke Yoshimura b a Graduate School of Pure and Applied Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8571,Japan b Center for Computational Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8577, Japan c Faculty of Engineering, Information and Systems, University of Tsukuba, Tsukuba, Ibaraki 305-8573, Japan
E-mail: [email protected] , [email protected] , [email protected] , [email protected] Abstract:
We analyze the chiral phase transition of the Nambu–Jona-Lasinio model in the coldand dense region on the lattice, developing the Grassmann version of the anisotropic ten-sor renormalization group algorithm. The model is formulated with the Kogut–Susskindfermion action. We use the chiral condensate as an order parameter to investigate therestoration of the chiral symmetry. The first-order chiral phase transition is clearly ob-served in the dense region at vanishing temperature with µ/T ∼ O (10 ) on a large volumeof V = 1024 . We also present the results for the equation of state. a r X i v : . [ h e p - l a t ] D ec ontents The phase structure and the equation of state for QCD at finite temperature and densityare essential ingredients to understand the evolution and the current state of the universequantitatively. Although the lattice QCD simulation has been expected to be an idealtool to investigate the non-perturbative aspects of QCD, it has not been successful toreveal the nature of QCD at finite density. This is due to the sign problem caused by theintroduction of the chemical potential in the lattice QCD simulations based on the MonteCarlo algorithm, see, e.g. , Ref. [1].The tensor renormalization group (TRG) method, which was originally proposed byLevin and Nave to study two-dimensional (2 d ) classical spin models in 2007 [2], has severalsuperior features over the Monte Carlo method. Firstly, the TRG method, and alsoother tensor network methods, are free from the sign problem. This virtue was confirmedby successful application of these methods to various 2 d quantum field theories whichcontain the sign problem [6, 8–13]. Moreover, the authors have successfully employed theanisotropic TRG (ATRG) algorithm to analyze the Bose condensation in the 4 d complex φ theory at finite density [16]. Secondly, the computational cost depends on the system size In this paper the TRG method or the TRG approach refers to not only the original numerical algorithmproposed by Levin and Nave but also its extensions [3–7]. The TRG study of the sign problem was also performed by introducing the complex coupling or thechemical potential to the 2 d classical O (2) spin model [14, 15]. – 1 –nly logarithmically. This is an essential ingredient to enable us to take the thermodynamiclimit at vanishing temperature. Thirdly, the Grassmann version of the TRG method,which was developed by some of the authors [6, 7, 17], allows direct manipulation of theGrassmann variables. Fourthly, we can obtain the partition function or the path-integralitself. In the thermodynamic limit, the pressure is directly related to the thermodynamicpotential so that the equation of state can be easily obtained with the TRG method.In this paper we investigate the phase structure of the Nambu–Jona-Lasinio (NJL)model [18, 19] at finite temperature T and chemical potential µ on the lattice, developingthe Grassmann version of the ATRG algorithm. The Lagrangian of the NJL model in thecontinuum is defined as follows: L = ¯ ψ ( x ) γ ν ∂ ν ψ ( x ) − g (cid:8) ( ¯ ψ ( x ) ψ ( x )) + ( ¯ ψ ( x )i γ ψ ( x )) (cid:9) , (1.1)which has the U(1) chiral symmetry with ψ ( x ) → e i αγ ψ ( x ) and ¯ ψ ( x ) → ¯ ψ ( x )e i αγ . Thisis an effective theory of QCD which describes the dynamical chiral symmetry breaking:once the strength of the coupling constant g exceeds a certain critical value the systemgenerates a non-trivial vacuum with (cid:104) ¯ ψ ( x ) ψ ( x ) (cid:105) (cid:54) = 0. The chiral phase structure of theNJL model on the T - µ plane is discussed by some analytical methods, e.g. , the mean-fieldapproximation (MFA) [20] and the functional renormalization group (FRG) [21]. Figure 1shows a schematic view of the expected phase structure, whose characteristic feature isthe first-order chiral phase transition in the dense region at very low temperature [22].This phase transition is our primary target to investigate, employing the chiral condensate (cid:104) ¯ ψ ( x ) ψ ( x ) (cid:105) as an order parameter. Since the chiral symmetry plays a crucial role in thisstudy, we use the Kogut–Susskind fermion to formulate the NJL model on the lattice.The analysis of the phase structure with the TRG method would help us understand thethermodynamic properties of dense QCD.This paper is organized as follows. In Sec. 2 we explain the formulation of the latticeNJL model with the Kogut–Susskind fermion and the algorithmic details of the GrassmannATRG (GATRG). In Sec. 3, we compare the numerical results with the analytical ones inthe heavy dense limit as a benchmark before numerical results for the chiral condensateand the equation of state are presented. Section 4 is devoted to summary and outlook. We use the Kogut–Susskind fermion to formulate the NJL model on the lattice. FollowingRefs. [23, 24], we define the model at finite chemical potential µ as S = 12 a (cid:88) n ∈ Λ 4 (cid:88) ν =1 η ν ( n ) (cid:104) e µaδ ν, ¯ χ ( n ) χ ( n + ˆ ν ) − e − µaδ ν, ¯ χ ( n + ˆ ν ) χ ( n ) (cid:105) + ma (cid:88) n ∈ Λ ¯ χ ( n ) χ ( n ) − g a (cid:88) n ∈ Λ 4 (cid:88) ν =1 ¯ χ ( n ) χ ( n ) ¯ χ ( n + ˆ ν ) χ ( n + ˆ ν ) , (2.1)– 2 – 𝜇 Critical end point ത𝜓𝜓 ≠ 0 ത𝜓𝜓 = 0
Figure 1 . Schematic view of expected phase diagram of the NJL model on the T - µ plane. Solidand broken curves represent the first- and second-order phase transitions, respectively. Closed circledenotes the critical end point (CEP) where the first-order phase transition line terminates. where n = ( n , n , n , n )( ∈ Z ) specifies a position in the lattice Λ, with the lattice spacing a . χ ( n ) and ¯ χ ( n ) are Grassmann-valued fields without the Dirac structure. Since theydescribe the Kogut–Susskind fermions, χ ( n ) and ¯ χ ( n ) are single-component Grassmannvariables. η ν ( n ) is the staggered sign function defined by η ν ( n ) = ( − n + ··· + n ν − with η ( n ) = 1. The partition function is defined in the ordinal manner: Z = (cid:90) (cid:32) (cid:89) n ∈ Λ d χ ( n )d ¯ χ ( n ) (cid:33) e − S . (2.2)For vanishing mass m , Eq. (2.1) is invariant under the following continuous chiral trans-formation: χ ( n ) → e i α(cid:15) ( n ) χ ( n ) , (2.3)¯ χ ( n ) → ¯ χ ( n )e i α(cid:15) ( n ) (2.4)with α ∈ R and (cid:15) ( n ) = ( − n + n + n + n . We introduce the tensor network representation for Eq. (2.2) in a similar way with Refs. [10,11]. Hereafter, we set a = 1 for simplicity. Firstly, we expand the local Boltzmann weights See Ref. [25] for a different TRG approach with the Kogut-Susskind fermion, where the TRG procedureis applied to the Schwinger model after integrating out the fermion fields analytically. – 3 –n the following manners to decompose the nearest-neighbor interactions:exp (cid:20) − e µδ ν, η ν ( n ) ¯ χ ( n ) χ ( n + ˆ ν ) (cid:21) = (cid:88) i ν, ( n )=0 (cid:90)(cid:32) e µ δ ν, √ η ν ( n ) ¯ χ ( n )dΦ ν ( n ) · e µ δ ν, √ χ ( n + ˆ ν )d ¯Φ ν ( n + ˆ ν ) · ¯Φ ν ( n + ˆ ν )Φ ν ( n ) (cid:33) i ν, ( n ) , (2.5)exp (cid:20) e − µδ ν, η ν ( n ) ¯ χ ( n + ˆ ν ) χ ( n ) (cid:21) = (cid:88) i ν, ( n )=0 (cid:90)(cid:32) e − µ δ ν, √ η ν ( n ) χ ( n )dΨ ν ( n ) · e − µ δ ν, √ χ ( n + ˆ ν )d ¯Ψ ν ( n + ˆ ν ) · ¯Ψ ν ( n + ˆ ν )Ψ ν ( n ) (cid:33) i ν, ( n ) , (2.6)e g ¯ χ ( n ) χ ( n )¯ χ ( n +ˆ ν ) χ ( n +ˆ ν ) = (cid:88) i ν, ( n )=0 ( √ g ¯ χ ( n ) χ ( n ) · √ g ¯ χ ( n + ˆ ν ) χ ( n + ˆ ν )) i ν, ( n ) . (2.7)Secondly, integrating out χ and ¯ χ at each lattice site n , we define T n ; i ( n ) i ( n ) i ( n ) i ( n ) i ( n − ˆ4) i ( n − ˆ1) i ( n − ˆ2) i ( n − ˆ3) = (cid:90) d χ d ¯ χ e − m ¯ χχ (cid:89) ν =1 (cid:32) e µ δ ν, √ η ν ( n ) ¯ χ dΦ ν ( n ) (cid:33) i ν, ( n ) (cid:32) e µ δ ν, √ χ d ¯Φ ν ( n ) (cid:33) i ν, ( n − ˆ ν ) (cid:32) e − µ δ ν, √ η ν ( n ) χ dΨ ν ( n ) (cid:33) i ν, ( n ) (cid:32) e − µ δ ν, √ χ d ¯Ψ ν ( n ) (cid:33) i ν, ( n − ˆ ν ) ( √ g ¯ χχ ) i ν, ( n ) ( √ g ¯ χχ ) i ν, ( n − ˆ ν ) (cid:0) ¯Φ ν ( n + ˆ ν )Φ ν ( n ) (cid:1) i ν, ( n ) (cid:0) ¯Ψ ν ( n + ˆ ν )Ψ ν ( n ) (cid:1) i ν, ( n ) . (2.8)This serves as a change of variables from χ, ¯ χ to the integer-valued fields i ν = ( i ν,p ) p =1 , , and alternative Grassmann variables Φ ν , Ψ ν . Renaming x = i , y = i , z = i , t = i ,Eq. (2.2) is expressed in the form, Z = (cid:88) { t,x,y,z } (cid:90) (cid:89) n ∈ Λ T n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) , (2.9)which is the tensor network representation of this model. In current construction, T n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) is factorized as T n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) = I n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) S n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) G n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) . (2.10) In Eq. (2.9), we omit arguments in tensor indices and introduce shorthand notations such as x (cid:48) = x ( n − ˆ1) , y (cid:48) = y ( n − ˆ2) , z (cid:48) = z ( n − ˆ3) , t (cid:48) = t ( n − ˆ4) – 4 – n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) denotes the contributions from the integration over χ ( n ) and ¯ χ ( n ). A straight-forward calculation shows I n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) =( − n ( y + y + z + z + t + t )+ n ( z + z + t + t )+ n ( t + t ) (cid:18) √ (cid:19) t + t + x + x + y + y + z + z + t (cid:48) + t (cid:48) + x (cid:48) + x (cid:48) + y (cid:48) + y (cid:48) + z (cid:48) + z (cid:48) √ g t + x + y + z + t (cid:48) + x (cid:48) + y (cid:48) + z (cid:48) e µ ( t − t + t (cid:48) − t (cid:48) ) (cid:2) − m ¯∆ txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) , ∆ txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) , + ¯∆ txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) , ∆ txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) , (cid:3) , (2.11)where ¯∆ txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) ,q = δ t + t + x + x + y + y + z + z + t (cid:48) + t (cid:48) + x (cid:48) + x (cid:48) + y (cid:48) + y (cid:48) + z (cid:48) + z (cid:48) ,q , (2.12)∆ txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) ,q = δ t + t + x + x + y + y + z + z + t (cid:48) + t (cid:48) + x (cid:48) + x (cid:48) + y (cid:48) + y (cid:48) + z (cid:48) + z (cid:48) ,q , (2.13)with q = 0 ,
1. ¯∆ , ∆ are derived from ¯ χ -, χ -integration, respectively. The second line inEq. (2.11) comes from the staggered sign factor η ν ( n ). Consequently, I n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) doesdepend on n ∈ Λ. Eq. (2.11) tells us that this tensor network is uniform in t -direction, buthas some periodic structure in x -, y -, z -directions. This periodicity corresponds to the parityof the spatial lattice site n = ( n , n , n ). A graphical representation of Eq. (2.9) is shownin Fig. 2 (A). As a result of Eqs. (2.5) and (2.6), some Grassmann variables are allowedto exist in T n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) . These Grassmann variables have been denoted by G n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) in Eq. (2.10). Some sign can arise reflecting on how we have arranged these Grassmannvariables in G n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) and we have set this sign S n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) in Eq. (2.10). We nowassume that G n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) =dΦ t dΨ t dΦ x dΨ x dΦ y dΨ y dΦ z dΨ z d ¯Ψ t (cid:48) d ¯Φ t (cid:48) d ¯Ψ x (cid:48) d ¯Φ x (cid:48) d ¯Ψ y (cid:48) d ¯Φ y (cid:48) d ¯Ψ z (cid:48) d ¯Φ z (cid:48) (cid:0) ¯Φ ( n + ˆ4)Φ ( n ) (cid:1) t (cid:0) ¯Ψ ( n + ˆ4)Ψ ( n ) (cid:1) t (cid:0) ¯Φ ( n + ˆ1)Φ ( n ) (cid:1) x (cid:0) ¯Ψ ( n + ˆ1)Ψ ( n ) (cid:1) x (cid:0) ¯Φ ( n + ˆ2)Φ ( n ) (cid:1) y (cid:0) ¯Ψ ( n + ˆ2)Ψ ( n ) (cid:1) y (cid:0) ¯Φ ( n + ˆ3)Φ ( n ) (cid:1) z (cid:0) ¯Ψ ( n + ˆ3)Ψ ( n ) (cid:1) z , (2.14)where all the Grassmann measures depend on n and their arguments are omitted. Accord-ing to this arrangement, S n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) is given by S n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) = ( − t ( t + x + y + z )+ x ( x + y + z )+ y ( y + z )+ z z ( − t (cid:48) ( t (cid:48) + x (cid:48) + y (cid:48) + z (cid:48) )+ x (cid:48) ( x (cid:48) + y (cid:48) + z (cid:48) )+ y (cid:48) ( y (cid:48) + z (cid:48) )+ z (cid:48) z (cid:48) ( − ( t + t + x + x + y + y + z + z )( t (cid:48) + x (cid:48) + y (cid:48) + z (cid:48) ) . (2.15) We now formulate Grassmann ATRG (GATRG) algorithm to coarse grain the tensor net-work defined by Eq. (2.9). The basic idea is that we combine the ATRG procedure to– 5 – 𝑧 𝑦 𝑥 𝑡 GATRG
GATRGGATRG (A) (B)(C)(D) 𝑧 𝑦 𝑥𝑡 𝑡 𝑧 𝑦 𝑥𝑥𝑦𝑧𝑡 Figure 2 . Schematic illustration of the coarse-graining with the GATRG. Coordinate axes in Λare shown as dotted lines with arrows. Different numbers are assigned to specify different tensors.(A) Initial tensor network in Eq. (2.9). Eight types of tensor are located at a spatial unit cube. Aperiodic structure is explicitly shown on x - y plane. The tensor network is uniform in t -direction. (B)The first coarse-graining along z -direction reduces types of tensor from eight to four. The tensornetwork becomes uniform in t -, z -directions. (C) The second coarse-graining along y -directionmakes the structure uniform except in x -direction. (D) Totally uniform tensor network is obtainedby the third coarse-graining along x -direction. In the following coarse-graining steps, the structureis invariant. compress I n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) S n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) in Eq. (2.10) with the Grassmann HOTRG (GHOTRG)procedure to deal with G n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) in Eq. (2.10).Let us consider the coarse-graining along z -direction. Our aim is to construct a kind ofblock-spin transformation, T n +ˆ3 · T n (cid:55)→ T (cid:48) . We start from constructing the transformation G n +ˆ3 · G n (cid:55)→ G (cid:48) , employing the GHOTRG procedure demonstrated in Refs. [7, 17]. Astraightforward extension to the 4 d system gives us G (cid:48) ˜ t ˜ x ˜ yz ˜ t (cid:48) ˜ x (cid:48) ˜ y (cid:48) z (cid:48) = σ n +ˆ3 ,n d η ˜ t d ξ ˜ x d θ ˜ y dΦ z dΨ z d¯ η ˜ t (cid:48) d ¯ ξ ˜ x (cid:48) d¯ θ ˜ y (cid:48) d ¯Ψ z (cid:48) d ¯Φ z (cid:48) (¯ ηη ) ˜ t ( ¯ ξξ ) ˜ x (¯ θθ ) ˜ y (cid:0) ¯Φ ( n + ˆ3)Φ ( n ) (cid:1) z (cid:0) ¯Ψ ( n + ˆ3)Ψ ( n ) (cid:1) z . (2.16)One can confirm this expression by integrating out some of the Grassmann variables in G n +ˆ3 · G n with the help of some shifting technique. From now on, we call the exponent ofGrassmann number as the fermion index. In Eq. (2.16), new fermion indices, ˜ t, ˜ x, ˜ y, ˜ t (cid:48) , ˜ x (cid:48) , ˜ y (cid:48) ,are introduced with new Grassmann variables, η, ξ, θ, ¯ η, ¯ ξ, ¯ θ . The resulting sign factor See Ref. [7] for more details. – 6 – n +ˆ3 ,n in Eq. (2.16) is σ n +ˆ3 ,n = ( − z ( n )+ z ( n ) ( − [ t ( n +ˆ3)+ t ( n +ˆ3)][ t ( n )+ t ( n )]+[ t ( n +ˆ3)+ t ( n +ˆ3)+ x ( n +ˆ3)+ x ( n +ˆ3)][ x ( n )+ x ( n )] ( − [ t ( n +ˆ3)+ t ( n +ˆ3)+ x ( n +ˆ3)+ x ( n +ˆ3)+ y ( n +ˆ3)+ y ( n +ˆ3)][ y ( n )+ y ( n )] ( − [ y (cid:48) ( n +ˆ3)+ y (cid:48) ( n +ˆ3)][ y (cid:48) ( n )+ y (cid:48) ( n )]+[ x (cid:48) ( n +ˆ3)+ x (cid:48) ( n +ˆ3)+ y (cid:48) ( n +ˆ3)+ y (cid:48) ( n +ˆ3)][ x (cid:48) ( n )+ x (cid:48) ( n )] ( − [ t (cid:48) ( n +ˆ3)+ t (cid:48) ( n +ˆ3)+ x (cid:48) ( n +ˆ3)+ x (cid:48) ( n +ˆ3)+ y (cid:48) ( n +ˆ3)+ y (cid:48) ( n +ˆ3)][ t (cid:48) ( n )+ t (cid:48) ( n )] . (2.17)We now move on to the ATRG procedure for I n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) S n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) . In the following,we set T n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) = I n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) S n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) . Incorporating with σ n +ˆ3 ,n in Eq. (2.16), atotal block-spin transformation, T n +ˆ3 · T n (cid:55)→ T (cid:48) , is accomplished. In other words, what wehave to formulate is a block-spin transformation, σ n +ˆ3 ,n · T n +ˆ3 · T n (cid:55)→ T (cid:48) . The sign factor σ n +ˆ3 ,n consists of two types of index; ones are z , z living on the coarse-graining directionand the others are t , t , x , x , y , y living on the other directions. We deal with them ina separate way. Let σ (CG) n = ( − z ( n )+ z ( n ) and σ (NCG) n +ˆ3 ,n = σ n +ˆ3 ,n /σ (CG) n . Then σ (CG) n isincluded in the swapping bond part of the ATRG. This is quite natural because one has tocontract the tensors with respect to the index z in the swapping step of the ATRG. On theother hand, σ (NCG) n +ˆ3 ,n is handled both in finding squeezers and in contracting these squeezersand local tensors. Modifying these procedures in the ATRG, the GATRG is formulated.The coarse-graining along z -direction is followed by the series of coarse-graining along y -, x -, and t -directions. Fig. 2 shows a schematic picture of the first three times of coarse-graining. Though the original tensor network in Eq. (2.9) consists of several types of localtensor, the GATRG reduces them under a sequential coarse-graining process and we obtaina uniform tensor network in all directions.As a final supplement, we comment on another implementation for GATRG. It is alsopossible to coarse-grain G n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) following the philosophy of the ATRG. G n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) can be decomposed by introducing some extra Grassmann variables in the swapping bondpart, based on the same idea in the Grassmann TRG [6, 10]. This implementation mustreproduce the same result with the GATRG explained above if no finite bond dimension isintroduced. However, the results obtained with the finite bond calculation are possibly dif-ferent because these two GATRGs assume non-identical cost functions in optimization. Wehave numerically confirmed that this another GATRG also works, applying it to evaluatethe tensor network in Eq. (2.9). We have also found that the deviation between resultingthermodynamic potentials obtained by two types of GATRG tends to be smaller as thebond dimension is increased. Eq. (2.11) reveals that T takes a finite value if and only if its Grassmann parity is even.This feature can be understood as a kind of Z symmetry, which enables us to introduce theblock diagonal representation for some tensors treated in the GATRG algorithm. We carry See Refs.[4, 26] for more details about the swapping bond part and the squeezers used in the ATRGalgorithm. – 7 –ut the singular value decomposition (SVD) in the swapping bond part and the higher-order SVD in finding squeezers under the block diagonal representation for correspondingtensors. This blocking technique is of essential importance because it naturally defines thefermion indices introduced in Eq. (2.16). For instance, in order to find the squeezers incoarse-graining along ν -direction, we need to indirectly carry out the following SVD, Q ν ν ν ν ≈ D (cid:88) k =1 U ν ν ,k s k V † k,ν ν , (2.18)with D the bond dimension in the GATRG. In the block diagonal representation, this isexpressed by (cid:34) Q (even) Q (odd) (cid:35) ≈ (cid:34) U (even) U (odd) (cid:35) (cid:34) s (even) s (odd) (cid:35) (cid:34) V (even) † V (odd) † (cid:35) . (2.19)Let ˜ k be the fermion index for k . Then we assign ˜ k = 0(1) when s k comes from the matrix s (even) ( s (odd) ). In addition, the information of the fermion index ˜ k can be encoded in theordering of k : k = 1 , · · · , d (cid:124) (cid:123)(cid:122) (cid:125) ˜ k =0 , d + 1 , · · · , D (cid:124) (cid:123)(cid:122) (cid:125) ˜ k =1 . This means that D largest singular values in s consist of d singular values in s (even) and D − d singular values in s (odd) . A similar ordering trick is also available in the initial tensornetwork representation. Each index in the initial tensor T n ; txyzt (cid:48) x (cid:48) y (cid:48) z (cid:48) is composed of threeintegers, say i = ( i , i , i ), running from (0 , ,
0) to (1 , , i , i and i havebeen introduced via Eqs. (2.5), (2.6) and (2.7), respectively. Then the following mappingis practically useful:(0 , , (cid:55)→ , (1 , , (cid:55)→ , (0 , , (cid:55)→ , (1 , , (cid:55)→ , (1 , , (cid:55)→ , (0 , , (cid:55)→ , (1 , , (cid:55)→ , (0 , , (cid:55)→ . As we have seen in Eq. (2.14), the third component, say i , does not affect the Grassmannparity in the initial tensor. As a result, the parity of the sum of the first two components, i and i , corresponds to the Grassmann parity. Then the fermion index can be encodedin the ordering as 1 , , , (cid:124) (cid:123)(cid:122) (cid:125) (fermion index)=0 , , , , (cid:124) (cid:123)(cid:122) (cid:125) (fermion index)=1 . Thanks to this trick, Eq. (2.17) is simplified with the corresponding fermion indices as σ n +ˆ3 ,n = ( − ˜ z ( n )+˜ t ( n +ˆ3)˜ t ( n )+[˜ t ( n +ˆ3)+˜ x ( n +ˆ3)]˜ x ( n )+[˜ t ( n +ˆ3)+˜ x ( n +ˆ3)+˜ y ( n +ˆ3)]˜ y ( n ) ( − ˜ y (cid:48) ( n +ˆ3)˜ y (cid:48) ( n )+[˜ x (cid:48) ( n +ˆ3)+˜ y (cid:48) ( n +ˆ3)]˜ x (cid:48) ( n )+[˜ t (cid:48) ( n +ˆ3)+˜ x (cid:48) ( n +ˆ3)+˜ y (cid:48) ( n +ˆ3)]˜ t (cid:48) ( n ) . (2.20) Another ordering trick is demonstrated in Ref. [27]. – 8 –he last technique to be mentioned is the parallel computation, which reduces theexecution time of the GATRG. The essence of this technique in the ATRG is demonstratedin Ref. [28]; the computational cost per process of tensor contraction is reduced from O ( D )to O ( D ). As in Refs. [16, 28], we employ the randomized SVD (RSVD) in the swappingbond part. The accuracy of the RSVD is controlled by the oversampling parameter p and q iterations of QR decomposition. Under the block diagonal representation, we apply theRSVD with p = 4 D and q = D to each block matrix. We choose a large value of g = 32 for the four-fermi coupling in Eq. (2.1), because theFRG analysis in Ref. [21] indicates the vanishing phase transition for smaller g . Thepartition function of Eq. (2.9) is evaluated, using the GATRG algorithm on a lattice up tothe volume V = L ( L = 2 m , m ∈ N ). We assume the periodic boundary conditions for x -, y -, z -directions and the anti-periodic boundary condition for t -direction.Before investigating the restoration of the chiral symmetry at vanishing fermion mass,we check the efficiency of the GATRG algorithm by benchmarking with the NJL modelin the heavy dense limit, which is defined as m → ∞ and µ → ∞ , keeping e µ /m fixed.The heavy dense limit gives us an opportunity to compare numerical results with the exactanalytical ones. In the heavy dense limit, the number density (cid:104) n (cid:105) and the fermion condensate (cid:104) ¯ χ ( n ) χ ( n ) (cid:105) at vanishing temperature can be derived analytically as (cid:104) n (cid:105) = Θ( µ − µ c ) , (3.1) (cid:104) ¯ χ ( n ) χ ( n ) (cid:105) = 1 m Θ( µ c − µ ) , (3.2)where Θ denotes the step function and µ c = ln(2 m ) [29].Figures 3 and 4 show the numerical results for (cid:104) n (cid:105) and (cid:104) ¯ χ ( n ) χ ( n ) (cid:105) obtained by theGATRG algorithm choosing m = 10 with D = 30. The number density is calculated bythe numerical derivative of the thermodynamic potential in terms of the chemical potential: (cid:104) n (cid:105) = 1 V ∂ ln Z ( µ ) ∂µ ≈ V ln Z ( µ + ∆ µ ) − ln Z ( µ )∆ µ . (3.3)In the vicinity of µ c , we have set ∆ µ = 4 . × − . The fermion condensate is also obtainedvia the numerical derivative of the thermodynamic potential in terms of m : (cid:104) ¯ χ ( n ) χ ( n ) (cid:105)| m =10 = 1 V ln Z ( m + ∆ m ) − ln Z ( m )∆ m (cid:12)(cid:12)(cid:12)(cid:12) m =10 (3.4)– 9 –ith ∆ m = 1. Since there is little difference between the L = 128 and 1024 results, the L = 1024 lattice is sufficiently large to estimate the thermodynamic limit at vanishing tem-perature. The numerical results well reproduce the analytical ones, including the locationof µ c = ln(2 m ) = 9 . (cid:104) n (cid:105) and (cid:104) ¯ χ ( n ) χ ( n ) (cid:105) in the heavy dense limit. Note thatthe results quickly converge with respect to D ; the difference between ln Z ( D = 25) andln Z ( D = 30) has been already suppressed less than 2 . × − % in the vicinity of µ c . m -0.20.00.20.40.60.81.01.2 < n > Heavy dense limit at T = 0 L = 128 L = 1024 Figure 3 . Number density at m = 10 and g = 32 on 128 and 1024 lattices as a function of µ with D = 30. ∆ µ = 4 . × − in the vicinity of µ c . Green line denotes the step function inEq. (3.1). Having confirmed the efficiency of the GATRG algorithm in the heavy dense limit, letus turn to the calculation with the light fermion masses. We first check the convergencebehavior of the thermodynamic potential by defining the quantity δ = (cid:12)(cid:12)(cid:12)(cid:12) ln Z ( D ) − ln Z ( D = 55)ln Z ( D = 55) (cid:12)(cid:12)(cid:12)(cid:12) (3.5)on V = 1024 . In Fig. 5, we plot the D dependence of δ at µ = 2 . µ = 4 .
0, which is in the dense region with the restored chiralsymmetry, as we will see below. Although both of them are in the cold and dense regioncharacterized with µ/T ∼ O (10 ), where the Monte Carlo simulation should be severelyhindered by the sign problem, good convergent behaviors are observed in the GATRGcalculation; near the transition point, δ is reduced to about 10 − up to D = 55 and betterconvergence behavior, δ (cid:46) − , is observed at µ = 4 .
0. Hereafter we present the resultsat D = 55. – 10 – .70 9.75 9.80 9.85 9.90 9.95 10.00 m -2 · -5 · -5 · -5 · -5 · -5 · -4 · -4 < c-c > Heavy dense limit at T = 0 L = 128 L = 1024 Figure 4 . Fermion condensate at m = 10 and g = 32 on 128 and 1024 lattices as a function of µ with D = 30. Green line denotes the step function in Eq. (3.2).
25 30 35 40 45 50 55 D -7 -6 -5 -4 -3 -2 d m = 2.875 m = 4.0 Figure 5 . Convergence behavior of thermodynamic potential as a function of D on V = 1024 with m = 0 . We investigate the chiral phase transition employing the chiral condensate (cid:104) ¯ χ ( n ) χ ( n ) (cid:105) ,as an order parameter, which is defined by (cid:104) ¯ χ ( n ) χ ( n ) (cid:105) = lim m → lim V →∞ V ∂∂m ln Z, (3.6)– 11 –n the cold region. We calculate (cid:104) ¯ χ ( n ) χ ( n ) (cid:105) with the numerical derivative of thermodynamicpotential and the chiral extrapolation with the corresponding results at finite mass in thethermodynamic limit. In this study, the partial derivative in Eq. (3.6) is numericallyevaluated via ∂∂m ln Z ≈ ln Z ( m + ∆ m ) − ln Z ( m )∆ m , (3.7)with ∆ m = 0 .
01. In Fig. 6, we plot the µ dependence of the chiral condensate at m = 0 . L = 1024 lattice. The signals show slight fluctuations as a functionof µ around the transition point. Away from the transition point, we have found littleresponse in (cid:104) ¯ χ ( n ) χ ( n ) (cid:105) to changes in mass. Figure 7 presents the results in the chirallimit obtained by the chiral extrapolation with the data at m = 0 .
01 and 0 .
02 on twovolumes of L = 128 and 1024 . It is hard to find the difference between the L = 128and 1024 results. This allows us to consider the L = 1024 result to be essentially in thethermodynamic limit. We observe the discontinuity from a finite value to zero for thechiral condensate at µ c = 3 . ± . D is more essential than addingthe data points at different fermion masses in order to increase the numerical accuracyaround µ c found in Fig. 7. m -0.020.000.020.040.060.080.10 < c-c > m = 0.01 m = 0.02 Figure 6 . Chiral condensate at m = 0 .
01 and 0 .
02 on 1024 lattice as a function of µ with D = 55. It is possible to evaluate the chiral condensate with the impurity tensor method [17, 30]. Since Eq. (2.9)consists of eight types of tensor, there are eight configurations of an impurity tensor. Consequently, thecomputational cost is eight times larger than that of coarse-graining Eq. (2.9). One can also evaluate thenumber density discussed below with the impurity tensor method, which requires four times larger costthan that to coarse-grain Eq. (2.9) . – 12 – .0 1.0 2.0 3.0 4.0 5.0 m -0.020.000.020.040.060.080.10 < c-c > L = 128 L = 1024 Figure 7 . Chiral condensate extrapolated in the chiral limit as a function of µ with D = 55 on128 and 1024 lattices. The equation of state is a relation between the pressure and the particle number density.Here we present both results as functions of µ , respectively. In the thermodynamic limit,the pressure P is directly obtained from the thermodynamic potential: P = ln ZV , (3.8)where the vast homogeneous system is assumed. In Fig. 8, we plot the µ dependence of thepressure at m = 0 .
01. We find a kink behavior at µ c = 3 . ± . m = 0 .
02 result shows little differencefrom the m = 0 .
01 one.Fig. 9 shows the µ dependence of the particle number density (cid:104) n (cid:105) obtained by Eq. (3.3).We observe an abrupt jump from (cid:104) n (cid:105) = 0 to (cid:104) n (cid:105) = 1 at µ c = 2 . ± . µ c compared tothe cases of chiral condensate and pressure is attributed to the definition of the numericalderivative in Eq. (3.3). We have investigated the restoration of the chiral symmetry of the NJL model in thedense region at very low temperature, employing the Kogut–Susskind fermion action onthe extremely large lattice of V = 1024 , which is in the thermodynamic limit at zerotemperature, essentially. The first-order phase transition is clearly observed using the– 13 – .0 1.0 2.0 3.0 4.0 5.0 m P L = 128 L = 1024 Figure 8 . Pressure at m = 0 .
01 as a function of µ on 128 and 1024 lattices. m -0.20.00.20.40.60.81.01.2 < n > L = 128 L = 1024 Figure 9 . Particle number density at m = 0 .
01 as a function of µ on 128 and 1024 lattices. chiral condensate as an order parameter. At the critical chemical potential, we also findthe jump of the number density.This is the third successful application of the TRG method to the 4 d lattice theories,following the Ising model [31] and the complex φ theory at finite density [16]. This study isalso the first application to the 4 d fermionic system. As a next step, it would be interestingto determine the critical end point of this model.– 14 – cknowledgments Numerical calculation for the present work was carried out with the Oakforest-PACS (OFP)and the Cygnus computers under the Interdisciplinary Computational Science Program ofCenter for Computational Sciences, University of Tsukuba. This work is supported in partby Grants-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports,Science and Technology (MEXT) (No. 20H00148).
References [1] P. de Forcrand,
Simulating QCD at finite density , PoS
LAT2009 (2009) 010, [ ].[2] M. Levin and C. P. Nave,
Tensor renormalization group approach to two-dimensionalclassical lattice models , Phys. Rev. Lett. (2007) 120601, [ cond-mat/0611687 ].[3] Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang and T. Xiang, Coarse-grainingrenormalization by higher-order singular value decomposition , Phys. Rev. B (Jul, 2012)045139.[4] D. Adachi, T. Okubo and S. Todo, Anisotropic Tensor Renormalization Group , Phys. Rev. B (2020) 054432, [ ].[5] D. Kadoh and K. Nakayama,
Renormalization group on a triad network , .[6] Y. Shimizu and Y. Kuramashi, Grassmann tensor renormalization group approach toone-flavor lattice Schwinger model , Phys. Rev.
D90 (2014) 014508, [ ].[7] R. Sakai, S. Takeda and Y. Yoshimura,
Higher order tensor renormalization group forrelativistic fermion systems , PTEP (2017) 063B07, [ ].[8] Y. Shimizu and Y. Kuramashi,
Critical behavior of the lattice Schwinger model with atopological term at θ = π using the Grassmann tensor renormalization group , Phys. Rev.
D90 (2014) 074503, [ ].[9] Y. Shimizu and Y. Kuramashi,
Berezinskii-Kosterlitz-Thouless transition in lattice Schwingermodel with one flavor of Wilson fermion , Phys. Rev.
D97 (2018) 034502, [ ].[10] S. Takeda and Y. Yoshimura,
Grassmann tensor renormalization group for the one-flavorlattice Gross-Neveu model with finite chemical potential , PTEP (2015) 043B01,[ ].[11] D. Kadoh, Y. Kuramashi, Y. Nakamura, R. Sakai, S. Takeda and Y. Yoshimura,
Tensornetwork formulation for two-dimensional lattice N = 1 Wess-Zumino model , JHEP (2018) 141, [ ].[12] D. Kadoh, Y. Kuramashi, Y. Nakamura, R. Sakai, S. Takeda and Y. Yoshimura, Investigation of complex φ theory at finite density in two dimensions using TRG , JHEP (2020) 161, [ ].[13] Y. Kuramashi and Y. Yoshimura, Tensor renormalization group study of two-dimensionalU(1) lattice gauge theory with a θ term , JHEP (2020) 089, [ ].[14] A. Denbleyker, Y. Liu, Y. Meurice, M. P. Qin, T. Xiang, Z. Y. Xie et al., Controlling SignProblems in Spin Models Using Tensor Renormalization , Phys. Rev.
D89 (2014) 016008,[ ]. – 15 –
15] L.-P. Yang, Y. Liu, H. Zou, Z. Xie and Y. Meurice,
Fine structure of the entanglemententropy in the O(2) model , Phys. Rev. E (2016) 012138, [ ].[16] S. Akiyama, D. Kadoh, Y. Kuramashi, T. Yamashita and Y. Yoshimura, Tensorrenormalization group approach to four-dimensional complex φ theory at finite density , JHEP (2020) 177, [ ].[17] Y. Yoshimura, Y. Kuramashi, Y. Nakamura, S. Takeda and R. Sakai, Calculation offermionic Green functions with Grassmann higher-order tensor renormalization group , Phys.Rev.
D97 (2018) 054511, [ ].[18] Y. Nambu and G. Jona-Lasinio,
Dynamical Model of Elementary Particles Based on anAnalogy with Superconductivity. I , Phys. Rev. (1961) 345–358.[19] Y. Nambu and G. Jona-Lasinio,
Dynamical Model of Elementary Particles Based on anAnalogy with Superconductivity. II , Phys. Rev. (1961) 246–254.[20] M. Buballa,
NJL model analysis of quark matter at large density , other thesis, 2005.10.1016/j.physrep.2004.11.004.[21] K.-I. Aoki, S.-I. Kumamoto and M. Yamada,
Phase structure of NJL model with weakrenormalization group , Nucl. Phys. B (2018) 105–131, [ ].[22] M. Asakawa and K. Yazaki,
Chiral Restoration at Finite Density and Temperature , Nucl.Phys. A (1989) 668–684.[23] I.-H. Lee and R. E. Shrock,
Chiral Symmetry Breaking Phase Transition in Lattice GaugeHiggs Theories With Fermions , Phys. Rev. Lett. (1987) 14.[24] S. Booth, R. Kenway and B. Pendleton, The Phase Diagram of the Gauge InvariantNambu-Jona-Lasinio Model , Phys. Lett. B (1989) 115–120.[25] N. Butt, S. Catterall, Y. Meurice, R. Sakai and J. Unmuth-Yockey,
Tensor networkformulation of the massless schwinger model with staggered fermions , Phys. Rev. D (May, 2020) 094509.[26] H. Oba,
Cost Reduction of Swapping Bonds Part in Anisotropic Tensor RenormalizationGroup , PTEP (2020) 013B02, [ ].[27] S. Akiyama and D. Kadoh,
More about the Grassmann tensor renormalization group , .[28] S. Akiyama, Y. Kuramashi, T. Yamashita and Y. Yoshimura, Phase transition offour-dimensional Ising model with tensor network scheme , PoS
LATTICE2019 (2020) 138.[29] J. M. Pawlowski and C. Zielinski,
Thirring model at finite density in 2+1 dimensions withstochastic quantization , Phys. Rev. D (2013) 094509, [ ].[30] S. Morita and N. Kawashima, Calculation of higher-order moments by higher-order tensorrenormalization group , Computer Physics Communications (2019) 65 – 71.[31] S. Akiyama, Y. Kuramashi, T. Yamashita and Y. Yoshimura,
Phase transition offour-dimensional Ising model with higher-order tensor renormalization group , Phys. Rev.
D100 (2019) 054510, [ ].].