text:
# Local Magnetohydrodynamical Models of Layered Accretion Disks

## 1 Introduction

One of the most promising mechanisms for angular momentum transport in accretion disks is magnetohydrodynamical (MHD) turbulence driven by the magnetorotational instability (MRI) (Balbus & Hawley 1998). For the instability to be present, the ionization fraction $`x=n_e/n_H`$ (where $`n_e`$ and $`n_H`$ are the number density of electrons and neutral hydrogen respectively) must be high enough that the neutral-ion collision frequency is greater than the local epicyclic frequency (Blaes & Balbus 1994). This condition can be met at very small values of $`x`$. For example, at $`r=1`$ AU in a disk with a density of $`n=10^{13}`$ cm-3, linear instability requires $`x>10^{11}`$. The fact that such low values of $`x`$ are sufficient to keep the gas and magnetic field well coupled implies that non-ideal MHD effects are important only in a few circumstances, such as in protostellar disks (see, e.g. Stone et al. 2001), or disks in dwarf novae systems (Gammie & Menou 1998).

The nonlinear evolution of the MRI in non-ideal MHD has been studied in several different regimes. If the ions and electrons are well coupled to the magnetic field, but the neutrals drift relative to the ionized species, the plasma is in the ambipolar diffusion regime. Hawley & Stone (1997) used two-fluid simulations to study the saturation of the MRI in this case. On the other hand, if collisions are so frequent that not only are the ions and neutrals well coupled, but also electrical currents are damped, then the plasma is in the Ohmic dissipation regime (Jin 1996; Sano & Miyama 1999). Saturation of the MRI in this case has been studied by Sano, Inutsuka, & Miyama (1998), and Fleming, Stone, & Hawley (2000, hereafter FSH), these authors find sustained MHD turbulence requires the magnetic Reynolds number $`Re_M=C_s^2/\eta \mathrm{\Omega }10^4`$ (where $`C_s`$ is the sound speed, $`\mathrm{\Omega }`$ the rotatonal frequency, and $`\eta `$ the coefficent of resistivity) for poloidal fields with no net flux. Finally, if the positive and negative charge carriers can drift relative to one another, the plasma is in the Hall regime. The linear properties of the MRI are modified in unexpected ways in Hall MHD (Wardle 1999; Balbus & Turquem 2001); the nonlinear saturation of the MRI in Hall MHD has recently been studied by Sano & Stone (2002a; 2002b).

Previous studies of the saturation of the MRI in non-ideal MHD have all been designed to study the local physics of the instability, and have therefore used a local approximation termed the shearing box (Hawley, Gammie, & Balbus 1995, hereafter HGB) with homogeneous initial conditions. In this case, the computational domain represents a small region deep in the interior of the disk. However, there are a variety of important questions that can only be addressed by including the vertical structure of the disk. For example, for the case of ideal MHD, local simulations of vertically stratified disks have been used to study buoyancy of MRI amplified field (Brandenburg et al. 1995; Stone et al. 1996, hereafter SHGB) and the formation of strongly magnetized coronae above weakly magnetized disks (Miller & Stone 1999).

There are reasons to expect even more interesting vertical structure in weakly ionized disks. If the disk is too cold for thermal ionization of rare earth elements to be important ($`T800`$ K (Umebayashi 1983)), then the only sources of ionization will be nonthermal, for example irradiation of the disk by cosmic rays or high-energy photons. At the interstellar ionization rate due to cosmic rays, $`\zeta 10^{17}\mathrm{s}^1`$, the surface layers will be kept sufficiently ionized that they will be well coupled to the magnetic field (Gammie 1996). (The presence of small dust grains can significantly change this picture, however, and lead to much lower ionization fractions, Wardle & Ng 1999; Sano et al. 2000. This will be discussed further in §2.2.) However, the cosmic ray flux is attenuated towards the disk midplane, meaning the ionization fraction drops significantly, and non-ideal MHD effects become important. This has led Gammie (1996) to propose a “layered” model for weakly ionized disks, in which the surface layers are well coupled to the magnetic field and therefore turbulent (the “active layers”), while the midplane is decoupled and queiscent (the “dead zone”). More recently, Fromang et al. (2002) have explored in more detail the ionization structure of protostellar disks, they find the extent of the dead zone is extremely sensitive to the mass accretion rate, the critical magnetic Reynolds number below which the MRI is suppressed, and the abundance of metals in the disk.

There are a number of interesting implications to Gammie’s model. In particular, since the accreting layers of the disk have a fixed depth in surface density, the mass accretion rate increases with the radius of the disk. Thus, steady-state solutions are not possible for layered disks, leading Gammie to suggest that layered accretion could account for FU Orionis outbursts in T Tauri disks as accumulated mass in the dead zone is flushed out by, e.g. gravitational instabilities. This idea has been explored further by Armitage, Livio, & Pringle (2001) by solving the one-dimensional evolution equation for the disk surface density. In addition, the layered disk model has important implications for dust coagulation and planet formation in protoplanetary disks, as the dead zone may provide a sheltered environment for both processes.

In this paper, we present three-dimensional numerical MHD calculations designed to explore the local dynamics of layered disks; that is vertically stratified accretion disks in which the ionization fraction (and therefore resistivity) varies with vertical height. We find the formation of active layers and a dead zone occurs naturally in such a model, with the critical magnetic Reynolds number found by FSH providing a good prediction of the demarcation between areas of high and low angular momentum transport. We are also able to measure the rate of mixing between the turbulent layers and the midplane, and how much heating there is in the midplane due to diffusion of magnetic energy from the active layers. Perhaps our most important result is that we find significant angular momentum transport in the dead layer by non-axisymmetric density waves driven by turbulent motions in the active layers. As discussed in §5, these results have interesting implications for the structure and evolution of layered disks.

The paper is organized as follows. In §2 we describe our numerical methods, initial conditions, and the profile of the resistivity adopted here. In §3 we discuss the results of simulations with initially vertical fields with zero net flux. In §4 we discuss the results of initially azimuthal fields with zero-net flux. Our conclusions are given in §5.

## 2 Method

### 2.1 Equations and Initial Conditions

We adopt the shearing box approximation (HGB) to evolve a local patch of the disk in a Cartesian frame that is in corotation at a radius $`R_0`$ and angular frequency $`\mathrm{\Omega }`$. In this frame, $`x,y,`$ and $`z`$ correspond to the radial, toroidal, and vertical directions, respectively. In these coordinates the equations of isothermal MHD may be written as

$$\frac{\rho }{t}+(\rho 𝐯)=0$$
(1)
$$\frac{𝐯}{t}+𝐯𝐯=\frac{1}{\rho }\left(P+\frac{B^2}{8\pi }\right)+\frac{𝐁𝐁}{4\pi \rho }2𝛀\times 𝐯+3\mathrm{\Omega }^2x\widehat{𝐱}\mathrm{\Omega }^2z\widehat{𝐳}$$
(2)
$$\frac{𝐁}{t}=\times [(𝐯\times 𝐁)\eta (z)𝐉],$$
(3)
$$P=c_s^2\rho $$
(4)
where $`c_s`$ is the sound speed, $`𝐉=(c/4\pi )\times 𝐁`$ is the current density, and $`\eta (z)`$ is the Ohmic resistivity (which is assumed to be a fixed function of vertical height), and the other symbols have their usual meaning. The method used to compute the vertical profile of the resistivity $`\eta (z)`$ described in the next section. The term $`\mathrm{\Omega }^2z\widehat{𝐳}`$ in the momentum equation (eq. 2) is the vertical component of the gravitational field of the central object. As in SHGB, we use an isothermal equation of state (equation 4) to prevent an increase in the scale height of the disk due to heating associated with turbulence in the nonlinear regime of the MRI.

All of our simulations begin with the disk in a state of hydrostatic equilibrium. Thus, the density profile is Gaussian,

$$\rho =\rho _0e^{z^2/H^2}$$
(5)
where $`H^2=2c_s^2/\mathrm{\Omega }^2`$ is the scale height of the disk. In order to facilitate comparison with SHGB we adopt the same parameters and initial conditions as that work. Thus, we set $`\mathrm{\Omega }=10^3`$, $`H=1`$, and $`\rho _0=1`$, leading to $`c_s=7.071\times 10^4`$ and $`P_0=5.0\times 10^7`$.

Our computational domain spans $`H/2xH/2`$, $`0y2\pi H`$, and $`2Hz2H`$. Our standard numerical resolution is $`60\times 128\times 256`$ grid points (equivalent to the high resolution runs in SHGB), except that in order to better resolve the vertical structure of the disk, we have used twice the number of grid points in $`z`$.

We study two types of initial magnetic field configurations: either initially purely vertical or purely azimuthal fields. In both cases the field strength is proportional to $`\mathrm{sin}2\pi x`$, so there is zero net magnetic flux through our computational domain. The magnetic field amplitude is parametrized by $`\beta `$, the ratio of thermal to magnetic pressure. At the midplane of the disk, $`\beta (0)=P(0)8\pi /B^2`$, all of the simulations presented here begin with $`\beta (0)=400`$. The exponential decrease in pressure and density reduces the plasma parameter $`\beta `$ as well. From a local analysis in ideal MHD, the critical wavelength of the MRI $`\lambda _c=9.18\beta ^{1/2}`$ (HGB), therefore we expect that the most unstable wavelengths will increase with height. For azimuthal fields, the field strength is given a gaussian profile in the vertical direction so that $`\beta `$ is initially uniform throughout the volume.

We employ periodic boundary conditions in the azimuthal ($`y`$) direction, and either periodic (for initially vertical field runs) or reflecting (for initially azimuthal field runs) boundary conditions in the vertical ($`z`$) direction. In the radial ($`x`$) direction, shearing-periodic boundary conditions are implemented on all runs (HGB). Outflow boundary conditions that allow mass, momentum, and energy to flow out of the computational volume would provide a more realistic model for accretion disks; however, utilizing outflow conditions along the $`z`$ direction is problematic. As described in SHGB, when strong ($`\beta <1`$) tangled magnetic fields are advected through free boundaries, field loops are broken. As the ends of the loop relax, a large Lorentz force will act on the fluid near the boundary. This disturbance is capable of affecting the entire disk. Because strong, highly tangled magnetic fields are characteristic of MHD turbulence, we cannot use outflow conditions to study the nonlinear evolution of the MRI unless the boundary is sufficiently removed from the disk midplane. For example, local isothermal simulations performed by Miller & Stone (2000) found that at five scale heights above the disk, a magnetically dominated coronal region which is stable to the MRI is formed, thus outflow conditions can be applied there. Since our simulations do not exceed two scale heights above the disk, we must adopt either periodic or reflective boundary conditions. SHGB found that the choice of vertical boundary conditions did not significantly alter the vertical structure of the disk.

The initial velocity field is a Keplerian shear flow represented by $`v_y=3\mathrm{\Omega }x/2`$. The MRI is seeded with small amplitude, spatially uncorrelated fluctuations in the velocity field distributed randomly throughout the simulation volume.

### 2.2 Vertical Profile of Resistivity

Our simulations are motivated primarily by the expected structure of disks around young stellar objects. We consider disk radii at which thermal ionization processes are ineffective. For protostellar disks, these regions lie at a radius

>

0.1AU

>

absent0.1𝐴𝑈\mathrel{\vbox{\offinterlineskip\hbox{$>$}
\kern 1.29167pt\hbox{$\sim$}}}0.1AU (Gammie 1996). We examine the case where the electron fraction is constant with time and cosmic rays are the only source of ionization, although similar results are expected if ionization occurs via X-rays instead (Igea & Glassgold 1999).

The resistivity is related to the electron fraction $`x`$ by

$$\eta =6.5\times 10^3x^1\mathrm{cm}^2\mathrm{s}^1$$
(6)
(Hayashi 1981), and assuming overall charge neutrality, the electron density is determined by

$$\frac{dn_e}{dt}=\zeta n_H\beta x^2n_H^2=0,$$
(7)
which yields

$$x=\left(\frac{\zeta }{\beta n_H}\right)^{1/2}.$$
(8)
where $`\beta `$ is the recombination coefficient. A standard gas phase recombination rate due to dissociative recombination is $`\beta 8.7\times 10^6T^{1/2}`$ cm3 s-1 (Glassgold, Lucas, & Omont 1986). At the surface of the disk we assume the ionization rate due to cosmic rays is the interstellar value, $`\zeta _010^{17}s^1`$ (Spitzer & Tomasko 1968). (We neglect the effect that magnetized outflows, if present, may have on the flux of cosmic rays that reach the disk surface.) However, the flux of cosmic rays is exponentially attenuated within the disk, with a decay length of $`\mathrm{\Sigma }_{CR}=100\mathrm{g}\mathrm{cm}^2`$. Thus, the ionization rate within the disk, $`\zeta (z)`$, is

$$\zeta (z)\zeta _0e^{\mathrm{\Sigma }(z)/\mathrm{\Sigma }_{CR}},$$
(9)
where the column density $`\mathrm{\Sigma }(z)=_z^{\mathrm{}}\rho 𝑑z^{}`$. Combining equations (6), (8), and (9) and normalizing with respect to disk scale height yields

$$\eta (z)=\eta _0\mathrm{exp}(z^2/2)\mathrm{exp}\left(\frac{\mathrm{\Sigma }_0}{\mathrm{\Sigma }_{CR}}\frac{1}{2\sqrt{pi}}_z^{\mathrm{}}e^{z^2}𝑑z^{}\right)$$
(10)
for an isothermal disk, where $`\mathrm{\Sigma }_0`$ is the total column density through the disk.

