Stability Analysis of a Bose-Einstein Condensate Trapped in a Generic Potential

TL;DR

This paper analyzes the stability of a Bose-Einstein condensate in a spherical power law potential using dynamical systems, revealing that stability depends critically on the potential's exponent.

Contribution

It provides a dynamical systems framework for understanding the stability of Bose-Einstein condensates in generic power law traps, highlighting the impact of the potential's exponent.

Findings

Stability for exponents ≤ 2

Instability increases for exponents > 2

Large exponents lead to particle loss from the condensate

Abstract

We investigate the dynamical behavior of the Gross-Pitaevskii equation for a Bose-Einstein condensate trapped in a spherical power law potential restricted to the repulsive case, from the dynamical system formalism point of view. A five-dimensional dynamical system is found (due the symmetry of the Gross-Pitaevskii equation interacting with a potential), where the Thomas-Fermi approximation constrains the parameter space of solutions. We show that for values of the power law exponent equal or smaller than 2 the system seems to be stable. However, when the corresponding exponent is bigger than 2, the instability of the system grows when the power law exponent grows, indicating that large values of the aforementioned parameter can be related to a loss in the number of particles from the condensed state. This fact can be used also to show that the stability conditions of the condensate are highly sensitive to the exponent associated with the external potential.

Figures

Tables

Table 1. Table 1: Critical points for the autonomous system ( 3.3 ) in where exist solutions related to the stability of the Gross–Pitaevskii equation

Label	Critical point
	($x_c,y_c,z_c,w_c,v_c$)
A	($0,-{\overline{\mu }}^2,0,0,0$)
B	($0,-{\overline{\mu }}^2,\pm \sqrt{\frac{2}{\kappa }}\overline{\mu },0,0$)

Table 2. Table 2: The corresponding eigenvalues related to the critical points.

Label	Eigenvalues
	($m_1,m_2,m_3,m_4,m_5$)
A	$(\frac{\sqrt{s^2-s}}{\overline{r}},-\frac{\sqrt{s^2-s}}{\overline{r}},\sqrt{2\overline{\mu }}i,\sqrt{2}\overline{\mu }i,2\sqrt{2}\overline{\mu }i/\overline{r})$
B	$(\frac{\sqrt{s^2-s}}{\overline{r}},-\frac{\sqrt{s^2-s}}{\overline{r}},\sqrt{2\overline{\mu }},\sqrt{2}\overline{\mu },2\sqrt{2}\overline{\mu }/\overline{r})$

Equations

$$i\hslash {\partial }_t\psi (r,t)$$

$$V(r)=A{\left(\frac{r}{a_0}\right)}^s.$$

$$\bar{t}$$

$$\bar{\psi }(\bar{r},\bar{t})$$

$$i{\partial }_{\bar{t}}\bar{\psi }(\bar{r},\bar{t})$$

$$\bar{\psi }(\bar{r},\bar{t})=e^{-i\bar{\mu }\bar{t}}\bar{\varphi }(\bar{r}),$$

$${\bar{\nabla }}^2\bar{\varphi }(\bar{r})=2\left[\bar{V}(\bar{r})-\bar{\mu }\right]\bar{\varphi }(\bar{r})+2\beta {\bar{\varphi }}^3(\bar{r}).$$

$$\frac{d^2\bar{\varphi }(\bar{r})}{d{\bar{r}}^2}+\frac{2}{\bar{r}}\frac{d\bar{\varphi }(\bar{r})}{d\bar{r}}=2\left[\bar{V}(\bar{r})-\bar{\mu }\right]\bar{\varphi }(\bar{r})+2\beta {\bar{\varphi }}^3(\bar{r}).$$

$$n(\bar{r})=|\bar{\varphi }(\bar{r}){|}^2=\frac{\bar{\mu }-\bar{V}(\bar{r})}{\beta },$$

$${\bar{r}}_e=(2\bar{\mu }{)}^{1/s}.$$

$$\begin{gathered} x\equiv {\bar{\varphi }}^{\prime },y\equiv \frac{{\bar{r}}^s}{2}-\bar{\mu },z\equiv \bar{\varphi },w\equiv \frac{s}{2}\left(\frac{{\bar{r}}^s}{\bar{r}}\right), \\ v=\frac{2}{\bar{r}}\bar{\varphi }, \end{gathered}$$

$$F(y,z):=\kappa |z{|}^2=-2y,$$

$$x^{\prime }$$

$$y^{\prime }$$

$$z^{\prime }$$

$$w^{\prime }$$

$$v^{\prime }$$

