# Ground state and vortex states of bosons in an anisotropic trap: A variational approach

## I Introduction

Observation of Bose-Einstein condensation in cooled and trapped dilute gases of alkali atoms  and spin polarized atomic hydrogen has generated a renewed theoretical interest in understanding such systems. In a meanfield approach, which is valid in the limit $`\rho a^3<<1`$, where $`\rho `$ is density of atoms and $`a`$ is the $`s`$-wave scattering length, ground state and vortex states of these systems can be described by Gross-Pitaevskii (GP) equation. Various numerical procedures and approximate analytical methods have been used to solve the GP equation. Among these variational scheme proposed by Baym and Pethick  to explain the experimental observations of Ref. is particularly appealing. In this approach the trial wave function was taken to be of the form of ground state of the trap potential (modeled by an anisotropic harmonic oscillator potential). Thus the wave function is represented by a three dimensional Gaussian with axial and transverse frequencies as variational parameters. This form of wave function, however, is valid only when the number of atoms in the trap is very small. As the number increases, the repulsive interaction between the atoms tends to expand the condensate and flatten the density profile in the central region of the trap where the density is maximum. Of these two effects, only expansion of the condensate can be described adequately by the Gaussian trial wave function. On the other hand, the Thomas-Fermi approximation  provides a wave function which is valid when the number of atoms is very large which is the case with the recent experiments . It is of interest, however, to have a form of wave function which apart from showing these limiting behaviours, is also valid in the intermediate region. In a variational approach, Fetter  has proposed such a wave function for the bosons in an isotropic trap.

In this paper we propose a simple form of wave function for the ground state of bosons confined in an anisotropic trap. The wave function is valid for a very wide range of particle numbers. When the number is small it tends to mimic a Gaussian, and in the opposite limit it resembles the Thomas-Fermi wave function. However for large number of atoms the wave function differs from the Thomas-Fermi wave function in the surface region, a desirable feature as noted in Ref.. The trial wave function has got an additional parameter comapared to the ones used in Ref.. This parameter takes care of the flattening of the density in the central region of the condensate. Thus, it provides a better lower bound for the ground state energy than the Gaussian trial wave function. We also compare the results obtained by our trial wave function with those obtained by other numerical procedures such as the minimization of energy functional by steepest descent method  and the integration of the nonlinear Schrodinger equation . These comparisons show the form chosen by us to be highly accurate for obtaining a host of physical qauantities of interest. In addition to providing acurate results for a very wide range of the particle numbers, the method is also very fast from the computational point of view. Further, using this wave function, the physical observables can be expressed analytically in terms of three variational parameters which are obtained by minimizing the GP energy functional. Since we have a simple analytical form for the energy in terms of three variational parameters the procedure of minimization is very simple. The novel achievement of this method lies in verfication of the predicted behaviour of the aspect ratio which is very important quantity from the experimental point of view. It could not be ascertained before because convergence of the aspect ratio to the highly repulsive limit is very slow . Since our method can handle even very large number of atoms in the trap without any difficulty we could verify the behaviour of the aspect ratio.

Based on physical reasoning we generalise the variational form to descibe the vortex states also. As is the case with the ground state, we find good agreement with the existing results, with considerably less computational effort.

Condensation has also been observed in <sup>7</sup>Li where interatomic interaction is attractive which is characterized by the negative s-wave scattering length. As the number of atoms in the trap inceases the condensate shrinks and nonuniformity in the central region increases. After a critical number of atoms in the trap, the condensate collapses. This situation is also very well descibed by our trial wave function. The parameter which accounts for flattening of the density profile in the case of repulsive interactions also takes care of the increase in density gradient in an effective way . The critical number of atoms for the case of <sup>7</sup>Li is in close agreement with that reported in Ref..

The paper is organized as follows. Section II contains the description of the variational scheme employed in the paper. It also contains the analytic expressions for the observables of interest. Results obtained from the variational procedure and their comparison with the existing ones are presented in Section III. Section IV contains summary and concluding remarks.

## II Variational Method

### A Ground state