Our expression for the resistive profile (equation 10) has two free parameters $`\eta _0`$ (the resistivity at the midplane) and the ratio $`\mathrm{\Sigma }_0/\mathrm{\Sigma }_{CR}`$. We choose the former so that the magnetic Reynolds number near the midplane $`Re_M10^3`$, this is below the critical value found by FSH for sustained turbulence and ensures the MRI will be damped in the midplane. The latter we choose to give a large dynamic range in the resitivity. Setting $`\mathrm{\Sigma }_0/\mathrm{\Sigma }_{CR}30`$ yields an $`Re_M`$ that varies by four orders of magnitude over two scale heights. Note this implies the disk surface density is 3000 gm cm-3, about a factor of two larger than the value for the minimum solar nebula. The range of values chosen for $`Re_M`$ was limited by the Courant condition for numerical stability. Choosing values lower than $`Re_M=100`$ would have placed a substantial computational burden upon our simulation. Thus it is worth noting that our midplane values for $`Re_M`$ are somewhat higher than those expected in protostellar disks at radii near 1AU. For example, with a number density of $`10^{13}cm^3`$, T=100 K, and a cosmic ray ionization rate $`\zeta =10^{17}s^1`$, the fractional ionization is $`x=10^{12}`$. This yields a $`Re_M1`$. However, as our simulations will demonstrate, the values of $`Re_M`$ considered in this work are sufficient to stabilize the disk. Moreover, once the critical magnetic Reynolds number which stabilizes the disk is reached, we do not expect the dynamics (such as the Reynolds stress) in the dead zone to be significantly different at lower $`Re_M`$. However, we might expect that decreasing the value of $`Re_M`$ at the midplane for our simulations would result in an increase in the size of the dead zone and therefore a thinner active layer.

We have neglected effects of dust grains in this analysis. For a given ionization rate, the presence of small dust grains can greatly decrease the ionization fraction since the surface of grains can act as a catalyst for recombination reactions. If dust grains with an interstellar size distribution are mixed throughout the disk, it is unlikely that even the surface layers of the disk will be sufficiently ionized to be unstable to the MRI (Sano et al. 2000). However, dust grains may grow in size through coagulation, thereby decreasing their surface area and reducing their effect on the recombination rate. Alternatively, in the abscence of turbulence which keeps the grains well mixed vertically, they may sediment toward the midplane leaving the upper layers of the disk dust free. By neglecting the effects of grains, we restrict the applicability of our results to protostellar disks which are in the late stages of evolution, where grains have either settled to the midplane or grown too large to affect the recombination rate. The calculation of the self-consistent vertical structure of the disk in the presence of small grains requires following both the non-ideal MHD of the gas, as well as the dynamics of the grains driven by aerodynamic forces (to account for gravitational settling and turbulent mixing); such a calculation is beyond the scope of the current paper.

## 3 Vertical Fields with Zero-Net Flux

Table 1 summarizes the parameters and results of our numerical simulations. We first discuss the results from vertical field simulations, which are labeled by the prefix Z (toroidal field simulations, which are discussed in the next section, are labelled with the prefix Y). At the numerical resolution used here the ratio of the critical wavelength $`\lambda _c`$ of the MRI at the midplane (assuming ideal MHD) to the grid spacing $`\mathrm{\Delta }z`$ ($`\mathrm{\Delta }y`$) for poloidal (toroidal) fields is 30, indicating it should be captured if present. (With a finite resistivity, the critical wavelength increases, thus our estimate is a lower bound on the number of grid points per unstable wavelength.) Columns (4), (5), and (6) give the time- and volume-averaged values of magnetic field energy, Maxwell stress, and Reynolds stress in both the active layers and the dead zone, all normalized with respect to the midplane pressure. These averages are taken from late times, when the active layer has reached saturation. The location of the boundaries between the active layers and dead zone is measured from the calculations themselves; as discussed below the boundaries are generally located where $`Re_M`$ is equal to the critical value found by FSH.

### 3.1 A Model with a small dead zone

We first discuss run Z1, comparing it to the equivalent ideal MHD run Z4 (run Z4 is identical to run IZ2 in SHGB but with a higher resolution). The magnetic Reynolds number in Z1 varies from $`1000`$ at the midplane to $`5.6\times 10^6`$ at $`z/H=2`$. Figure 1 is a “spacetime” plot of horizontally averaged quantities $`F(z,t)`$ versus vertical height and time. The plot is created by horizontally averaging each four-dimensional quantity $`f(x,y,z,t)`$ via

$$F(z,t)=\frac{f(x,y,z,t)𝑑x𝑑y}{𝑑x𝑑y}.$$
(11)
The figure can be compared directly to figure 3 in SHGB.

The linear phase of the MRI occurs over the first three orbits of evolution. Suprisingly, the plot of magnetic energy shows the earliest growth occurs near the midplane, where the resistivity is highest and the growth rates of the MRI should be lowest. In fact, this growth is not associated with the MRI at all, but occurs as a consequence of reconnection in the initial field topology. In the highly resistive midplane, the oppositely directed vertical field components in the initial conditions reconnect, forming radial fields that are subsequently sheared by the differential rotation of the disk. This shear produces an azimuthal magnetic field which grows linearly with time. It is this azimuthal field that produces an increase in the magnetic energy. This process does not occur in ideal MHD runs.

The shear amplified field is soon overwhelmed by the exponentially growing MRI. The initial saturation of the MRI occurs around orbit 3; it is evident as a strong peak in the magnetic energy and other variables at this time. The highest field energy occurs at $`|z/H|<1.5`$ as a consequence of disk stratification. The decrease in density with increasing $`|z|`$ means that the Alfvén speed, and also critical wavelength $`\lambda _c`$, rise with increasing $`|z|`$. The increase in $`\lambda _c`$ lowers the peak magnetic field strength at saturation. In fact, some of the long wavelength modes near the top and bottom of the disk may exceed the length of our simulation volume and thus will not be present. Near $`z=0`$ (minimum $`Re_M`$), the field energy is reduced, indicating the shorter wavelength modes of the MRI have been damped by resistivity; only long wavelength modes develop in the midplane region. Note the maximum magnetic energy at orbit 3 has discrete peaks at various vertical positions. These are associated with the channel solution that occurs for vertical fields (Hawley & Balbus 1991; SHGB).

Beyond the linear phase (after orbit 3) we see the breakdown of the channel solution via parasitic instabilities (Goodman & Xu 1994) into MHD turbulence. In an ideal MHD run (cf. figure 3 of SHGB), this turbulence fills the entire volume of the disk. However, the variation of the resistivity with height causes a non-uniform structure in the disk: MHD turbulence is confined to the outer layers alone. Thus, as anticipated by Gammie (1996), the disk may be divided into two regions: two active layers near the surfaces of the disk (characterized by low resistivity and MHD turbulence), and a dead zone near the midplane (characterized by high resistivity and a quiescent flow). The dead zone forms spontaneasouly after saturation, by orbit 15 it is well developed. The boundary between zones is located at approximately $`|z|=z_B=0.4`$. At $`|z|=z_B`$ the magnetic Reynolds number $`Re_M2\times 10^4`$. This is in good agreement with the critical $`Re_M`$ obtained by FSH for zero mean poloidal fields in uniform disks. This dead zone persists throughout the remainder of the simulation. Given the location of the boundaries, and since the vertical profile of density is little changed from the initial state (eq. 5) (see below), the column density in the active layers $`\mathrm{\Sigma }_a`$ compared to the dead zone $`\mathrm{\Sigma }_d`$ follows directly from

$$\frac{\mathrm{\Sigma }_a}{\mathrm{\Sigma }_d}=\frac{2_{z_B}^{\mathrm{}}\mathrm{exp}(z/H)^2}{_{z_B}^{z_B}\mathrm{exp}(z/H)^2}0.75$$
(12)
The two layer structure in the disk is evident in the plots of the Maxwell and Reynolds stress. Angular momentum transport via magnetic fields is confined to the active zone. Interestingly, although the Reynolds stress is reduced in the dead zone, it is not eliminated there. Figure 1 shows that there are irregular increases in Reynolds stress within the dead zone. Vertical striations in the plot of the Reynolds stress are associated with spiral density waves in the disk (SHGB). A comparison with the field energy reveals that the fluctuations in the Reynolds stress correspond with significant field energy increases in the active layers. Strong Lorentz forces in the active zone may be driving these fluctations well into the midplane region. We will examine this question in greater detail below.

The time-averaged vertical structure of the disk is shown in Figure 2, which plots the density, plasma $`\beta `$ parameter, magnetic and kinetic energy, and Maxwell and Reynolds stress averaged over orbits 15-45. The density retains its initial gaussian profile with the same scale height throughout the nonlinear evolution of the MRI, indicating thermal pressure is still the dominant force which balances gravity. Within the dead zone the effects of magnetic diffusion are clearly evident. At the midplane MHD turbulence has been suppressed resulting in an extremely weak ($`\beta =1000`$) field there. This is a significant decline in the field energy as compared to the outer regions; $`\beta `$ in the midplane is $``$ 100 times greater than in the active layers of the disk. The maximum value of $`\beta `$ in the dead zone is $``$ 22 times greater than the peak value in the ideal run Z4. However, in the dead zone of run Z1, the field rises by an order of magnitude as the demarcation line, $`z=z_B`$, is approached. By $`z=1.0`$, $`\beta `$ has declined to a value comparable to that at the corresponding point of the ideal run Z4. This steep increase in field energy is consistent with the high gradient in resistivity present in the dead zone.

Tables 1 and 2 list volume averaged values for a number of quantities in run Z1 in both the active layers and the dead zone. From the tables it is evident that the average magnetic energy is three times larger in the active layer than in the dead zone. The magnetic energy in the active layer is confined mostly to the azimuthal component, with $`B_y^27.5B_x^217.5B_z^2`$. Similar results are obtained in the ideal run BZ4 and in SGHB. Azimuthal field also dominates all other components in the dead zone as well. In all cases the background shear favors the growth of the azimuthal component. The amplitude of each of the field components in the active layers of run Z1 is significantly lower (by a factor of at least two) than in the corresponding locations in the ideal MHD runs. This is mostly due to the finite resistivity in the active layers themselves, rather than due to diffusion into the dead zone. This is shown by the magnitude of the Poynting flux

$$F_P=_s\eta (𝐉\times 𝐁)\mathrm{𝐝𝐀}+_sB^2𝐯\mathrm{𝐝𝐀}_s(𝐯𝐁)𝐁\mathrm{𝐝𝐀}$$
(13)
measured at the boundaries between the active layers and dead zone in run Z1. The average magnetic energy flux out of the active regions in one orbit is only $`0.3\%`$ of the total magnetic energy in those regions. We conclude the lower field energy in those regions is a result of Ohmic dissipation within the active zone itself.

From Table 1, in the active regions the time- and volume-averaged Maxwell stress exceeds the Reynolds stress by a factor of $`4`$; consistent with the ideal MHD results (run Z4 and SHGB). However, as with field energy saturation amplitudes, the stress amplitudes in the active regions are reduced by about 50% compared to their values in the same volume of an ideal MHD simulation due to the finite resistivity in these layers. Within the dead zone, the Maxwell stress drops rapidly (see Figure 2). The time- and volume-averaged Maxwell stress in the dead zone is about 1/4 the value in the active layers. Interestingly, the same effect is not observed with the Reynolds stress. The average Reynolds stress in the dead zone is $`60\%`$ of the average Reynolds stress in the active layers, and $`34\%`$ of the average Reynolds stress for the ideal run. Since the drop in Reynolds stress in the dead zone is not as significant as the drop in Maxwell stress, the Reynolds stress is now the dominant mechanism by which transport may occur. Similar dominance of Reynolds stress over the Maxwell stress at low $`Re_M`$ was observed in the simulations performed by FSH.

When scaled by the thermal pressure at the midplane, the total stress in the active layers in run Z1 is $`\alpha =W_{xy}/P_0=0.006`$, whereas in the region within $`|z|<0.4`$, $`\alpha =0.002`$ The amplitude of the total stress at the midplane (which is dominated by the Reynolds stress) is much larger than expected from the uniform resistivity runs of FSH. For example, the minimum Reynolds stress (at the midplane) is about 50 times larger than the Reynolds stress for a uniform resitivity run with the same $`Re_M`$ from FSH. To test whether turbulence in the active layers is driving velocity fluctuations in the dead zone and producing an enhanced Reynolds stress there, we have measured the flux of kinetic energy,

$$F_K=_S(\rho v^2)𝐯\mathrm{𝐝𝐀},$$
(14)
across the surfaces $`|z|=0.5`$. We find a net flux into the dead zone which transports $`5\%`$ of the average total kinetic energy into the active layers in one orbit. This greatly exceeds the magnitude of the Poynting flux into the dead zone.

Overshooting of turbulent eddies results in significant mas mixing between the active and the dead layers of the disk. The average mixing time (defined as the amount of time taken to transport the entire mass of the dead zone across the active zone boundary) was 3.8 orbits for this simulation. We find later that this time increases as the size of the active layer is decreased.

### 3.2 A Model with a Larger Dead Zone

Run Z2 is identical to Z1 except the resistivity at the midplane is increased so that $`Re_M=100`$ at $`z=0`$ and increases to $`5.6\times 10^5`$ at $`|z|=2`$. This simulation ran for 60 orbits. The qualitative evolution of this model is identical to run Z1: a queiscent region is formed by 15 orbits at the midplane of the disk which is surrounded by active, turbulent layers. Once again the boundary between the active and dead layers is the value of $`z`$ where $`Re_M2\times 10^4`$, in this model this occurs within $`|z|<0.9`$. Thus, the column density of the active layers compared to the dead zone is $`\mathrm{\Sigma }_a/\mathrm{\Sigma }_d0.25`$, three times smaller than in run Z1.

