\documentstyle[preprint,aps]{revtex}
%%%%%%%%%%%%%%%%%%
%\documentstyle[prb,aps]{revtex}
%\documentstyle[prl,aps,psfig]{revtex}
%%%%%%%%%%%%%%%%%%
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
\newcommand{\beqa}{\begin{eqnarray}}
\newcommand{\eeqa}{\end{eqnarray}}
\newcommand{\ba}{\begin{array}}
\newcommand{\ea}{\end{array}}
\begin{document}
\draft
%%%%%%%%%%%%%%%%%%
%\twocolumn[\hsize\textwidth\columnwidth\hsize\csname
%@twocolumnfalse\endcsname
%%%%%%%%%%%%%%%%%%
\widetext
\title{Time-Dependent Variational Approach \\
to Bose-Einstein Condensation}
\author{Luca Salasnich}
\address{Istituto Nazionale per la Fisica della Materia,
Unit\`a di Milano \\
Dipartimento di Fisica, Universit\`a di Milano \\
Via Celoria 16, 20133 Milano, Italy}
\maketitle
\begin{abstract}
We discuss the mean-field approximation for a trapped
weakly-interacting Bose-Einstein
condensate (BEC) and its connection with the exact many-body problem
by deriving the Gross-Pitaevskii action of the condensate.
The mechanics of the BEC in a harmonic potential
is studied by using a variational
approach with time-dependent Gaussian trial wave-functions.
In particular,
we analyze the static configurations, the stability
and the collective oscillations for both ground-state and vortices.
\end{abstract}
\vskip 1. truecm
\pacs{PACS numbers: 32.80.Pj; 67.65.+z; 51.10.+y}
%%%%%%%%%%%%%%%%%%
%]
%%%%%%%%%%%%%%%%%%
\narrowtext
\newpage
\section{Introduction}
Few years ago, the Bose-Einstein condensation
has been observed with alkali-metal vapors
confined in magnetic traps at very low temperatures.$^{1-3}$
By now, the total number of groups that have achived Bose-Einstein
condensation is over a dozen.
The confining potential is accurately described by a harmonic trap,
and theoretical studies of the Bose-Einstein condensate (BEC) in
harmonic traps have been performed for the ground state,
elementary excitations and vortices (for a review see Ref. 4).
\par
In the first part of this paper, we discuss the time-dependent
Gross-Pitaevskii (GP) equation,$^{5,6}$ which describes the macroscopic
wave-function of a weakly-interacting Bose condensate at zero temperature.
We deduce the GP action and the GP equation by imposing
the least action principle to the quantum many-body action of the
system within the mean-field approximation.$^{7}$
Such approach is mathematically equivalent to the field theory one,
but it avoids the tricky issues related with spontaneous symmetry breaking.
\par
The rest of the paper is devoted to the study of BEC
by using Gaussian trial wave-functions
to minimize the GP action.
Minimizing the GP action with a Gaussian trial function
makes it possible to obtain approximate solutions for the condensate
wave function, which can be treated, to a large extent, analytically.
This method has been used by several authors
to study the BEC static properties
and the dynamics near the ground-state.$^{8-14}$
Here we extend its application also to vortex states.
The presence of vortex states is a signature of the
macroscopic phase coherence of the system. Moreover,
vortices are important to characterize the superfluid properties
of Bose systems.$^{7}$
\section{Variational Principle and Mean-Field approximation}
Let us consider a $N$-body quantum system with
Hamiltonian ${\hat H}$. The exact time-dependent Schr\"odinger equation
can be obtained by imposing the quantum least action principle to the action
\beq
S= \int dt <\psi (t) | i\hbar
{\partial \over \partial t} - {\hat H} |\psi (t) > \; ,
\eeq
where $\psi$ is the many-body wave-function of the system and
\beq
<\psi (t) |{\hat A} |\psi (t) >
= \int \; d^3{\bf r}_1 ... d^3{\bf r}_N\;
\psi^*({\bf r}_1,...,{\bf r}_N,t) {\hat A}
\psi ({\bf r}_1,...,{\bf r}_N,t) \; ,
\eeq
for any quantum operator ${\hat A}$.
Looking for stationary points of $S$ with respect
to variation of the conjugate wave-function $\psi^*$ gives
\beq
i\hbar {\partial \over \partial t}\psi = {\hat H}\psi \; ,
\eeq
which is the many-body time-dependent Schr\"odinger equation.
\par
As is well known, except for integrable systems,
it is impossible to obtain the exact solution
of the many-body Schr\"odinger equation and some approximation must be used.
\par
In the mean-field approximation the total wave-function
is assumed to be composed of independent particles, i.e. it can be
written as a product of single-particle wave-functions $\phi_j$.
In the case of identical fermions, $\psi$ must be antisymmetrized.$^{7}$
By looking for stationary action with respect to variation of a
particular single-particle conjugate wave-function $\phi_j^*$ one finds
a time-dependent Hartree-Fock equation for each $\phi_j$:
\beq
i\hbar {\partial \over \partial t}\phi_j = {\delta \over \delta \phi_j^*}
<\psi | {\hat H}| \psi > = {\hat h} \phi_j \; ,
\eeq
where ${\hat h}$ is a one-body operator.
In general, the one-body operator ${\hat h}$ is nonlinear. Thus
the Hartree-Fock equations are non-linear integro-differential
equations. Note that for alkali-metal vapors, which are quite
diluite, the independent-particle approximation gives
reliable results.$^{4}$
Obviously, it is possible to go beyond the independent-particle
approximation by including correlations in the many-body
wave-function.$^{8}$ This is particularly useful to reproduce
the properties of the superfluid $^4$He, that is a strongly
interacting system.$^{9,10}$
\section{Gross-Pitaevskii Action}
In this section we discuss the zero-temperature mean-field approximation
for a system of trapped weakly-interacting
bosons in the same quantum state, i.e. a Bose-Einstein condensate.$^{7}$
In this case the Hartree-Fock equations reduce to only one equation,
the Gross-Pitaevskii equation,$^{5,6}$ which describes the dynamics
of the condensate. As previously discussed,
this equation is intensively studied
because of the recent experimental achievement of Bose-Einstein
condensation for atomic gasses
in magnetic traps at very low temperatures.$^{1-3}$
\par
The Hamiltonian operator of a system of $N$ identical
bosons of mass $m$ is given by
\beq
{\hat H}=\sum_{i=1}^N \Big( -{\hbar^2\over 2 m} \nabla_i^2
+ V_0({\bf r}_i) \Big) +
{1\over 2} \sum_{ij=1}^N V({\bf r}_i,{\bf r}_j) \; ,
\eeq
where $V_0({\bf r})$ is an external potential and $V({\bf r},{\bf r}')$
is the interaction potential.
In the mean-field approximation the totally symmetric
many-particle wave-function of the Bose-Einstein condensate reads
\beq
\psi({\bf r}_1,...,{\bf r}_N,t) = \phi({\bf r}_1,t) \; ...
\; \phi({\bf r}_N,t) \; ,
\eeq
where $\phi ({\bf r},t)$ is the single-particle wave-function.
Note that such factorization of the total wave-function is exact
in the case of a non-interacting condensate.
The quantum action of the system is then simply given by
\beq
S_{GP}=N \int dt <\phi (t) | i\hbar
{\partial \over \partial t} - {\hat h_s} |\phi (t) > \; ,
\eeq
where
\beq
{\hat h_s}= -{\hbar^2\over 2 m} \nabla^2
+ V_0({\bf r}) + {1\over 2}(N-1)
\int d^3{\bf r}' |\phi ({\bf r}')|^2 V({\bf r},{\bf r}') \; .
\eeq
We call $S_{GP}$ the Gross-Pitaevskii (GP) action of the Bose condensate.
By using the quantum least action principle
we get the Euler-Lagrange equation
\beq
i\hbar {\partial \over \partial t}\phi ({\bf r},t)=
\Big[ -{\hbar^2\over 2m} \nabla^2
+ V_0({\bf r}) + (N-1)
\int d^3{\bf r}' V({\bf r},{\bf r}') |\phi ({\bf r}',t)|^2
\Big] \phi ({\bf r},t) \; ,
\eeq
which is an integro-differential nonlinear Schr\"odinger equation.
Such equation and the effect of a finite-range interaction have
been analyzed only by few authors (see, for example, Ref. 14-16).
In fact, at low energies, it is possible to substitute
the true interaction with a pseudo-potential$^{7}$
\beq
V({\bf r},{\bf r}') = B \; \delta^3 ({\bf r}-{\bf r}') \; ,
\eeq
where $B={4\pi \hbar^2 a_s/m}$ is the scattering amplitude and $a_s$
the s-wave scattering length. In this way one obtains
the so-called time-dependent GP equation$^{5,6}$
\beq
i\hbar {\partial \over \partial t}\phi ({\bf r},t)=
\Big[ -{\hbar^2\over 2m} \nabla^2
+ V_0({\bf r}) + B (N-1) |\phi ({\bf r},t)|^2 \Big] \phi ({\bf r},t) \; ,
\eeq
that is the starting point of many calculations.$^{4}$
Note that the GP equation is accurate to describe the condensate
of weakly-interacting bosons only
near zero temperature, where thermal excitations can be neglected.
At finite temperature, one must consider a generalized GP equation plus
the Bogoliubov equations, which describe the quasi-particle
elementary excitations of the condensate.$^{7}$
\section{Trial Wave-Function of the Condensate}
\par
In this section we use a Gaussian trial wave-function to
describe the ground-state of the Bose condensate and to minimize
the GP action.
\par
We consider a triaxially asymmetric harmonic trapping potential of
the form
\begin{equation}
V_0\left({\bf r}\right)=
{1\over 2}m\omega_0^2\left(\omega_1^2 x^2+
\omega_2^2 y^2+\omega_3^2 z^2\right)\;,
\end{equation}
where $\omega_i$ ($i=1,2,3$) are adimensional constants proportional
to the spring constants of the potential along the three axes.
In the rest of this section we write the lengths in units
of the characteristic length of the
trap $a_0=\sqrt{\hbar /(m\omega_0)}$,
the time in units $\omega_0^{-1}$, the action in units $\hbar$ and
the energy in units $\hbar \omega_0$.
By using the delta pseudo-potential to model
the interaction between particles, the GP action
is explicitly given by
\beq
S_{GP} = N \int dt \; d^3{\bf r}\;
\phi^*({\bf r},t)
\Big[
i {\partial \over \partial t} +{1\over 2} \nabla^2
- V_0({\bf r}) - {1 \over 2} B (N-1) |\phi ({\bf r})|^2
\Big] \phi ({\bf r},t) \; ,
\eeq
where $B=4\pi a_s/a_0$ in our units.
\par
We want minimize this GP action $S_{GP}$ by choosing
an appropriate single-particle trial wave-function.
A good choice is the following
\beq
\phi\left({\bf r},t\right)=
\Big[ {1\over \pi^3 {\sigma}_1^2(t) {\sigma}_2^2(t)
{\sigma}_3^2(t)} \Big]^{1/4}
\prod_{i=1,2,3}
\exp\left\{-\frac{x_i^2}{2 {\sigma}_i^2(t)}
+i \beta_i(t) x_i^2 \right\}\;,
\eeq
with $(x_1,x_2,x_3)\equiv(x,y,z)$. $\sigma_i$ and $\beta_i$ are
the time-dependent variational parameters. The $\sigma_i$ are
the widths of the condensate in the three axial directions.
Note that, in order to describe the time evolution of the variational
function, the phase factor $i \beta_i(t) x_i^2$ is needed.$^{11,12}$
The choice of a Gaussian shape for the condensate
is well justified in the limit of weak interatomic coupling,
because the exact ground-state of the linear Schr\"odinger
equation with harmonic potential is a Gaussian. Moreover,
for the description of the collective
dynamics of Bose-Einstein condensates,
it has already been shown that the variational technique based
on Gaussian trial functions leads to
reliable results even in the large condensate number limit.$^{15-17}$
\par
By inserting the trial wave-function,
after spatial integration, the GP action becomes
\beq
S_{GP}= N \int dt \; L({\dot \beta_i},\beta_i, \sigma_i) \; ,
\eeq
where the effective Lagrangian is given by
\beq
L=-{1\over 2} \Big[ \sum_{i=1}^3
\big( {\dot \beta_i}\sigma_i^2 + {1\over 2 \sigma_i^2} +
2 \sigma_i^2 \beta_i^2 + {1\over 2} \omega_i^2 \sigma_i^2 \big)
+ {{\tilde g} \over \sigma_1\sigma_2\sigma_3} \Big] \; ,
\eeq
with $g=(2/\pi)^{1/2}(a_s/a_0)$ and ${\tilde g}=g(N-1)$.
Again, by imposing the least action principle, we find
six Euler-Lagrange equations
\beq
\beta_i = -{{\dot \sigma_i}\over 2 \sigma_i}
\; , \;\;\;\;\;\;\;\;
{\ddot \sigma_i}+\omega_i^2 \sigma_i = {1\over \sigma_i^3} +
{{\tilde g} \over \sigma_1 \sigma_2 \sigma_3 \sigma_i} \; ,
\eeq
with $i=1,2,3$. Note that the time dependence of
$\beta_i$ is fully determined by that of $\sigma_i$. Moreover,
for ${\tilde g}=0$ the six equations are exact solutions
of the time-dependent GP equation.
The last three differential equations are the classical
equations of motion of a system with coordinates $\sigma_i$ and
total energy per particle given by
\beq
{E\over N} = {1\over 2}
\Big[ \sum_{i=1}^3{1\over 2}\dot{\sigma}_i^2 \Big] +
U(\sigma_1,\sigma_2,\sigma_3) \; ,
\eeq
where the potential energy is given by
\beq
U(\sigma_1,\sigma_2,\sigma_3)= {1\over 2} \Big[
\sum_{i=1}^3
\big({1\over 2}\omega_i^2\sigma_i^2 + {1\over 2\sigma_i^2}\big)
+{{\tilde g}\over \sigma_1\sigma_2\sigma_3} \Big] \; .
\eeq
We use the previous equations to calculate the ground-state energy
and the low-energy excitations of the condensate.
\par
The ground-state energy per particle of the Bose
condensate is simply $E^{(0)}/N=U({\vec \sigma}^*)$, where
${\vec \sigma^* }=(\sigma_1^*,\sigma_2^*,\sigma_3^*)$
is the the minimum of the effective potential energy $U$.
This equilibrium point is identified by the following equations:
\beq
\omega_1^2\sigma_1^4 -{\tilde g}\frac{\sigma_1}{\sigma_2\sigma_3} =1\; ,
\;\;\;
\omega_2^2\sigma_2^4 -{\tilde g}\frac{\sigma_2}{\sigma_1\sigma_3} =1\; ,
\;\;\;
\omega_3^2\sigma_3^4 -{\tilde g}\frac{\sigma_3}{\sigma_1\sigma_2} =1\; .
\end{equation}
\par
The low-energy collective excitations of the condensate are
the small oscillations of variables $\sigma_i$'s around the equilibrium
point. The calculation of the normal mode frequencies for the motion of the
condensate is reduced to an eigenvalue problem for the
Hessian matrix $\Lambda$,
given by
\begin{equation}
\Lambda_{ij}=
\left.\frac{\partial^2U}{\partial\sigma_i\partial\sigma_j}
\right|_{{\vec \sigma}={\vec \sigma}^*}\; .
\end{equation}
\par
In the case ${\tilde g}=0$, the stationary solution is simply
$\sigma_i^*=\omega_i^{1/2}$ $(i=1,2,3)$ and
the energy per particle of the non-interacting condensate is given by
${E^{(0)}/N}=(\omega_1+\omega_2 + \omega_3)/2$.
The Hessian matrix $\Lambda$ of the potential energy $U$ is diagonal
and the frequency modes are simply $\Omega_i=2 \omega_i$ $(i=1,2,3)$.
\par
The general solution for stationary state
and eigenvalue problem of the collective mode frequencies can
be obtained numerically,$^{17}$ but analytical results are available
in the Thomas-Fermi limit.
In such limit, i.e. when ${\tilde g}\gg 1$,
the right-hand side in the equations of the equilibrium point,
related to the kinetic pressure, can be
neglected and the following expression for the
equilibrium point is obtained
\beq
\sigma_1^*=\Big({{\tilde g} \; \omega_2 \; \omega_3
\over \omega_1^4} \Big)^{1/5} \; ,
\;\;\;\;\;
\sigma_2^*=\Big({{\tilde g} \; \omega_1 \; \omega_3
\over \omega_2^4} \Big)^{1/5} \; ,
\;\;\;\;\;
\sigma_3^*=\Big({{\tilde g} \; \omega_1 \; \omega_2
\over \omega_3^4} \Big)^{1/5} \; .
\eeq
The energy per particle of the ground-state results
\beq
{E^{(0)}\over N} = {5\over 4}
\Big( {\tilde g} \; \omega_1 \; \omega_2 \; \omega_3 \Big)^{2/5} \; ,
\eeq
and the Hessian matrix $\Lambda$ is given by
\beq
\Lambda =
\left( \begin{array}{ccc}
3 \; \omega_1^2 & \omega_1 \omega_2 & \omega_1 \omega_3 \\
\omega_1 \omega_2 & 3 \; \omega_2^2 & \omega_2 \omega_3 \\
\omega_1 \omega_3 & \omega_2 \omega_3 & 3 \; \omega_3^2 \\
\end{array} \right) \; .
\eeq
Note that ${\tilde g}$ does not appear in the matrix.
The eigenfrequencies $\Omega$
of the collective motion are found as the solutions of the equation$^{14}$
\beq
\Omega^6
-3\left(\omega_1^2+\omega_2^2+\omega_3^2\right)\Omega^4+
8\left(\omega_1^2\; \omega_2^2+\omega_1^2\; \omega_3^2
+\omega_2^2\; \omega_3^2 \right)\Omega^2
-20\; \omega_1^2\; \omega_2^2\; \omega_3^2=0 \; .
\eeq
This equation has been also obtained by Stringari$^{4}$ using a hydrodynamical
approach, which is effective in the Thomas-Fermi limit.
\par
For an axially symmetric trap, where
$\omega_1=\omega_2=\omega_{\bot} \neq \omega_3$,
the previous equations give
$\sigma_1^* =\sigma_2^* =({\tilde g}\omega_3/\omega_{\bot}^3)^{1/5}$,
$\sigma_3^*=({\tilde g}\omega_{\bot}/\omega_3^3)^{1/5}$,
$E^{(0)}/N=(5/4){\tilde g}^{2/5} \omega_{\bot}^{4/5}\omega_3^{2/5}$, and
\beq
\Omega_{1,2}=\sqrt{2\omega_{\bot}^2 + {3\over 2}\omega_3^2 \pm
{1\over 2}\sqrt{16 \omega_{\bot}^4 + 9\omega_3^4 -
16\omega_{\bot}^2\omega_3^2}} \; ,
\;\;\;\;
\Omega_3 =\sqrt{2} \; \omega_{\bot} \; .
\eeq
Because of the axial symmetry of the trap, the angular momentum
along the $z$ axis is conserved and one can thus label the modes
by angular ($l$) and azimuthal ($m_z$) quantum numbers.
One finds$^{4,11}$ $m_z=0$ for the coupled monopole ($l=0$)
and quadrupole ($l=2$) modes
$\Omega_{1,2}$, and $|m_z|=2$ for the quadrupole ($l=2$)
mode $\Omega_3$. Note that the experimental results obtained
on sodium vapors at MIT are in excellent agreement with these
theoretical values.$^{18}$
\par
For an isotropic harmonic trap with frequency $\omega$, we obtain
$\sigma_i^* =({\tilde g}/\omega^2)^{1/5}$ ($i=1,2,3$),
$E^{(0)}/N=(5/4){\tilde g}^{2/5} \omega^{6/5}$, and
$\Omega_{1,2}=\sqrt{5}\; \omega$, $\Omega_3 = \sqrt{2} \; \omega$.
Here, differently from the axially symmetric case, the frequencies
do not depend on $m_z$. In particular,
$\sqrt{5}\omega$ is the monopole oscillation, also called
breathing mode, characterized by radial quantum number
$n_r=1$ and angular quantum number $l=0$.
Instead $\sqrt{2}\omega$ is the quadrupole ($l=2$) surface ($n_r=0$)
oscillation.
\par
It is worth noting that the variational formalism
makes it possible to investigate the case of negative scattering length.
In this case, the large $N$ limit of the Thomas-Fermi approximation
does not apply because
the condensate number cannot exceed a critical value.
Let us consider again the isotropic trap with frequency
$\omega$. For the ground-state, the number of particles
satisfies the following equation
\beq
N = 1+ {1\over g} (\omega^2 \sigma_0^5-\sigma_0) \; ,
\eeq
where $\sigma_0 =\sigma_i^*$ ($i=1,2,3$) is mean radius of the
condensate, that is the same in the three axial directions.
It follows that, for $g>0$, with $N=1$ we have $\sigma_0 =1$
and with $N\to \infty$ then $\sigma_0\to \infty$.
Instead, for $g<0$, the mean radius $\sigma_0$ decreases
by increasing the number of bosons
$N$ to a critical radius $\sigma_0^c =(5\omega^2)^{-1/4}$
with a critical number of bosons given by
$N^c=1+|g|^{-1}\omega^{-1/2}(5^{-1/4}-5^{-5/4})$.
Moreover, one finds a simple expression for the monopole ($l=0$)
frequency of surface oscillation of the condensate, namely$^{16}$
\beq
\Omega = \Big[5\omega^2 -\sigma_0^{-4}\Big]^{1/2} \; .
\eeq
We have seen that $\sigma_0=1$ for $N=1$ or $g=0$. Moreover,
when $N>1$, it is $\sigma_0>0$ for $g>0$
and $\sigma_0^c<\sigma < 1$ for $g<0$.
Regarding the frequency $\Omega$,
one easily verifies that $\Omega \to \sqrt{5}\omega$ for
$\sigma_0 \to \infty$ and $\Omega \to 0$ for $\sigma_0\to \sigma_0^c$.
Note that, for the triaxial trap, at the critical point
the two highest frequencies diverge and the
lowest one falls to zero.$^{17}$
\par
It is important to observe that it is possible to include
in the Gaussian trial wave-function other three parameters
describing the position of the center of mass of the condensate.
Due to the harmonic confinement, the motion of the center
of mass (dipole mode) is periodic and the frequencies of
oscillation are equal to those of the harmonic trap.$^{11,12}$
\section{Trial Wave-Function of Vortex States}
\par
Let us consider states having a vortex line
along the $z$ axis and all bosons flowing around it with
quantized circulation. The observation of these states would be
a signature of macroscopic phase coherence of trapped BEC.
The time-dependent variational approach can be used
to describe also these vortex states of the condensate in a
axially symmetric trap
\beq
V_0({\bf r})={1\over 2} \left( \omega_{\bot}^2 r^2 +
\omega_3^2 z^2 \right) \; ,
\eeq
where $r=\sqrt{x^2+y^2}$, $\theta=$arctg$(y/x)$ and $z$ are the
cylindrical coordinates.
\par
The trial wave function of the vortex states of the condensate
can be chosen as
\beq
\phi_k \left({\bf r},t \right)=
\Big[{1\over k! \pi^3 {\sigma_{\bot}}^{4k+4}(t)
{\sigma}_3^2(t) }\Big]^{1/4}
r^k e^{i k \theta}
\exp\left\{-\frac{r^2}{2 \sigma_{\bot}^2(t)}-{z^2\over 2 \sigma_3^2(t)}
+i \beta_{\bot}(t) r^2 + i \beta_3(t) z^2 \right\}\;,
\eeq
where $k$ is the vortex quantum number and
$\sigma_{\bot}$, $\sigma_3$, $\beta_{\bot}$ and $\beta_3$ are
the time-dependent variational parameters.
This trial wave-function describes exactly the
vortex states of a non-interacting gas of bosons in the harmonic trap.
\par
By following the same procedure of the previous section,
namely by inserting the trial wave-function in the GP action,
and after spatial integration, one gets
\beq
S_{GP}= N \int dt \; L({\dot \beta_{\bot}},{\dot \beta_3},
\beta_{\bot} ,\beta_3, \sigma_{\bot}, \sigma_3) \; ,
\eeq
where the effective Lagrangian is given by
$$
L=-{1\over 2} \Big[ (k+1)\big( {\dot \beta_{\bot}}\sigma_{\bot}^2 +
{1\over 2 \sigma_{\bot}^2} + 2 \sigma_{\bot}^2 \beta^2 + {1\over 2}
\omega_{\bot}^2 \sigma^2 \big)+
$$
\beq
+\big( {\dot \beta_3}\sigma_3^2 + {1\over 2 \sigma_3^2} +
2 \sigma_3^2 \beta_3^2 + {1\over 2}
\omega_3^2 \sigma_3^2 \big)
+ {{\tilde g}(k) \over \sigma_{\bot}^2\sigma_3} \Big] \; ,
\eeq
with $g=(2/\pi)^{1/2}(a_s/a_0)$ and
${\tilde g}(k)=g(N-1)(2k)!/(2^{2k}k!^2)$.
The four Euler-Lagrange equations of the system are given by
\beq
\beta_{\bot} = -{{\dot \sigma_{\bot}}\over 2 \sigma_{\bot}}
\; , \;\;\;\;\;\;\;\;\;\;
\beta_3 = -{{\dot \sigma_3}\over 2 \sigma_3} \; ,
\eeq
\beq
(k+1){\ddot \sigma_{\bot}}+(k+1)\omega_{\bot}^2 \sigma_{\bot}
= {(k+1)\over \sigma_{\bot}^3} +
{{\tilde g}(k) \over \sigma_{\bot}^3 \sigma_3} \; ,
\;\;\; \;\;\;
{\ddot \sigma_3}+\omega_3^2 \sigma_3 = {1\over \sigma_3^3} +
{{\tilde g}(k) \over \sigma_{\bot}^2 \sigma_3^2} \; .
\eeq
The time dependence of $\beta_{\bot}$ and $\beta_3$
is fully determined by that of $\sigma_{\bot}$ and $\sigma_3$;
moreover, for ${\tilde g}=0$ (non-interacting condensed vortex)
these four equations are exact.
The last two differential equations correspond to the classical
equations of motion of a system with energy per particle
\beq
{E_k\over N} = {1\over 2} \Big[
(k+1)\dot{\sigma}_{\bot}^2 + {1\over 2}\dot{\sigma}_3^3 \Big]
+ U(\sigma_{\bot},\sigma_3) \; ,
\eeq
where the potential energy is given by
\beq
U(\sigma_{\bot},\sigma_3)= {1\over 2} \Big[
(k+1)\omega_{\bot}^2\sigma_{\bot}^2 +
{1\over 2}\omega_3^2\sigma_3^2 +
{(k+1)\over \sigma_{\bot}^2}+ {1\over 2\sigma_3^2}
+{{\tilde g}(k)\over \sigma_{\bot}^2 \sigma_3} \Big] \; .
\eeq
The equilibrium point $(\sigma_{\bot}^*,\sigma_3^*)$,
corresponding to the minimum
of the effective potential energy, is given by the following equations:
\beq
(k+1)\omega_{\bot}^2\sigma_{\bot}^4
-{\tilde g}(k){1\over \sigma_3} =(k+1)\; ,
\;\;\; \;\;\;\;
\omega_3^2\sigma_3^4 -{\tilde g}(k)\frac{\sigma_3}{\sigma_{\bot}^2} =1\; .
\eeq
\par
In the non-interacting case $({\tilde g}=0)$
the stationary solution is simply
$\sigma_{\bot}^*=\omega_{\bot}^{1/2}$ and
$\sigma_3^*=\omega_3^{1/2}$.
The stationary energy of the non-interacting vortex state
is exactly given by
\beq
{E_k^{(0)}\over N} = (k+1)\omega_{\bot}+{1\over 2}\omega_3 \; .
\eeq
It is important to observe that the energy grows only
linearly with the vortex quantum number $k$, in contrast
to the homogenous case$^{7}$ where the vortex energy grows as $k^2$.
In addition, the two ($m_z=0$) frequency modes are
$\Omega_1=2 \omega_{\bot}$ and $\Omega_2=2 \omega_3$
and do not depend on the vortex quantum number $k$.
\par
In the Thomas-Fermi limit (${\tilde g}>>1$), the analytic
expression of the equilibrium point reads
\beq
\sigma_{\bot}^*=\left[ {\tilde g}(k) \; \omega_3
\over (k+1)^{3/2}\omega_{\bot}^3\right]^{1/5} \; , \;\;\; \;\;\;
\sigma_3^*=\left[{{\tilde g}(k) \; (k+1)\; \omega_{\bot}^2
\over \omega_3^4} \right]^{1/5} \; .
\eeq
The stationary energy of the vortex state results
\beq
{E_k^{(0)}\over N}= {5\over 4} \Big( (k+1) \; {\tilde g}(k)
\; \omega_{\bot}^2 \;
\omega_3 \Big)^{2/5} \; .
\eeq
>From this formula it is easy to calculate the critical frequency
${\bar \Omega}_c$ at which a vortex can be produced.
One has to compare the energy of a vortex state in frame rotating
with angular frequency ${\bar \Omega}$, that is
$E^{(0)}_k-{\bar \Omega} L_z$, with
the energy $E^0$ of the ground state (with no vortices).$^{7}$
Since the angular momentum per particle is $\hbar k$, the critical
frequency is given by $\hbar {\bar \Omega}_c =(E_k^{(0)}/N - E^{(0)}/N)/k$.
By using our previous formula, we get
\beq
{\bar \Omega}_c = {5\over 4} k^{-1} \omega_{\bot}^{4/5} \omega_3^{2/5}
\Big[ {\tilde g}(k)^{2/5} (k+1)^{2/5} - {\tilde g}(0)^{2/5} \Big] \; .
\eeq
In the presence of a vortex,
for frequencies smaller than the critical frequency ${\bar \Omega}_c$,
the system is not thermodynamically
stable (global stability), but can be locally stable (metastability)
against small deformations of the system. Rohksar$^{19}$
has shown that the vortex state is locally unstable
for ${\bar \Omega}=0$. Actually, in a harmonic trap,
there is local instability for $0\leq {\bar \Omega} < {\bar \Omega}_s$, where
the frequency ${\bar \Omega}_s$ (${\bar \Omega}_s < {\bar \Omega}_c$)
can be estimated numerically.$^{20}$
On the other hand, we have recently shown that,
in toroidal traps, vortices can be locally stable also
in absence of external forced rotation.$^{21}$
\par
It is important to observe that, in our present treatment,
the condensate has a unique vortex line, which is fixed along
the $z$ axis. Butts and Rokhsar$^{22}$ have numerically shown that
increasing the rotational frequency ${\bar \Omega}$ well beyond
${\bar \Omega}_c$ will favour the creation of more complex vortical
configurations, associated with the occurrence of $2$ or more vortex lines.
\par
To calculate the collective oscillations
in the Thomas-Fermi limit, we follow the same
procedure of the previous section. The Hessian matrix $\Lambda$ of
the potential energy $U$ at the equilibrium point is given by
\beq
\Lambda =
\left( \begin{array}{cc}
8\; (k+1)\; \omega_{\bot}^2 & 2\; (k+1)^{1/2}\; \omega_{\bot}\; \omega_3 \\
2\; (k+1)^{1/2}\; \omega_{\bot}\; \omega_3 & 3\; \omega_3^2 \\
\end{array} \right) \; .
\eeq
Again the ${\tilde g}(k)$ dependence drops out.
Moreover, the mass matrix of the kinetic energy is not the identity
matrix but instead
\beq
M=
\left( \begin{array}{cc}
2 \; (k+1) & 0 \\
0 & 1 \\
\end{array} \right) \; .
\eeq
The eigenfrequencies $\Omega$
of the collective motion are found as the solutions of the equation
$|\Lambda - \Omega^2 M| = 0$, which gives the $m_z=0$ frequencies
\beq
\Omega_{1,2}=\sqrt{2\omega_{\bot}^2 + {3\over 2}\omega_3^2 \pm
{1\over 2}\sqrt{16 \omega_{\bot}^4 + 9\omega_3^4 -
16\omega_{\bot}^2\omega_3^2}} \; .
\eeq
This is an interesting result: The $m_z=0$ frequencies
of a harmonically trapped condensate do not depend on
the vortex quantum number $k$. Note that Zambelli and Stringari$^{23}$
have shown, by using a sum-rule approach, that the
$m_z\neq 0$ collective frequencies do instead depend on $k$.
Nevertheless, all the modes are $k$ sensitive because of
the amplitudes, which are given by $^{1/2}=(k+1)^{1/2}\sigma_{\bot}$
and $^{1/2}=\sigma_3/\sqrt{2}$.
\par
Finally, we study an isotropic trap with
frequency $\omega$. We have $\sigma^*_{\bot}=\sigma^*_3$
only for ${\tilde g}=0$. In order to obtain some analytic
results also in the case of negative scattering length,
we consider the vortex trial wave-function with
a unique variational parameter $\sigma$, i.e.
we put $\sigma_{\bot}=\sigma_3=\sigma$.
For the ground-state, the number of particles
satisfies the equation $N = 1+g(k)^{-1}
(\omega^2\sigma_0^5-\sigma_0) (1+{2/3}k)$,
where $\sigma_0$ is mean radius of the condensate in the
stationary configuration and $g(k)=g(2k)!/(2^{2k}k!^2)$.
It follows that for $g<0$, the critical number of bosons
for the collapse of the wave-function is given
by $N^c=1+|g(k)|^{-1}\omega^{-1/2}(5^{-1/4}-5^{-5/4})(1+2/3k)$.
Thus, compared with the condensed state without vortices,
the critical number of trapped atoms in the vortex state increases.
This variational result is in full agreement
with the $k=1$ numerical results of Dalfovo and Stringari.$^{24}$
\section{Conclusions}
We have studied a Bose-Einstein condensate trapped
in a harmonic potential by using a variational method based on
the minimization of the Gross-Pitaevskii action with a Gaussian ansatz
for the shape of the macroscopic wave function of the
condensate. Our analysis has been applied to the whole range of
coupling strengths between particles (attractive and
repulsive interaction) to study both the ground-state
and vortex states. In particular, for vortex states,
we have determined the critical rotational
frequency at which a single vortex becomes globally stable.
Moreover, we have calculated the collective oscillations
and the critical number of bosons for the collapse of the condensate
in the case of attractive interaction.
Nowadays it is experimentally possible to trap and condense a very large
number of alkali atoms (more than $10^{9}$) and to obtain
a less diluite system. Thus, the study
of the time-dependent variational approach
without the independent-particle approximation
begins to appear an interesting problem also for Bose gases
in an external potential (the case of a superfluid film containing a vortex
has been discussed in Ref. 9 and 10).
This will be one of our future projects.
\section*{Acknowledgments}
This work has been supported by INFM through a Research Advanced Project
on Bose-Einstein Condensation. The author thanks L. Reatto,
A. Parola and D. Pini for useful discussions.
\newpage
\section*{References}
\begin{description}
\item{\ 1.} M.H. Anderson {\it et al.}, Science {\bf 269}, 198 (1995).
\item{\ 2.} C.C. Bradley {\it et al.},
Phys. Rev. Lett. {\bf 75}, 1687 (1995).
\item{\ 3.} K.B. Davis {\it et al.},
Phys. Rev. Lett. {\bf 75}, 3969 (1995).
\item{\ 4.} F. Dalfovo, S. Giorgini, L.P. Pitaevskii, and S. Stringari,
Rev. Mod. Phys. {\bf 71}, N.3 (1999).
\item{\ 5.} E.P. Gross, Nuovo Cimento {\bf 20}, 454 (1961).
\item{\ 6.} L.P. Pitaevskii, Zh. Eksp. Teor. Fiz. {\bf 40},
646 (1961) [English Transl. Sov. Phys. JETP {\bf 13}, 451 (1961)].
\item{\ 7.} A. Fetter and J. Walecka, {\it Quantum Theory
of Many Particle Systems} (McGraw-Hill, New York, 1971).
\item{\ 8.} G.D. Mahan, {\it Many-Particle Physics}
(Plenum Press, New York, 1990).
\item{\ 9.} Q. Niu, P. Ao, and D.J. Thouless,
Phys. Rev. Lett. {\bf 72}, 1706 (1994).
\item{\ 10.} D. Demircan, P. Ao, and Q. Niu,
Phys. Rev. B {\bf 54}, 10027 (1996).
\item{\ 11.} V.M. P\'erez-Garcia, H. Michinel, J.I. Cirac, M. Lewenstein,
and P. Zoller, Phys. Rev. Lett. {\bf 77}, 5320 (1996)
\item{\ 12.} V.M. P\'erez-Garcia, H. Michinel, J.I. Cirac, M. Lewenstein,
and P. Zoller, Phys. Rev. A {\bf 56} 1424 (1997).
\item{\ 13.} L. Salasnich, Mod. Phys. Lett. B {\bf 11}, 1249 (1997).
\item{\ 14.} L. Salasnich, Mod. Phys. Lett. B {\bf 12}, 649 (1998).
\item{\ 15.} A. Parola, L. Salasnich, and L. Reatto,
Phys. Rev. A {\bf 57}, R3180 (1998).
\item{\ 16.} L. Reatto, A. Parola, and L. Salasnich,
J. Low Temp. Phys. {\bf 113}, 195 (1998).
\item{\ 17.} E. Cerboneschi, R. Mannella, E. Arimondo, and L. Salasnich,
Phys. Lett. A {\bf 249}, 495 (1998).
\item{\ 18.} D.M. Stamper-Kurn {\it et al.}, cond-mat/9801262.
\item{\ 19.} D.S. Rokhsar, Phys. Rev. Lett. {\bf 79}, 2164 (1997).
\item{\ 20.} A.A. Svidzinsky and A.L. Fetter, cond-mat/9811348.
\item{\ 21.} L. Salasnich, A. Parola, and L. Reatto,
Phys. Rev. A {\bf 59}, 2990 (1999).
\item{\ 22.} D. A. Butts and R. Rokhsar, Nature {\bf 397}, 327 (1999).
\item{\ 23.} Zambelli and Stringari,
Phys. Rev. Lett. {\bf 81}, 1754 (1998).
\item{\ 24.} F. Dalfovo and S. Stringari,
Phys. Rev. A {\bf 53}, 2477 (1996).
\end{description}
\end{document}