text:
# Hamiltonian lattice QCD at finite density: equation of state in the strong coupling limit

## I Introduction

Lattice gauge theory is currently the only known method of solving Quantum Chromodynamics (QCD) from first principles. It has developed sufficiently to be able to calculate a broad range of low and intermediate energy hadronic observables from ground state hadron masses to pion–nucleon scattering lengths. In addition, lattice studies of QCD at finite temperature ($`T`$), especially its spectacular success in demonstrating the deconfinement phase transition, have been invaluable and continue to play an important role in the search for the quark–gluon plasma.

However, as is well known, one of the outstanding problems in lattice gauge theory is the consistent implimentation of chemical potential in numerical simulations M. Creutz (2001). Progress in lattice QCD calculations at finite chemical potential ($`\mu `$) with dynamical quarks has been hindered by the presence of the complex fermion determinant which renders standard Monte Carlo techniques useless. In fact, currently there is only one numerical method of simulating finite density QCD with three colors at zero temperature. This method is based on the Monomer–Dimer–Polymer algorithm developed by Karsch and Mütter F. Karsch and K.–H. Mütter (1989). However its applicability is limited to the strong coupling limit, and furthermore a recent study R. Aloisio, V. Azcoiti, G. Di Carlo, A. Galante and A.F. Grillo (2000) indicates that this algorithm might not be reliable for studying the chiral phase transition at finite density. Therefore even a qualitative description of finite density lattice QCD is welcome.

One method of studying finite density QCD on the lattice is to invoke the strong coupling approximation where analytical methods are applicable. Although far from the realistic continuum limit, the strong coupling approximation has played an important role in the development of QCD lattice gauge theory from its very inception. In the renowned paper by Wilson K.G. Wilson (1974) this approximation was invoked to demonstrate quark confinement on the Euclidean space–time lattice. Soon thereafter Kogut and Susskind J. Kogut and L. Susskind (1975) formulated the Hamiltonian lattice gauge theory and concluded that in the strong coupling limit the quark dynamics is best described by a collection of non–Abelian electric flux tubes with quarks attached at their ends. This was followed by the work of Baluni and Willemsen who used a variant of the Kogut–Susskind formalism to demonstrate quantitatively that dynamical chiral symmetry breaking indeed takes places in lattice QCD at strong coupling V. Baluni and J.F. Willemsen (1976). Finally, calculations by Kogut, Pearson and Shigemitsu J. Kogut, R.B. Pearson and J. Shigemitsu (1979) and by Creutz M. Creutz (1979) suggesting the absence of a phase transition between the strong and weak coupling regimes of QCD motivated numerous studies using the strong coupling approximation.

Strong coupling QCD at finite $`T`$ and/or $`\mu `$ has previously been studied analytically both in the Euclidean P.H. Damgaard, D. Hochberg and N. Kawamoto (1985); E.–M. Ilgenfritz and J. Kripfganz (1985); N. Bilić, K. Demeterfi and B. Petersson (1992) and in the Hamiltonian A. Patel (1984); C.P. Van Den Doel (1984); A. Le Yaouanc et al. (1988); E.B. Gregory, S.–H. Guo, H. Kröger and X.–Q. Luo (2000) formulations. One of the main objectives of these studies was to investigate the nature of chiral phase transition at finite temperature and density. In each study this was accomplished by constructing some effective action or Hamiltonian for strongly coupled lattice QCD using Kogut–Susskind fermions. These effective descriptions involve the introduction of composite meson and baryon fields which are treated in the mean field approximation.111The Monomer–Dimer–Polymer algorithm F. Karsch and K.–H. Mütter (1989) also uses composite meson and baryons fields. The consensus is that at zero or low temperatures strong coupling QCD undergoes a first order chiral phase transition from the broken symmetry phase below a critical chemical potential $`\mu _\text{C}`$ to a chirally symmetric phase above $`\mu _\text{C}`$. The only exception is the work by Le Yaouanc et al. A. Le Yaouanc et al. (1988) which does not involve effective composite fields but is equivalent to the mean field approximation. In this case the chiral phase transition was found to be of second order.

In this paper we present the equation of state of Hamiltonian lattice QCD at finite density in the strong coupling limit using both Kogut–Susskind and Wilson fermions. As in previous studies we begin with an effective theory by using a Hamiltonian describing the ground state of strongly coupled QCD. However, our approach differs from earlier works in that we do not introduce composite fields but explictly construct a solution to the field equations of motion corresponding to the effective Hamiltonian. This solution exactly diagonalizes the Hamiltonian to second order in field operators for all densities and is used to calculate the vacuum energy density from which we obtain the equation of state. We find that up to and beyond the chiral symmetry restoration density the quark Fermi sea can have negative pressure indicating its mechanical instability. Our result is in qualitative agreement with those obtained using continuum effective QCD models M. Buballa (1996); M. Alford, K. Rajagopal and F. Wilczek (1998) and should be verifiable by future lattice simulationns of strongly coupled QCD at finite density.

Our approach admits to a systematic extension to finite temperature and to the description of bound states. In fact, we first introduce temperature and chemical potential simultaneously into our formalism and then take the limit of vanishing temperature to examine the consequences. Description of bound states is accomplished by interpreting our solution within the context of the N–quantum approach (NQA) to quantum field theory O.W. Greenberg (1965, 1994) which we shall discuss in the concluding section. In the same section we propose how the NQA may be combined with the present approach to study the nature of the deconfinement phase transition.

In Section II we introduce our effective Hamiltonian for the ground state of strong coupling QCD using Wilson fermions, and discuss the condition under which it can be extended to finite $`T`$ and $`\mu `$. Our ansatz for the lattice quark field at finite $`T`$ and $`\mu `$ is presented in Section III. The equation of motion at finite $`\mu `$ is then derived in Section IV and used to diagonalize the effective Hamiltonian to second order in field operators and to evaluate the vacuum energy density. In the same section we determine the unknown quantities in our ansatz by deriving coupled equations for the dynamical quark mass and the total chemical potential and solving them self–consistently. Having constructed a solution for the quark field we present in Section V the equation of state of Hamiltonian lattice QCD at finite density in the strong coupling approximation. We summarize our results in Section VI and discuss how our approach may be extended to incorporate temperature and to describe the deconfinement phase transition. A review of the properties of free lattice Wilson fermions using the Hamiltonian formulation is given in the Appendix.

## II The effective Hamiltonian

We begin by introducing Smit’s effective Hamiltonian J. Smit (1980) describing the ground state of strongly coupled QCD. This state is the one in which no links are excited by the electric flux. It is also infinitely degenerate since various color singlet states may be created at each lattice site without increasing the ground state energy. This degeneracy is lifted by the propagation of quarks on the lattice. The simplest type of such a propagation involves a quark exciting a flux link and an anti–quark deexciting the same link and corresponds to the propagation of a meson. Smit obtained an effective Hamiltonian describing this propagation using second order perturbation theory involving only the quark field $`\mathrm{\Psi }`$ with a nearest neighbour interaction. The Hamiltonian is effective because it only acts on the space of states with no excited links. Nevertheless, it serves our purpose since the main quantity of interest in this work is the vacuum energy density which is obtained by diagonalizing Smit’s Hamiltonian.

In the Hamiltonian formulation of lattice field theory J. Kogut and L. Susskind (1975) only the spatial coordinates are discreticized while the temporal coordinate remain continuous. We adopt the notation of Smit J. Smit (1980) where the discrete sums over the spatial and momentum coordinates are given by

$$\underset{\stackrel{}{x}}{}=a^3\underset{\stackrel{}{m}}{}$$
(1)
and

$$\underset{\stackrel{}{p}}{}=\frac{1}{V}\underset{\stackrel{}{n}}{}$$
(2)
with $`\stackrel{}{x}=a\stackrel{}{m}=a(m_1,m_2,m_3)`$ and $`\stackrel{}{p}=\frac{\pi }{La}\stackrel{}{n}=\frac{\pi }{La}(n_1,n_2,n_3)`$. Here $`a`$ is the lattice spacing and $`L`$ defines the number of unit lattice cells with $`m_l,n_l=0,\pm 1,\pm 2,\mathrm{},\pm L`$. With this notation the volume $`V`$ is given by $`V=(2La)^3`$. Henceforth we shall work in lattice units where $`a=1`$ so that $`\pi p_l+\pi `$, and use Wilson fermions.

The effective Hamiltonian derived by Smit J. Smit (1980) using the temporal gauge is