The time-averaged vertical structure of the disk in run Z2 is shown in Figure 3. Here all quantities have been averaged from orbit 30 to 60, after saturation has been reached in the active layers. As with the other simulations, the gaussian density profile remains constant. The lower $`Re_M`$ in the dead zone has had a significant effect on the field energy there; the maximum value of $`\beta `$ has grown by a factor of 8 over the peak for Z1. The total range in $`\beta `$ from the dead zone to the active layers is now over 5 orders of magnitude. From the plot of field energy the saturation amplitude in the active layers has been reduced to $`20\%`$ of its value in Z1 due to the higher resistivity. Within the active layers, the azimuthal component of magnetic field energy dominates by an even greater amount than in Z1; ordering is $`B_y^210B_x^227.5B_z^2`$, averaging from orbit 30 to 60. This is consistent with FSH, where it was found that increasing resistivity resulted in greater dominance of the azimuthal field. In figure 3 we see that the field energy drops off by 2 orders of magnitude between the active zone and the center of the dead zone. This is an order of magnitude greater than the drop off observed for Z1. The saturation amplitude of the kinetic energy are lower in all regions of the disk for Z2 as compared to Z1.

An interesting feature of the evolution of run Z2 is the excitation of vertical oscillations in the disk produced by asymmetrical magnetic pressures in the active layers. The two active layers act as independent dynamos, with local amplification of magnetic field balanced by resistive dissipation. There are significant temporal fluctuations in the total magnetic energy in the active layers associated with the turbulence. Since the two active layers are decoupled, these fluctuations are not correlated, leading to small differences in the magnetic pressure gradients at the boundary of the dead zone. These gradients excite vertical oscillations of the disk midplane with an amplitude of $`0.1H`$ and a period of 4 orbits. The oscillations begin around orbit 10 and continue for the remainder of the calculation with no significant change in amplitude. They do not appear to affect the internal dynamics of the active and dead zones, but might have interesting implications for the production of winds and outflows from the disk.

The time- and volume-averaged Maxwell stress in the active layers dominates over the Reynolds stress by $`3.5`$. When scaled by the thermal pressure, this leads to a total stress of $`\alpha =W_{r\varphi }/P_0=0.00126`$ in the active regions, and $`\alpha =0.00026`$ in the dead zone ($`|z|<0.9`$). These represent reductions of $`79\%`$ and $`91\%`$ as compared to the active and dead zones of Z1. The higher average resistivity in these active layers as compared to Z1 reduces the saturated stress amplitudes. For example, the Maxwell stress in the active regions is $`18\%`$ of the value obtained in the active layers of Z1, and from $`z=1`$ to $`z=0`$ the Maxwell stress drops by approximately 4 orders of magnitude. This represents a far greater drop than occurred in the lower resitive run. However, much like Z1, the Reynolds stress does not suffer a similar decline in the midplane. From Figure 3, we see that the average Reynolds stress within the dead zone does not vary substantialy from its value in the active layer. As with Z1, the velocity fluctuations which give rise to the Reynolds stress are supported by transport of kinetic energy into this region from the active layers. The flux of kinetic energy into the dead zone during one orbit is $`12\%`$ of the average kinetic energy within the active layers. This transport ratio is larger than that which was obtained for Z1, due to the depressed saturation amplitudes in Z2. On the other hand, the average Poynting flux across the interface between the active and dead layers ($`|z|=0.9`$) is insignificant, with the net average value being only $`0.15\%`$ of the average total magnetic energy in the active regions.

The very small amplitude of the Reynolds stress might be a result of purely numerical effects rather than being driven by turbulence in the active layers. To test this, we have restarted run Z2 at orbit 45 with the Lorentz forces removed from the equations of motion (i.e. a pure hydrodynamics simulation). This restart ran for 10 orbits. The effect is immediate. Turbulent flow in the active layers, and the Reynolds stress it produces in the dead zone, quickly decays away on a timescale of only a few orbits. This result is in agreement with the Lorentz force off run performed by HGB. By orbit 55 (ten orbits after the restart), the Reynolds stress in the dead zone is 10 times smaller than when Z2 is continued with MHD. Thus we conclude that the Reynolds stress in the dead layer requires MHD turbulence in the active regions.

Another way to quantify the effect of increasing resistivity on the dynamics of an accretion disk is to measure the mixing time across an arbitrary boundary within the disk, defined as the time required for an amount of mass equal to the total mass within a volume to flow into that volume. The mass flux into the volume $`|z|<0.5`$ was measured for runs Z1, Z2, and Z4. (Note that these fluxes are measured at the same location, although the boundary of the dead zone and active layers is different in each case.) Consistent with our measurements of kinetic energy above, we find that increasing the resistivity decreases the mass flux into the volume $`|z|<0.5`$. The average mixing time measured for the ideal run, Z4, was 3.8 orbits. For the resistive run Z1, this time increases to 5.8 orbits, and for the highly resitive run Z2, it is 6.2 orbits. These results suggest that there is a decreasing level of interaction between the dead and active layers as resistivity is increased. Note that even though we have measured the mass flux at a location well within the dead zone of run Z2, the mixing time is still relatively short.

Finally, we have also run a model identical to runs Z1 and Z2, but with a very large resistivity at the midplane so that $`Re_M=10`$ there; this run is labelled Z3 in Table 1. In this case the magnetic Reynolds is so small throughout the entire volume of the disk that MHD turbulence cannot be sustained. There is transitory growth of the MRI in the early stages of evolution, however once the MRI saturates the field dies away and the disk becomes quiescent. For this reason the model was only run to 20 orbits; the time-averaged values of the Maxwell and Reynolds stress during this peroid are given in Table 1. In comparison to runs Z1 and Z2, the saturation amplitude of both quantities at this $`Re_M`$ is very small. Based on the results of FSH, sustained turbulence and transport can only be supported at such low values of $`Re_M`$ if the magnetic field has a net poloidal flux; our results are consistent with this finding.

## 4 Azimuthal Fields with Zero-Net Flux

To examine whether the saturated state of the MRI in a vertically stratified disk depends on the initial field geometry, we have run two simulations with a purely toroidal field, labeled Y1 and Y2 in Table 1. In run Y1, the profile of the resistivity was identical to that used in run Z1, that is $`Re_M=1000`$ at the midplane, to $`Re_M=5.6\times 10^6`$ at $`|z|=2`$. The resistivity was reduced in run Y2, so that $`Re_M=5700`$, at the midplane, to $`Re_M=2.8\times 10^7`$ at $`|z|=2`$. These simulations were performed at our standard resolution of $`60\times 128\times 256`$, and used reflective boundary conditions at the surfaces $`z=2`$ and $`z=2`$.

Figure 4 is a spacetime plot of the horizontally averaged magnetic energy, Maxwell stress, Reynolds stress, and kinetic energy for run Y1. As expected, the linear growth rate of the MRI on purely azimuthal fields is lower than for vertical fields (Balbus & Hawley 1992), thus saturation does not occur until orbit 20. Also note that shear amplification of radial fields generated by reconnection which is visible at early time in run Z1 (see Figure 1) does not occur here.

Beyond orbit 15, a dead zone ranging from $`z=1`$ to $`z=1`$ is clearly evident within the disk. Despite the same profile of resistivity, this is approximately twice as large as the dead zone observed in Z1. The boundary for the dead zone occurs at $`Re_M3\times 10^5`$. This is an order of magnitude larger than the value of $`Re_M`$ that marks the boundary for the dead zone in BZ1. A similar change in critical $`Re_M`$ with initial field geometery was observed in FSH; an order of magnitude increase in the critical magnetic Reynolds number for instability for a uniform azimuthal compared to a uniform poloidal field geometry. The column density in the active layers compared to the dead zone is $`\mathrm{\Sigma }_a/\mathrm{\Sigma }_d0.187`$, five times smaller than in run Z1.

The volume-averaged properties of Y1 are more similar to run Z2 (which has roughly the same ration of column density in the dead zone to active layers) than run Z1, which has the same resistivity profile. This is because azimuthal fields are more strongly affected by magnetic diffusion than poloidal fields. For example, the magnetic energy in the active zone is approximately 10 times greater than in the dead zone. The ordering of field energy components in the active zone of Y1 is $`B_y^210B_x^233B_z^2`$, which is also similar to that in the corresponding region of Z2.

Figure 5 plots the vertical profiles of the density, plasma parameter $`\beta `$, magnetic and kinetic energy densities, and Maxwell and Reynolds stress averaged from orbits 30 to 45 for Y1. Note this is deep in the nonlinear regime, well after saturation of the MRI. We see that the saturation amplitude of field energy in the active zone is $`2\times 10^3`$ (normalized by $`P_0`$) for Y1, which is $`13\%`$ of its value for Z1. Additionaly, the Maxwell stress in the active layers of Y1 is also only $`13\%`$ of its value for Z1, leading to $`\alpha 8\times 10^4`$ in the active layers. In the center of the dead zone ($`|z|<0.5`$), $`\alpha 5.7\times 10^5`$, a value two orders of magnitude lower than found for Z1. Finally, we find that the kinetic energy in the active region is also reduced by about an order of magnitude as compared to Z1, due to the decreased levels of MHD turbulence in the active layers.

The qualitative behavior of the Maxwell and Reynolds stress in these azimuthal runs is the same as in the poloidal runs; while the Maxwell stress suffers a substantial decline in the dead zone, the Reynolds stress is only slightly lower there than in the active region. In the active layers, Maxwell dominates over Reynolds stress by a factor of $``$ 4. This is the same ratio as in the poloidal field runs. However, within the dead zone, the steep decline of the Maxwell stress causes it to fall below the Reynolds stress. The closer we get to the midplane, the greater the Reynolds stress dominates. For example, within the volume enclosed by $`|z|<0.5`$, the average Reynolds stress is 4.7 times larger than the Maxwell stress. For the volume bounded by $`|z|<0.2`$, the Reynolds stress is two orders of magnitude greater than that observed in Z1 for $`|z|<0.2`$.

Run Y2 has a much lower resistivity in the midplane than Y1 discussed above, so that $`Re_M=5700`$ there. The properties of this run are given in Table 1. We find the production of only a tenous dead zone in the early stages of the evolution, which disappears once the active layers are fully saturated. Shortly after the onset of the MRI, we find that within the region $`|z|<0.5`$, average Maxwell stress is $`10\%`$ of its value in the outer regions. However, as the simulation progresses, the difference between the Maxwell stress in the inner and outer regions decreases. By orbit 40, the average Maxwell stress within $`|z|<0.5`$ has climbed to $`60\%`$ of its value in the outer regions. In addition, we find that at later stages in the simulation the Maxwell stress never falls below the Reynolds stress, even at the midplane. Correspondingly we find that angular momentum transport rates are high in the central regions of the disk during the last 20 orbits of the simulation. Clearly, the existence of a dead zone depends critically on whether or not the magnetic Reynolds number in the disk midplane is small.

## 5 Summary

Using numerical MHD simulations, we have studied the evolution of the MRI in stratified accretion disks including resitivity that varies with vertical height. The vertical profile of resistivity is computed assuming that the dominant source of ionization is external irradiation by cosmic rays (or perhaps X-rays), so that the ionization fraction drops dramatically below a penetration depth of $`\mathrm{\Sigma }_{CR}=100`$ g cm-2. This profile is appropriate for cold disks around young stellar objects, or in dwarf nova systems.

The evolution of two different initial field geometries were examined: a zero net mean poloidal field and a zero net mean azimuthal field. Additionaly, several different column densities for the disk were studied (models with larger column densities have a much higher resistivity at the midplane). In each case, we find the development of MHD turbulence is the less-resistive surface layers of the disk with a queiscent zone in the highly resistive midplane, consistent with the model of layered accretion put forward by Gammie (1996). The boundary between the active layers and the dead zone occurs at the location where the magnetic Reynolds number $`Re_M`$ is equal to the critical value found by FSH for MHD turbulence in the nonlinear regime of the MRI. For initially poloidal fields with zero net flux this critical value is of order $`10^4`$, for initially toroidal fields it is ten times larger. This criterion is therefore a useful means by which the scale of the dead zone may be determined in stratified disks. Since the critical value is smaller for fields with net flux, the size of the dead zone may be correspondingly smaller for that field geometry.

The active regions are charcterized by MHD turbulence and outward angular momentum transport. In these regions, as with ideal simulations, Maxwell stress is the significant source of angular momentum transport, dominating Reynolds stress by a factor of $`4`$. However, the saturation levels for field energy and Maxwell stress are less than those obtained in ideal MHD models since Ohmic dissipation in the active zone is never completely negligable.

Our simulations reveal some unexpected properties of the layered accretion disk model. For example, although the MRI is suppressed in the dead zone the Reynolds stress fails to completely vanish there even when $`Re_M`$ is far below the critical value. Measurement of the flux of magnetic and kinetic energy from the active layers into the dead zone supports the idea that velocity fluctuations in the dead zone are driven by turbulence in the active layers. These fluctuations drive non-axisymetric density waves which transport angular momentum (provide a non-zero Reynolds stress) in the dead zone. Thus, there is a minimum level of transport in the dead zone in the presence of active layers.

