The Staggered Fermion for the Gross-Neveu Model at Non-zero Temperature and Density

The 2+1d Gross-Neveu model with finite density and finite temperature are studied by the staggered fermion discretization. The kinetic part of this staggered fermion in momentum space is used to build the relation between the staggered fermion and Wilson-like fermion. In the large Nf limit (the number Nf of staggered fermion flavors), the chiral condensate and fermion density are solved from the gap equation in momentum space, and thus the phase diagram of fermion coupling, temperature and chemical potential are obtained. Moreover, an analytic formula for the inverse of the staggered fermion matrix are given explicitly, which can be calculated easily by parallelization. The generalization to the 1+1d and 3+1d cases are also considered.


I. INTRODUCTION
The chiral phase transition in QCD from the hadronic phase at low temperature T (low density µ B ) to the quarkgluon plasma phase at high temperature (high density) has been studied intensively in the last decade. Although the relative firm statements for the phase structure can be made in two limit cases: finite T with small baryon density µ B ≪ T and asymptotically high density µ B ≫ Λ QCD , the phase structures at the intermediate baryon density are not clear. For a recent review of QCD with finite density, see Ref. [1][2] [3].
Since the chiral symmetry breaking and restoration are intrinsically non-perturbative, the number of techniques are limited and most results comes from the lattice QCD. Unfortunately the lattice QCD at finite density suffers from the notorious sign problem, especially for the intermediate or large baryon density. For some simpler quantum field models, e.g., the dense two-color QCD [4], the sign problem can be avoided. The recent progress of the sign problem in lattice field models can refer to [5] and references therein. This paper address a simplest four-fermion model with Z 2 symmetry: Gross-Neveu model at non-zero temperature and density [6] [7][8] [9] [10]. The 2+1d Gross-Neveu model has an interesting continuum limit and there is a critical coupling indicating the threshold for the symmetry breaking at zero temperature and density. Although the 2+1d Gross-Neveu model is not renormalisable in the weak coupling expansion, it is renormalisable in 1/N f expansion [6], where N f is the number of flavors of fermions.
The symmetry breaking of Gross-Neveu model for the 1+1d case has been studied extensively [11] [12][13] [14]. Recently the domain wall fermion was used to study the symmetry breaking for 2+1d Gross-Neveu model at zero temperature and density [15].
Compared with the Wilson fermion, the staggered fermion are more adequate for studying spontaneous chiral symmetry breaking. Another advantage of the staggered fermion is due to the reduced computational cost since the Dirac matrices have been replaced by the staggered phase factor. The reconstruction of the Wilson-like fermion from the staggered fermion is rather technique and thus need a careful explanation of the physical fermions for lattice QCD [16] and for Gross-Neveu model [8].
In this paper we revisit the staggered fermion for the 1+1d, 2+1d and 3+1d Gross-Neveu model at non-zero temperature and finite density. The gap equation, which is based on the large N f limit, are solved in the momentum space. Moreover we derive an explicit formula for the inverse matrix of the staggered fermion matrix, which is easy to be implemented by parallelization and thus make the large scale calculation of the gap equation feasible.
The arrangement of the paper is as follows. The continuum 2+1d Gross-Neveu model at finite density and non-zero temperature is introduced in section II. In section III the 2+1d staggered fermion is shown and non-dimensional quantities are introduced. The kinetic part of staggered fermion in the momentum space is given in section IV, where the trace of the inverse matrix and elements of inverse matrix are given explicitly in momentum space. In section V the results in section IV are generalized to the 1+1d and 3+1d staggered fermion matrix. Gap equation are given in section VI, where the chiral condensate and fermion density are calculated. The simulation results in the large N f limit are obtained in section VII. Finally the conclusion are given in section VIII.