$`H_{\text{eff}}`$ $`=`$ $`M_0{\displaystyle \underset{\stackrel{}{x}}{}}\left(\mathrm{\Psi }_{a\alpha }^{}\right)_\rho (\stackrel{}{x})\left(\gamma _0\right)_{\rho \nu }\left(\mathrm{\Psi }_{a\alpha }\right)_\nu (\stackrel{}{x})`$ (3)
$`{\displaystyle \frac{K}{2N_\text{c}}}{\displaystyle \underset{\stackrel{}{x}}{}}{\displaystyle \underset{l}{}}`$
$`[\left(\mathrm{\Psi }_{a\alpha }^{}\right)_\rho (\stackrel{}{x})\left(\mathrm{\Sigma }_l\right)_{\rho \nu }\left(\mathrm{\Psi }_{b\alpha }\right)_\nu (\stackrel{}{x}+\widehat{a}_l)\left(\mathrm{\Psi }_{b\beta }^{}\right)_\gamma (\stackrel{}{x}+\widehat{a}_l)\left(\mathrm{\Sigma }_l\right)_{\gamma \delta }^{}\left(\mathrm{\Psi }_{a\beta }\right)_\delta (\stackrel{}{x})`$
$`+\left(\mathrm{\Psi }_{b\beta }^{}\right)_\gamma (\stackrel{}{x}+\widehat{a}_l)\left(\mathrm{\Sigma }_l\right)_{\gamma \delta }^{}\left(\mathrm{\Psi }_{a\beta }\right)_\delta (\stackrel{}{x})\left(\mathrm{\Psi }_{a\alpha }^{}\right)_\rho (\stackrel{}{x})\left(\mathrm{\Sigma }_l\right)_{\rho \nu }\left(\mathrm{\Psi }_{b\alpha }\right)_\nu (\stackrel{}{x}+\widehat{a}_l)]`$
where $`\mathrm{\Sigma }_li(\gamma _0\gamma _lir\gamma _0)`$ and $`\widehat{a}_l`$ is a unit vector along the positive $`l`$–axis. We denote color, flavor and Dirac indices by $`(ab)`$, $`(\alpha \beta )`$ and $`(\rho \nu \gamma \delta )`$, respectively. Summation convention for repeated indices is implied. The three parameters in this Hamiltonian are the Wilson parameter $`r`$ which takes on values between 0 and 1, the current quark mass $`M_0`$ and the effective coupling constant $`K=2N_\text{c}/(N_\text{c}^21)\mathrm{\hspace{0.17em}1}/g^2`$ where $`g`$ is the QCD coupling constant. $`N_\text{c}`$ is the number of colors. When $`r=0`$ the quark fields become Kogut–Susskind fermions.

Smit’s Hamiltonian is valid to order $`𝒪(1/g^2)`$ in the strong coupling expansion. The $`𝒪(1/g^2)`$ corrections involve products of quark bilinears which describe meson propagation mentioned above and are known as “meson terms”. For $`N_\text{c}=3`$, contributions from the subsequent order in the $`1/g^2`$ expansion would consist of products of terms which are trilinear in the quark fields called “baryon terms”. These meson and baryon terms appear in the strong coupling expansions of both Euclidean and Hamiltonian lattice QCD and are the motivations for introducing effective composite meson and baryon fields. In this work we do not take the baryon terms into account but our formalism presented here is also applicable if such terms were present in the effective Hamiltonian.

In the absence of the current quark mass and the Wilson parameter ($`M_0=r=0`$), $`H_{\text{eff}}`$ posseses a $`U(4N_\text{f})`$ symmetry with $`N_\text{f}`$ being the number of flavors. This symmetry is spontaneously broken to $`U(2N_\text{f})U(2N_\text{f})`$ accompanied by the appearance of $`8N_\text{f}^2`$ Goldstone bosons J. Smit (1980). A finite current quark mass also breaks the original $`U(4N_\text{f})`$ symmetry, albeit explicitly, to $`U(2N_\text{f})U(2N_\text{f})`$. Introduction of the Wilson term explicitly breaks the latter symmetry further down to $`U(N_\text{f})`$ thereby solving the fermion doubling problem.

We shall work exclusively in momentum space. Our convention for the Fourier transform from configuration to momentum space is $`\mathrm{\Psi }(\stackrel{}{x})=_\stackrel{}{p}\mathrm{\Psi }(\stackrel{}{p})e^{\stackrel{}{p}\stackrel{}{x}}`$, which implies that the volume $`V`$ is given by $`V=_\stackrel{}{x}=\delta _{\stackrel{}{p},\stackrel{}{p}}.`$ Then the charge conjugation symmetric form of Smit’s Hamiltonian in momentum space is given by

$`H_{\text{eff}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\stackrel{}{p}}{}}M_0\left(\gamma _0\right)_{\rho \nu }[\left(\mathrm{\Psi }_{a\alpha }^{}\right)_\rho (\stackrel{}{p}),\left(\mathrm{\Psi }_{a\alpha }\right)_\nu (\stackrel{}{p})]_{}`$ (4)
$`{\displaystyle \frac{K}{8N_\text{c}}}{\displaystyle \underset{\stackrel{}{p}_1,\mathrm{},\stackrel{}{p}_4}{}}{\displaystyle \underset{l}{}}\delta _{\stackrel{}{p}_1+\mathrm{}+\stackrel{}{p}_4,\stackrel{}{0}}\left[e^{i((\stackrel{}{p}_1+\stackrel{}{p}_2)\widehat{n}_l)}+e^{i((\stackrel{}{p}_3+\stackrel{}{p}_4)\widehat{n}_l)}\right]`$
$`\left[\left(\mathrm{\Psi }_{a\alpha }^{}\right)_\rho (\stackrel{}{p}_1)\left(\mathrm{\Sigma }_l\right)_{\rho \nu }\left(\mathrm{\Psi }_{b\alpha }\right)_\nu (\stackrel{}{p}_2)\left(\mathrm{\Psi }_{a\alpha }\right)_\nu (\stackrel{}{p}_1)\left(\mathrm{\Sigma }_l\right)_{\rho \nu }^{}\left(\mathrm{\Psi }_{b\alpha }^{}\right)_\rho (\stackrel{}{p}_2)\right]`$
$`\left[\left(\mathrm{\Psi }_{b\beta }^{}\right)_\gamma (\stackrel{}{p}_3)\left(\mathrm{\Sigma }_l\right)_{\gamma \delta }^{}\left(\mathrm{\Psi }_{a\beta }\right)_\delta (\stackrel{}{p}_4)\left(\mathrm{\Psi }_{b\beta }\right)_\delta (\stackrel{}{p}_3)\left(\mathrm{\Sigma }_l\right)_{\gamma \delta }\left(\mathrm{\Psi }_{a\beta }^{}\right)_\gamma (\stackrel{}{p}_4)\right]`$
This effective Hamiltonian is the starting point of the present investigation. Our method for obtaining the equation of state consists of extending $`H_{\text{eff}}`$ to finite $`\mu `$ and constructing a quark field operator $`\mathrm{\Psi }`$ which diagonalizes the Hamiltonian to second order in field operators for all densities. Once this solution has been found it can be used to evaluate the vacuum energy density from which we obtain the pressure of the many body system.

However before extending $`H_{\text{eff}}`$ to finite $`T`$ and/or $`\mu `$, it is necessary to impose a condition on these external parameters so that all links would remain in their ground states. In the strong coupling limit the amount of energy required to excite one color electric flux link is

$$E=\frac{1}{2N_\text{c}}(N_\text{c}^21)g^2=\frac{1}{K}$$
(5)
Therefore an extension of $`H_{\text{eff}}`$ to finite $`T`$ and/or $`\mu `$ will be valid as long as $`T,\mu <1/K`$ A. Le Yaouanc et al. (1988)222Note that in A. Le Yaouanc et al. (1988) $`E`$ has been approximated by $`EN_\text{c}g^2`$. since the Hamiltonian only acts on the space of states with no excited links. We shall see that this condition is satisfied in the present work.

The effective Hamiltonian is extended to finite $`T`$ and $`\mu `$ in two steps. The first one is to make the following trivial replacement of the current quark mass term in $`H_{\text{eff}}`$

$$M_0\left(\gamma _0\right)_{\rho \nu }M_0\left(\gamma _0\right)_{\rho \nu }\mu _0\delta _{\rho \nu }$$
(6)
where $`\mu _0`$ is the quark chemical potential. Note that $`\mu _0`$ should not be identified with the total chemical potential $`\mu _{\text{tot}}`$ of the interacting many body system. As we shall see below, the $`𝒪(1/g^2)`$ interaction terms in $`H_{\text{eff}}`$ will induce a correction to $`\mu _0`$ which in general is momentum dependent. We shall therefore refer to $`\mu _0`$ as the ”bare” quark chemical potential and treat it as an input parameter. The second step is to introduce an ansatz for the quark field at finite $`T`$ and $`\mu `$.

## III An ansatz for finite $`T`$ and $`\mu `$

We proceed by presenting our ansatz for the $`\mathrm{\Psi }`$ field in $`H_{\text{eff}}`$ at any temperature and chemical potential. The special case of this ansatz for free space was given in Y. Umino (2000). It has the same structure as the free lattice Dirac field and obeys the free lattice Dirac equation with a mass which is interpreted as the dynamical quark mass. This mass is the only unknown quantity in the free space ansatz and is determined by solving a gap equation. It was shown in Y. Umino (2000) that this ansatz exactly diagonalizes $`H_{\text{eff}}`$ to second order in field operators. Properties of free lattice Dirac fields using Wilson fermions are given in the Appendix.

