Restoration of chiral symmetry in cold and dense Nambu--Jona-Lasinio model with tensor renormalization group

We analyze the chiral phase transition of the Nambu--Jona-Lasinio model in the cold and dense region on the lattice developing the Grassmann version of the anisotropic tensor renormalization group algorithm. The model is formulated with the Kogut--Susskind fermion action. We use the chiral condensate as an order parameter to investigate the restoration of the chiral symmetry. The first-order chiral phase transition is clearly observed in the dense region at vanishing temperature with $\mu/T\sim O(10^3)$ on a large volume of $V=1024^4$. We also present the results for the equation of state.


Introduction
The phase structure and the equation of state for QCD at finite temperature and density are essential ingredients to understand the evolution and the current state of the universe quantitatively. Although the lattice QCD simulation has been expected to be an ideal tool to investigate the nonperturbative aspects of QCD, it has not been successful to reveal the nature of QCD at finite density. This is due to the sign problem caused by the introduction of the chemical potential in lattice QCD simulations based on the Monte Carlo algorithm, see, e.g., Ref. [1].
The tensor renormalization group (TRG) method, which was originally proposed by Levin and Nave to study two-dimensional (2d) classical spin models in 2007 [2], has several superior features over the Monte Carlo method 1 . Firstly, the TRG method is free from the sign problem. This virtue was confirmed by successful application of the TRG method to various 2d quantum field theories which contain the sign problem [6,[8][9][10][11][12][13]. Moreover, the authors have successfully employed the anisotropic TRG (ATRG) algorithm to analyze the Bose condensation in the 4d complex φ 4 theory at finite density [14]. Secondly, the computational cost depends on the system size only logarithmically. This is an essential ingredient to enable us to take the thermodynamic limit at vanishing temperature. Thirdly, the Grassmann version of the TRG method, which was developed by some of the authors [6,7,15], allows direct manipulation of the Grassmann variables. Fourthly, we can obtain 1 In this paper the TRG method or the TRG approach refers to not only the original numerical algorithm proposed by Levin and Nave but also its extensions [3][4][5][6][7]. the partition function or the path-integral itself. In the thermodynamic limit the pressure is directly related to the thermodynamic potential 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 [16,17] at finite temperature T and chemical potential µ on the lattice developing the Grassmann version of the ATRG algorithm. The Lagrangian of the NJL model in the continuum is defined as follows: which has the U(1) chiral symmetry with ψ(x) → e iαγ 5 ψ(x) andψ(x) →ψ(x)e iαγ 5 . This is an effective theory of QCD which describes the dynamical chiral symmetry breaking: once the strength of the coupling constant g 0 exceeds a certain critical value the system generates a non-trivial vacuum with ψ (x)ψ(x) = 0. The chiral phase structure of the NJL model on the T -µ plane is discussed by some analytical methods, e.g., the mean-field approximation (MFA) [18] and the functional renormalization group (FRG) [19]. Figure 1 shows a schematic view of the expected phase structure, whose characteristic feature is the first-order chiral phase transition in the dense region at very low temperature [20]. This phase transition is our primary target to investigate, employing the chiral condensate ψ (x)ψ(x) as an order parameter. Since the chiral symmetry plays a crucial role in this study, 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 the thermodynamic properties of dense QCD. This paper is organized as follows. In Sec. 2 we explain the formulation of the lattice NJL model with the Kogut-Susskind fermion and the algorithmic details of the Grassmann ATRG (GATRG). Numerical results for the chiral condensate and the equation of state are presented in Sec. 3. Section 4 is devoted to summary and outlook.

Tensor network representation
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 in the following manners to decompose the nearest-neighbor interactions: This serves as a change of variables from χ,χ to the integer-valued fields i ν = (i ν,p ) p=1,2,3 and alternative Grassmann variables Φ ν , Ψ ν . Renaming which is the tensor network representation of this model 2 . In current construction, T n;txyzt x y z is factorized as T n;txyzt x y z = I n;txyzt x y z S n;txyzt x y z G n;txyzt x y z .

Procedure of the algorithm
We now formulate Grassmann ATRG (GATRG) algorithm to coarse-grain the tensor network defined by Eq. (2.9). The basic idea is that we combine the ATRG procedure to compress I n;txyzt x y z S n;txyzt x y z in Eq. (2.10) with the Grassmann HOTRG (GHOTRG) procedure to coarse-grain G n;txyzt x y z in Eq. (2.10). Let us consider the coarse-graining along z-direction. Our aim is to construct a kind of block-spin transformation, T n+3 · T n → T . We start from constructing the transformation G n+3 · G n → G , employing the GHOTRG procedure demonstrated in Refs. [7,15]. A straightforward extension to the 4d system gives us The first coarse-graining along z-direction reduces numbers of tensor from eight to four. The tensor network becomes uniform in t-, z-directions. (C) The second coarse-graining along y-direction makes the structure uniform except in x-direction. (D) Totally uniform tensor network is obtained by the third coarse-graining along x-direction. In the following coarse-graining steps, the structure is invariant.
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 3 . From now on, we call the expotent of Grassmann number as the fermion index. In Eq. (2.16), new fermion indices,t,x,ỹ,t ,x ,ỹ , are introduced with new Grassmann variables, η, ξ, θ,η,ξ,θ. The resulting sign factor σ n+3,n in Eq. (2.16) is (2.17) We now move on to the ATRG procedure for I n;txyzt x y z S n;txyzt x y z . In the following, we set T n;txyzt x y z = I n;txyzt x y z S n;txyzt x y z . Incorporating with σ n+3,n in Eq. (2.16), a total block-spin transformation, T n+3 · T n → T , is accomplished. In other words, what we have to formulate is a block-spin transformation, σ n+3,n · T n+3 · T n → T . The sign factor σ n+3,n consists of two types of index; ones are z 1 , z 2 living on the coarse-graining direction and the others are t 1 , t 2 , x 1 , x 2 , y 1 , y 2 living on the other directions. We deal with them in a separate way. Let σ (CG) n = (−1) z 1 (n)+z 2 (n) and σ (NCG) n+3,n = σ n+3,n /σ (CG) n . Then σ (CG) n is included in swapping bond part in the ATRG 4 . This is quite natural because one has to contract the tensors with respect to the index z in swapping step in the ATRG. On the other hand, σ (NCG) n+3,n is handled both in finding squeezers and in contracting these squeezers and local tensors. Modifying these procedures in the ATRG, the GATRG is finally 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 coarsegraining. Though the original tensor network in Eq. (2.9) consists of several types of local tensor, the GATRG reduces them under a sequential coarse-graining process and finally we obtain a uniform tensor network in all directions.
As a final supplement, we comment on another implementation for GATRG. It is also possible to coarse-grain G n;txyzt x y z following the philosophy of the ATRG. G n;txyzt x y z can be decomposed by introducing some extra Grassmann variables in swapping bond part, based on the same idea in the Grassmann TRG [6,10]. This implementation must reproduce the same result with the GATRG explained above if no finite bond dimension is introduced. However, the results obtained with the finite bond calculation are possibly different because these two GATRGs assume non-identical cost functions in optimization. We have numerically confirmed that this another GATRG also works, applying it to evaluate the tensor network in Eq. (2.9). We have also found that the deviation between resulting thermodynamic potentials obtained by two types of GATRG tends to be smaller as the bond dimension is increased.

Some techniques
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 2 symmetry, which enables us to introduce the block diagonal representation for some tensors treated in the GATRG algorithm. We carry out the singular value decomposition (SVD) in swapping bond part and the higherorder SVD in finding squeezers under the block diagonal representation for corresponding tensors. This blocking technique is of essential importance because it naturally defines the fermion indices introduced in Eq. (2.16). For instance, in order to find the squeezers in coarse-graining along ν-direction, we need to indirectly carry out the following SVD, with D the bond dimension in the GATRG. In the block diagonal representation, this is expressed by Letk be the fermion index for k. Then we assignk = 0(1) when s k comes from the matrix s (even) (s (odd) ). In addition, the information of the fermion indexk can be encoded in the ordering of k: 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 tensor network representation 5 . Each index in the initial tensor T n;txyzt x y z is composed of three integers, say i = (i 1 , i 2 , i 3 ), running from (0, 0, 0) to (1, 1, 1). Note that i 1 , i 2 and i 3 have been introduced via Eqs. (2.5), (2.6) and (2.7), respectively. Then the following mapping is practically useful: (1, 0, 0) → 5, (0, 1, 0) → 6, (1, 0, 1) → 7, (0, 1, 1) → 8.
As we have seen in Eq. (2.14), the third component, say i 3 , does not affect the Grassmann parity in the initial tensor. As a result, the parity of the sum of the first two components, i 1 and i 2 , corresponds to the Grassmann parity. Then the fermion index can be encoded in the ordering as .
Thanks to this trick, Eq. (2.17) is simplified with the corresponding fermion indices as σ n+3,n = (−1)z (n)+t(n+3)t(n)+[t(n+3)+x(n+3)]x(n)+[t(n+3)+x(n+3)+ỹ(n+3)]ỹ(n) The last technique to be mentioned is the parallel computation, which reduces the execution time of the GATRG. The essence of this technique in the ATRG is demonstrated in Ref. [25]; the computational cost per process of tensor contraction is reduced from O(D 9 ) to O(D 8 ). As in Refs. [14,25], we employ the randomized SVD (RSVD) in swapping bond 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 the RSVD with p = 4D and q = D to each block matrix.

Setup
We choose a large value of g 0 = 32 for the four-fermi coupling in Eq. (2.1), because the FRG analysis in Ref. [19] indicates the vanishing phase transition for smaller g 0 . The partition function of Eq. (2.9) is evaluated using the GATRG algorithm on a lattice up to the volume V = L 4 (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. To investigate the convergence behavior of the thermodynamic potential, we define the quantity,

Chiral phase transition
We investigate the chiral phase transition employing the chiral condensate χ(n)χ(n) , as an order parameter, which is defined by in the cold region. We calculate χ(n)χ(n) with the numerical derivative of thermodynamic potential and the chiral extrapolation with the corresponding results at finite mass in the thermodynamic limit 6 . In this study, the partial derivative in Eq. (3.2) is numerically evaluated via with ∆m = 0.01. In Fig. 4 we plot the µ dependence of the chiral condensate at m = 0.01 and 0.02 on the L 4 = 1024 4 lattice. The signals show slight fluctuations as a function of µ around the transition point. Away from the transition point, we have found little response in χ(n)χ(n) to changes in mass. Figure 5 presents the results in the chiral limit obtained by the chiral extrapolation with those at m = 0.01 and 0.02 on two volumes of L 4 = 128 4 and 1024 4 . It is hard to find the difference between the L = 128 and 1024 results. This allows us to consider the L = 1024 result to be essentially in the thermodynamic limit. We observe the discontinuity from a finite value to zero for the chiral condensate at µ c = 3.0625 ± 0.0625, which is a clear indication of the first-order phase transition.

Equation of state
Equation of state is a relation between the pressure and the particle number density. Here we presents both results as a function of µ, respectively. In the thermodynamic limit, the pressure P is directly obtained from the thermodynamic potential: where the vast homogeneous system is assumed. In Fig. 6 we plot the µ dependence of the pressure at m = 0.01. We find a kink behavior at µ c = 3.0625 ± 0.0625, where the chiral condensate shows the discontinuity. Note that the m = 0.02 result shows little difference from the m = 0.01 one. The particle number density is obtained by the numerical derivative of pressure in terms of the chemical potential: The µ dependence of n is shown in Fig. 7. We observe an abrupt jump from n = 0 to n = 1 at µ c = 2.9375 ± 0.0625. This is another indication of the first-order phase transition. The small shift of µ c compared to the chiral condensate case is attributed to the definition of the numerical derivative in Eq. (3.5).   with the impurity tensor, following the ideas in Refs. [15,26]. This is also the case with the evaluation of number density discussed below.

Summary and outlook
We have investigated the restoration of the chiral symmetry of the NJL model in the dense region at very low temperature employing the Kogut-Susskind fermion action on the extremely large lattice of V = 1024 4 , which is essentially in the thermodynamic limit at zero temperature. The first-order phase transition is clearly observed using the chiral condensate as an order parameter. At the critical chemical potential we also find the jump of the number density. This is the third successful application of the TRG method to the 4d lattice theories following the Ising model [27] and the complex φ 4 theory at finite density [14]. This study is also the first application to the 4d fermionic system. As a next step it would be interesting to determine the critical end point.