Our limited numerical resolution allows us to explore only a small range in the size of the dead zone compared to the active layer; we have performed simulations in which the ratio of the column densities of the active to dead layers $`\mathrm{\Sigma }_a/\mathrm{\Sigma }_d`$ ranges from 0.187 to 0.75. However, over this range the Reynolds stress in the midplane does not seem to depend on the size of the dead zone but rather the amplitude of the turbulence in the active layers; we find the Reynolds stress in the midplane is always about 10% of the Maxwell stress in the active layers. For models with very small active layers finite Ohmic dissipation lowers the saturation amplitude of the MRI, in this case the Reynolds stress in the midplane scaled by the gas pressure is $`10^4`$, whereas in the case of $`\mathrm{\Sigma }_a/\mathrm{\Sigma }_d0.75`$ this ratio was $`10^3`$. Angular momentum transport in the midplane via driven spiral density waves could have important implications for the global evolution of layered accretion disks, for example the timescale for mass accumulation in the dead zone (Gammie 1996) may be changed.

We find significant mass mixing can occur between the active layers and the dead zone via overshooting of turbulent eddies. The average mixing time (defined as the time to transport the entire mass of the dead zone across the boundary with the active layers) varied from 3.8 to 6.2 orbits as the size of the active layers decreased. This mixing may strongly affect dust settling times in the dead zone.

We also find that turbulence in the active layers can drive global, vertical oscillations of the disk when the active layers are large. As is known from ideal MHD simulations, magnetic field fluctuations in the tubulece driven by the MRI can be large. Since each surface layer acts independently in a layered disk, magnetic pressure fluctuations in the two layers can be unbalanced, and provide a small amplitude vertical forcing which drives low-amplitude oscillations.

The flux of energy from the active layers to the dead zone is important for determining the temperature structure in the disk. In particular, if the dead zone is strongly heated by an energy flux from the active layers, it may become sufficiently ionized that it will couple to the magnetic field and become active if the temperature rises above $`1000`$ K. We find the Poynting flux from the active layers to the dead zone to be negligable in all the models considered here. For example, in run Z1 the the Poynting flux was only $`0.3\%`$ of the field energy in the active regions.

Finally, it is important to note that we have ignored the effects of dust in all the models computed here. Sano et al. (2000) have shown that even a tiny fraction of small dust grains can significantly alter the ionization fraction in the disk, and therefore the size of the dead zone and active layers (see also Fromang et al. 2002). However, the effect of dust depends strongly on both the size and spatial distribution of grains. In particular, gravitational settling toward the midplane, or coagulation into larger grains, can reduce the effect of dust. Since gravitational settling depends on the degree to which the disk is turbulent, which in turn depends on the ionization fraction and therefore dust distribution in the disk, it is clear that self-consistent MHD models of dusty disks will require following the motion dust grains in concert with the gas dynamics.

We thank Charles Gammie for insightful comments, especially regarding the calculation of the vertical profile of resistivity. This work was supported by the NSF and by the NASA Origins of Solar Systems program.


format: text

text:
# Hamiltonian lattice quantum chromodynamics at finite density with Wilson fermions

## I Introduction

Quantum Chromodynamics (QCD) is the fundamental theory of strong interactions. It is a SU(3) gauge theory of quarks and gluons. Precise determination of the QCD phase diagram on temperature $`T`$ and chemical potential $`\mu `$ plane will provide valuable information for the experimental search for quark-gluon plasma (QGP). The ultimate goal of machines like the Relativistic Heavy Ion Collider (RHIC) at BNL and the Large Hadron Collider (LHC) at CERN is to create the QGP phase, and replay the birth and evolution of the Universe. Such a new state of matter may also exist in the core of neutron stars or quark stars at low temperature $`T`$ and large chemical potential $`\mu `$. Lattice gauge theory (LGT), proposed by Wilson is a first principle non-perturbative method for QCD. Although it is the most reliable technique for investigating phase transitions in QCD, it is not free of problems: complex action at finite chemical potential and species doubling with naive fermions.

In Lagrangian formulation of LGT at finite chemical potential, the success is limited to SU(2) gauge theory, while in the physical SU(3) case, complex action spoils numerical simulations with importance sampling. Even though much effort has recently been made for SU(3) LGT, and some very interesting information on the phase diagram at large $`T`$ and small $`\mu `$ has been obtained, it is still extremely difficult to do simulations at large chemical potential. QCD at large $`\mu `$ is of particular importance for neutron star or quark star physics. Hamiltonian formulation of LGT doesn’t encounter the notorious “complex action problem”. Recently, we proposed a Hamiltonian approach to LGT with naive fermions at finite chemical potential and solve it in the strong coupling regime. We predicted that at zero temperature, there is a first order chiral phase transition at critical chemical potential $`\mu _C=m_{dyn}^{(0)}=M_N^{(0)}/3`$, with $`m_{dyn}^{(0)}`$ and $`M_N^{(0)}`$ being the dynamical mass of quark and nucleon mass at $`\mu =0`$ respectively. (We expect this is also true for Kogut-Susskind fermions.) By solving the gap and Bethe-Salpeter equations, the authors of Ref.  obtained the critical point the same as ours; but they concluded that the chiral transition is of second order, different from ours. Our order of transition is consistent with other lattice simulation results.

Wilson’s approach to lattice fermions has been extensively used in hadron spectrum calculations as well as in QCD at finite temperature. It avoids the species doubling and preserves the flavor symmetry, but it explicitly breaks the chiral symmetry, one of the most important symmetries of the original theory. Non-perturbative fine-tuning of the bare fermion mass has to be done, in order to define the chiral limit.

In this paper, we study Hamiltonian lattice QCD with Wilson fermions at finite chemical potential. We derive the effective Hamiltonian in the strong coupling regime and diagonalize it by Bogoliubov transformation. The vacuum energy, chiral condensate, and masses of pseudo-scalar, vector meson and nucleons are computed. In the non-perturbatively defined chiral limit, we obtain reasonable results for the critical point and some physical quantities in the large $`N_c`$ limit, with $`N_c`$ the number of colors.

To our knowledge, the only existing literature about the same system ($`r0`$ and $`\mu 0`$) is Ref., where the author used a very different approach: the solution to the gap equation. In contrary to the conventional predictions, the author found that even at $`\mu =0`$, there is a critical value for the effective four fermion coupling $`K`$, below which dynamical mass of quark vanishes. He introduced the concept of total chemical potential and found the transition order depends on the input parameters $`K`$ and $`r`$ as well as the momentum. In contrast, we find that in the chiral limit, dynamical mass of quark doesn’t vanish for all values of $`K`$ if $`\mu <\mu _C`$ (the chiral-symmetry broken phase); and at $`\mu =\mu _C`$, our order of chiral phase transition doesn’t depend on the input parameter.

The rest of the paper is organized as follows. In Sec.II, we derive the effective Hamiltonian at finite chemical potential. In Sec.III, we present the results for the vacuum energy, chiral condensate, and hadron masses. In Sec.IV, we estimate the critical chemical potential at zero temperature. The results are summarized in Sec.V.

## II Effective Hamiltonian in the strong coupling regime

### A The $`\mu =0`$ case

We begin with QCD Hamiltonian***We use a representation of $`\gamma `$ matrices described in D. Lurie, Particles and Fields, Interscience Publishers, John Wiley $`\&`$ Sons, Inc. (1968). with Wilson fermions at chemical potential $`\mu =0`$ on 1 dimensional continuum time and 3 dimensional spatial discretized lattice,

$`H`$ $`=`$ $`M{\displaystyle \underset{x}{}}\overline{\psi }(x)\psi (x)+{\displaystyle \frac{1}{2a}}{\displaystyle \underset{x}{}}{\displaystyle \underset{k=\pm 1}{\overset{\pm d}{}}}\overline{\psi }(x)\gamma _kU(x,k)\psi (x+\widehat{k})`$ (1)
$``$ $`{\displaystyle \frac{r}{2a}}{\displaystyle \underset{x}{}}{\displaystyle \underset{k=\pm 1}{\overset{\pm d}{}}}\overline{\psi }(x)U(x,k)\psi (x+\widehat{k})+{\displaystyle \frac{g^2}{2a}}{\displaystyle \underset{x}{}}{\displaystyle \underset{j=1}{\overset{d}{}}}E_j^\alpha (x)E_j^\alpha (x)`$ (2)
$``$ $`{\displaystyle \frac{1}{ag^2}}{\displaystyle \underset{p}{}}\mathrm{Tr}\left(U_p+U_p^+2\right),`$ (3)
where

$`M=m+{\displaystyle \frac{rd}{a}},`$ (4)
$`d=3`$ is the spatial dimension and $`m`$, $`a`$, $`r`$ and $`g`$ are respectively the bare fermion mass, spatial lattice spacing, Wilson parameter, and bare coupling constant. $`U(x,k)`$ is the gauge link variable at site $`x`$ and direction $`\widehat{k}`$. Fermion field $`\psi `$ on each lattice site carries spin (Dirac), color and flavor indexes; here and in the following, whenever there is a summation sign “$`_x`$”, summations over these indexes are implied. The convention $`\gamma _k=\gamma _k`$ is used. $`E_j^\alpha (x)`$ is the color-electric field at site $`x`$ and direction $`j`$, and summation over $`\alpha =1,2,\mathrm{},8`$ is implied. $`U_p`$ is the product of gauge link variables around an elementary spatial plaquette, and it represents the color magnetic interactions. In the continuum limit $`a0`$, Eq.(3) approaches to the continuum QCD Hamiltonian in the temporal gauge $`A_4=0`$.

The effective Hamiltonian, obtained by strong coupling expansion up to the second order, is

$`H_{eff}`$ $`=`$ $`M{\displaystyle \underset{x}{}}\overline{\psi }_f(x)\psi _f(x){\displaystyle \frac{K(r^2+1)d}{a}}{\displaystyle \underset{x}{}}\psi _f^{}(x)\psi _f(x)`$ (5)
$`+`$ $`{\displaystyle \frac{K}{8aN_c}}{\displaystyle \underset{x}{}}{\displaystyle \underset{k=\pm j}{}}((r^2+1)\psi _{f_1}^{}(x)\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\psi _{f_1}(x+\widehat{k})`$ (6)
$`+`$ $`(r^21)\psi _{f_1}^{}(x)\gamma _4\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\gamma _4\psi _{f_1}(x+\widehat{k})`$ (7)
$``$ $`(r^21)\psi _{f_1}^{}(x)\gamma _5\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\gamma _5\psi _{f_1}(x+\widehat{k})`$ (8)
$`+`$ $`(r^2+1)\psi _{f_1}^{}(x)\gamma _4\gamma _5\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\gamma _4\gamma _5\psi _{f_1}(x+\widehat{k})`$ (9)
$`+`$ $`\left(r^2+(12\delta _{|k|,j})\right)\psi _{f_1}^{}(x)\gamma _4\gamma _j\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\gamma _4\gamma _j\psi _{f_1}(x+\widehat{k})`$ (10)
$``$ $`\left(r^2(12\delta _{|k|,j})\right)\psi _{f_1}^{}(x)\gamma _j\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\gamma _j\psi _{f_1}(x+\widehat{k})`$ (11)
$``$ $`\left(r^2+(12\delta _{|k|,j})\right)\psi _{f_1}^{}(x)\gamma _4\sigma _j\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\gamma _4\sigma _j\psi _{f_1}(x+\widehat{k})`$ (12)
$``$ $`(r^2(12\delta _{|k|,j}))\psi _{f_1}^{}(x)\sigma _j\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\sigma _j\psi _{f_1}(x+\widehat{k})).`$ (13)
This effective Hamiltonian is equivalent to that in Ref. , though different representations of $`\gamma `$ matrices are used. Here $`\sigma _j=ϵ_{jj_1j_2}\gamma _{j_1}\gamma _{j_2}`$. The flavor indexes are explicitly written. This effective Hamiltonian describes the nearest-neighbor four fermion interactions, with

$`K={\displaystyle \frac{1}{g^2C_N}}`$ (14)
being the effective four fermion coupling constant. Here $`C_N=(N_c^21)/(2N_c)`$ is the Casimir invariant of the SU$`(N_c)`$ gauge group.

### B The $`\mu 0`$ case

In the continuum, the grand canonical partition function of QCD at finite temperature $`T`$ and chemical potential $`\mu `$ is

$`Z=\mathrm{Tr}\mathrm{e}^{\beta \left(H\mu N\right)},\beta =(k_BT)^1,`$ (15)
where $`k_B`$ is the Boltzmann constant and $`N`$ is particle number operator

$`N={\displaystyle \underset{x}{}}\psi ^{}(x)\psi (x).`$ (16)
According to Eq.(15) and following the procedure in Sec.II A, the role of the Hamiltonian at strong coupling is now played by

$`H_{eff}^\mu =H_{eff}\mu N.`$ (17)
In this Hamiltonian, there are three input parameters: $`r`$, $`m`$ and $`\mu `$. Suppose we study the phase structure of the system in the chiral limit. Such a limit can be reached by fine-tuning the bare quark mass $`m`$ so that the pion becomes massless. In such a case, there are only two free parameters left: $`r`$ and $`\mu `$.

The vacuum energy is the expectation value of $`H\mu N`$ in its ground state $`|\mathrm{\Omega }`$, and also the expectation value of $`H_{eff}^\mu `$ in its ground state $`|\mathrm{\Omega }_{eff}`$, given by

$`E_\mathrm{\Omega }`$ $`=`$ $`\mathrm{\Omega }|H\mu N|\mathrm{\Omega }=\mathrm{\Omega }_{eff}|H_{eff}^\mu |\mathrm{\Omega }_{eff}.`$ (18)
## III Physical quantities at $`\mu 0`$ and $`T=0`$

