\thanks[cor]{Corresponding author; Email: viollier@physci.uct.ac.za}
\begin{abstract}
The dynamics of a self-gravitating cold Fermi gas is described
using the analogy with an interacting self-gravitating
Bose condensate having the same Thomas-Fermi limit.
The dissipationless formation of a heavy neutrino star
through gravitational collapse
and ejection of matter
is demonstrated numerically.
Such neutrino stars offer an alternative to
black holes for the supermassive
compact dark objects at the centers of galaxies.
\end{abstract}
\begin{keyword}
neutrino stars \sep boson stars \sep
self-gravitating systems \sep galactic centers
\PACS
02.60.Cb \sep 95.30.Lz \sep 95.35.+d \sep 98.35.Jk
\end{keyword}
\end{frontmatter}
Supermassive neutrino stars, in which self-gravity is balanced by the
degeneracy pressure of the cold fermions, have been a subject of
speculation for more than three decades \cite{1}. Originally, these
objects were proposed as models for dark matter in galactic halos and
clusters of galaxies, with neutrino masses in the $\sim$ eV
range. More recently, however, degenerate superstars composed of
weakly interacting fermions in the $\sim$ 10 keV range, have been
suggested as an alternative to the supermassive black holes that are
purported to exist at the centers of galaxies \cite{2,3,4,5,6}.
In fact, it has been shown \cite{4} that such degenerate fermion stars
could explain the whole range of supermassive compact dark objects
which have been observed so far, with masses ranging from $10^{6}$ to
$3 \times 10^{9}$$M_{\odot}$, merely assuming that a weakly
interacting quasi-stable fermion of mass $m_{f} \simeq$ 15 keV exists
in nature.
As an example, the most massive compact dark object ever
observed, is located at the center of M87, with a mass $M \simeq 3.2
\times 10^{9} M_{\odot}$ \cite{7}.
Interpreting this object as a relativistic
fermion star at the Op\-pen\-hei\-mer-Volkoff \cite{8} limit, yields a
fermion mass of $m_{f} \simeq 15$ keV and
a fermion star radius of $R = 4.45 R_{\mathrm{S}}
\simeq 1.5$ light-days \cite{3,4}, where
$R_{\mathrm{S}}$ is the Schwarzschild
radius.
In this case there is little difference between the fermion
star and black hole scenarios,
because the radius of the last stable orbit around a
Schwarzschild black hole is $R=3R_{\mathrm{S}}$ anyway.
Extrapolating this down to the compact dark object at the center of
our galaxy \cite{9}, which,
having a
mass $M \simeq 2.6 \times 10^{6} M_{\odot}$,
is at the lower limit of the mass range
of the observed compact dark objects,
we obtain, using
the same fermion mass,
a radius
$R\simeq 20$ light-days $\simeq 5\times10^{4} R_{\mathrm{S}}$
\cite{2}.
As the potential inside such a nonrelativistic fermion star is shallow,
the spectrum of radiation emitted by accreting baryonic matter is cut
off for frequencies larger than $10^{13}$ Hz \cite{3,5}, as is
actually
observed in the spectrum of the
enigmatic radio source Sgr A$^{*}$ at the
galactic center \cite{10}.
A fermion star
with radius $R$
\raisebox{-1.5mm}{$\stackrel{\textstyle <}{\sim}$}
20 light-days and mass
$M \simeq 2.6 \times 10^6 M_{\odot}$
is also consistent
\cite{6} with the observed motion of stars within a projected distance
of 6 to 30 light-days from Sgr A$^{*}$ \cite{9}.
Of course, it is well-known that 15 keV lies squarely in the
cosmologically forbidden mass range for stable active neutrinos $\nu$
\cite{11}.
The existence of such a massive active neutrino is also
disfavoured by the Super-Kamiokande data
\cite{fuk}.
However,
as shown by Shi and Fuller
\cite{12} for an initial lepton asymmetry of $\sim
10^{-3}$, a sterile neutrino $\nu_{s}$ of mass $m_{s} \sim 10$ keV
may be resonantly produced in the early universe
with near closure density, i.e., $\Omega \simeq 1$.
The resulting energy spectrum is not thermal
but rather cut off so that it
approximates a degenerate Fermi gas.
Sterile neutrinos in this mass range are also
constrained
by astrophysical bounds on
the radiative decay
$\nu_{s} \rightarrow \nu \gamma$ \cite{13}.
However, the allowed parameter space includes $m_{s} \simeq 15$ keV,
contributing $\Omega_{s} \simeq 0.3$
to the critical density,
as favoured by the BOOMERANG data
\cite{14}.
The statics of degenerate fermion stars is well understood, being the
Op\-pen\-hei\-mer-Volkoff equation in the relativistic case \cite{8}, or the
Lan\'{e}-Emden equation with polytropic index $n = 3/2$ in the
nonrelativistic limit \cite{15}. Alternatively,
one may understand these as the Thomas-Fermi theory applied
to self-gravitating fermion systems.
The extension of the Thomas-Fermi theory
to finite temperature \cite{16,23} has been used
to show that, at a certain
critical temperature, weakly interacting massive fermionic matter
undergoes a first-order gravitational phase transition from a diffuse
to a clustered state, i.e. a nearly degenerate fermion star.
However, such
studies do not bear on the crucial dynamical question of
whether the fermion star can form through gravitational collapse of
density fluctuations in an orthodox cosmological setting.
Indeed,
since collisional damping is negligible, one would expect that only a
virialized cloud results \cite{11}.
$N$-body simulations of the collisionless Boltzmann or Vlasov equation
evidence a rather different picture: the collapse is followed by a
series of bounces, with matter expelled at each bounce, leaving behind
a condensed object \cite{18}.
By Liouville's theorem, the Vlasov
equation describes an incompressible fluid in phase-space so that it
respects a form of the exclusion principle.
Hence, these $N$-body
simulations are effectively fermion simulations.
What transpires is
that gravity, being attractive, self-organises the phase-space fluid
into a high-density/momentum core at the expense of other
low-density/momentum regions as seen in the evolution of the spherical
Vlasov equation \cite{19}.
Much the same behaviour is observed in the formation of mini-boson
stars through so-called gravitational cooling \cite{20}.
Such a boson
star is stable by balancing the
uncertainty and gravitational pressures.
A similar mechanism works in the presence of a quartic self-interaction
\cite{21}, which dominates over uncertainty pressure, resulting in an
equilibrium radius $R \gg 1/m$, where $m$ is the boson mass
\cite{22}.
Hence, we have a universal description of the physics underlying
the formation process: once the collapse proceeds far enough,
uncertainty, interaction or degeneracy pressure results in a bounce,
in which
the outgoing shock wave carries away the binding energy.
The virial counter argument mentioned
above is bypassed, because the ejected matter invalidates
its assumption that there is no flow through the boundary.
The purpose of
this letter is to verify this picture for the formation of a
fermion star from a cold gravitationally unstable configuration.
The dynamical
Thomas-Fermi theory, developed long ago by Bloch for the electron gas
\cite{23}, amounts to Euler's equations for irrotational flow,
subject to an equation of state $P = P(\rho)$.
The problem is that,
due to the presence of shocks
and instabilities,
a naive integration of the Bloch equations is
precluded
in the gravitational case.
The usual remedy is to introduce some small numerical
viscosity. However, it seems imprudent to draw conclusions based on
introducing an {\em ad-hoc}
dissipation into what is fundamentally a dissipationless
system. Here we take another, literally conservative
approach.
We base our approach on the equivalence of a degenerate
non-interacting fermion star
with a certain type of self-interacting cold boson star.
We shall demonstrate this equivalence using scaling arguments
which are similar in spirit
to the analysis of self-interacting boson stars
by Colpi, Shapiro and Wassermann
\cite{22}.
Consider a complex scalar field $\psi$ with a
repulsive Lagrangean interaction term $U(|\psi|^2)$.
For simplicity, we use the Newtonian approximation,
but our analysis may be easily extended in a general relativistic
context.
In the Newtonian limit, a self-interacting boson star is governed by
the Gross-Pitaevskii-like equations
%.............
\begin{equation}
i \frac{\partial \psi}{\partial t} = \left[ - \frac{\Delta}{2 m} +
\frac{1}{m} \frac{\d U}{\d |\psi|^2}
+ m \varphi \right] \, \psi,
\label{dpsi}
\end{equation}
\begin{equation}
\Delta \varphi = 4 \pi G \rho,
\label{delta}
\end{equation}
\begin{equation}
\rho=m^2|\psi|^2.
\label{rhopsi}
\end{equation}
The pressure tensor is given by
\begin{equation}
\label{eq102}
P_{ij}=
\mathrm{Re}\,
\frac{\partial\psi}{\partial x_i}
\frac{\partial\psi^*}{\partial x_j}
+\delta_{ij}\left(|\psi|^2\frac{\d U}{\d |\psi|^2} -U
-\frac{\Delta|\psi|^2}{4}\right).
\end{equation}
Following Colpi et al. \cite{22}, we introduce
a dimensionless parameter $\Lambda$ and define
\begin{equation}
\label{eq104}
R_{*} \equiv \frac{\Lambda}{m} \equiv
\lambda\frac{m_{\mathrm{Pl}}}{m^2};
\;\;\;
M_{*} = R_{*} m^2_{\mathrm{Pl}}\,,
\end{equation}
as the length and mass scales, respectively.
Here $m_{\mathrm{Pl}} \equiv 1/\sqrt{G}$ denotes the Planck mass
and the arbitrary constant $\lambda$
will be fixed later using physical arguments.
%EQ 8
The substitution
\begin{equation}
\psi = \frac{m}{\sqrt{4\pi}\lambda} \Psi
\end{equation}
yields the coupled dimensionless equations
%EQ 9
\begin{equation}
\frac{i}{\Lambda} \; \frac{\partial \Psi}{\partial t} =
\left[ - \frac{\Delta}{2 \Lambda^{2}} + \varphi +
V(|\Psi|^2)
\right] \Psi,
\label{eq009a}
\end{equation}
\begin{equation}
\label{eq009b}
\Delta \varphi = |\Psi|^2 ,
\end{equation}
where we have introduced the dimensionless potential
\begin{equation}
V(|\Psi|^2)=
\frac{4\pi\lambda^2}{m^4} \frac{\d U}{\d |\Psi|^2}.
\end{equation}
The pressure tensor (\ref{eq102}), written in terms of
dimensionless variables, reads
\begin{equation}
\label{eq109}
P_{ij}=
\frac{m^4}{4\pi\lambda^2}\frac{1}{\Lambda^2}
\left(
\mathrm{Re}\,
\frac{\partial\Psi}{\partial x_i}
\frac{\partial\Psi^*}{\partial x_j}
-\delta_{ij}\frac{\Delta|\Psi|^2}{4} \right)
+\delta_{ij}\left(|\Psi|^2\frac{\d U}{\d |\Psi|^2} -U\right).
\end{equation}
A static, spherically symmetric solution,
usually referred to as boson star, is obtained
by the ansatz
%
\begin{equation}
\Psi =\e^{-{\mathrm i}\epsilon R_* t}\Phi(r),
\label{ansatz}
\end{equation}
where $\Phi(r)$ is a real function.
It is clear from the definition (\ref{eq104})
that $\Lambda \gg 1$, as long as $m/m_{\mathrm{Pl}}\ll \lambda$
which is a reasonable assumption.
Hence, for sufficiently large $\Lambda$,
Eq. (\ref{eq009a}) with (\ref{ansatz}) degenerates to
\begin{equation}
\label{eq105}
\frac{\epsilon}{m} -\varphi -
V(\Phi^2)=0\, .
\end{equation}
This equation is exact in the limit
$\Lambda\rightarrow \infty$, which is just
the Thomas-Fermi limit \cite{24}.
In this limit
the derivative terms in Eq. (\ref{eq109}) vanish
so that the pressure tensor becomes diagonal,
with all the components equal to
$P\equiv P_{ii}$.
We obtain an effective equation of state given by
\begin{equation}
\label{eq106a}
\rho = \frac{m^4}{4\pi\lambda^2}\Phi^2,
\end{equation}
\begin{equation}
\label{eq106b}
P=
\rho V -U.
\end{equation}
It may be easily shown that the last equation
combined with (\ref{eq105}) yields the equation
of hydrostatic
equilibrium
\begin{equation}
\label{eq107}
\frac{dP}{\rho}=-d\varphi.
\end{equation}
For a particular $U(|\psi|^2)$, the equation of
state in the form $P=P(\rho)$ is obtained by eliminating $\Phi^2$ from
(\ref{eq106a}) and
(\ref{eq106b}).
On the other hand,
if the equation of state,
e.g., of the polytropic type,
is given,
we can determine the potential $U$ and its
derivative $V$
by integrating (\ref{eq107}).
For a general polytropic equation of state
%EQ 5
\begin{equation}
P (\rho) = K \rho^{1+1/n},
\end{equation}
Eq. (\ref{eq106b}) yields
%EQ 6
\begin{equation}
V=(n + 1) K \rho^{1/n} \; .
\label{eq6}
\end{equation}
This, together with Eq. (\ref{eq106a}), gives the potential in terms
of $\Phi^2$.
As we still have the freedom to choose
the parameter $\lambda$ conveniently,
we fix $\lambda$ so that
the potential $V$ in (\ref{eq009a}) takes a simple form
\begin{equation}
\label{eq108}
V=|\Psi|^{2/n} .
\end{equation}
The polytropic equation of state with $n=3/2$, together with
the equation of hydrostatic equilibrium (\ref{eq107}) and
Poisson's equation (\ref{delta}), describes a degenerate
fermion star.
Hence, we have demonstrated that a
degenerate fermion star
is equivalent to a self-interacting boson star
in the limit $\Lambda\rightarrow \infty$.
Moreover, it may be shown numerically that,
even for moderate values of $\Lambda$ of the order of 10 to 20,
the static solutions
are almost degenerate and
are quite well
approximated by
the static solution for an infinite $\Lambda$.
This has been demonstrated
for a quartic self-interaction \cite{22}
(note that our $\Lambda$ corresponds
to $\sqrt{\Lambda}$ in reference \cite{22}).
Next we proceed to solving
Eqs. (\ref{eq009a}),
(\ref{eq009b})
and
(\ref{eq108})
numerically
for large, but finite value of
$\Lambda$.
For weakly interacting degenerate fermions with
the polytropic index $n = 3/2$,
the length and mass scales are given by
\begin{equation}
R_{*} = \left( \frac{9 \pi^{2}}{32 g_{f}^{2}} \right)^{1/4} \;
\frac{m_{\mathrm{Pl}}}{m_{f}^{2}} = 0.2325
\left( \frac{\mbox{keV}}{m_{f}} \right)^{2} \sqrt{\frac{2}
{g_{f}}} \; \mathrm{lyr} \; ,
\end{equation}
\begin{equation}
M_{*} = 1.489 \left( \frac{\mathrm{keV}}{m_{f}} \right)^{2}
\sqrt{ \frac{2}{g_{f}}} \times 10^{12} \; M_{\odot} \; ,
\end{equation}
where $g_{f}$ is the spin degeneracy and $m_f$ the
mass of the fermion.
As $M_*$ is of the order of the
Oppenheimer-Volkoff limit,
the validity of the Newtonian
approximation in the static case requires
\begin{equation}
M=(4 \pi)^{-1} \int \d^{3}r \, | \Psi |^{2} \ll 1 .
\end{equation}
By construction,
once the overall scale is specified,
a large but finite $\Lambda$ allows us to
simulate the fermionic problem as a bosonic one,
as long as their
Thomas-Fermi limits coincide,
while providing an explicitly energy conserving
way of controlling the shocks
and instabilities.
The basic regulating mechanism is the
kinetic part of (\ref{eq009a})
which penalises density spikes.
Of course,
$\Lambda$ must be so large that this term does not change
the static scaling relationship \cite{15}
\begin{equation}
M R^{\left( \frac{3-n}{n-1} \right)} = C_{n},
\label{eq012}
\end{equation}
arising from the polytropic equation of state.
Our criterion is that
the ratio of kinetic and pressure contributions to the static
energy functional should be small.
In particular for
a Gaussian $\Psi
= \alpha\exp\left[ - (r/\beta)^{2}\right]$
we find
%EQ 13
\begin{equation}
\frac{1}{2\Lambda^2}
\frac{\int \d^3r |\nabla \Psi|^2}{
\int \d^3r |\Psi|^{2+2/n}}
= \left( 1 + \frac{1}{n} \right)^{3/2}
\frac{3}{2 \Lambda^{2} \beta^{2} \alpha^{2/n}} \; .
\end{equation}
For $n = 3/2$ this is
independent of the size of $\beta$ for a given mass,
yielding the weak condition $M \gg 0.91 \Lambda^{-3}$.
When the mass violates this inequality
the self-interaction $V$ becomes irrelevant and
there is a crossover to mini-boson star behaviour
$M R=$ const.
\begin{figure}[h]
\centering
\epsfig{file=fig1.eps,height=8.0cm,width=12.5cm,angle=270}
\caption{
Combined
contour-density plot for the evolution of $|r \Psi |^{2}$
from the initial configuration
described in the text.
Light contour lines denote levels
from $10^{-5}$ to $10^{-4}$ while black lines denote
levels below $10^{-5}$.
%The lightest gray area correspond to a threshold value of
%$10^{-6}$, while the darkest area corresponds to levels
%at or exceeding $10^{-5}$.
Gravitational collapse is followed by ejection of
excess matter, leaving a fermion star at the center.}
\label{fig1}
\end{figure}
In Fig.\ \ref{fig1} we display the evolution of $|r \Psi |^{2}$ for
the spherical collapse of a mass $M=0.008$, initially in the form of
the above
Gaussian with $\beta = 100 = \Lambda$.
A fairly coarse sampling is
used to make this contour plot, as this allows one to see the basic
features of matter ejection without the complicated fine detail hiding
it. Eqs. (\ref{eq009a}) and
(\ref{eq009b})
are first written in terms of $\eta = r
\Psi$, and then solved assuming spherical symmetry using the method of
lines \cite{25} for $r = 0$ to 1000. The region is divided into 4000
intervals,
and a fourth order polynomial is fitted in each interval, with
the function continuous up to its second derivatives across each
interval.
The time step is adjusted dynamically to achieve the desired
accuracy of $10^{-6}$. The boundary conditions at $r = 0$ and $r =
1000$ are reflective. To prevent matter from being artificially
reflected by the boundary at $r = 1000$, we have
introduced an
$r$-dependent imaginary part to the potential,
or `sponge' \cite{21},
from $r =$ 700 to 1000, which removes the ejected fermionic matter.
The mass in the region remains
constant to within the numerical integration
error until $t = 5 \times 10^{4}$, which is a quarter of the
simulation period,
when the first ejected matter reaches the
`sponge'.
The expected features of bounce and ejection, leaving a
condensed core, are evident.
We obtain identical results using
the Crank-Nicholson method.
\begin{figure}
\centering
\epsfig{file=fig2.ps,height=9.0cm}
\caption{$| r \Psi |^{2}$ (solid line) and $\xi$ (dashed line)
versus $r$ at $t = 1.22 \times 10^{4}$.}
\label{fig2}
\end{figure}
\begin{figure}
\centering
\epsfig{file=fig3.ps,height=9.0cm}
\caption{ $| r \Psi |^{2}$ versus $r$ on the time slice
$t=2\times 10^5$
in Figure 1.}
\label{fig4}
\end{figure}
A pertinent question is whether the Thomas-Fermi
approximation
is satisfied during the ejection process.
Similar to WKB,
the Thomas-Fermi limit
is obtained when the scale of variation of density
is small compared to that of the potential $V$, i.e.,
\begin{equation}
\xi^2 \equiv
\frac{ |\nabla\rho|^2}{
\rho^2\, 8\, \d U/\d |\psi|^2}=
\frac{|\nabla |\Psi||^2}{
2\Lambda |\Psi|^{2+2/n}}
\ll 1.
\end{equation}
In Fig.\ \ref{fig2} we
exhibit $|r \Psi|^{2}$ and $\xi$ versus $r$ at time $1.22
\times 10^{4}$, where ejecta first develop.
Evidently, the Thomas-Fermi
condition, is well satisfied, except near sharp density minima.
This is expected since in a pure Thomas-Fermi system evolving under
Bloch's equations, the sharp density gradient would evolve into a shock,
which violates the Thomas-Fermi condition.
In Fig. \ref{fig4} we show $|r \Psi|^{2}$ on the time slice
$t=2\times 10^5$ of
our $\Lambda=100$ simulation. The core mass $M = 0.0057$ and radius
$R = 28$ are roughly commensurate with (\ref{eq012}) and \cite{15}
$C_{3/2} = 132.3843$. Of course, at this point there is still evolution
and oscillation of the core, similar to the mini-boson star case
\cite{20}, with $R$ varying between about 25 and 31. In general one can
expect a persistent seismology with a period $T$ of the order of the
radius divided by the average speed of sound, $T \sim 51/M$, which
agrees with the simulations. The seismology includes
higher-order modes as well.
In summary, using a bosonic representation of the
dynamical Thomas-Fermi
theory for a self-gravitating gas, we have shown that nonrelativistic,
degenerate and weakly interacting
fermionic matter will form supermassive
fermion stars through gravitational collapse accompanied by ejection.
For a fermion mass of
$m_{f} \simeq 15$ keV,
and a final mass $M\simeq 2.6\times 10^6 M_{\odot}$
such a superstar is consistent
with the observations of the compact dark object at
the center of our galaxy.
A similar demonstration for the formation
of such a star near the Oppenheimer-Volkoff limit, and
the question of cosmology with
degenerate dark matter, requires a general
relativistic extension which is under
development and will be reported
elsewhere.
%\newpage
%________________________
%{\bf Acknowledgement}
\begin{ack}
The authors wish to thank Duncan Elliott ($\dag$)
%\footnote{Deceased 20 July 2000}
for many useful discussions regarding the simulations.
This
research is in part supported by the Foundation of Fundamental
Research (FFR) grant number PHY99-01241 and the Research Committee of
the University of Cape Town. The work of N.B. is supported in part by
the Ministry of Science and Technology of the Republic of Croatia
under Contract No. 00980102.
\end{ack}
\begin{thebibliography}{99}
\bibitem{1} M.A. Markov, Phys. Lett. {\bf 10}, 122 (1964);
R. Cowsik and J. McClelland, Astrophys. J. {\bf 180}, 7 (1973);
R. Ruffini, Lett. Nuovo Cim. {\bf 29}, 161 (1980).
\bibitem{2} R.D. Viollier, D. Trautmann and G.B. Tupper,
Phys. Lett. {\bf B306}, 79 (1993).
\bibitem{3} N. Bili\'{c}, D. Tsiklauri and R.D. Viollier,
Prog. Part. Nucl. Phys. {\bf 40}, 17 (1998);
N. Bili\'{c} and R.D. Viollier,
Nucl. Phys. (Proc. Suppl.) {\bf B66}, 256 (1998).
\bibitem{4} N. Bili\'{c}, F. Munyaneza and R.D. Viollier,
Phys. Rev. {\bf D59}, 024003 (1999).
\bibitem{5} D. Tsiklauri and R.D. Viollier,
Astropart. Phys. {\bf 12},
199 (1999).
\bibitem{6} F. Munyaneza, D. Tsiklauri and R.D. Viollier,
Astrophys. J. {\bf 509}, L105 (1998);
ibid. {\bf 526}, 744 (1999);
F. Munyaneza and R.D. Viollier,
astro-ph/0103466.
\bibitem{7} F. Macchetto et al.,
Astrophys. J. {\bf 489}, 579 (1997).
\bibitem{8} J.R. Oppenheimer and G.M. Volkoff,
Phys. Rev. {\bf 55}, 374 (1939).
\bibitem{9} A. Eckart and R. Genzel,
Mon. N.R. Astron. Soc. {\bf 284}, 576 (1997);
A.M. Ghez, B.L. Klein, M. Morris and E.E. Becklin,
Astrophys. J. {\bf 509}, 678 (1998).
\bibitem{10} R. Mahadevan, Nature {\bf 394}, 651 (1998).
\bibitem{11} E.W. Kolb and M.S. Turner, The Early Universe
(Addison-Wesley, San Francisco, 1989).
\bibitem{fuk} S. Fukuda et al.,
Phys. Rev. Lett. {\bf 85}, 3999 (2000).
\bibitem{12} X. Shi and G.M. Fuller,
Phys. Rev. Lett. {\bf 82}, 2832 (1999);
K. Abazajian, G.M. Fuller, and M. Patel,
astro-ph/0101524.
\bibitem{13} M. Drees and D. Wright, hep-ph/0006274.
\bibitem{14} P. de Bernardis et al.,
Nature {\bf 404}, 955 (2000).
\bibitem{15} S. Chandrasekhar, Stellar Structure
(University of Chicago Press, Chicago, 1939).
\bibitem{16} N. Bili\'{c} and R.D. Viollier,
Phys. Lett. {\bf B408}, 75 (1997);
Gen. Rel. Grav. {\bf 31},
1105 (1999); Eur. Phys. J. {\bf C11}, 173 (1999).
\bibitem{17} W. Thirring, Z. Physik {\bf 235}, 339 (1970);
P. Hertel, H. Narnhofer and W. Thirring, Comm. Math. Phys.
{\bf 28}, 159 (1972);
J. Messer, J. Math. Phys. {\bf 22}, 2910 (1981).
\bibitem{18} J. Binney and S. Tremaine, Galactic Dynamics
(Princeton University Press, Princeton, New Jersey, 1987),
and references cited therein.
\bibitem{19} R.N. Henrikson and L.M. Widrow,
Phys. Rev. Lett. {\bf 78}, 3426 (1997).
\bibitem{20} E. Seidel and W.-M. Suen,
Phys. Rev. Lett. {\bf 72}, 2516
(1994).
\bibitem{21} J. Balakrishna, E. Seidel and W.-M. Suen,
Phys. Rev. {\bf D58}, 104004 (1998).
\bibitem{22} M. Colpi, S.L. Shapiro and I. Wasserman,
Phys. Rev. Lett. {\bf 57}, 2485 (1986).
\bibitem{23} F. Bloch, Z. Physik {\bf 81}, 363 (1933).
\bibitem{24} In the context of Bose condensates see:
A.S. Parkins and D.F. Walls
Phys. Rep. {\bf 303}, 1 (1998).
\bibitem{25} W. Schiesser,
The Numerical Method of Lines (Academic Press,
San Diego, 1991).
\end{thebibliography}
\end{document}