Temporarily dropping color and flavor indices the free space ansatz given in Y. Umino (2000) is

$$\mathrm{\Psi }_\nu ^{\text{Free}}(t,\stackrel{}{p})=b(\stackrel{}{p})\xi _\nu (\stackrel{}{p})e^{i\omega (\stackrel{}{p})t}+d^{}(\stackrel{}{p})\eta _\nu (\stackrel{}{p})e^{+i\omega (\stackrel{}{p})t}$$
(7)
with $`\nu `$ denoting the Dirac index. The annihilation operators for particles $`b`$ and anti–particles $`d`$ annihilate an interacting vacuum state $`|𝒢_0`$, and obey the free fermion anti–commutation relations. The properties of the lattice spinors $`\xi `$ and $`\eta `$ are given in the Appendix. The free lattice Dirac equation fixes the excitation energy $`\omega (\stackrel{}{p})`$ to be

$$\omega (\stackrel{}{p})=\left(\underset{l}{}\mathrm{sin}^2(\stackrel{}{p}\widehat{n}_l)+M^2(\stackrel{}{p})\right)^{1/2}$$
(8)
where $`M(\stackrel{}{p})`$ is the dynamical quark mass.

In order to extend Eq. (7) to finite $`T`$ and $`\mu `$ we observe that the annihilation operators $`b`$ and $`d`$ in $`\mathrm{\Psi }^{\text{Free}}`$ no longer annihilate the interacting vacuum state at finite $`T`$ and $`\mu `$ denoted as $`|𝒢(T,\mu )`$. To construct operators that annihilate $`|𝒢(T,\mu )`$ we apply a generalized thermal Bogoliubov transformation to the $`b`$ and $`d`$ operators following the formalism of thermal field dynamics H. Umezawa (1992)

$`b(\stackrel{}{p})`$ $`=`$ $`\alpha _pB(\stackrel{}{p})\beta _p\stackrel{~}{B}^{}(\stackrel{}{p})`$ (9a)
$`d(\stackrel{}{p})`$ $`=`$ $`\gamma _pD(\stackrel{}{p})\delta _p\stackrel{~}{D}^{}(\stackrel{}{p})`$ (9b)
The thermal field operators $`B`$ and $`\stackrel{~}{B}^{}`$ annihilate a quasi–particle and create a quasi–hole at finite $`T`$ and $`\mu `$, respectively, while $`D`$ and $`\stackrel{~}{D}^{}`$ are the annihilation operator for a quasi–anti–particle and creation opertor for a quasi–anti–hole, respectively.

These thermal annihilation operators annihilate the interacting thermal vacuum state for each $`T`$ and $`\mu `$.

$$B(\stackrel{}{p})|𝒢(T,\mu )=\stackrel{~}{B}(\stackrel{}{p})|𝒢(T,\mu )=D(\stackrel{}{p})|𝒢(T,\mu )=\stackrel{~}{D}(\stackrel{}{p})|𝒢(T,\mu )=0$$
(10)
We note that the thermal doubling of the Hilbert space accompanying the thermal Bogoliubov transformation is implicit in Eq. (10) where a vacuum state which is annihilated by operators $`B`$, $`\stackrel{~}{B}`$, $`D`$ and $`\stackrel{~}{D}`$ is defined. In addition, since we shall be working only in the space of quantum field operators it is not necessary to specify the structure of the thermal vacuum $`|𝒢(T,\mu )`$.

The thermal operators also satisfy the Fermion anti–commutation relations

$$[B^{}(\stackrel{}{p}),B(\stackrel{}{q})]_+=[D^{}(\stackrel{}{p}),D(\stackrel{}{q})]_+=[\stackrel{~}{B}^{}(\stackrel{}{p}),\stackrel{~}{B}(\stackrel{}{q})]_+=[\stackrel{~}{D}^{}(\stackrel{}{p}),\stackrel{~}{D}(\stackrel{}{q})]_+=\delta _{\stackrel{}{p},\stackrel{}{q}}$$
(11)
with vanishing anti–commutators for the remaining combinations. The coefficients of the transformation are

$`\alpha _p`$ $`=`$ $`\sqrt{1n_p^{}}`$ (12a)
$`\beta _p`$ $`=`$ $`\sqrt{n_p^{}}`$ (12b)
$`\gamma _p`$ $`=`$ $`\sqrt{1n_p^+}`$ (12c)
$`\delta _p`$ $`=`$ $`\sqrt{n_p^+}`$ (12d)
where

$$n_p^\pm =\frac{1}{e^{(\omega (\stackrel{}{p})\pm \mu )/T}+1}$$
(13)
are the Fermi distribution functions for particles ($`n_p^{}`$) and anti–particles ($`n_p^+`$). We stress that the chemical potential appearing in the Fermi distribution functions is the total chemical potential of the interacting many body system. The coefficients are chosen so that $`n_p^\pm `$ are given by

$`n_p^{}`$ $`=`$ $`𝒢(T,\mu )|b^{}(\stackrel{}{p})b(\stackrel{}{p})|𝒢(T,\mu )`$ (14a)
$`n_p^+`$ $`=`$ $`𝒢(T,\mu )|d^{}(\stackrel{}{p})d(\stackrel{}{p})|𝒢(T,\mu )`$ (14b)
Hence in this approach temperature and chemical potential are introduced simultaneously through the coefficients of the thermal Bogoliubov transformation and are treated on an equal footing.

After applying the Bogoliubov transformation to Eq. (7) our ansatz at finite $`T`$ and $`\mu `$ becomes

$`\mathrm{\Psi }_\nu (t,\stackrel{}{p})`$ $`=`$ $`\left[\alpha _pB(\stackrel{}{p})\beta _p\stackrel{~}{B}^{}(\stackrel{}{p})\right]\xi _\nu (\stackrel{}{p})e^{i[\omega (\stackrel{}{p})\mu _{\text{tot}}]t}`$ (15)
$`+\left[\gamma _pD^{}(\stackrel{}{p})\delta _p\stackrel{~}{D}(\stackrel{}{p})\right]\eta _\nu (\stackrel{}{p})e^{+i[\omega (\stackrel{}{p})+\mu _{\text{tot}}]t}`$
and satisfies the equation of motion corresponding to the free lattice Dirac Hamiltonian at finite chemical potential given by

$$H^0=\frac{1}{2}\underset{\stackrel{}{p}}{}\left[\underset{l}{}\text{sin}(\stackrel{}{p}\widehat{n}_l)(\gamma _0\gamma _l)_{\eta \nu }+M(\stackrel{}{p})(\gamma _0)_{\eta \nu }\mu _{\text{tot}}\delta _{\eta \nu }\right][\mathrm{\Psi }_\eta ^{}(t,\stackrel{}{p}),\mathrm{\Psi }_\nu (t,\stackrel{}{p})]_{}$$
(16)
The spinors $`\xi `$ and $`\eta `$ in Eq. (15) obey the same properties as in free space and the excitation energy $`\omega (\stackrel{}{p})`$ has the same form as in Eq. (8). The unknown quantities in our ansatz Eq. (15) are the dynamical quark mass $`M(\stackrel{}{p})`$ and the total chemical potential $`\mu _{\text{tot}}`$ which will be determined in the following section.

In this work we shall take the $`T0`$ limit which amounts to setting $`\gamma _p=1`$ and $`\delta _p=0`$ in the Bogoliubov transformation Eq. (9) thereby suppressing the excitation of anti–holes. In this limit $`\beta _p^2`$ becomes the Heaviside function $`\beta _p^2=\theta (\mu _{\text{tot}}\omega (\stackrel{}{p}))`$ defining the Fermi momentum $`\stackrel{}{p}_F`$ through the relation

$$\mu _{\text{tot}}=\left(\underset{l}{}\mathrm{sin}^2(\stackrel{}{p}_F\widehat{n}_l)+M^2(\stackrel{}{p}_F)\right)^{1/2}$$
(17)
Note that we define chemical potential such that $`\mu _{\text{tot}}M(\stackrel{}{p}_F)`$ which differs from the conventional definition of chemical potential used in lattice calculations where $`\mu 0`$.

One of the simplest quantities to calculate using the ansatz of Eq. (15) in the $`T0`$ limit is the quark number density $`n`$ given by