### A Meson masses

One way to compute the masses of mesons as well as the contributions of the mesons to the vacuum energy, is to bosonize the effective Hamiltonian Eq. (13). We introduce the following operators

$`\mathrm{\Pi }_{f_1f_2}(x)`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{\overline{v}}}}\psi _{f_1}^{}(x)(1\gamma _4)\gamma _5\psi _{f_2}(x),`$ (19)
$`\mathrm{\Pi }_{f2f1}^{}(x)`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{\overline{v}}}}\psi _{f_2}^{}(x)(1+\gamma _4)\gamma _5\psi _{f_1}(x),`$ (20)
$`V_{jf_1f_2}(x)`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{\overline{v}}}}\psi _{f_1}^{}(x)(1\gamma _4)\gamma _j\psi _{f_2}(x),`$ (21)
$`V_{jf_2f_1}^{}(x)`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{\overline{v}}}}\psi _{f_2}^{}(x)(1+\gamma _4)\gamma _j\psi _{f_1}(x).`$ (22)
$`j`$ stands for the positive spatial direction, and $`\overline{v}`$ and $`v^{}`$ denote respectively expectation value of $`\overline{\psi }\psi `$ and $`\psi ^{}\psi `$ in the vacuum state $`|\mathrm{\Omega }_{eff}`$ of $`H_{eff}`$, i.e.,

$`\overline{v}`$ $`=`$ $`\overline{\psi }(x)\psi (x)_{eff}={\displaystyle \frac{1}{N_fN_s}}\mathrm{\Omega }_{eff}|{\displaystyle \underset{x}{}}\overline{\psi }(x)\psi (x)|\mathrm{\Omega }_{eff},`$ (23)
$`v^{}`$ $`=`$ $`\psi ^{}(x)\psi (x)_{eff}={\displaystyle \frac{1}{N_fN_s}}\mathrm{\Omega }_{eff}|{\displaystyle \underset{x}{}}\psi ^{}(x)\psi (x)|\mathrm{\Omega }_{eff}.`$ (24)
Here $`N_s`$ is the total number of lattice sites and $`N_f`$ the number of flavors. It is shown in Appendix A that under the linearization prescription, operators defined in Eq. (22) satisfy the canonical commutation relations for bosons and then the effective Hamiltonian $`H_{eff}^\mu `$ in Eq.(17) can be expressed in terms of these operators in the following way

$`H_{eff}^\mu `$ $``$ $`H_{Linear}^\mu =E_\mathrm{\Omega }^{(0)}+H_\mathrm{\Pi }+H_V,`$ (25)
where

$`E_\mathrm{\Omega }^{(0)}`$ $`=`$ $`N_fN_s[M\overline{v}({\displaystyle \frac{Kd\left(1+r^2\right)}{a}}+\mu )v^{}+{\displaystyle \frac{Kdr^2}{4aN_c}}(v_2^{}+\overline{v}_2)`$ (26)
$``$ $`{\displaystyle \frac{Kd}{4aN_c}}(\overline{v}_2v_2^{}+v_{2\sigma _j}^{}+\overline{v}_{2\sigma _j})],`$ (27)
and

$`H_\mathrm{\Pi }`$ $`=`$ $`\left(2M{\displaystyle \frac{Kd(1r^2)}{aN_c}}\overline{v}\right){\displaystyle \underset{x,f_1,f_2}{}}\mathrm{\Pi }_{f_2f_1}^{}(x)\mathrm{\Pi }_{f_1f_2}(x)`$ (28)
$`+`$ $`{\displaystyle \frac{Kr^2}{4aN_c}}\overline{v}{\displaystyle \underset{x,f_1,f_2,k}{}}\left(\mathrm{\Pi }_{f_1f_2}^{}(x)\mathrm{\Pi }_{f_2f_1}(x+\widehat{k})+\mathrm{\Pi }_{f_1f_2}(x)\mathrm{\Pi }_{f_2f_1}^{}(x+\widehat{k})\right)`$ (29)
$``$ $`{\displaystyle \frac{K}{4aN_c}}\overline{v}{\displaystyle \underset{x,f_1,f_2,k}{}}\left(\mathrm{\Pi }_{f_1f_2}^{}(x)\mathrm{\Pi }_{f_2f_1}^{}(x+\widehat{k})+\mathrm{\Pi }_{f_1f_2}(x)\mathrm{\Pi }_{f_2f_1}(x+\widehat{k})\right),`$ (30)
$`H_V`$ $`=`$ $`\left(2M{\displaystyle \frac{Kd(1r^2)}{aN_c}}\overline{v}\right){\displaystyle \underset{x,f_1,f_2,j}{}}V_{jf_2f_1}^{}(x)V_{jf_1f_2}(x)`$ (31)
$`+`$ $`{\displaystyle \frac{Kr^2}{4aN_c}}\overline{v}{\displaystyle \underset{x,f_1,f_2,k,j}{}}\left(V_{jf_1f_2}^{}(x)V_{jf_2f_1}(x+\widehat{k})+V_{jf_1f_2}(x)V_{jf_2f_1}^{}(x+\widehat{k})\right)`$ (32)
$``$ $`{\displaystyle \frac{K}{4aN_c}}\overline{v}{\displaystyle \underset{x,f_1,f_2,k,j}{}}\left(V_{jf_1f_2}^{}(x)V_{jf_2f_1}^{}(x+\widehat{k})+V_{jf_1f_2}(x)V_{jf_2f_1}(x+\widehat{k})\right)(12\delta _{|k|,j}).`$ (33)
$`v_2^{}`$, $`\overline{v}_2`$, $`v_{2\sigma _j}^{}`$ and $`\overline{v}_{2\sigma _j}`$ are expectation value of four fermion operators

$`\overline{v}_2`$ $`=`$ $`{\displaystyle \frac{1}{2dN_fN_s}}\mathrm{\Omega }_{eff}|{\displaystyle \underset{x,f_1,f_2,k}{}}\overline{\psi }_{f_1}(x)\psi _{f_2}(x)\overline{\psi }_{f_2}(x+\widehat{k})\psi _{f_1}(x+\widehat{k})|\mathrm{\Omega }_{eff},`$ (34)
$`v_2^{}`$ $`=`$ $`{\displaystyle \frac{1}{2dN_fN_s}}\mathrm{\Omega }_{eff}|{\displaystyle \underset{x,f_1,f_2,k}{}}\psi _{f_1}^{}(x)\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\psi _{f_1}(x+\widehat{k})|\mathrm{\Omega }_{eff},`$ (35)
$`\overline{v}_{2\sigma _j}`$ $`=`$ $`{\displaystyle \frac{1}{2dN_fN_s}}\mathrm{\Omega }_{eff}|{\displaystyle \underset{x,f_1,f_2,k,j}{}}\left(r^2+\left(12\delta _{|k|,j}\right)\right)\overline{\psi }_{f_1}(x)\sigma _j\psi _{f_2}(x)\overline{\psi }_{f_2}(x+\widehat{k})\sigma _j\psi _{f_1}(x+\widehat{k})|\mathrm{\Omega }_{eff},`$ (36)
$`v_{2\sigma _j}^{}`$ $`=`$ $`{\displaystyle \frac{1}{2dN_fN_s}}\mathrm{\Omega }_{eff}|{\displaystyle \underset{x,f_1,f_2,k,j}{}}\left(r^2\left(12\delta _{|k|,j}\right)\right)\psi _{f_1}^{}(x)\sigma _j\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\sigma _j\psi _{f_1}(x+\widehat{k})|\mathrm{\Omega }_{eff}.`$ (37)
(38)
After a Fourier transformation

$`\mathrm{\Pi }_{f_1f_2}(x)={\displaystyle \underset{p}{}}e^{ipx}\stackrel{~}{\mathrm{\Pi }}_{f_1f_2}(p),`$ (39)
$`H_\mathrm{\Pi }`$ in Eq.(33) becomes

$`H_\mathrm{\Pi }`$ $`=`$ $`\left(2M{\displaystyle \frac{Kd}{aN_c}}(1r^2)\right){\displaystyle \underset{p,f_1,f_2}{}}\stackrel{~}{\mathrm{\Pi }}_{f_1f_2}^{}(p)\stackrel{~}{\mathrm{\Pi }}_{f_2f_1}(p)`$ (40)
$`+`$ $`{\displaystyle \frac{Kr^2}{2aN_c}}\overline{v}{\displaystyle \underset{f_1f_2}{}}{\displaystyle \underset{p}{}}\left(\stackrel{~}{\mathrm{\Pi }}_{f_1f_2}^{}(p)\stackrel{~}{\mathrm{\Pi }}_{f_2f_1}(p)+\stackrel{~}{\mathrm{\Pi }}(p)_{f_1f_2}\stackrel{~}{\mathrm{\Pi }}_{f_2f_1}^{}(p)\right){\displaystyle \underset{j=1}{\overset{d}{}}}\mathrm{cos}p_ja`$ (41)
$``$ $`{\displaystyle \frac{K}{2aN_c}}\overline{v}{\displaystyle \underset{p,f_1,f_2}{}}\left(\stackrel{~}{\mathrm{\Pi }}_{f_1f_2}^{}(p)\stackrel{~}{\mathrm{\Pi }}_{f_2f_1}^{}(p)+\stackrel{~}{\mathrm{\Pi }}_{f_1f_2}(p)\stackrel{~}{\mathrm{\Pi }}_{f_2f_1}(p)\right){\displaystyle \underset{j=1}{\overset{d}{}}}\mathrm{cos}p_ja.`$ (42)
The Bogoliubov transformation

$`\stackrel{~}{\mathrm{\Pi }}(p)`$ $``$ $`\stackrel{~}{\mathrm{\Pi }}(p)\mathrm{cosh}u_p+\stackrel{~}{\mathrm{\Pi }}^{}(p)\mathrm{sinh}u_p,`$ (43)
$`\stackrel{~}{\mathrm{\Pi }}^{}(p)`$ $``$ $`\stackrel{~}{\mathrm{\Pi }}^{}(p)\mathrm{cosh}u_p+\stackrel{~}{\mathrm{\Pi }}(p)\mathrm{sinh}u_p,`$ (44)
diagonalizes $`H_\mathrm{\Pi }`$ if

$`\mathrm{tanh}2u_p`$ $`=`$ $`{\displaystyle \frac{2G_2}{G_1}}{\displaystyle \underset{l=1}{\overset{d}{}}}\mathrm{cos}p_la,`$ (45)
$`G_1`$ $`=`$ $`2M{\displaystyle \frac{Kd}{aN_c}}\overline{v}(1r^2)+{\displaystyle \frac{Kr^2}{aN_c}}\overline{v}{\displaystyle \underset{l=1}{\overset{d}{}}}\mathrm{cos}p_la,`$ (46)
$`G_2`$ $`=`$ $`{\displaystyle \frac{K}{2aN_c}}\overline{v}.`$ (47)
The resulting $`H_\mathrm{\Pi }`$ is

$`H_\mathrm{\Pi }`$ $`=`$ $`G_1{\displaystyle \underset{p,f_1,f_2}{}}(1\mathrm{tanh}^22u_p)^{\frac{1}{2}}\stackrel{~}{\mathrm{\Pi }}_{f_1f_2}^{}(p)\stackrel{~}{\mathrm{\Pi }}_{f_2f_1}(p)`$ (48)
$``$ $`{\displaystyle \frac{G_1}{2}}N_f^2{\displaystyle \underset{p}{}}\left(\left(1(1\mathrm{tanh}^22u_p)^{\frac{1}{2}}\right)+{\displaystyle \frac{2G_2r^2}{G_1}}{\displaystyle \underset{l=1}{\overset{d}{}}}\mathrm{cos}p_la\right).`$ (49)
Now $`\stackrel{~}{\mathrm{\Pi }}_{f_1f_2}^{}(p)`$ stands for the pseudo-scalar creation operator in the momentum space. According to Eq. (49), the difference between the pseudo-scalar meson energy and vacuum energy is

$`E_\mathrm{\Pi }`$ $`=`$ $`G_1\left(1\mathrm{tanh}^22u_p\right)^{\frac{1}{2}}`$ (50)
$`=`$ $`\left(2M{\displaystyle \frac{Kd}{aN_c}}\overline{v}(1r^2){\displaystyle \frac{K}{aN_c}}(1r^2)\overline{v}{\displaystyle \underset{l=1}{\overset{d}{}}}\mathrm{cos}p_la\right)^{\frac{1}{2}}`$ (51)
$`\times `$ $`\left(2M{\displaystyle \frac{Kd}{aN_c}}\overline{v}(1r^2)+{\displaystyle \frac{K}{aN_c}}(1+r^2)\overline{v}{\displaystyle \underset{l=1}{\overset{d}{}}}\mathrm{cos}p_la\right)^{\frac{1}{2}},`$ (52)
which gives the pseudo-scalar mass when $`p_j=0`$. The pseudo-scalar mass square is

$`M_\mathrm{\Pi }^2=E_\mathrm{\Pi }^2|_{p_j=0}`$ $`=`$ $`4\left(M+{\displaystyle \frac{Kdr^2}{aN_c}}\overline{v}\right)\left(M+{\displaystyle \frac{Kdr^2}{aN_c}}\overline{v}{\displaystyle \frac{Kd}{aN_c}}\overline{v}\right).`$ (53)
In order to define the chiral limit, one has to fine tune $`MM_{chiral}`$ so that the pion becomes massless. From Eq. (53), we get

