text:
# New Velocity Distribution for Cold Dark Matter in the Context of the Eddington Theory

## 1 Introduction

In recent years the consideration of exotic dark matter has become necessary in order to close the Universe (Jungman et al., 1996). The COBE data (Smoot et al., 1992) suggest that CDM (Cold Dark Matter) component is at least is at least $`60\%`$ (Gawser et al., 1988; Gross et al., 1998) of the total mass. On the other hand evidence from two different teams, the High-z Supernova Search Team (Riess et al., 1998) and the Supernova Cosmology Project (Somerville et al., 2002) , (Perlmutter et al, 1999; Perlmutter et al., 1997) suggests that the Universe may be dominated by the cosmological constant $`\mathrm{\Lambda }`$. Thus the situation can be adequately described by a baryonic component $`\mathrm{\Omega }_B=0.1`$ along with the exotic components $`\mathrm{\Omega }_{CDM}=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.6`$. In a more detailed $`\mathrm{\Lambda }CDM`$ analysis (Primack, 2002) one finds:

$$\mathrm{\Omega }_b=0.040\pm 0.002,\mathrm{\Omega }_m=\mathrm{\Omega }_{CDM}=0.33\pm 0.035,\mathrm{\Omega }_{HDM}0.05,\mathrm{\Omega }_\mathrm{\Lambda }=0.73\pm 0.08$$

(see also (Turner, 1990) and (Einasto, 2001)). Since the non exotic component cannot exceed $`40\%`$ of the CDM (Jungman et al., 1996)(Alcock et al., 1995), there is room for the exotic WIMP’s (Weakly Interacting Massive Particles). In fact the DAMA experiment (Bernabei et al., 1996) has claimed the observation of one signal in direct detection of a WIMP, which with better statistics has subsequently been interpreted as a modulation signal (Bernabei et al, 1998; Bernabei et al., 1999).

In the most favored scenario of supersymmetry the LSP can be simply described as a Majorana fermion, a linear combination of the neutral components of the gauginos and Higgsinos (Jungman et al., 1996; Vergados, 1990; Gomez and Vergados, 2001; Gomez et al., 2000a, b; Gomez and Vergados, 2001; Anrowit and Nath, 1995, 1996; Bottino et al., 1997; Bednyakov et al., 1994).

Since this particle is expected to be very massive, $`m_\chi 30GeV`$, and extremely non relativistic with average kinetic energy $`T100KeV`$, it can be directly detected (Vergados, 1996; Spira et al., 1995; Kosmas and Vergados, 1997) mainly via the recoiling of a nucleus (A,Z) in elastic scattering.

In order to compute the event rate one needs the following ingredients:

1) An effective Lagrangian at the elementary particle (quark) level obtained in the framework of supersymmetry as described , e.g., in (Jungman et al., 1996; Anrowit and Nath, 1995, 1996; Bottino et al., 1997; Bednyakov et al., 1994).

2) A procedure in going from the quark to the nucleon level, i.e. a quark model for the nucleon. The results depend crucially on the content of the nucleon in quarks other than u and d. This is particularly true for the scalar couplings as well as the isoscalar axial coupling (Drees and Noijiri, 1993b, c; Cheng, 1988, 1989).

3) Compute the relevant nuclear matrix elements (Ressel et al., 1993; Divari et al., 2000) using as reliable as possible many body nuclear wave functions. The situation is a bit simpler in the case of the scalar coupling, in which case one only needs the nuclear form factor.

4) The LSP density and velocity distribution. Among other other things, since the detection rates are expected to be very small, the velocity distribution is crucial in exploiting the characteristic experimental signatures provided by the reaction, namely: a) the modulation of the event rates due to the earth’s revolution around the sun (Vergados, 1998, 1999)-(Vergados, 2000) and b) the correlation of the rates with the Sun’s direction of motion in directional experiments, i.e. experiments in which the direction of the recoiling nucleus is observed (Vergados, 1990; Buckland et al., 2000). To obtain the right density and velocity distributions is the purpose of the present paper.

In the past various velocity distributions have been considered. The most popular one is the isothermal Maxwell-Boltzmann velocity distribution with $`<\upsilon ^2>=(3/2)\upsilon _0^2`$ where $`\upsilon _0`$ is the velocity of the sun around the galaxy, i.e. $`220km/s`$. Extensions of this M-B distribution were also considered, in particular those that were axially symmetric with enhanced dispersion in the galactocentric direction (Drucker et al., 1986; Vergados, 2000). In such distributions an upper cutoff $`\upsilon _{esc}=2.84\upsilon _0`$ was introduced by hand.

Non isothermal models have also been considered. Among those one should mention the late infall of dark matter into the galaxy, i.e caustic rings (Sikivie, 1999, 1998; Vergados, 2001a; Green, 2001; Gelmini and Gondolo, 2001), and dark matter orbiting the Sun (Copi et al., 1999).

The correct approach in our view is to consider the Eddington approach (Eddington, 1916), i.e. to obtain both the density and the velocity distribution from a mass distribution, which depends both on the velocity and the gravitational potential. This approach has been extensively studied by Merritt (Merritt, 1985) and recently applied to dark matter by Ullio and Kamionkowski(Ullio and Kamiokowski, 2001)

## 2 Density Profiles

As we have seen in the introduction the matter distribution can be given as follows

$$dM=2\pi f(\mathrm{\Phi }(r),\upsilon _r,\upsilon _t)dxdydz\upsilon _td\upsilon _td\upsilon _r$$
(1)
where the function f the distribution function, which depends on $`r`$ through the potential $`\mathrm{\Phi }(r)`$ and the tangential and radial velocities $`\upsilon _t`$ and $`\upsilon _r`$ (we assume axial symmetry in velocity space). Thus the density of matter $`\rho `$ satisfies the equation

$$d\rho =2\pi f(\mathrm{\Phi }(r),\upsilon _r,\upsilon _t)\upsilon _td\upsilon _td\upsilon _r$$
(2)
It is more convenient instead of the velocities to use the total energy $`E`$ and the angular momentum $`J`$ via the equations

$$J=\upsilon _tr,2E=\upsilon _r^2+\frac{J^2}{r^2}+2\mathrm{\Phi }(r)$$
(3)
The use of these variables, which are constants of motion, is very useful, when one wants to study equilibrium states. We thus find

$$\rho =\frac{2\pi }{r^2}\frac{f(E,J)J}{\sqrt{2(E\mathrm{\Phi }(r))J^2/r^2}}𝑑J𝑑E$$
(4)
The limits of integration for $`E`$ are from $`\mathrm{\Phi }`$ to $`0`$ and for $`J`$ from $`0`$ to $`[2r^2(E\mathrm{\Phi }(r))]^{1/2}`$.

Following Eddington we will choose a distribution function of the form

$$f(E,J)=K_\lambda (2E)^\lambda $$
(5)
($`E`$ is negative for a bound system), where $`\lambda `$ is a parameter, which will depend on the type of matter, and $`K_\lambda `$ is a normalization constant, which will be related to the density at some point. With this choice of the distribution function it is quite straightforward to find the relationship between the density $`\rho `$ and the potential. The result is

$$\rho =K_\lambda 2^{\lambda +3/2}\pi |\mathrm{\Phi }(r)|^{\lambda +3/2}\beta (\lambda +1,3/2)$$
(6)
with

$$\beta (a,b)=\frac{\mathrm{\Gamma }(a)\mathrm{\Gamma }(b)}{\mathrm{\Gamma }(a+b)}$$
(7)
Eq. (6) can be cast in the simple form:

$$\rho =\rho _\lambda (0)(\chi (x))^{\lambda +3/2}$$
(8)
with

$$\chi (x)=\frac{\mathrm{\Phi }(x)}{\mathrm{\Phi }_0},x=r/r_s$$
(9)
with $`\mathrm{\Phi }_0=\mathrm{\Phi }(0)`$ and $`r_s`$ the galactic radius (position of the sun). The constant $`K_\lambda `$ is related to the density at the origin via the relation:

$$K_\lambda =\frac{\rho _\lambda (0)}{\pi (2|\mathrm{\Phi }_0|)^{\lambda +3/2}\beta (\lambda +1,3/2)}$$
(10)
The above distribution function is isotropic. Following Merritt (Merritt, 1985) we can introduce an anisotropy by modifying the distribution function as follows:

$$f(E,J)=K_\lambda (2E)^\lambda (1\pm \frac{J^2}{r_a^2})$$
(11)
Instead of the parameter $`r_a`$ we find it convenient for our applications later on (see below) to adopt the more recent conventions and write the above equation as follows:

$$f(E,J)=K_\lambda (2E)^\lambda [1+\alpha _s\frac{J^2}{(r_s\upsilon _m)^2}]$$
(12)
where $`\upsilon _m`$ is the maximum velocity allowed by the potential, to be specified below, and $`\alpha _s`$ the asymmetry parameter. Proceeding as above we find that this induces a correction to the density of the form:

$$\mathrm{\Delta }\rho =K_\lambda \frac{4}{3}2^{\lambda +3/2}\pi |\mathrm{\Phi }(r)|^{\lambda +5/2}\beta (\lambda +1,5/2)\frac{\alpha _s}{\upsilon _m^2}x^2$$
(13)
Combining Eqs (6) and (13) we get