$`n`$ $`=`$ $`{\displaystyle \frac{1}{2VN_\text{f}N_\text{c}}}\overline{\mathrm{\Psi }}\gamma _0\mathrm{\Psi }={\displaystyle \frac{1}{2VN_\text{f}N_\text{c}}}{\displaystyle \frac{1}{2}}{\displaystyle \underset{\stackrel{}{p}}{}}:[\left(\overline{\mathrm{\Psi }}_{a,\alpha }^{}\right)_\rho (\stackrel{}{p}),\left(\mathrm{\Psi }_{a,\alpha }\right)_\nu (\stackrel{}{p})]_{}:(\gamma _0)_{\rho \nu }`$ (18a)
$`=`$ $`{\displaystyle \underset{\stackrel{}{p}}{}}\beta _p^2={\displaystyle \underset{\stackrel{}{p}}{}}\theta (\mu _{\text{tot}}\omega (\stackrel{}{p}))`$ (18b)
where the symbol : : denotes normal ordering with respect to the vacuum at zero temperature $`|𝒢(T=0,\mu )`$. Therefore, above a sufficiently large value of $`\mu _{\text{tot}}`$ the quark number density becomes a constant which with the present normalization will equal unity. This saturation effect is purely a lattice artifact originating from the $`\text{sin}^2(\stackrel{}{p}\widehat{n}_l)`$ term in $`\omega (\stackrel{}{p})`$.

Another quantity that may be readily calculated using the $`T0`$ ansatz is the chiral condensate. It is found to be proprotional to the dynamical quark mass

$$\frac{1}{2VN_\text{f}N_\text{c}}\overline{\mathrm{\Psi }}\mathrm{\Psi }=\underset{\stackrel{}{p}}{}\alpha _p^2\frac{M(\stackrel{}{p})}{\omega (\stackrel{}{p})}$$
(19)
Below we shall derive a gap equation for $`M(\stackrel{}{p})`$ and show that for a given physically reasonable set of parameters there exists a critical chemical potential above which $`M(\stackrel{}{p})=0`$. Thus the chiral condensate may be identified as being the order parameter for the chiral phase transition at finite density.

## IV Applications of the Equation of Motion

### IV.1 The equation of motion

We now calculate the equation of motion corresponding to $`H_{\text{eff}}`$ with our ansatz for finite $`\mu `$ using two light flavors. The result is used to show that our ansatz exactly diagonalizes the Hamiltonian to second order in field operators for all densities and to calculate the vacuum energy density. In addition, by analyzing the Dirac structure of the equation of motion we derive coupled equations for the dynamical quark mass and the total chemical potential. They are solved to lowest order in the $`1/N_\text{c}`$ expansion thereby completing our construction of a solution to the lattice field theory defined by $`H_{\text{eff}}`$.

The equation of motion for $`H_{\text{eff}}`$ is obtained by exploiting the fact that our ansatz also satisfies the equation of motion corresponding to the free lattice Dirac Hamiltonian $`H^0`$ given in Eq. (16). We therefore have the relation

$$:[\left(\mathrm{\Psi }_{a\alpha }\right)_\nu (t,\stackrel{}{q}),H_{\text{eff}}]_{}:=:[\left(\mathrm{\Psi }_{a\alpha }\right)_\nu (t,\stackrel{}{q}),H^0]_{}:$$
(20)
which plays a crucial role in our construction of a solution for the quark field $`\mathrm{\Psi }`$. Evaluating both sides of Eq. (20) and equating terms which are linear in the field operators we obtain the equation of motion for $`\mathrm{\Psi }`$

$`\left[{\displaystyle \underset{l}{}}\text{sin}(\stackrel{}{q}\widehat{n}_l)(\gamma _0\gamma _l)_{\rho \delta }+M(\stackrel{}{q})(\gamma _0)_{\rho \delta }\mu _{\text{tot}}\delta _{\rho \delta }\right]\left(\mathrm{\Psi }_{a\alpha }\right)_\delta (t,\stackrel{}{q})=`$ (21)
$`\{M_0\left(\gamma _0\right)_{\rho \delta }\mu _0\delta _{\rho \delta }`$
$`+{\displaystyle \frac{1}{N_\text{c}}}K{\displaystyle \underset{\stackrel{}{p}}{}}{\displaystyle \underset{l}{}}\alpha _p^2\mathrm{\Lambda }_{\nu \gamma }^+(\stackrel{}{p})`$
$`[\mathrm{cos}(\stackrel{}{p}\stackrel{}{q})\widehat{n}_l(\left(\mathrm{\Sigma }_l\right)_{\gamma \nu }\left(\mathrm{\Sigma }_l\right)_{\rho \delta }^{}+\left(\mathrm{\Sigma }_l\right)_{\rho \nu }^{}\left(\mathrm{\Sigma }_l\right)_{\gamma \delta })`$
$`+\mathrm{cos}(\stackrel{}{p}+\stackrel{}{q})\widehat{n}_l(\left(\mathrm{\Sigma }_l\right)_{\gamma \nu }^{}\left(\mathrm{\Sigma }_l\right)_{\rho \delta }^{}+\left(\mathrm{\Sigma }_l\right)_{\rho \nu }\left(\mathrm{\Sigma }_l\right)_{\gamma \delta })]`$
$`{\displaystyle \frac{1}{N_\text{c}}}{\displaystyle \frac{K}{4}}{\displaystyle \underset{\stackrel{}{p}}{}}{\displaystyle \underset{l}{}}\left[2\alpha _p^2\mathrm{\Lambda }_{\nu \gamma }^+(\stackrel{}{p})\delta _{\nu \gamma }\right]`$
$`[N_\text{c}(\left(\mathrm{\Sigma }_l\right)_{\rho \nu }\left(\mathrm{\Sigma }_l\right)_{\gamma \delta }^{}+\left(\mathrm{\Sigma }_l\right)_{\rho \nu }^{}\left(\mathrm{\Sigma }_l\right)_{\gamma \delta })`$
$`+\mathrm{cos}(\stackrel{}{p}+\stackrel{}{q})\widehat{n}_l(\left(\mathrm{\Sigma }_l\right)_{\rho \nu }^{}\left(\mathrm{\Sigma }_l\right)_{\gamma \delta }^{}+\left(\mathrm{\Sigma }_l\right)_{\rho \nu }\left(\mathrm{\Sigma }_l\right)_{\gamma \delta })]\left\}\right(\mathrm{\Psi }_{a\alpha })_\delta (t,\stackrel{}{q})`$
with $`\mathrm{\Lambda }^+(\stackrel{}{p})\xi (\stackrel{}{p})\xi ^{}(\stackrel{}{p})`$ being the positive energy projection operator defined in Eq. (48a).

### IV.2 Diagonalization of $`H_{\text{eff}}`$ and the vacuum energy density

We shall now show that our $`T0`$ ansatz exactly diagonalizes the effective Hamiltonian to second order in field operators. The diagonalization procedure involves only algebraic substitutions and does not require any approximations. The quantity of interest here is the off–diagonal Hamiltonian which, to second order in field operators, is found to be

$`H_{\text{off}}|𝒢(0,\mu )`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{q}}{}}\{\alpha _q\xi _\rho ^{}(\stackrel{}{q})[M_0\left(\gamma _0\right)_{\rho \delta }\mu _0\delta _{\rho \delta }]`$ (22)
$`+{\displaystyle \frac{1}{N_\text{c}}}K{\displaystyle \underset{\stackrel{}{p}}{}}{\displaystyle \underset{l}{}}\alpha _p^2\alpha _q\mathrm{\Lambda }_{\nu \rho }^+(\stackrel{}{p})`$
$`\xi _\gamma ^{}(\stackrel{}{q})[\mathrm{cos}(\stackrel{}{p}\stackrel{}{q})\widehat{n}_l(\left(\mathrm{\Sigma }_l\right)_{\rho \nu }\left(\mathrm{\Sigma }_l\right)_{\gamma \delta }^{}+\left(\mathrm{\Sigma }_l\right)_{\rho \nu }^{}\left(\mathrm{\Sigma }_l\right)_{\gamma \delta })`$
$`+\mathrm{cos}(\stackrel{}{p}+\stackrel{}{q})\widehat{n}_l(\left(\mathrm{\Sigma }_l\right)_{\rho \nu }^{}\left(\mathrm{\Sigma }_l\right)_{\gamma \delta }^{}+\left(\mathrm{\Sigma }_l\right)_{\rho \nu }\left(\mathrm{\Sigma }_l\right)_{\gamma \delta })]`$
$`{\displaystyle \frac{1}{N_\text{c}}}{\displaystyle \frac{K}{4}}{\displaystyle \underset{\stackrel{}{p}}{}}{\displaystyle \underset{l}{}}\alpha _q\left[2\alpha _p^2\mathrm{\Lambda }_{\nu \gamma }^+(\stackrel{}{p})\delta _{\nu \gamma }\right]`$
$`\xi _\rho ^{}(\stackrel{}{q})[N_\text{c}(\left(\mathrm{\Sigma }_l\right)_{\rho \nu }\left(\mathrm{\Sigma }_l\right)_{\gamma \delta }^{}+\left(\mathrm{\Sigma }_l\right)_{\rho \nu }^{}\left(\mathrm{\Sigma }_l\right)_{\gamma \delta })`$
$`+\mathrm{cos}(\stackrel{}{p}+\stackrel{}{q})\widehat{n}_l(\left(\mathrm{\Sigma }_l\right)_{\rho \nu }^{}\left(\mathrm{\Sigma }_l\right)_{\gamma \delta }^{}+\left(\mathrm{\Sigma }_l\right)_{\rho \nu }\left(\mathrm{\Sigma }_l\right)_{\gamma \delta })]\}\eta _\delta (\stackrel{}{q})`$
$`B_{\alpha ,a}^{}(\stackrel{}{q})D_{\alpha ,a}^{}(\stackrel{}{q})|𝒢(0,\mu )`$
We see from Eq. (22) that the elementary excitations of the effective Hamiltonian are color singlet (quasi) quark–anti–quark excitations coupled to zero total three momentum. They correspond to the meson propagation on the lattice responsible for lifting the degeneracy of the ground state of strongly coupled QCD.

With the use of the equation of motion for the $`\mathrm{\Psi }`$ field Eq. (21), the equation of motion for the $`\eta `$ spinor Eq. (46) and the orthonormality condition for the $`\xi `$ and $`\eta `$ spinors Eq. (38), we can show that

$`H_{\text{off}}|𝒢(0,\mu )`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{q}}{}}\left\{\alpha _q\xi _\rho ^{}(\stackrel{}{q})\left[{\displaystyle \underset{l}{}}\text{sin}(\stackrel{}{q}\widehat{n}_l)(\gamma _0\gamma _l)_{\rho \delta }M(\stackrel{}{q})(\gamma _0)_{\rho \delta }+\mu _{\text{tot}}\delta _{\rho \delta }\right]\eta _\delta (\stackrel{}{q})\right\}`$ (23)
$`B_{\alpha ,a}^{}(\stackrel{}{q})D_{\alpha ,a}^{}(\stackrel{}{q})|𝒢(0,\mu )`$
$`=`$ $`{\displaystyle \underset{\stackrel{}{q}}{}}\left\{\alpha _q\xi _\rho ^{}(\stackrel{}{q})\left[\omega (\stackrel{}{q})+\mu _{\text{tot}}\right]\eta _\rho (\stackrel{}{q})\right\}B_{\alpha ,a}^{}(\stackrel{}{q})D_{\alpha ,a}^{}(\stackrel{}{q})|𝒢(0,\mu )`$
$`=`$ $`0`$
Note that this result is valid for any dynamical quark mass and total chemical potential. Therefore in the $`T0`$ limit our ansatz shown in Eq. (15) exactly diagonalizes the effective Hamiltonian to second order in field operators for all densities.