$`M_{chiral}={\displaystyle \frac{Kdr^2}{aN_c}}\overline{v}.`$ (54)
In this limit, the pseudo-scalar mass square behaves as $`M_\mathrm{\Pi }^2MM_{chiral}`$, which is the PCAC relation.

The $`H_V`$ sector in Eq. (33) can be considered in a similar way. After a Fourier transformation

$`V_j(x)={\displaystyle \underset{p}{}}e^{ipx}\stackrel{~}{V}_j(p)`$ (55)
and a Bogoliubov transformation

$`\stackrel{~}{V}_j(p)`$ $``$ $`\stackrel{~}{V}_j(p)\mathrm{cosh}w_p^{(j)}+\stackrel{~}{V}_j^{}(p)\mathrm{sinh}w_p^{(j)},`$ (56)
$`\stackrel{~}{V}_j^{}(p)`$ $``$ $`\stackrel{~}{V}_j^{}(p)\mathrm{cosh}w_j(p)+\stackrel{~}{V}_j(p)\mathrm{sinh}w_p^{(j)},`$ (57)
$`H_V`$ becomes a diagonalized one

$`H_V`$ $`=`$ $`G_1{\displaystyle \underset{p,j,f_1,f_2}{}}\left(1\mathrm{tanh}^22w_p^{(j)}\right)^{\frac{1}{2}}\stackrel{~}{V}_{jf_1f_2}^{}(p)V_{jf_2f_1}(p)`$ (58)
$``$ $`{\displaystyle \frac{G_1}{2}}N_f^2{\displaystyle \underset{p,j}{}}\left(\left(1\left(1\mathrm{tanh}^22w_p^{(j)}\right)^{\frac{1}{2}}\right)+{\displaystyle \frac{2G_2}{G_1}}r^2{\displaystyle \underset{l=1}{\overset{d}{}}}\mathrm{cos}p_la\right),`$ (59)
if

$`\mathrm{tanh}2w_p^{(j)}={\displaystyle \frac{2G_2}{G_1}}({\displaystyle \underset{l=1}{\overset{d}{}}}\mathrm{cos}p_la2\mathrm{cos}p_ja).`$ (60)
Now $`\stackrel{~}{V}_{jf_1f_2}^{}(p)`$ stands for the vector creation operator in the momentum space. The vector mass is

$`M_V=G_1\left(1\mathrm{tanh}^22w_0^{(j)}\right)^{\frac{1}{2}}^{MM_c}{\displaystyle \frac{2K\sqrt{d1}}{aN_c}}\overline{v}.`$ (61)
According to Eqs. (25), (49) and (59), the vacuum energy reads

$`E_\mathrm{\Omega }`$ $`=`$ $`\mathrm{\Omega }|H^\mu |\mathrm{\Omega }`$ (62)
$`=`$ $`E_\mathrm{\Omega }^{(0)}{\displaystyle \frac{G_1}{2}}N_f^2{\displaystyle \underset{p}{}}\left(\left(1(1\mathrm{tanh}^22u_p)^{\frac{1}{2}}\right)+{\displaystyle \frac{2G_2r^2}{G_1}}{\displaystyle \underset{l=1}{\overset{d}{}}}\mathrm{cos}p_la\right)`$ (63)
$``$ $`{\displaystyle \frac{G_1}{2}}N_f^2{\displaystyle \underset{p,j}{}}\left(\left(1(1\mathrm{tanh}^22w_p^{(j)})^{\frac{1}{2}}\right)+{\displaystyle \frac{2G_2}{G_1}}r^2{\displaystyle \underset{l=1}{\overset{d}{}}}\mathrm{cos}p_la\right).`$ (64)
This also gives the thermodynamic potential (grand potential) at $`T=0`$.

As shown in Ref. , at $`\mu =0`$ the results above are consistent with those of Smit in Ref.  where the $`1/N_c`$ expansion was used.

### B Results in the large $`N_c`$ limit

As shown in Appendix B, in the chiral limit, the dominant contributions to the vacuum energy for large $`N_c`$ is

$`E_\mathrm{\Omega }`$ $``$ $`N_fN_s\left[M_{chiral}\overline{v}\left({\displaystyle \frac{Kd\left(1+r^2\right)}{a}}+\mu \right)v^{}+{\displaystyle \frac{Kdr^2}{4aN_c}}\left(\left(v^{}\right)^2+\overline{v}^2\right){\displaystyle \frac{Kd}{4aN_c}}\left(\overline{v}^2\left(v^{}\right)^2\right)\right].`$ (65)
As shown Appendix C, under the mean-field approximation, i.e., by Wick-contracting a pair of fermion fields in the four fermion terms in Eq. (13), one can obtain a bilinear Hamiltonian in the large $`N_c`$ limit

$`H_{eff}H_{MFA}=A{\displaystyle \underset{x}{}}\overline{\psi }(x)\psi (x)+B{\displaystyle \underset{x}{}}\psi ^{}(x)\psi (x)+C,`$ (67)
where

$`A`$ $`=`$ $`M_{chiral}{\displaystyle \frac{Kd}{2aN_c}}(1r^2)\overline{v},`$ (68)
$`B`$ $`=`$ $`{\displaystyle \frac{Kd(1+r^2)}{a}}\left({\displaystyle \frac{v^{}}{2N_c}}1\right),`$ (69)
$`C`$ $`=`$ $`{\displaystyle \frac{Kd}{4aN_c}}\left((1+r^2)v_{}^{}{}_{}{}^{2}(1r^2)\overline{v}^2\right)N_sN_f.`$ (70)
The coefficient $`A`$ plays the role of dynamical mass of quark. According to Eq. (18) and Eq. (67), the vacuum energy in presence of the chemical potential $`\mu `$ is now $`E_\mathrm{\Omega }=\mathrm{\Omega }_{eff}|H_{MFA}\mu N|\mathrm{\Omega }_{eff}`$, which agrees with Eq. (LABEL:vacuum\_energy\_large\_Nc), derived from the large $`N_c`$ limit of a bosonized Hamiltonian Eq. (25).

In the large $`N_c`$ limit, the fermion field $`\psi `$ can be expressed as

$`\psi (x)=\left(\begin{array}{c}\xi (x)\\ \eta ^{}(x)\end{array}\right).`$ (73)
The 2-spinors $`\xi `$ and $`\eta ^{}`$ are the annihilation operator of positive energy fermion and creation operator of negative energy fermion respectively. Let us define the state $`|n_p,\overline{n}_p`$ in the momentum space by

$`\xi _p|0_p,\overline{n}_p=0,\xi _p^{}|0_p,\overline{n}_p=|1_p,\overline{n}_p,\xi _p|1_p,\overline{n}_p=|0_p,\overline{n}_p,\xi _p^{}|1_p,\overline{n}_p=0,`$ (74)
$`\eta _p|n_p,0_p=0,\eta _p^{}|n_p,0_p=|n_p,1_p,\eta _p|n_p,1_p=|n_p,0_p,\eta _p^{}|n_p,1_p=0.`$ (75)
The numbers $`n_p`$ and $`\overline{n}_p`$ take the values 0 or 1 due to the Pauli principle. By definition, the up and down components of the fermion field are decoupled in such a state $`|n_p,\overline{n}_p`$. For the vacuum state of $`H_{eff}^\mu `$, we make an ansatz

$`|\mathrm{\Omega }_{eff}={\displaystyle \underset{n_p,\overline{n}_p,p}{}}f_{n_p,\overline{n}_p}|n_p,\overline{n}_p.`$ (76)
In the leading $`N_c`$ limit, the dominant contributions to the chiral condensate and quark number density are

$`\overline{\psi }\psi `$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }|\underset{x}{}\overline{\psi }(x)\psi (x)|\mathrm{\Omega }}{N_fN_s}}\overline{v}`$ (77)
$`=`$ $`{\displaystyle \frac{1}{N_fN_s}}{\displaystyle \underset{n_p,\overline{n}_p,p}{}}C_{n_p,\overline{n}_p}n_p,\overline{n}_p|\overline{\psi }\psi |n_p,\overline{n}_p=2N_c(n+\overline{n}1),`$ (78)
$`n_q`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }|\underset{x}{}\psi ^{}(x)\psi (x)|\mathrm{\Omega }}{2N_cN_fN_s}}1{\displaystyle \frac{v^{}}{2N_c}}1`$ (79)
$`=`$ $`{\displaystyle \frac{1}{2N_cN_fN_s}}{\displaystyle \underset{n_p,\overline{n}_p,p}{}}C_{n_p,\overline{n}_p}n_p,\overline{n}_p|\psi ^{}\psi |n_p,\overline{n}_p1=n\overline{n}.`$ (80)
Here we denote $`C_{n_p,\overline{n}_p}=f_{n_p,\overline{n}_p}^2`$. Using Eq. (LABEL:vacuum\_energy\_large\_Nc) and Eq. (80), in the large $`N_c`$ limit we obtain the normalized vacuum energy

$`ϵ_\mathrm{\Omega }`$ $`=`$ $`{\displaystyle \frac{E_\mathrm{\Omega }}{2N_cN_fN_s}}`$ (81)
$`=`$ $`M_{chiral}\left(n+\overline{n}1\right){\displaystyle \frac{Kdr^2}{a}}\left(n\overline{n}+1\right)`$ (82)
$`+`$ $`{\displaystyle \frac{Kdr^2}{a}}\left(n^2+\overline{n}^2+12\overline{n}\right)+{\displaystyle \frac{Kd}{a}}\left(n+\overline{n}2n\overline{n}1\right)\mu \left(n\overline{n}+1\right),`$ (83)
where the quark number $`n`$ and anti-quark number $`\overline{n}`$

$`n`$ $`=`$ $`n={\displaystyle \underset{n_p,\overline{n}_p,p}{}}C_{n_p,\overline{n}_p}n_p,`$ (84)
$`\overline{n}`$ $`=`$ $`\overline{n}={\displaystyle \underset{n_p,\overline{n}_p,p}{}}C_{n_p,\overline{n}_p}\overline{n}_p`$ (85)
are constrained in the range of $`[0,1]`$ and determined by minimizing the vacuum energy.

For a generic nucleon operator $`O_{Nucl}`$ consisting of three quarks, the thermo mass is

$`M_{Nucl}=\mathrm{\Omega }_{eff}|O_{Nucl}H_{eff}^\mu O_{Nucl}^{}|\mathrm{\Omega }_{eff}E_\mathrm{\Omega }.`$ (86)
Under the mean-field approximation in the large $`N_c`$ limit, it becomes (see Appendix D)

$`M_{Nucl}=3(A+B)3\mu .`$ (87)
## IV Phase structure at $`T=0`$ and $`\mu 0`$

We now consider in the larger $`N_c`$ limit, the critical behavior of the system at $`T=0`$ and $`\mu 0`$. The ground state of the system corresponds to the lowest value of the vacuum energy. At some given inputs of Wilson parameter $`r`$ and chemical potential $`\mu `$, we can find the value of $`n`$ and $`\overline{n}`$ when $`ϵ_\mathrm{\Omega }`$ Eq. (83) is minimized. The result is

$`n`$ $`=`$ $`\mathrm{\Theta }(\mu \mu _C),`$ (88)
$`\overline{n}`$ $`=`$ $`0.`$ (89)
Here $`\mathrm{\Theta }(\mu \mu _C)`$ is the step function: it is 0 for $`\mu <\mu _C`$ and 1 for $`\mu >\mu _C`$, where $`\mu _C`$ is the critical chemical potential

$`\mu _C={\displaystyle \frac{Kd}{a}}(1+2r^2).`$ (90)
Substituting Eq. (89) into Eq. (80), we obtain the chiral condensate and quark number density

$`\overline{\psi }\psi `$ $`=`$ $`2N_c\left(\mathrm{\Theta }\left(\mu \mu _C\right)1\right),`$ (91)
$`n_q`$ $`=`$ $`\mathrm{\Theta }(\mu \mu _C).`$ (92)
There is clearly is a first order chiral phase transition. For $`\mu <\mu _C`$, the system is in the confinement phase with chiral-symmetry breaking. For $`\mu >\mu _C`$, chiral symmetry is restored.

According to Eqs. (70), (87) and (89), the thermo mass of the nucleon is

$`M_{Nucl}=M_{Nucl}^{(0)}3\mu ,`$ (93)
where

$`M_{Nucl}^{(0)}=3m_{dyn}^{(0)}`$ (94)
is the nucleon mass at $`\mu =0`$, and

$`m_{dyn}^{(0)}={\displaystyle \frac{Kd(1+r^2)}{a}}`$ (95)
is the dynamical mass of quark at $`\mu =0`$. From Eq. (93), one sees that the nucleon thermo mass vanishes at $`\mu =M_{Nucl}^{(0)}/3=m_{dyn}^{(0)}`$, before the chiral phase transition takes place. This is not surprising because Wilson fermions break explicitly the chiral symmetry. The value of $`\mu _C`$ should coincide with $`M_{Nucl}^{(0)}/3`$ when $`r`$ is very small, i.e., as in the case of naive or Kogut-Susskind fermions .

## V Discussions

In the preceding sections, we have investigated (d+1)-dimensional Hamiltonian lattice QCD at finite density with Wilson fermions in the strong coupling regime. We compute the vacuum energy, meson and nucleon masses, chiral condensate and quark number density. At finite chemical potential, there is an interplay between the bare fermion mass in the chiral limit and the chiral condensate, which has to be determined self-consistently. The critical behavior of the system in the large $`N_c`$ limit is considered: a first order chiral phase transition is found at $`\mu =\mu _C`$; the nucleon thermo mass vanishes before the chiral transition takes place, which is due to the explicit breakdown of chiral symmetry by Wilson fermions.