Bose Einstein condensation in experiments with cooled and trapped atoms can be described within the framework of the GP theory. Validity of such a description has been analysed by Stenholm. In a situation where the trap can be modeled by an anisotropic harmonic oscillator potential characterized by the two angular frequencies $`\omega _{}^0`$ and $`\omega _z^0`$, and the interatomic ineractions can be replaced by an effective pseudo-potential involving $`s`$-wave scattering length $`a`$, the ground state energy for condensed bosons of mass $`m`$ is given by the Gross-Pitaevskii functional

$$\frac{E_1[\psi ]}{N}=𝑑𝐫_1\frac{1}{2}\left[|_1\psi _1(𝐫_\mathrm{𝟏})|^2+\left(x_1^2+y_1^2+\lambda _0^2z_1^2\right)|\psi (𝐫_1)|^2+\frac{u_1}{2}|\psi _1(𝐫_1)|^4\right].$$
(1)   
Here we have used the length scale $`a_{}=\sqrt{\mathrm{}/m\omega _{}^0}`$ and the energy scale $`\mathrm{}\omega _{}^0`$ provided by the trap potential to express Eq.(1) in terms of the dimensionless variables. $`\psi _1(𝐫_1)`$ is the condensate wave function which satisfies the normalization condition

$$𝑑𝐫_1|\psi _1(𝐫_1)|^2=1.$$
(2)   
$`\lambda _0=\omega _z^0/\omega _{}^0`$ is the anisotropy parameter of the trap, $`N`$ is the total number of atoms in the condensate and $`u_1=8\pi aN/a_{}`$. The exact form of the wave function can be determined by minimizing the energy functional in Eq.(1) with the normalization constraint of Eq.(2). Such a minimization results in the nonlinear Schrodinger equation

$$\left[_1^2+\left(x_1^2+y_1^2+\lambda _0^2z_1^2\right)+u_1|\psi _1(𝐫_1)|^2\right]\psi _1(𝐫_1)=2\mu _1\psi _1(𝐫_1),$$
(3)   
where, $`\mu _1`$ is the chemical potential. It is not possible to find exact analytic solution to Eq.(3). Consequently various numerical techniques have been developed to study the ground state property of such systems within the framework of the GP theory. These techniques involve either the direct numerical minimization of Eq.(1) with the constraint of Eq.(2)  or numerical integration of Eq.(3) or its time dependent version. Another approach is to use the variational method which has been extensively used in different branches of physics. The main advantage of this method is that with a suitable guess for the form of the wave function it is possible to save a lot of computational effort and time. In addition, it may also provide physical insights which generally get obscured in the complicated computational procedures. The first study of this kind was done by Baym and Pethick  in light of the experimental observations in <sup>87</sup>Rb . They took the trial wave function for the ground state as

$$\psi (𝐫)=N^{1/2}\omega _{}^{1/2}\omega _z^{1/4}\left(\frac{m}{\pi \mathrm{}}\right)^{3/4}e^{m(\omega _{}r_{}^2+\omega _zz^2)/2\mathrm{}}$$
(4)   
with effective frequencies, $`\omega _{}`$ and $`\omega _z`$, treated as variational parameters. However, the wave function above brings out only the qualitative features of the condensate e.g. expansion of the condensate in different directions, shifts in the angular frequencies and the scaling behaviour of energy with the number of atoms in the trap. Further, this form of the wave function is valid only for small number of atoms in the trap (see Fig.2 below). We now propose a variational form of the wave function and demonstrate its applicability and utility in providing accurate description of the condensate for a wide range of the particle numbers. The form of the trial wave function we choose is

$$\psi _1(𝐫_1)=\sqrt{\frac{p}{2\pi \mathrm{\Gamma }(\frac{3}{2p})}}\lambda ^{\frac{1}{4}}\left(\frac{\omega _{}}{\omega _{}^0}\right)^{\frac{3}{4}}e^{\frac{1}{2}\left(\omega _{}/\omega _{}^0\right)^p\left(r_1^2+\lambda z_1^2\right)^p},$$
(5)   
where, $`\lambda ,\omega _{}`$ and $`p`$ are the variational parameters which are obtained by minimizing the energy $`E_1`$ in Eq.(1) with respect to these parameters. It is easily verified that the wave function satisfies the normalization condition of Eq.(2). The expression of the ground state energy $`E_1`$ in terms of $`\lambda `$, $`\omega _{}`$ and $`p`$ is