Having diagonalized the second order Hamiltonian we can proceed to evaluate the vacuum energy density. Using Eq. (21) once more we find

$`{\displaystyle \frac{1}{V}}𝒢(0,\mu )|H_{\text{eff}}|𝒢(0,\mu )`$ $`=`$ $`N_cN_f{\displaystyle \underset{\stackrel{}{p}}{}}\left[\alpha _p^2M_0\text{Tr}\left(\mathrm{\Lambda }^+(\stackrel{}{p})\gamma _0\right)+2\beta _p^2\mu _0\right]`$ (24)
$`K{\displaystyle \underset{\stackrel{}{p},\stackrel{}{q}}{}}{\displaystyle \underset{l}{}}\alpha _p^2\alpha _q^2\mathrm{\Lambda }_{\nu \rho }^+(\stackrel{}{p})\mathrm{\Lambda }_{\delta \gamma }^+(\stackrel{}{q})`$
$`[\mathrm{cos}(\stackrel{}{p}\stackrel{}{q})\widehat{n}_l(\left(\mathrm{\Sigma }_l\right)_{\rho \nu }\left(\mathrm{\Sigma }_l\right)_{\gamma \delta }^{}+\left(\mathrm{\Sigma }_l\right)_{\rho \nu }^{}\left(\mathrm{\Sigma }_l\right)_{\gamma \delta })`$
$`+\mathrm{cos}(\stackrel{}{p}+\stackrel{}{q})\widehat{n}_l(\left(\mathrm{\Sigma }_l\right)_{\rho \nu }^{}\left(\mathrm{\Sigma }_l\right)_{\gamma \delta }^{}+\left(\mathrm{\Sigma }_l\right)_{\rho \nu }\left(\mathrm{\Sigma }_l\right)_{\gamma \delta })]`$
$`+{\displaystyle \frac{K}{2}}{\displaystyle \underset{\stackrel{}{p},\stackrel{}{q}}{}}{\displaystyle \underset{l}{}}\alpha _p^2\left[\alpha _q^2\mathrm{\Lambda }_{\nu \gamma }^+(\stackrel{}{q})\delta _{\nu \gamma }\right]\mathrm{\Lambda }_{\delta \rho }^+(\stackrel{}{p})`$
$`[N_\text{c}(\left(\mathrm{\Sigma }_l\right)_{\rho \nu }\left(\mathrm{\Sigma }_l\right)_{\gamma \delta }^{}+\left(\mathrm{\Sigma }_l\right)_{\rho \nu }^{}\left(\mathrm{\Sigma }_l\right)_{\gamma \delta })`$
$`+\mathrm{cos}(\stackrel{}{p}+\stackrel{}{q})\widehat{n}_l(\left(\mathrm{\Sigma }_l\right)_{\rho \nu }^{}\left(\mathrm{\Sigma }_l\right)_{\gamma \delta }^{}+\left(\mathrm{\Sigma }_l\right)_{\rho \nu }\left(\mathrm{\Sigma }_l\right)_{\gamma \delta })]`$
$`=`$ $`2N_c{\displaystyle \underset{\stackrel{}{p},\stackrel{}{q}}{}}\{\alpha _p^2[{\displaystyle \frac{3}{2}}K(1+r^2)+\omega (\stackrel{}{p})+{\displaystyle \frac{M(\stackrel{}{p})}{\omega (\stackrel{}{p})}}M_0`$
$`{\displaystyle \frac{1}{N_c}}{\displaystyle \frac{K}{2}}(1r^2)\mathrm{cos}(\stackrel{}{p}+\stackrel{}{q})\widehat{n}_l\mu _{\text{tot}}]+(1+\beta _p^2)\mu _0\}`$
For free space the difference of the vacuum energy densities in the Wigner–Weyl ($`M(\stackrel{}{q})=0`$) and Nambu–Goldstone ($`M(\stackrel{}{q})0`$) phases of the theory is positive

$$\mathrm{\Delta }E=\frac{1}{V}𝒢|H_{\mathrm{eff}}|𝒢|_{M(\stackrel{}{q})=0}\frac{1}{V}𝒢|H_{\mathrm{eff}}|𝒢|_{M(\stackrel{}{q})0}>0$$
(25)
Numerically we find that Eq. (25) also holds for finite chemical potential. Therefore the true ground state of our interacting many body system is in the phase with broken chiral symmetry.

### IV.3 Dynamical quark mass and $`\mu _{\text{tot}}`$

We now derive the equations for the dynamical quark mass and the total chemical potential and solve them to determine our solution Eq. (15) for each density at zero temperature. To accomplish this we explicitly evaluate the right hand side of Eq. (21) to reveal its Dirac structure. The result may be cast in the following compact form

$`\left[{\displaystyle \underset{l}{}}\mathrm{sin}(\stackrel{}{q}\widehat{n}_l)(\gamma _0\gamma _l)_{\nu \delta }+M(\stackrel{}{q})(\gamma _0)_{\nu \delta }\mu _{\mathrm{tot}}\delta _{\nu \delta }\right]\left(\mathrm{\Psi }_{a\alpha }\right)_\delta (t,\stackrel{}{q})=`$ (26)
$`\left[A(\stackrel{}{q})(\gamma _0\gamma _l)_{\nu \delta }+B(\stackrel{}{q})(\gamma _0)_{\nu \delta }+C(\stackrel{}{q})\delta _{\nu \delta }\right]\left(\mathrm{\Psi }_{a\alpha }\right)_\delta (t,\stackrel{}{q})`$
The equations for $`M(\stackrel{}{p})`$ and $`\mu _{\mathrm{tot}}`$ are obtained by equating the coefficents of the $`\gamma _0`$ operator and the Kronecker delta function, respectively.

The gap equation determining $`M(\stackrel{}{p})`$ is given by the coefficient $`B(\stackrel{}{q})`$

$`M(\stackrel{}{q})`$ $`=`$ $`B(\stackrel{}{q})`$ (27)
$`=`$ $`M_0+{\displaystyle \frac{3}{2}}K(1r^2){\displaystyle \underset{\stackrel{}{p}}{}}\left(1\beta _p^2\right){\displaystyle \frac{M(\stackrel{}{p})}{\omega (\stackrel{}{p})}}`$
$`+{\displaystyle \frac{K}{N_\text{c}}}{\displaystyle \underset{\stackrel{}{p},l}{}}\left(1\beta _p^2\right){\displaystyle \frac{M(\stackrel{}{p})}{\omega (\stackrel{}{p})}}`$
$`\left\{8r^2\mathrm{cos}(\stackrel{}{p}\widehat{n}_l)\mathrm{cos}(\stackrel{}{q}\widehat{n}_l){\displaystyle \frac{1}{2}}(1+r^2)\mathrm{cos}(\stackrel{}{p}+\stackrel{}{q})\widehat{n}_l\right\}`$
The structure of this gap equation is very similar to the one in free space ($`\beta _p^2=0`$) found in Y. Umino (2000). The dynamical quark mass is a constant to lowest order in $`N_\text{c}`$ but becomes momentum dependent once $`1/N_\text{c}`$ correction is taken into account.