We have not yet specified the nature of the chiral-symmetric phase for $`\mu >\mu _C`$. Is it a QGP phase or a color-superconducting phase? Up to now, there has been no first principle investigation of such a phase in SU(3) gauge theory. The answer to this question might be very important to our understanding of the formation of the neutron star or quark star.

We also know that the strong coupling regime is far from the continuum limit. One has to develop a new numerical method to study the continuum physics. The Monte Carlo Hamiltonian method developed recently might eventually be useful for such a purpose. We hope to discuss these interesting issues in the future.

###### Acknowledgements.

We thank V. Azcoiti, S.H. Guo, H. Kröger, V. Laliena, and M. Lombardo for useful discussions. X.Q.L. is supported by the National Science Fund for Distinguished Young Scholars (19825117), Key Project of National Science Foundation (10235040), Guangdong Natural Science Foundation (990212), National and Guangdong Ministries of Education, and Hong Kong Foundation of the Zhongshan University Advanced Research Center. Y.Z.F. is partly supported by the National Natural Science Foundation (10235040 and 10275098).

## A BOSONIZATION AND LINEARIZATION

In Sect. III, to bosonize the Hamiltonian Eq. (13), we introduce the operators $`\mathrm{\Pi }`$, $`\mathrm{\Pi }^{}`$, $`V`$ and $`V^{}`$, defined by Eq. (22). In Eq. (13), there are terms having direct correspondence to these operators. In addition, there are also terms irrelevant to these operators.

In order that the operators in Eq. (22) represent the appropriate mesons, they must satisfy the commutation relations for boson operators. However, a direct calculation shows

$`[\mathrm{\Pi }_{f_1,f_2}(x),\mathrm{\Pi }_{f_2,f_1}^{}(x)]`$ $`=`$ $`{\displaystyle \frac{1}{2\overline{v}}}\left(\overline{\psi }_{f_1}(x)\psi _{f_1}(x)+\overline{\psi }_{f_2}(x)\psi _{f_2}(x)\right),`$ (A1)
$`[V_{j,f_1,f_2}(x),V_{l,f_2,f_1}^{}(x)]`$ $`=`$ $`\delta _{jl}{\displaystyle \frac{1}{2\overline{v}}}\left(\overline{\psi }_{f_1}(x)\psi _{f_1}(x)+\overline{\psi }_{f_2}(x)\psi _{f_2}(x)\right)`$ (A3)
$`{\displaystyle \frac{1}{4\overline{v}}}\left(\psi _{f_1}^{}(x)(\gamma _j\gamma _l\gamma _l\gamma _j)\psi _{f_1}(x)+\psi _{f_2}^{}(x)(\gamma _j\gamma _l\gamma _l\gamma _j)\psi _{f_2}(x)\right),`$
which are not consistent with the commutation relations between the annihilation and creation operators for bosons. Similar situation also appears in quantum theory of magnetization and many-particle systems, where a linearization prescription is used to simplify the theory. Using such a procedure, i.e., with the fermion bilinears in the r.h.s. of the commutation relations replaced by the vacuum expectation value $`\overline{\psi }(x)\psi (x)\overline{v}`$ and for degenerate quarks, Eq. (A3) becomes

$`[\mathrm{\Pi }(x),\mathrm{\Pi }^{}(x)]`$ $``$ $`{\displaystyle \frac{1}{\overline{v}}}{\displaystyle \frac{\mathrm{\Omega }_{eff}|\underset{x}{}\overline{\psi }\psi |\mathrm{\Omega }_{eff}}{N_fN_s}}=1,`$ (A4)
$`[V_j(x),V_l^{}(x)]`$ $``$ $`\delta _{jl}{\displaystyle \frac{\mathrm{\Omega }_{eff}|\underset{x}{}\overline{\psi }\psi |\mathrm{\Omega }_{eff}}{N_fN_s\overline{v}}}{\displaystyle \frac{\mathrm{\Omega }_{eff}|\underset{x}{}\psi ^{}(x)(\gamma _j\gamma _l\gamma _l\gamma _j)\psi (x)|\mathrm{\Omega }_{eff}}{2N_fN_s\overline{v}}}=\delta _{jl},`$ (A5)
which are the correct commutation relations for bosons. As shown in Refs. , the linearization prescription leads to results consistent with Ref. , where systematic $`1/N_c`$ expansion was used.

The four fermion terms in Eq. (13), which have direct correspondence to the boson operators, are

$`{\displaystyle \underset{x,k}{}}\left(\psi _{f_1}^{}(x)\gamma _5\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\gamma _5\psi _{f_1}(x+\widehat{k})+\psi _{f_1}^{}(x)\gamma _4\gamma _5\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\gamma _4\gamma _5\psi _{f_1}(x+\widehat{k})\right)`$ (A6)
$`=`$ $`2\overline{v}{\displaystyle \underset{x,k}{}}\left(\mathrm{\Pi }_{f_1f_2}^{}(x)\mathrm{\Pi }_{f_2f_1}^{}(x+\widehat{k})+\mathrm{\Pi }_{f_1f_2}(x)\mathrm{\Pi }_{f_2f_1}(x+\widehat{k})\right),`$ (A8)
$`{\displaystyle \underset{x,k}{}}\left(\psi _{f_1}^{}(x)\gamma _5\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\gamma _5\psi _{f_1}(x+\widehat{k})\psi _{f_1}^{}(x)\gamma _4\gamma _5\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\gamma _4\gamma _5\psi _{f_1}(x+\widehat{k})\right)`$
$`=`$ $`2\overline{v}{\displaystyle \underset{x,k}{}}\left(\mathrm{\Pi }_{f_1f_2}^{}(x)\mathrm{\Pi }_{f_2f_1}(x+\widehat{k})+\mathrm{\Pi }_{f_1f_2}(x)\mathrm{\Pi }_{f_2f_1}^{}(x+\widehat{k})\right),`$ (A11)
$`{\displaystyle \underset{x,k,j}{}}\left(\psi _{f_1}^{}(x)\gamma _4\gamma _j\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\gamma _4\gamma _j\psi _{f_1}(x+\widehat{k})+\psi _{f_1}^{}(x)\gamma _j\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\gamma _j\psi _{f_1}(x+\widehat{k})\right)`$
$`\times \left(12\delta _{|k|,j}\right)`$
$`=`$ $`2\overline{v}{\displaystyle \underset{x,k,j}{}}\left(V_{jf_1f_2}^{}(x)V_{jf_2f_1}^{}(x+\widehat{k})+V_{jf_1f_2}(x)V_{jf_2f_1}(x+\widehat{k})\right)\left(12\delta _{|k|,j}\right),`$ (A13)
$`{\displaystyle \underset{x,k,j}{}}\left(\psi _{f_1}^{}(x)\gamma _4\gamma _j\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\gamma _4\gamma _j\psi _{f_1}(x+\widehat{k})\psi _{f_1}^{}(x)\gamma _j\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\gamma _j\psi _{f_1}(x+\widehat{k})\right)`$
$`=`$ $`2\overline{v}{\displaystyle \underset{x,k,j}{}}\left(V_{jf_1f_2}^{}(x)V_{jf_2f_1}(x+\widehat{k})+V_{jf_1f_2}(x)V_{jf_2f_1}^{}(x+\widehat{k})\right).`$ (A14)
For the bilinear $`_x\overline{\psi }(x)\psi (x)`$ in Eq. (13), a direct calculation shows

$`[{\displaystyle \frac{1}{2}}{\displaystyle \underset{x^{}}{}}\overline{\psi }(x^{})\psi (x^{}),\mathrm{\Pi }_{f_1f_2}(x)]`$ $`=`$ $`{\displaystyle \frac{1}{4\sqrt{\overline{v}}}}(\psi _{f_1}^{}(x)\gamma _4(1\gamma _4)\gamma _5\psi _{f_2}(x)`$ (A16)
$`\psi _{f_1}^{}(x)\gamma _4(1\gamma _4)\gamma _5\gamma _4\psi _{f_2}(x))`$
$`=`$ $`{\displaystyle \frac{1}{2\sqrt{\overline{v}}}}\psi _{f_1}^{}(x)(1\gamma _4)\gamma _5\psi _{f_2}(x)=\mathrm{\Pi }_{f_1f_2}(x).`$ (A17)
Similarly, one can also show

$`[{\displaystyle \frac{1}{2}}{\displaystyle \underset{x^{}}{}}\overline{\psi }(x^{})\psi (x^{}),\mathrm{\Pi }_{f_2f_1}^{}(x)]`$ $`=`$ $`\mathrm{\Pi }_{f_2f_1}^{}(x),`$ (A18)
$`[{\displaystyle \frac{1}{2}}{\displaystyle \underset{x^{}}{}}\overline{\psi }(x^{})\psi (x^{}),V_{jf_1f_2}(x)]`$ $`=`$ $`V_{jf_1f_2}(x),`$ (A19)
$`[{\displaystyle \frac{1}{2}}{\displaystyle \underset{x^{}}{}}\overline{\psi }(x^{})\psi (x^{}),V_{jf_2f_1}^{}(x)]`$ $`=`$ $`V_{jf_2f_1}^{}(x),`$ (A20)
Eq. (A17) and Eq. (A20) imply that the bilinear $`_x\overline{\psi }(x)\psi (x)`$ can be bosonized as

$`{\displaystyle \underset{x}{}}\overline{\psi }(x)\psi (x){\displaystyle \underset{x}{}}\left(\overline{v}+2\mathrm{\Pi }_{f_2f_1}^{}(x)\mathrm{\Pi }_{f_1f_2}(x)+2{\displaystyle \underset{j}{}}V_{jf_2f_1}^{}(x)V_{jf_1f_2}(x)\right).`$ (A21)
For the four fermion term $`_x\overline{\psi }_{f_1}(x)\psi _{f_2}(x)\overline{\psi }_{f_2}(x+\widehat{k})\psi _{f_1}(x+\widehat{k})`$ in Eq. (13), we have

$`[{\displaystyle \frac{1}{2}}{\displaystyle \underset{x^{}}{}}\overline{\psi }_{f_1^{}}(x^{})\psi _{f_2^{}}(x^{})\overline{\psi }_{f_2^{}}(x^{}+\widehat{k})\psi _{f_1^{}}(x^{}+\widehat{k}),\mathrm{\Pi }_{f_1f_2}(x)]`$ (A22)
$`=`$ $`{\displaystyle \underset{x^{}}{}}\overline{\psi }_{f_1^{}}(x^{})\psi _{f_2^{}}(x^{})[{\displaystyle \frac{1}{2}}\overline{\psi }_{f_2^{}}(x^{}+\widehat{k})\psi _{f_1^{}}(x^{}+\widehat{k}),\mathrm{\Pi }_{f_1f_2}(x)]`$ (A24)
$`+{\displaystyle \underset{x^{}}{}}[{\displaystyle \frac{1}{2}}\overline{\psi }_{f_1^{}}(x^{})\psi _{f_2^{}}(x^{}),\mathrm{\Pi }_{f_1f_2}(x)]\overline{\psi }_{f_2^{}}(x^{}+\widehat{k})\psi _{f_1^{}}(x^{}+\widehat{k})`$
$``$ $`2\overline{v}\mathrm{\Pi }_{f_1f_2}(x),`$ (A25)
where in the r.h.s., Eq. (A17) and the linearization prescription have been used. In addition, we also have

$`[{\displaystyle \frac{1}{2}}{\displaystyle \underset{x^{}}{}}\overline{\psi }_{f_1^{}}(x^{})\psi _{f_2^{}}(x^{})\overline{\psi }_{f_2^{}}(x^{}+\widehat{k})\psi _{f_1^{}}(x^{}+\widehat{k}),\mathrm{\Pi }_{f_2f_1}^{}(x)]`$ $``$ $`2\overline{v}\mathrm{\Pi }_{f_2f_1}^{}(x),`$ (A26)
$`[{\displaystyle \frac{1}{2}}{\displaystyle \underset{x^{}}{}}\overline{\psi }_{f_1^{}}(x^{})\psi _{f_2^{}}(x^{})\overline{\psi }_{f_2^{}}(x^{}+\widehat{k})\psi _{f_1^{}}(x^{}+\widehat{k}),V_{jf_1f_2}(x)]`$ $``$ $`2\overline{v}V_{jf_1f_2}(x),`$ (A27)
$`[{\displaystyle \frac{1}{2}}{\displaystyle \underset{x^{}}{}}\overline{\psi }_{f_1^{}}(x^{})\psi _{f_2^{}}(x^{})\overline{\psi }_{f_2^{}}(x^{}+\widehat{k})\psi _{f_1^{}}(x^{}+\widehat{k}),V_{jf_2f_1}^{}(x)]`$ $``$ $`2\overline{v}V_{jf_2f_1}^{}(x).`$ (A28)
Eq. (A25) and Eq. (A28) also imply that the four fermion term $`_x\overline{\psi }_{f_1}(x)\psi _{f_2}(x)\overline{\psi }_{f_2}(x+\widehat{k})\psi _{f_1}(x+\widehat{k})`$ can be bosonized as

$`{\displaystyle \underset{x}{}}\overline{\psi }_{f_1}(x)\psi _{f_2}(x)\overline{\psi }_{f_2}(x+\widehat{k})\psi _{f_1}(x+\widehat{k})`$ (A29)
$``$ $`{\displaystyle \underset{x}{}}\left(\overline{v}_2+4\overline{v}\mathrm{\Pi }_{f_2f_1}^{}(x)\mathrm{\Pi }_{f_1f_2}(x)+4\overline{v}{\displaystyle \underset{j}{}}V_{jf_2f_1}^{}(x)V_{jf_1f_2}(x)\right).`$ (A30)
In the linear prescription, the four fermion operators in Eq. (13), which are irrelevant to the operators in Eq. (22) are replaced by their vacuum expectation values. For example,