$`E_1`$ $`=`$ $`{\displaystyle \frac{1}{12}}{\displaystyle \frac{\omega _{}}{\omega _{}^0}}\left(1+{\displaystyle \frac{\lambda }{2}}\right){\displaystyle \frac{\mathrm{\Gamma }\left(\frac{1}{2p}\right)}{\mathrm{\Gamma }\left(\frac{3}{2p}\right)}}(1+2p)+{\displaystyle \frac{1}{3}}{\displaystyle \frac{\omega _{}^0}{\omega _{}}}\left(1+{\displaystyle \frac{\lambda _0^2}{2\lambda }}\right){\displaystyle \frac{\mathrm{\Gamma }\left(\frac{5}{2p}\right)}{\mathrm{\Gamma }\left(\frac{3}{2p}\right)}}`$ (7)   
$`+N{\displaystyle \frac{a}{a_{}}}\left({\displaystyle \frac{\omega _{}}{\omega _{}^0}}\right)^{\frac{3}{2}}{\displaystyle \frac{\sqrt{\lambda }p}{\mathrm{\Gamma }\left(\frac{3}{2p}\right)}}\left({\displaystyle \frac{1}{2}}\right)^{\frac{3}{2p}}.`$   
For a particular value of $`N`$ the parameters $`\omega _{}`$, $`\lambda `$ and $`p`$ are obtained by minimizing the energy above using standard numerical routines. We have used Mathematica for this and it takes a few seconds of the real time on 166 MHz Pentium-1 computer to get the answer. Next we discuss how different physical observables can be obtained in terms of the parameters of the wave function. Aspect ratio which characterizes the anisotropy of the velocity distribution of the condensate is defined as $`\sqrt{<p_z^2>/<p_x^2>}`$. This can be easily shown to be

$$\sqrt{\frac{<p_z^2>}{<p_x^2>}}=\sqrt{\frac{<x_1^2>}{<z_1^2>}}=\sqrt{\lambda }.$$
(8)   
Width of the condensate in the transervse direction is given by

$$<x_1^2>=\frac{\omega _{}^0\mathrm{\Gamma }\left(\frac{5}{2p}\right)}{3\omega _{}\mathrm{\Gamma }\left(\frac{3}{2p}\right)}$$
(9)   
and the width of the the momentum distribution in this direction is given by

$$<p_x^2>=\frac{N\mathrm{}m\omega _{}\mathrm{\Gamma }\left(\frac{1}{2p}\right)(1+2p)}{12\mathrm{\Gamma }\left(\frac{3}{2p}\right)}.$$
(10)   
The peak density of the condensate is $`N\psi _1^2(0)/a_{}^3`$. Life time of the condensate is related to the density distribution. The loss rate due to the two body loss rate and the three body loss rate  is given by

$`R(N)`$ $`=`$ $`\alpha {\displaystyle 𝑑𝐫|\psi (𝐫)|^4}+L{\displaystyle 𝑑𝐫|\psi (𝐫)|^6}`$ (11)   
$`=`$ $`{\displaystyle \frac{\alpha N^2\sqrt{\lambda }(\omega _{}/\omega _{}^0)^{3/2}p}{2\pi 2^{3/2p}a_{}^3\mathrm{\Gamma }(3/2p)}}+{\displaystyle \frac{LN^3\lambda (\omega _{}/\omega _{}^0)^3p^2}{4\pi ^23^{3/2p}a_{}^6\mathrm{\Gamma }^2(3/2p)}},`$ (12)   
where, $`\alpha `$ is the two-body dipolar relaxation loss rate coefficient and $`L`$ is the three-body recombination loss rate coefficient.

### B Vortex States

We consider here the states having a vortex line along the $`z`$ axis. Wave function of such a state can be written as

$$\mathrm{\Psi }(𝐫)=\psi (𝐫)e^{ı\kappa \varphi }$$
(13)   
where $`\kappa `$ is an integer denoting the quantum of circulation. Subsituting the complex wave function $`\mathrm{\Psi }`$ in place of $`\psi `$ in Eq.(1) we get Gross-Pitaevskii functional for the vortex states in terms of the scaled variables

