format: text

text:
# 1 Introduction

## 1 Introduction

Interest in QCD at low temperature with a non-zero background baryon density is currently at a high level. Understanding the ground state of any strongly interacting many-body system requires non-perturbative techniques, and lattice QCD is usually considered one of the most systematic and successful. Unfortunately, attempts to develop Monte Carlo simulation techniques for systems with a non-vanishing baryon chemical potential $`\mu `$ have so far met with little success ; understandably, therefore, there has been comparatively little parallel development in relating correlation functions measured using standard lattice field theory techniques to the phenomenology of dense systems. Exceptions are recent attempts to measure diquark condensates $`qq`$ in Two Color QCD (TCQCD) and Nambu – Jona-Lasinio (NJL) models (see eg. ) inspired by the possibility that the QCD ground state at large $`\mu `$ and small $`T`$ is a color superconductor . However, conceptually simpler issues such as the restoration of chiral symmetry and the spectrum of excitations in a dense interacting medium continue to present interesting challenges. One obvious example is whether and how masses and widths of mesonic bound states such as the $`\rho `$ meson are modified in a dense medium . A broader question is how a Fermi surface manifests itself in Euclidean simulations. Recall that color superconductivity phenomenology is based on the BCS mechanism for electronic superconductivity, which in turn is based on a small, though non-perturbative, deformation of a sharp Fermi surface.

As part of the learning process we have studied the simplest non-trivial model simulable with $`\mu 0`$ using standard algorithms, namely the Gross-Neveu (GN) model in 2+1 dimensions. Its Lagrangian density in Euclidean metric is written in terms of $`4N_f`$-component spinors $`\psi ,\overline{\psi }`$ as

$$=\overline{\psi }(/+m)\psi \frac{g^2}{2N_f}(\overline{\psi }\psi )^2.$$
(1)
In the chiral limit $`m=0`$ the model has a global Z2 symmetry $`\psi \gamma _5\psi `$, $`\overline{\psi }\overline{\psi }\gamma _5`$. At tree level the model’s content is $`N_f`$ fermion flavors each with parity-invariant mass $`m`$, interacting via a four-fermion contact term. However at leading order in an expansion in $`1/N_f`$ it can be shown that for sufficiently large coupling strength $`g^2`$ the Z2 symmetry is spontaneously broken by a condensate $`\overline{\psi }\psi 0`$ leading to a dynamically generated fermion mass gap $`\mathrm{\Sigma }_0=g^2\overline{\psi }\psi m`$. The critical coupling $`g_c^2`$ at which the gap $`\mathrm{\Sigma }_0/\mathrm{\Lambda }_{UV}0`$ defines an ultra-violet stable fixed point of the renormalisation group at which an interacting continuum limit may be taken. This picture has been verified both at next-to-leading order in $`1/N_f`$ and by Monte Carlo simulations with finite $`N_f`$ (eg. ).

The thermodynamic phase structure of the GN model is also known to leading order in $`1/N_f`$ . At $`T=0`$, $`\overline{\psi }\psi `$ remains constant as $`\mu `$ is increased until the Z2 symmetry is restored in a first-order transition at a critical $`\mu _c=\mathrm{\Sigma }_0`$. At the same point the baryon density $`nN_f^1\overline{\psi }\gamma _0\psi `$ jumps from zero to $`\mu _c^2/2\pi `$ and then continues to rise quadratically as the Fermi disk of relativistic “quark matter” grows. This point is hence often called the onset. Simulations with finite $`N_f`$ have again verified this picture with , and importantly have showed that the first order nature of the chiral transition at $`\mu _c\stackrel{<}{}\mathrm{\Sigma }_0`$ persists at small but non-zero $`T`$, implying the existence of a tricritical point in the $`(\mu ,T)`$ plane . With currently available precision there is however no evidence for a “nuclear matter” phase with both $`\overline{\psi }\psi 0,n0`$.

The high-$`\mu `$ phase of a related model, the 2+1$`d`$ NJL model, which has a global SU(2)$`{}_{L}{}^{}`$SU(2)R chiral symmetry, has also been extensively studied . The main issue of interest has been whether diquark condensation takes place at the Fermi surface in the isopsin-symmetric degenerate matter formed for $`\mu >\mu _c`$, leading to a spontaneous breakdown of baryon number symmetry whose physical manifestation is superfluidity. Diquark pairing is favoured in this model because the formulation with staggered lattice fermions admits a scalar isoscalar $`qq`$ operator localised on a single lattice site . It should be noted that in the current Z2 model a local $`O(N_f)`$ invariant diquark operator is forbidden by the Pauli exclusion principle. Although we have not studied the diquark sector in the current work, we will see in Sec. 3.1 below that the two models show dramatic differences in the spectrum of spin-$`\frac{1}{2}`$ excitations around the Fermi surface.

In Sec. 2.1 we will use the technology of the $`1/N_f`$ expansion in the so-called “hard dense loop” approach to calculate the long-range inter-particle interaction in medium, ie. for $`\mu >\mu _c`$. The main result is that the Debye mass is infinite, but the plasma frequency $`M_\sigma `$ is finite and actually vanishes as the onset $`\mu \mu _{c+}`$ is approached. The main new feature emerging for $`\mu >\mu _c`$ is the existence of a new physical scale, the Fermi momentum $`k_F`$. Since the lowest energy excitations of the ground state have $`|\stackrel{}{k}|k_F`$, measurements of Euclidean timeslice correlators with $`\stackrel{}{k}0`$ are mandatory. The resulting energy $`E(\stackrel{}{k})`$ is known as the dispersion relation. The main numerical results of the paper are thus dispersion relations in various channels of interest. In the spin-$`\frac{1}{2}`$ fermion channel, the low-energy excitations about the Fermi surface, which may be hole- or particle-like, are generically known as quasiparticles. There is a well-developed framework for describing quasiparticles and their interactions known as Fermi liquid theory. In Sec. 2.2 we will present theoretical predictions of Fermi liquid properties to leading non-trivial order in $`1/N_f`$. In Sec. 2.3 we will describe another consequence of the new scale; the wavefunction describing inter-particle spatial correlations, which decays monotonically in the vacuum, now oscillates with spatial frequency related to $`k_F`$.

Sec. 3 describes the results of our numerical study. As mentioned above, the fermion dispersion relation is presented in Sec. 3.1, and mesonic (ie. $`q\overline{q}`$) channels are similarly investigated in Secs. 3.2 and 3.3. We will find that the predictions of the $`1/N_f`$ approach are verified with reasonable accuracy, and that in-medium effects are clearly visible. In most meson channels the lightest excitations are particle-hole pairs with effectively zero energy whose temporal decay is hence algebraic rather than exponential, but which nonetheless show a non-trivial dependence on the ratio $`|\stackrel{}{k}|/\mu `$. An interesting phenomenon not described within the $`1/N_f`$ approach is a collective excitation of the whole system (ie. not a quasiparticle) corresponding to a distortion in the shape of the Fermi surface which propagates as zero sound. Tentative evidence for an excitation of this form is presented in Sec. 3.3. Finally in Sec. 3.4 we present results for wavefunctions in mesonic channels which clearly show the oscillatory behaviour described above. After a brief summary, various technical details of the calculation of relevant correlators in free-field theory are given in two appendices.

## 2 Theoretical Results

In most theoretical approaches to the GN model it is convenient to encode interactions by the introduction of an auxiliary scalar field $`\sigma `$ related to the fermions via the equation of motion $`\sigma =\frac{g}{\sqrt{N_f}}\overline{\psi }\psi `$. Although at tree level $`\sigma `$ is non-propagating, the leading order $`1/N_f`$ expansion at $`\mu =0`$ predicts that it acquires dynamics through quantum corrections due to virtual $`q\overline{q}`$ pairs, resulting in a propagator $`D_\sigma (k^2)`$ of the form