$$x^{\prime }=2{\alpha }_{1}z+\kappa z^3,z^{\prime }=x,v^{\prime }$$

$$\displaystyle x^{\prime}=2yz+\kappa z^{3}=2\Big{(}\frac{\bar{r}}{2}-\bar{\mu}\Big{)}z+\kappa z^{3},$$

$$y^{\prime }=1/2,z^{\prime }=x,v^{\prime }=\frac{2}{\bar{r}}(x-z),$$

$$\displaystyle x^{\prime}=2yz+\kappa z^{3}=2\Big{(}\frac{\bar{r}^{2}}{2}-\bar{\mu}\Big{)}z+\kappa z^{3},$$

$$y^{\prime }=w,z^{\prime }=x,w^{\prime }=1,v^{\prime }=\frac{2}{\bar{r}}(x-z)$$

$$\displaystyle x^{\prime}=2yz+\kappa z^{3}=2\Big{(}\frac{\bar{r}^{3}}{2}-\bar{\mu}\Big{)}z+\kappa z^{3},$$

$$y^{\prime }=w=\frac{3}{2}{\bar{r}}^2,z^{\prime }=x,w^{\prime }=\frac{6}{{\bar{r}}^2}(y+\bar{\mu })$$

$$v^{\prime }=v+\frac{2}{\bar{r}}{z}^{\prime }.$$

$$\displaystyle x^{\prime}=2yz+\kappa z^{3}=\Big{(}\bar{r}^{4}-2\bar{\mu}\Big{)}z+\kappa z^{3},$$

$$y^{\prime }=w=2{\bar{r}}^3,z^{\prime }=x,$$

$$w^{\prime }=\frac{12}{{\bar{r}}^2}(y+\bar{\mu })=6{\bar{r}}^2,v^{\prime }=v+\frac{2}{\bar{r}}{z}^{\prime }$$

$${ϵ}_p=\frac{p^2}{2m}+V(r)+2U_{0}n(r),$$

$$n(r,p)=\frac{1}{e^{\beta ({ϵ}_p-\mu )}-1},$$

$$N=\frac{1}{(2\pi \hslash {)}^3}\int d^{3}r{d}^{3}pn(r,p),$$

$$n(r)=\int d^{3}pn(r,p),$$

$$n(r)=\lambda^{-3}g_{3/2}\Bigl{(}e^{\beta[\mu-V(r)-2U_{0}n(r)]}\Bigr{)}$$

$$g_{\nu }(z)=\frac{1}{\Gamma (\nu )}{\int }_0^{\infty }\frac{x^{\nu -1}dx}{z^{-1}{e}^x-1}.$$

$$n(r)$$

$$Z=e^{\frac{1}{\kappa T}[\mu -V(r)]},$$

Full text

1 Introduction

The phenomenon of Bose–Einstein condensation Bose ; Bose1 ; Bose2 is one of the most remarkable N–body quantum phenomena that can be observed at macroscopic scales. Since its laboratory observation Bose-Einstein condensation of dilute atomic gases has stimulated an enormous amount of related work Anderson ; A3 ; A4 ; A5 . We found mathematical issues related to Bose–Einstein condensates Lieb , several theoretical and heuristic aspects Dalfovo ; S , also as test tools in gravitational physics, or even to model dark matter in the universe. The study of its associated thermodynamic properties is naturally also a very pertinent aspect Dalfovo ; Pethick ; ueda ; ref ; ref1 ; ref2 . The experimental understanding of Bose–Einstein condensates can be achieved in several alkali atoms which allows the use of various atomic systems, for instance, 87Rb, 23Na, which have repulsive interactions within the condensate Pethick ; TA . Conversely, 7 Li and 85Rb that have attractive interactions Pethick ; TA . The low energy properties and specifically the stability of the condensate properties can be characterized through the s–wave scattering length. These conditions of the condensate trapped in magnetic traps has been extensively analyzed from the numerical and analytical point of view TA . It is seems that the condensate was predicted to be (meta)stable in magnetic traps, when the number of atoms is below some critical number K . In other words, the possible values of the s–wave scattering length are intimately linked to the stability of the system under these circumstances. However, there is at least one more parameter linked to the stability of the condensate, the specific form of the trapping potential. The form of the trapping potential is strongly related to the macroscopic behavior of the condensate, and clearly the properties of the condensed cloud must be highly sensitive to its shape dm . Then, it is interesting to analyze the stability of the condensate when one varies the shape of the trap besides the value of the scattering length. In this work, we will assume that the value of the scattering length is positive, i.e., we assume from the very beginning that our system is a repulsively condensate. However, there are in the literature several works related to the stability analysis when the value of the scattering length is negative 111For works related to the stability of attractively interacting condensates, see refs. ATRA ; ATRA1 ; ATRA2 ; ATRA3 ; ATRA4 . We must also mention that experimental realizations of Bose Einstein condensates in higher power law potential traps has been extensively investigated, see for instance refs. EXP ; EXP1 ; EXP2 ; EXP3 .