$$\frac{E_1[\psi ]}{N}=𝑑𝐫_1\frac{1}{2}\left[|_1\psi _1(𝐫_\mathrm{𝟏})|^2+\left(\kappa ^2r_1^2+r_1^2+\lambda _0^2z_1^2\right)|\psi (𝐫_1)|^2+\frac{u_1}{2}|\psi _1(𝐫_1)|^4\right].$$
(14)   
The corressponding nonlinear Schrodinger equation is

$$\left[_1^2+\kappa ^2r_1^2+r_1^2+\lambda _0^2z_1^2+u_1|\psi _1(𝐫_1)|^2\right]\psi _1(𝐫_1)=2\mu _1\psi _1(𝐫_1).$$
(15)   
We assume the trial wave function of the form

$$\psi _1(𝐫_1)=Ar_1^qe^{\frac{1}{2}\left(\omega _{}/\omega _{}^0\right)^p\left(r_1^2+\lambda z_1^2\right)^p}$$
(16)   
where $`q`$ is an additional variational parameter. This particular form of the wave function is motivated by the following considerations.

1. Presence of the centrifugal term $`\kappa ^2/r_1^2`$ forces the wave function to vanish on the $`z`$ axis.
2. It has been shown that for a weakly interacting Bose gas  the wave function corresponding to $`k`$th quantum circulation behaves as

$$\psi r_1^k$$
(17)   
near the $`z`$ axis.

Proportionality factor in Eq.(16) is determined by the normalization condition (Eq.2)

$$A^2=\frac{\sqrt{\lambda }p\mathrm{\Gamma }\left(\frac{3}{2}+q\right)}{\pi ^{3/2}\mathrm{\Gamma }(1+q)\mathrm{\Gamma }\left(\frac{3+2q}{2p}\right)}\left(\frac{\omega _{}}{\omega _{}^0}\right)^{(3+2q)/2}.$$
(18)   
For a vortex line descibed by the wave function in Eq.(16) the density peaks at

$$r_1=\sqrt{\frac{\omega _{}^0}{\omega _{}}}\left(\frac{q}{p}\right)^{\frac{1}{2p}},$$
(19)   
and the peak density is given by

$$\rho _{max}=\frac{N}{a_{}^3}A^2\left(\frac{\omega _{}^0}{\omega _{}}\right)^q\left(\frac{q}{p}\right)^{\frac{q}{p}}e^{q/p}.$$
(20)   
It is also sraightforward to get the analytic expression for $`E_1`$ in terms of the variational parameters $`\omega `$, $`\lambda `$, $`p`$ and $`q`$ which, in turn, are obtained by minimization of $`E_1`$. The kinetic energy is given by

$$(E_1/N)_{kin}=\frac{\omega _{}(1+2q)\left[(1+2p)(1+\lambda /2)+q(2p+2q+\lambda )\right]\mathrm{\Gamma }\left(\frac{1+2q}{2p}\right)}{4(3+2q)\omega _{}^0\mathrm{\Gamma }\left(\frac{3+2q}{2p}\right)}.$$
(21)   
Energy corresponding to the rotational motion is given by

$$(E_1/N)_{rot}=\frac{\kappa ^2\omega _{}(1+2q)\mathrm{\Gamma }\left(\frac{1+2q}{2p}\right)}{4\omega _{}^0q\mathrm{\Gamma }\left(\frac{3+2q}{2p}\right)}.$$
(22)   
The oscillator energy is given by

$$(E_1/N)_{HO}=\frac{\omega _{}^0\lambda \left(1+q+\frac{\lambda _0^2}{2\lambda }\right)\mathrm{\Gamma }\left(\frac{5+2q}{2p}\right)}{\omega _{}\left(3+2q\right)\mathrm{\Gamma }\left(\frac{3+2q}{2p}\right)}.$$
(23)   
The interaction energy is given by

$$(E_1/N)_{pot}=2\left(\frac{\omega _{}}{\omega _{}^0}\right)^{3/2}\frac{a\sqrt{\lambda }Np(1+2q)^2\mathrm{\Gamma }(2q)\mathrm{\Gamma }^2\left(\frac{1}{2}+q\right)\mathrm{\Gamma }\left(\frac{3+4q}{2p}\right)}{2^{(3+4q)/2p}\pi ^{1/2}a_{}q(1+4q)\mathrm{\Gamma }^2(q)\mathrm{\Gamma }\left(\frac{1}{2}+2q\right)\mathrm{\Gamma }^2\left(\frac{3+2q}{2p}\right)}.$$
(24)   
It is easy to verify that ground state is obtained by setting $`\kappa `$ and $`q`$ equal to zero. Once we have the energy of the states with and without vortices we can calculate the critical angular velocity for the formation of the vortex state. In the unit of $`\omega _{}^0`$ it is given by 