Similarly, the total chemical potential is given by the coefficient $`C(\stackrel{}{q})`$

$`\mu _{\mathrm{tot}}`$ $`=`$ $`C(\stackrel{}{q})`$ (28)
$`=`$ $`\mu _0+{\displaystyle \frac{1}{4}}{\displaystyle \frac{K}{N_\text{c}}}{\displaystyle \underset{l}{}}{\displaystyle \underset{\stackrel{}{p}}{}}\beta _p^2\left[2N_\text{c}\left(1+r^2\right)2\left(1r^2\right)\mathrm{cos}\left(\stackrel{}{p}+\stackrel{}{q}\right)\right]`$
Thus $`\mu _{\text{tot}}`$ is a sum of the bare chemical potential $`\mu _0`$ and an interaction induced chemical potential which is proportional to the effective coupling constant $`K`$. Furthermore, the latter contribution to $`\mu _{\text{tot}}`$ is momentum dependent and this dependence is a $`1/N_\text{c}`$ correction just as in the case of the gap equation. It should be noted that the above shifting of the bare chemical potential by the interaction is not a new effect. For example, in the well–known and well–studied Nambu–Jona–Lasinio model Y. Nambu and G. Jona–Lasinio (1961) at finite $`T`$ and $`\mu `$ the interaction induces a contribution to the total chemical potential which is proportional to the number density M. Asakawa and K. Yazaki (1989); S.P. Klevansky (1992).

The two equations Eqs. (27) and (28) are coupled and therefore solutions for $`M`$ and $`\mu _{\text{tot}}`$ must be found self–consistently for each value of the input parameter $`\mu _0`$. We shall solve the coupled equations to $`𝒪(N_\text{c}^0)`$ which is equivalent to invoking the mean field approximation. At this order in $`N_\text{c}`$ both the dynamical mass and the total chemical potential are momentum independent. It is also the same order in the $`1/N_\text{c}`$ expansion used to obtain results in all previous studies of strongly coupled lattice QCD. All results are presented using $`M_0`$ = 0 and $`N_\text{c}=3`$.

We first discuss the solutions to the gap equation in free space. In Figure 1 we show dynamical quark masses in free space $`M(\text{Free})`$ as functions of the coupling constant $`K`$ for Wilson parameters $`r`$ = 0.00, 0.50 and 0.75. The dynamical quark masses were obtained in a straightfoward manner by solving the free space gap equation

$$M(\text{Free})=M_0+\frac{3}{2}K(1r^2)\underset{\stackrel{}{p}}{}\frac{M(\text{Free})}{\omega (\stackrel{}{p})}$$
(29)
The figure shows that for each value of $`r`$ there exists a critical coupling constant $`K_\text{C}>0`$ above which the theory is in the broken symmetry phase. This is also true for the $`r=0`$ case corresponding to the use of Kogut–Susskind fermions. In this case the symmetry breaking takes place only for $`K0.732`$.

The dependence of the dynamical mass, and consequently of the chiral condensate through Eq. (19), on the coupling constant is qualitatively different from the results obtained previously using the same effective Hamiltonian J. Smit (1980); A. Le Yaouanc, L. Oliver, O. Pène and J.–C. Raynal (1986). In both J. Smit (1980) and A. Le Yaouanc, L. Oliver, O. Pène and J.–C. Raynal (1986) $`q\overline{q}`$ pair condensation occurs for any value of $`K>0`$. We find that the attraction between a quark and an anti–quark must be sufficiently large enough for a $`q\overline{q}`$ pair to condensate in the vacuum. Thus our approach provides a mechanism for chiral symmetry breaking which other approaches do not. In addition, our results are consistent with the works by Finger and Mandula J. Finger and J. Mandula (1982) and by Amer, Le Yaouanc, Oliver, Pene and Raynal A. Amer, A. Le Yaouanc, L. Oliver, O. Pène and J.–C. Raynal (1983) who have shown that in QCD in the Coulomb gauge $`q\overline{q}`$ condensation takes place only above a critical coupling constant.

Examples of finite $`\mu `$ solutions to the coupled equations Eqs. (27) and (28) to lowest order in $`N_\text{c}`$ are shown in Figures 2 and 3. In Figure 2 we show dynamical masses as functions of the bare chemical potential $`\mu _0`$ to highlight the importance of solving the coupled equations consistently. The figure shows the dynamical quark mass obtained by solving only Eqs. (27) with $`\mu _{\text{tot}}=\mu _0`$ as well the mass obtained by solving the coupled equations consistently. Using $`r=0`$ and $`K=0.9`$ a first order phase transition is observed in both cases, but the values of the critical $`\mu _0`$ are 0.825 when $`\mu _{\text{tot}}=\mu _0`$ and 0.785 when the two equations are solved self–consistently. Therefore the critical chemical potential will be overestimated if interaction induced corrections to the bare chemical potential are ignored.

In Figure 3 we present the dynamical mass as a function of the total chemical potential $`\mu _{\text{tot}}`$ for two values of $`K`$ obtained with $`r=0.25`$. From the figure we see that the phase transition can be either first or second order depending on the value of the coupling constant. When $`K`$ = 0.9 we find a second order phase transition with a critical chemical potential of $`(\mu _{\text{tot}})_\text{C}0.716`$, while if the coupling constant is increased to $`K`$ = 1.0 the phase transition becomes first order with a larger critical chemical potential of $`(\mu _{\text{tot}})_\text{C}0.871`$. This increase in the critical chemical potential with $`K`$ has also been observed in N. Bilić, K. Demeterfi and B. Petersson (1992). Furthermore, we find that when $`K=0.9`$ lattice saturation sets in above $`(\mu _{\text{tot}})_\text{C}`$ at around $`\mu _{\text{tot}}0.898`$ while this effect takes place immediately above $`(\mu _{\text{tot}})_\text{C}`$ for $`K=1.0`$. These values of chemical potentials are smaller than the energy $`E=1/K`$ required to excite one color electric flux link as given in Eq. (5). Therefore with a reasonable set of parameters it is possible to extend Smit’s effective Hamiltonian to finite density as was first pointed out in A. Le Yaouanc et al. (1988).

Having solved the self–consistency equations for the dynamical quark mass and the total chemical potential to lowest order in $`N_\text{c}`$ we have constructed a mean field solution for the quark field appearing in the effective Hamiltonian Eq. (4). In Figure 4 we show the quark number density obtained with this solution as a function of $`\mu _{\text{tot}}`$ for $`K`$ = 0.9 and 1.0. In both cases the number density is a monotonically increasing function of $`\mu _{\mathrm{tot}}`$ in the broken symmetry phase. When $`K=0.9`$ there is a jump in the number density at the phase transition point at $`(\mu _{\text{tot}})_\text{C}0.716`$ from $`n0.037`$ to $`n0.070`$. Beyond this point the number density continues to increase monotonically until when the lattice saturation sets in at $`\mu _{\text{tot}}0.898`$. This behaviour of the number density is qualitatively the same as the one obtained numerically using the Monomer–Dimer–Polymer algorithm as can be seen from a comparison with Figure 5 of F. Karsch and K.–H. Mütter (1989). For $`K=1.0`$ the lattice saturation takes place at the phase transition point at $`(\mu _{\text{tot}})_\text{C}0.871`$ and beyond this point the number density remains a constant at $`n=1`$. Noting that the number density at the phase transition point is $`n0.013`$, the number density for $`K=1.0`$ may be approximated by a Heaviside function of the form $`n=\theta (\mu _{\text{tot}}(\mu _{\text{tot}})_\text{C})`$. This is exactly the result obtained in Eq. (2.47) of E.B. Gregory, S.–H. Guo, H. Kröger and X.–Q. Luo (2000) where a different effective Hamiltonian was used to study strongly coupled lattice QCD at finite density. Furthermore the functional form of the number density found in E.B. Gregory, S.–H. Guo, H. Kröger and X.–Q. Luo (2000) is independent of the strength of the interaction. Therefore the results presented in E.B. Gregory, S.–H. Guo, H. Kröger and X.–Q. Luo (2000) for $`\mu \mu _\text{C}`$ represent nothing but those obtained in the lattice saturation limit.

## V Equation of State