$$\rho _\lambda (r)=\rho _\lambda (0)\psi _\lambda (x)$$
(14)
$$\psi _\lambda (x)=(\chi (x))^{\lambda +3/2}[1+\frac{4}{3}a\frac{\beta (\lambda +1,5/2)}{\beta (\lambda +1,3/2)}x^2\chi (x)]$$
(15)
with $`a=\alpha _s|\mathrm{\Phi }_0|/\upsilon _m^2`$.

Since the scale of the potential appears only via the parameter $`a`$ one, in principle, could have two $`a`$ parameters, one for Matter ($`a_m`$) and one for dark matter ($`a_{dm}`$). In the present work we will assume that they are equal.We remind the reader that $`\alpha _s`$ is the asymmetry parameter to be treated phenomenologically.

## 3 Allowed Density Functions

The central question is to specify density $`\psi _\lambda (x)`$, the potential $`\chi (x)`$ and the mass density distribution entering Eq. (1). One can adopt one of two procedures:

1) Start out with a given density, obtained, e.g., phenomenologically, and find the potential by solving Poisson’s equation. This way one obtains, at least parametrically, e.g. with the radial coordinate as a parameter, the proper relation between the density and the potential. This approach has been adopted by a lot of researchers, see e.g. (Widrow, 2000), (Evans, 1994), (Henrriksen and Widrow, 1994), (Ullio and Kamiokowski, 2001). From this one can obtain the mass density function $`f(\mathrm{\Phi },\upsilon )`$. This approach leads, in general, to non analytic, i.e. complicated, velocity distribution, hard to implement in dark matter calculations (Ullio and Kamiokowski, 2001), especially for not spherically symmetric mass distributions.

2) Use the above simple relation between the density and the potential and solve the differential equation resulting from Poisson’s equation and thus obtain each one of them. One hopes that this way one will get a density distribution, which describes adequately dark matter distribution globally, via fitting the rotational curves, and in our vicinity. This approach can easily deal with asymmetric density distributions. The simplicity of the relation between the density and the potential thus yields the bonus that the obtained velocity distribution is analytic and can easily implemented in obtaining the event rates for direct dark matter detections.It can also be extended to include more than one power in $`\lambda `$. In this paper we will follow this approach, and leave it for the future to use semi-analytic mass distributions obtained in the approach outlined in the previous paragraph.

For a spherically symmetric potential Poisson’s equation leads to a differential equation of the type:

$$x\chi ^{^{\prime \prime }}(x)+2\chi ^{^{}}(x)=\mathrm{\Lambda }x(\chi (x))^n[1+\frac{2a}{n+1}x^2\chi (x)]$$
(16)
with $`n=\lambda +1/2`$ and

$$\mathrm{\Lambda }=\frac{4\pi G_Nr_s^2\rho _\lambda (0)}{\mathrm{\Phi }_0}$$
(17)
with $`G_N`$ Newton’s gravitational constant. The dimensionless quantity $`\mathrm{\Lambda }`$ is assumed to be positive (attractive potential). Since the asymmetry parameter is assumed to be small, the first term in the right hand side dominates. Introducing the variable $`\xi =\sqrt{\mathrm{\Lambda }}x`$ and $`\chi (x)=u(x\sqrt{\mathrm{\Lambda }})`$ we arrive at:

$$\xi u^{^{\prime \prime }}(\xi )+2u^{^{}}(\xi )=\xi (u(\xi ))^n[1+\frac{2b}{n+1}\xi ^2u(\xi )],b=\frac{a}{\mathrm{\Lambda }}$$
(18)
The last equation is a non linear differential equation, which must be solved for $`\xi >0`$, with the conditions $`u(0)=1`$ and $`u^{^{}}(0)=0`$. We also demand that the solution remains positive, i.e. the solution drops from unity to zero. For $`\alpha _s=0`$ the density is a rapidly decreasing function of $`\xi `$.

The above equation can be solved analytically in the special case $`a=0`$ and only when $`\lambda =1/2`$ (n=1) or $`\lambda =7/2`$ (n=5). In these cases the solutions are:

$$\lambda =1/2u(\xi )=\frac{sin\xi }{\xi }$$
(19)
and

$$\lambda =7/2u(\xi )=\frac{1}{(1+\xi ^2/3)^{1/2}}$$
(20)
The case $`\lambda =1/2`$ may not have physical meaning, but it will serve as an illustrative example for the realistic cases to be discussed below.

We should stress that the last solution associated with $`\lambda =7/2`$, admissible over all space, can be considered as providing an adequate description of ordinary matter111The designation of $`\lambda =7/2`$ as representing ordinary matter is due to the work of Eddington (Eddington, 1916), where he points out that this value leads to the Plummer distribution. Furthermore this distribution provides an adequate description of the galaxy (of ordinary matter) he was concerned with. A discussion of some of the thermodynamical aspects of the $`\lambda =7/2`$ distribution is also given in this paper.. In the absence of asymmetry no constraint is imposed on the parameter $`\mathrm{\Lambda }`$ by the positivity and finiteness condition of the solution. The solution given by Eq. (19) is constrained in the range $`0\xi \pi `$. If we demand that the range of the potential for dark matter extends far out in the halo, e.g. up to $`x_r=r/r_s=20`$, we must have $`\mathrm{\Lambda }=[\pi /x_r]^2=0.025`$. This gives a useful constraint between $`\mathrm{\Phi }_0`$ and $`\rho _\lambda (0)`$.

The differential equation (18) can, in general, be solved only numerically. For $`\lambda =1/2`$ we find that there is a singularity at $`\xi =12.43`$, while the first zero of the potential occurs at $`\xi =4.35`$. This leads to the constraint $`\mathrm{\Lambda }=(4.35/20)^2=0.047`$. For $`\lambda =1`$ the corresponding value is $`\mathrm{\Lambda }=0.065`$.

The above results are significantly modified when the asymmetry parameter is turned on. One now can see that, depending on the parameters $`\mathrm{\Lambda }`$ and $`\alpha _s`$, the shape of the density is significantly modified.

We will examine the special case of $`\lambda =7/2`$, which corresponds to the Eddington solution for ordinary matter. As we have seen in the absence of asymmetry the solution can be found exactly. In the presence of asymmetry we distinguish two cases:

1) Positive values of $`\alpha _s`$.

For small values of $`b`$ the solution has no roots. Beyond a critical value of $`b`$ the solution attains the value zero. Then further out for still larger $`\xi `$ it becomes negative. At some point it becomes singular. The extracted values of $`\mathrm{\Lambda }`$ range between $`1.5`$ and 0.12.

2) Negative values of $`\alpha _s`$.

In this case when the absolute value of $`b`$ becomes sufficiently large the solution goes through zero. Again a value of $`\mathrm{\Lambda }`$ can be extracted. The situation is quite unstable. Generally speaking for arbitrary values of $`b`$ the solution blows up to infinity without going through zero. One cannot extract values of $`\mathrm{\Lambda }`$ in this case.

In any case, these observations are intended as illustrations of what is expected for dark matter where the exact solution cannot be obtained. In the case of ordinary matter we will not consider $`\alpha _s`$ different from zero, since the range of the solution is small ( the density is assumed to vanish for $`rr_s`$).

Let us return to Eq. (12) and express it in terms of the velocity in our vicinity. We obtain:

$$f(\upsilon )=N[2\mathrm{\Phi }(r_s)\upsilon ^2]^\lambda (1+\alpha _s\frac{\upsilon _t^2}{\upsilon _m^2})$$
(21)
where $`N`$ is a normalization constant, which depends on $`\lambda ,\alpha _s`$ and $`\upsilon _m`$. From this we see that the maximum velocity in our vicinity , $`\upsilon _m`$, depends on the value of the potential, i.e.

$$\upsilon _m^2=2\mathrm{\Phi }(r_s)=2\mathrm{\Phi }_0u(\sqrt{\mathrm{\Lambda }})$$
(22)
This means that

$$b=\frac{a}{\mathrm{\Lambda }}=\frac{\alpha _s}{2\mathrm{\Lambda }u(\sqrt{\mathrm{\Lambda }})}$$
(23)
Our strategy is clear:

1) Vary $`b`$ so that the solution has its first zero at some value $`\xi _0`$ and is monotonically decreasing up to that point.

2) Determine $`\mathrm{\Lambda }`$ so that $`\xi _0`$ corresponds to a range $`x_r`$ of the potential $`\chi (x)`$.

3) From Eq. (23) obtain an acceptable (self-consistent) value of $`\alpha _s`$, $`\alpha _s=2\mathrm{\Lambda }bu(\sqrt{\mathrm{\Lambda }})`$.

4) The parameter $`\upsilon _m`$ is obtained via the relation:

$$\upsilon _m^2==2\mathrm{\Phi }_0u(\sqrt{\lambda })=8\pi G_Nr_s^2\frac{u(\sqrt{\mathrm{\Lambda }})\rho _\lambda (0)}{\mathrm{\Lambda }}$$
(24)
or equivalently