$$\mathrm{\Omega }_c=\kappa ^1\left[(E_1/N)_\kappa (E_1/N)_0\right].$$
(25)   
To demonstrate the applicability of this method we have performed calculations for <sup>87</sup>Rb and <sup>7</sup>Li. The $`s`$ wave scattering length is positive for <sup>87</sup>Rb. It is negative for <sup>7</sup>Li. Consequently the interatomic interaction is repulsive in the former and attractive in the latter. We now present the results and their comparison with the existing calculations.

## III Results

### A Positive scattering length: <sup>87</sup>Rb

In this section we report calculations on <sup>87</sup>Rb. We perform our calculations by employing the experimental numbers for the asymmetry parameter of the trap, the axial frequency and the $`s`$ wave scattering length corresponding to the experimental situation of Ref. and the subsequent theoretical calculations . Accordingly, $`\lambda ^0=\omega _z^0/\omega _{}^0=\sqrt{8}`$; $`\omega _z^0/2\pi `$ is 220 Hz; $`a`$ is 100$`a_0`$, where $`a_0`$ is the Bohr radius. The corresponding characteristic length is $`a_{}=1.222\times 10^4`$ cm and the ratio $`a/a_{}`$ is $`4.33\times 10^4`$.

First, we obtain the energy $`E_1`$ (Eq.7) by minimizing it with respect to the variational parameters $`\omega _{}`$, $`\lambda `$ and $`p`$ for various values of the particle number $`N`$. We present the results in Fig.1. It is evident that the results obtained by us are in close agreement with the results in Ref.(see also Table 1 below). As pointed out above, these agree well with the results obtained by using the Gaussian trial wave function when $`N`$ is small and with those obtained by using the Thomas-Fermi approximation when $`N`$ is very large. These comparisons clearly establish the validity of our wave function for a very wide range of the particle numbers.

Next, we compare the proposed wave function with the Gaussian trial wave function of Eq.(4) and also the one given by the Thomas-Fermi approximation in Fig.2. It is clear that when the number of atoms in the trap is small the proposed wavefunction has more resemblence with the Gaussian wavefunction (Fig.2a). As $`N`$ is increased the wave function tends to flatten in the central region and the resulting form is a mixture of the two limiting forms i.e. the Gaussian and the Thomas-Fermi wave function. In the central region it is close to the latter while it resembles the former away from the trap centre (Fig.2b). When $`N`$ is very large the resemblence is more with the Thomas-Fermi wave function (Fig.2c). However we note that the wave function vanishes smoothly far away from the centre of the trap. As mentioned above this a desirable feature which is crucial for the calculation of some relevent physical observables e.g. the aspect ratio . It is clear that our wave function not only provides a better lower bound for the energy but also shows the correct and the desirable limiting behaviour.

Results of calculation of various quantities e.g. chemical potential, total energy , kinetic energy, potential energy, interaction energy, average size of the condensate in the tranverse direction and in the axial direction have been presented in Table 1. The close agreement with the results of Ref. is evident. We have also looked at the variation of the peak density and the total loss rate of the atoms with the particle number $`N`$. We find them to be in good agreement with the result of Ref..

We now present the behaviour of the aspect ratio which is a very important quantity from the experimental point of view . As mentioned in Ref.  it is equal to $`\sqrt{\lambda _0}`$ in the non-interacting limit and tends to $`\lambda _0`$ in the highly repulsive limit, which is the case when $`N`$ is very large. However the convergence to the repulsive limit is very slow  and therefore this behaviour has not been seen explicitly so far. On the other hand, with a variational wave function, calculations can be performed for any $`N`$ with equal ease. Consequently we have been able to verify the predicted behaviour of the aspect ratio. It is seen from Fig.3 that it requires calculations up to $`N10^6`$ to see the aforesaid behaviour.