$$D_\sigma (k^2)\{\begin{array}{cc}(k^2+4\mathrm{\Sigma }_0^2)^1,\hfill & k\mathrm{\Sigma }_0\text{;}\hfill \\ (k^2)^{\frac{1}{2}},\hfill & k\mathrm{\Sigma }_0\text{.}\hfill \end{array}$$
(2)
For soft momenta $`\sigma `$ thus resmbles an orthodox scalar boson formed as a weakly-bound $`q\overline{q}`$ state. In the Born approximation $`\sigma `$ exchange generates an interaction between the fermions of range $`O(\mathrm{\Sigma }_0^1)`$ and strength $`O(N_f^1)`$. By contrast, as $`k\mathrm{}`$ $`D_\sigma `$ is a harder function of $`k`$, though still softer than the original tree-level prediction $`D_\sigma k^0`$. The crossover from IR to UV is encoded as a branch cut in the complex $`k`$ plane, which means that the decay of $`D_\sigma `$ in Euclidean time is not well-described by a simple pole. In the next subsection we will extend the leading order calculation of $`D_\sigma `$ to $`\mu >\mu _c`$, restricting our attention to soft momenta $`k\mu `$. Even with this restriction the resulting dispersion relation is surprisingly complicated. The calculation also gives access to the forward scattering amplitude between quasiparticles in the medium, which in turn enables contact with the relativistic theory of Fermi liquids . It is hence possible to find quantitative relations between the Fermi liquid parameters $`k_F`$, $`\epsilon _F`$ and $`\beta _F`$, respectively the Fermi momentum, energy and velocity.

### 2.1 The $`\sigma `$ Propagator in Medium

Here we examine the auxiliary scalar propagator $`D_\sigma `$ in the presence of a non-zero baryon density. To leading order in the $`1/N_f`$ expansion (eg. ), we have

$`D_\sigma ^1(k)`$ $`=`$ $`1\mathrm{\Pi }(k;\mu )`$ (3)
$`=`$ $`g^2\left[{\displaystyle \frac{1}{\mathrm{\Sigma }_0}}{\displaystyle _{q,\mu =0}}\text{tr}{\displaystyle \frac{1}{iq/+\mathrm{\Sigma }_0}}+{\displaystyle _{q,\mu >\mu _c}}\text{tr}{\displaystyle \frac{1}{iq/+\mu \gamma _0}}{\displaystyle \frac{1}{i(q/k/)+\mu \gamma _0}}\right],`$
where $`\mathrm{\Pi }`$ is the virtual fermion – anti-fermion vacuum polarisation bubble, we have used the gap equation at $`\mu =0`$ to express $`1/g^2`$ in terms of the zero density gap $`\mathrm{\Sigma }_0`$, and assumed that the gap vanishes in the integral defining $`\mathrm{\Pi }`$ for $`\mu >\mu _c`$. It is technically more convenient to evaluate the integrals for temperature $`T>0`$ and for arbitrary spatial dimension $`d(1,3)`$; after tracing over Dirac indices we find

$`D_\sigma ^1(k_0,\stackrel{}{k})=4g^2T{\displaystyle \underset{n}{}}{\displaystyle }{\displaystyle \frac{d^d\stackrel{}{q}}{(2\pi )^d}}\{{\displaystyle \frac{1}{\omega _n^2+\stackrel{}{q}^{\mathrm{\hspace{0.17em}2}}+\mathrm{\Sigma }_0^2}}{\displaystyle \frac{1}{(\omega _ni\mu )^2+\stackrel{}{q}^{\mathrm{\hspace{0.17em}2}}}}`$
$`{\displaystyle \frac{k_0(\omega _nk_0i\mu )+\stackrel{}{k}.(\stackrel{}{q}\stackrel{}{k})}{[(\omega _ni\mu )^2+\stackrel{}{q}^{\mathrm{\hspace{0.17em}2}}][(\omega _nk_0i\mu )^2+(\stackrel{}{q}\stackrel{}{k})^2]}}\},`$ (4)
where the sum is over discrete Matsubara modes $`\omega _n=(2n1)\pi T`$.

The first two terms of (4) are individually divergent, but their sum is finite. In the limit $`T0`$ (see eg. ) their contribution is

$$\frac{4g^2}{(4\pi )^{\frac{d}{2}}\mathrm{\Gamma }(\frac{d}{2})}\left\{_0^\mathrm{\Lambda }𝑑q\frac{q^{d1}}{\sqrt{q^2+\mathrm{\Sigma }_0^2}}_\mu ^\mathrm{\Lambda }𝑑qq^{d2}\right\}=\frac{4g^2(\mu ^{d1}\mu _c^{d1})}{(4\pi )^{\frac{d}{2}}\mathrm{\Gamma }(\frac{d}{2})(d1)},$$
(5)
the critical chemical potential in $`d+1`$ dimensions in the large-$`N_f`$ limit being

$$\mu _c=\left(\frac{(1d)\mathrm{\Gamma }(\frac{d}{2})\mathrm{\Gamma }(\frac{1}{2}\frac{d}{2})}{2\sqrt{\pi }}\right)^{\frac{1}{d1}}\mathrm{\Sigma }_0.$$
(6)
The frequency sum for the momentum-dependent term may be evaluated using the formulæ(eg. )

$`T`$ $`{\displaystyle \underset{n}{}}\mathrm{\Delta }^+(\omega _n,E_1)\mathrm{\Delta }^{}(\omega \omega _n,E_2)={\displaystyle \underset{s_i=\pm 1}{}}{\displaystyle \frac{1f_+(s_1E_1)f_{}(s_2E_2)}{4s_1s_2E_1E_2(i\omega s_1E_1s_2E_2)}};`$
$`iT`$ $`{\displaystyle \underset{n}{}}\omega _n\mathrm{\Delta }^+(\omega _n,E_1)\mathrm{\Delta }^{}(\omega \omega _n,E_2)={\displaystyle \underset{s_i=\pm 1}{}}{\displaystyle \frac{(\mu s_1E_1)(1f_+(s_1E_1)f_{}(s_2E_2))}{4s_1s_2E_1E_2(i\omega s_1E_1s_2E_2)}},`$
where $`\mathrm{\Delta }^\pm (\omega ,E)((\omega i\mu )^2+E^2)^1`$, the Fermi-Dirac function $`f_\pm (x)(\mathrm{exp}((x\mu )/T)+1)^1`$, and $`E_i`$ is taken positive throughout. Using the identity $`1f_\pm (x)f_{}(x)=0`$, and noting that for $`\mu >0`$ $`lim_{T0}f_+(E)=1\theta (E\mu )`$, $`lim_{T0}f_{}(E)=0`$, we find in this limit

$`g^2{\displaystyle }{\displaystyle \frac{d^d\stackrel{}{q}}{(2\pi )^d}}{\displaystyle \frac{1}{E_1E_2}}\{[k_0^2+`$ $`\stackrel{}{k}.(\stackrel{}{q}\frac{1}{2}\stackrel{}{k})iE_1k_0]({\displaystyle \frac{\theta (E_1\mu )}{ik_0E_1E_2}}+`$
$`{\displaystyle \frac{\theta (E_2\mu )\theta (E_1\mu )}{ik_0E_1+E_2}})`$ $``$ $`[k_0^2+\stackrel{}{k}.(\stackrel{}{q}\frac{1}{2}\stackrel{}{k})+iE_1k_0]{\displaystyle \frac{\theta (E_2\mu )}{ik_0+E_1+E_2}}\}`$ (8)
where the spatial loop momentum has been rerouted by $`\stackrel{}{q}\stackrel{}{q}+\frac{1}{2}\stackrel{}{k}`$ so that $`E_{1,2}=|\stackrel{}{q}\pm \frac{1}{2}\stackrel{}{k}|`$. We now restrict our attention to $`k_0,|\stackrel{}{k}|\mu `$, thus probing the inter-particle interaction for soft momentum transfer in the so-called hard dense loop (HDL) approximation 111For asymptotically large $`k\mu ,\mathrm{\Sigma }_0`$ the form of the propagator reverts to $`D_\sigma ^1(k)k^{d1}`$ .. The first and third terms together yield

$`{\displaystyle \frac{4g^2}{(4\pi )^{\frac{d}{2}}\mathrm{\Gamma }(\frac{d}{2})}}{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\theta }{2\pi }}{\displaystyle _{\mu \frac{|\stackrel{}{k}|}{2}\mathrm{cos}\theta +\mathrm{}}^{\mathrm{}}}𝑑qq^{d2}{\displaystyle \frac{k_0^2+|\stackrel{}{k}|^2+ik_0|\stackrel{}{k}|\mathrm{cos}\theta +O(k^3\mathrm{cos}\theta /q,k^4/q^2)}{k_0^2+4q^2}}`$
$`={\displaystyle \frac{g^2\mu ^{d3}}{(4\pi )^{\frac{d}{2}}\mathrm{\Gamma }(\frac{d}{2})}}\left({\displaystyle \frac{k_0^2+|\stackrel{}{k}|^2}{3d}}+i{\displaystyle \frac{k_0}{4\mu }}|\stackrel{}{k}|^2+O(k^4/\mu ^2)\right)`$ (9)
The integrand of the second term is non-zero only in a crescent-shaped region near the Fermi surface; we find

$`g^2{\displaystyle \frac{|\stackrel{}{k}|}{(4\pi )^{\frac{d}{2}}\mathrm{\Gamma }(\frac{d}{2})}}{\displaystyle _{\frac{\pi }{2}}^{\frac{\pi }{2}}}{\displaystyle \frac{d\theta }{2\pi }}{\displaystyle \frac{\mathrm{cos}\theta }{k_0^2+|\stackrel{}{k}|^2\mathrm{cos}^2\theta }}{\displaystyle _{\mu \frac{|\stackrel{}{k}|}{2}\mathrm{cos}\theta +\mathrm{}}^{\mu +\frac{|\stackrel{}{k}|}{2}\mathrm{cos}\theta +\mathrm{}}}𝑑qq^{d3}`$
$`\times \left(2(k_0^2+|\stackrel{}{k}|^2)+i{\displaystyle \frac{k_0}{q}}|\stackrel{}{k}|^2\mathrm{cos}\theta \mathrm{sin}^2\theta +O(k^4/q^2)\right)`$ $`=`$ (10)
$`{\displaystyle \frac{g^2\mu ^{d3}}{(4\pi )^{\frac{d}{2}}\mathrm{\Gamma }(\frac{d}{2})}}{\displaystyle \frac{|\stackrel{}{k}|^2}{k_0+(k_0^2+|\stackrel{}{k}|^2)^{\frac{1}{2}}}}[(k_0^2+|\stackrel{}{k}|^2)^{\frac{1}{2}}+i{\displaystyle \frac{k_0|\stackrel{}{k}|^2}{4\mu [k_0+(k_0^2+|\stackrel{}{k}|^2)^{\frac{1}{2}}]}}`$ $`+`$ $`O\left({\displaystyle \frac{k^3}{\mu ^2}}\right)]`$
Eqns. (5), (9) and (10) combined yield

$`D_\sigma ^1(k_0,\stackrel{}{k})={\displaystyle \frac{g^2\mu ^{d3}}{(4\pi )^{\frac{d}{2}}\mathrm{\Gamma }(\frac{d}{2})(3d)}}[4{\displaystyle \frac{(3d)}{(d1)}}\mu ^{3d}(\mu ^{d1}\mu _c^{d1})+`$ (11)
$`k_0^2+|\stackrel{}{k}|^2(1{\displaystyle \frac{(3d)(k_0^2+|\stackrel{}{k}|^2)^{\frac{1}{2}}}{k_0+(k_0^2+|\stackrel{}{k}|^2)^{\frac{1}{2}}}})+i{\displaystyle \frac{k_0|\stackrel{}{k}|^2}{4\mu }}(1+{\displaystyle \frac{|\stackrel{}{k}|^2}{[k_0+(k_0^2+|\stackrel{}{k}|^2)^{\frac{1}{2}}]^2}})+O\left({\displaystyle \frac{k^4}{\mu ^2}}\right)].`$
The first observation to make is that in the static limit $`k_0=0`$ the momentum dependence of (11) vanishes for the physical value $`d=2`$, implying complete screening of the static potential due to $`\sigma `$ exchange, or in other words infinite Debye mass. By contrast, for $`\stackrel{}{k}=\stackrel{}{0}`$ $`D_\sigma `$ is proportional to a conventional zero-momentum boson propagator with mass

$$M_\sigma ^2=4\frac{(3d)}{(d1)}\mu ^{3d}(\mu ^{d1}\mu _c^{d1})$$
(12)
or for the physical $`d=2`$

$$M_\sigma =2\sqrt{\mu (\mu \mu _c)}.$$
(13)
Just above the transition, therefore, there is a tightly-bound state whose decay at rest into a particle - anti-particle pair is Pauli-blocked. The free-field result $`M_\sigma =2\mu `$ is only approached as $`\mu \mathrm{}`$. $`M_\sigma `$ may also be interpreted as the plasma frequency of the system.

For states in motion, things are more complicated. For $`|\stackrel{}{k}||k_0|`$ in $`d=2`$ we have

$$D_\sigma ^1(k_0,\stackrel{}{k})M_\sigma ^2+k_0^2+\frac{1}{2}|\stackrel{}{k}|^2+i\frac{k_0|\stackrel{}{k}|^2}{4\mu }+O\left(\frac{|\stackrel{}{k}|^4}{k_0^2}\right).$$
(14)
Eqn. (14) yields two pure imaginary poles in the complex $`k_0`$ plane leading to solutions with differing decay lengths in forward and backward directions in Euclidean time. In principle this may be probed by a Euclidean lattice simulation. However, linear response theory dictates that on analytic continuation to Minkowski metric, the physically relevant solution giving the energy $`E(\stackrel{}{k})`$ of a collective excitation of the system is given by poles of the retarded propagator , which corresponds to the solution of (14) with negative imaginary part:

$$E(\stackrel{}{k})=\left(M_\sigma ^2+\frac{1}{2}|\stackrel{}{k}|^2\right)^{\frac{1}{2}}+\frac{|\stackrel{}{k}|^2}{8\mu }M_\sigma +|\stackrel{}{k}|^2\left(\frac{1}{4M_\sigma }+\frac{1}{8\mu }\right).$$
(15)
In the limit $`\mu \mathrm{}`$ this becomes $`EM_\sigma +|\stackrel{}{k}|^2/2M_\sigma `$, appropriate for a free non-relativistic particle of mass $`M_\sigma `$. Another limit we may consider is $`\mu \mu _{c+}`$, in which case $`M_\sigma ^2|\stackrel{}{k}|^2k_0^2`$ in (11), yielding the dispersion relation of a standard relativistic boson:

$$E^2(\stackrel{}{k})M_\sigma ^2+|\stackrel{}{k}|^2.$$
(16)
Dispersion curves obtained from a numerical solution of (11) for three different $`M_\sigma `$ are shown in Fig. 1.

### 2.2 Fermi Liquid Properties

Given what we know about the fermion - fermion interaction in medium from the calculation of $`D_\sigma `$, it is possible to make quantitative statements in the context of the Fermi liquid description of degenerate fermions, first developed for non-relativistic systems such as 3He by Landau  (see also ). The generalisation to relativistic systems has subsequently been given by Baym and Chin . The essential physical idea is that the dominant low-energy excitations in the neighbourhood of the Fermi surface are long-lived quasiparticles having energy $`\epsilon _\stackrel{}{k}`$, width $`(\epsilon _\stackrel{}{k}\mu )^2`$ and equilibrium distribution $`n_\stackrel{}{k}`$ related by the Fermi-Dirac function

$$n_\stackrel{}{k}=f_+(\epsilon _\stackrel{}{k}).$$
(17)
For $`T=0`$ we expect $`\epsilon _\stackrel{}{k}`$ to have the form

$$\epsilon _\stackrel{}{k}\mu +\beta _F(|\stackrel{}{k}|k_F)$$
(18)
where $`k_F`$, $`\beta _F|\stackrel{}{}_\stackrel{}{k}\epsilon _{|\stackrel{}{k}|=k_F}|`$ are respectively the Fermi momentum and Fermi velocity, which in the absence of further information are phenomenological parameters. The quasiparticle energy is implicitly determined by the equation for its variation under a small departure from equilibrium;

$$\delta \epsilon _\stackrel{}{k}=\frac{d^dk^{}}{(2\pi )^d}_{\stackrel{}{k},\stackrel{}{k}^{}}\delta n_\stackrel{}{k}^{},$$
(19)
where $`_{\stackrel{}{k},\stackrel{}{k}^{}}=_{\stackrel{}{k}^{},\stackrel{}{k}}`$ is the Fermi liquid interaction, which is in turn related to the two particle forward scattering amplitude

$$_{\stackrel{}{k},\sigma ,\stackrel{}{k}^{},\sigma ^{}}=\frac{1}{4\epsilon _\stackrel{}{k}\epsilon _\stackrel{}{k}^{}}_{\stackrel{}{k},\sigma ,\stackrel{}{k}^{},\sigma ^{}},$$
(20)
where we have retained explicit spin states $`\sigma ,\sigma ^{}`$ and $``$ is the Lorentz invariant matrix element. Landau’s theory becomes predictive if this amplitude is calculable in medium.

To leading non-trivial order in $`1/N_f`$ there are both direct and exchange contributions to $``$, as shown in Fig. 2, but in the chiral limit only the exchange term is non-vanishing. After averaging over spin states we find

$$_{\stackrel{}{k},\stackrel{}{k}^{}}=\frac{g^2}{4N_f}\left[1\frac{\stackrel{}{k}.\stackrel{}{k}^{}}{\epsilon _\stackrel{}{k}\epsilon _\stackrel{}{k}^{}}\right]D_\sigma (\epsilon _\stackrel{}{k}\epsilon _\stackrel{}{k}^{},\stackrel{}{k}\stackrel{}{k}^{}).$$
(21)
Specialising to $`d=2`$ and using the result (11) for $`D_\sigma `$ we find that for $`\stackrel{}{k},\stackrel{}{k}^{}`$ on the Fermi surface separated by angle $`\theta `$,

$$_{\stackrel{}{k},\stackrel{}{k}^{}}=\frac{\pi \mu }{N_fM_\sigma ^2(\mu )}(1\mathrm{cos}\theta )$$
(22)
with $`M_\sigma ^2=4\mu (\mu \mu _c)`$. This particularly simple form arises due to the vanishing spatial momentum dependence of $`D_\sigma `$ for zero energy exchange since the $`\sigma `$ has infinite Debye mass in medium. The factor of $`N_f^1`$ ensures that $``$ can be treated systematically as a weak correction to the leading order case of non-interacting particles; note that this would not hold if direct interactions were involved, since that would result in a compensating “degeneracy factor” of $`N_f`$ in the numerator.

It is possible to derive relations between the Fermi liquid parameters. For instance, by requiring consistency of the above framework under Lorentz boosts the following relation may be derived :

$$\epsilon _\stackrel{}{k}\stackrel{}{}_\stackrel{}{k}\epsilon _\stackrel{}{k}=\stackrel{}{k}+\frac{d^2k^{}}{(2\pi )^2}_{\stackrel{}{k},\stackrel{}{k}^{}}\epsilon _\stackrel{}{k}^{}\stackrel{}{}_\stackrel{}{k}^{}n_\stackrel{}{k}^{}.$$
(23)
Dotting this equation with the unit vector $`\widehat{\stackrel{}{k}}`$, and noting that at $`T=0`$, $`n_\stackrel{}{k}^{}=\theta (\mu \epsilon _\stackrel{}{k}^{})`$ and $`\beta _F\delta (\epsilon _\stackrel{}{k}^{}\mu )=\delta (|\stackrel{}{k}^{}|k_F)`$ we arrive at the following relation at the Fermi surface:

$$\mu \beta _F+\mu k_F𝔤\frac{\pi \mu }{N_fM_\sigma ^2(\mu )}_0^{2\pi }\frac{d\theta }{4\pi ^2}(1\mathrm{cos}\theta )\mathrm{cos}\theta =k_F,$$
(24)
ie.

$$\beta _F=\frac{k_F}{\mu }+𝔤\frac{k_F\mu }{4N_fM_\sigma ^2(\mu )}.$$
(25)
Here $`𝔤`$ is the degeneracy of the states participating in the exchange interaction; for the Z2 GN model $`𝔤=2`$ counts the number of spin states, whereas $`𝔤=4`$ for the isospin-symmetric matter described by the NJL model when $`\pi `$ exchange is also taken into account.

Similarly, the compressibility follows from differentiating the relation $`\mu =\epsilon _{k_F}`$ with respect to $`n`$:

$$\frac{\mu }{n}=\frac{\epsilon _k}{k_F}\frac{k_F}{n}+\overline{}=\beta _F\frac{2\pi }{𝔤k_F}+\frac{\pi \mu }{N_fM_\sigma ^2(\mu )},$$
(26)
where we have used the equality of particle and quasiparticle densities $`n=\frac{d^2\stackrel{}{k}}{(2\pi )^2}n_\stackrel{}{k}`$ and the bar denotes the average over the circle. Using relation (25) we then derive the first sound velocity $`\beta _1`$:

$$\beta _1^2\frac{n}{\mu }\frac{\mu }{n}=\frac{k_F^2}{2\mu ^2}+\frac{3𝔤k_F^2}{8N_fM_\sigma ^2(\mu )}.$$
(27)
To find the relation between $`\epsilon _{k_F}=\mu `$ and $`k_F=(4\pi n/𝒢)^{\frac{1}{2}}`$ we need to integrate (26):

$$_{k_{F_c}}^{k_F}k_F𝑑k_F=_{\mu _c}^\mu \frac{\mu d\mu }{1+\frac{3𝔤}{16N_f}\frac{\mu }{\mu \mu _c}}.$$
(28)
The integration limits specify the conditions at onset; note that $`k_{F_c}=\mu _c+O(N_f^1)=\mathrm{\Sigma }_0+O(N_f^1)`$ with the corrections given by the two-loop gap equation. The solution of (28) expressed consistently to $`O(N_f^1)`$ is

$$k_F^2=\mu ^2+(k_{F_c}^2\mu _c^2)\frac{3𝔤}{16N_f}\left[\mu ^2\mu _c^2+2\mu _c(\mu \mu _c)+2\mu _c^2\mathrm{ln}\left(\frac{\mu \mu _c}{3𝒢\mu _c/16N_f}\right)\right].$$
(29)
We can combine this with (25) to get an absolute prediction for $`\beta _F`$. For $`\mu \mu _c`$ we find

$$\frac{k_F}{\mu }=1\frac{3𝔤}{32N_f};\beta _F=1\frac{𝔤}{32N_f},$$
(30)
which, happily, is consistent with the causality requirement $`\beta _F1`$. In fact, the full solution (29) violates this condition for $`\mu \mu _cO(\mu _c)`$, although a precise statement is impossible without knowledge of $`k_{F_c}`$ to $`O(N_f^1)`$. This suggests that the HDL approximation leading to (11) is insufficient at such small $`\mu `$.

### 2.3 Friedel Oscillations

Another characteristic of a Fermi surface can be exposed by considering spatial correlations between $`q\overline{q}`$ pairs. A convenient way to do this is to consider correlation functions of the form

$`C(\stackrel{}{y};x_0)`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{x}}{}}\overline{\psi }\mathrm{\Gamma }\psi (\stackrel{}{0},0)\overline{\psi }\mathrm{\Gamma }\psi (\stackrel{}{x}+\stackrel{}{y},x_0)`$ (31)
$`=`$ $`\text{tr}{\displaystyle _p}{\displaystyle _q}\mathrm{\Gamma }{\displaystyle \frac{e^{ipx}}{ip/+\mu \gamma _0+M}}\mathrm{\Gamma }{\displaystyle \frac{e^{iqx}e^{i\stackrel{}{q}.\stackrel{}{y}}}{iq/+\mu \gamma _0+M}}`$ (32)
where in $`2+1`$ dimensions $`_p=d^3p/(2\pi )^3`$ and the matrix $`\mathrm{\Gamma }`$ determines the quantum numbers of the channel. $`M`$ is the fermion mass, which we will take to have the leading order form $`M(\mu )=\mathrm{\Sigma }_0\theta (\mu _c\mu )`$. Performing the Dirac trace and the sum over $`\stackrel{}{x}`$ yields

$$C(\stackrel{}{y};x_0)=\frac{4}{(2\pi )^4}𝑑q_0𝑑p_0d^2\stackrel{}{p}\frac{(\pm (p_0i\mu )(q_0i\mu )\pm \stackrel{}{p}^{\mathrm{\hspace{0.17em}2}}+M^2)e^{i(p_0q_0)x_0}e^{i\stackrel{}{p}.\stackrel{}{y}}}{[(p_0i\mu )^2+\stackrel{}{p}^{\mathrm{\hspace{0.17em}2}}+M^2][(q_0i\mu )^2+\stackrel{}{p}^{\mathrm{\hspace{0.17em}2}}+M^2]}$$
(33)
where the $`\pm `$ signs correspond to pseudoscalar and scalar channels respectively.

First consider the chirally broken phase $`\mu <\mu _c=\mathrm{\Sigma }_0`$; the $`p_0`$ ($`q_0`$) integral may be performed using Cauchy’s theorem by completing the contour in the upper (lower) half-plane, picking up residues at the poles at $`i\mu \pm i\sqrt{\stackrel{}{p}^{\mathrm{\hspace{0.17em}2}}+M^2}`$. Therefore

$`C(\stackrel{}{y};x_0)`$ $`=`$ $`{\displaystyle \frac{2}{(2\pi )^2}}{\displaystyle d^2\stackrel{}{p}e^{2x_0\sqrt{\stackrel{}{p}^{\mathrm{\hspace{0.17em}2}}+M^2}}e^{i\stackrel{}{p}.\stackrel{}{y}}\times \{\begin{array}{cc}1,\hfill & \text{pseudoscalar}\hfill \\ \frac{\stackrel{}{p}^{\mathrm{\hspace{0.17em}2}}}{\stackrel{}{p}^{\mathrm{\hspace{0.17em}2}}+M^2},\hfill & \text{scalar}\hfill \end{array}}`$ (34)
$`=`$ $`2{\displaystyle _0^{\mathrm{}}}p𝑑pJ_0(py)e^{2x_0\sqrt{p^2+M^2}}.`$ (35)
where $`J`$ is the Bessel function and in the second line we have specialised to the pseudoscalar case. For large $`x_0`$, the final integral over $`p`$ can now be approximated using Laplace’s method of asymptotic expansion, by expanding the integrand in powers of $`p/M`$, yielding

$$\underset{x_0\mathrm{}}{lim}C(\stackrel{}{y};x_0)\frac{M}{x_0}e^{2Mx_0}\mathrm{exp}\left(\frac{|\stackrel{}{y}|^2M}{4x_0}\right).$$
(36)
The general profile of the wavefunction is a gaussian with width increasing as $`\sqrt{x_0}`$.

In the chirally restored phase $`\mu >\mu _c`$, $`M`$ vanishes. For $`|\stackrel{}{p}|<\mu `$, both poles of the integrand of (33) lie in the upper half-plane, making the $`q_0`$ integral vanish. The effect is to raise the lower limit of the integral of (35) from 0 to $`\mu `$. The asymptotic expansion is now made by expanding the integrand in powers of $`(\frac{p}{\mu }1)`$, with the resulting form

$$\underset{x_00}{lim}C(\stackrel{}{y};x_0)\frac{\mu }{x_0}e^{2\mu x_0}J_0(\mu |\stackrel{}{y}|).$$
(37)
The wavefunction profile no longer changes with $`x_0`$, but instead oscillates with a spatial frequency determined by $`\mu `$, which to this order may be identified with the Fermi momentum $`k_F`$.

The oscillations observed in $`C(\stackrel{}{y};x_0)`$ are characteristic of a sharp Fermi surface and are reminiscent of oscillations in either the density-density correlation function, or the screened inter-particle potential, in degenerate systems generically known as Friedel Oscillations . The resemblance is only qualitative, however, because while the large-$`y`$ form of (37) is $`Cy^{\frac{1}{2}}\mathrm{sin}(k_Fy\frac{\pi }{4})`$, the form of eg. the screened potential in high-density QED is $`y^3\mathrm{cos}(2k_Fy)`$, and even $`y^4\mathrm{sin}(2k_Fy)`$ in the relativistic limit $`m_e\mu `$ . The disparity is due in part to the lower dimensionality, and partly because the function $`C(\stackrel{}{y};x_0)`$ probes correlations between fermion field variables rather than scalar densities. In Appendix A we present a similar asymptotic analysis of the density-density correlator $`C_{nn}(\stackrel{}{y};x_0)=_\stackrel{}{x}j_0(\stackrel{}{0},0)j_0^{}(\stackrel{}{x},x_0)j_0(\stackrel{}{0},0)j_0^{}(\stackrel{}{x}+\stackrel{}{y},x_0)`$, where $`j_0(x)`$ is the local baryon density $`\overline{\psi }\gamma _0\psi (x)`$, and find a decay with $`y`$ more closely resembling the forms in the literature. In particular, because $`C_{nn}`$ describes correlations between particle-hole pairs at the Fermi surface, the asymptotic spatial frequency of the oscillations is $`2k_F`$ rather than $`k_F`$. Because of the similar origin, however, we will refer to oscillations of the form (37) as Friedel oscillations in what follows.

## 3 Numerical Results

In this section we present results from a numerical simulation of the Z2 GN model using staggered lattice fermions. The lattice formulation and simulation algorithm are described in detail in Refs. . We used a $`32^2\times 48`$ lattice and throughout chose a bare fermion mass $`m=0`$ and $`N=2`$ lattice flavors. Unless otherwise stated the coupling constant used was $`a/g^2=0.75`$, corresponding to a physical fermion mass at $`\mu =0`$ of $`\mathrm{\Sigma }_0a0.17`$. At each value of $`\mu `$ between 2000 and 4000 hybrid Monte Carlo trajectories of mean length 1.0 were taken. All results have been taken in the chirally-restored phase with $`\mu >\mu _c`$, whose value is estimated to be $`\mu _ca0.16(1)`$. Values for the baryon density per staggered fermion flavor $`na^2`$ for values of chemical potential $`\mu a[0.2,0.8]`$ are given in Table 1, showing that the density rises steadily as $`\mu `$ is increased. The reported densities are roughly twice the free-field value $`n=\mu ^2/2\pi `$, consistent with fermion doubling; $`N`$ flavors of staggered lattice fermion describe $`N_f=2N`$ continuum flavors in 2+1 dimensions, so that $`N_f=4`$ in our simulation. Henceforth we will quote results in units in which the lattice spacing $`a=1`$.

### 3.1 Fermion Dispersion Relation

The fermion timeslice propagator at fixed spatial momentum $`\stackrel{}{k}`$ is defined by $`𝒢(\stackrel{}{k},x_0)=_\stackrel{}{x}𝒢(\stackrel{}{0},0;\stackrel{}{x},x_0)e^{i\stackrel{}{k}.\stackrel{}{x}}`$, where $`𝒢(x,y)\psi (x)\overline{\psi }(y)`$. The decay in Euclidean time $`x_0`$ then yields the fermion dispersion relation $`E(\stackrel{}{k})`$. For a degenerate system with Fermi momentum $`k_F`$, excitations with $`|\stackrel{}{k}|<k_F`$ are hole-like and those with $`|\stackrel{}{k}|>k_F`$ particle-like; the lowest energy excitations are those in the neighbourhood of the Fermi surface with $`|\stackrel{}{k}|k_F`$. Fig. 3 shows $`E(\stackrel{}{k})`$ measured for $`\stackrel{}{k}`$ parallel to the $`x`$-axis for four values of $`\mu >\mu _c`$, using the procedure specified in . In order to obtain a smooth curve we plot $`E(|\stackrel{}{k}|)`$ as negative for $`|\stackrel{}{k}|<k_F`$. Note also that $`|\stackrel{}{k}|`$ is plotted as a fraction of $`\pi `$ – since the staggered fermion action is only invariant under translations of length $`2a`$ the maximum momentum which can be probed is $`\pi /2`$.

It is obvious that there is no mass gap at the Fermi surface characteristic of a BCS instability, since for certain combinations of $`\mu `$ and $`\stackrel{}{k}`$ our data have $`E(|\stackrel{}{k}|)0`$. Indeed, it is possible to fit a smooth function

$$E(|\stackrel{}{k}|)=E_0+D\mathrm{sinh}^1(\mathrm{sin}|\stackrel{}{k}|)$$
(38)
to the data, which corrects the Fermi liquid form (18) for lattice discretisation effects, and permits the extraction of an effective Fermi momentum $`K_F`$ and velocity $`\beta _F`$ given by

$$K_F=\mathrm{sinh}^1(\mathrm{sin}k_F)=\frac{E_0}{D};\beta _F=D\frac{\mathrm{cosh}E_0}{\mathrm{cosh}K_F}.$$
(39)
The fitted values of $`K_F,\beta _F`$ are given in Table 2. The departures from free-field values are small and of similar magnitude to the predictions (30) of the $`1/N_f`$ expansion $`K_F/\mu 0.953`$, $`\beta _F0.984`$. Also shown is the ratio $`K_F/\mu \beta _F`$, which apart from $`\mu =0.2`$ shows no significant deviation from its free-field value 1. This constrasts sharply with Monte Carlo studies of the 2+1$`d`$ NJL model, for which $`𝔤=4`$, where values in the range 1.3 - 1.5 were observed . This suggests that the quasiparticle spectrum observed in that model is determined by physics outside the scope of the $`1/N_f`$ expansion, and as such the current result provides indirect support for the exotic gapless superfluid scenario proposed in .

### 3.2 The $`\sigma `$ Dispersion Relation

Just as in the continuum $`1/N_f`$ expansion, the lattice formulation of the GN model encodes the interactions between the fermions via an auxiliary bosonic scalar field $`\sigma `$ whose value is given by the equation of motion $`\sigma =\frac{g}{\sqrt{N_f}}\overline{\psi }\psi `$ (in this paper we define $`\sigma `$ with dimension $`\frac{3}{2}`$). For $`\mu =0`$ and $`k\mathrm{\Sigma }_0`$ the auxiliary propagator (2) resembles that of a standard boson of mass $`2\mathrm{\Sigma }_0`$. Since this mass coincides with the threshold for an $`q\overline{q}`$ continuum, $`D_\sigma `$ is dominated by a branch cut rather than an isolated pole in the complex $`k`$-plane making standard lattice spectroscopy techniques difficult to apply . A recent study using the Maximum Entropy Method has shown evidence for a non-zero binding energy in the $`\sigma `$ channel for finite $`N_f`$ .

For $`\mu >\mu _c`$, however, the calculation of Sec. 2.1 suggests that the $`\sigma `$ is more tightly bound. We have studied this by measuring the auxiliary timeslice propagator $`𝒟(\stackrel{}{k},|x_0|)=_{\stackrel{}{x},x_0=\pm |x_0|}\sigma (\stackrel{}{0},0)\sigma (\stackrel{}{x},x_0)e^{i\stackrel{}{k}.\stackrel{}{x}}`$, once again choosing $`\stackrel{}{k}`$ along the $`x`$-axis222Note that as discussed below eqn. (14) $`𝒟`$ may not in fact be time-reversal symmetric for $`\stackrel{}{k}\stackrel{}{0}`$.. Fig. 4 plots $`𝒟(x_0)`$ for $`\stackrel{}{k}=\stackrel{}{0}`$, showing that with the exception of $`\mu =0.2`$, simple pole fits of the form $`𝒟(x_0)e^{M_\sigma x_0}`$ describe the data reasonably well before it descends into noise. Values of $`M_\sigma `$ extracted from fits to timeslices $`x_0[1,5]`$ are plotted in Fig. 5. Also shown is the leading order prediction $`M_\sigma =2\sqrt{\mu (\mu \mu _c)}`$, with a value 0.16(1) assumed for $`\mu _c`$. The numerical results fall 20-30% lower than the leading order prediction, which is consistent with a correction of $`O(N_f^1)`$.

Results from simple pole fits to $`𝒟(\stackrel{}{k},|x_0|)`$ for $`|\stackrel{}{k}|\frac{\pi }{2}`$ yielding the $`\sigma `$ dispersion relation $`E(k)`$ are plotted in Fig. 6 333Since under the global Z2 $`\sigma (x)(1)^{x_0+x_1+x_2}\sigma (x)`$, in the symmetric phase the propagator obeys $`𝒟(k_i)=𝒟(\pi k_i)`$.. Only for the lowest value $`\mu =0.2`$ is there any marked change in $`E`$ as $`k`$ is increased, in which case $`Ek`$. For larger $`\mu `$ the dispersion is flatter. There is qualitative agreement with the leading order $`1/N_f`$ result shown in Fig. 1, although an improved understanding of discretisation effects close to $`|\stackrel{}{k}|=\frac{\pi }{2}`$ would be needed to make the comparison more quantitative.

Nonetheless, Figs. 4 and 6 demonstrate unambiguously that in-medium effects are observed in the lattice simulation.

### 3.3 Mesonic Dispersion Relations and Zero Sound

Next we investigate mesonic (ie. $`q\overline{q}`$) states. We have studied mesonic correlation functions $`𝒞_\mathrm{\Gamma }(\stackrel{}{k},x_0)=_\stackrel{}{x}j_\mathrm{\Gamma }(\stackrel{}{0},0)j_\mathrm{\Gamma }^{}(\stackrel{}{x},x_0)e^{i\stackrel{}{k}.\stackrel{}{x}}`$ where the bilinears $`j_\mathrm{\Gamma }(x)`$ are defined with scalar, pseudoscalar or vector quantum numbers; in terms of staggered fermion fields $`\chi ,\overline{\chi }`$ the operators are

$$j_{11}(x)=\overline{\chi }_x\chi _x;j_{\gamma _5}(x)=\epsilon _x\overline{\chi }_x\chi _x;j_{\gamma _i}(x)=\frac{\eta _{ix}}{2}[\overline{\chi }_x\chi _{x+\widehat{ı}}+\overline{\chi }_{x+\widehat{ı}}\chi _x],$$
(40)
where $`\eta _{1x}=(1)^{x_0}`$, $`\eta _{2x}=(1)^{x_0+x_1}`$ and $`\epsilon _x=(1)^{x_0+x_1+x_2}`$. As before we study the timeslice correlators in each channel as a function of spatial momentum $`\stackrel{}{k}\widehat{x}`$.

As expected for $`\mu >\mu _c`$ where the Z2 global symmetry is restored, the results in scalar and pseudoscalar channels coincide. The results on odd timeslices, averaged over $`k`$ and $`\pi k`$, are shown for four different $`\mu >\mu _c`$ in Fig. 7. Note the difference in the vertical logarithmic scales between the upper two panels $`\mu =0.2,0.4`$ and the lower two $`\mu =0.6,0.8`$. For each $`\mu `$, $`𝒞_{\gamma _5}(\stackrel{}{k}=\stackrel{}{0})`$ is almost flat, indicating a massless state.

Appendix B sketches the calculation of $`𝒞_\mathrm{\Gamma }(\stackrel{}{k},x_0)`$ as a function of $`\mu `$ in free field theory. As $`x_0\mathrm{}`$ the function is dominated by a continuum of particle-hole pairs at or near the Fermi surface, which effectively cost zero energy to excite. The generic result is that for $`|\stackrel{}{k}|2\mu `$ the decay is algebraic, with (59,66)

$$𝒞_\mathrm{\Gamma }(\stackrel{}{k},x_0)\{\begin{array}{cc}x_0^2,\hfill & |\stackrel{}{k}|\mu \text{;}\hfill \\ x_0^{\frac{3}{2}},\hfill & |\stackrel{}{k}|2\mu \text{.}\hfill \end{array}$$
(41)
Only once $`|\stackrel{}{k}|>2\mu `$ does it become kinematically impossible to excite a pair with zero energy, resulting in exponential decay:

$$𝒞_\mathrm{\Gamma }(\stackrel{}{k},x_0)x_0^{\frac{3}{2}}\mathrm{exp}\left((|\stackrel{}{k}|2\mu )x_0\right).$$
(42)
The sequence of plots in Fig. 7 is in qualitative agreement with these findings. For each $`\mu `$ there is a particular value of $`|\stackrel{}{k}|`$, highlighted with a solid line in the plots, for which the temporal falloff is particularly slow, corresponding to $`|\stackrel{}{k}|2\mu `$ (eg. for $`\mu =0.4`$ the slow falloff occurs for $`|\stackrel{}{k}|=\frac{\pi }{4}0.785`$). For $`|\stackrel{}{k}|`$ larger than this value the decay is much steeper, although only for $`\mu =0.2`$ does it resemble an exponential form. Because of the technical difficulties in treating correlators with power-law decays in a finite volume (see eg. ) we have made no attempt to fit the numerical data for $`𝒞_{\gamma _5}(|\stackrel{}{k}|,x_0)`$ to a functional form. It can be observed, however, that the overall magnitude of $`𝒞_{\gamma _5}`$ increases with $`\mu `$ as the size of the Fermi surface and hence the number of participating particle-hole states grows. There is, however, no obvious correspondence with the prefactor $`|\stackrel{}{k}|/\mu `$ predicted in (59), probably because the approximation $`|\stackrel{}{k}|\mu `$ is not valid over the accessible momentum range.

The situation is more interesting in the vector channel, corresponding to the quantum numbers of the $`\rho `$ meson, because for $`|\stackrel{}{k}|>0`$ we can distinguish between $`𝒞_\gamma _{}(\stackrel{}{k})`$, in which the component of the vector is parallel to $`\stackrel{}{k}`$, and $`𝒞_\gamma _{}(\stackrel{}{k})`$. In Fig. 8 we plot the correlator for several $`|\stackrel{}{k}|`$ values at $`\mu =0.6`$ in each case. For $`𝒞_\gamma _{}`$, just as for $`𝒞_{\gamma _5}`$, there is no difference between the data for $`|\stackrel{}{k}|`$ and $`\pi |\stackrel{}{k}|`$. For $`𝒞_\gamma _{}`$, however, because of the point-split nature of the operator $`j_{\gamma _i}`$ this symmetry no longer holds; in this case we choose to plot data with $`|\stackrel{}{k}|[\frac{\pi }{2},\pi ]`$ (the data for $`|\stackrel{}{k}|\frac{\pi }{2}`$ are very similar to those for $`𝒞_\gamma _{}`$). Each plot covers two decades on the vertical axis; however both the magnitude and the shape of the curves is very different. For $`𝒞_\gamma _{}`$ the curves are qualitatively very similar to those of Fig. 7 at $`\mu =0.6`$, with a distinguished momentum $`|\stackrel{}{k}|=\frac{3\pi }{8}`$. The correlator $`𝒞_\gamma _{}`$ is much smaller in magnitude, and is consistent with exponential rather than algebraic decay. The lines are fits of the form $`𝒞_\gamma _{}(|\stackrel{}{k}|,x_0)=A(e^{Ex_0}+e^{E(L_tx_0)})`$. The resulting $`E(|\stackrel{}{k}|)`$ is shown in Fig. 9. For small $`|\stackrel{}{k}|`$ it resembles that of a massless pole, ie. $`E=\beta _0|\stackrel{}{k}|`$, with velocity $`\beta _00.5`$. It should be mentioned, however, that although $`𝒞_\gamma _{}`$ and $`𝒞_\gamma _{}`$ still differed, no evidence for a massless pole was seen in the data at $`\mu =0.8`$.

Light states in the $`\rho `$ channel in medium have long been of interest . One recently proposed theoretical scenario is that the $`\rho `$ should be viewed as an “almost” Goldstone boson due to a spontaneous “induced” breaking of Lorentz symmetry over and above the explicit breaking due to $`\mu 0`$ . In this case the $`\rho `$ mass is predicted to scale as $`m_\rho ^22\mu n`$. Another, proposed in the context of Two Color QCD (TCQCD) in which mesons and baryons fall in the same multiplets and are degenerate at $`\mu =0`$, is that for sufficiently large $`\mu `$ there is a condensation of vector diquark states leading to the spontaneous breaking of rotational invariance and vector Goldstone states with both $`E|\stackrel{}{k}|`$ and $`E|\stackrel{}{k}|^2`$ . In-medium modification of the $`\rho `$ correlator has recently been observed in lattice Monte Carlo simulations of TCQCD .

Since our results suggest a massless state and there is no Fermi surface in TCQCD we favour an alternative explanation, namely that the pole signals a collective excitation of the system in which the shape of the Fermi surface is distorted, leading to a rotationally non-invariant disturbance propagating at velocity $`\beta _0`$ known as zero sound. The conditions for zero sound to propagate are that there must be interactions, ie. $`_{\stackrel{}{k},\stackrel{}{k}^{}}0`$, and that inter-particle collisons are negligible, ensured if $`T0`$ . Indeed, zero sound solutions can in principle be found self-consistently using the theoretical framework of Sec. 2.2, although the specific form (22) for the Fermi liquid interaction entails the solution of an integral equation of the second kind, beyond the scope of this paper. A potential problem for this picture is that the extracted velocity $`\beta _00.5`$ is less than $`\beta _F1`$, implying that quasiparticles would experience damping via Čerenkov radiation of zero sound. Clearly further work exploring the systematic effects of varying $`\mu `$, $`\mathrm{\Sigma }_0`$, $`\stackrel{}{k}`$ and volume will be needed for a more complete understanding to emerge.

### 3.4 Mesonic Wavefunctions and Friedel Oscillations

We now switch attention to spatial correlations between particles, probed via the wavefunction $`\mathrm{\Psi }(\stackrel{}{x})`$ defined by

$$\mathrm{\Psi }_\mathrm{\Gamma }(\stackrel{}{x})=\underset{x_0\mathrm{}}{lim}\mathrm{\Psi }_\mathrm{\Gamma }^1(\stackrel{}{x}=\stackrel{}{0})\underset{\stackrel{}{y}}{}𝒢_q(\stackrel{}{0},0;\stackrel{}{y},x_0)\mathrm{\Gamma }𝒢_{\overline{q}}(\stackrel{}{0},0;\stackrel{}{y}+\stackrel{}{x},x_0)\mathrm{\Gamma }$$
(43)
where as ususal $`\mathrm{\Gamma }`$ projects out the quantum numbers of the channel of interest. In the GN and related models $`\mathrm{\Psi }(\stackrel{}{x})`$ is technically much easier to define and measure than in QCD-like theories where they are gauge-dependent. Meson wavefunctions in the 2+1$`d`$ GN model have been studied at $`T,\mu =0`$ in , where further technical details are given.

Fig. 10 shows $`\mathrm{\Psi }(\stackrel{}{x})`$ measured along the $`x`$-axis in the scalar channel at a constant $`\mu =0.5`$ at various values of the the coupling $`1/g^2=0.65,0.70,0.75`$. The oscillatory behaviour described in Sec. 2.3 is clearly seen, and is plainly not a discretisation artifact since its form is stable as $`a(g)`$ varies. Fig. 11 shows that there are no significant differences among the various mesonic channels implying that in contrast to the situation at $`\mu =0`$ , effects due to eg. $`\sigma `$ exchange are very hard to detect.

In the next four figures we plot the scalar wavefunctions at $`1/g^2=0.75`$ for a sequence of $`\mu `$ values 0.2, 0.4, 0.6, and 0.8. The simulation data are connected by solid lines. As $`\mu `$, and hence $`k_F`$, increases the oscillations decrease in wavelength in accordance with theoretical expectation, and provide a graphic illustration of the presence of a sharp Fermi surface. The dashed lines show measurements taken with the same lattice parameters but with the interaction switched off. The disparity with the interacting theory is small, though increasing with $`|\stackrel{}{x}|`$, showing that the free field description of the oscillations is qualitatively correct. A possible explanation is the infinite value predicted for the Debye mass in Sec. 2.1, implying that interactions between static quarks at non-zero $`|\stackrel{}{x}|`$ are completely screened. Also shown is the theoretical form (37) $`\mathrm{\Psi }(x)=J_0(k_Fx)`$ with $`k_F\mu `$, showing good agreement with the data for small $`x\stackrel{<}{}2k_F^1`$. Unfortunately it appears hard to obtain more quantitative information, such as an independent fit for the Fermi momentum $`k_F`$, because $`|J_0(kx)|`$ decays only as $`x^{\frac{1}{2}}`$. This means that fits should not only include the backwards-propagating signal $`J_0(k_F(L_sx))`$ but also image contributions $`J_0(k_F(nL_sx))`$  – our attempts to find a satisfactory fit were unsuccessful. The figures therefore simply show both the “forwards” $`J_0(\mu x)`$ and “forwards-and-backwards” $`J_0(\mu x)+J_0(\mu (L_sx))`$ forms, showing that neither gives a satisfactory description of the data over the full range of $`x`$ and $`\mu `$. It is also possible that more sophisticated fits taking proper account of discretisation effects are needed; at $`\mu =0.8`$ the baryon density $`n=0.31`$ (see Table 1), a significant fraction of its saturation value of 1, so that the Fermi surface may be distorted.

## 4 Summary

The GN and related models remain the only interacting field theories both simulable by standard lattice methods at $`\mu 0`$ and displaying a Fermi surface, thought to be an essential feature of dense quark matter. In this paper we have attempted to develop an understanding of how orthodox lattice observables are altered in such conditions, and how they should be interpreted using the language of many-body physics. The main achievements have been:

* An analytic calculation of the auxiliary boson propagator $`D_\sigma `$ in medium to leading order in $`N_f^1`$ in the HDL approximation. The branch cut in the complex-$`k`$ plane is modified to become an isolated pole. Physically, the Debye mass is infinite, but the plasma frequency finite and vanishing as onset is approached. The dispersion relation $`E(k)`$ is non-trivial. Lattice simulation has verified these results with acceptable precision (ie. to within $`O(N_f^1`$)), and have unambiguously identified in-medium effects in this channel.
* A systematic calculation of Fermi liquid parameters up to $`O(N_f^1)`$, which is consistent with causality constraints and has in turn been verified by lattice simulation. This contrasts with existing studies of the quasiparticle dispersion relation in the 2+1$`d`$ NJL model with $`\mu >\mu _c`$, which are not consistent with the large-$`N_f`$ predictions.
* A study of meson correlation functions at $`\stackrel{}{k}0`$; instead of showing exponential falloff with Euclidean time these generically decay algebraically, signalling the presence of massless particle-hole excitations for all $`\stackrel{}{k}`$. There is qualitative agreement between the results of free field theory and lattice simulation, which both show a non-trivial dependence on $`|\stackrel{}{k}|/\mu `$. There is also tentative evidence for a massless pole in the vector channel, which is possibly a manifestation of zero sound.
* The first observation of oscillatory behaviour in mesonic wavefunctions, which resemble the Friedel oscillations familiar in many-body physics. Figs. 12-15 offer a graphic confirmation, if one is still needed, of the presence of a sharp Fermi surface in this model.

It would be interesting to apply the same analysis to the lattice NJL model at $`\mu 0`$, which has been studied in both 2+1$`d`$  and 3+1$`d`$ . Independent of the issue of whether the ground state has a non-vanishing diquark condensate, these models are more realistic because they contain physical pions, whose behaviour at the quark-hadron transition is thought to change dramatically. Indeed, the result of Sec. 2.2 that the plasma frequency (pole mass) in the chirally symmetric phase decreases as $`\mu \mu _{c+}`$ while the Debye (screening) mass diverges is reminiscent of recent results obtained by Son and Stephanov in the chirally broken phase  using an effective theory and scaling arguments. It would be straightforward to repeat our analysis for the pion dispersion relation. More generally, it is worth noting that should an algorithm effective at $`\mu 0`$ ever emerge our results in the meson sector are in principle accessible as gauge invariant correlators in a QCD simulation. The $`\sigma `$ channel in the current model is dominated by disconnected $`q\overline{q}`$ bubbles which are technically hard to compute if an auxiliary field is absent, but in the lightest channel in QCD, namely the pion, the relevant diagrams have connected quark lines. The main difference in the baryon sector, of course, is that gauge invariant quasiparticles now have $`qqq`$ quantum numbers, and at least below the transition should be non-relativistic at the Fermi surface.

## Acknowledgements

SJH was supported by a PPARC Senior Research Fellowship, CGS by the Leverhulme Trust and DOE grant DE-FG02-96ER40945, and JBK by NSF grant NSF-PHY-0102409. We are grateful to Kurt Langfeld for discussions in the early stages of this work.

## Appendix A

In this appendix we outline how oscillatory behaviour as a function of spatial separation $`\stackrel{}{y}`$ resembling Friedel oscillations appears in the large Euclidean time $`x_0\mathrm{}`$ limit of the density-density correlator $`C_{nn}(\stackrel{}{y};x_0)=_\stackrel{}{x}j_0(\stackrel{}{0},0)j_0^{}(\stackrel{}{x},x_0)j_0(\stackrel{}{0},0)j_0^{}(\stackrel{}{x}+\stackrel{}{y},x_0)`$, where $`j_0(x)=\overline{\psi }\gamma _0\psi (x)`$. The relevant expression is

$$\underset{\stackrel{}{x}}{}\text{tr}_p_q_k_{\mathrm{}}\left\{\gamma _0\frac{e^{ipx}}{ip/+\mu \gamma _0}\gamma _0\frac{e^{iqx}}{iq/+\mu \gamma _0}\gamma _0\frac{e^{ikx}e^{i\stackrel{}{k}.\stackrel{}{y}}}{ik/+\mu \gamma _0}\gamma _0\frac{e^{i\mathrm{}x}e^{i\stackrel{}{\mathrm{}}.\stackrel{}{y}}}{i\mathrm{}/+\mu \gamma _0}\right\},$$
(44)
where we have assumed massless quark propagation in the high density phase. The trace yields three terms of which the first is

$$_p_q\frac{[\stackrel{}{p}.\stackrel{}{q}(p_0i\mu )(q_0i\mu )]e^{i(pq)x}}{[(p_0i\mu )^2+\stackrel{}{p}^2][(q_0i\mu )^2+\stackrel{}{q}^2]}\times _k_{\mathrm{}}\frac{[\stackrel{}{k}.\stackrel{}{\mathrm{}}(k_0i\mu )(\mathrm{}_0i\mu )]e^{i(k\mathrm{})x}e^{i(\stackrel{}{k}\stackrel{}{\mathrm{}}).\stackrel{}{y}}}{[(k_0i\mu )^2+\stackrel{}{k}^2][(\mathrm{}_0i\mu )^2+\stackrel{}{\mathrm{}}^2]}$$
(45)
The integrals over timelike momenta are performed by completing the contour in the upper half-plane for $`𝑑p_0`$, $`𝑑k_0`$ and in the lower half-plane for $`𝑑q_0`$, $`𝑑\mathrm{}_0`$, noting the pole positions at eg. $`p_0=i\mu \pm i|\stackrel{}{p}|`$. In fact, only the pole with negative imaginary part contributes to $`𝑑q_0`$ and $`𝑑\mathrm{}_0`$, and it is straightforward to show that the dominant term for large $`x_0`$ in the resulting expression comes from the poles at $`p_0=i\mu i|\stackrel{}{p}|`$ and $`k_0=i\mu i|\stackrel{}{k}|`$ since all others are suppressed by powers of $`e^{\mu x_0}`$. The resulting contribution from all three trace terms is

$`{\displaystyle \underset{\stackrel{}{x}}{}}{\displaystyle \frac{1}{4}}{\displaystyle _0^\mu }{\displaystyle \frac{d^dp}{(2\pi )^d}}{\displaystyle _\mu ^{\mathrm{}}}{\displaystyle \frac{d^dq}{(2\pi )^d}}{\displaystyle _0^\mu }{\displaystyle \frac{d^dk}{(2\pi )^d}}{\displaystyle _\mu ^{\mathrm{}}}{\displaystyle \frac{d^d\mathrm{}}{(2\pi )^d}}e^{i(\stackrel{}{p}\stackrel{}{q}+\stackrel{}{k}\stackrel{}{\mathrm{}}).\stackrel{}{x}}e^{(|\stackrel{}{q}||\stackrel{}{p}|+|\stackrel{}{\mathrm{}}||\stackrel{}{k}|)x_0}e^{i(\stackrel{}{k}\stackrel{}{\mathrm{}}).\stackrel{}{y}}\times `$
$`\left\{(1+\mathrm{cos}\theta _{pq})(1+\mathrm{cos}\theta _k\mathrm{})(1\mathrm{cos}\theta _{pk})(1\mathrm{cos}\theta _q\mathrm{})+(1+\mathrm{cos}\theta _p\mathrm{})(1+\mathrm{cos}\theta _{qk})\right\},`$ (46)
where $`\theta _{pq}`$ denotes the angle between vectors $`\stackrel{}{p}`$ and $`\stackrel{}{q}`$, etc.

Now, the sum over the timeslice $`_\stackrel{}{x}`$ results in an overall momentum-conserving delta-function. Since, however, the integrand is dominated by momenta in the immediate vicinity of the Fermi surface, a physically reasonable approximation, which preserves all angles under the remaining $`d^dpd^dq`$, is given by

$$(2\pi )^d\delta ^d(\stackrel{}{p}\stackrel{}{q}+\stackrel{}{k}\stackrel{}{\mathrm{}})(2\pi )^d\delta ^d(\stackrel{}{p}+\stackrel{}{k})(2\pi )^d\delta ^d(\stackrel{}{q}+\stackrel{}{\mathrm{}}).$$
(47)
The expression (46) reduces to

$$_0^\mu \frac{d^dp}{(2\pi )^d}_\mu ^{\mathrm{}}\frac{d^dq}{(2\pi )^d}e^{2(|\stackrel{}{q}||\stackrel{}{p}|)x_0}e^{i(\stackrel{}{p}\stackrel{}{q}).\stackrel{}{y}}\frac{1}{2}(\mathrm{cos}^2\theta _{pq}1).$$
(48)
Denote the angle between $`\stackrel{}{p}`$ ($`\stackrel{}{q}`$) and $`\stackrel{}{y}`$ by $`\theta _1`$ ($`\theta _2`$). Then

$$\mathrm{cos}\theta _{pq}=\mathrm{cos}\theta _1\mathrm{cos}\theta _2+\mathrm{sin}\theta _1\mathrm{sin}\theta _2\mathrm{cos}\mathrm{\Theta }_{pq},$$
(49)
where $`\mathrm{\Theta }_{pq}`$ is the angle between the vectors projected onto the subspace $`\mathrm{\Omega }^{d1}`$ orthogonal to $`\stackrel{}{y}`$. The angular integrals over this subspace may be performed using

$$𝑑\mathrm{\Omega }^{d1}=\frac{2\pi ^{(d1)/2}}{\mathrm{\Gamma }(\frac{d1}{2})};𝑑\mathrm{\Omega }^{d1}\mathrm{cos}^2\mathrm{\Theta }=\frac{1}{d1}\frac{2\pi ^{(d1)/2}}{\mathrm{\Gamma }(\frac{d1}{2})},$$

to yield

$`{\displaystyle \frac{1}{(2\pi )^{2d}}}{\displaystyle \frac{2\pi ^{d1}}{\mathrm{\Gamma }^2(\frac{d1}{2})}}{\displaystyle _0^\mu }p^{d1}dp{\displaystyle _\mu ^{\mathrm{}}}q^{d1}dqe^{2(qp)x_0}{\displaystyle _0^\pi }d\theta _1{\displaystyle _0^\pi }d\theta _2e^{ipy\mathrm{cos}\theta _1}e^{iqy\mathrm{cos}\theta _2}\times `$
$`\mathrm{sin}^{d2}\theta _1\mathrm{sin}^{d2}\theta _2\left\{\mathrm{cos}^2\theta _1\mathrm{cos}^2\theta _2+{\displaystyle \frac{1}{d1}}\mathrm{sin}^2\theta _1\mathrm{sin}^2\theta _21\right\}.`$ (50)
The remaining angular integrals are readily performed to yield

$`{\displaystyle \frac{(d1)}{2(2\pi )^d}}`$ $`{\displaystyle _0^\mu }p^{d1}dp{\displaystyle _\mu ^{\mathrm{}}}q^{d1}dqe^{2(qp)x_0}\times `$ (51)
$`\left({\displaystyle \frac{1}{\sqrt{pq}y}}\right)^{d2}`$ $`\left\{{\displaystyle \frac{1}{py}}J_{\frac{d}{2}}(py)J_{\frac{d2}{2}}(qy){\displaystyle \frac{1}{qy}}J_{\frac{d}{2}}(qy)J_{\frac{d2}{2}}(py)+{\displaystyle \frac{d}{pqy^2}}J_{\frac{d}{2}}(py)J_{\frac{d}{2}}(qy)\right\}.`$
The integrals over vector lengths $`p`$ and $`q`$ are estimated using Laplace’s method;

$$\underset{x_0\mathrm{}}{lim}_0^\mu 𝑑pp^{d1}pe^{2px_0}f(py)_\mu ^{\mathrm{}}𝑑qq^{d1}e^{2qx_0}g(qy)\frac{\mu ^{2d2}}{4x_0^2}f(\mu y)g(\mu y)\left(1+O(x_0^1)\right),$$
(52)
to yield the final result

$$\frac{(d1)}{8(2\pi )^d}\frac{\mu ^{2d2}}{x_0^2(\mu y)^{d2}}\left[\frac{2}{\mu y}J_{\frac{d2}{2}}(\mu y)J_{\frac{d}{2}}(\mu y)+\frac{d}{(\mu y)^2}J_{\frac{d}{2}}^2(\mu y)\right].$$
(53)
One may now use the asymptotic form of the Bessel function to deduce the long-distance behaviour of $`C_{nn}(|\stackrel{}{y}|;x_0)`$:

$$C_{nn}(y;x_0)\{\begin{array}{cc}\frac{1}{16\pi ^3x_0^2}\frac{\mathrm{cos}2\mu y}{y^2},\hfill & d=2\text{;}\hfill \\ \frac{\mu }{16\pi ^4x_0^2}\frac{\mathrm{sin}2\mu y}{y^3},\hfill & d=3\text{.}\hfill \end{array}$$
(54)
## Appendix B

In this section we analyse the behaviour of mesonic correlators as a function of spatial momentum Euclidean time in the high density phase. For simplicity we fix $`d=2`$. The general expression is

$$𝒞_\mathrm{\Gamma }(\stackrel{}{k};x_0)=\underset{\stackrel{}{x}}{}_p_q\frac{\text{tr}\{\mathrm{\Gamma }(i\stackrel{}{p}.\stackrel{}{\gamma }i(p_0i\mu )\gamma _0)\mathrm{\Gamma }(i\stackrel{}{q}.\stackrel{}{\gamma }i(q_0i\mu )\gamma _0)\}}{[(p_0i\mu )^2+|\stackrel{}{p}|^2][(q_0i\mu )^2+|\stackrel{}{q}|^2]}e^{i(pq)x}e^{i\stackrel{}{k}.\stackrel{}{x}},$$
(55)
where we will consider $`\mathrm{\Gamma }\{11,\gamma _5,\gamma _{},\gamma _{},\gamma _o\}`$, $``$ and $``$ denoting components parallel and perpendicular to $`\stackrel{}{k}`$. On performing the trace, the timeslice sum $`_\stackrel{}{x}`$ and the timelike momentum integrals, we are left with, eg.

$$𝒞_{\gamma _5}(\stackrel{}{k};x_0)=\frac{d^2\stackrel{}{p}}{(2\pi )^2}\theta (\mu |\stackrel{}{p}|)\theta (|\stackrel{}{p}+\stackrel{}{k}|\mu )\frac{|\stackrel{}{p}||\stackrel{}{p}+\stackrel{}{k}|\stackrel{}{p}.(\stackrel{}{p}+\stackrel{}{k})}{|\stackrel{}{p}||\stackrel{}{p}+\stackrel{}{k}|}e^{(|\stackrel{}{p}+\stackrel{}{k}||\stackrel{}{p}|)x_0}$$
(56)
where we have assumed $`x_0>0`$ and once again focussed on the contribution from the dominant pole $`p_0=i(\mu |\stackrel{}{p}|)`$. We can evaluate the remaining integral in two limits:

$`\left|\stackrel{}{k}\right|\mu `$:

In this case the integration region is a narrow crescent in the vicinity of the Fermi surface, and the integrand and limits can be expanded in powers of $`|\stackrel{}{k}|/|\stackrel{}{p}|`$. The general form is then

$$𝒞_\mathrm{\Gamma }(\stackrel{}{k};x_0)=2\mu |\stackrel{}{k}|_{\frac{\pi }{2}}^{\frac{\pi }{2}}\frac{d\theta }{2\pi }h_\mathrm{\Gamma }(\theta )\mathrm{exp}(|\stackrel{}{k}|x_0\mathrm{cos}\theta )\left(1+O\left(\frac{|\stackrel{}{k}|}{\mu }\right)\right)$$
(57)
where the leading contribution for each channel is given by

$`h_{\gamma _5}=h_{11}={\displaystyle \frac{|\stackrel{}{k}|^2}{4\mu ^2}}\mathrm{cos}\theta \mathrm{sin}^2\theta ;h_\gamma _{}=\mathrm{cos}\theta \mathrm{sin}^2\theta ;h_{\gamma _0}=\mathrm{cos}\theta ;`$
$`h_\gamma _{}=\mathrm{cos}^3\theta +{\displaystyle \frac{3|\stackrel{}{k}|}{2\mu }}\mathrm{cos}^2\theta \mathrm{sin}^2\theta +{\displaystyle \frac{|\stackrel{}{k}|^2}{4\mu ^2}}\mathrm{cos}\theta \mathrm{sin}^2\theta (1+2\mathrm{sin}^2\theta ).`$ (58)
The angular integral leads to the following results for the leading behaviour in $`|\stackrel{}{k}|/\mu `$:

$`𝒞_{\gamma _5,11}(\stackrel{}{k};x_0)`$ $`=`$ $`\pm {\displaystyle \frac{|\stackrel{}{k}|^3}{2\pi \mu }}\left[{\displaystyle \frac{1}{3}}+{\displaystyle \frac{\pi }{2}}\left({\displaystyle \frac{𝐋_2(|\stackrel{}{k}|x_0)I_2(|\stackrel{}{k}|x_0)}{|\stackrel{}{k}|x_0}}\right)\right]\pm {\displaystyle \frac{|\stackrel{}{k}|}{2\pi \mu x_0^2}}+O(x_0^4)`$ (59)
$`𝒞_\gamma _{}(\stackrel{}{k};x_0)`$ $`=`$ $`{\displaystyle \frac{4\mu ^2}{|\stackrel{}{k}|^2}}𝒞_{\gamma _5}(\stackrel{}{k};x_0){\displaystyle \frac{2\mu }{\pi |\stackrel{}{k}|x_0^2}}+O(x_0^4)`$ (60)
$`𝒞_\gamma _{}(\stackrel{}{k};x_0)`$ $`=`$ $`{\displaystyle \frac{11|\stackrel{}{k}|^3}{30\pi \mu }}+{\displaystyle \frac{3|\stackrel{}{k}|^3}{2\mu }}\left({\displaystyle \frac{𝐋_3(|\stackrel{}{k}|x_0)I_3(|\stackrel{}{k}|x_0)}{(|\stackrel{}{k}|x_0)^2}}\right)`$ (61)
$`+|\stackrel{}{k}|^2({\displaystyle \frac{|\stackrel{}{k}|}{4\mu }}`$ $``$ $`{\displaystyle \frac{3}{2}}{\displaystyle \frac{d}{d\beta }}\left)\right({\displaystyle \frac{𝐋_2(\beta )I_2(\beta )}{\beta }})_{\beta =|\stackrel{}{k}|x_0}+\mu |\stackrel{}{k}|{\displaystyle \frac{d^2}{d\beta ^2}}(𝐋_1(\beta )I_1(\beta ))_{\beta =|\stackrel{}{k}|x_0}`$
$``$ $`{\displaystyle \frac{3|\stackrel{}{k}|}{2\pi \mu x_0^2}}+{\displaystyle \frac{6}{\pi |\stackrel{}{k}|x_0^3}}+{\displaystyle \frac{12\mu }{\pi |\stackrel{}{k}|^3x_0^4}}+O((\mu |\stackrel{}{k}|x_0^4)^1,(|\stackrel{}{k}|^3x_0^5)^1,\mu (|\stackrel{}{k}|^5x_0^6)^1)`$
$`𝒞_{\gamma _0}(\stackrel{}{k};x_0)`$ $`=`$ $`{\displaystyle \frac{2\mu |\stackrel{}{k}|}{\pi }}\left[1+{\displaystyle \frac{\pi }{2}}\left(𝐋_1(|\stackrel{}{k}|x_0)I_1(|\stackrel{}{k}|x_0)\right)\right]{\displaystyle \frac{2\mu }{\pi |\stackrel{}{k}|x_0^2}}+O(x_0^4)`$ (62)
where $`I`$ and L are respectively modified Bessel and Struve functions. Fortunately their difference is a tabulated function , which enables the asymptotic form of the second equality in eqns. (59 \- 62) to be determined. The decay in the timelike direction is generically $`x_0^2`$, the interesting exception being for $`\mathrm{\Gamma }=\gamma _{}`$ where there are significant corrections at shorter distances.

$`\left|\stackrel{}{k}\right|2\mu `$:

In this case the $`\theta `$ functions in (56) yield unity for all $`|\stackrel{}{p}|\mu `$, simplifying the integration limits. We will restrict analysis to the case $`\mathrm{\Gamma }=\gamma _5`$:

$`{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\theta }{2\pi }}{\displaystyle _0^\mu }p𝑑p\left[1{\displaystyle \frac{p^2+p|\stackrel{}{k}|\mathrm{cos}\theta }{p^2(1+\frac{2|\stackrel{}{k}|}{p}\mathrm{cos}\theta +\frac{|\stackrel{}{k}|^2}{p^2})^{\frac{1}{2}}}}\right]e^{px_0\left[\sqrt{1+{\scriptscriptstyle \frac{2|\stackrel{}{k}|}{p}}\mathrm{cos}\theta +{\scriptscriptstyle \frac{|\stackrel{}{k}|^2}{p^2}}}1\right]}=`$
$`{\displaystyle \frac{\mu ^2}{\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑u{\displaystyle _0^1}𝑑vv\left[1{\displaystyle \frac{1+\frac{v+\kappa }{v\kappa }u^2}{[(1+u^2)(1+(\frac{v+\kappa }{v\kappa })^2u^2)]^{\frac{1}{2}}}}\right]{\displaystyle \frac{e^{\mu x_0\left[\sqrt{{\scriptscriptstyle \frac{(v\kappa )^2+(v+\kappa )^2u^2}{1+u^2}}}v\right]}}{1+u^2}},`$
where in the second line we have used the substitutions $`u=\mathrm{cot}\frac{1}{2}\theta `$, $`v=p/\mu `$, and $`\kappa =|\stackrel{}{k}|/\mu 2`$. To focus on the large-$`x_0`$ behaviour we now employ Laplace’s method, expanding the integrand in powers of $`u^2`$ and noting that since $`v\kappa <0`$ the negative branch of the square roots should be chosen. The dominant contribution is

$$\frac{2\mu ^2}{\pi }_0^1𝑑vve^{\mu x_0(\kappa 2v)}_{\mathrm{}}^{\mathrm{}}𝑑ue^{\mu x_0\frac{2\kappa v}{\kappa v}u^2}=\frac{\mu ^2}{\sqrt{\pi }}\left(\frac{2}{\kappa \mu x_0}\right)^{\frac{1}{2}}e^{\kappa \mu x_0}_0^1𝑑vv^{\frac{1}{2}}(\kappa v)^{\frac{1}{2}}e^{2\mu x_0v}.$$
(64)
The final integral is estimated by replacing the factor $`(\kappa v)^{\frac{1}{2}}`$, which is non-singular within the integration range, by $`(\kappa 1)^{\frac{1}{2}}`$; we obtain the lower bound

$$𝒞_{\gamma _5}(\kappa ;x_0)\frac{\mu ^2}{\sqrt{\pi }}\left(\frac{2(\kappa 1)}{\kappa \mu x_0}\right)^{\frac{1}{2}}e^{\kappa \mu x_0}(2\mu x_0)^{\frac{3}{2}}\gamma (\frac{3}{2},2\mu x_0)$$
(65)
where $`\gamma `$ is the incomplete Gamma function, whose asymptotic form for large $`\mu x_0`$ yields the final answer

$$𝒞_{\gamma _5}(\stackrel{}{k};x_0)\frac{\mu ^2}{2\sqrt{\pi }}\frac{e^{(|\stackrel{}{k}|2\mu )x_0}}{(\mu x_0)^{\frac{3}{2}}}\left(1+O(\mu x_0)^1\right).$$
(66)
For $`|\stackrel{}{k}|>2\mu `$ the correlator thus decays exponentially; however for the critical value $`|\stackrel{}{k}|=2\mu `$ the decay is algebraic, and moreover slower than that for $`|\stackrel{}{k}|\mu `$ (59).