Additionally, the dynamical system formalism libro0 is a mathematical tool that is implemented in order to rewrite the evolution equations as plane-autonomous system and to analyze the stability of the equation under consideration. The results of the latter can give some information about the physics of the system behaviour which can be hidded in the corresponding critical points of the system.

Therefore, we will compare the stability conditions of the condensate with the stability of the Gross–Pitaevskii equation for a Bose–Einstein condensate trapped in a spherical power law potential. As far as we know, this mathematical tool has never been used in this context. Thus, the analysis of stability conditions of the Gross–Pitaevskii equation as an autonomous system over the condensate in a generic potential can be helpful in order to find the most favorable scenarios in which the system actually acquires more stable configurations, when repulsive interactions are present.

The outline of the paper is as follows: in Section 2 we describe the main properties of the Gross–Pitaevskii equation and of the corresponding spherical power law potential. Also, some issues regarded the Thomas–Fermi approximation are revised. Section 3 is dedicated to the dynamical analysis of the cases in which a Bose–Einstein condensate is endowed in a power law potential. We first show that there are four cases of interest, related to the values of the corresponding exponent $s$. To have a complete picture of the solution for the Gross–Pitaevskii equation in a Thomas–Fermi approximation, these equations are written as an autonomous system, and we study its critical points and general trajectories in the phase space of the resulting dynamical variables. In Section 4 we analyze some thermodynamical properties of the system and its relation with the results obtained in dynamical analysis. Finally, Section 5 is devoted to final comments about the general properties of the Bose–Einstein condensate trapped in a spherical power law potential.

2 The spherical power law potential and the Gross–Pitaevskii equation

Let us start with the time dependent Gross–Pitaevskii equation in 3–dimensions

[TABLE]

where $m$ is the atomic mass and $U_{0}=4\pi\hbar^{2}a_{s}/m$ describes the interaction between the particles in the condensate with $a_{s}$ the $s$–wave scattering length, which in this approach is assumed to be positive. The validity of the Gross–Pitaevskii equation is based on the condition that the $s$-wave scattering length be much smaller than the average distance between atoms and that the number of atoms in the condensate be much larger than one. The Gross–Pitaevskii equation can be used, at temperatures below the condensation temperature, to explore the macroscopic behavior of the system, characterized by variations of the order parameter $\psi(r,t)$ over distances larger than the mean distance between atoms.

The trapping potential in Eq. (2.1) is given by:

[TABLE]

The potential (2.2) is the so–called spherical power law potential where, for simplicity, we assume spherical symmetry (see for instance Refs. zobay ; zobay1 ; zobay2 and references therein for more general potentials). In the generic potential (2.2) we have $A=\hbar\omega_{0}/2$ and $a_{0}=\sqrt{\hbar/m\omega_{0}}$, which are characteristic scales of energy and length associated with the trap eli ; eli1 , with $\omega_{0}$ the corresponding frequency and $\hbar$ the reduced Planck’s constant. The parameter $s$ depicts the exponent of the spherical power law potential, which in our approach will play a relevant role in the dynamics. Notice that if we set $s=2$ in the spherical power law potential (2.2), we recover the more common harmonic trap. These type of potentials are known as generic spherical potentials. It is quite interesting to analyze the behavior of the condensate when the shape of the trap varies as a power of the radial coordinate for several reasons. This class of trapping potentials could be useful in order to cool the system with an adiabatic procedure, by changing the shape of the trap Dalfovo . The relevant thermodynamical properties associated with the system, e.g., the condensation temperature, the number of particles within the condensate, the corresponding density, the entropy, etc., explicitly exhibits a sensitive trap–dependence, i.e., the macroscopic behavior of the condensate is highly sensitive to the shape of the external trapping potential. In this sense, it is relevant to analyze the stability of the system when this spherical power law potential is present.

In order to obtain a dimensionless version of the Gross–Pitaevskii equation (2.1), let us define, as usual, the following dimensionless variables in function of the characteristic scales associated with the system M.M ; M.M1 ; M.M2

[TABLE]

By using the dimensionless variables defined above, the Gross–Pitaevskii equation in 3–dimensions becomes