In Fig.4 we show the wave function of the vortex state corresponding to $`\kappa =1`$ for $`N=5000`$. The atoms are pushed away from the $`z`$ axis. Peak density is 7.155$`\times 10^{13}`$ which occurs at $`r_1=1.611`$. The position of the peak moves further away from the $`z`$ axis as $`N`$ is increased. For $`N=10000`$ the peak occurs at $`r_1=1.844`$ while it is at $`r_1=2.161`$ for $`N=20000`$. Value of the peak density remains much the same for $`N=10000`$ and 20000( 9.524$`\times 10^{13}`$ and 12.26$`\times 10^{13}`$, respectively). We also compare the energy of the vortex state with different values of $`\kappa `$ with that obtained in the Thomas-Fermi approximation . For calculation in the Thomas-Fermi approximation we have neglected the change in the chemical potential due to the presence of a vortex line. Results are presented in Table 2. The Thomas-Fermi approximation yields lower value of the energy which is expected as it ignores the kinetic energy contribution from the core and the surface regions of the vortex state. This approximation is valid only for very large $`N`$. We do find the results of the two calculations to be in close agreement for large $`N`$. For $`N=10^6`$ the two results differ only by $`1\%`$. We have also calculated the critical angular velocity for the vortex state with $`\kappa =1`$ for various values of the particle number $`N`$. We find that it decreases rapidly in the begining (up to $`N2000`$). Thereafter it varies slowly. The result shows good quantitative agreement with those in the Ref.. For example for $`N=2000`$ the critical angular velocity is 52$`\%`$ of the noninteracting value given by the transverse angular frequency $`\omega _{}^0`$ of the trap, in comparison with 49.33$`\%`$ in Ref.. For $`N>5000`$ it is less than 43$`\%`$ of the noninteracting value, which compares well with the figure of 40$`\%`$ in Ref.. The critical angular velocity increases with $`\kappa `$. For $`N=10000`$ we find $`\mathrm{\Omega }_c/2\pi `$ equal to 30.57, 38.5, and 45.42 Hz for $`\kappa =`$1, 2, and 3 respectively. These figures are 26, 35, and 41 Hz respectively in Ref..

### B Negative scattering length: <sup>7</sup>Li

In this section we report calculations on <sup>7</sup>Li. These atoms interact via an attractive interaction and consequently the scattering length in this case is a negative quantity. Numerical values of the parameters used in the calculations correspond to the experimental situation of Ref. and the subsequent theoretical calculations . Accordingly, the asymmetry parameter of the trap is $`\lambda ^0=\omega _z^0/\omega _{}^0=0.72`$. The axial frequency $`\omega _z^0/2\pi `$ is taken to be 117 Hz. The $`s`$-wave scattering length $`a`$ is $`27a_0`$ . The corresponding characteristic length is $`a_{}=2.972\times 10^4`$ cm and the ratio $`a/a_{}`$ is $`4.33\times 10^4`$.

We find the value of the critcal number $`N_c=1270`$ beyond which the ground state collapses because of the attractive interaction. This is in good agreement with the figure of $`N_c1300`$ in Ref. and the experimental observation (see the second paper of Ref.) that $`N1300`$. Wave functions for $`N=`$ 500 and 1270 are shown in Fig.5. For $`N=`$ 500 there is hardly any difference between the proposed wave function and the Gaussian trial wave function. However, the difference is rather significant for $`N=`$ 1270 as is evident from the figure. We plot the aspect ratio for various values of $`NN_c`$ in Fig.6. At $`NN_c`$ the aspect ratio tends to 1. Since for a wave function of the form given by Eq.(5) the aspect ratio also gives ratio of spatial widths in the transverse and the axial directions, the condensate tends to be isotropic for $`NN_c`$. This becomes further evident in Table 3 where we have listed the results for various quantities of interest. As reported in Ref. the variation in the various quantities is smooth from $`N`$ = 1 to 1000. However, we can also note the sharp variation as we reach the critical number. This behaviour is consistent with that reported in Ref.. We have also calculated the peak density and loss rate for $`NN_c`$. Once again we find a sharp increase near $`NN_c`$. These results also match well with those of Ref..