II. THE GROSS-NEVEU MODEL
The Gross-Neveu (GN) model for interacting fermions in 2+1d is defined by the continuum Euclidiean Lagrangian density at finite density where ∂ / = 2 ν=0 γ ν ∂ ν ,μ is the chemical potential,m the bare mass, ψ andψ are an N f -flavor 4 component spinor fields. Here we choose the Gamma matrices where σ i (i = 1, 2, 3) are the Pauli matrices. The Gamma matrices satisfies There is a discrete Z 2 symmetry ψ → γ 5 ψ,ψ → −ψγ 5 , which is broken by the mass term but not the interaction. Introducing the bosonic filed σ, the interaction between fermions is decoupled, The dimension of quantities for the 2+1d GN model is as follows The partition function for this model is where ≡ β 0 dx 0 L 0 dx 1 dx 2 with the inverse temperature β = 1/T and the space size L.ψ and ψ are antiperiodic in x 0 direction, and are periodic in x 1 and x 2 directions. We want to calculate the chiral condensate for one flavor where V = βL 2 is the volume of 2+1d system. In the second equality we used Since the Lagrangian density is translation invariant, ψ (x)ψ(x) and σ(x) does not depend on x. This model in the large N f limit can be solved exactly [8] in the chiral limitm = 0, which is based on the saddle approximation (gap equation) in (6) where in the third equality we write the trace of operator in momentum space and the summation over k = (k 0 , k 1 , k 2 ) k 0 = (2n − 1)πT, k ν = 2n ν π/L, n, n ν ∈ Z, ν = 1, 2

III. THE STAGGERED FERMION
The staggered fermion discretization of the action L is with staggered phase factor η x,0 = 1, η x,1 = (−1) x0/a , η x,2 = (−1) (x0+x1)/a . aN x = L, a t N t = β = 1/T . The boundary condition for ψ andψ are accounted for by the sign s 1 and s 2 Here φ is defined on lattice x by σ(x) where [x,x] denotes 8 dual latticesx which is neighbour to x. The auxiliary field on dual lattice for two dimensional GN model was first studied in Ref. [18]. According to (5), the non-dimensional quantities are introduced by and thus the action in (9) can be rewritten as The partition function for the Gross-Neveu model with N f flavors: where ψ i andψ i denote the Grassmann fields of flavors i = 0, · · · , N f − 1 at the sites x, σ is the real field define at the dual lattice sitesx. The action is where The derivative of this matrix D with respect to the chemical potential and bare mass are rather simple The real matrix D(µ, m) satisfies the following symmetry where ε x = (−1) x0+x1+x2 is the parity of site x. By integrating the Grassmann fields, the partition function in (14) can be rewritten as with the effective action and The computational results, e.g., non-dimensional chiral condensate and fermion density, depend on the nondimensional quantities The physical dimensional quantities can be recovered from the non-dimensional ones by introducing lattice size a according to (12)(13). For notation simplicity, we set a t = a and thus a 1 = 1 in the following discussion.