[TABLE]

where we have defined $\beta=U_{0}/\hbar\omega_{0}a^{3}_{0}$. As was mentioned above, we restrict our analysis to positive values of the corresponding self–interaction potential since positive scattering lengths provides repulsive interactions among the particles.

Consequently, the time independent dimensionless Gross–Pitaevskii equation can be obtained by written the dimensionless order parameter $\bar{\psi}(\bar{r},\bar{t})$ as follows:

[TABLE]

where we have defined the dimensionless chemical potential as $\bar{\mu}=\mu/\hbar\omega_{0}$. Inserting (2.5) into (2.4), this leads to the following equation

[TABLE]

where the Laplacian operator is given by $\bar{\nabla}^{2}=\bar{r}^{-2}\partial_{\bar{r}}\bar{r}^{2}\partial_{\bar{r}}$. Consequently Eq. (2.6) reads

[TABLE]

An analytical and formal solution associated with our system can be obtained by neglecting the kinetic energy term from the very beginning in Eq. (2.6), this is the so–called Thomas–Fermi approximation Pethick , with the result

[TABLE]

while $n(r)=0$ outside of this region 222 We are able to estimate the spatial extend $\bar{r}_{e}$ of the cloud in the dimensionless variables, through the condition $\bar{\mu}=\bar{V}(\bar{r})$, which in our case is given by

$\bar{r}_{e}=(2\bar{\mu})^{1/s}.$

(2.9)

This spatial extend fixes the regime of validity of the Thomas–Fermi approximation and consequently, also fixes the validity of the analysis followed in this paper.. Here, $n(\bar{r})=|\bar{\phi}(\bar{r})|^{2}$ depicts the density of particles within the condensate in the Thomas–Fermi approximation. Therefore, the Thomas–Fermi approach is valid for systems at very low temperatures $T<T_{c}$, where $T_{c}$ is the condensation temperature, for weakly interacting systems, for sufficient large clouds, and when the kinetic energy is negligible with respect to the potential and the interaction ones. Additionally, the Thomas–Fermi approximation is valid near to the minima of the potential. Finally, we must add that the Thomas-Fermi approximation fails for trapped condensates near the edge of the cloud, due to the divergent behavior of the kinetic energy (i.e. the total kinetic energy per unit area diverges on the boundary of the system) DAL1 .

3 Dynamical system structure

The first step to approach to our goal consists in study the evolution of the Gross–Pitaevskii equation (2.6) in which the condensate is endowed with the spherical power law potential (2.2). As in the standard case, it is possible to perform a dynamical study of the model so that its relevant solutions are easily unveiled Copeland:1997et ; EscamillaRivera:2010py ; Boehmer:2014vea . We will see that the Thomas–Fermi approximation leads to a constraint equation on the parameter space whose solution defines a subspace where stable solutions of the system are possible. In order to construct a dynamical system for a Bose–Einstein condensate trapped in the spherical power law potential (2.2), we first introduce a set of conveniently chosen variables which may allow rewriting the evolution equation as an autonomous system subject to the constraint arising from the Thomas–Fermi approximation. We choose the following variables

[TABLE]

where the prime denote derivatives with respect to $\bar{r}$. Let us clarify that these dynamical coordinates correspond to the phase space. In these variables the Thomas–Fermi approximation constraint is implemented by setting

[TABLE]

where $\kappa\equiv 2\beta$. We shall restrict ourselves to the parameter space delimited by this hyperbolic region, which contains the solutions in this approximation for the Gross–Pitaevskii equation. Combining (2.7) and the variables (3.1), the equations of motion read

[TABLE]

The stability of the Gross–Pitaevskii equation through its critical points is investigated by using linear perturbations (See Appendix A). The results of the analysis of the dynamical system (3.3), their critical points and its stability properties, are summarized in Table 1. At this point we are interested in where the trapping potential, (i.e., for different values of the exponent $s$), can set a stable (or unstable) scenario. To proceed with these ideas, we follow the cases related to the Bose–Einstein condensate physics. Notice that if we set $s=0$ in the spherical power law potential (2.2), we obtain a free gas. For $s=1$ a linear trap and for $s=2$ the usual harmonic oscillator. Additionally, when $s\rightarrow\infty$, we obtain the box potential.

Let us now establish the stability of the Gross–Pitaevskii equation conditions for some potentials of interest.

• •

Case $s=0$, i.e $\bar{V}=\text{const}$. The autonomous system is reduced to the following three dimensional system:

[TABLE]

