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

Daming Li

School of Mathematical Sciences, Shanghai Jiao Tong University, Shanghai, China.

DOI: 10.4236/jmp.2021.1213105                

Abstract

The 2 + 1d Gross-Neveu model with finite density and finite temperature is 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 N_f limit (the number N_f 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 is obtained. Moreover, an analytic formula for the inverse of the staggered fermion matrix is given explicitly, which can be calculated easily by parallelization. The generalization to the 1 + 1d and 3 + 1d cases is also considered.

Keywords

Gross-Neveu Model, Phase Diagram, Staggered Fermion, Gap Equation

1. Introduction

The chiral phase transition in quantum chromodynamics (QCD) from the hadronic phase at low temperature T (low density ${\mu }_B$ ) to the quark-gluon 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 ${\mu }_B\ll T$ and asymptotically high density ${\mu }_B\gg {\Lambda }_{\text{QCD}}$, the phase structures at the intermediate baryon density are not clear. For a recent and review and related work of QCD with finite density, see Ref. [1] - [9].

Since the chiral symmetry breaking and restoration are intrinsically non-perturbative, the number of techniques is limited and most results come 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 [10], the sign problem can be avoided. The recent progress of the sign problem in lattice field models can refer to [11] and references therein. In the last decades, the tensor network becomes very popular in condensed matter physics and high energy physics, especial for lower dimension models, since probability is not used and thus it is free of sign problem [12] [13] [14] [15].

This paper addresses a simplest four-fermion model with $Z_2$ symmetry: Gross-Neveu model at non-zero temperature and density [16] [17] [18] [19] [20]. 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 [16], 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 [21] - [29]. Recently, 2 + 1d Gross-Neveu model is used to study the inhomogeneous phases [30] and the symmetry breaking [31].

Compared with the Wilson fermion, the staggered fermion is 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, thus needing a careful explanation of the physical fermions for lattice QCD [32] and for Gross-Neveu model [18].

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, is 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 2. In Section 3, 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 4, where the trace of the inverse matrix and elements of inverse matrix are given explicitly in momentum space. In Section 5, the results in Section 4 are generalized to the 1 + 1d and 3 + 1d staggered fermion. The gap equation is given in Section 6, where the chiral condensate and fermion density are calculated. The simulation results in the large $N_f$ limit are obtained in Section 7. Finally, the conclusion is given in Section 8.

2. The Gross-Neveu Model

The Gross-Neveu model for interacting fermions in 2 + 1d is defined by the continuum Euclidiean Lagrangian density at finite density

$\mathcal{L}=\bar{\psi }\left(\cancel{\partial }+\tilde{\mu }{\gamma }_0+\tilde{m}\right)\psi -\frac{{\tilde{g}}^2}{2N_f}{\left(\bar{\psi }\psi \right)}^2$ (1)

where $\cancel{\partial }={\displaystyle {\sum }_{\nu =0}^2}{\gamma }_{\nu }{\partial }_{\nu }$, $\tilde{\mu }$ is the chemical potential, $\tilde{m}$ the bare mass, $\psi$ and $\bar{\psi }$ are an $N_f$ -flavor 4 component spinor fields. Here we choose the Gamma matrices

${\gamma }_{\nu }=\left(\begin{array}{cc}{\sigma }_{\nu +1}& 0\\ 0& -{\sigma }_{\nu +1}\end{array}\right),\nu =0,1,2$ (2)

${\gamma }_3=\left(\begin{array}{cc}& -i{\mathbb{I}}_2\\ i{\mathbb{I}}_2& \end{array}\right),{\gamma }_5={\gamma }_0{\gamma }_1{\gamma }_2{\gamma }_3=\left(\begin{array}{cc}& {\mathbb{I}}_2\\ {\mathbb{I}}_2& \end{array}\right)$ (3)

where ${\sigma }_i\left(i=1,2,3\right)$ are the Pauli matrices. The Gamma matrices satisfies

${\gamma }_{\mu }{\gamma }_{\nu }+{\gamma }_{\nu }{\gamma }_{\mu }={\delta }_{\mu \nu }2{\mathbb{I}}_4,\mu ,\nu =0,1,2,3,5$

There is a discrete $Z_2$ symmetry $\psi \to {\gamma }_5\psi$, $\bar{\psi }\to -\bar{\psi }{\gamma }_5$, which is broken by the mass term but not the interaction. Introducing the bosonic field $\sigma$, the interaction between fermions is decoupled with the Lagrangian density,

$L=\bar{\psi }\left(\cancel{\partial }+\tilde{\mu }{\gamma }_0+\tilde{m}+\sigma \right)\psi +\frac{N_f}{2{\tilde{g}}^2}{\sigma }^2$ (4)

The dimension of quantities for the 2 + 1d Gross-Neveu model is as follows

$\left[\bar{\psi }\right]=\left[\psi \right]=\left[\tilde{\mu }\right]=\left[\tilde{m}\right]=\left[\sigma \right]={\text{length}}^{-1}\mathrm{,}\left[\tilde{g}\right]={\text{length}}^{1/2}$ (5)

The partition function for this model is

$\begin{array}{c}Z={\displaystyle \int \text{d}\bar{\psi }\text{d}\psi \text{d}\sigma {\text{e}}^{-{\displaystyle \int \mathcal{L}}}}\\ ={\displaystyle \int \text{d}\sigma {\text{e}}^{-{\displaystyle \int }\frac{N_f}{2{\tilde{g}}^2}{\sigma }^2}{\left[\det \left(\cancel{\partial }+\tilde{\mu }{\gamma }_0+\tilde{m}+\sigma \right)\right]}^{N_f}}\\ ={\displaystyle \int \text{d}\sigma }\exp \left(-{\displaystyle \int }\frac{N_f}{2{\tilde{g}}^2}{\sigma }^2+N_f\ln \left[\det \left(\cancel{\partial }+\tilde{\mu }{\gamma }_0+\tilde{m}+\sigma \right)\right]\right)\end{array}$ (6)

where ${\displaystyle \int }\equiv {\displaystyle {\int }_0^{\beta }}\text{d}{x}_0{\displaystyle {\int }_0^L}\text{d}{x}_1\text{d}{x}_2$ with the inverse temperature $\beta =1/T$ and the space size L. $\bar{\psi }$ and $\psi$ 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

$\frac{1}{N_{f}V}\frac{\partial \ln Z}{\partial \tilde{m}}=\langle -\frac{1}{V}{\displaystyle \int {\bar{\psi }}_i{\psi }_i}\rangle =\frac{1}{{\tilde{g}}^2}\langle \frac{1}{V}{\displaystyle \int \sigma }\rangle \equiv \frac{1}{{\tilde{g}}^2}\Sigma$ (7)

where $V=\beta L^2$ is the volume of 2 + 1d system. In the second equality we used

$0={\displaystyle \int \text{d}\bar{\psi }\text{d}\psi \text{d}\sigma \frac{\delta }{\delta \sigma \left(x\right)}{\text{e}}^{-{\displaystyle \int \mathcal{L}}}}={\displaystyle \int \text{d}\bar{\psi }\text{d}\psi \text{d}\sigma {\text{e}}^{-{\displaystyle \int \mathcal{L}}}\left(-1\right)\left(\bar{\psi }\psi +\frac{N}{{\tilde{g}}^2}\sigma \right)(\; x\; )}$

Since the Lagrangian density is translation invariant, $\langle \bar{\psi }\left(x\right)\psi \left(x\right)\rangle$ and $\langle \sigma \left(x\right)\rangle$ does not depend on x. This model in the large $N_f$ limit can be solved exactly [18] in the chiral limit $\tilde{m}=0$, which is based on the saddle approximation (gap equation) in (6)

$\begin{array}{c}0=-\frac{V}{{\tilde{g}}^2}\Sigma +\frac{\text{d}}{\text{d}\Sigma }\ln \left[\det \left(\cancel{\partial }+\tilde{\mu }{\gamma }_0+\tilde{m}+\Sigma \right)\right]\\ =-\frac{V}{{\tilde{g}}^2}\Sigma +\text{Tr}{\left(\cancel{\partial }+\tilde{\mu }{\gamma }_0+\tilde{m}+\Sigma \right)}^{-1}\\ =-\frac{V}{{\tilde{g}}^2}\Sigma +{\displaystyle \underset{k}{\sum }}\text{tr}{\left(i\cancel{k}+\tilde{\mu }{\gamma }_0+\tilde{m}+\Sigma \right)}^{-1}\\ =-\frac{V}{{\tilde{g}}^2}\Sigma +4\left(\tilde{m}+\Sigma \right){\displaystyle \underset{k}{\sum }}{\left({\left(k_0-i\tilde{\mu }\right)}^2+{\displaystyle \underset{\nu =\mathrm{1,2}}{\sum }}{k}_{\nu }^2+{\left(\tilde{m}+\Sigma \right)}^2\right)}^{-1}\end{array}$ (8)

where in the third equality we write the trace of operator in momentum space and the summation over $k=\left(k_0\mathrm{,}{k}_1\mathrm{,}{k}_2\right)$

$k_0=\left(2n-1\right)\pi T\mathrm{,}{k}_{\nu }=2n_{\nu }\pi /L\mathrm{,}n\mathrm{,}{n}_{\nu }\in Z\mathrm{,}\nu =\mathrm{1,2}$

3. The Staggered Fermion

The staggered fermion discretization of the action ${\displaystyle \int \mathcal{L}}$ is

$\begin{array}{c}S=a^2{a}_t{\displaystyle \underset{x,y}{\sum }}\bar{\psi }\left(x\right)\left({\displaystyle \underset{\alpha =1,2}{\sum }}\frac{{\eta }_{x,\alpha }}{2a}\left({\delta }_{x+\hat{\alpha },y}-{\delta }_{x,y+\hat{\alpha }}\right)\right)\psi \left(y\right)\\ +a^2{a}_t{\displaystyle \underset{x,y}{\sum }}\bar{\psi }\left(x\right)\left(\frac{{\eta }_{x,0}}{2a_t}\left({\text{e}}^{a_t\tilde{\mu }}{s}_x^1{\delta }_{x+\hat{0},y}-{\text{e}}^{-a_t\tilde{\mu }}{s}_x^2{\delta }_{x,y+\hat{0}}\right)\right)\psi \left(y\right)\\ +a^2{a}_t{\displaystyle \underset{x}{\sum }}\left(\tilde{m}+\varphi \left(x\right)\right)\bar{\psi }\left(x\right)\psi \left(x\right)+a^2{a}_t\frac{N_f}{2{\tilde{g}}^2}{\displaystyle \underset{\tilde{x}}{\sum }}\sigma {\left(\tilde{x}\right)}^2\end{array}$ (9)

with staggered phase factor ${\eta }_{x,0}=1$, ${\eta }_{x\mathrm{,1}}={\left(-1\right)}^{x_0/a}$, ${\eta }_{x\mathrm{,2}}={\left(-1\right)}^{\left(x_0+x_1\right)/a}$. $a{N}_x=L$, $a_t{N}_t=\beta =1/T$. The boundary condition for $\psi$ and $\bar{\psi }$ are accounted for by the sign $s^1$ and $s^2$

$s_x^1=\{\begin{array}{ll}-1\hfill & \text{if}{x}_0=N_t-1\hfill \\ 1\hfill & \text{Otherwise}\hfill \end{array}\mathrm{,}{s}_x^2=\{\begin{array}{ll}-1\hfill & \text{if}{x}_0=0\hfill \\ 1\hfill & \text{Otherwise}\hfill \end{array}$ (10)

Here $\varphi$ is defined on lattice $x$ by $\sigma (\; x\; ˜\; )$

$\varphi \left(x\right)=\frac{1}{8}{\displaystyle \underset{\left[x\mathrm{,}\tilde{x}\right]}{\sum }}\sigma \left(\tilde{x}\right)\iff \sigma \left(\tilde{x}\right)=\frac{1}{8}{\displaystyle \underset{\left[x\mathrm{,}\tilde{x}\right]}{\sum }}\varphi \left(x\right)$ (11)

where $\left[x\mathrm{,}\tilde{x}\right]$ denotes 8 dual lattices $\tilde{x}$ which is neighbour to $x$. The auxiliary field on dual lattice for two dimensional Gross-Neveu model was first studied in Ref. [33].

According to (5), the non-dimensional quantities are introduced by

$a\sigma \to \sigma \mathrm{,}a\varphi \to \varphi \mathrm{,}a\bar{\psi }\to \bar{\psi }\mathrm{,}a\psi \to \psi$ (12)

$a\tilde{\mu }=\mu ,a\tilde{m}=m,a^{-1/2}\tilde{g}=g,x/a\to x,a_1=a_t/a$ (13)

and thus the action in (9) can be rewritten as

$\begin{array}{c}S=a_1{\displaystyle \underset{x,y}{\sum }}\bar{\psi }\left(x\right)\left({\displaystyle \underset{\alpha =1,2}{\sum }}\frac{{\eta }_{x,\alpha }}{2}\left({\delta }_{x+\hat{\alpha },y}-{\delta }_{x,y+\hat{\alpha }}\right)\right)\psi \left(y\right)\\ +{\displaystyle \underset{x,y}{\sum }}\bar{\psi }\left(x\right)\left(\frac{{\eta }_{x,0}}{2}\left({\text{e}}^{a_1\mu }{s}_x^1{\delta }_{x+\hat{0},y}-{\text{e}}^{-a_1\mu }{s}_x^2{\delta }_{x,y+\hat{0}}\right)\right)\psi \left(y\right)\\ +a_1{\displaystyle \underset{x}{\sum }}\left(m+\varphi \left(x\right)\right)\bar{\psi }\left(x\right)\psi \left(x\right)+a_1\frac{N_f}{2g^2}{\displaystyle \underset{\tilde{x}}{\sum }}\sigma {\left(\tilde{x}\right)}^2\end{array}$

The partition function for the Gross-Neveu model with $N_f$ flavors is:

$Z={\displaystyle \int {\displaystyle \underset{i}{\prod }\text{d}{\bar{\psi }}_i\text{d}{\psi }_i\text{d}\sigma {\text{e}}^{-S}}}$ (14)

where ${\psi }_i$ and ${\bar{\psi }}_i$ denote the Grassmann fields of flavors $i=0,\cdots ,N_f-1$ at the sites $x$, $\sigma$ is the real field defined at the dual lattice sites $\tilde{x}$. The action is

$S={\displaystyle \underset{i,x,y}{\sum }}{\bar{\psi }}_i\left(x\right)D_{x,y}{\psi }_i\left(y\right)+{\displaystyle \underset{i,x}{\sum }}{a}_1\varphi \left(x\right){\bar{\psi }}_i\left(x\right){\psi }_i\left(x\right)+\frac{a_1{N}_f}{2g^2}{\displaystyle \underset{\tilde{x}}{\sum }}\sigma {\left(\tilde{x}\right)}^2$ (15)

where

$D_{x,y}=\{\begin{array}{ll}{a}_1\frac{{\eta }_{x,\alpha }}{2}\hfill & \text{if}y=x+\hat{\alpha },\alpha =1,2\hfill \\ -a_1\frac{{\eta }_{x,\alpha }}{2}\hfill & \text{if}y=x-\hat{\alpha },\alpha =1,2\hfill \\ \frac{{\eta }_{x,0}}{2}{\text{e}}^{a_1\mu }{s}_x^1\hfill & \text{if}y=x+\hat{0}\hfill \\ -\frac{{\eta }_{x,0}}{2}{\text{e}}^{-a_1\mu }{s}_x^2\hfill & \text{if}y=x-\hat{0}\hfill \\ a_{1}m\hfill & \text{if}y=x\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$ (16)

The derivative of this matrix D with respect to the chemical potential and bare mass are rather simple

$\frac{\partial D_{x,y}}{\partial \left(a_1\mu \right)}=\frac{{\text{e}}^{a_1\mu }}{2}{s}_x^1{\delta }_{x+\hat{0},y}+\frac{{\text{e}}^{-a_1\mu }}{2}{s}_x^2{\delta }_{x,y+\hat{0}},\frac{\partial D_{x,y}}{\partial \left(a_{1}m\right)}={\delta }_{x,y}$

The real matrix $D\left(\mu \mathrm{,}m\right)$ satisfies the following symmetry

$D{\left(\mu \mathrm{,}m\right)}_{x\mathrm{,}y}=-D{\left(-\mu \mathrm{,}-m\right)}_{y\mathrm{,}x}$

${\epsilon }_{x}D{\left(\mu \mathrm{,}m\right)}_{x\mathrm{,}y}{\epsilon }_y=-D{\left(\mu \mathrm{,}-m\right)}_{x\mathrm{,}y}=D{\left(-\mu \mathrm{,}m\right)}_{y\mathrm{,}x}$

where ${\epsilon }_x={\left(-1\right)}^{x_0+x_1+x_2}$ is the parity of site $x$.

By integrating the Grassmann fields, the partition function in (14) can be rewritten as

$Z={\displaystyle \int {\displaystyle \underset{\tilde{x}}{\prod }\text{d}\sigma \left(\tilde{x}\right){\text{e}}^{-S_{\text{eff}}}}}$ (17)

with the effective action

$S_{\text{eff}}=\frac{a_1{N}_f}{2g^2}{\displaystyle \underset{\tilde{x}}{\sum }}{\sigma }^2\left(\tilde{x}\right)-N_f\ln \det D\left[\varphi \right]$ (18)

and

${\left(D\left[\varphi \right]\right)}_{x\mathrm{,}y}=D_{x\mathrm{,}y}+a_1\varphi \left(x\right){\delta }_{x\mathrm{,}y}$ (19)

The computational results, e.g., non-dimensional chiral condensate and fermion density, depend on the non-dimensional quantities

$\left(N_f\mathrm{,}g\mathrm{,}\mu \mathrm{,}m\mathrm{,}{N}_x\mathrm{,}{N}_t\right)$

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.

4. Staggered Fermion in Momentum Space

The kinetic part in (15) in one flavor is ${\displaystyle {\sum }_{x,y}}\bar{\chi }\left(x\right)D_{x,y}\chi \left(y\right)$ where

$D_{x,y}=\{\begin{array}{ll}\frac{{\eta }_{x,\alpha }}{2}\hfill & \text{if}y=x+\hat{\alpha },\alpha =1,2\hfill \\ -\frac{{\eta }_{x,\alpha }}{2}\hfill & \text{if}y=x-\hat{\alpha },\alpha =1,2\hfill \\ \frac{{\eta }_{x,0}}{2}{\text{e}}^{\mu }{s}_x^1\hfill & \text{if}y=x+\hat{0}\hfill \\ -\frac{{\eta }_{x,0}}{2}{\text{e}}^{-\mu }{s}_x^2\hfill & \text{if}y=x-\hat{0}\hfill \\ m\hfill & \text{if}y=x\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$ (20)

$\bar{\chi }$ and $\chi$ are the Grassmann fields defined on lattices. A Wilson-like fermion can be obtained from the stagger fermion ${\displaystyle {\sum }_{x\mathrm{,}y}}\bar{\chi }\left(x\right)D_{x\mathrm{,}y}\chi \left(y\right)$ [18].

Assume that $N_x$ and $N_t$ are even integers. Let $Y=\left(Y_0\mathrm{,}{Y}_1\mathrm{,}{Y}_2\right)$ denotes a site on a lattice of twice the spacing of the original, and $A=\left(A_0,A_1,A_2\right),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

$\chi \left(x\right)=\chi \left(2Y+A\right)=\chi \left(A,Y\right)$

A shift along $\mu$ direction can be represented by

$\begin{array}{c}\chi \left(x+\hat{\mu }\right)=\chi \left(2Y+A+\hat{\mu }\right)=\chi \left(2\left(Y+\hat{\mu }\right)+A-\hat{\mu }\right)\\ ={\displaystyle \underset{A^{\prime }}{\sum }}\left({\delta }_{A+\hat{\mu }\mathrm{,}{A}^{\prime }}\chi \left(A^{\prime }\mathrm{,}Y\right)+{\delta }_{A-\hat{\mu }\mathrm{,}{A}^{\prime }}\chi \left(A^{\prime }\mathrm{,}Y+\hat{\mu }\right)\right)\end{array}$ (21)

Similarly,

$\chi \left(x-\hat{\mu }\right)={\displaystyle \underset{A^{\prime }}{\sum }}\left({\delta }_{A-\hat{\mu },A^{\prime }}\chi \left(A^{\prime },Y\right)+{\delta }_{A+\hat{\mu },A^{\prime }}\chi \left(A^{\prime },Y-\hat{\mu }\right)\right)$ (22)

$\chi \left(x\right)$ is defined on the fine lattice sites $x$ with lattice size $a=1$

$\left\{x=\left(x_0,x_1,x_2\right),0\le x_0<N_t,0\le x_1,x_2<N_x\right\}$ (23)

while $\chi \left(A\mathrm{,}\cdot \right)$ on the coarse lattice sites Y with lattice size $2a=2$

$\left\{2Y=2\left(Y_0,Y_1,Y_2\right),0\le Y_0<N_t/2,0\le Y_1,Y_2<N_x/2\right\}$ (24)

A unitary transformation of $\chi \left(A\mathrm{,}\cdot \right)$ is defined by [34]

$u^{\alpha a}\left(Y\right)=\frac{1}{4\sqrt{2}}{\displaystyle \underset{A}{\sum }}{\Gamma }_A^{\alpha a}\chi \left(A,Y\right),d^{\alpha a}\left(Y\right)=\frac{1}{4\sqrt{2}}{\displaystyle \underset{A}{\sum }}{B}_A^{\alpha a}\chi \left(A,Y\right)$ (25)

${\bar{u}}^{\alpha a}\left(Y\right)=\frac{1}{4\sqrt{2}}{\displaystyle \underset{A}{\sum }}\bar{\chi }\left(A,Y\right){\Gamma }_A^{\ast \alpha a},{\bar{d}}^{\alpha a}\left(Y\right)=\frac{1}{4\sqrt{2}}{\displaystyle \underset{A}{\sum }}\bar{\chi }\left(A,Y\right)B_A^{\ast \alpha a}$ (26)

where $2\times 2$ matrices ${\Gamma }_A$ and $B_A$ is given by

${\Gamma }_A={\sigma }_1^{A_0}{\sigma }_2^{A_1}{\sigma }_3^{A_2},B_A={\left(-{\sigma }_1\right)}^{A_0}{\left(-{\sigma }_2\right)}^{A_1}{\left(-{\sigma }_3\right)}^{A_2}$ (27)

${\Gamma }_A$ and $B_A$ satisfies the following properties (The indices $\alpha \mathrm{,}{\alpha }^{\prime }\mathrm{,}\beta \mathrm{,}a\mathrm{,}{a}^{\prime }$ and b always run from 1 to 2)

${\Gamma }_{A\pm \hat{\mu }}={\eta }_{\mu }\left(A\right){\sigma }_{\mu +1}{\Gamma }_A\mathrm{,}{B}_{A\pm \hat{\mu }}={\eta }_{\mu }\left(A\right)\left(-{\sigma }_{\mu +1}\right)B_A\mathrm{,}\mu =\mathrm{0,1,2}$ (28)

$\text{Tr}\left({\Gamma }_A^{†}{\Gamma }_{A^{\prime }}+B_A^{†}{B}_{A^{\prime }}\right)=4{\delta }_{A{A}^{\prime }}$ (29)

$\begin{array}{l}{\displaystyle \underset{A}{\sum }}{\Gamma }_A^{\alpha a}{\Gamma }_A^{\ast \beta b}={\displaystyle \underset{A}{\sum }}{B}_A^{\alpha a}{B}_A^{\ast \beta b}=4{\delta }_{\alpha \beta }{\delta }_{ab},\\ {\displaystyle \underset{A}{\sum }}{\Gamma }_A^{\alpha a}{B}_A^{\ast \beta b}={\displaystyle \underset{A}{\sum }}{B}_A^{\alpha a}{\Gamma }_A^{\ast \beta b}=0\end{array}$ (30)

${\displaystyle \underset{A,A_{\mu }=1}{\sum }}{\Gamma }_A^{\alpha a}{\left({\Gamma }_A^{*}\right)}^{{\alpha }^{\prime }{a}^{\prime }}={\displaystyle \underset{A,A_{\mu }=0}{\sum }}{\Gamma }_A^{\alpha a}{\left({\Gamma }_A^{*}\right)}^{{\alpha }^{\prime }{a}^{\prime }},\mu =0,1,2$ (31)

Equation (31) is also valid if $\Gamma$ is replaced by B.

${\displaystyle \underset{A,A_{\mu }=1}{\sum }}{\Gamma }_A^{\alpha a}{\left({\sigma }_{\mu +1}^{\ast }{B}_A^{\ast }\right)}^{{\alpha }^{\prime }{a}^{\prime }}=-2{\sigma }_{\mu +1}^{\ast a{a}^{\prime }}{\delta }_{\alpha {\alpha }^{\prime }}$ (32)

${\displaystyle \underset{A,A_{\mu }=0}{\sum }}{\Gamma }_A^{\alpha a}{\left({\sigma }_{\mu +1}^{\ast }{B}_A^{\ast }\right)}^{{\alpha }^{\prime }{a}^{\prime }}=2{\sigma }_{\mu +1}^{\ast a{a}^{\prime }}{\delta }_{\alpha {\alpha }^{\prime }}$ (33)

See Appendix A for these properties.

Using (29), the inverse transformation of (25) and (26) are

$\chi \left(A,Y\right)=\sqrt{2}{\displaystyle \underset{\alpha ,a}{\sum }}\left[{\Gamma }_A^{\ast \alpha a}{u}^{\alpha a}\left(Y\right)+B_A^{\ast \alpha a}{d}^{\alpha a}\left(Y\right)\right]$ (34)

$\bar{\chi }\left(A,Y\right)=\sqrt{2}{\displaystyle \underset{\alpha ,a}{\sum }}\left[{\bar{u}}^{\alpha a}\left(Y\right){\Gamma }_A^{\alpha a}+{\bar{d}}^{\alpha a}\left(Y\right)B_A^{\alpha a}\right]$ (35)

Let us introduce the two Dirac fields with 4 components ( $a=1,2$ )

$q^a\left(Y\right)=\left(\begin{array}{c}{q}_1^a\left(Y\right)\\ q_2^a\left(Y\right)\end{array}\right)=\left(\begin{array}{c}{u}^{\alpha a}\left(Y\right)\\ d^{\alpha a}\left(Y\right)\end{array}\right),{\bar{q}}^a\left(Y\right)=\left({\bar{q}}_1^a\left(Y\right),{\bar{q}}_2^a\left(Y\right)\right)=\left({\bar{u}}^{\alpha a}\left(Y\right),{\bar{d}}^{\alpha a}(\; Y\; )\right)$

From the properties (30), it is easy to show that

$\begin{array}{l}{\displaystyle \underset{x}{\sum }}\bar{\chi }\left(x\right)\chi \left(x\right)\\ ={\displaystyle \underset{A,Y}{\sum }}\sqrt{2}{\displaystyle \underset{\alpha ,a}{\sum }}\left({\bar{u}}^{\alpha a}\left(Y\right){\Gamma }_A^{\alpha a}+{\bar{d}}^{\alpha a}\left(Y\right)B_A^{\alpha a}\right)\sqrt{2}{\displaystyle \underset{{\alpha }^{\prime },a^{\prime }}{\sum }}\left[{\Gamma }_A^{\ast {\alpha }^{\prime }{a}^{\prime }}{u}^{{\alpha }^{\prime }{a}^{\prime }}\left(Y\right)+B_A^{\ast {\alpha }^{\prime }{a}^{\prime }}{d}^{{\alpha }^{\prime }{a}^{\prime }}\left(Y\right)\right]\\ =8{\displaystyle \underset{Y}{\sum }}{\displaystyle \underset{\alpha ,a}{\sum }}\left({\bar{u}}^{\alpha a}\left(Y\right)u^{\alpha a}\left(Y\right)+{\bar{d}}^{\alpha a}\left(Y\right)d^{\alpha a}\left(Y\right)\right)\\ =8{\displaystyle \underset{Y,a}{\sum }}{\bar{q}}^a\left(Y\right)q^a\left(Y\right)=8{\displaystyle \underset{Y}{\sum }}\bar{q}\left(Y\right)q\left(Y\right)=8{\displaystyle \underset{k}{\sum }}\bar{q}\left(k\right)q(\; k\; )\end{array}$

where in the last equality the inner produce between $\bar{q}$ and q is given in momentum space corresponding to the coarse lattice with lattice size 2

$k=2\pi \left(\frac{m_0+\frac{1}{2}}{N_t}\mathrm{,}\frac{m_1}{N_x}\mathrm{,}\frac{m_2}{N_x}\right),0\le m_0<N_t/2,0\le m_1,m_2<N_x/2$ (36)

For any fixed $\mu =0,1,2$,

$\begin{array}{l}\frac{1}{2}{\displaystyle \underset{x}{\sum }}{\eta }_{\mu }\left(x\right)\bar{\chi }\left(x\right)\left(\chi \left(x+\hat{\mu }\right)-\chi \left(x-\hat{\mu }\right)\right)\\ =\frac{1}{2}{\displaystyle \underset{A,A^{\prime },Y}{\sum }}{\eta }_{\mu }\left(A\right)\bar{\chi }\left(A,Y\right)({\delta }_{A+\hat{\mu },A^{\prime }}\left(\chi \left(A^{\prime },Y\right)-\chi \left(A^{\prime },Y-\hat{\mu }\right)\right)\\ +{\delta }_{A-\hat{\mu },A^{\prime }}\left(\chi \left(A^{\prime },Y+\hat{\mu }\right)-\chi \left(A^{\prime },Y\right)\right))\\ =\frac{1}{2}{\displaystyle \underset{A,A^{\prime },Y}{\sum }}{\eta }_{\mu }\left(A\right)\bar{\chi }\left(A,Y\right)(\frac{{\delta }_{A+\hat{\mu },A^{\prime }}+{\delta }_{A-\hat{\mu },A^{\prime }}}{2}{\partial }_{\mu }\chi \left(A^{\prime },Y\right)\\ +\frac{{\delta }_{A-\hat{\mu },A^{\prime }}-{\delta }_{A+\hat{\mu },A^{\prime }}}{2}{\partial }_{\mu }^2\chi \left(A^{\prime },Y\right))\end{array}$

$\begin{array}{l}=\frac{1}{2}{\displaystyle \underset{A,A^{\prime },Y}{\sum }}{\eta }_{\mu }\left(A\right)\sqrt{2}{\displaystyle \underset{\alpha ,a}{\sum }}\left({\bar{u}}^{\alpha a}\left(Y\right){\Gamma }_A^{\alpha a}+{\bar{d}}^{\alpha a}\left(Y\right)B_A^{\alpha a}\right)\\ \times \{\frac{{\delta }_{A+\hat{\mu },A^{\prime }}+{\delta }_{A-\hat{\mu },A^{\prime }}}{2}\sqrt{2}{\displaystyle \underset{{\alpha }^{\prime },a^{\prime }}{\sum }}\left({\Gamma }_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{\partial }_{\mu }{u}^{{\alpha }^{\prime }{a}^{\prime }}\left(Y\right)+B_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{\partial }_{\mu }{d}^{{\alpha }^{\prime }{a}^{\prime }}\left(Y\right)\right)\\ +\frac{{\delta }_{A-\hat{\mu },A^{\prime }}-{\delta }_{A+\hat{\mu },A^{\prime }}}{2}\sqrt{2}{\displaystyle \underset{{\alpha }^{\prime },a^{\prime }}{\sum }}\left({\Gamma }_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{\partial }_{\mu }^2{u}^{{\alpha }^{\prime }{a}^{\prime }}\left(Y\right)+B_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{\partial }_{\mu }^2{d}^{{\alpha }^{\prime }{a}^{\prime }}\left(Y\right)\right)\}\end{array}$

where in the second equality (21) and (22) are used. According to the properties of ${\Gamma }_A$ and $B_A$ in (30) (31) (32) and (33)

$\begin{array}{l}\frac{1}{2}{\displaystyle \underset{x}{\sum }}{\eta }_{\mu }\left(x\right)\bar{\chi }\left(x\right)\left(\chi \left(x+\hat{\mu }\right)-\chi \left(x-\hat{\mu }\right)\right)\\ =2{\displaystyle \underset{Y}{\sum }}({\bar{u}}^{\alpha a}\left(Y\right){\left({\sigma }_{\mu +1}\right)}^{\alpha {\alpha }^{\prime }}{\delta }_{a{a}^{\prime }}{\partial }_{\mu }{u}^{{\alpha }^{\prime }{a}^{\prime }}\left(Y\right)+{\bar{d}}^{\alpha a}\left(Y\right){\left(-{\sigma }_{\mu +1}\right)}^{\alpha {\alpha }^{\prime }}{\delta }_{a{a}^{\prime }}{\partial }_{\mu }{d}^{{\alpha }^{\prime }{a}^{\prime }}\left(Y\right)\\ +{\bar{u}}^{\alpha a}\left(Y\right){\left({\sigma }_{\mu +1}^{*}\right)}^{a{a}^{\prime }}{\delta }_{\alpha {\alpha }^{\prime }}{\partial }_{\mu }^2{d}^{{\alpha }^{\prime }{a}^{\prime }}\left(Y\right)+{\bar{d}}^{\alpha a}\left(Y\right){\left(-{\sigma }_{\mu +1}^{*}\right)}^{a{a}^{\prime }}{\delta }_{\alpha {\alpha }^{\prime }}{\partial }_{\mu }^2{u}^{{\alpha }^{\prime }{a}^{\prime }}\left(Y\right))\\ =2{\displaystyle \underset{Y}{\sum }}\left[\bar{q}\left(Y\right)\left({\gamma }_{\mu }\otimes {\mathbb{I}}_2\right){\partial }_{\mu }q\left(Y\right)+\bar{q}\left(Y\right)\left(i{\gamma }_3\otimes {\sigma }_{\mu +1}^{*}\right){\partial }_{\mu }^{2}q\left(Y\right)\right]\\ =8{\displaystyle \underset{Y}{\sum }}\left[\bar{q}\left(Y\right)\left({\gamma }_{\mu }\otimes {\mathbb{I}}_2\right)\frac{{\partial }_{\mu }q\left(Y\right)}{4}+\bar{q}\left(Y\right)\left(i{\gamma }_3\otimes {\sigma }_{\mu +1}^{*}\right)\frac{{\partial }_{\mu }^{2}q\left(Y\right)}{4}\right]\\ =8{\displaystyle \underset{k}{\sum }}\left[\bar{q}\left(k\right)\left({\gamma }_{\mu }\otimes {\mathbb{I}}_2\right)\frac{i}{2}\sin \left(2k_{\mu }\right)q\left(k\right)+\bar{q}\left(k\right)\left(i{\gamma }_3\otimes {\sigma }_{\mu +1}^{*}\right)\frac{1}{2}\left[\cos \left(2k_{\mu }\right)-1\right]q\left(k\right)\right]\end{array}$ (37)

where we used the notations

${\partial }_{\mu }q\left(Y\right)=q\left(Y+\hat{\mu }\right)-q\left(Y-\hat{\mu }\right)$

${\partial }_{\mu }^{2}q\left(Y\right)=q\left(Y+\hat{\mu }\right)-2q\left(Y\right)+q\left(Y-\hat{\mu }\right)$

and the summation over $k$ is taken for all modes in (36). Similarly, we have (see Appendix B)

$\begin{array}{l}\frac{1}{2}{\displaystyle \underset{x}{\sum }}\bar{\chi }\left(x\right)\left(\chi \left(x+\hat{0}\right)+\chi \left(x-\hat{0}\right)\right)\\ =8{\displaystyle \underset{k}{\sum }}[\bar{q}\left(k\right)\left(i{\gamma }_3\otimes {\sigma }_1^{*}\right)i{2}^{-1}\sin \left(2k_0\right)q\left(k\right)\\ +\bar{q}\left(k\right)\left({\gamma }_0\otimes {\mathbb{I}}_2\right)2^{-1}\left[\cos \left(2k_0\right)+1\right]q\left(k\right)]\\ \equiv 8{\displaystyle \underset{k}{\sum }}\bar{q}\left(k\right)A_{+}\left(k\right)q\left(k\right)\end{array}$ (38)

Using

$\begin{array}{l}\frac{1}{2}{\displaystyle \underset{x}{\sum }}\bar{\chi }\left(x\right)\left({\text{e}}^{\mu }\chi \left(x+\hat{0}\right)-{\text{e}}^{-\mu }\chi \left(x-\hat{0}\right)\right)\\ =\cosh \mu \left[\frac{1}{2}{\displaystyle \underset{x}{\sum }}\bar{\chi }\left(x\right)\left(\chi \left(x+\hat{0}\right)-\chi \left(x-\hat{0}\right)\right)\right]\\ +\sinh \mu \left[\frac{1}{2}{\displaystyle \underset{x}{\sum }}\bar{\chi }\left(x\right)\left(\chi \left(x+\hat{0}\right)+\chi \left(x-\hat{0}\right)\right)\right]\end{array}$

and (37) (38), the kinetic part ${\displaystyle {\sum }_{x\mathrm{,}y}}\bar{\chi }\left(x\right)D_{x\mathrm{,}y}\chi \left(y\right)$ can be rewritten as in the momentum space

${\displaystyle \underset{x,y}{\sum }}\bar{\chi }\left(x\right)D_{x,y}\chi \left(y\right)=8{\displaystyle \underset{k}{\sum }}\bar{q}\left(k\right)D\left(k\right)q\left(k\right)$ (39)

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

$\begin{array}{c}D\left(k\right)=m+{\displaystyle \underset{\mu =\mathrm{1,2}}{\sum }}\frac{i}{2}\left\{\left({\gamma }_{\mu }\otimes {\mathbb{I}}_2\right)\sin \left(2k_{\mu }\right)+\left({\gamma }_3\otimes {\sigma }_{\mu +1}^{\mathrm{<mo>\; *\; </mo>}}\right)\left[\cos \left(2k_{\mu }\right)-1\right]\right\}\\ +\frac{1}{2}\{\left({\gamma }_0\otimes {\mathbb{I}}_2\right)\left(i\cosh \mu \sin \left(2k_0\right)+\sinh \mu \left[\cos \left(2k_0\right)+1\right]\right)\\ +\left({\gamma }_3\otimes {\sigma }_1^{\mathrm{<mo>\; *\; </mo>}}\right)\left(i\cosh \mu \left[\cos \left(2k_0\right)-1\right]-\sinh \mu \sin \left(2k_0\right)\right)\}\\ \equiv m+{\displaystyle \underset{\mu =\mathrm{0,1,2}}{\sum }}\left({\gamma }_{\mu }\otimes {\mathbb{I}}_2\right)a_{\mu }+{\displaystyle \underset{c=\mathrm{1,2,3}}{\sum }}\left({\gamma }_3\otimes {\sigma }_c^{\mathrm{<mo>\; *\; </mo>}}\right)b_c\end{array}$ (40)

where $a_{\mu }$ and $b_c$ depends on k. The inverse matrix of $D\left(k\right)$ is

$D{\left(k\right)}^{-1}=\frac{1}{N\left(k\right)}\left[m-{\displaystyle \underset{\mu =0,1,2}{\sum }}\left({\gamma }_{\mu }\otimes {\mathbb{I}}_2\right)a_{\mu }-{\displaystyle \underset{c=1,2,3}{\sum }}\left({\gamma }_3\otimes {\sigma }_c^{*}\right)b_c\right]$ (41)

where

$\begin{array}{c}N\left(k\right)=m^2+\frac{1}{4}{\displaystyle \underset{\mu =0,1,2}{\sum }}{\left(\sin 2k_{\mu }\right)}^2+\frac{1}{4}{\displaystyle \underset{\mu =0,1,2}{\sum }}{\left(1-\cos 2k_{\mu }\right)}^2\\ -{\sinh }^2\mu \cos 2k_0-i\cosh \mu \sinh \mu \sin 2k_0\end{array}$ (42)

We can calculate the trace of inverse matrix D in (20) from (39)

$\begin{array}{c}{\displaystyle \underset{x}{\sum }}{D}_{x\mathrm{,}x}^{-1}=-\frac{{\displaystyle \int {\text{e}}^{-\bar{\chi }D\chi }{\displaystyle {\sum }_x\bar{\chi }\left(x\right)\chi \left(x\right)}}}{{\displaystyle \int {\text{e}}^{-\bar{\chi }D\chi }}}\\ =-\frac{{\displaystyle \int {\text{e}}^{-{\displaystyle {\sum }_k\bar{q}\left(k\right)8D\left(k\right)q\left(k\right)}}8{\displaystyle {\sum }_k\bar{q}\left(k\right)q\left(k\right)}}}{{\displaystyle \int {\text{e}}^{-{\displaystyle {\sum }_k\bar{q}\left(k\right)8D\left(k\right)q\left(k\right)}}}}\\ =-8{\displaystyle \underset{k}{\sum }}\frac{{\displaystyle \int {\text{e}}^{-\bar{q}\left(k\right)8D\left(k\right)q\left(k\right)}\bar{q}\left(k\right)q\left(k\right)}}{{\displaystyle \int {\text{e}}^{-\bar{q}\left(k\right)8D\left(k\right)q\left(k\right)}}}\\ =8{\displaystyle \underset{k}{\sum }}\text{tr}\left[{\left(8D\left(k\right)\right)}^{-1}\right]\\ ={\displaystyle \underset{k}{\sum }}\text{tr}\left[D{\left(k\right)}^{-1}\right]={\displaystyle \underset{k}{\sum }}\frac{8m}{N\left(k\right)}\end{array}$ (43)

where the summation over $k$ is given by (36). Note that the right hand side of (43) is real since ${\displaystyle {\sum }_{k_0}}\sin 2k_0/{\left|N\left(k\right)\right|}^2=0$ for any $k_1$ and $k_2$ modes in (36). Similarly,

${\displaystyle \underset{x}{\sum }}\left(D_{x+\hat{0},x}^{-1}{s}_x^1+D_{x-\hat{0},x}^{-1}{s}_x^2\right)=8{\displaystyle \underset{k}{\sum }}\frac{b_1\sin 2k_0-a_0\left(\cos 2k_0+1\right)}{N\left(k\right)}$ (44)

and

${\displaystyle \underset{x}{\sum }}\left(D_{x+\hat{0},x}^{-1}{s}_x^1-D_{x-\hat{0},x}^{-1}{s}_x^2\right)=\left(-8i\right){\displaystyle \underset{k}{\sum }}\frac{a_0\sin 2k_0+b_1\left(\cos 2k_0-1\right)}{N\left(k\right)}$ (45)

The inverse matrix of D in (20) is

$\begin{array}{l}{D}_{x^{\prime },x}^{-1}=\frac{1}{4}{\displaystyle \underset{\alpha ,a,{\alpha }^{\prime },a^{\prime }}{\sum }}[{\Gamma }_A^{\alpha a}{\Gamma }_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{D}_{\left(Y^{\prime }{a}^{\prime }{\alpha }^{\prime }1;Ya\alpha 1\right)}^{-1}+{\Gamma }_A^{\alpha a}{B}_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{D}_{\left(Y^{\prime }{a}^{\prime }{\alpha }^{\prime }2;Ya\alpha 1\right)}^{-1}\\ +B_A^{\alpha a}{\Gamma }_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{D}_{\left(Y^{\prime }{a}^{\prime }{\alpha }^{\prime }1;Ya\alpha 2\right)}^{-1}+B_A^{\alpha a}{B}_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{D}_{\left(Y^{\prime }{a}^{\prime }{\alpha }^{\prime }2;Ya\alpha 2\right)}^{-1}]\end{array}$ (46)

See Appendix C for the derivation of (44)-(46).

Since D is diagonal in momentum space, the inverse matrix in the $\bar{q}q$ basis is

$\begin{array}{l}{D}_{Y^{\prime };Y}^{-1}=\frac{1}{N_t/2{\left(N_x/2\right)}^2}{\displaystyle \underset{k}{\sum }}{\text{e}}^{ik\cdot 2\left(Y^{\prime }-Y\right)}{D}^{-1}\left(k\right)\\ =\frac{1}{N_t/2{\left(N_x/2\right)}^2}{\displaystyle \underset{k}{\sum }}{\text{e}}^{ik\cdot 2\left(Y^{\prime }-Y\right)}\frac{1}{N\left(k\right)}\left[m-{\displaystyle \underset{\mu =0,1,2}{\sum }}\left({\gamma }_{\mu }\otimes {\mathbb{I}}_2\right)a_{\mu }-{\displaystyle \underset{c=1,2,3}{\sum }}\left({\gamma }_3\otimes {\sigma }_c^{*}\right)b_c\right]\\ \equiv m\left({\mathbb{I}}_4\otimes {\mathbb{I}}_2\right)\tilde{1}\left(Y^{\prime }-Y\right)-{\displaystyle \underset{\mu =0,1,2}{\sum }}\left({\gamma }_{\mu }\otimes {\mathbb{I}}_2\right){\tilde{a}}_{\mu }\left(Y^{\prime }-Y\right)-{\displaystyle \underset{c=1,2,3}{\sum }}\left({\gamma }_3\otimes {\sigma }_c^{*}\right){\tilde{b}}_c\left(Y^{\prime }-Y\right)\end{array}$

where the notation with tilde denotes the inverse Fourier transformation, e.g.,

$\begin{array}{c}{\tilde{a}}_{\mu }\left(Y\right)=\frac{1}{\frac{N_t}{2}{\left(\frac{N_x}{2}\right)}^2}{\displaystyle \underset{k}{\sum }}{\text{e}}^{ik\cdot 2Y}\frac{a_{\mu }\left(k\right)}{N\left(k\right)}\\ ={\text{e}}^{i\frac{2\pi Y_0}{N_t}}\frac{1}{\frac{N_t}{2}{\left(\frac{N_x}{2}\right)}^2}{\displaystyle \underset{m_0=0}{\overset{\frac{N_t}{2}-1}{\sum }}}{\displaystyle \underset{m_1=0}{\overset{\frac{N_x}{2}-1}{\sum }}}{\displaystyle \underset{m_2=0}{\overset{\frac{N_x}{2}-1}{\sum }}}{\text{e}}^{i2\pi \left(\frac{m_0{Y}_0}{N_t/2}+\frac{m_1{Y}_1}{N_x/2}+\frac{m_2{Y}_2}{N_x/2}\right)}\frac{a_{\mu }\left(m_0,m_1,m_2\right)}{N\left(m_0,m_1,m_2\right)}\end{array}$

for $\left|Y_0\right|\le \frac{N_t}{2}-1$, $\left|Y_1\right|\mathrm{,}\left|Y_2\right|\le \frac{N_x}{2}-1$. We first use the fast Fourier transformation to calculate ${\tilde{a}}_{\mu }\left(Y\right)\exp \left(-i\frac{2\pi Y_0}{N_t}\right)$ and thus ${\tilde{a}}_{\mu }\left(Y\right)$ for $0\le Y_0\le \frac{N_t}{2}-1$, $0\le Y_1\mathrm{,}{Y}_2\le \frac{N_x}{2}-1$. Then ${\tilde{a}}_{\mu }\left(Y\right)$ for $\left|Y_0\right|\le \frac{N_t}{2}-1$, $\left|Y_1\right|\mathrm{,}\left|Y_2\right|\le \frac{N_x}{2}-1$ can be obtained since it is anti-periodic in $Y_0$ direction and periodic in $Y_1$ and $Y_2$ direction.

Each term in $D_{Y^{\prime };Y}^{-1}$ has a tensor product $A\otimes B$ between $4\times 4$ matrix $A={\left(A_{ij}\right)}_{i,j=1,2}$ with $2\times 2$ matrix $A_{ij}$ and $2\times 2$ matrix B. The indices of $D_{\left(Y^{\prime }{a}^{\prime }{\alpha }^{\prime }i\mathrm{<mo>\; ;\; </mo>}Ya\alpha j\right)}^{-1}$ of the inverse matrix $D_{Y^{\prime };Y}^{-1}$ in (46) is related to ${\left(A_{ij}\right)}_{{\alpha }^{\prime }\alpha }{B}_{a^{\prime }a}$. The analytic formula for the inverse matrix of the staggered fermion is the main contribution of this paper. Compared to the computational complexity $O\left({\left(N_t{N}_x^2\right)}^3\right)$ of the usual inverse matrix, the computational cost is $O\left(16{\left(N_t{N}_x^2\right)}^2\right)$ since each element of the inverse matrix needs the summation over $\alpha ,a,{\alpha }^{\prime },a^{\prime }=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)

${\displaystyle \underset{x}{\sum }}{D}_{x,x}^{-1}={\displaystyle \underset{\alpha ,a}{\sum }}\left[D_{\left(Ya\alpha 1;Ya\alpha 1\right)}^{-1}+D_{\left(Ya\alpha 2;Ya\alpha 2\right)}^{-1}\right]={\displaystyle \underset{k}{\sum }}\frac{8m}{N(\; k\; )}$

5. The 1 + 1d and 3 + 1d Staggered Fermion

The staggered fermion matrix in (20) can be generalized to the 1 + 1d and 3 + 1d case, where $\alpha$ is 1 for the 1 + 1d case and $\alpha$ run from 1 to 3 for the 3 + 1d case.

For the 1 + 1d case, the $2\times 2$ matrices ${\gamma }_{\mu }$ are defined to be

${\gamma }_{\mu }={\sigma }_{\mu },\mu =1,2,{\gamma }_5=i{\gamma }_1{\gamma }_2,{\gamma }_{\mu }{\gamma }_{\nu }+{\gamma }_{\nu }{\gamma }_{\mu }={\delta }_{\mu \nu }2{\mathbb{I}}_2,\mu ,\nu =1,2,5$

The unitary transformation in (25) and (26) are modified to be

${\psi }^{\alpha a}\left(Y\right)=\frac{1}{2}{\displaystyle \underset{A}{\sum }}{\Gamma }_A^{\alpha a}\chi \left(A,Y\right),{\bar{\psi }}^{\alpha a}\left(Y\right)=\frac{1}{2}{\displaystyle \underset{A}{\sum }}\bar{\chi }\left(A,Y\right){\Gamma }_A^{\ast \alpha a}$

The kinetic part ${\displaystyle {\sum }_{x\mathrm{,}y}}\bar{\chi }\left(x\right)D_{x\mathrm{,}y}\chi \left(y\right)$ can be written as

${\displaystyle \underset{x,y}{\sum }}\bar{\chi }\left(x\right)D_{x,y}\chi \left(y\right)={\displaystyle \underset{k}{\sum }}\bar{\psi }\left(k\right)D\left(k\right)\psi \left(k\right)$ (47)

where the summation is taken over all modes

$k=2\pi \left(\frac{m_0+\frac{1}{2}}{N_t},\frac{m_1}{N_x}\right),0\le m_0<N_t/2,0\le m_1<N_x/2$ (48)

The fermion matrix in momentum space is diagonal

$\begin{array}{c}D\left(k\right)=2m+{\displaystyle \underset{\mu =1}{\sum }}\left\{\left({\gamma }_{\mu +1}\otimes {\mathbb{I}}_4\right)i\sin \left(2k_{\mu }\right)+\left({\gamma }_5\otimes {\gamma }_{\mu +1}^{*}{\gamma }_5^{*}\right)\left[\cos \left(2k_{\mu }\right)-1\right]\right\}\\ +\{\left({\gamma }_1\otimes {\mathbb{I}}_4\right)\left(i\cosh \mu \sin \left(2k_0\right)+\sinh \mu \left[\cos \left(2k_0\right)+1\right]\right)\\ +\left({\gamma }_5\otimes {\gamma }_1^{*}{\gamma }_5^{*}\right)\left(\cosh \mu \left[\cos \left(2k_0\right)-1\right]+i\sinh \mu \sin \left(2k_0\right)\right)\}\\ \equiv 2m+{\displaystyle \underset{\mu =0,1}{\sum }}\left({\gamma }_{\mu +1}\otimes {\mathbb{I}}_4\right)a_{\mu }+{\displaystyle \underset{\mu =0,1}{\sum }}\left({\gamma }_5\otimes {\gamma }_{\mu +1}^{*}{\gamma }_5^{*}\right)b_{\mu }\end{array}$ (49)

with its inverse

$D{\left(k\right)}^{-1}=\frac{1}{N\left(k\right)}\left[2m-{\displaystyle \underset{\mu =0,1}{\sum }}\left({\gamma }_{\mu +1}\otimes {\mathbb{I}}_4\right)a_{\mu }-{\displaystyle \underset{\mu =0,1}{\sum }}\left({\gamma }_5\otimes {\gamma }_{\mu +1}^{*}{\gamma }_5^{*}\right)b_{\mu }\right]$ (50)

where

$\begin{array}{c}N\left(k\right)=4m^2+{\displaystyle \underset{\mu =1}{\sum }}{\left(\sin 2k_{\mu }\right)}^2-{\left(i\cosh \mu \sin 2k_0+\sinh \mu \left(\cos 2k_0+1\right)\right)}^2\\ +{\displaystyle \underset{\mu =1}{\sum }}{\left(1-\cos 2k_{\mu }\right)}^2+{\left(\cosh \mu \left(\cos 2k_0-1\right)+i\sinh \mu \sin 2k_0\right)}^2\end{array}$ (51)

The trace of the inverse matrix is

${\displaystyle \underset{x}{\sum }}{D}_{x\mathrm{,}x}^{-1}={\displaystyle \underset{k}{\sum }}\frac{16m}{N\left(k\right)}$ (52)

The inverse matrix of D can be calculated

$D_{x^{\prime }\mathrm{,}x}^{-1}={\displaystyle \underset{\alpha \mathrm{,}a\mathrm{,}{\alpha }^{\prime },a^{\prime }}{\sum }}{\Gamma }_A^{\alpha a}{\Gamma }_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{D}_{\left(Y^{\prime }{a}^{\prime }{\alpha }^{\prime }\mathrm{<mo>\; ;\; </mo>}Ya\alpha \right)}^{-1}$ (53)

where

$D_{Y^{\prime };Y}^{-1}=\frac{1}{N_t/2\left(N_x/2\right)}{\displaystyle \underset{k}{\sum }}{\text{e}}^{ik\cdot 2\left(Y^{\prime }-Y\right)}{D}^{-1}\left(k\right)$ (54)

For the 3 + 1d case, the $4\times 4$ matrices ${\gamma }_{\mu }$ are defined to be

${\Gamma }_A={\gamma }_1^{A_0}{\gamma }_2^{A_1}{\gamma }_3^{A_2}{\gamma }_4^{A_3},\mu =1,2,3,4,{\gamma }_5={\gamma }_1{\gamma }_2{\gamma }_3{\gamma }_4$

${\gamma }_{\mu }{\gamma }_{\nu }+{\gamma }_{\nu }{\gamma }_{\mu }={\delta }_{\mu \nu }2{\mathbb{I}}_2,\mu ,\nu =1,2,3,4,5$

The unitary transformation in (25) and (26) are modified to be

${\psi }^{\alpha a}\left(Y\right)=\frac{1}{2\sqrt{2}}{\displaystyle \underset{A}{\sum }}{\Gamma }_A^{\alpha a}\chi \left(A\mathrm{,}Y\right)\mathrm{,}{\bar{\psi }}^{\alpha a}\left(Y\right)=\frac{1}{2\sqrt{2}}{\displaystyle \underset{A}{\sum }}\bar{\chi }\left(A\mathrm{,}Y\right){\Gamma }_A^{\ast \alpha a}$

The kinetic part can also be written as (47) where the summation is taken for all modes

$k=2\pi \left(\frac{m_0+\frac{1}{2}}{N_t},\frac{m_1}{N_x},\frac{m_2}{N_x},\frac{m_3}{N_x}\right),0\le m_0<N_t/2,0\le m_1,m_2,m_3<N_x/2$

Equations (49) - (51) are still valid except that $\mu$ runs from 1 to 3. Equations (52) - (54) are modified to be

${\displaystyle \underset{x}{\sum }}{D}_{x,x}^{-1}={\displaystyle \underset{k}{\sum }}\frac{64m}{N\left(k\right)}$ (55)

$D_{x^{\prime },x}^{-1}=\frac{1}{2}{\displaystyle \underset{\alpha ,a,{\alpha }^{\prime },a^{\prime }}{\sum }}{\Gamma }_A^{\alpha a}{\Gamma }_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{D}_{\left(Y^{\prime }{a}^{\prime }{\alpha }^{\prime }\mathrm{<mo>\; ;\; </mo>}Ya\alpha \right)}^{-1}$ (56)

$D_{Y^{\prime }\mathrm{<mo>\; ;\; </mo>}Y}^{-1}=\frac{1}{N_t/2{\left(N_x/2\right)}^3}{\displaystyle \underset{k}{\sum }}{\text{e}}^{ik\cdot 2\left(Y^{\prime }-Y\right)}{D}^{-1}\left(k\right)$ (57)

respectively. We have checked the formula (46), (53), (56) for the inverse matrices by Matlab.

6. 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\to \infty$,

$\frac{\Sigma }{g^2}=\frac{1}{N_t{N}_x^2}{\displaystyle \underset{x}{\sum }}{D}_{x,x}^{-1}$ (58)

Here D is defined in (20) where $m$ is replaced by $m+\Sigma$. The right hand side of (58) can be calculated from (42), (43) where $m$ is replaced by $m+\Sigma$. The first derivative of ${\Sigma }^2$ with respect to $\mu$ can be computed from the gap equation (For simplicity, we assume that $m=0$ )

$\frac{\partial {\Sigma }^2}{\partial \mu }=\frac{{\displaystyle {\sum }_k\left(\sinh 2\mu \cos 2k_0+i\cosh 2\mu \sin 2k_0\right)N{\left(k\right)}^{-2}}}{{\displaystyle {\sum }_{k}N{\left(k\right)}^{-2}}}$ (59)

If the average $\Sigma$ of $\sigma$ has been calculated from the gap equation, the free energy density in the large $N_f$ limit is

$\ln Z=-N_t{N}_x^2\frac{{\Sigma }^2}{2g^2}+\ln \det D$

where $\ln \det D={\displaystyle {\prod }_k}\det D\left(k\right)$ up to a constant. The other thermodynamic quantities can be calculated. For example, the fermion density can be analytically calculated

$\begin{array}{c}\frac{1}{N_t{N}_x^2}\frac{\partial \ln Z}{\partial \mu }=-\frac{1}{2g^2}\frac{\partial {\Sigma }^2}{\partial \mu }+\frac{1}{N_t{N}_x^2}\left({\text{e}}^{\mu }{\displaystyle \underset{x}{\sum }}{D}_{x+\hat{0},x}^{-1}{s}_x^1+{\text{e}}^{-\mu }{\displaystyle \underset{x}{\sum }}{D}_{x-\hat{0},x}^{-1}{s}_x^2\right)\\ =-\frac{1}{2g^2}\frac{\partial {\Sigma }^2}{\partial \mu }+\frac{1}{N_t{N}_x^2}(\cosh \mu {\displaystyle \underset{x}{\sum }}\left(D_{x+\hat{0},x}^{-1}{s}_x^1+D_{x-\hat{0},x}^{-1}{s}_x^2\right)\\ +\sinh \mu {\displaystyle \underset{x}{\sum }}\left(D_{x+\hat{0},x}^{-1}{s}_x^1-D_{x-\hat{0},x}^{-1}{s}_x^2\right))\end{array}$ (60)

where $\frac{\partial {\Sigma }^2}{\partial \mu }$, and two sums over $x$ in (60) are given in (59), (44) and (45), respectively. The $N\left(k\right)$ for each mode $k$ in (44), (45), (59) is given by (42) with the replacement of $m$ by $m+\Sigma$ (Here for simplicity we assume that $m=0$ ) and $\Sigma$ is solved from the gap equation (58).

7. Simulation Results

7.1. Large Volume Limit

Let us consider the large volume limit for the non-interacting 2 + 1d Gross-Neveu model. The partition function $Z={\displaystyle \int \text{d}\bar{\chi }\text{d}\chi {\text{e}}^{-\bar{\chi }D\chi }}=\det D$, where the stagger fermion matrix D is given by (20). The ratio of the non-dimensional chiral condensate $a^2\langle \bar{\psi }\psi \rangle$ and non-dimensional mass $m=a\tilde{m}$ is

$\frac{a^2\langle \bar{\psi }\psi \rangle }{a\tilde{m}}=\frac{\langle \bar{\chi }\chi \rangle }{a\tilde{m}}=\frac{\langle {\displaystyle {\sum }_x\bar{\chi }\left(x\right)\chi \left(x\right)}\rangle }{a\tilde{m}\left(N_t{N}_x^2\right)}=\frac{{\displaystyle {\sum }_x{D}_{x,x}^{-1}}}{a\tilde{m}\left(N_t{N}_x^2\right)}=\frac{8}{N_t{N}_x^2}{\displaystyle \underset{k}{\sum }}\frac{1}{N\left(k\right)}$ (61)

where in the last equality we used Equation (43) where $N\left(k\right)$, depending on $m$ and $\mu$, is given by (61). Note that there are $N_t{N}_x^2/8$ modes $k$ in (61). The ratio of the non-dimensional fermion density $a^3\rho$ and ${\left(a\tilde{\mu }\right)}^3$

$\begin{array}{c}\frac{a^3\rho }{{\left(a\tilde{\mu }\right)}^3}=\frac{1}{{\tilde{\mu }}^3}\left(\frac{1/\beta }{\beta L^2}\right)\frac{\partial \ln Z}{\partial \tilde{\mu }}=\frac{1}{N_t\beta L^2{\tilde{\mu }}^3}\frac{\partial \ln Z}{\partial \mu }\\ =\frac{1}{N_t\beta L^2{\tilde{\mu }}^3}[\frac{\cosh \mu }{2}8{\displaystyle \underset{k}{\sum }}\frac{b_1\sin 2k_0-a_0\left(\cos 2k_0+1\right)}{N\left(k\right)}\\ +\frac{\sinh \mu }{2}\left(-8i\right){\displaystyle \underset{k}{\sum }}\frac{a_0\sin 2k_0+b_1\left(\cos 2k_0-1\right)}{N\left(k\right)}]\end{array}$ (62)

where in the last equality we used (44) and (45).

We consider the case $L=\beta$, $a=a_t$ and thus $N_x=N_t\equiv N$. We fix $\tilde{\mu }L$ and $\tilde{m}L$ and then calculate $\frac{a^2\langle \bar{\psi }\psi \rangle }{a\tilde{m}}$ and $\frac{a^3\rho }{{\left(a\tilde{\mu }\right)}^3}$ in the large $N$ limit for fixed lattice size $a$. In fact $\frac{a^2\langle \bar{\psi }\psi \rangle }{a\tilde{m}}$ and $\frac{a^3\rho }{{\left(a\tilde{\mu }\right)}^3}$ does not depend on the lattice size $a$ since the non-dimensional mass $m=a\tilde{m}=\frac{\tilde{m}L}{N}$ and non-dimensional chemical potential $\mu =a\tilde{\mu }=\frac{\tilde{\mu }L}{N}$ does not depends on lattice size $a$. Figure 1 shows the dependence of $\frac{a^2\langle \bar{\psi }\psi \rangle }{a\tilde{m}}$ on $N$ with fixed $\tilde{\mu }L,\tilde{m}L=0,1$. The linear fitting with respect to 1/N shows that the large $N$ limit of $\frac{a^2\langle \bar{\psi }\psi \rangle }{a\tilde{m}}$ is close to 1.008 for all four cases, this is because $m=1/N$ and $\mu =1/N$ both vanish for large $N$ limit. Figure 2 shows the dependence of $\frac{a^3\rho }{{\left(a\tilde{\mu }\right)}^3}$ on $N$, where $\tilde{\mu }L=1$ and $\tilde{m}L=0,1$. The large $N$ limit is close to 1.9271 for $m=0$ and 1.9234 for $m=0.1/N$, respectively.

7.2. Phase Diagram

The phase diagram of the 2 + 1d Gross-Neveu model in the large $N_f$ limit is well known [16] [17] [18]. In this limit the phase diagram of $\left(g^{-2}\mathrm{,}\mu \mathrm{,}T\right)$ is based on the calculation of $\Sigma$. Basically for $T=0$ and $\mu =0$, there is a critical coupling $g_c^{-2}$ such that the chiral symmetry is broken $\Sigma >0$ if the coupling is

Figure 1. The dependence of $\frac{a^2\langle \bar{\psi }\psi \rangle }{a\tilde{m}}$ on $N$, $N=4,8,16,32,64,128,256,512$. (1) $m=1/N,\mu =1/N$ with fitting $-0.9563/N+1.009$, (2) $m=1/N,\mu =0$ with fitting $-0.6051/N+1.008$, (3) $m=0,\mu =1/N$ with fitting $-0.7904/N+1.008$, (4) $m=0,\mu =0$ with fitting $-0.3224/N+1.007$.

Figure 2. The dependence of $\frac{a^3\rho }{{\left(a\tilde{\mu }\right)}^4}$ on $N$, $N=8,16,32,64,128,256$. (1) $m=0,\mu =1/N$ with fitting $14.4370/N^2-0.4345/N+1.9271$, (2) $m=0.1/N,\mu =0$ with fitting $14.4288/N^2-0.4343/N+1.9234$.

strong enough $g^{-2}<g_c^{-2}$. This critical coupling depends on the regularization of the continuum model. For the lattice regularization in this paper, $g_c^{-2}\sim a^{-1}$ where $a$ is the lattice size. For fixed coupling $g^{-2}<g_c^{-2}$ which is not far away from the critical coupling (Otherwise, the continuum limit $a\to 0$ cannot be taken), denote ${\Sigma }_0$ be the value of $\Sigma$ at this coupling $g^{-2}$ with vanishing temperature T and chemical potential $\mu$. The gap Equation (8), which is solved exactly in the chiral limit in Ref. [18], shows that there exists a critical temperature

$T_c=\frac{{\Sigma }_0}{2\ln 2}$ such that the chiral symmetry is broken if $T<T_c$ at this coupling

$g^{-2}$ and $\mu =0$. Moreover, there is another critical chemical potential ${\mu }_c={\Sigma }_0$ such that this symmetry is broken only if $\mu <{\mu }_c$ at this coupling $g^{-2}$ and $T=0$. The mean field results predict that the first order transition only exists at $T=0$ and $\mu ={\mu }_c$ for this coupling $g^{-2}$.

For the 2 + 1d Gross-Neveu model, we first study the dependence of $\Sigma$ on the coupling g and temperature $T=1/N_t$ with vanishing chemical potential $\mu =0$. Figure 3 is the phase diagram of $\left(N_t\mathrm{,}1/g^2\right)$ for $m=0$ and $N_x=36$. We always choose $N_x=36$ to ensure the thermodynamic limit is achieved: the simulation results change very small for larger $N_x$. The marks + separate the symmetry phase $\Sigma =0$ (above marks) and the chiral symmetry broken phase $\Sigma >0$ (below marks). For fixed temperature T there is a critical coupling $g_c^{-2}$ such that $\Sigma$ decreases to zero if $1/g^2$ is increasing to $1/g_c^2$. Figure 3 shows that $1/g_c^2$ 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\left(g\right)$ such that $\Sigma$ is increasing from zero if T is decreasing from $T_c$.

Figure 4 shows the dependence of $\Sigma$ on $N_t$ for the different coupling $1/g^2$. For small $1/g^2$, e.g., $1/g^2=0.65$, $\Sigma$ changes small with the temperature.

Figure 3. Phase diagram of $\left(N_t\mathrm{,}1/g^2\right)$ for $\mu =0$, $m=0$, $N_x=36$. Below the marks + is the broken phase $\Sigma >0$.

Figure 4. $\Sigma$ versus $N_t$, $\mu =0$, $m=0$, $N_x=36$. $1/g^2=0.65,0.70,0.75,0.80,0.83,0.85,0.90,0.95,1.00$ from top to bottom.

For these range of parameters, it is in the deep chiral symmetry broken phase and we cannot obtain the chiral symmetry phase $\Sigma =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 $\frac{1}{8}$ and $\frac{1}{10}$ in lattice unit. The symmetry phase and broken phase are realized for $T>T_c\left(g\right)$ and $T<T_c\left(g\right)$, respectively.

Figure 5 shows the dependence of $\Sigma$ on $1/g^2$ at different temperature. $\Sigma$ drops continuously to 0 if $1/g^2$ is increasing to $1/g_c^2\left(T\right)$ from below, which show that the transition at the critical coupling constant $g_c\left(T\right)$ is second order. At very low temperature $T=1/N_t=1/36$, $g_c\left(T\right)$ is close to 1, which is consistent with those obtained in [19]. This is because in the limit of $N_t\mathrm{,}{N}_x\to \infty$, the gap equation at $\Sigma =0$ is reduced to

$\frac{1}{g^2}=\frac{1}{N_t{N}_x^2}{\displaystyle \underset{k}{\sum }}\frac{8}{N\left(k\right)}\approx \frac{8}{{\pi }^3}{\displaystyle {\int }_0^{\pi /2}}\text{d}k\frac{1}{{\displaystyle \underset{\mu =0,1,2}{\sum }}{\left(\cos k_{\mu }\right)}^2}=1$

The critical temperature $T_c=\frac{{\Sigma }_0}{2\ln 2}$ at the coupling $g^{-2}$ and $\mu =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_c^{-2}\approx 1$. We also choose $N_t=36$ such that it is very close to zero temperature, the value of $\Sigma$ at the zero temperature and vanishing chemical potential is ${\Sigma }_0=0.0944$. To calculate the critical temperature at this coupling, we calculate $\Sigma$ at $N_t=8,\cdots ,36$ and found that $\Sigma$ 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=\frac{{\Sigma }_0}{2\ln 2}=\frac{0.0944}{2\ln 2}=0.0680$.

Now let us study the effect of chemical potential on the chiral condensate $\Sigma$. Figure 6 shows the dependence of $\Sigma$ on the chemical potential at the different

Figure 5. $\Sigma$ versus $1/g^2$ for different $N_t$. $\mu =0$, $m=0$, $N_x=36$.

Figure 6. $\Sigma$ versus $\mu$, $m=0$, $g=1.19525$ ( $1/g^2=0.70$ ), $N_x=36$.

temperature $T=1/N_t$. $\Sigma$ drops sharply near ${\mu }_c\approx 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, ${\mu }_c={\Sigma }_0$ where ${\Sigma }_0$ is the $\Sigma$ with $\mu =0$ [18]. For the temperature $T=1/16$, ${\Sigma }_0\approx 0.47$ is slightly larger than ${\mu }_c\approx 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 [18]. Figure 7 shows the dependence of $\Sigma$ on $\mu$ for a larger $1/g^2=0.80$. Compared with Figure 6, $\Sigma$ at $\mu =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 ${\mu }_c$ is decreased if the temperature is increased. At zero temperature, the mean field exact result show the critical chemical potential ${\mu }_c$ is just the value of ${\Sigma }_0$ at the vanishing chemical potential. This is exactly recovered in Figure 7 where ${\mu }_c=0.32$ for $g^{-2}=0.80$ with $N_t=16$.

Figure 8 shows the dependence of $\Sigma$ and fermion density on the chemical potential at $1/g^2=0.7$. At low temperature $N_t=16$, $\Sigma$ drops rapidly near the critical chemical potential ${\mu }_c\approx 0.45$, and the fermion density increase very fast, which suggest $\Sigma$ and fermion density are not continuous at ${\mu }_c$ at zero temperature and thus they can be regarded as the order parameters.

For the 3 + 1d Gross-Neveu model, we also calculate the dependence of $\Sigma$ on the coupling and chemical potential at different temperature. Figure 9 shows the value of $\Sigma$ depending on the coupling for the vanishing chemical potential.

Figure 7. $\Sigma$ versus $\mu$, $m=0$, $g=1.1180$ ( $1/g^2=0.80$ ), $N_x=36$.

Figure 8. $\Sigma$ and fermion density vs $\mu$, $m=0$, $g=1.19525$ ( $1/g^2=0.70$ ), $N_x=36$.

Figure 9. $\Sigma$ versus $1/g^2$ for different $N_t$. $\mu =0$, $m=0$, $N_x=36$.

Compared to Figure 5 for the 2 + 1d model, the critical coupling becomes smaller. Moreover, the dependence of $\Sigma$ on the temperature is less sensitive. Figure 10 shows the dependence of $\Sigma$ 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 larger for the 2 + 1d model than those for the 3 + 1d model.

8. Conclusions

The staggered fermion for the Gross-Neveu model at finite density and temperature is 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 is found, which can be implemented by parallelization. This formula can also be generalized to the other space dimensions. For the odd space dimension,

Figure 10. $\Sigma$ versus $\mu$, $m=0$, ( $1/g^2=0.58$ ), $N_x=36$. Left (3 + 1d), Right (2 + 1d).

the orthogonal transformation was found [33]. The key point to find the explicit formula for the inverse matrix is to use the properties of ${\Gamma }_A$ and $B_A$ as shown in Section 4. These properties for the even number of space dimension are simpler, 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 compare the chiral condensate for the 2 + 1d and 3 + 1d case in the same range of parameters, showing that the reason for symmetry breaking and restoration can be explained by the suitable choice of the coupling, temperature and chemical potential.

Acknowledgements

Daming Li was supported by the National Science Foundation of China (No. 11271258, 11971309).

Appendix A. Proof of Properties of ${\Gamma }_A$ and $B_A$

The notations for ${\left\{A_i\right\}}_{i=0}^2$ in (27) is a little awkward. I replace $A_0$, $A_1$ and $A_2$ in (27) by $A_1$, $A_2$ and $A_3$, respectively. Thus

${\Gamma }_A={\sigma }_1^{A_1}{\sigma }_2^{A_2}{\sigma }_3^{A_3},B_A={\left(-{\sigma }_1\right)}^{A_1}{\left(-{\sigma }_2\right)}^{A_2}{\left(-{\sigma }_3\right)}^{A_3}={\left(-1\right)}^{A_1+A_2+A_3}{\Gamma }_A$ (A1)

The three Pauli matrices

$\begin{array}{l}{\left({\sigma }_1\right)}^{\alpha \beta }={\left(-1\right)}^{\beta }{\epsilon }_{\alpha \beta },\\ {\left({\sigma }_2\right)}^{\alpha \beta }=\left(-i\right){\epsilon }_{\alpha \beta },\\ {\left({\sigma }_3\right)}^{\alpha \beta }={\left(-1\right)}^{\beta -1}{\delta }_{\alpha \beta },\alpha ,\beta =1,2\end{array}$

satisfies the completeness relation

${\delta }_{\alpha a}{\delta }_{\beta b}+{\displaystyle \underset{\mu =1,2,3}{\sum }}{\sigma }_{\mu }^{\alpha a}{\sigma }_{\mu }^{*\beta b}=2{\delta }_{\alpha \beta }{\delta }_{ab}$ (A2)

We first have

$\begin{array}{c}{\displaystyle \underset{A,A_3=0}{\sum }}{\Gamma }_A^{\alpha a}{\Gamma }_A^{*\beta b}={\displaystyle \underset{A_1,A_2}{\sum }}{\left({\sigma }_1^{A_1}{\sigma }_2^{A_2}\right)}^{\alpha a}{\left({\sigma }_1^{\ast A_1}{\sigma }_2^{\ast A_2}\right)}^{\beta b}\\ ={\displaystyle \underset{A_1,A_2}{\sum }}{\left({\sigma }_1^{A_1}\right)}^{\alpha \gamma }{\left({\sigma }_2^{A_2}\right)}^{\gamma a}{\left({\sigma }_1^{\ast A_1}\right)}^{\beta {\gamma }^{\prime }}{\left({\sigma }_2^{\ast A_2}\right)}^{{\gamma }^{\prime }b}\\ ={\displaystyle \underset{A_1}{\sum }}{\left({\sigma }_1^{A_1}\right)}^{\alpha \gamma }{\left({\sigma }_1^{\ast A_1}\right)}^{\beta {\gamma }^{\prime }}{\displaystyle \underset{A_2}{\sum }}{\left({\sigma }_2^{A_2}\right)}^{\gamma a}{(}^{{\sigma }_2^{\ast A_2}}\\ =\left({\delta }_{\alpha \gamma }{\delta }_{\beta {\gamma }^{\prime }}+{\sigma }_1^{\alpha \gamma }{\sigma }_1^{\ast \beta {\gamma }^{\prime }}\right)\left({\delta }_{\gamma a}{\delta }_{{\gamma }^{\prime }b}+{\sigma }_2^{\gamma a}{\sigma }_2^{\ast {\gamma }^{\prime }b}\right)\\ ={\delta }_{\alpha a}{\delta }_{\beta b}+{\sigma }_2^{\alpha a}{\sigma }_2^{\ast \beta b}+{\sigma }_1^{\alpha a}{\sigma }_1^{\ast \beta b}+{\left({\sigma }_1{\sigma }_2\right)}^{\alpha a}{\left({\sigma }_1^{\ast }{\sigma }_2^{\ast }\right)}^{\beta b}\\ ={\delta }_{\alpha a}{\delta }_{\beta b}+{\sigma }_2^{\alpha a}{\sigma }_2^{\ast \beta b}+{\sigma }_1^{\alpha a}{\sigma }_1^{\ast \beta b}+{\sigma }_3^{\alpha a}{\sigma }_3^{\ast \beta b}\end{array}$

$=2{\delta }_{\alpha \beta }{\delta }_{ab}$ by (A2) (A3)

which is also valid if $\left(\mathrm{1,2}\right)$ is replaced by $\left(\mathrm{1,3}\right)$ or $\left(\mathrm{2,3}\right)$. Secondly,

$\begin{array}{c}{\displaystyle \underset{A}{\sum }}{\Gamma }_A^{\alpha a}{\Gamma }_A^{\ast \beta b}={\displaystyle \underset{A_1,A_2,A_3}{\sum }}{\left({\sigma }_1^{A_1}{\sigma }_2^{A_2}{\sigma }_3^{A_3}\right)}^{\alpha a}{\left({\sigma }_1^{\ast A_1}{\sigma }_2^{\ast A_2}{\sigma }_3^{\ast A_3}\right)}^{\beta b}\\ ={\displaystyle \underset{A_1,A_2,A_3}{\sum }}{\left({\sigma }_1^{A_1}{\sigma }_2^{A_2}\right)}^{\alpha t}{\left({\sigma }_3^{A_3}\right)}^{ta}{\left({\sigma }_1^{\ast A_1}{\sigma }_2^{\ast A_2}\right)}^{\beta t^{\prime }}{\left({\sigma }_3^{\ast A_3}\right)}^{t^{\prime }b}\\ =2{\delta }_{\alpha \beta }{\delta }_{t{t}^{\prime }}\left({\delta }_{ta}{\delta }_{t^{\prime }b}+{\left(-1\right)}^{a+b}{\delta }_{ta}{\delta }_{t^{\prime }b}\right)\\ =4{\delta }_{\alpha \beta }{\delta }_{ab}\end{array}$

Inserting $B_A={\left(-1\right)}^{A_1+A_2+A_3}{\Gamma }_A$ in the above equality, we have ${\displaystyle {\sum }_A}{B}_A^{\alpha a}{B}_A^{\ast \beta b}=4{\delta }_{\alpha \beta }{\delta }_{ab}$.

$\begin{array}{c}{\displaystyle \underset{A}{\sum }}{\Gamma }_A^{\alpha a}{B}_A^{\ast \beta b}={\displaystyle \underset{A_1,A_2,A_3}{\sum }}{\left(-1\right)}^{A_1+A_2+A_3}{\left({\sigma }_1^{A_1}{\sigma }_2^{A_2}{\sigma }_3^{A_3}\right)}^{\alpha a}{\left({\sigma }_1^{\ast A_1}{\sigma }_2^{\ast A_2}{\sigma }_3^{\ast A_3}\right)}^{\beta b}\\ ={\displaystyle \underset{A_1,A_2,A_3}{\sum }}{\left(-1\right)}^{A_1+A_2}{\left(-1\right)}^{A_3}{\left({\sigma }_1^{A_1}{\sigma }_2^{A_2}\right)}^{\alpha t}{\left({\sigma }_3^{A_3}\right)}^{ta}{\left({\sigma }_1^{\ast A_1}{\sigma }_2^{\ast A_2}\right)}^{\beta t^{\prime }}{\left({\sigma }_3^{\ast A_3}\right)}^{t^{\prime }b}\\ =\left({\delta }_{\alpha t}{\delta }_{\beta t^{\prime }}-{\sigma }_2^{\alpha t}{\sigma }_2^{\ast \beta t^{\prime }}-{\sigma }_1^{\alpha t}{\sigma }_1^{\ast \beta t^{\prime }}+{\sigma }_3^{\alpha t}{\sigma }_3^{\ast \beta t^{\prime }}\right)\left({\delta }_{ta}{\delta }_{t^{\prime }b}-{\left(-1\right)}^{a+b}{\delta }_{ta}{\delta }_{t^{\prime }b}\right)\\ =\left({\delta }_{\alpha a}{\delta }_{\beta b}-{\sigma }_2^{\alpha a}{\sigma }_2^{\ast \beta b}-{\sigma }_1^{\alpha a}{\sigma }_1^{\ast \beta b}+{\sigma }_3^{\alpha a}{\sigma }_3^{\ast \beta b}\right)\left(1-{\left(-1\right)}^{a+b}\right)=0\end{array}$

where in the last equality we used

$\begin{array}{l}{\delta }_{\alpha a}{\delta }_{\beta b}-{\sigma }_2^{\alpha a}{\sigma }_2^{\ast \beta b}-{\sigma }_1^{\alpha a}{\sigma }_1^{\ast \beta b}+{\sigma }_3^{\alpha a}{\sigma }_3^{\ast \beta b}\\ ={\delta }_{\alpha a}{\delta }_{\beta b}-{\epsilon }_{\alpha a}{\epsilon }_{\beta b}-{\left(-1\right)}^{a+b}{\epsilon }_{\alpha a}{\epsilon }_{\beta b}+{\left(-1\right)}^{a+b}{\delta }_{\alpha a}{\delta }_{\beta b}\\ =\left(1+{\left(-1\right)}^{a+b}\right)\left({\delta }_{\alpha a}{\delta }_{\beta b}-{\epsilon }_{\alpha a}{\epsilon }_{\beta b}\right)=\mathrm{0,}\text{if}a\ne b\end{array}$

To prove that

${\displaystyle \underset{A,A_{\mu }=1}{\sum }}{\Gamma }_A^{\alpha a}{\left({\sigma }_{\mu }^{\ast }{B}_A^{\ast }\right)}^{{\alpha }^{\prime }{a}^{\prime }}=-2{\left({\sigma }_{\mu }^{\ast }\right)}^{a{a}^{\prime }}{\delta }_{\alpha {\alpha }^{\prime }}$ (A4)

we want to prove that

${\displaystyle \underset{A,A_{\mu }=1}{\sum }}{\sigma }_{\mu }^{\ast ba}{\Gamma }_A^{\alpha a}{\left({\sigma }_{\mu }^{\ast }{B}_A^{\ast }\right)}^{{\alpha }^{\prime }{a}^{\prime }}=-2{\sigma }_{\mu }^{\ast ba}{\left({\sigma }_{\mu }^{\ast }\right)}^{a{a}^{\prime }}{\delta }_{\alpha {\alpha }^{\prime }}$

i.e.,

${\displaystyle \underset{A,A_{\mu }=1}{\sum }}{\left({\Gamma }_A{\sigma }_{\mu }\right)}^{\alpha b}{\left({\sigma }_{\mu }^{\ast }{B}_A^{\ast }\right)}^{{\alpha }^{\prime }{a}^{\prime }}=-2{\delta }_{b{a}^{\prime }}{\delta }_{\alpha {\alpha }^{\prime }}$

This is obvious since the left hand side is

$\begin{array}{l}{\displaystyle \underset{A\mathrm{,}{A}_{\mu }=1}{\sum }}{\left({\sigma }_1^{A_1}\cdots {\sigma }_{\mu -1}^{A_{\mu -1}}{\sigma }_{\mu +1}^{A_{\mu +1}}\cdots {\sigma }_3^{A_3}\right)}^{\alpha b}{\left(-1\right)}^{A_{\mu +1}+\cdots +A_3}\\ \times {\left({\sigma }_1^{\ast A_1}\cdots {\sigma }_{\mu -1}^{\ast A_{\mu -1}}{\sigma }_{\mu +1}^{\ast A_{\mu +1}}\cdots {\sigma }_3^{\ast A_3}\right)}^{{\alpha }^{\prime }{a}^{\prime }}{\left(-1\right)}^{A_1+\cdots +A_{\mu -1}}{\left(-1\right)}^{A_1+A_2+A_3}\\ ={\displaystyle \underset{A\mathrm{,}{A}_{\mu }=1}{\sum }}{\left(-1\right)}^{A_{\mu }}{\left({\sigma }_1^{A_1}\cdots {\sigma }_{\mu -1}^{A_{\mu -1}}{\sigma }_{\mu +1}^{A_{\mu +1}}\cdots {\sigma }_3^{A_3}\right)}^{\alpha b}{\left({\sigma }_1^{\ast A_1}\cdots {\sigma }_{\mu -1}^{\ast A_{\mu -1}}{\sigma }_{\mu +1}^{\ast A_{\mu +1}}\cdots {\sigma }_3^{\ast A_3}\right)}^{{\alpha }^{\prime }{a}^{\prime }}\end{array}$

$=-2{\delta }_{b{a}^{\prime }}{\delta }_{\alpha {\alpha }^{\prime }}$, by (A3) if $\mu =3$ (A5)

Similarly, (A4) is also valid if $A_{\mu }=1$ and -2 are replaced by $A_{\mu }=0$ and +2, respectively. This is because the sign ${\left(-1\right)}^{A_{\mu }}=-1$ in (A5) is replaced by ${\left(-1\right)}^{A_{\mu }}=+1$. Obviously,

${\Gamma }_{A\pm \hat{\mu }}={\eta }_{\mu }\left(A\right){\sigma }_{\mu }{\Gamma }_A,\mu =1,2,3$

For example, $\mu =2$,

${\Gamma }_{A\pm \hat{2}}={\sigma }_1^{A_1}{\sigma }_2^{A_2\pm 1}{\sigma }_3^{A_3}={\sigma }_1^{A_1}{\sigma }_2^{A_2+1}{\sigma }_3^{A_3}={\eta }_2\left(A\right){\sigma }_2{\Gamma }_A,{\eta }_2\left(A\right)={\left(-1\right)}^{A_1}$

Finally, we have

$\frac{1}{4}\text{Tr}\left({\Gamma }_A^{†}{\Gamma }_{A^{\prime }}+B_A^{†}{B}_{A^{\prime }}\right)={\delta }_{A{A}^{\prime }}$

since the left hand side is

$\begin{array}{l}\frac{1}{4}\text{Tr}\left[{\left(\begin{array}{cc}{\Gamma }_A& \\ & B_A\end{array}\right)}^{†}\left(\begin{array}{cc}{\Gamma }_{A^{\prime }}& \\ & B_{A^{\prime }}\end{array}\right)\right]\\ =\frac{1}{4}\text{Tr}\left(\left({\gamma }_3^{A_3}{\gamma }_2^{A_2}{\gamma }_1^{A_1}\right)\left({\gamma }_1^{{A^{\prime }}_1}{\gamma }_2^{{A^{\prime }}_2}{\gamma }_3^{{A^{\prime }}_3}\right)\right)\\ =\frac{1}{4}{\left(-1\right)}^{\left(A_1+{A^{\prime }}_1\right)\left(A_2+A_3\right)+\left(A_2+{A^{\prime }}_2\right)A_3}\text{Tr}\left({\gamma }_1^{A_1+{A^{\prime }}_1}{\gamma }_2^{A_2+{A^{\prime }}_2}{\gamma }_3^{A_3+{A^{\prime }}_3}\right)={\delta }_{A{A}^{\prime }}\end{array}$

where we used

${\gamma }_{\mu }^i{\gamma }_{\nu }^j={\left(-1\right)}^{ij}{\gamma }_{\nu }^j{\gamma }_{\mu }^i\mathrm{,}\mu \ne \nu \mathrm{,}i\mathrm{,}j=\mathrm{0,1,2}$

$\text{Tr}\left({\gamma }_{\mu }\right)=\mathrm{0,}\text{Tr}\left({\gamma }_{\mu }{\gamma }_{\nu }\right)=\mathrm{0,}\mu \ne \nu \mathrm{,}\text{Tr}\left({\gamma }_1{\gamma }_2{\gamma }_3\right)=0$

Here the we define ${\gamma }_{\mu }=\left(\begin{array}{cc}{\sigma }_{\mu }& 0\\ 0& -{\sigma }_{\mu }\end{array}\right)\left(\mu =1,2,3\right)$.

Appendix B: The Derivation of (38)

The derivation of (38) is similar to the calculation of

$\frac{1}{2}{\displaystyle {\sum }_x}{\eta }_{\mu }\left(x\right)\bar{\chi }\left(x\right)\left(\chi \left(x+\hat{\mu }\right)-\chi \left(x-\hat{\mu }\right)\right)$.

$\begin{array}{l}\frac{1}{2}{\displaystyle \underset{x}{\sum }}\bar{\chi }\left(x\right)\left(\chi \left(x+\hat{0}\right)+\chi \left(x-\hat{0}\right)\right)\\ =\frac{1}{2}{\displaystyle \underset{A\mathrm{,}{A}^{\prime }\mathrm{,}Y}{\sum }}\bar{\chi }\left(A\mathrm{,}Y\right)({\delta }_{A+\hat{0}\mathrm{,}{A}^{\prime }}\left(\chi \left(A^{\prime }\mathrm{,}Y\right)+\chi \left(A^{\prime }\mathrm{,}Y-\hat{0}\right)\right)\\ +{\delta }_{A-\hat{0}\mathrm{,}{A}^{\prime }}\left(\chi \left(A^{\prime }\mathrm{,}Y+\hat{0}\right)+\chi \left(A^{\prime }\mathrm{,}Y\right)\right))\\ =\frac{1}{2}{\displaystyle \underset{A\mathrm{,}{A}^{\prime }\mathrm{,}Y}{\sum }}\bar{\chi }\left(A\mathrm{,}Y\right)\left(\frac{{\delta }_{A-\hat{0}\mathrm{,}{A}^{\prime }}-{\delta }_{A+\hat{0}\mathrm{,}{A}^{\prime }}}{2}{\partial }_0\chi \left(A^{\prime }\mathrm{,}Y\right)+\frac{{\delta }_{A-\hat{0}\mathrm{,}{A}^{\prime }}+{\delta }_{A+\hat{0}\mathrm{,}{A}^{\prime }}}{2}\delta \chi \left(A^{\prime }\mathrm{,}Y\right)\right)\end{array}$

$\begin{array}{l}=\frac{1}{2}{\displaystyle \underset{A\mathrm{,}{A}^{\prime }\mathrm{,}Y}{\sum }}\sqrt{2}{\displaystyle \underset{\alpha \mathrm{,}a}{\sum }}\left({\bar{u}}^{\alpha a}\left(Y\right){\Gamma }_A^{\alpha a}+{\bar{d}}^{\alpha a}\left(Y\right)B_A^{\alpha a}\right)\\ \times \{\frac{{\delta }_{A-\hat{0}\mathrm{,}{A}^{\prime }}-{\delta }_{A+\hat{0}\mathrm{,}{A}^{\prime }}}{2}\sqrt{2}{\displaystyle \underset{{\alpha }^{\prime },a^{\prime }}{\sum }}\left({\Gamma }_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{\partial }_0{u}^{{\alpha }^{\prime }{a}^{\prime }}\left(Y\right)+B_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{\partial }_0{d}^{{\alpha }^{\prime }{a}^{\prime }}\left(Y\right)\right)\\ +\frac{{\delta }_{A-\hat{0}\mathrm{,}{A}^{\prime }}+{\delta }_{A+\hat{0}\mathrm{,}{A}^{\prime }}}{2}\sqrt{2}{\displaystyle \underset{{\alpha }^{\prime },a^{\prime }}{\sum }}\left({\Gamma }_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}\delta u^{{\alpha }^{\prime }{a}^{\prime }}\left(Y\right)+B_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}\delta d^{{\alpha }^{\prime }{a}^{\prime }}\left(Y\right)\right)\}\end{array}$

$\begin{array}{l}=2{\displaystyle \underset{Y}{\sum }}({\bar{u}}^{\alpha a}\left(Y\right){\left({\sigma }_1^{\mathrm{<mo>\; *\; </mo>}}\right)}^{a{a}^{\prime }}{\partial }_0{d}^{\alpha a^{\prime }}\left(Y\right)+{\bar{d}}^{\alpha a}\left(Y\right){\left(-{\sigma }_1^{\mathrm{<mo>\; *\; </mo>}}\right)}^{a{a}^{\prime }}{\partial }_0{u}^{\alpha a^{\prime }}\left(Y\right)\\ +{\bar{u}}^{\alpha a}\left(Y\right){\left({\sigma }_1\right)}^{\alpha {\alpha }^{\prime }}\delta u^{{\alpha }^{\prime }a}\left(Y\right)+{\bar{d}}^{\alpha a}\left(Y\right){\left(-{\sigma }_1\right)}^{\alpha {\alpha }^{\prime }}\delta d^{{\alpha }^{\prime }a}\left(Y\right))\\ =8{\displaystyle \underset{Y}{\sum }}\left[\bar{q}\left(Y\right)\left(i{\gamma }_3\otimes {\sigma }_1^{\mathrm{<mo>\; *\; </mo>}}\right)\frac{{\partial }_{0}q\left(Y\right)}{4}+\bar{q}\left(Y\right)\left({\gamma }_0\otimes {\mathbb{I}}_2\right)\frac{\delta q\left(Y\right)}{4}\right]\\ =8{\displaystyle \underset{k}{\sum }}[\bar{q}\left(k\right)\left(i{\gamma }_3\otimes {\sigma }_1^{\mathrm{<mo>\; *\; </mo>}}\right)i{2}^{-1}\sin \left(2k_0\right)q\left(k\right)\\ +\bar{q}\left(k\right)\left({\gamma }_0\otimes {\mathbb{I}}_2\right)2^{-1}\left[\cos \left(2k_0\right)+1\right]q\left(k\right)]\end{array}$

where

$\delta q\left(Y\right)=q\left(Y+\hat{0}\right)+2q\left(Y\right)+q\left(Y-\hat{0}\right)$

In the fourth equality, we used the formula like

$\begin{array}{l}{\displaystyle \underset{A,A^{\prime }}{\sum }}{\Gamma }_A^{\alpha a}{B}_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}\left({\delta }_{A-\hat{0},A^{\prime }}-{\delta }_{A+\hat{0},A^{\prime }}\right)\\ ={\displaystyle \underset{A,A_0=1}{\sum }}{\Gamma }_A^{\alpha a}{B}_{A-\hat{0}}^{\ast {\alpha }^{\prime }{a}^{\prime }}-{\displaystyle \underset{A,A_0=0}{\sum }}{\Gamma }_A^{\alpha a}{B}_{A+\hat{0}}^{\ast {\alpha }^{\prime }{a}^{\prime }}\\ ={\displaystyle \underset{A,A_0=1}{\sum }}{\Gamma }_A^{\alpha a}{\left(-{\sigma }_1{B}_A\right)}^{\ast {\alpha }^{\prime }{a}^{\prime }}-{\displaystyle \underset{A,A_0=0}{\sum }}{\Gamma }_A^{\alpha a}{\left(-{\sigma }_1{B}_A\right)}^{\ast {\alpha }^{\prime }{a}^{\prime }}\\ =4{\sigma }_1^{\ast a{a}^{\prime }}{\delta }_{\alpha {\alpha }^{\prime }}\end{array}$

Appendix C. The Derivation of (44)-(46)

First,

$\begin{array}{l}{\displaystyle \underset{x}{\sum }}\left(D_{x+\hat{0}\mathrm{,}x}^{-1}{s}_x^1+D_{x-\hat{0}\mathrm{,}x}^{-1}{s}_x^2\right)\\ =-\frac{{\displaystyle \int {\text{e}}^{-\bar{\chi }D\chi }}{\displaystyle {\sum }_x\bar{\chi }\left(x\right)}\left[\chi \left(x+\hat{0}\right)+\chi \left(x-\hat{0}\right)\right]}{{\displaystyle \int {\text{e}}^{-\bar{\chi }D\chi }}}\end{array}$

$=-\frac{{\displaystyle \int {\text{e}}^{-{\displaystyle {\sum }_k\bar{q}\left(k\right)8D\left(k\right)q\left(k\right)}}}16{\displaystyle {\sum }_k\bar{q}\left(k\right)A_{+}\left(k\right)q\left(k\right)}}{{\displaystyle \int {\text{e}}^{-{\displaystyle {\sum }_k\bar{q}\left(k\right)8D\left(k\right)q\left(k\right)}}}}$ by (38) (39)

$\begin{array}{l}=-16{\displaystyle \underset{k}{\sum }}\frac{{\displaystyle \int {\text{e}}^{-\bar{q}\left(k\right)8D\left(k\right)q\left(k\right)}\bar{q}\left(k\right)A_{+}\left(k\right)q\left(k\right)}}{{\displaystyle \int {\text{e}}^{-\bar{q}\left(k\right)8D\left(k\right)q\left(k\right)}}}\\ =16{\displaystyle \underset{k}{\sum }}\text{tr}\left[{\left(8D\left(k\right)\right)}^{-1}{A}_{+}\left(k\right)\right]\\ =2{\displaystyle \underset{k}{\sum }}\text{tr}\left[D{\left(k\right)}^{-1}{A}_{+}\left(k\right)\right]\end{array}$

$\begin{array}{l}={\displaystyle \underset{k}{\sum }}\frac{2}{N\left(k\right)}\text{tr}\{\left[m-{\displaystyle \underset{\mu =\mathrm{0,1,2}}{\sum }}\left({\gamma }_{\mu }\otimes {\mathbb{I}}_2\right)a_{\mu }-{\displaystyle \underset{c=\mathrm{1,2,3}}{\sum }}\left({\gamma }_3\otimes {\sigma }_c^{\mathrm{<mo>\; *\; </mo>}}\right)b_c\right]\\ \times \left[\left(i{\gamma }_3\otimes {\sigma }_1^{\mathrm{<mo>\; *\; </mo>}}\right)i{2}^{-1}\sin \left(2k_0\right)+\left({\gamma }_0\otimes {\mathbb{I}}_2\right)2^{-1}\left[\cos \left(2k_0\right)+1\right]\right]\}\end{array}$ by (41)

$\begin{array}{l}={\displaystyle \underset{k}{\sum }}\frac{2}{N\left(k\right)}\text{tr}\left\{\left({\mathbb{I}}_4\otimes {\mathbb{I}}_2\right)\left(b_1{2}^{-1}\sin 2k_0-a_0{2}^{-1}\left(\cos 2k_0+1\right)\right)\right\}\\ =8{\displaystyle \underset{k}{\sum }}\frac{b_1\sin 2k_0-a_0\left(\cos 2k_0+1\right)}{N(\; k\; )}\end{array}$

Similarly,

$\begin{array}{l}{\displaystyle \underset{x}{\sum }}\left(D_{x+\hat{0}\mathrm{,}x}^{-1}{s}_x^1-D_{x-\hat{0}\mathrm{,}x}^{-1}{s}_x^2\right)\\ ={\displaystyle \underset{k}{\sum }}\frac{2}{N\left(k\right)}\text{tr}\{\left[m-{\displaystyle \underset{\mu =\mathrm{0,1,2}}{\sum }}\left({\gamma }_{\mu }\otimes {\mathbb{I}}_2\right)a_{\mu }-{\displaystyle \underset{c=\mathrm{1,2,3}}{\sum }}\left({\gamma }_3\otimes {\sigma }_c^{\mathrm{<mo>\; *\; </mo>}}\right)b_c\right]\\ \begin{array}{c}\\ \end{array}\times \left[\left({\gamma }_0\otimes {\mathbb{I}}_2\right)i{2}^{-1}\sin \left(2k_0\right)+\left(i{\gamma }_3\otimes {\sigma }_1^{\mathrm{<mo>\; *\; </mo>}}\right)2^{-1}\left[\cos \left(2k_0\right)-1\right]\right]\}\\ ={\displaystyle \underset{k}{\sum }}\frac{2}{N\left(k\right)}\text{tr}\left\{\left({\mathbb{I}}_4\otimes {\mathbb{I}}_2\right)\left(-a_{0}i{2}^{-1}\sin 2k_0-b_{1}i{2}^{-1}\left(\cos 2k_0-1\right)\right)\right\}\\ =\left(-8i\right){\displaystyle \underset{k}{\sum }}\frac{a_0\sin 2k_0+b_1\left(\cos 2k_0-1\right)}{N(\; k\; )}\end{array}$

The inverse matrix of D in (20) can be calculated as follows

$\begin{array}{l}{D}_{x^{\prime }\mathrm{,}x}^{-1}=-\frac{{\displaystyle \int {\text{e}}^{-\bar{\chi }D\chi }}\bar{\chi }\left(x\right)\chi \left(x^{\prime }\right)}{{\displaystyle \int {\text{e}}^{-\bar{\chi }D\chi }}}\\ =-2{\displaystyle \underset{\alpha \mathrm{,}a\mathrm{,}{\alpha }^{\prime }\mathrm{,}{a}^{\prime }}{\sum }}\frac{{\displaystyle \int {\text{e}}^{-\bar{q}8Dq}}\left[{\bar{u}}^{\alpha a}\left(Y\right){\Gamma }_A^{\alpha a}+{\bar{d}}^{\alpha a}\left(Y\right)B_A^{\alpha a}\right]\left[{\Gamma }_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{u}^{{\alpha }^{\prime }{a}^{\prime }}\left(Y^{\prime }\right)+B_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{d}^{{\alpha }^{\prime }{a}^{\prime }}\left(Y^{\prime }\right)\right]}{{\displaystyle \int {\text{e}}^{-\bar{q}8Dq}}}\\ =-2{\displaystyle \underset{\alpha \mathrm{,}a\mathrm{,}{\alpha }^{\prime }\mathrm{,}{a}^{\prime }}{\sum }}[{\Gamma }_A^{\alpha a}{\Gamma }_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}\frac{{\displaystyle \int {\text{e}}^{-\bar{q}8Dq}}{\bar{q}}_1^{\alpha a}\left(Y\right)q_1^{{\alpha }^{\prime }{a}^{\prime }}\left(Y^{\prime }\right)}{{\displaystyle \int {\text{e}}^{-\bar{q}8Dq}}}+{\Gamma }_A^{\alpha a}{B}_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}\frac{{\displaystyle \int {\text{e}}^{-\bar{q}8Dq}}{\bar{q}}_1^{\alpha a}\left(Y\right)q_2^{{\alpha }^{\prime }{a}^{\prime }}\left(Y^{\prime }\right)}{{\displaystyle \int {\text{e}}^{-\bar{q}8Dq}}}\end{array}$

$\begin{array}{l}+B_A^{\alpha a}{\Gamma }_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}\frac{{\displaystyle \int {\text{e}}^{-\bar{q}8Dq}}{\bar{q}}_2^{\alpha a}\left(Y\right)q_1^{{\alpha }^{\prime }{a}^{\prime }}\left(Y^{\prime }\right)}{{\displaystyle \int {\text{e}}^{-\bar{q}8Dq}}}+B_A^{\alpha a}{B}_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}\frac{{\displaystyle \int {\text{e}}^{-\bar{q}8Dq}}{\bar{q}}_2^{\alpha a}\left(Y\right)q_2^{{\alpha }^{\prime }{a}^{\prime }}\left(Y^{\prime }\right)}{{\displaystyle \int {\text{e}}^{-\bar{q}8Dq}}}]\\ =\frac{1}{4}{\displaystyle \underset{\alpha \mathrm{,}a\mathrm{,}{\alpha }^{\prime }\mathrm{,}{a}^{\prime }}{\sum }}[{\Gamma }_A^{\alpha a}{\Gamma }_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{D}_{\left(Y^{\prime }{a}^{\prime }{\alpha }^{\prime }\mathrm{1\; <mo>\; ;\; </mo>}Ya\alpha 1\right)}^{-1}+{\Gamma }_A^{\alpha a}{B}_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{D}_{\left(Y^{\prime }{a}^{\prime }{\alpha }^{\prime }\mathrm{2\; <mo>\; ;\; </mo>}Ya\alpha 1\right)}^{-1}\\ +B_A^{\alpha a}{\Gamma }_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{D}_{\left(Y^{\prime }{a}^{\prime }{\alpha }^{\prime }\mathrm{1\; <mo>\; ;\; </mo>}Ya\alpha 2\right)}^{-1}+B_A^{\alpha a}{B}_{A^{\prime }}^{\ast {\alpha }^{\prime }{a}^{\prime }}{D}_{\left(Y^{\prime }{a}^{\prime }{\alpha }^{\prime }\mathrm{2\; <mo>\; ;\; </mo>}Ya\alpha 2\right)}^{-1}]\end{array}$

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	de Forcrand, P. (2009) PoS, LAT2009, 10. https://doi.org/10.22323/1.091.0010
[2]	Fukushima, K. and Hatsuda, T. (2011) Reports on Progress in Physics, 74, Article ID: 014001. https://doi.org/10.1088/0034-4885/74/1/014001
[3]	Bazavov, A., et al. (2012) Physical Review D, 85, Article ID: 054503. https://doi.org/10.1103/PhysRevD.85.054503
[4]	Bhattacharya, T., et al. (2014) Physical Review Letters, 113, Article ID: 082001. https://doi.org/10.1103/PhysRevLett.113.082001
[5]	Schmidt, C. and Sharma, S. (2017) Journal of Physics G: Nuclear and Particle Physics, 44, Article ID: 104002. https://doi.org/10.1088/1361-6471/aa824a
[6]	Bazavov, A., et al. (2019) Physics Letters B, 795, 15. https://doi.org/10.1016/j.physletb.2019.05.013
[7]	Fischer, C.S. (2019) Progress in Particle and Nuclear Physics, 105, 1. https://doi.org/10.1016/j.ppnp.2019.01.002
[8]	Fu, W.J., Pawlowski, J.M. and Rennecke, F. (2020) Physical Review D, 101, Article ID: 054032. https://doi.org/10.1103/PhysRevD.101.054032
[9]	Guenther, J.N. (2021) The European Physical Journal A, 57, 136. https://doi.org/10.1140/epja/s10050-021-00354-6
[10]	Cotter, S., Giudice, P., Hands, S. and Skullerud, J.-I. (2013) Physical Review D, 87, Article ID: 034507. https://doi.org/10.1103/PhysRevD.87.034507
[11]	Li, D.M. (2016) Physical Review D, 94, Article ID: 114501. https://doi.org/10.1103/PhysRevD.94.114501
[12]	Banuls, M.C. and Cichy, K. (2020) Reports on Progress in Physics, 83, Article ID: 024401. https://doi.org/10.1088/1361-6633/ab6311
[13]	Hauschild, J. and Pollmann, F. (2018) SciPost Physics Lecture Notes, 5, 1.
[14]	Verstraete, F. and Ignacio Cirac, J. (2010) Physical Review Letters, 104, Article ID: 190405. https://doi.org/10.1103/PhysRevLett.104.190405
[15]	Haegeman, J., Ignacio Cirac, J., Osborne, T.J. and Verstraete, F. (2013) Physical Review B, 88, Article ID: 085118. https://doi.org/10.1103/PhysRevB.88.085118
[16]	Rosenstein, B., Warr, B.J. and Park, S.H. (1991) Physics Reports, 205, 59. https://doi.org/10.1016/0370-1573(91)90129-A
[17]	Rosenstein, B., Warr, B.J. and Park, S.H. (1989) Physical Review D, 39, 3088. https://doi.org/10.1103/PhysRevD.39.3088
[18]	Hands, S., Kocic, A. and Kogut, J.B. (1993) Nuclear Physics B, 390, 355. https://doi.org/10.1016/0550-3213(93)90460-7
[19]	Hands, S., Kocić, A. and Kogut, J.B. (1993) Annals of Physics, 224, 29. https://doi.org/10.1006/aphy.1993.1039
[20]	Kogut, J.B. and Strouthos, C.G. (2001) Physical Review D, 63, Article ID: 054502. https://doi.org/10.1103/PhysRevD.63.054502
[21]	Wolff, U. (1985) Physics Letters B, 157, 303. https://doi.org/10.1016/0370-2693(85)90671-9
[22]	Thies, M. and Urlichs, K. (2003) Physical Review D, 67, Article ID: 125015. https://doi.org/10.1103/PhysRevD.67.125015
[23]	Schnetz, O., Thies, M. and Urlichs, K. (2004) Annals of Physics (Amsterdam), 314, 425. https://doi.org/10.1016/j.aop.2004.06.009
[24]	Thies, M. (2006) Journal of Physics A, 39, 12707. https://doi.org/10.1088/0305-4470/39/41/S04
[25]	de Forcrand, P. and Wenger, U. (2006) PoS, LAT2009, 152.
[26]	Wagner, M. (2007) Physical Review D, 76, Article ID: 076002. https://doi.org/10.1103/PhysRevD.76.076002
[27]	Lenz, J., Pannullo, L., Wagner, M., Wellegehausen, B. and Wipf, A. (2020) Physical Review D, 101, Article ID: 094512. https://doi.org/10.1103/PhysRevD.101.094512
[28]	Lenz, J.J., Pannullo, L., Wagner, M., Wellegehausen, B. and Wipf, A. (2020) Physical Review D, 102, Article ID: 114501. https://doi.org/10.1103/PhysRevD.102.114501
[29]	Castorina, P., Mazza, M. and Zappala, D. (2003) Physics Letters B, 567, 31-38. https://doi.org/10.1016/j.physletb.2003.06.005
[30]	Buballa, M., Kurth, L., Wagner, M. and Winstel, M. (2021) Physical Review D, 103, Article ID: 034503. https://doi.org/10.1103/PhysRevD.103.034503
[31]	Hands, S. (2016) JHEP, 11, 15. https://doi.org/10.1007/JHEP11(2016)015
[32]	Rothe, H.J. (2005) Lattice Gauge Theory, an Introduction. World Scientific Lecture Notes in Physics Vol. 74, World Scientific Publishing, Singapore. https://doi.org/10.1142/5674
[33]	Cohen, Y., Elittzur, S. and Rabinovici, E. (1983) Nuclear Physics B, 220, 102. https://doi.org/10.1016/0550-3213(83)90136-0
[34]	Burden, C.J. and Burkitt, A.N. (1987) Europhysics Letters, 3, 545. https://doi.org/10.1209/0295-5075/3/5/006