It is possible to have very large number of particles ($`N>>1300`$) in the vortex states even when the interatomic interaction is attractive. We have considered the vortex states with $`\kappa =`$ 1, 2 and 3. The particle number is 3500, 6000 and 8000 respectively. Peak densities for these states are $`1.266\times 10^{13}`$, $`2.239\times 10^{13}`$ and $`2.744\times 10^{13}`$ which occur at $`r_1=`$ 0.922, 1.257 and 1.571 respectively. Although the particle number is quite different in the three cases, the peak densities are not very different. Also they remain less than the peak density $`3.984\times 10^{13}`$ which corressponds to $`\kappa =0`$ and $`N=1270`$. These observations are consistent with those in Ref.. Stability of the vortex states for attractive interaction can be physically explained as the interplay between the restoring force and and the centrifugal force. The restoring force tries to attract the particles to the centre while the centrifugal force tries to push them out. The net effect is that the peak density does not change much even when there is siginificant variation in the particle number. Since the interparticle interaction depends on the density, for low dnesities it does not cause the collapse of the condensate.

In case of attractive interaction it takes more energy to create a vortex state than that requiread in the noninteracting case. Consequenty the critical angular velocity is greater than unity. For $`\kappa =1`$ and $`N=1000`$ we find $`\mathrm{\Omega }_c=1.119`$ which compares very well with $`\mathrm{\Omega }_c=1.2`$ reported in Ref..

## IV Conclusion

We have proposed a variational scheme to describe the ground state and vortex states of weakly interacting atomic gases confined by harmonic traps within the framework of the meanfield theory of Gross and Pitaevskii . It is based on a judicious choice of the form of trial wave function for the ground state which has a simple functional form and at the same time is valid for a wide range of the particle numbers. When the number is small it tends towards a Gaussian and in the opposite limit it resembles the Thomas-Fermi wave function. However, for large $`N`$ it provides a better description of the surface region than the Thomas-Fermi wave function. In the intermediate regime it combines the feature of both in an effective way. In the central region of the trap, where density is high, it matches with the Thomas-Fermi wave function. Away from the centre of the trap, where density is low, it matches the Gaussian trial wave function. We easily generalize the wave function for the vortex states. We have demonstrated the applicability of our method by performing calculations of various physical quantities for the experimental situations of Ref. and Ref.. We find our results to be in good agreement with the existing results. The method is semi-analytic and consequently computationally easy to implement. As our method poses no additional problems even for very large particle numbers we have been able to verify, for the first time, the predicted behaviour of the aspect ratio. The formalism is quite general involving only the scaled $`s`$-wave scattering length and the asymmetry parameter of the trap. We therefore believe that it will be useful in analysing a variety of experiments. In addition it may serve as a very good starting point for the theories where quantum fluctuations play an important role. Generalization of our method to the time dependent case  is straightforward and it will be reported in a future publication.

## V acknowledgment

We thank Dr. S. C. Mehendale and Dr. M. K. Harbola for helpful discussions and critical reading of the manuscript.

Figure Captions

* Ground state energy per atom for <sup>87</sup>Rb as a function of $`N`$. Energy is in the units of $`\mathrm{}\omega _{}^0`$. The solid line is the result of our variational calculation. The dashed line is the result obtained by using the Gaussian trial wave function of Ref., while the dotted line is the result obtained by using the Thomas-Fermi approximation. The filled circles are the results of Ref..
* Ground state wave function for <sup>87</sup>Rb along the $`x`$ axis for different values of $`N`$. The solid line is the result of our variational calculation. The dotted line is the result obtained by using the Gaussian trial wave function of Ref., while the dashed line is the result obtained by using the Thomas-Fermi approximation.(a) $`N=100`$. (b)$`N=1000`$. (c)$`N=50000`$.
* Aspect ratio in <sup>87</sup>Rb as a function of $`N`$. The lower and upper horizontal lines correspond to $`\sqrt{\lambda _0}`$ and $`\lambda _0`$ respectively.
* Vortex-state wave function of 5000 <sup>87</sup>Rb atoms along the $`x`$ axis for $`\kappa =1`$.
* Ground state wave function for <sup>7</sup>Li along the $`x`$ axis for $`N=500`$ (lower curves) and $`N=1270`$ (upper curves). The solid lines are the results of our variational calculation. The dotted line are the results obtained by using the Gaussian trial wave function of Ref..
* Aspect ratio in <sup>7</sup>Li as a function of $`N`$.