where $\alpha_{1}=1/2-\bar{\mu}$. The real critical point for this system is $(x_{c}=0,z_{c}=0,v_{c}=0)$. The constraint is fixed by $|z|^{2}=-2\alpha_{1}/\kappa$.

The trajectories in Figure 1 for the case $\xi_{+}$ of $m_{+}$ start at the critical point then diverge to infinite-distance out. As a result, the critical point for the case results in a saddle point, which is always unstable. When $s=0$ the potential becomes in a constant that can be re-absorbed in the dimensionless chemical potential. This type of potential may not be of physical relevance. However, we have included this case for mathematical reasons. In this scenario, the system is unstable, according to this approach, which makes sense since there is no confinement region. This behavior is showed in Figure 1.

• •

Case $s=1$, i.e $\bar{V}=\bar{r}/2$, which corresponds to the linear trapping potential. The autonomous system is reduced to the following set of equations:

[TABLE]

with two critical points at $(x_{c_{1}}=0,y_{c_{1}}=y,z_{c_{1}}=0,v_{c_{1}}=0)$ and $(x_{c_{2}}=0,y_{c_{2}}=-\kappa z^{2}/2,z_{c_{2}}=z,v_{c_{2}}=2z_{c_{2}}/\bar{r})$ The constraint is fixed by the curve $|z|^{2}=(2\bar{\mu}-\bar{r})/\kappa$. The instability happens only along the eigenvectors corresponding to $\pm\sqrt{2y}$ and $\pm\sqrt{2\kappa}z$, which are unstable (saddle) points. The first eigenvalues are null, and then the system is indifferent under perturbations along the constraint curve for this case. Also we notice that for a vanishing potential $\bar{\mu}\rightarrow 0$, the fourth dynamical variable is $v=z/y$, therefore fixing a new minimal constraint in one dimension according to (3.2) as $F(z):=4\beta z^{3}=-v_{(s=1,\bar{\mu}\rightarrow 0)}$. The phase space of solutions for this case is reduced to a one-dimension in comparison with Eq.(3.2). Inside this constriction, only the unstable solutions with eigenvalues $\pm\sqrt{2\kappa}z$ are allowed. Additionally, according to the spatial extend for this case $\bar{r}_{e}=2\bar{\mu}$, the above constraint indicates that when $\bar{\mu}\rightarrow 0$ then $\bar{r}_{e}$ decreases with respect to the $\bar{r}_{e}$ obtained from Eq. (3.2).

• •

Case $s=2$, i.e $\bar{V}=\bar{r}^{2}/2$, which corresponds to the more common harmonic oscillator potential. The autonomous system is reduced in this case to the system:

[TABLE]

with two critical points at $(x_{c_{1}}=0,y_{c_{1}}=y,z_{c_{1}}=0,w_{c_{1}}=0,v_{c_{1}}=0)$ and $(x_{c_{2}}=0,y_{c_{2}}=-(\kappa z^{2})/2,z_{c_{2}}=z,w_{c_{2}}=0,v_{c_{2}}=2z_{c}/\bar{r})$. The constraint is fixed by the curve $|z|^{2}=(2\bar{\mu}-\bar{r}^{2})/\kappa$. For the origin critical point we observe that only the non-vanishes eigenvalues, which points in $z$-direction and $v$-direction gives the information about the stability of this point. In such case, when $\bar{r}>2\bar{\mu}$ ($\bar{r}<2\bar{\mu}$) the point will be unstable (stable). For the case $s=2$, or the usual harmonic trap, the system seems to be stable in all cases. In other words, note that in our approach the Thomas–Fermi approximation fixes the permisible values for the spatial extend of the cloud as $\bar{r}_{e}=(2\bar{\mu})^{1/s}$, this means that for all $\bar{r}$ of interest $\bar{r}<(2\bar{\mu})^{1/s}$, and consequently $\bar{r}<2\bar{\mu}$ is always satisfied, leading to a stable scenario.

• •

Case $s=3$, i.e $\bar{V}=\bar{r}^{3}/2$, which corresponds to a cubic potential. The autonomous system is reduced to the system:

[TABLE]

