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 μ/T ∼ O(103) on a large volume of V = 10244. 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 non-perturbative 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 the 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, and also other tensor network methods, are free from the sign problem. This virtue was confirmed by successful application of these methods to various 2d quantum field theories which contain the sign problem [6,[8][9][10][11][12][13]. 2 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 [16]. Secondly, the computational cost depends on the system size 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]. 2 The TRG study of the sign problem was also performed by introducing the complex coupling or the chemical potential to the 2d classical O(2) spin model [14,15].

JHEP01(2021)121
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,17], allows direct manipulation of the Grassmann variables. Fourthly, we can obtain 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 [18,19] 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) [20] and the functional renormalization group (FRG) [21]. 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 [22]. 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 section 2 we explain the formulation of the lattice NJL model with the Kogut-Susskind fermion and the algorithmic details of the Grassmann ATRG (GATRG). In section 3, we compare the numerical results with the analytical ones in the heavy dense limit as a benchmark before numerical results for the chiral condensate and the equation of state are presented. 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]. 3 Hereafter, we set a = 1 for simplicity. Firstly, we expand the local Boltzmann weights

JHEP01(2021)121
in the following manners to decompose the nearest-neighbor interactions: exp − e µδ ν, 4 2 η ν (n)χ(n)χ(n +ν) Secondly, integrating out χ andχ at each lattice site n, we define This serves as a change of variables from χ,χ to the integer-valued fields i ν = (i ν,p ) p=1,2,3 and alternative Grassmann variables which is the tensor network representation of this model. 4 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 . (2.10)

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 JHEP01(2021)121 compress I n;txyzt x y z S n;txyzt x y z in eq. (2.10) with the Grassmann HOTRG (GHOTRG) procedure to deal with 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,17]. A straightforward extension to the 4d system gives us 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 σ 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 formulated.
The coarse-graining along z-direction is followed by the series of coarse-graining along y-, x-, and t-directions. Figure 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 local tensor, the GATRG reduces them under a sequential coarse-graining process and 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 the 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 JHEP01(2021)121 out the singular value decomposition (SVD) in the 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: . (1, 0, 0) → 5, (0, 1, 0) → 6, (1, 0, 1) → 7, (0, 1, 1) → 8.

This means that
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 (2.20)

JHEP01(2021)121
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. [28]; the computational cost per process of tensor contraction is reduced from O(D 9 ) to O(D 8 ). As in refs. [16,28], we employ the randomized SVD (RSVD) in the 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. [21] 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.
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 model in 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 exact analytical ones.

Heavy dense limit as a benchmark
In the heavy dense limit, the number density n and the fermion condensate χ(n)χ(n) at vanishing temperature can be derived analytically as where Θ denotes the step function and µ c = ln(2m) [29]. Figures 3 and 4 show the numerical results for n and χ(n)χ(n) obtained by the GATRG algorithm choosing m = 10 4 with D = 30. The number density is calculated by the numerical derivative of the thermodynamic potential in terms of the chemical potential: In the vicinity of µ c , we have set ∆µ = 4.0 × 10 −3 . The fermion condensate is also obtained via the numerical derivative of the thermodynamic potential in terms of m:  of µ c = ln(2m) = 9.903, both for n and χ(n)χ(n) in the heavy dense limit. Note that the results quickly converge with respect to D; the difference between ln Z(D = 25) and ln Z(D = 30) has been already suppressed less than 2.1 × 10 −3 % in the vicinity of µ c .

Chiral phase transition
Having confirmed the efficiency of the GATRG algorithm in the heavy dense limit, let us turn to the calculation with the light fermion masses. We first check the convergence 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. 8 In this study, the partial derivative in eq.

evaluated via
with ∆m = 0.01. In figure 6, 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 7 presents the results in the chiral limit obtained by the chiral extrapolation with the data 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. Note that enlarging the bond dimension D is more essential than adding the data points at different fermion masses in order to increase the numerical accuracy around µ c found in figure 7.

Equation of state
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 Z V , (3.8) where the vast homogeneous system is assumed. In figure 8, we plot the µ dependence of the pressure at m = 0.01. We find a kink behavior at µ c = 3. condensate shows the discontinuity. Note that the m = 0.02 result shows little difference from the m = 0.01 one. Figure 9 shows the µ dependence of the particle number density n obtained by eq. (3.3). 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 cases of chiral condensate and pressure is attributed to the definition of the numerical derivative in eq.

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 in the thermodynamic limit at zero temperature, essentially. 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 [31] and the complex φ 4 theory at finite density [16]. 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 of this model.