We are now in a position to determine the equation of state by numerically evaluating the thermodynamic potential density using the mean field solution determined above in the vacuum energy density Eq. (24). In Figure 5 we plot pressure as a function of $`\mu _{\text{tot}}`$ for $`K`$ = 0.8 and 0.9. The value of the Wilson parameter is $`r=0.0`$ so that the results have been obtained using Kogut–Susskind fermions. For both values of $`K`$ we find that the pressure of the quark Fermi sea is negative and monotonically decreasing in the broken symmetry phase. For $`K`$ = 0.8 the pressure remains negative but increases in the symmetry restored phase, at least until the lattice saturation point, and has a cusp where the two phases meet. Unfortunately, for $`K=0.9`$ we can not make a definite quantitative statement concerning the behaviour of the pressure in the symmetry restored phase due to lattice saturation, except to mention that there is a discontinuity when going from one phase to another. We find qualitatively similar results when Wilson fermions are used to calculate the pressure as shown in Figure 6. The parameter used in this figure are $`r`$ = 0.25 and $`K`$ = 0.9 and 1.0. We may therefore conclude, at least at the mean field level, that up to and beyond the chiral symmetry restoration point the quark Fermi sea can have negative pressure and therefore can be mechanically unstable with an imaginary speed of sound.

Our conclusion regarding the (strongly coupled) quark matter stability at finite density is consistent with similar studies using effective continuum models of QCD. In the Nambu–Jona–Lasinio model M. Buballa (1996) and the instanton induced ’t Hooft interaction model M. Alford, K. Rajagopal and F. Wilczek (1998), mean field calculations show that cold and dense quark matter may be unstable in the phase with spontaneously broken chiral symmetry, but can become stable in the symmetry restored phase at high enough density. In particular, the result for the pressure obtained in M. Alford, K. Rajagopal and F. Wilczek (1998) is qualitatively the same as the one shown in Figures 5 and 6 as can be seen by comparing the figures with Figure 1 of M. Alford, K. Rajagopal and F. Wilczek (1998). The possibility of unstable quark mattter lead the authors of M. Buballa (1996) and M. Alford, K. Rajagopal and F. Wilczek (1998) to speculate the formation of nucleon droplets, reminiscent of the MIT bag model, in the broken symmetry phase. We shall not indulge on such a speculation here since we are working in an artificial strong coupling regime. Nevertheless, our results concerning the negative pressure is certainly verifiable in future lattice simulations of finite density QCD at strong coupling.

## VI Conclusion and Outlook

In this work we studied the equation of state of two flavored Hamiltonian lattice QCD in the strong coupling limit at finite density using both Kogut–Susskind and Wilson fermions. Starting from an effective lattice Hamiltonian for the ground state of the strongly coupled QCD, we constructed a mean field solution which exactly diagonalizes the Hamiltonian to second order in field operators for all densities. This solution obeys the free lattice Dirac equation with a dynamical quark mass and total chemical potential which are determined by solving a coupled set of equations obtained from the equation of motion. From the gap equation determining the dynamical quark mass we find that at the mean field level the order of the chiral phase transition can be either first or second order depending on the values of input parameters.

The equation of state was obtained by evaluating the thermodynamic potential density from the vacuum energy density using our solution. We find that the pressure of the strongly interacting many body system may be negative in the broken symmetry phase indicating the mechanical instability of our quark Fermi sea. There are indications of this instability beyond the phase transition point although no definite conclusions could be reached for very high densities due to lattice saturation. Nevertheless this behaviour of the pressure was found both for the case of Kogut–Susskind and Wilson fermions and seems, at least at the mean field level, to be robust. In addition, our result concerning negative pressure is in qualitative agreement with studies using continuum effective QCD models, and therefore should certainly be verified by future lattice simulations of strongly coupled QCD at finite density.

To include temperature into our formalism we simply repeat our calculations using the ansatz given in Eq. (15) at non–zero $`T`$. Preliminary calculations indicated that, in addition to particle–anti–particle excitations, the elementary excitations would now involve particle–hole, anti–particle–anti–hole and hole–anti–hole excitations. Because of these additional types of excitations our ansatz would no longer be able to exactly diagonalize the second order Hamiltonian. In fact, a simple exercise would show that at finite $`T`$ even the free lattice Dirac Hamiltonian Eq. (16) is not diagonal due to particle–hole and anti–particle–anti–hole excitations.

We now turn our attention to the possibility of studying the nature of the confinement–deconfinement phase transition. Our solution presented in this work is non–confining and therefore it would be hopeless to use it to study this important phase transition. What is lacking in our formalism is the description of bound states. However, our solution presented here is by no means unique or complete and it can be systematically improved to include all the bound states allowed by the effective Hamiltonian. This is accomplished by interpreting our solution within the context of the N–quantum approach (NQA) to quantum field theory O.W. Greenberg (1965, 1994).

NQA is a method to solve field equations of motion by expanding the interacting Heisenberg fields in terms of asymptotic fields obeying the free field equations of motion. Here the on–shell masses can but need not equal the physical masses of the fields. This expansion is known as the Haag expansion R. Haag (1955) and our ansatz presented here is nothing but the first term in this expansion. Note that because we are working in the Hamiltonian (Kogut–Susskind) formulation of lattice field theory the time variable is continuous and therefore we can introduce and work with the concept of asymptotic fields. The second order terms in the Haag expansion would consist of a product of fermionic quark fields and bosonic elementary color singlet $`\overline{q}q`$ bound state fields. The coefficient of each of the second order terms are interpreted as creation amplitudes for the bound states and are known as Haag amplitudes.

Supressing color and flavor indices for simplicity, our extended ansatz for the Heisenberg quark field $`\mathrm{\Psi }`$ to second order in the Haag expansion in free space will have the following structure

$`\mathrm{\Psi }_\nu ^{\text{Free}}(\stackrel{}{q})`$ $`=`$ $`\mathrm{\Psi }_\nu ^0(\stackrel{}{q})+{\displaystyle \underset{i}{}}{\displaystyle d^3kd^3b\delta ^3(\stackrel{}{q}+\stackrel{}{k}\stackrel{}{b})f_{\nu \rho }^{(i)}(\stackrel{}{k},\stackrel{}{b})}:\mathrm{\Psi }_\rho ^0(\stackrel{}{k})B_{(i)}^0(\stackrel{}{b}):`$ (30)
where the symbol $`::`$ denotes normal ordering. In Eq. (30), $`B`$ is the elementary bosonic field while the superscript $`0`$ indicates that the fields obey their corresponding free field equations of motion. The Haag amplitudes are denoted by $`f_{\nu \rho }^{(i)}`$ with the sum over the index $`i`$ running through all the possible bound states allowed by the Hamiltonian. These states are the color singlet $`q\overline{q}`$ elementary excitations identified in this work. The basic idea of NQA is to use the field equations of motion and derive integral equations for the Haag amplitudes and solve them to obatin a solution to the equation of motion.

In order to solve for the Haag amplitudes it is necessary to calculate the mass and the coupling constant for each of the bound states. This has been accomplished successfully at finite $`T`$ and $`\mu `$ for the two flavored ‘t Hooft interaction model Y. Umino and V. Vento (2000). In addition to bound state masses and coupling constants it is also possible to determine the widths of these states. This quantity is the key to studying the confinement–deconfinement phase transition within our formalism. In the confined phase the bound states will have vanishing widths while in the deconfined phase we expect to see unbound resonant states with finite widths. Hence we propose to use the widths of the $`q\overline{q}`$ states as an order parameter to study the nature of the deconfinement phase transition within strong coupling QCD.

###### Acknowledgements.

Part of this work was completed while I was a Junior Visiting Scientist at ECT in Trento, Italy. I thank the Center for its hospitality and generous support.

*

## Appendix A Properties of free Wilson fermions

In this Appendix we present the properties of free Wilson fermions on the lattice in the Hamiltonian formulation J. Kogut and L. Susskind (1975). The free lattice Dirac Hamiltonian is given by

$`H^0`$ $`=`$ $`{\displaystyle \frac{1}{2i}}{\displaystyle \underset{\stackrel{}{x},l}{}}\left[\mathrm{\Psi }^{}(\stackrel{}{x})\gamma _0\gamma _l\mathrm{\Psi }(\stackrel{}{x}+\widehat{n}_l)\mathrm{\Psi }^{}(\stackrel{}{x}+\widehat{n}_l)\gamma _0\gamma _l\mathrm{\Psi }(\stackrel{}{x})\right]+M{\displaystyle \underset{\stackrel{}{x}}{}}\mathrm{\Psi }^{}(\stackrel{}{x})\gamma _0\mathrm{\Psi }(\stackrel{}{x})`$ (31)
$`{\displaystyle \frac{r}{2}}{\displaystyle \underset{\stackrel{}{x},l}{}}\left[\mathrm{\Psi }^{}(\stackrel{}{x})\gamma _0\mathrm{\Psi }(\stackrel{}{x}+\widehat{n}_l)+\mathrm{\Psi }^{}(\stackrel{}{x}+\widehat{n}_l)\gamma _0\mathrm{\Psi }(\stackrel{}{x})\right]`$
where the third term is the Wilson term. For $`r=0`$ there is an eightfold fermion multiplicity which is removed when $`r0`$. At each lattice site the free Dirac field in configuration space is given by