with three critical points at $(x_{c_{1}}=0,y_{c_{1}}=-\bar{\mu},z_{c_{1}}=0,w_{c_{1}}=0,v_{c_{1}}=0)$ and $(x_{c_{(2,3)}}=0,y_{c_{(2,3)}}=-\bar{\mu},z_{c_{(2,3)}}=\pm\sqrt{\frac{2\bar{\mu}}{\kappa}},w_{c_{(2,3)}}=0,v_{(2,3)}=0)$. The constraint is fixed by the curve $|z|^{2}=2\bar{\mu}/\kappa-8z^{3}/v^{3}$. To have a Bose–Einstein condensate the values of $\bar{\mu}>0$, i.e the first energy state is a non-vanishing constant. Then, for the first critical point the trajectories stay in an elliptical traces. However, with each revolution, their distance from $c_{1}$ grow/decay exponentially according to the term $e^{m_{+}\bar{r}}$ (see (A.4)). Therefore, the 5-dimensional phase portrait shows trajectories that attractor to the critical point (when $m_{c_{1}}>0$). Or trajectories that loiter to the attractor converge to the critical point (when $m_{c_{1}}<0$). As a result, when $\bar{\mu}>0$, $c_{1}$ is a attractor point, which it is asymptotically stable if $m_{c_{1}}<0$ and it is unstable if $m_{c_{1}}>0$.

• •

Case $s=4$, i.e $\bar{V}=\bar{r}^{4}/2$, which corresponds to the quartic potential. The autonomous system is now reduced to the following system of equations:

[TABLE]

with four critical points at $(x_{c_{(1,2)}}=0,y_{c_{(1,2)}}=\pm\bar{\mu},\\ z_{c_{(1,2)}}=0,w_{c_{(}1,2)}=0)$ and $(x_{c_{(3,4)}}=0,y_{c_{(3,4)}}=-\bar{\mu},z_{c_{(3,4)}}\\ =\pm\sqrt{\frac{2\bar{\mu}}{\kappa}},w_{(3,4)}=0,v_{(3,4)}=0)$. The constraint is fixed by the curve $|z|^{2}=2\bar{\mu}/\kappa-8z^{4}/v^{4}$. As in the case $s=3$, for the first critical point the trajectories are elliptical. Therefore, we have the same 5-dimensional phase portrait with spiral trajectories when $\bar{\mu}>0$, $c_{1}$ is a spiral point, which it is asymptotically stable if $m_{c_{1}}<0$ and it is unstable if $m_{c_{1}}>0$. Notice that $m_{c_{1}}$ is always positive according to Table (2) and then is unstable. In comparison the latter case, we have a critical point for $\bar{\mu}<0$, which means $\bar{r}_{e}=-2\bar{\mu}^{1/4}$, which is a mathematical solution but is unphysical since $\bar{\mu}$ is always positive for temperatures smaller that the condensation temperature $T_{c}$.

Furthermore, for the cubic potential $s=3$ and the quartic potential $s=4$ the numerical stability analysis shows that these systems present an instability due to the presence of attractors trajectories. If we interpreted these trajectories as the velocity field of the particles in the condensate, we can immediately observe that when $s$ grows the instability increases. These facts can be interpreted as follows, when $s>2$ grows, the particles forming the condensate tends to leave the ground state to form a sea of excited particles, which interacts with the condensate affecting its stability.

Finally, let us remark that due to the corresponding limitations associated with the Thomas–Fermi approach used in the above analysis, the deduced instabilities are valid in a certain region near to the minima of the potential, i.e., for the region $\bar{r}<(2\bar{\mu})^{1/s}$. Thus, as a consequence we can not extrapolate our results to the entire condensate.

4 Condensation temperature and condensate fraction

The second step to approach to our goal consists of an analysis of some of the thermodynamic properties of our bosonic gas trapped in the spherical power law potential Eq. (2.2) in order to study its relation to the dynamical structure analyzed in the previous section. The properties of a condensate trapped in a spherical power law potential have been extensively analyzed, see for instance, zobay ; eli ; H ; L ; sal and references therein.

To this aim, let us propose a particularly simple Hartree–Fock type spectrum, in the semi–classical approximation. This approach basically consists of the assumption that the constituents of the gas behave like non–interacting bosons moving in a self–consistent mean field, valid when the semiclassical energy spectrum $\epsilon_{p}$ is bigger than the corresponding chemical potential $\mu$, for dilute gases Dalfovo ; Pethick

[TABLE]

where $p$ is the momentum and $m$ is the mass of a single particle respectively. The term $2U_{0}n(\vec{r})$ is a mean field generated by the interactions with the other constituents of the bosonic gas, being $n(r)$ the spatial density of the cloud Dalfovo .

In the semiclassical approximation, the single–particle phase–space distribution may be written as Dalfovo ; Pethick

[TABLE]

where $\beta=1/\kappa T$, $\kappa$ is the Boltzmann constant, $T$ is the temperature, and $\mu$ is the chemical potential. The number of particles in the 3–dimensional space obeys the normalization condition Dalfovo ; Pethick ,