$$\upsilon _m=\upsilon _0\sqrt{\eta _\chi }h_\upsilon (n,a,\mathrm{\Lambda }),h_\upsilon (n,a,\mathrm{\Lambda })=[\frac{u(\sqrt{\mathrm{\Lambda }})}{\mathrm{\Lambda }((u(\sqrt{\mathrm{\Lambda }}))^n[1+\frac{2a}{n+1}u(\sqrt{\mathrm{\Lambda }})]}]^{1/2},n=\lambda +3/2$$
(25)
where $`\upsilon _0`$ is our velocity of rotation around the center of our galaxy (220 $`km/s`$). In obtaining the last equation we used our solution to relate the density $`\rho _\lambda (0)`$ to the density of dark matter in our vicinity, $`\rho _0`$, via the equation:

$$\rho _\lambda (0)=\frac{\rho _0}{((u(\sqrt{\mathrm{\Lambda }}))^n[1+\frac{2a}{n+1}u(\sqrt{\mathrm{\Lambda }})]},n=\lambda +3/2$$
(26)
The parameter $`\eta _\chi `$ is defined by $`\eta _\chi =\rho _0/(0.3GeV/cm^3)`$. $`\eta _\chi `$ is normally taken to be unity (Jungman et al., 1996).

It is clear that in our treatment for a given $`\lambda `$ we have three independent parameters, $`b`$, $`x_r`$ and $`\rho _0`$. The obtained results for some interesting typical cases of $`\alpha _s`$ are shown in Table 1. The value of the parameter $`h_\upsilon `$ has the meaning of $`y_{esc}`$ (see section 6) and it is rewarding that it is not far from the value of 2.84 deduced phenomenologically (Drucker et al., 1986).

The potentials obtained for the various such parameters are shown in Figs 1 \- 4.

It is clear that the shapes of these potentials are similar.

The densities obtained with the same set of parameters are shown in Figs 5 \- 8.

We observe that, even though the potentials are similar, the densities differ substantially. The asymmetry parameter has a big effect on the density, especially for $`\alpha _s`$ greater than zero, when it changes substantially around $`x=8`$. For $`\alpha _s`$ negative there is small reduction of the density in the region of interest, but, as we have already mentioned, only special values of $`b`$ are acceptable. A detailed rigorous mathematical study of the behavior of the solutions for $`b`$ negative is currently under study and it will appear elsewhere.

## 4 Rotational Velocities

We are now in position to obtain the rotational velocity curves. The rotational velocity is given by

$$\upsilon _\lambda (x)=\frac{\upsilon _0}{\sqrt{2}}\sqrt{g_\lambda (x)}$$
(27)
with

$$g_\lambda (x)=\frac{1}{x}_0^x(x^{^{}})^2𝑑x^{^{}}\psi _\lambda (x^{^{}})$$
(28)
We see that the scale of the rotational velocity is set by $`\frac{\upsilon _0}{\sqrt{2}}`$ in agreement with the data, see, e.g., the recent review (Jungman et al., 1996). For dark matter the above integrals can only be performed numerically. The obtained rotational velocities are shown in Figs 9 \- 11.

For ordinary matter ($`\lambda =7/2)`$ we make the choice of $`\mathrm{\Lambda }=3`$ to get the familiar solution

$$\chi (x)=\frac{1}{(1+x^2)^{1/2}},\psi _{7/2}=\psi (x)=\frac{1}{(1+x^2)^{5/2}}$$
(29)
(We will drop the index $`7/2`$ and use the index $`m`$, for matter, when necessary). Thus the rotational velocity due to ordinary matter is now given by

$$\upsilon _m(x)=\frac{\upsilon _0}{\sqrt{2}}\sqrt{g_m(x)}$$
(30)
For $`x<1`$ we write

$$g_m(x)=g(x)$$
(31)
while for $`x>1`$ we have:

$$g_m(x)=\frac{1}{x}g(1)$$
(32)
We now distinguish two cases:

i) Spherical galaxy.

Using the above density we find:

$$g(x)=\frac{1}{3}\frac{x^2}{(1+x^2)^{3/2}}$$
(33)
ii) Spiral galaxy (disk).

In this case we will use the same density as above. In the Eddington theory the equation, which relates the potential and the density is no longer of the above simple form. One finds:

$$\frac{\mathrm{\Phi }}{\mathrm{\Phi }_0}=\frac{\mathrm{\Lambda }}{12}[2(\frac{\rho }{\rho (0)})^5\mathrm{ln}[\frac{1(\rho /\rho (0))^5}{1+(\rho /\rho _0)^5}]]$$
(34)
This complexity is of no concern to us since we do not need to use this complicated equation. In the case of ordinary matter we do not need to obtain the velocity distribution. Thus with the assumed density we find:

$$g(x)=\frac{1}{3x}[1\frac{1}{(1+x^2)^{3/2}}]$$
(35)
The above functions are plotted in Fig. 12. In the graphs we show not only the physically interesting case, in which the density vanishes outside the radius of the galaxy, but, for illustrative purposes, the case of infinite extent of the same density function . We see that, since the density falls very fast, it makes very little difference, which form we use.

By comparing the above graphs of dark matter with those of ordinary matter we see that in the Eddington theory ordinary matter can make a significant contribution to the rotational velocities, only if the density at the origin is much bigger compared to that of dark matter. As we have mentioned above for ordinary matter there is no constraint between $`\mathrm{\Phi }_0`$ and $`\rho _0`$. So we simply rescale ordinary matter by a factor $`C_{mdm}`$ and write:

$$g_{matter}(x)=C_{mdm}g_m(x)$$
(36)
the scaling factor can be determined by a comparison to the experimental rotational curves.

Thus the rotational velocity due to both matter and dark matter is given by

$$\upsilon _\lambda (x)=\frac{\upsilon _0}{\sqrt{2}}\sqrt{g_\lambda (x)+g_{matter}(x)}$$
(37)
The obtained rotational velocities for two values of $`C_{mdm}`$ are shown in Figs 13 \- 16

The larger value of $`C_{mdm}`$ seems to be in better agreement with the data, see, e.g., the recent review (Jungman et al., 1996).

It is clear from these curves that the presence of asymmetry has a dramatic effect on the rotational velocities. Thus, in the context of the Eddington theory as discussed here, a large asymmetry is excluded by the data on rotational velocities. A detailed fitting of our input parameters to the rotational curves will not be attempted here.

## 5 Velocity distribution with respect to the galactic center

The above density via the Eddington formula leads to a velocity distribution of the form:

$$F(\upsilon ,r)[2\mathrm{\Phi }(r)\upsilon ^2]^\lambda (1+\alpha _s\frac{\upsilon _t^2}{\upsilon _m^2})$$
(38)
The above velocities and the distance $`r`$ are defined with respect to the center of the matter distribution, i.e to the center of our galaxy. We note in the context of the Eddington theory the velocity distribution cannot be Maxwellian. For a given distance it goes to zero at the boundaries of the corresponding ellipsoid. It is customary to consider the value of the above distribution in our vicinity, $`r=r_s`$. This way it reduces to the product of the local density and the velocity distribution. The latter is given by Eq. (21) where $`N`$ is a normalization constant, which depends on $`\lambda ,\alpha _s`$ and $`\upsilon _m`$. The above notation was introduced to make the last equation coincide with the standard expression when the function f is chosen to be Maxwellian, i.e.

$`Exp(\frac{\upsilon ^2\alpha _s\upsilon _t^2}{\upsilon _0^2})`$ $`(1+\alpha _s\frac{\upsilon _t^2}{\upsilon _0^2})`$ $`Exp(\frac{\upsilon ^2}{\upsilon _0^2})`$

for sufficiently small $`\alpha _s`$. In this limit we see that $`\alpha _s`$ coincides with the parameter $`\lambda `$ of Vergados (Vergados, 2000) et al Drucker et al. (1986) (in the present work $`\lambda `$ is used for another purpose).

It is straightforward to find that the normalization factor $`\stackrel{~}{N}`$ is given by

$$\stackrel{~}{N}^1(\lambda ,\alpha _s,\upsilon _m)=2\pi \upsilon _m^{2\lambda +3}\beta (\lambda +1,3/2)[1+\frac{4}{3}\alpha _s\frac{\beta (\lambda +1,5/2}{\beta (\lambda +1,3/2})]$$
(39)
with $`\upsilon _m`$ given by Eq. 22. In the special case of dark matter ($`\lambda =1/2)`$ it becomes

$$\stackrel{~}{N}^1(1/2,\alpha _s,\upsilon _m)=\frac{\pi ^2}{4}\upsilon _m^4[1+\frac{1}{3}\alpha _s]$$
(40)
From the above formulas we see that the velocity of dark matter with respect to the galactic center ranges from $`0`$ to a maximum speed $`\upsilon _m=|2\mathrm{\Phi }(r_s)|^{1/2}`$. Since the escape velocity is determined by the potential $`\mathrm{\Phi }(r_s)`$ due to all kinds of matter, the velocity $`\upsilon _m`$ is not simply related to $`\upsilon _{esc}`$. Thus, since the distribution function must remain positive, if $`\alpha _s<0`$ its absolute value is bounded. This imposes a constraint, since in the traditional analysis with only axially symmetric Gaussian distribution it leads to negative $`\alpha _s`$, i.e. enhanced dispersion in the galactocentric direction, a phenomenologically preferred result (Drucker et al., 1986). The data on rotational curves may provide an additional constraint on the negative values of $`\alpha _s`$. For positive values of $`\alpha _s`$ the constraint coming from the rotational curves is, as we have seen, more stringent.

From then on one proceeds in the usual way to obtain the velocity distribution with respect to the laboratory.

## 6 Velocity distribution with respect to the laboratory

For this transformation one needs the velocity of the sun around the galaxy $`\upsilon _0=220Km/s`$, a fraction of the escape velocity, which is $`\upsilon _{esc}=625Km/s=2.84\upsilon _0`$ (Drucker et al., 1986).

It is convenient to choose coordinate system with its polar z-axis in the the direction of the disc’s rotation, i.e. in the direction of the motion of the the sun, the x-axis is in the outward radial direction and the y-axis is perpendicular to the plane of the galaxy. Since the axis of the ecliptic (Kosmas and Vergados, 1997). lies very close to the $`y,z`$ plane the velocity of the earth around the the center of the galaxy is given by

$$𝝊_E=\upsilon _0\widehat{𝐳}+𝝊_1=\upsilon _0\widehat{𝐳}+\upsilon _1(\mathrm{sin}\alpha \widehat{𝐱}\mathrm{cos}\alpha \mathrm{cos}\gamma \widehat{𝐲}+\mathrm{cos}\alpha \mathrm{sin}\gamma \widehat{𝐳})$$
(41)
where $`\alpha `$ is the phase of the earth’s orbital motion, $`\alpha =2\pi (tt_1)/T_E`$, where $`t`$ is the time of observation and $`t_1`$ is around second of June and $`T_E=1`$ year. The magnitude of the Earth’s velocity is much smaller than that of the sun, i.e. $`\delta _1=2\upsilon _1/\upsilon _0=0.27`$ The velocity of the earth around its own axis is even smaller and it is usually neglected.

One can now express the above distribution in the laboratory frame by writing $`𝝊^{^{}}=𝝊+𝝊_E`$ , where the prime indicates the velocity with respect to the center of the galaxy. Indicating by $`y`$ the velocity of the LSP (Lightest Supersymmetric Particle) in units of $`\upsilon _0`$, i.e. by defining $`y=\upsilon /\upsilon _0`$, we find

$$f(y,\theta ,\varphi )=N(\lambda ,\alpha _s,y_m)[y_m^2Y(y,\theta ,\varphi )]^\lambda [1+\alpha _s(Y(y,\theta )(ysin\theta \mathrm{cos}\varphi \frac{\delta _1}{2}\mathrm{sin}\alpha )^2)]$$
(42)
with

$$N(\lambda ,\alpha _s,y_m)=\upsilon _0^{2\lambda +3}\stackrel{~}{N}(\lambda ,\alpha _s,\upsilon _m)$$
(43)
where $`\stackrel{~}{N}(\lambda ,\alpha _s,\upsilon _m)`$ is given by Eq. 39, $`y_m=\upsilon _m/\upsilon _0`$ and

$`Y(y,\theta ,\varphi )`$ $`=`$ $`1+{\displaystyle \frac{\delta _1^2}{4}}+y^2+2y\mathrm{cos}\theta +\delta _1[y\mathrm{cos}\theta \mathrm{cos}\alpha \mathrm{sin}\gamma `$ (44)
$``$ $`y\mathrm{sin}\theta \mathrm{sin}\varphi \mathrm{cos}\alpha \mathrm{cos}\gamma +y\mathrm{sin}\theta \mathrm{cos}\varphi \mathrm{sin}\alpha ]`$
with In the conventional axially symmetric Gaussian velocity distribution this function is given by

$`f(y,\theta ,\varphi )`$ $`=`$ $`{\displaystyle \frac{N(\alpha _s,y_{esc})}{\pi \sqrt{\pi }}}Exp[(\alpha _s+1)Y(y,\theta ,\varphi )`$
$`+`$ $`\alpha _s(ysin\theta \mathrm{cos}\varphi {\displaystyle \frac{\delta _1}{2}}\mathrm{sin}\alpha )^2)]`$
In this case $`0yy_{esc}`$, but the upper cutoff is introduced here artificially. The normalization here is defined so that $`N(\alpha _s=0,y_{esc}\mathrm{})=1`$ (Vergados, 2000)

The detection rate in direct dark matter experiments is obtained by convoluting the the relevant cross section with the above velocity distribution. If the dark matter candidate is the LSP, the $`\alpha `$-dependence of the above distribution , present only when $`\delta _10`$, gives rise to the modulation effect, i.e. the dependence of the rate on the Earth’s motion. This signal can be used to discriminate against background.

### 6.1 The non directional rate

The non directional differential event rate is given by:

$$\frac{dR}{du}=\overline{R}\sqrt{\frac{2}{3}}T(u),T(u)=a^2|F(u)|^2\mathrm{\Psi }(a\sqrt{u})$$
(46)
for the coherent mode and

$$\frac{dR}{du}=\overline{R}\sqrt{\frac{2}{3}}T_{spin}(u),T_{spin}=a^2|F_{11}u)|\mathrm{\Psi }(a\sqrt{u})$$
(47)
where $`\overline{R}`$ and $`\overline{R}_{spin}`$ are the total rates for the coherent and the spin contributions associated with an average LSP velocity, $`\sqrt{<\upsilon ^2>}=\sqrt{(3/2)}\upsilon _0`$. These parameters, which carry the dependence on the SUSY parameters, are the most important ones, but they are not of interest in our present calculation. $`F(u)`$ is the form factor, entering the coherent scattering and $`F_{11}(u)`$ is the spin response function entering via the axial current (Vergados, 2001b). The function $`\mathrm{\Psi }`$ depends on the LSP distribution velocity employed and is a function of the energy $`Q`$ transferred to the nucleus

$$u=\frac{Q}{Q_0},Q_0=4.1\times 10^4A^{4/3}KeV$$
(48)
where $`A`$ is the nuclear mass number and the parameter $`a`$ is given by:

$$a=[\sqrt{2}\mu _rb\upsilon _0]^1$$
(49)
where $`\mu _r`$ is the reduced mass of the LSP-nucleus system and $`b`$ is the (harmonic oscillator) size parameter.

The function, which is basic to us, $`\mathrm{\Psi }`$, is given by

$$\mathrm{\Psi }(x)=_0^x𝑑y_0^\pi \mathrm{sin}\theta d\theta _0^{2\pi }𝑑\varphi yf(y,\theta ,\varphi )$$
(50)
with $`0xy_m1+(\delta _1/2)\mathrm{cos}\alpha \mathrm{sin}\beta `$.

The total rate is given by:

$$R=_{u_{min}}^{u_{max}}\frac{dR}{du}𝑑u$$
(51)
where $`u_{min}`$ is determined from the cutoff energy of the detector and $`u_{max}=(y_m/a)^2`$

From now on we will specialize our results to the case $`\lambda =1/2`$, which yields an adequate description of the rotation curves in the case of dark matter. In the case of this velocity distribution, unlike the Gaussian one, one cannot approximate the distribution by a power series in $`\delta _1`$. The reason is that there may be threshold problems, when the argument of the square root approaches zero. To simplify matters we will still make use of the fact that the velocity of the Earth around the Sun is much smaller than the velocity of the Sun around the galaxy, $`\delta _1<<1`$. So, if we expand the previous expression into a Fourier series with respect to the phase of the Earth, $`\alpha `$, only the lowest terms will become important. In other words to leading order in $`\delta _1`$ it can be put in the form:

$$R=\overline{R}[R_0+(R_1\mathrm{cos}\alpha \mathrm{sin}\gamma R_2\mathrm{cos}\alpha \mathrm{cos}\gamma +R_3\mathrm{sin}\alpha )\delta _1/2]$$
(52)
It turns out that the expansion coefficients $`R_2`$ and $`R_3`$ are zero. We can thus conveniently fit the rate with the formula:

$$R=\overline{R}t[1+h\mathrm{cos}\alpha ]$$
(53)
where $`h`$ is the modulation amplitude (the difference between the maximum and the minimum is equal to $`2|h|`$).

In the case of no modulation, $`\delta _1=0`$, the angular integrals can be done analytically to yield:

$$\mathrm{\Psi }(x)=2\pi N(1/2,\alpha _s,y_m)[\frac{(1+\alpha _sy_m^2)}{3}J_3(x,y_m)\frac{\alpha _s}{20}J_5(x,ym)]$$
(54)
with $`0xy_m1`$ and

$$J_n(x,y_m)=J_{int}(n,y_m,x1)J_{int}(n,y_m,x+1)+2J_{int}(n,y_m,1)$$
(55)
$$J_{int}(n,y_m,y)=_0^y[y_m^2z^2]^{n/2}𝑑z$$
(56)
The above integral can be done analytically to yield:

$$J_{int}(n,y_m,y)=y(y_m)^n{}_{2}{}^{}F_{1}^{}(\frac{1}{2},\frac{n}{2},\frac{3}{2},\frac{y^2}{y_m^2})$$
(57)
where $`{}_{2}{}^{}F_{1}^{}`$ is the usual hypergeometric function. For the cases of interest to us here the hypergeometric function can be simplified to yield:

$$J_{int}(3,y_m,y)=\frac{1}{8}[2y(y_m^2y^2)^{3/2}+3y_m^2y(y_m^2y^2y^2)^{1/2}+3y_m^4\mathrm{sin}^1(y/y_m)]$$
(58)
$`J_{int}(5,y_m,y)`$ $`=`$ $`{\displaystyle \frac{1}{24}}[4y(y_m^2y^2)^{5/2}+5y_m^2y(y_m^2y^2)^{3/2}`$
$`+`$ $`15y_m^4y(y_m^2y^2)^{1/2}+15y_m^6\mathrm{sin}^1(y/y_m)]`$
$`J_{int}(7,y_m,y)`$ $`=`$ $`{\displaystyle \frac{1}{192}}[24y(y_m^2y^2)^{7/2}+28y_m^2y(y_m^2y^2)^{5/2}+35y_m^4y(y_m^2)^{1/2}`$
$`+`$ $`105y_m^6y(y_m^2)^{1/2}+105y_m^8\mathrm{sin}^1(y/y_m)]`$
The corresponding expressions for the Gaussian expressions for $`\alpha _s`$ cannot be done analytically. For the symmetric case, $`\alpha _s=0`$, one finds that:

$$\mathrm{\Psi }(x)=\frac{1}{2}[erf(x1)erf(x+1)]+2erf(1)$$
(61)
with $`erf(x)`$ the error function:

$$erf(y)=\frac{2}{\pi }_0^ye^{t^2}𝑑t$$
(62)
The above functions $`\mathrm{\Psi }`$ are plotted in Figs 17 and 18.

From Fig. 18 we see that both the value of the function $`\mathrm{\Psi }(x)`$ as well as its range depend crucially on the parameters of the model. In the absence of asymmetry the Gaussian model, with an appropriate escape velocity put in by hand, yields results, which are almost four times larger than those of the Eddington theory. In the presence of asymmetry our results in the Eddington theory are substantially larger. This is due to the larger peak value attained as well as the larger allowed range of x. On the other hand we know that in the Gaussian model the total rate is not significantly affected by the asymmetry (Vergados, 1999). The reason for the strong dependence of the results in the present model on the asymmetry is not, of course, the asymmetry per se, but the fact that all parameters change with it. In particular the value of $`y_m`$ (see table 1). In the Gaussian model the introduction of asymmetry did not affect the velocity distribution in any other way, e.g. it did not affect the cutoff value of the velocity distribution. The large upper values of $`x`$ in the function $`\mathrm{\Psi }(x)`$, allowed in the present model for large $`y_m`$, are, of course, somewhat controlled by the nuclear form factor. In order to get a better feeling for such an effect on the rate, we also plot the function $`T(u)`$, which is proportional to the differential non directional rate. For illustration purposes we have chosen to present results for the popular target $`{}_{}{}^{127}I`$ and a typical LSP mass of $`100GeV`$ (see Fig. 19). The results presented are for the coherent mechanism, but we expect very small changes when the spin contribution is considered. It is clear that the introduction of asymmetry has a profound effect on the rate. We have seen, however, that the large positive values of the $`\alpha _s`$ can, in principle, be eliminated from the data on the rotational curves.

Effects like those discussed above may be more pronounced in the case of modulation, not studied in this work.

### 6.2 The directional rate

The directional differential rate $`(dR/du)_{dir}`$ is proportional to

$`T_{dir}(u)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}\overline{R}a^2|F(u)|^2\mathrm{\Psi }_{dir}(a\sqrt{u})`$
$`T_{dir}^{spin}`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}\overline{R}_{spin}a^2|F_{11}u)|\mathrm{\Psi }_{dir}(a\sqrt{u})`$ (63)
with

$$\mathrm{\Psi }_{dir}(x)=_0^x𝑑y_0^\pi sin\theta 𝑑\theta _0^{2\pi }𝑑\varphi yf(y,\theta ,\varphi )XH(X)$$
(64)
with $`H(X)`$ the well known Heaviside function and $`X`$ is given by

$$X=\mathrm{cos}\mathrm{\Theta }\mathrm{cos}\theta +\mathrm{sin}\mathrm{\Theta }sin\theta [\mathrm{sin}\mathrm{\Phi }\mathrm{sin}\varphi +\mathrm{cos}\mathrm{\Phi }\mathrm{cos}\varphi ]$$
(65)
where $`\mathrm{\Theta }`$ and $`\mathrm{\Phi }`$ describe the direction of observation $`\widehat{e}`$

$$\widehat{e}=\mathrm{sin}\mathrm{\Theta }(\mathrm{cos}\mathrm{\Phi }\widehat{e}_x+\mathrm{sin}\mathrm{\Phi }\widehat{e}_y)+\mathrm{cos}\mathrm{\Theta }\widehat{e}_z$$

(There should be no confusion of the angle $`\mathrm{\Phi }`$ used here with the potential $`\mathrm{\Phi }(r)`$ used earlier). Note the presence of the factor of $`1/(2\pi )`$, since the azimuthal integration of the recoiling nucleus is not present and we intend to use the same nucleon cross-section both in the directional and the non directional case.

The total rate is proportional to:

$$R_{dir}=_{u_{min}}^{u_{max}}(\frac{dR}{du})_{dir}𝑑u$$
(66)
Taking again the lowest Fourier components of the obtained rate as a function of the phase of the Earth we get an expression similar to Eq. 52. Thus, to leading order in $`\delta _1`$, we can fit the total rate by an expression of the form:

$$R_{dir}=\frac{1}{2\pi }\overline{R}t_{dir}[1+(h_1h_2)\mathrm{cos}\alpha +h_3\mathrm{sin}\alpha ]$$
(67)
The parameters $`t_{dir}`$ and $`h_i`$ $`i=1,2,3`$ are obviously functions of the direction of observation, i.e. $`\mathrm{\Theta }`$ and $`\mathrm{\Phi }`$. If one observes in the direction of the Sun’s velocity $`h_2=h_3=0`$. Similarly if one observes in a plane perpendicular to the Sun’s velocity $`h_1=0`$. Instead of $`t_{dir}`$ it is best to use the ratio:

$$\kappa =2\pi \frac{R_{dir}}{R}=\frac{t_{dir}}{t}$$
(68)
The parameter $`\kappa `$ is essentially independent of the LSP mass, the nuclear parameters and the asymmetry parameter $`\alpha _s`$. But it depends strongly on the direction of observation and is expected to correlate strongly with the angle between $`\widehat{e}`$ and the Sun’s direction of motion. This correlation provides a an experimental signature perhaps better the modulation with the Earth’s motion in non directional experiments.

For $`\delta _10`$ the above integrals over $`y,\theta ,\varphi `$, especially in the directional case, can only be done numerically. Such results will appear elsewhere (Braun et al., 2002).

## 7 Conclusions

In the present paper we studied the density and velocity distributions of cold dark matter in the context of the Eddington theory, considering not only symmetric but axially symmetric distributions as well. In our approach we used standard simple distribution functions of the energy and angular momentum. This lead us to simple relations between the density $`\rho `$ and the potential $`\mathrm{\Phi }`$. $`\rho `$ and $`\mathrm{\Phi }`$ were obtained by solving numerically poisson’s equation in a suitable region of space. This procedure allowed us to determine the maximum permitted dark matter velocity. Then we were able to study both the rotational curves as well as the velocity distribution of dark matter in our vicinity. We should mention that this distribution is not Maxwellian and has an upper velocity cutoff built into into it and not put in by hand as in the traditional treatment with Gaussian distribution.

Our results depend on a minimum set of parameters, which were treated as free (see table 1).

We saw that in the context of this theory the predicted rotational velocities , see Figs 13 \- 16, depending on the input parameters, can vary significantly in the presence of the asymmetry. This comes mainly from the factor $`x^2`$ in the second term of Eq. 15, when the parameter $`a`$ is reasonably large. By comparing our results to the observed rotational curves (Jungman et al., 1996), one may constrain the parameters of the model. The best choice seems to be the case with a small asymmetry parameter $`\alpha _s`$.

We have also made a preliminary study of the effect of the new velocity distribution on the direct detection rates for cold dark matter. For illustration purposes we have selected the case $`\lambda =1/2`$. We have seen that, in the context of the Eddington approach, the total rates, unlike the case of the Gaussian distribution (Vergados, 1999, 2000), sensitively depend on the asymmetry parameter $`\alpha _s`$. This is due to the fact that, when the asymmetry changes, the upper value of the velocity distribution also changes.

It is thus not surprising that, in the Eddington theory, the total (non directional and non modulated) event rates for direct LSP detection maybe substantially different from those of the phenomenological Gaussian distributions (compare Figs 17 and 18). The strong dependence of the rate on the asymmetry parameter remains even after the nuclear form factor has been incorporated (see Fig. 19). Results of more detailed calculations of the event rates will appear elsewhere.

The dependence of the directional and/or modulated rates on the velocity distribution is currently under study. We expect this dependence to be more pronounced than on the total rates.

## 8 acknowledgments

J.D.V. would like to thank the Physics Department of UNISA and Professor S. Sofianos for their hospitality. D.O appreciates the hospitality provided by the University of Ioannina.


format: text

text:
# Quark Condensates and Momentum-Dependent Quark Masses in a Nonlocal Nambu–Jona-Lasinio Model

## I INTRODUCTION

The Nambu–Jona-Lasinio model has been extensively studied for several decades \[1-3\]. In recent years there has been strong interest in the study of quark matter at high densities, using the NJL model and related models. In particular, one finds color superconductivity under certain conditions and it has been suggested that some compact stars might be made of superconducting quark matter \[4-10\]. It is our belief that in such studies one should use a model which reproduces, as well as possible, know features of QCD. In this work we wish to generalize the SU(3)-flavor version of the NJL model to be consistent with the type of momentum-dependent quark masses found in lattice simulations of QCD . For example, in Figs. 1 and 2 we show some of the results obtained in Ref. . It may be seen from these figures that, in Euclidean space, the quark mass goes over to the current mass for Euclidean momentum $`k`$ 2 GeV. On the other hand, the NJL model, in the standard analysis , gives rise to a constant value for the constituent quark mass. As we will see, the form of the nonlocality used here is different from that used in Refs. \[6, 12-16\]. For example, in reference  the $`q\overline{q}`$ vertex is modified by a form factor which depends on the relative momentum of the quark and antiquark, while, in another scheme, a form factor is associated with each quark line appearing in a diagram. In the latter procedure no further regularization is needed. However, for the problem considered in this work, the nonlocal models that appear in the literature are of limited applicability, since, in the limit of zero current quark mass, the quark self-energy is proportional to the form factors used to define the nonlocality . In contrast, in the current work, we calculate the form of the quark self-energy after introducing a regulator and a nonlocal quark interaction. We stress that the procedure used here differs from any that appears in the literature.

The organization of our work is as follows. In Section II we review the standard analysis for the condensates and “gap equation” of the SU(3)-flavor NJL model . We then go on to review a procedure for including the contribution to the quark self-energy due to the addition of a confining interaction for the case of the SU(2)-flavor model. In Section III we introduce a nonlocal interaction in the SU(3)-flavor NJL model while maintaining the separable nature of the interaction. We also describe the approximation used for the ’t Hooft interaction in the nonlocal model. In Section IV and V we present some of the results of our numerical calculations. In Section VI we discuss the density dependence of the quark condensate and the momentum dependent quark mass. Finally, Section VII contains some additional discussion and conclusions.

## II The quark self-energy in the NJL model

In order to best introduce the nonlocal model, we will first review the calculation of the quark self-energy in the local SU(3)-flavor NJL model and then proceed to add a confinement interaction, as was done in an earlier study of the quark self-energy . We consider the generalization to a model with nonlocal short-range and ’t Hooft interactions in the next section. The Lagrangian of the model is

$`=`$ $`\overline{q}(i\text{ / }m^0)q+{\displaystyle \frac{G_S}{2}}{\displaystyle \underset{i=0}{\overset{8}{}}}[(\overline{q}\lambda ^iq)^2+(\overline{q}i\gamma _5\lambda ^iq)^2]`$ (2.1)
$`+{\displaystyle \frac{G_D}{2}}\{det[\overline{q}(1+\gamma _5)q]+det[\overline{q}(1\gamma _5)q]\}.`$
Here $`m^0`$ is the matrix of quark current masses, $`m^0=\text{diag }(m_u^0,m_d^0,m_s^0)`$ and the $`\lambda ^i(i=1,\mathrm{},8)`$ are the Gell-Mann matrices. Further, $`\lambda _0=\sqrt{2/3}𝟙`$, with 1being the unit matrix in the flavor space.

The quark propagator is written as

$$iS(k)=\frac{i}{\text{ / }k\mathrm{\Sigma }(k)+iϵ},$$
(2.2)
with

$$\mathrm{\Sigma }(k)=A(k^2)+B(k^2)\text{ / }k.$$
(2.3)
We may define

$$M_u(k^2)=\frac{A_u(k^2)}{1B_u(k^2)},$$
(2.4)
and

$$Z_u(k^2)=\frac{1}{1B_u(k^2)},$$
(2.5)
with similar definitions for $`M_d(k^2),M_s(k^2),Z_d(k^2)`$ and $`Z_s(k^2)`$.

In the absence of a confinement model, we have $`B(k^2)=0,A_u(k^2)=A_d(k^2)=m_u`$ and $`A_s(k^2)=m_s`$, where $`m_u`$ and and $`m_s`$ are constants. (Here, we have take $`m_u^0=m_d^0`$.) In this case we have

$$m_u=m_u^02G_S\overline{u}uG_D\overline{d}d\overline{s}s,$$
(2.6)
$$m_d=m_d^02G_S\overline{d}dG_D\overline{u}u\overline{s}s,$$
(2.7)
$$m_s=m_s^02G_S\overline{s}sG_D\overline{u}u\overline{d}d.$$
(2.8)
These equations are depicted in Fig. 3a, where the last term represents the ’t Hooft interaction. The up quark vacuum condensate is given by

$`\overline{u}u`$ $`=`$ $`N_ci{\displaystyle \frac{d^4k}{(2\pi )^4}\text{Tr}\frac{C(k^2)}{\text{ / }km_u+iϵ}},`$ (2.9)
$`=`$ $`4N_ci{\displaystyle \frac{d^4k}{(2\pi )^4}\frac{m_uC(k^2)}{k^2m_u^2+iϵ}}.`$ (2.10)
Here, $`C(k^2)`$ is a function needed to regulate the integral. In this work we will use the Pauli-Villars procedure and evaluate the integral in Euclidean space, as was done in Ref. . In the general case we may write

$$\overline{u}u=4N_ci\frac{d^4k}{(2\pi )^4}\frac{Z_u(k^2)M_u(k^2)C(k^2)}{k^2M_u^2(k^2)+iϵ},$$
(2.11)
where $`Z_u(k^2)`$ and $`M_u^2(k^2)`$ were defined in Eqs. (2.4) and (2.5).

In the appendix of Ref.  we considered Lorentz-vector confinement, with

$$\overline{V^c}(k_Ek_E^{})=\gamma ^\mu (1)\gamma _\mu (2)V^c(k_Ek_E^{})$$
(2.12)
and

$$V^c(k_Ek_E^{})=8\pi \kappa \left\{\frac{1}{[(k_Ek_E^{})^2+\mu ^2]^2}\frac{4\mu ^2}{[(k_Ek_E^{})^2+\mu ^2]^3}\right\},$$
(2.13)
where $`k_E^\mu k_E^\mu `$ denotes the Euclidean-space momentum transfer. Here, $`\mu `$ is a small parameter introduced to soften the momentum-space singularities. Note that the form

$$V^c(\stackrel{}{k}\stackrel{}{k}^{})=8\pi \kappa \left\{\frac{1}{[(\stackrel{}{k}\stackrel{}{k}^{})^2+\mu ^2]^2}\frac{4\mu ^2}{[(\stackrel{}{k}\stackrel{}{k}^{})^2+\mu ^2]^3}\right\},$$
(2.14)
represents the Fourier transform of $`V^c(r)=\kappa re^{\mu r}`$, so that for small $`\mu ,V^c(r)`$ approximates a linear potential over the relevant range of *r*. In Fig. 3b we show the equation for the self-energy when the confining field is included.

In Ref.  we obtained the following coupled equations in the case of the SU(2)-flavor model

$$A(k^2)=i\frac{d^4k^{}}{(2\pi )^4}\frac{[4V^c(kk^{})+4N_cn_fG_S]A(k^{\mathrm{\hspace{0.17em}2}})}{k^{\mathrm{\hspace{0.17em}2}}[1B(k^{\mathrm{\hspace{0.17em}2}})]^2A^2(k^{\mathrm{\hspace{0.17em}2}})+iϵ},$$
(2.15)
$$k^2B(k^2)=i\frac{d^4k^{}}{(2\pi )^4}\frac{2(kk^{})[1B(k^{\mathrm{\hspace{0.17em}2}})]V^c(kk^{})}{k^{\mathrm{\hspace{0.17em}2}}[1B(k^{\mathrm{\hspace{0.17em}2}})]^2A^2(k^{\mathrm{\hspace{0.17em}2}})+iϵ}.$$
(2.16)
These equations were solved after passing to Euclidean space and including a Pauli-Villars regulator of the form

$$C(k_E^{\mathrm{\hspace{0.17em}2}})=\frac{2\mathrm{\Lambda }^4}{[k_E^{\mathrm{\hspace{0.17em}2}}+A^2(k_E^{\mathrm{\hspace{0.17em}2}})+\mathrm{\Lambda }^2][k_E^{\mathrm{\hspace{0.17em}2}}+A^2(k_E^{\mathrm{\hspace{0.17em}2}})+2\mathrm{\Lambda }^2]}$$
(2.17)
in Euclidean space. Note that the form

$$\stackrel{~}{C}(k_E^{\mathrm{\hspace{0.17em}2}})=\frac{2\mathrm{\Lambda }^4}{[k_E^{\mathrm{\hspace{0.17em}2}}\left(1B(k_E^{\mathrm{\hspace{0.17em}2}})\right)^2+A^2(k_E^{\mathrm{\hspace{0.17em}2}})+\mathrm{\Lambda }^2][k_E^{\mathrm{\hspace{0.17em}2}}\left(1B(k_E^{\mathrm{\hspace{0.17em}2}})\right)^2+A^2(k_E^{\mathrm{\hspace{0.17em}2}})+2\mathrm{\Lambda }^2]}$$
(2.18)
may also be used. (In Eqs. (2.15) and (2.16) $`V^c`$ appears with a sign opposite to that given in Ref. , since, in that work, we used a negative value of $`\kappa `$. For the present work we use a positive value of $`\kappa `$ to be consistent with all of our other publications.)

## III A nonlocal NJL model

In this section we describe the procedures we use to create a nonlocal version of our generalized NJL model. We consider the second term of Eq. (2.1) and make the replacement

$`{\displaystyle \frac{G_S}{2}}`$ $`{\displaystyle \underset{i=0}{\overset{8}{}}}[(\overline{q}(x)\lambda _iq(x))^2+(\overline{q}(x)i\gamma _5\lambda _iq(x))^2]`$
$`{\displaystyle \frac{G_S}{2}}{\displaystyle \underset{i=0}{\overset{8}{}}}\{[\overline{q}(x)\lambda ^if(x)q(x)\overline{q}(y)\lambda ^if(y)q(y)`$
$`+[\overline{q}(x)i\gamma _5\lambda ^if(x)q(x)\overline{q}(y)i\gamma _5\lambda ^if(y)q(y)]\}.`$
This replacement corresponds to the use of a separable interaction $`V(xy)=G_Sf(x)f(y)`$.

A related modification may be made for the ’t Hooft interaction. It is useful, however, to describe these modifications as they affect momentum-space calculations. With reference to Fig. 4, we replace $`G_S`$ by $`f(k_1k_2)G_Sf(k_3k_4).`$ In the evaluation of the second term of Fig. 2 we need $`G_S(kk^{})=f(kk^{})G_Sf(kk^{})`$ and choose to write

$`G_S(kk^{})`$ $`=`$ $`\mathrm{exp}[(kk^{})^{\mathrm{\hspace{0.17em}2}n}/2\beta ]G_S\mathrm{exp}[(kk^{})^{\mathrm{\hspace{0.17em}2}n}/2\beta ]`$
$`=`$ $`G_S\mathrm{exp}[(kk^{})^{\mathrm{\hspace{0.17em}2}n}/\beta ].`$
In this work we take $`n=4`$ and $`\beta =20`$ GeV8. In Fig. 5 we exhibit the function $`F(k^2)=\mathrm{exp}[k^{\mathrm{\hspace{0.17em}2}n}/\beta ]`$. It is clear that many other functions may be chosen.

We now rewrite Eqs. (2.15) and (2.16) for the up, down and strange quarks. For example, for the SU(3)-flavor case

$`A_u(k^2)`$ $`=`$ $`m_u^0+i{\displaystyle \frac{d^4k^{}}{(2\pi )^4}\frac{[4V^c(kk^{})+8N_cG_S(kk^{})]A_u(k^{\mathrm{\hspace{0.17em}2}})}{k^{\mathrm{\hspace{0.17em}2}}[1B_u(k^{\mathrm{\hspace{0.17em}2}})]^2A_u^2(k^{\mathrm{\hspace{0.17em}2}})+iϵ}},`$ (3.3)
$`k^2B_u(k^2)`$ $`=`$ $`i{\displaystyle \frac{d^4k^{}}{(2\pi )^4}\frac{2(kk^{})[1B_u(k^{\mathrm{\hspace{0.17em}2}})]V^c(kk^{})}{k^{\mathrm{\hspace{0.17em}2}}[1B_u(k^{\mathrm{\hspace{0.17em}2}})]^2A_u^2(k^{\mathrm{\hspace{0.17em}2}})+iϵ}},`$ (3.4)
with similar equations for $`A_d(k^2),B_d(k^2),`$ etc. Again, these equations are solved after passing to Euclidean space and introducing regulator functions: $`C_u(k^2),C_d(k^2)\text{ and }C_s(k^2)`$. In Euclidean space we write

$`C_u(k^2)={\displaystyle \frac{2\mathrm{\Lambda }^4}{[k^2+A_u^2(k^2)+\mathrm{\Lambda }^2][k^2+A_u^2(k^2)+2\mathrm{\Lambda }^2]}},`$ (3.5)
etc. Note that without the ’t Hooft interaction the equations for the up, down and strange quarks are uncoupled.

Our treatment of the ’t Hooft interaction is based upon a generalization of the last term in Eqs. (2.6)-(2.8). With reference to the third term on the right in Fig. 6, we introduce a correlation between the quark of momentum *k* and the quark of momentum $`k^{}`$. We also include a correlation between the quark of momentum *k* and that of momentum $`k^{}^{}`$. (Therefore, our procedure does not introduce a correlation between the two quarks in the separate condensates. At this stage of the development of our model that seems to be a reasonable approximation and avoids having to define a three-quark correlation function, $`f(kk^{},k^{}k^{\prime \prime },kk^{\prime \prime }).`$) For example, we generalize the term $`G_D\overline{d}d\overline{s}s`$ to obtain the contribution to $`A_u(k)`$:

$`A_u^t(k)=G_D\left[4N_ci{\displaystyle \frac{d^4k^{}}{(2\pi )^4}\frac{C_d(k^{\mathrm{\hspace{0.17em}2}})Z_u(k^{\mathrm{\hspace{0.17em}2}})M_u(k^{\mathrm{\hspace{0.17em}2}})f^2(kk^{})}{k^{\mathrm{\hspace{0.17em}2}}[1B_u(k^{\mathrm{\hspace{0.17em}2}})]^2A_u^2(k^{\mathrm{\hspace{0.17em}2}})}}\right]`$ (3.6)
$`\times \left[4N_ci{\displaystyle \frac{d^4k^{\prime \prime }}{(2\pi )^4}\frac{C_s(k^{\prime \prime }{}_{}{}^{2})Z_s(k^{\prime \prime }{}_{}{}^{2})M_s(k^{\prime \prime }{}_{}{}^{2})f^2(kk^{\prime \prime })}{k^{\prime \prime }{}_{}{}^{2}[1B_s(k^{\prime \prime }{}_{}{}^{2})]_{}^{2}A_s^2(k^{\prime \prime }{}_{}{}^{2})}}\right]`$
in Minkowski space. This expression is then evaluated in Euclidean space and added to the right-hand side of Eq. (3.3). In a similar fashion, we calculate $`A_d^t(k)`$ and $`A_s^t(k)`$.

## IV Numerical Results: Condensates and Constituent Mass Values

There is a good deal of flexibility in choosing the regulators $`C_u(k^2)`$, $`C_d(k^2)`$, and $`C_s(k^2)`$. Also, various forms could be chosen for the correlation functions, $`f(kk^{})`$. Previously, in our Minkowski-space studies of the $`\eta `$ mesons we used $`G_S=11.84`$ GeV-2 and $`G_D200`$ GeV-5 . However, in that work we used a Gaussian regulator in Minkowski space so that a direct comparison with the present study can not be made. On the other hand, we do not expect to find radically different parameters, if the constituent masses in the two calculations are similar. For example, in our earlier work, in which $`m_u`$ and $`m_s`$ were parameters, we used $`m_u=0.364`$ GeV, which can be compared to the value of $`M_u(0)`$ calculated here. We have also used either $`m_s=0.565`$ GeV \[19-23\] or $`m_s=0.585`$ GeV , values which may be compared to $`M_s(0)`$.

To proceed, we take $`\mathrm{\Lambda }=1.0`$ GeV, $`G_S=13.30`$ GeV-2, $`\kappa =0.055`$ GeV2, $`m_u^0=0.0055`$ GeV, $`m_s^0=0.130`$ GeV, $`\mu =0.010`$ GeV and $`\beta =20.0`$ GeV8. We then consider values of $`G_D=0`$, $`G_D=20G_S`$, $`G_D=30G_S`$ and $`G_D=40G_S`$. The results of our calculations are given in Table I. Recall that the function $`F(k)`$ does not appear in our expression for the condensates. The calculation of the condensates includes the Pauli-Villars regulators, $`C_u(k^2)`$, $`C_d(k^2)`$ and $`C_s(k^2)`$, however. \[See Eqs. (2.9)-(2.11).\] In our calculation of the properties of the $`\eta `$ mesons  we had $`G_D/G_S1518`$, since we used $`G_D`$ values in the range $`180`$ GeV-5$`G_S220`$ GeV-5 in that work.

In order to specify a value of $`G_D`$ for this work, we note that a calculation based upon chiral perturbation theory yields $`\overline{s}s/\overline{u}u=1.689`$ . Inspection of Table I suggests that the values of $`G_D`$, other than $`G_D=0`$, given in Table I are acceptable. For $`G_D=266`$ GeV-5 we have $`M_u(0)=0.377`$ GeV and $`M_s(0)=0.555`$ GeV, which are reasonably close to the phenomenological parameters $`m_u=0.364`$ GeV and $`m_s=0.565`$ GeV used in our earlier work \[19-23\].

It is worth noting that, in standard application of the SU(3)-flavor NJL model, one finds $`\overline{s}s/\overline{u}u1.1`$ , so that the results shown in Table I are encouraging, given that the value for $`\overline{s}s/\overline{u}u`$ obtained using chiral perturbation theory is about 1.7 , as noted above.

## V numerical results:
momentum dependence of the constituent quark masses

In Fig. 7 we show $`M_u(k)`$, where *k* is the magnitude of the Euclidean momentum. The dashed line exhibits the result without the confining interaction ($`\kappa =0`$). It is interesting to see that inclusion of confinement improves the shape of the curve when we compare our results to the lattice results shown in Figs. 1 and 2. We note that $`M_u(k)`$ goes over to $`m_u^0=0.0055`$ GeV for large *k*. In Fig. 8 we show $`B_u(k)`$. (Recall that $`Z_u(k)=[1B_u(k)]^1`$.) We remark that $`B_u(k)=0`$ when $`\kappa =0`$. In Figs. 9 and 10 we show $`M_s(k)`$ and $`B_s(k)`$, respectively. As expected, we find that $`M_s(k)`$ goes over to $`m_s^0=0.130`$ GeV when *k* is large.

## VI density dependence of the quark mass and quark condensate

As stated earlier, the behavior of the NJL model for finite values of the baryon density is an extensively explored topic , with particular recent emphasis on color superconductivity \[4-10\]. In this section we explore the behavior of our model at finite baryon density. (It should be noted that the nature of the phase transition describing chiral symmetry restoration at finite density is quite model dependent. For example, the inclusion of current quark masses can change a strong first-order transition to a smooth second-order transition .)

A comprehensive study of the thermodynamics of the three-flavor NJL model has been reported in Ref. . There it is found that the up, down and strange quark masses are essentially constant up to the density where a first-order phase transition appears. At that point, the up and down quark masses drop from a value of about 380 MeV to about 30 MeV. That behavior differs from the behavior expected at low density. For example, we have the well-known relation between the value of the condensate and the baryon density of nuclear matter

$`{\displaystyle \frac{\overline{q}q_\rho }{\overline{q}q_0}}=1{\displaystyle \frac{\sigma _N\rho _B}{f_\pi ^2m_\pi ^2}},`$ (6.1)
where $`\sigma _N`$ is the pion-nucleon sigma term. This relation is valid to first-order in the density. It may be derived, in the case of nuclear matter, by writing

$`\overline{q}q_\rho =\overline{q}q_0+N|\overline{q}q|N\rho _B`$ (6.2)
and making use of the definition of the pion-nucleon sigma term, $`\sigma _N`$, and the Gell-Mann–Oakes–Renner relation. If we put $`\sigma _N=0.045`$ GeV, we have $`\overline{q}q_\rho /\overline{q}q_0=10.273\rho _B`$, where $`\rho _B`$ is in GeV3 units. For nuclear matter $`\rho _B=`$( 0.109 GeV)3, so we see that the condensate is reduced by about 35%, if we evaluate Eq. (6.1) at nuclear matter density. We can check whether the density dependence given by Eq. (6.1) is reproduced in our model, since it should not matter whether the scalar density of the background matter is generated by quarks in nucleons or by the presence of free quarks. In the former case, we may write, for the baryon density,

$`\rho _B=4{\displaystyle ^{k_F}}{\displaystyle \frac{d^3k}{(2\pi )^3}}`$ (6.3)
where the factor of 4 is arises from the product of the spin and isospin factors. In the case of quarks, we have

$`\rho _B=4N_c({\displaystyle \frac{1}{3}}){\displaystyle ^{k_F}}{\displaystyle \frac{d^3k}{(2\pi )^3}}`$ (6.4)
where, in this case, the factor of 4 again arises from the spin and isospin factor. (Both up and down quarks are present in equal numbers.) The color factor, $`N_c=3`$, is cancelled by the baryon number of 1/3 of each quark.

We need to modify the equations for the quark self-energy to take into account the presence of the Fermi seas of up and down quarks whose Fermi momentum is $`k_F`$. We take one Fermi sea to be composed of on-mass-shell up quarks with constituent mass $`M_u(0)`$. The following term is then added to the equation for $`A_u(k)`$.

$`A_u^{(\rho )}(k)=(2G_S)N_c2{\displaystyle ^{k_F}}{\displaystyle \frac{d^3k^{}}{(2\pi )^3}}{\displaystyle \frac{M_u(0)}{E_u(k)}}f^2(kk^{})`$ (6.5)
where $`E_u(k)=\left[\stackrel{}{k}^2+M_u^2(0)\right]^{\frac{1}{2}}`$. The second factor of 2 in Eq. (6.5) reflects the spin degeneracy. We note that $`M_u(0)`$ is density-dependent and could be written as $`M_u(0,\rho )`$ in keeping with the labelling of Fig. 11, where $`M_u(k,\rho )`$ is used. Also, if $`f(kk^{})=1`$, $`A_u^{(\rho )}(k)`$ would then represent $`2G_S\rho _S`$, where $`\rho _S`$ is the scalar density associated with the up-quark Fermi sea.

In Fig. 11 we show $`M(k,\rho )`$ calculated for four values of $`\rho `$ and in Fig. 12 we show *M*(0) as a function of $`k_F^3`$. Since the quarks in the Fermi sea are taken to be on-mass-shell, we would, in principle, require $`M_u(k,\rho )`$ in Minkowski space. However, since $`k_F`$ = 0.268 GeV for the case of nuclear matter, only a very modest extrapolation of the curves shown in Fig. 11 is needed for the densities considered in this work. In Fig. 13 we show the value of the up quark condensate as a function of $`k_F^3`$. (Note that $`k_F^3=19.2\times 10^3`$ GeV3 represents the density of nuclear matter.) It is seen that, for small values of the density, the density dependence of the condensate reproduces what is expected from Eq. (6.1). (If we extrapolate the curve using a linear approximation, the condensate is reduced by about 30% at nuclear matter density.)

## VII discussion

We have remarked earlier in this work that the large values of the condensate ratio $`\overline{s}s/\overline{u}u`$ seen in Table I play a role in obtaining a good fit to the mixing angles of the $`\eta `$(947) and $`\eta `$(958) mesons . To understand this remark we note that the effective singlet-octet coupling constants for pseudoscalar states are

$`G_{00}^P`$ $`=`$ $`G_S{\displaystyle \frac{2}{3}}(\alpha +\beta +\gamma ){\displaystyle \frac{G_D}{2}},`$ (7.1)
$`G_{88}^P`$ $`=`$ $`G_S{\displaystyle \frac{1}{3}}(\gamma 2\alpha 2\beta ){\displaystyle \frac{G_D}{2}},`$ (7.2)
and

$`G_{08}^P={\displaystyle \frac{\sqrt{2}}{6}}(2\gamma \alpha \beta ){\displaystyle \frac{G_D}{2}},`$ (7.3)
where $`\alpha =\overline{u}u`$, $`\beta =\overline{d}d`$ and $`\gamma =\overline{s}s`$. We take $`\alpha =\beta `$, so that

$`G_{08}^P={\displaystyle \frac{\sqrt{2}}{3}}(\gamma \alpha ){\displaystyle \frac{G_D}{2}}.`$ (7.4)
If $`\gamma =1.7\alpha `$, the result for $`G_{08}^P`$ is six times larger than when $`\gamma =1.1\alpha `$.

In addition to the effects of $`G_{08}^P`$, singlet-octet mixing is induced by the quantity

$`E_{08}(k)={\displaystyle \frac{2\sqrt{2}}{3}}[E_u(k)E_s(k)],`$ (7.5)
where $`E_u(k)=\left[\stackrel{}{k}^2+m_u^2\right]^{\frac{1}{2}}`$, etc. It is found that, since $`G_{08}^P`$ and $`E_{08}(k)`$ tend to cancel in our formalism, the significant singlet-octet mixing generated by $`E_{08}(k)`$ is reduced by the values of $`G_{08}^P`$ obtained for the larger value of the ratio $`\overline{s}s/\overline{u}u`$, with the result that we reproduce the values of the mixing angles found in other studies that make use of experimental data to obtain values for the mixing angles .

In our earlier work, which was carried out in Minkowski space, the values of $`m_u=m_d`$ and $`m_s`$ were taken as parameters. Inspection of our figures which exhibit values of $`M_u(k)`$ and $`M_s(k)`$ suggests that an extrapolation into Minkowski space may be made if $`k^2`$ is not too large. The fact that $`M_u(0)`$ and $`M_s(0)`$ are close to our phenomenological parameters for $`G_D=266`$ GeV-5 is encouraging and suggests that some support for our choice of quark mass parameters may be found in our Euclidean-space analysis.

The full consequences of separating the specification of the nonlocality of the quark interaction from the choice of the regulator of the theory should be explored more fully. Although that feature of our model introduces greater flexibility, that comes with the disadvantage of having to introduce other parameters in the model. We have made only limited variation of the form of the nonlocality and the regulator. For further applications it may be of interest to explore a more comprehensive parameter variation. It is also necessary to extend the calculations reported in Figs. 11–13 to larger values of the density than those considered here. That step will require more complex methods for solving our nonlinear equations for the self-energy than the simple iteration scheme we have used thus far.

Our work may be compared to that of Alkofer, Watson and Weigel  who have solved the Schwinger-Dyson equation using a gluon propagator whose low-momentum behavior is enhanced by a Gaussian function. (That modification requires the introduction of two phenomenological parameters .) The behavior found for $`A(k)`$ and $`B(k)`$ in Euclidean space is similar to that obtained in this work. (See Fig. 1 of Ref. .) Those authors also solve the Bethe-Salpeter equation to obtain the properties of various $`q\overline{q}`$ mesons with generally satisfactory results. It is of interest to note that the Minkowski space solution for $`A(k)`$ and $`B(k)`$ is such that the quark can go on-mass-shell. That feature may be related to our work \[18-23\] in which we use on-mass-shell quarks with masses $`m_u=m_d=0.364`$ GeV and $`m_s=0.565`$ GeV (or 0.585 GeV ) when solving the Bethe-Salpeter equation in our study of $`q\overline{q}`$ mesons.

References