$`\psi ^{}(x)\psi (x)`$ $``$ $`v^{},`$ (A31)
$`{\displaystyle \underset{x,k}{}}\psi _{f_1}^{}(x)\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\psi _{f_1}(x+\widehat{k})`$ $``$ $`2dN_fN_sv_2^{},`$ (A32)
$`{\displaystyle \underset{x,k,j}{}}\left(r^2+\left(12\delta _{|k|,j}\right)\right)\overline{\psi }_{f_1}(x)\sigma _j\psi _{f_2}(x)\overline{\psi }_{f_2}(x+\widehat{k})\sigma _j\psi _{f_1}(x+\widehat{k})`$ $``$ $`2dN_fN_s\overline{v}_{2\sigma _j},`$ (A33)
$`{\displaystyle \underset{x,k,j}{}}\left(r^2\left(12\delta _{|k|,j}\right)\right)\psi _{f_1}^{}(x)\sigma _j\psi _{f_2}(x)\psi _{f_2}^{}(x+\widehat{k})\sigma _j\psi _{f_1}(x+\widehat{k})`$ $``$ $`2dN_fN_sv_{2\sigma _j}^{}.`$ (A34)
Collecting the above results, we obtain Eqs. (25), (27) and (33).

## B VACUUM ENERGY IN THE LARGE $`N_c`$ LIMIT

According to Eq. (38),

$`\overline{v}_2`$ $`=`$ $`{\displaystyle \frac{1}{2dN_fN_s}}{\displaystyle \underset{x,k,f_1,f_2}{}}(\overline{\psi }_{f_1,c_1}(x)\psi _{f_2,c_1}(x)_{eff}\overline{\psi }_{f_2,c_2}(x+\widehat{k})\psi _{f_1,c_2}(x+\widehat{k})_{eff}`$ (B1)
$`+`$ $`\psi _{f_2}(x)\overline{\psi }_{f_2,c_2}(x+\widehat{k})_{eff}\psi _{f_1}(x+\widehat{k})\overline{\psi }_{f_1,c_1}(x)_{eff}),`$ (B2)
$`v_2^{}`$ $`=`$ $`{\displaystyle \frac{1}{2dN_fN_s}}{\displaystyle \underset{x,k,f_1,f_2}{}}(\psi _{f_1,c_1}^{}(x)\psi _{f_2,c_1}(x)_{eff}\psi _{f_2,c_2}^{}(x+\widehat{k})\psi _{f_1,c_2}(x+\widehat{k})_{eff}`$ (B3)
$`+`$ $`\psi _{f_2}(x)\psi _{f_2,c_2}^{}(x+\widehat{k})_{eff}\psi _{f_1}(x+\widehat{k})\psi _{f_1,c_1}^{}(x)_{eff}),`$ (B4)
where the color and flavor indexes for the fermion fields are explicitly specified and summation over the repeated color index is implied. $`\mathrm{}_{eff}`$ stands for the expectation value taken in vacuum state $`|\mathrm{\Omega }_{eff}`$. Because

$`\overline{\psi }_{f_1,c_1}(x)\psi _{f_2,c_1}(x)_{eff}`$ $`=`$ $`\overline{v}\delta _{f_1,f_2}N_c,`$ (B5)
$`\psi _{f_2,c_1}(x)\overline{\psi }_{f_2,c_2}(x+\widehat{k})_{eff}`$ $``$ $`\delta _{c_1,c_2},`$ (B6)
$`\psi _{f_1,c_1}^{}(x)\psi _{f_2,c_1}(x)_{eff}`$ $`=`$ $`v^{}\delta _{f_1,f_2}N_c,`$ (B7)
$`\psi _{f_2,c_1}(x)\psi _{f_2,c_2}^{}(x+\widehat{k})_{eff}`$ $``$ $`\delta _{c_1,c_2},`$ (B8)
terms from Wick-contraction of bilinear at different lattice sites are subdominant for large $`N_c`$. Therefore, we have

$`\overline{v}_2\overline{v}^2,v_2^{}\left(v^{}\right)^2,`$ (B9)
with $`\overline{v}`$ and $`v^{}`$ expectation value of fermion bilinears defined in Eq. (24). Expectation values of other four fermion operators in Eq. (38) can be considered in the similar way,

$`\overline{v}_{2\sigma _j}\left(\mathrm{Tr}\sigma _j\right)^2=0,v_{2\sigma _j}^{}\left(\mathrm{Tr}\sigma _j\right)^2=0`$ (B10)
in the large $`N_c`$ limit.

In the chiral limit Eq. (54), $`G_1`$ defined in Eq. (47) scales as $`\overline{v}K/N_cK`$ and $`G_1/G_2`$ doesn’t depend on $`N_c`$. The $`N_f^2`$ terms in Eq. (64) are subdominant because they are proportional to $`K`$, while other non-vanishing terms of $`E_\mathrm{\Omega }^{(0)}`$ in Eq. (64), defined in Eq. (27) are dominant because they behave as $`\overline{v}KKN_c`$ or $`v^{}KKN_c`$. Using this fact and substituting Eqs. (B9) and (B10) into Eq. (27), we obtain Eq. (LABEL:vacuum\_energy\_large\_Nc) from Eq. (64), i.e., the vacuum energy in the large $`N_c`$ limit.

## C MEAN-FIELD APPROXIMATION IN THE LARGE $`N_c`$ LIMIT

The mean-field approximation is a popular technique wildly used in quantum field theory with four fermion interactions (for example, to bilinearize the Nambu-Jona-Lasino model) and quantum theory of many-particle systems.

In order to show how to derive Eq. (67) from Eq. (13), let’s look at the first and last of the four fermion terms in Eq. (13) as an example. Under the mean-field approximation, i.e., replacing a pair of fermion fields by their vacuum expectation value, the first four fermion term in Eq. (13) becomes

$`{\displaystyle \underset{x,k,f_1,f_2}{}}\psi _{f_1,c_1}^{}(x)\psi _{f_2,c_1}(x)\psi _{f_2,c_2}^{}(x+\widehat{k})\psi _{f_1,c_2}(x+\widehat{k})`$ (C1)
$``$ $`{\displaystyle \underset{x,k,f_1,f_2}{}}(\psi _{f_1,c_1}^{}(x)\psi _{f_2,c_1}(x)_{eff}\psi _{f_2,c_2}^{}(x+\widehat{k})\psi _{f_1,c_2}(x+\widehat{k})`$ (C5)
$`+\psi _{f_1,c_1}^{}(x)\psi _{f_2,c_1}(x)\psi _{f_2,c_2}^{}(x+\widehat{k})\psi _{f_1,c_2}(x+\widehat{k})_{eff}`$
$`+\psi _{f_1,c_1}^{}(x)\psi _{f_2,c_1}(x)\psi _{f_2,c_2}^{}(x+\widehat{k})_{eff}\psi _{f_1,c_2}(x+\widehat{k})`$
$`+\psi _{f_2,c_1}(x)\psi _{f_2,c_2}^{}(x+\widehat{k})\psi _{f_1,c_1}^{}(x)\psi _{f_1,c_2}(x+\widehat{k})_{eff})2dN_fN_sv_2^{}.`$
Again, the color and flavor indexes for the fermion fields are explicitly specified and summation over the repeated color index is implied. According to Eq. (B8), for large $`N_c`$, Wick contractions of bilinear at different lattice sites give subdominant contribution, therefore

$`{\displaystyle \underset{x,k,f_1,f_2}{}}\psi _{f_1,c_1}^{}(x)\psi _{f_2,c_1}(x)\psi _{f_2,c_2}^{}(x+\widehat{k})\psi _{f_1,c_2}(x+\widehat{k})`$ (C6)
$``$ $`v^{}{\displaystyle \underset{x,k,f_1,f_2}{}}\delta _{f_1,f_2}\left(\psi _{f_2,c}^{}(x+\widehat{k})\psi _{f_1,c}(x+\widehat{k})+\psi _{f_1,c}^{}(x)\psi _{f_2,c}(x)\right)2dN_fN_s\left(v^{}\right)^2`$ (C7)
$`=`$ $`v^{}{\displaystyle \underset{x,k,f}{}}\left(\psi _{f,c}^{}(x)\psi _{f,c}(x)+\psi _{f,c}^{}(x+\widehat{k})\psi _{f,c}(x+\widehat{k})\right)2dN_fN_s\left(v^{}\right)^2`$ (C8)
$`=`$ $`4dv^{}{\displaystyle \underset{x}{}}\psi ^{}(x)\psi (x)2dN_fN_s\left(v^{}\right)^2.`$ (C9)
Under the mean-field approximation in the large $`N_c`$ limit, the dominant contribution to the last four fermion term in Eq. (13) is

$`{\displaystyle \underset{x,k,f_1,f_2}{}}\psi _{f_1,c_1}^{}(x)\sigma _j\psi _{f_2,c_1}(x)\psi _{f_2,c_2}^{}(x+\widehat{k})\sigma _j\psi _{f_1,c_2}(x+\widehat{k})`$ (C10)
$``$ $`{\displaystyle \underset{x,k,f_1,f_2}{}}(\psi _{f_1,c_1}^{}(x)\sigma _j\psi _{f_2,c_1}(x)_{eff}\psi _{f_2,c_2}^{}(x+\widehat{k})\sigma _j\psi _{f_1,c_2}(x+\widehat{k})`$ (C12)
$`+\psi _{f_1,c_1}^{}(x)\sigma _j\psi _{f_2,c_1}(x)\psi _{f_2,c_2}^{}(x+\widehat{k})\sigma _j\psi _{f_1,c_2}(x+\widehat{k})_{eff})2dN_fN_sv_{2\sigma _j}^{}`$
$``$ $`Tr(\sigma _j)=0.`$ (C13)
Treating other four fermion terms in the same way, one obtains Eq. (67), i.e., the bilinear Hamiltonian $`H_{MFA}`$.

## D NUCLEON MASS

Under the mean-field approximation in the large $`N_c`$ limit, Eq. (86) becomes

$`M_{Nucl}=\mathrm{\Omega }_{eff}|O_{Nucl}\left(H_{MFA}\mu N\right)O_{Nucl}^{}|\mathrm{\Omega }_{eff}E_\mathrm{\Omega }.`$ (D1)
Substituting the bilinear Hamiltonian Eq. (67), the fermion field Eq. (73), and the particle number operator Eq. (16) into Eq. (D1), we have

$`M_{Nucl}`$ $`=`$ $`\mathrm{\Omega }_{eff}|O_{Nucl}\left(A{\displaystyle \underset{x}{}}\overline{\psi }(x)\psi (x)+(B\mu ){\displaystyle \underset{x}{}}\psi ^{}(x)\psi (x)+C\right)O_{Nucl}^{}|\mathrm{\Omega }_{eff}E_\mathrm{\Omega }`$ (D2)
$`=`$ $`\mathrm{\Omega }_{eff}|O_{Nucl}\left(A{\displaystyle \underset{x}{}}(\xi ^{}\xi \eta \eta ^{})+(B\mu ){\displaystyle \underset{x}{}}(\xi ^{}\xi +\eta \eta ^{})+C\right)O_{Nucl}^{}|\mathrm{\Omega }_{eff}E_\mathrm{\Omega }`$ (D3)
$`=`$ $`\mathrm{\Omega }_{eff}|O_{Nucl}\left(A{\displaystyle \underset{x}{}}\xi ^{}\xi +(B\mu ){\displaystyle \underset{x}{}}\xi ^{}\xi \right)O_{Nucl}^{}|\mathrm{\Omega }_{eff},`$ (D4)
where we have used the fact that the terms with $`\eta ^{}\eta `$ and $`C`$ are canceled by $`E_\mathrm{\Omega }`$.

Let us take the proton as an example. One can write the operator $`O_{Nucl}^{}`$ explicitly as

$`O_{Nucl}^{}={\displaystyle \frac{1}{\sqrt{18N_s}}}{\displaystyle \underset{x}{}}ϵ_{c_1c_2c_3}\xi _{c_1,u,1}^{}(x)\left(\xi _{c_2,u,1}^{}(x)\xi _{c_3,d,2}^{}(x)\xi _{c_2,u,2}^{}(x)\xi _{c_3,d,1}^{}(x)\right),`$ (D5)
where $`c_1,c_2`$ and $`c_3`$ are the color indexes, $`u`$ and $`d`$ stand respectively for the u-quark and d-quark, and $`1`$ and $`2`$ are the spin up and down indexes. Using the anti-commutation relationship for fermions

$`\xi _{c,f,s}(x)\xi _{c^{},f^{},s^{}}^{}(x^{})=\delta _{c,c^{}}\delta _{f,f^{}}\delta _{s,s^{}}\delta _{x,x^{}}\xi _{c^{},f^{},s^{}}^{}(x^{})\xi _{c,f,s}(x),`$ (D6)
one obtains

$`{\displaystyle \underset{x}{}}\xi ^{}(x)\xi (x)O_{Nucl}^{}|\mathrm{\Omega }_{eff}=3O_{Nucl}^{}|\mathrm{\Omega }_{eff}.`$ (D7)
Eq. (87) is a consequence of Eq. (D4) and Eq. (D7).

At $`\mu =0`$, the nucleon mass equals to $`3m_{dyn}^{(0)}`$ (see Eq. (94)), which agrees with Refs. .