[TABLE]

where

[TABLE]

is the corresponding spatial density. Using expression (4.1), and integrating expression (4.2) over momentum, with the help of (4.4), we get the spatial distribution associated with our semi–classical spectrum (4.1)

[TABLE]

where \lambda=\Bigl{(}\frac{2\pi\hbar^{2}}{m\kappa T}\Bigr{)}^{1/2}, is the de Broglie thermal wavelength. The function $g_{\nu}(z)$ is the so–called Bose–Einstein function defined by Phatria

[TABLE]

By using the properties of the Bose–Einstein functions Phatria , we are able to expand expression (4.5) to first order in $U_{0}$, with the result

[TABLE]

where

[TABLE]

being $n_{0}(r)$ the space density distribution for the ideal case $U_{0}=0$,

[TABLE]

Integrating the normalization condition (4.3) and using expression (4.7) with the corresponding potential (2.2), this allows us to obtain an expression for the number of particles as a function of the chemical potential $\mu$, the temperature $T$, the coupling constant $U_{0}$

[TABLE]

with the following definition

[TABLE]

where $\Gamma(x)$ is the Euler function, $N_{0}$ the number of particles in the grown state and $z=\exp^{\frac{\mu}{\kappa T}}$ is the so–called fugacity.

Setting $U_{0}=0$ in Eq. (4.10), we are able to calculate the condensation temperature $T_{0}$ in the ideal case. To do this, we assume that the number of particles in the ground state in negligible, together with $\mu=0$ at the condensation temperature. Thus, we obtain

[TABLE]

where we have used that $g_{\nu}(z=1)=\zeta(\nu)$, being $\zeta(\nu)$ the Riemann Zeta Function.

Since the total number of particles is given by $N=N_{0}+N_{e}$, where $N_{0}$ is the number of particles in the ground state and $N_{e}$ is the number of particles in the exited states, the condensate fraction can be written as $N_{0}/N=1-N_{e}/N$, and according to (4.10) and (4.12), at temperatures $T$ below $T_{0}$ , the condensate fraction is given by:

[TABLE]

where as usual, we have used the fact that the number of particles that can be located in the excited states has a maximum

[TABLE]

In order to obtain the shift in the condensation temperature caused by interactions, let us expand expression (4.10) to first order in $T=T_{0}$, $\mu=0$ and, $U_{0}=0$. Recalling that $T_{0}$ is the condensation temperature given by expression (4.12) (see for instance Refs. eli ; sal ), with the result :

[TABLE]

where

[TABLE]

Where we have setting in Eq. (4.15) the condensation conditions in the case of weakly interacting systems, i.e., at the condensation temperature $T_{c}$ for large $N$, in the mean field approach the chemical potential takes the value Dalfovo ; Pethick

[TABLE]

In other words, $n_{0}(r=0)=\lambda_{c}^{-3}\zeta(3/2)$ in the large $N$ limit, which means that the critical density near the center of the spherical power law potential for any $s$ is the same as that of the uniform model. For instance, in the case $s=2$, i. e., the usual harmonic oscillator potencial the current high precision experiments for ${}^{39}_{19}K$, show that the shift in the condensation temperature with respect to the ideal result, caused by the interactions among the constituents of the gas is about $5\times 10^{-2}$ with a 1% of error RP . In other words, when repulsive interactions ($a>0$) are present, the condensation temperature (depending of the exponent $s$) is smaller compared with the ideal case, see Eq. (4.15).

Additionally, when interactions are present the condensate fraction can be obtained from eqns. (4.10), as

$N_{0}/N=1-N_{e}/N$, where, as in the non–interacting case the number of particles that can be located in the excited states can be expressed as

[TABLE]

We show in Fig. 2 the plot for some specific values of the exponent $s$. Moreover, we have used typical laboratory conditions, i. e., $N\sim 10^{5}$ particles, $a\sim 10^{-9}$ m, and $m\sim 10^{-26}$ kg Dalfovo , with a typical condensation temperature of order $T_{c}\sim 10^{-9}$ $K$.

Note that when interactions are present and for temperatures smaller than $T_{c}$ the condensed fraction is also smaller compared with the ideal case for $a>0$. From this analysis we infer that when $s$ grows, the number of particles in the condensate decreases, which is in agreement with the results given by the dynamical analysis, in which we interpreted the increasing of $s$ as a increase of instability. For instance, the analysis performed in Sec. (3) together with the results exposed in the present section, allow us to conclude that the cases $s=3$ and $s=4$ are less stable than the usual harmonic trap $s=2$. In this line of thinking, the case $s\rightarrow\infty$ which describes a box potential, is clearly, less stable than the aforementioned scenarios. Also, the spatial extend of the condensate deduced from the Thomas–Fermi approach decreases as the exponent $s$ grows.