IV. STAGGERED FERMION IN MOMENTUM SPACE
The kinetic part in (15) in one flavor is x,yχ (x)D x,y χ(y) where χ and χ are the Grassmann fields defined on lattices. A Wilson-like fermion can be obtained from the stagger fermion x,yχ (x)D x,y χ(y) [8]. Assume that N x and N t are even integers. Let Y = (Y 0 , Y 1 , Y 2 ) denotes a site on a lattice of twice the spacing of the original, and A = (A 0 , A 1 , A 2 ), A i = 0, 1 is a lattice vector, which ranges over the corners of the elementary cube associated with Y , so that each site on the original lattice x uniquely corresponds to A and Y : x = 2Y + A. Introducing notation A shift along µ direction can be represented by Similarly, χ(x) is defined on the fine lattice sites x with lattice size a = 1 while χ(A, ·) on the coarse lattice sites Y with lattice size 2a = 2 where 2 × 2 matrices Γ A and B A is given by Γ A and B A satisfies the following properties (The indices α, α ′ , β, a, a ′ and b always run from 1 to 2) Eq. (31) is also valid if Γ is replaced by B.
See A for these properties. Using (29), the inverse transformation of (25) and (26) are Let us introduce the two Dirac fields with 4 components (a = 1, 2) From the properties (30), it is easy to show that where in the last equality the inner produce betweenq and q is given in momentum space corresponding to the coarse lattice with lattice size 2 For any fixed µ = 0, 1, 2, where in the second equality (21) and (22) where we used the notations and the summation over k is taken for all modes in (36). Similarly, we have (B) and (37)(38), the kinetic part x,yχ (x)D x,y χ(y) can be rewritten as in the momentum space where the summation over k is taken for all momentum mode of coarse lattice according to (36), and the staggered matrix in the momentum space is diagonal where a µ and b c depends on k. The inverse matrix of D(k) is We can calculate the trace of inverse matrix D in (20) from (39) where the summation over k is given by (36). Note that the right hand side of (43) is real since k0 sin 2k 0 /|N (k)| 2 = 0 for any k 1 and k 2 modes in (36). Similarly, The inverse matrix of D in (20) is See C for the derivation of (44)(45)(46). Since D is diagonal in momentum space, the inverse matrix in theqq basis is where the notation with tilde denotes the inverse Fourier transformation, e.g., The analytic formula for the inverse matrix of the staggered fermion is the main contribution of this paper. Compared to the computational complexity O ((N t N 2 x ) 3 ) of the usual inverse matrix, the computational cost is O(16(N t N 2 x ) 2 ) since each element of the inverse matrix needs the summation over α, a, α ′ , a ′ = 1, 2. Moreover a parallel implementation can be realized easily for the formula (46).
The trace of the inverse matrix in (43) can be derived from (46) The staggered fermion matrix in (20) can be generalized to the 1+1d and 3+1d case, where α is 1 for the 1+1d case and α run from 1 to 3 for the 3+1d case.
For the 1+1d case, the 2 × 2 matrices γ µ are defined to be The unitary transformation in (25) and (26) are modified to be The kinetic part x,yχ (x)D x,y χ(y) can be written as where the summation is taken over all modes The fermion matrix in momentum space is diagonal with its inverse where N (k) = 4m 2 + µ=1 (sin 2k µ ) 2 − i cosh µ sin 2k 0 + sinh µ(cos 2k 0 + 1) The trace of the inverse matrix is The inverse matrix of D in can be calculated where For the 3+1d case, the 4 × 4 matrices γ µ are defined to be The unitary transformation in (25) and (26) are modified to be The kinetic part can also be written as (47) where the summation is taken for all modes Eq. (49)(50)(51) are still valid except that µ runs from 1 to 3. Eq. (52)(53)(54) are modified to be respectively. The details of staggered fermion matrix can be found in the supplement material. We have checked the formula (46)(53)(56) for the inverse matrices by Matlab.

VI. THE GAP EQUATION
The main contribution of the effective action (18) to the partition function can be obtained by the gap equation if N f → ∞, Here D is defined in (20) where m is replaced by m + Σ. The right hand side of (58) can be calculated from (42)(43) where m is replaced by m + Σ. The first derivative of Σ 2 with respect to µ can be computed from the gap equation (For simplicity, we assume that m = 0) If the average Σ of σ has been calculated from the gap equation, the free energy density in the large N f limit is where ln det D = k det D(k) up to a constant. The other thermodynamic quantities can be calculated. For example, the fermion density can be analytically calculated where ∂Σ 2 ∂µ , and two sums over x in (60) are given in (59) (44) and (45), respectively. The N (k) for each mode k in (44)(45)(59) is given by (42) with the replacement of m by m + Σ (Here for simplicity we assume that m = 0) and Σ is solved from the gap equation (58).

A. Large volume limit
Let us consider the large volume limit for the non-interacting theory. The partition function Z = dχdχe −χDχ = det D, where the stagger fermion matrix D is given by (20). The ratio of the non-dimensional chiral condensate a 2 ψ ψ and non-dimensional mass m = am is where in the last equality we used Equ. (43) where N (k), depending on m and µ, is given by (42). Note that there are N t N 2 x /8 modes k in (61). The ratio of the non-dimensional fermion density a 3 ρ and (aμ) 3 where in the last equality we used (44) and (45). We consider the case L = β, a = a t and thus N x = N t ≡ N . We fixμL andmL and then calculate a 2 ψ ψ am and a 3 ρ (aμ) 3 in the large N limit for fixed lattice size a. In fact a 2 ψ ψ am and a 3 ρ (aμ) 3 does not depend on the lattice size a since the non-dimensional mass m = am =m L N and non-dimensional chemical potential µ = aμ =μ L N does not depends on lattice size a. Figure 1 and shows the dependence of a 2 ψ ψ am on N with fixedμL,mL = 0, 1. The linear fitting with respect to 1/N shows that the large N limit of a 2 ψ ψ am is close to 1.008 for all four cases, this is because m = 1/N and µ = 1/N both vanish for large N limit. Figure 2 shows the dependence of a 3 ρ (aμ) 3 on N , whereμL = 1 andmL = 0, 1. The large N limit is close to 1.9271 for m = 0 and 1.9234 for m = 0.1/N , respectively.