$$\mathrm{\Psi }_\nu (t,\stackrel{}{x})=\underset{\stackrel{}{p}}{}\left[b(\stackrel{}{p})\xi _\nu (\stackrel{}{p})e^{i(\omega (\stackrel{}{p})t\stackrel{}{p}\stackrel{}{x})}+d^{}(\stackrel{}{p})\eta _\nu (\stackrel{}{p})e^{i(\omega (\stackrel{}{p})t\stackrel{}{p}\stackrel{}{x})}\right]$$
(32)
with $`\nu `$ denoting the Dirac index. The excitation energy $`\omega (\stackrel{}{p})`$ will be determined shortly. The annihilation operators $`b`$ and $`d`$ annihilate the non–interacting vacuum state $`|\mathrm{\hspace{0.17em}0}`$. For our purpose it is not necessary to know the structure of the spinors $`\xi `$ and $`\eta `$.

The only assumption that we shall make is that the creation and annihilation operators obey the free fermion anti–commutation relations

$$[b^{}(\stackrel{}{p}),b(\stackrel{}{q})]_+=[d^{}(\stackrel{}{p}),d(\stackrel{}{q})]_+=\delta _{\stackrel{}{p},\stackrel{}{q}}$$
(33)
Using this assumption we can recover the anti–commutation relations for the field operators

$$[\mathrm{\Psi }_\rho (t,\stackrel{}{x}),\mathrm{\Psi }_\nu ^{}(t,\stackrel{}{y})]_+=\delta _{\stackrel{}{x},\stackrel{}{y}}\delta _{\rho \nu }$$
(34)
provided that $`\xi `$ and $`\eta `$ satisfy the relation

$$\xi _\rho (\stackrel{}{p})\xi _\nu ^{}(\stackrel{}{p})+\eta _\rho (\stackrel{}{p})\eta _\nu ^{}(\stackrel{}{p})=\delta _{\rho \nu }$$
(35)
We normalize the spinors by demanding that the number density is given by

$$𝒩=\underset{\stackrel{}{x}}{}:\mathrm{\Psi }^{}(t,\stackrel{}{x})\mathrm{\Psi }(t,\stackrel{}{x}):=2\underset{\stackrel{}{p}}{}(b^{}(\stackrel{}{p})b(\stackrel{}{p})d^{}(\stackrel{}{p})d(\stackrel{}{p}))$$
(36)
where the symbol : : denotes normal ordering with respect to $`|\mathrm{\hspace{0.17em}0}`$ and the factor of 2 accounts for the spin degrees of freedom. Eq. (36) fixes the normalizations of $`\xi `$ and $`\eta `$ to be

$`\xi _\nu (\stackrel{}{p})\xi _\nu ^{}(\stackrel{}{p})`$ $`=`$ $`\eta _\nu (\stackrel{}{p})\eta _\nu ^{}(\stackrel{}{p})=2`$ (37)
$`\xi _\nu ^{}(\stackrel{}{p})\eta _\nu (\stackrel{}{p})`$ $`=`$ $`\eta _\nu ^{}(\stackrel{}{p})\xi _\nu (\stackrel{}{p})=0`$ (38)
which are consistent with Eq. (35).

In momentum space the charge conjugaton symmetric form of $`H^0`$ is

$$H^0=\frac{1}{2}\underset{\stackrel{}{p}}{}\left(\underset{l}{}\mathrm{sin}(\stackrel{}{p}\widehat{n}_l)\gamma _0\gamma _l+M(\stackrel{}{p})\gamma _0\right)_{\rho \nu }[\mathrm{\Psi }_\rho ^{}(t,\stackrel{}{p}),\mathrm{\Psi }_\nu (t,\stackrel{}{p})]_{}$$
(39)
where the momentum dependent mass term is given by

$$M(\stackrel{}{p})Mr\underset{l}{}\mathrm{cos}(\stackrel{}{p}\widehat{n}_l)$$
(40)
The free Dirac field now becomes

$$\mathrm{\Psi }_\nu (t,\stackrel{}{p})=b(\stackrel{}{p})\xi _\nu (\stackrel{}{p})e^{i\omega (\stackrel{}{p})t}+d^{}(\stackrel{}{p})\eta _\nu (\stackrel{}{p})e^{+i\omega (\stackrel{}{p})t}$$
(41)
which is used to derive the equation of motion corresponding to Eq. (39)

$`i\dot{\mathrm{\Psi }}(t,\stackrel{}{p})`$ $`=`$ $`:[\mathrm{\Psi }(t,\stackrel{}{p}),H^0]_{}:`$ (42)
$`=`$ $`\left({\displaystyle \underset{l}{}}\mathrm{sin}(\stackrel{}{p}\widehat{n}_l)\gamma _0\gamma _l+M(\stackrel{}{p})\gamma _0\right)\mathrm{\Psi }(t,\stackrel{}{p})`$ (43)
From Eq. (43) one obtains the excitation energy

$$\omega (\stackrel{}{p})=\left(\underset{l}{}\mathrm{sin}^2(\stackrel{}{p}\widehat{n}_l)+M^2(\stackrel{}{p})\right)^{1/2}$$
(44)
and the equations of motion for the $`\xi `$ and $`\eta `$ spinors

$`\omega (\stackrel{}{p})\xi (\stackrel{}{p})`$ $`=`$ $`\left({\displaystyle \underset{l}{}}\mathrm{sin}(\stackrel{}{p}\widehat{n}_l)\gamma _0\gamma _l+M(\stackrel{}{p})\gamma _0\right)\xi (\stackrel{}{p})`$ (45)
$`\omega (\stackrel{}{p})\eta (\stackrel{}{p})`$ $`=`$ $`\left({\displaystyle \underset{l}{}}\mathrm{sin}(\stackrel{}{p}\widehat{n}_l)\gamma _0\gamma _l+M(\stackrel{}{p})\gamma _0\right)\eta (\stackrel{}{p})`$ (46)
When $`r=0`$ these equations of motion are relativistic near the eight corners of the Brillouin zone denoted by $`\stackrel{}{\pi }_0=(0,0,0)`$, $`\stackrel{}{\pi }_\text{x}=(\pi ,0,0)`$, $`\stackrel{}{\pi }_\text{y}=(0,\pi ,0)`$, $`\stackrel{}{\pi }_\text{z}=(0,0,\pi )`$, $`\stackrel{}{\pi }_{\text{xy}}=(\pi ,\pi ,0)`$, $`\stackrel{}{\pi }_{\text{xz}}=(\pi ,0,\pi )`$, $`\stackrel{}{\pi }_{\text{yz}}=(0,\pi ,\pi )`$ and $`\stackrel{}{\pi }_{\text{xyz}}=(\pi ,\pi ,\pi )`$. The excitation energies near these values of momenta are equal which corresponds to the eightfold multiplicity mentioned above. This degeneracy is lifted when $`r0`$ due to the momentum dependent mass term Eq. (40). Using the equations of motion for $`\xi `$ and $`\eta `$ it is a simple excercise to show that the off–diagonal Hamiltonian vanishes and that the vacuum energy is given by

$$\mathrm{\hspace{0.17em}0}|H^0|\mathrm{\hspace{0.17em}0}=2V\underset{\stackrel{}{p}}{}\omega (\stackrel{}{p})$$
(47)
Finally, we construct positive and negative energy projection operators $`\mathrm{\Lambda }^+(\stackrel{}{p})`$ and $`\mathrm{\Lambda }^{}(\stackrel{}{p})`$ as follows

$`\mathrm{\Lambda }^+(\stackrel{}{p})`$ $``$ $`\xi (\stackrel{}{p})\xi ^{}(\stackrel{}{p})={\displaystyle \frac{1}{2}}\left[1+{\displaystyle \frac{1}{\omega (\stackrel{}{p})}}{\displaystyle \underset{l}{}}\mathrm{sin}(\stackrel{}{p}\widehat{n}_l)\gamma _0\gamma _l+{\displaystyle \frac{M(\stackrel{}{p})}{\omega (\stackrel{}{p})}}\gamma _0\right]`$ (48a)
$`\mathrm{\Lambda }^{}(\stackrel{}{p})`$ $``$ $`\eta (\stackrel{}{p})\eta ^{}(\stackrel{}{p})={\displaystyle \frac{1}{2}}\left[1{\displaystyle \frac{1}{\omega (\stackrel{}{p})}}{\displaystyle \underset{l}{}}\mathrm{sin}(\stackrel{}{p}\widehat{n}_l)\gamma _0\gamma _l{\displaystyle \frac{M(\stackrel{}{p})}{\omega (\stackrel{}{p})}}\gamma _0\right]`$ (48b)
Note that the projection operators obey the condition

$$\left[\mathrm{\Lambda }^+(\stackrel{}{p})+\mathrm{\Lambda }^{}(\stackrel{}{p})\right]_{\rho \nu }=\delta _{\rho \nu }$$
(49)
as is required by Eq. (35).