5 Discussion

In order to study the stability of the condensate when the exponent of the trap varies, we have analyzed the stability conditions of a Bose–Einstein condensate in a spherical power law potential using techniques from dynamical systems. Also, we have confronted these results with the analysis of the corresponding condensed fraction, showing that both predictions are compatible.

We observe that when the exponent $s$ grows, the system tends to be unstable, i.e., when the exponent fulfill the condition $s>2$. Basically, the perturbations deduced in the above sections are directly related to perturbations of the order parameter (see the system (A.1)), and consequently, related to fluctuations of the condensed state. This fact can be interpreted as follows: when $s$ grows, the number of particles in the ground state decreases and clearly the number of particles in the excited states increase. Also, the size of the cloud decreases as the exponent $s$ increases. Let us add, that according to our results the most stable scenarios for the condensate when a positive self–interaction parameter is assumed, are given for small enough values of the exponent $s$.

Let us mention here that according to our interpretation, i.e., that the instability means a loss of particles of the condensate when the exponent $s$ grows, does not means that the equilibrium condensate density in the instable cases, i.e. $s>2$ is zero, or at least does not allows to conclude this fact. As was analyzed in the previous sections, the loss of particles of the condensed state grows when $s$ grows. In other words, the so–called depletion increases with the exponent $s$. A more suitable interpretation is that the systems with bigger exponent $s$ are less stable than the systems with small exponent $s$. However, this fact do not implies that there are not a well defined equilibrium condensate densities for $s>2$. Moreover, as was mentioned previously in the manuscript, the limitations associated with the Thomas–Fermi approach used in the above analysis, allows us to deduce instabilities in a certain region, i.e., $\bar{r}<(2\bar{\mu})^{1/s}$. Thus, as a consequence we can not extrapolate completely our results to the entire condensate but can be used as an insight of the condensate stability.

In addition, the dynamical system formalism can extended to more general traps, even time–dependent traps and rotating traps in order to analyze the stability of the system under more general conditions, where it would be interesting to show if the solutions are independent of the initial conditions. Also, a numerical simulation can be used in order to match the results obtained in the present report with those eventually obtained in a numerical analysis. Clearly these topics deserve deeper investigation, which will be presented elsewhere.

6 Acknowledgments

E. C. and C. E-R acknowledges MCTP/UNACH for financial support. We thanks to P. Sloane for his opinion on the manuscript.

7 Authors contributions

All the authors were involved in the preparation of the manuscript and all authors contributed equally to the paper. All the authors have read and approved the final manuscript.

Appendix A Dynamical system background

The stability of a particular system of equations system can be studied through its critical points and performing linear perturbations around the critical values of the form $\mathbb{x}=\mathbb{x_{0}}+\mathbb{u}$, where $\mathbb{x}=(x,y,z)$ and $\mathbb{u}=(\delta x,\delta y,\delta z)$. The equations of motion (3.3) can be written as $\mathbb{x}^{\prime}=\mathbb{f}(\mathbb{x})$, which upon linearization reads

[TABLE]

where $\mathcal{M}$ is called the linearization matrix. The eigenvalues $m$ of $\mathcal{M}$ determine the stability of the critical points, whereas the eigenvectors determine the principal directions of the perturbations. In general, if $Re(m)<0$

($Re(m)>0$) the critical point is called stable (unstable). We should study the perturbations of the three dynamical variables $(x,y,z)$ in the phase space, but remember that they are not all independent because they are bound by the Thomas–Fermi approximation (2.8), and the same happens to their perturbations. This constraint defines a two dimensional surface (3.2) upon which lie all physically relevant phase space trajectories (solutions), and then we will be interested on perturbations lying also on the constraint surface. In other words, perturbations which are perpendicular to the Thomas–Fermi constraint surface should be excluded from the analysis.

To resolve Eq.(A.1) we list the eigenvalues of the stability matrix $\mathcal{M}$ for each of the critical point with the following ansatz at first order in the perturbations:

[TABLE]

where

[TABLE]

is the system to resolve which has solutions

[TABLE]

where $[m_{+},m{-}]$ are the eigenvalues of the system and $\xi$ the eigenvectors and $[u_{+},u_{-},v{+},v{-},w_{+},w_{-}]$ are arbitrary constants which are determined by the initial conditions.