B. Phase diagram
The phase diagram of the Gross-Neveu model in the large N f limit is well known [6][7] [8]. In this limit the phase diagram of (g −2 , µ, T ) is based on the calculation of Σ. Basically for T = 0 and µ = 0, there is a critical coupling g −2 c such that the chiral symmetry is broken Σ > 0 if the coupling is strong enough g −2 < g −2 c . This critical coupling depends on the regularization of the continuum model. For the lattice regularization in this paper, g −2 c ∼ a −1 where a is the lattice size. For fixed coupling g −2 < g −2 c which is not far away from the critical coupling (Otherwise, the continuum limit a → 0 can not be taken), denote Σ 0 be the value of Σ at this coupling g −2 with vanishing temperature T and chemical potential µ. The gap equation (8), which is solved exactly in the chiral limit in Ref. [8], shows that there exists a critical temperature T c = Σ0 2 ln 2 such that the chiral symmetry is broken if T < T c at this coupling g −2 and µ = 0. Moreover, there is another critical chemical potential µ c = Σ 0 such that this symmetry is broken only if µ < µ c at this coupling g −2 and T = 0. Moreover the mean field results predict that the first order transition only exists at T = 0 and µ c for this coupling.
We first study the dependence of Σ on the coupling g and temperature T = 1/N t without chemical potential µ = 0. Figure 3 is the phase diagram of (N t , 1/g 2 ) for m = 0 and N x = 36. The marks + separate the symmetry phase Σ = 0 (above marks) and the chiral symmetry broken phase Σ > 0 (below marks). For fixed temperature T there is a critical coupling g −2 c such that Σ decreases to zero if 1/g 2 is increasing to 1/g 2 c . Figure 3 shows that 1/g 2 c is a increasing function of N t = 1/T and it will close to 1 at very low temperature. On the other hand, if g −2 is fixed, there is a critical temperature T c = T c (g) such that Σ is increasing from zero if T is decreasing from T c . Figure 4 shows the dependence of Σ on N t for the different coupling 1/g 2 . For small 1/g 2 , for example, 1/g 2 = 0.65, Σ changes small with the temperature. For these range of parameters, it is in the deep chiral symmetry broken phase and we can not obtain the chiral symmetry phase Σ = 0 even at very high temperature. For a slightly larger 1/g 2 , for example, 1/g 2 = 0.90 (black dots in Figure 4), we can find a transition point T c , which is between 1 8 and 1 10 in lattice unit. The symmetry phase and broken phase are realized for T > T c (g) and T < T c (g), respectively. Figure 5 shows the dependence of Σ on 1/g 2 at different temperature. Σ drops continuously to 0 if 1/g 2 is increasing to 1/g 2 c (T ) from below, which show that the transition at the critical coupling constant g c (T ) is second order. At very low temperature T = 1/N t = 1/36, g c (T ) is close to 1, which is consistent with those obtained in [9]. This is because in the limit of N t , N x → ∞, the gap equation at Σ = 0 is reduced to The critical temperature T c = Σ0 2 ln 2 at the coupling g −2 and µ = 0 can be verified numerically. Here we choose N x = 36 and g −2 = 0.95 which is not too far away from the critical coupling g −2 c ≈ 1. We also choose N t = 36 such that it is very close to zero temperature, the value of Σ at the zero temperature and vanishing chemical potential is Σ 0 = 0.0944. To calculate the critical temperature at this coupling, we calculate Σ at N t = 8, · · · , 36 and found that Σ is zero if N t is between 14 and 16. Thus the critial temperature is between 1/16 = 0.0625 and 1/14 = 0.0667 which is very close to T c = Σ0 2 ln 2 = 0.0944 2 ln 2 = 0.0680. Now let us study the effect of chemical potential on the chiral condensate Σ. Figure 6 shows the dependence of Σ on the chemical potential at the different temperature T = 1/N t . Σ drops sharply near µ c ≈ 0.45 in the limit of zero temperature N t = 16, i.e., T = 1/16, which suggest a first order transition at the zero temperature. This first order transition at the zero temperature is verified by the analytical calculation, µ c = Σ 0 where Σ 0 is the Σ with µ = 0 [8]. For the temperature T = 1/16, Σ 0 ≈ 0.47 is slightly larger than µ c ≈ 0.45. If the temperature is raised, e.g., N t = 6, it is more difficult to find a critical chemical potential such that the chiral symmetry is restored. This is not caused by the smallness of N x = 36, since the our results is always obtained for N x = 36, which is very close to the thermodynamics limit, i.e., the result changes very small if N x is larger than 36. We also note that the transition at finite temperature is the second order, as explained in [8]. Figure 7 shows the dependence of Σ on µ for a larger 1/g 2 = 0.80. Compared with Figure 6, Σ at µ = 0 and the critical chemical potential in Figure 7 become smaller, and thus the figures in Figure 7 is obtained by moving those figures of Figure 6 in the left-down direction. For the same temperature, for example, N t = 16, it is more difficult to find the critical chemical potential in Figure 7 than those in Figure 6. Both Figure 6 and Figure 7 show that the critical chemical potential µ c is decreased if the temperature is increased. At zero temperature, the mean field exact result show the critical chemical potential µ c is just the value of Σ 0 at the vanishing chemical potential. This is exactly recovered in Figure 7 where µ c = 0.32 for g −2 = 0.80 with N t = 16. Figure 8 shows the dependence of Σ and fermion density on the chemical potential at 1/g 2 = 0.7. At low temperature N t = 16, Σ drops rapidly near the critical chemical potential µ c ≈ 0.45, and the fermion density increase very fast, which suggest Σ and fermion density are not continuous at µ c at zero temperature and thus they can be regarded as the order parameters.
For finite N f = 2, 6, 12, the hybrid Monte Carlo method was used to study the symmetry breaking and symmetry restoration [8] [10]. The Monte Carlo results for N f = 12 is consistent with those obtained from the correction of mean field theory to order 1/N 2 f . The staggered fermion of Gross-Neveu model with N f flavors is equivalent to the Wilson-like fermion (4 Dirac components) with 2N f flavors. To avoid the sign problem in the hybrid Monte Carlo method, N f must be even.
For the 3+1d Gross-Neveu model, we also calculate the dependence of Σ on the coupling and chemical potential at different temperature. Figure 9 shows the value of Σ depending on the coupling for the vanishing chemical potential. Compared to Figure 5 for the 2+1d model, the critical coupling becomes smaller. Moreover, the dependence of Σ on the temperature is less sensitive. Figure 10 shows the dependence of Σ on the chemical potential at the coupling 1/g 2 = 0.58 for the 2+1d and 3+1d Gross-Neveu model, the critical chemical potential is smaller for the 2+1d model than those for the 3+1d model.

VIII. CONCLUSION
The staggered fermion for the Gross-Neveu model at finite density and temperature are revisited. In the large N f limit, this model in 1+1d, 2+1d and 3+1d dimension can be easily solved in momentum space. Moreover an explicit formula for the inverse matrix for the 1+1d, 2+1d and 3+1d staggered fermion matrix are found, which can be implemented by parallelization. This formula can also generalized to the other space dimensions. For the odd space dimension, the orthogonal transformation were found [17]. The key point to find the explicit formula for the inverse matrix is to use the properties of Γ A and B A as shown in Section IV. These properties for the even number of space dimension is more simple, as shown in the supplement material. The dependence of chiral condensate and fermion density on the coupling, temperature and chemical potential are obtained by solving the gap equation. Our results for the 2+1d case reproduce the analytical results. We also compared the chiral condensate for the 2+1d and 3+1d case in the same range of parameters and shows that the reason for symmetry breaking and restoration can be explained by the suitable choice of the coupling, temperature and chemical potential.
which is also valid if (1, 2) is replaced by (1,3) or (2,3). Secondly, To prove that we want to prove that