Analytical solutions for the spin-1 Bose–Einstein condensate in a harmonic trap

Abstract

The homotopy analysis method and Galerkin spectral method are applied to find the analytical solutions for the Gross–Pitaevskii equations, a set of nonlinear Schrödinger equation used in simulation of spin-1 Bose–Einstein condensates trapped in a harmonic potential. We investigate the one-dimensional case and get the approximate analytical solutions successfully. Comparisons between the analytical solutions and the numerical solutions have been made. The results indicate that they are in agreement well with each other when the atomic interaction is weakly. We also find a class of exact solutions for the stationary states of the spin-1 system with harmonic potential for a special case.

Graphical abstract

Keywords

spin-1 Bose–Einstein condensate / Gross–Pitaevskii equation / homotopy analysis method / analytical solution

1 Introduction

In recent years, the experimental achievement of spinor condensates [1–5] offers new regimes to study various quantum phenomena that are generally absent in a single-component condensate. For the F=1 spinor condensate, such as ²³Na, ³⁹K, ⁸⁷Rb, there are three degrees of freedom represented by magnetic quantum number, m_{\mathrm{F}}=-1, 0, 1 due to the vectorial nature. As a result, the order parameter, or wave function, of the condensate characterizing the F=1 spinor condensate have three components, denoted by \Psi ={({\psi }_1,{\psi }_0,{\psi }_{-1})}^{\mathrm{T}}, where \mathrm{T} denotes the transpose. Theoretically, the spin-1 BEC was first studied by Ho [6], Ohmi and Machida [7] by generalizing the Gross–Pitaevskii equation (GPE) under the restriction of gauge and spin-rotation symmetry. The multi-component GPEs describing the evolution of the wave functions of the condensate can be obtained from the derivative of the energy functional, \mathrm{i}\mathrm{\hslash }\frac{\partial \Psi }{\partial t}=\frac{\delta E[\Psi ]}{\delta {\Psi }^{\ast }}. When the Zeeman terms are negligible, the energy functional E[\Psi ] is given by [6, 8]

\begin{array}{}\text{(1)}& \begin{array}{cc}E[\Psi ]=\hfill & {\displaystyle \int }\mathrm{d}\overrightarrow{r}\{{\displaystyle \sum _{\alpha =-1}^1}{\psi }_{\alpha }^{\ast }[-{\displaystyle \frac{{\mathrm{\hslash }}^2}{2m}}{\nabla }^2+U(\overrightarrow{r})]{\psi }_{\alpha }\hfill \\ \hfill & +{\displaystyle \frac{c_0}{2}}{\rho }^2+{\displaystyle \frac{c_2}{2}}{⟨\widehat{f}⟩}^2\}\hfill \end{array}\end{array}

where m is the atomic mass, \overrightarrow{r}=(x,y,z) the Cartesian coordinate vector, U(\overrightarrow{r}) the external potential, \mathrm{\hslash } the reduced Planck constant, \rho ={\sum }_{\alpha =-1}^1{|{\psi }_{\alpha }|}^2 the particle density. The coupling constants c_0,c_2 govern the nonlinear interaction between different spin components of the condensate and c_0=\frac{g_0+2g_2}{3}=\frac{4\pi {\mathrm{\hslash }}^2}{3m}(a_0+2a_2),c_2=\frac{g_2-g_0}{3}=\frac{4\pi {\mathrm{\hslash }}^2}{3m}(a_2-a_0). The suffixes 0, 2 denote the total spin (for F=1, total spin 0 or 2 are possible for bosons) and g_k=\frac{4\pi {\mathrm{\hslash }}^2{a}_k}{m}(k=0,2) characterize the strength of the interaction between two particles having a total spin of F and it is related to the corresponding s-wave scattering length a_s, which can be adjusted experimentally by optical Feshbach resonance [9]. The usual mean-field interaction c₀ is positive for repulsive interaction and negative for attractive interaction. However, the energy functional shows that c₀ must be nonnegative; otherwise, the system would collapse. The spin-exchange interaction c₂ is positive for anti-ferromagnetic (polar) interaction [10] and negative for ferromagnetic interaction [11]. ⟨\widehat{\mathcal{f}}⟩=(⟨{\widehat{f}}_x⟩,⟨{\widehat{f}}_y⟩,⟨{\widehat{f}}_z⟩) is the spin vector and ⟨{\widehat{f}}_{\alpha }⟩={\sum }_{k,l=-1}^1{\psi }_k^{\ast }{(f_{\alpha })}_{kl}{\psi }_l(\alpha =x,y,z) are the spin density, where f_x,f_y,f_z are the spin-1 matrices

\begin{array}{}\text{(2)}& \begin{array}{c}{f}_x={\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{ccc}0& 1& 0\\ 1& 0& 1\\ 0& 1& 0\end{array}\right),f_y={\displaystyle \frac{\mathrm{i}}{\sqrt{2}}}\left(\begin{array}{ccc}0& -1& 0\\ 1& 0& -1\\ 0& 1& 0\end{array}\right)\hfill \\ f_z=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& -1\end{array}\right)\hfill \end{array}\end{array}

From Eq. (1) the coupled field equations are obtained in the form

\begin{array}{}\text{(3)}& \begin{array}{cc}\mathrm{i}\mathrm{\hslash }{\displaystyle \frac{\partial {\psi }_1}{\partial t}}=\hfill & [-{\displaystyle \frac{{\mathrm{\hslash }}^2}{2m}}{\nabla }^2+U(\overrightarrow{r})+c_0\rho +c_2(|{\psi }_1{|}^2+|{\psi }_0{|}^2\hfill \\ \hfill & -|{\psi }_{-1}{|}^2)]\psi {}_1+c{}_2{\psi }_0^2{\psi }_{-1}^{\ast }\hfill \\ \mathrm{i}\mathrm{\hslash }{\displaystyle \frac{\partial {\psi }_0}{\partial t}}=\hfill & [-{\displaystyle \frac{{\mathrm{\hslash }}^2}{2m}}{\nabla }^2+U(\overrightarrow{r})+c_0\rho +c_2(|{\psi }_1{|}^2\hfill \\ \hfill & +|{\psi }_{-1}{|}^2)]\psi {}_0+2c{}_2\psi {}_1{\psi }_0^{\ast }\psi {}_{-1}\hfill \\ \mathrm{i}\mathrm{\hslash }{\displaystyle \frac{\partial {\psi }_{-1}}{\partial t}}=\hfill & [-{\displaystyle \frac{{\mathrm{\hslash }}^2}{2m}}{\nabla }^2+U(\overrightarrow{r})+c_0\rho +c_2(|{\psi }_0{|}^2\hfill \\ \hfill & +|{\psi }_{-1}{|}^2-|{\psi }_1{|}^2)]\psi {}_{-1}+c{}_2{\psi }_0^2{\psi }_1^{\ast }\hfill \end{array}\end{array}

with the normalization condition {\sum }_{\alpha =-1}^1\int {|{\psi }_{\alpha }|}^{2}d\overrightarrow{r}=N, where N is the total number of atoms in the condensate, and the asterisk denotes the complex conjugation.

In most cases, the confinement due to the magnetic trap is well-described by a harmonic oscillator potential

\begin{array}{}\text{(4)}& U(\overrightarrow{r})={\displaystyle \frac{1}{2}}m({\omega }_x^2{x}^2+{\omega }_y^2{y}^2+{\omega }_z^2{z}^2)\end{array}

with {\omega }_x,{\omega }_y and {\omega }_zbeing the trap frequencies in the x-, y-, and z-directions, respectively. In the disk-shaped condensates, i.e., {\omega }_x\approx {\omega }_y and {\omega }_z\gg {\omega }_x, the three-dimensional (3D) GPEs can be reduced to 2D GPEs. In the cigar-shaped condensates, i.e., {\omega }_y\gg {\omega }_x and {\omega }_z\gg {\omega }_x, the 3D GPEs can be reduced to 1D GPEs.

In this paper, we only consider the cigar-shaped condensates. By introducing the dimensionless variables t^{\prime }={\omega }_{x}t, x^{\prime }=x/a_{\mathrm{ho}}, {\psi }_{\alpha }^{\prime }=\sqrt{a_{\mathrm{ho}}^3/N}{\psi }_{\alpha }(\alpha =0,\pm 1), where a_{\mathrm{ho}}=\sqrt{\frac{\mathrm{\hslash }}{m{\omega }_x}} is the characterized length of harmonic oscillator, and omitting the primes, we get the dimensionless coupled GPEs from Eq. (3)

\begin{array}{}\text{(5)}& \begin{array}{cc}\mathrm{i}{\displaystyle \frac{\partial {\psi }_1}{\partial t}}=\hfill & {\widehat{H}}_0{\psi }_1+[{\beta }_{0}n+{\beta }_2(|{\psi }_1{|}^2+|{\psi }_0{|}^2\hfill \\ \hfill & -|{\psi }_{-1}{|}^2)]\psi {}_1+\beta {}_2{\psi }_0^2\psi {}_{-1}{}^{\ast }\hfill \\ \mathrm{i}{\displaystyle \frac{\partial {\psi }_0}{\partial t}}=\hfill & {\widehat{H}}_0{\psi }_0+[{\beta }_{0}n+{\beta }_2({|{\psi }_1|}^2+{|{\psi }_{-1}|}^2)]{\psi }_0\hfill \\ \hfill & +2{\beta }_2{\psi }_1{\psi }_0^{\ast }{\psi }_{-1}\hfill \\ \mathrm{i}{\displaystyle \frac{\partial {\psi }_{-1}}{\partial t}}=\hfill & {\widehat{H}}_0{\psi }_{-1}+[{\beta }_{0}n+{\beta }_2(|{\psi }_0{|}^2+|{\psi }_{-1}{|}^2\hfill \\ \hfill & -|{\psi }_1{|}^2)]\psi {}_{-1}+\beta {}_2{\psi }_0^2\psi {}_1{}^{\ast }\hfill \end{array}\end{array}

with the normalization condition

\begin{array}{}\text{(6)}& {\displaystyle \sum _{\alpha =-1}^1}{\displaystyle {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}}{|{\psi }_{\alpha }(x,t)|}^{2}dx=1\end{array}

where n={\sum }_{\alpha =-1}^1{|{\psi }_{\alpha }|}^2, {\gamma }_y={\omega }_y/{\omega }_x,{\gamma }_z={\omega }_z/{\omega }_x, {\beta }_0=\frac{c_{0}mN}{{\mathrm{\hslash }}^2{a}_{\mathrm{ho}}}\frac{\sqrt{{\gamma }_y{\gamma }_z}}{2\pi }=\frac{4\pi N(a_0+2a_2)}{3a_{\mathrm{ho}}}\frac{\sqrt{{\gamma }_y{\gamma }_z}}{2\pi }, {\beta }_2=\frac{c_{2}mN}{{\mathrm{\hslash }}^2{a}_{\mathrm{ho}}}\frac{\sqrt{{\gamma }_y{\gamma }_z}}{2\pi }=\frac{4\pi N(a_2-a_0)}{3a_{\mathrm{ho}}}\frac{\sqrt{{\gamma }_y{\gamma }_z}}{2\pi }, and

\begin{array}{}\text{(7)}& {\widehat{H}}_0=-{\displaystyle \frac{1}{2}}{\displaystyle \frac{{\partial }^2}{\partial x^2}}+{\displaystyle \frac{1}{2}}{x}^2\end{array}

We also define the dimensionless magnetization per particle

\begin{array}{}\text{(8)}& M={\displaystyle {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}}\left({|{\psi }_1|}^2-{|{\psi }_{-1}|}^2\right)dx,-1⩽M⩽1\end{array}

The dimensionless energy (measured by \mathrm{\hslash }{\omega }_x) per particle can be expressed as

\begin{array}{}\text{(9)}& \begin{array}{cc}E=\hfill & {\displaystyle {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}}\{{\displaystyle \sum _{\alpha =-1}^1}({\displaystyle \frac{1}{2}}|{\displaystyle \frac{\partial {\psi }_{\alpha }}{\partial x}}{|}^2+{\displaystyle \frac{1}{2}}{x}^2|{\psi }_{\alpha }{|}^2)+{\displaystyle \frac{{\beta }_0}{2}}{n}^2\hfill \\ \hfill & +{\beta }_2({\psi }_1^{\ast }{\psi }_0^2{\psi }_{-1}^{\ast }+{\psi }_1{\psi }_0^{\ast 2}{\psi }_{-1})+{\displaystyle \frac{{\beta }_2}{2}}[|{\psi }_1{|}^4+|{\psi }_{-1}{|}^4\hfill \\ \hfill & +2|{\psi }_0{|}^2(|{\psi }_1{|}^2+|{\psi }_{-1}{|}^2)-2|{\psi }_1{|}^2|{\psi }_{-1}{|}^2]\}\mathrm{d}x\hfill \end{array}\end{array}

It can be seen from above that all the physical parameters have been absorbed in the nonlinear parameter {\beta }_0 and {\beta }_2, which can be varied by tuning the trap frequency and they are proportional to the number of atoms. Calculations with the same {\beta }_0 and {\beta }_2 can correspond to results for different species, but in diverse experimental conditions. We make the transformation {\psi }_{\alpha }(x,t)=\widetilde{{\psi }_{\alpha }}(x){\mathrm{e}}^{-\mathrm{i}{\mu }_{\alpha }t} to seek the stationary-state solutions to Eq. (5). Without loss of generality, we assume that \widetilde{{\psi }_{\alpha }}(x)(\alpha =0,\pm 1) are all real-valued. For simplicity, we write the equations in matrix form

\begin{array}{}\text{(10)}& {\widehat{H}}_0\psi -\mu \psi +{\beta }_0{\mathcal{A}}_1(\psi )\psi +{\beta }_2{\mathcal{A}}_2(\psi )\psi =0\end{array}

where \psi ={({\psi }_1,{\psi }_0,{\psi }_{-1})}^{\mathrm{T}}(the tilde is already omitted), \mu =\mathrm{diag}({\mu }_1,{\mu }_0,{\mu }_{-1}), {\mu }_1,{\mu }_0,{\mu }_{-1} are the chemical potentials, which satisfy

\begin{array}{}\text{(11)}& {\mu }_0=({\mu }_1+{\mu }_{-1})/2\end{array}

{\mathcal{A}}_1(\psi )=\mathrm{diag}(n,n,n), n={\sum }_{\alpha =-1}^1{\psi }_{\alpha }^2, {\mathcal{A}}_2(\psi ) is defined as

\begin{array}{}\text{(12)}& {\mathcal{A}}_2(\psi )=\left(\begin{array}{ccc}{\psi }_1^2+{\psi }_0^2-{\psi }_{-1}^2& 0& {\psi }_0^2\\ 0& {\psi }_1^2+{\psi }_{-1}^2+2{\psi }_1{\psi }_{-1}& 0\\ {\psi }_0^2& 0& {\psi }_{-1}^2+{\psi }_0^2-{\psi }_1^2\end{array}\right)\end{array}

Because of the degeneracy of the spin-1 BEC system, the solution of Eq. (10) is probably not unique.

The theoretical studies of the condensate have been carried out in several papers since the achievement of it in experiments [12–23]. However, to our best knowledge, there is little works on seeking the analytical solutions of Eq. (10), which has a harmonic potential and the nonlinear interaction coefficients are constants. In fact, it is a rather difficult task to find the exact solutions of Eq. (10). In this paper, we apply the homotopy analysis method (HAM) [24, 25] and the Galerkin spectral method to find the approximated analytical solutions of Eq. (10). The HAM is an analytic approximation method presented by Liao in 1992 [26] for solving nonlinear problems. By now, it has been widely applied to solve different types of nonlinear problems in science and engineering [24, 25, 27–42]. We also obtain the numerical results for the ground state solutions to Eq. (5) by the CNGF method [20]. Comparisons between the analytical results and the numerical results also have been done. The results show that they are in agreement well with each other when the nonlinear interactions are not too strong. Furthermore, we successfully find a class of exact solutions to Eq. (10) under certain conditions, which implies that the HAM probably can be used to seek the exact solutions for some nonlinear evolution equations.

2 The Homotopy Analysis Method

The HAM is independent of any small/large physical parameters at all: one can always transfer a nonlinear problem into an infinite number of linear subproblems by means of the HAM. It also provides us a convenient way to guarantee the convergence of solution series so that it is valid even if nonlinearity becomes rather strong. For interested readers, we refer to Refs. [24, 25] and relative references for more details. As is done in the standard HAM, we construct a homotopy

\begin{array}{}\text{(13)}& \begin{array}{cc}(1-q)({\widehat{H}}_0\Phi -\Lambda \Phi )=\hfill & qp[{\widehat{H}}_0\Phi -\Lambda \Phi +{\beta }_0{\mathcal{A}}_1(\Phi )\Phi \hfill \\ \hfill & +{\beta }_2{\mathcal{A}}_2(\Phi )\Phi ]\hfill \end{array}\end{array}

where \Phi =\Phi (x,q)={(F_1,F_0,F_{-1})}^{\mathrm{T}}, \Lambda =\Lambda (q)=\mathrm{diag}({\lambda }_1,{\lambda }_0,{\lambda }_{-1}), F_{\alpha }=F_{\alpha }(x,q),{\lambda }_{\alpha }={\lambda }_{\alpha }(q)(\alpha =0,\pm 1), q\in [0,1] is the embedding parameter and p is an auxiliary parameter (an adjustable constant), called the convergence-control parameter. When q=0, Eq. (13) becomes

\begin{array}{}\text{(14)}& {\widehat{H}}_0\Phi (x,0)=\Lambda (0)\Phi (x,0)\end{array}

Eq. (14) is a set of linear Schrödinger equations with harmonic potential for stationary state, thus \Phi (x,0) can be easily obtained. When q=1, Eq. (13) suggests that \psi (x)=\Phi (x,1) and \mu =\Lambda (1). It can be seen from Eq. (13) that as the parameter q increases from 0 to 1, the solution \Phi (x,q) varies from \Phi (x,0) to \psi (x), so does the \Lambda (q) from {\Lambda }_0\equiv \Lambda (0), the initial guess of μ, to μ. If this continuous variation is smooth enough, the Maclaurin’s series with respect to q can be constructed for \Phi (x,q) and \Lambda (q) respectively, and further, if all the series are convergent at q=1, we then have

\begin{array}{}\text{(15)}& \begin{array}{c}\psi (x)={\Phi (x,q)|}_{q=1}={\Phi }_0(x)+{\displaystyle \sum _{k=1}^{+\mathrm{\infty }}}{\Phi }_k(x)\hfill \\ \mu ={\Lambda (q)|}_{q=1}={\Lambda }_0+{\displaystyle \sum _{k=1}^{+\mathrm{\infty }}}{\Lambda }_k\hfill \end{array}\end{array}

with

\begin{array}{}\text{(16)}& \begin{array}{cc}{\Phi }_k(x)\hfill & ={{\displaystyle \frac{1}{k!}}{\displaystyle \frac{{\partial }^k\Phi (x,q)}{\partial q^k}}|}_{q=0}\hfill \\ \hfill & ={(F_{1,k}(x),F_{0,k}(x),F_{-1,k}(x))}^{\mathrm{T}}\hfill \\ {\Lambda }_k=\hfill & {{\displaystyle \frac{1}{k!}}{\displaystyle \frac{{\mathrm{d}}^k\Lambda (q)}{\mathrm{d}{q}^k}}|}_{q=0}=\mathrm{diag}({\lambda }_{1,k},{\lambda }_{0,k},{\lambda }_{-1,k})\hfill \\ \hfill & k=0,1,2,\mathrm{\cdots }\hfill \end{array}\end{array}

being the k-th order deformation derivatives, where

\begin{array}{}\text{(17)}& \begin{array}{c}{F}_{\alpha ,k}(x)={{\displaystyle \frac{1}{k!}}{\displaystyle \frac{{\partial }^k{F}_{\alpha }(x,q)}{\partial q^k}}|}_{q=0}\hfill \\ {\lambda }_{\alpha ,k}={{\displaystyle \frac{1}{k!}}{\displaystyle \frac{{\mathrm{d}}^k{\lambda }_{\alpha }(q)}{\mathrm{d}{q}^k}}|}_{q=0},\alpha =0,\pm 1\hfill \end{array}\end{array}

Note that {\Phi }_k(x) is a matrix that has three rows and one column and {\Lambda }_k is a diagonal matrix that has three rows and three columns.

Differentiating Eq. (13) k times with respect to q then setting q=0 and finally dividing them by k!, we gain the k-th order deformation equation for {\Phi }_k(x)

\begin{array}{}\text{(18)}& \begin{array}{c}k=0:{\widehat{H}}_0{\Phi }_0={\Lambda }_0{\Phi }_0\hfill \\ k=1:{\widehat{H}}_0{\Phi }_1-{\Lambda }_0{\Phi }_1-{\Lambda }_1{\Phi }_0=p\left({\beta }_0{\mathcal{F}}_0+{\beta }_2{\mathcal{G}}_0\right)\hfill \\ k>1:{\widehat{H}}_0({\Phi }_k-{\Phi }_{k-1})+{\displaystyle \sum _{j=0}^{k-1}}{\Lambda }_j({\Phi }_{k-1-j}-{\Phi }_{k-j})-{\Lambda }_k{\Phi }_0\hfill \\ =p\left({\widehat{H}}_0{\Phi }_{k-1}-{\displaystyle \sum _{j=0}^{k-1}}{\Lambda }_j{\Phi }_{k-1-j}+{\beta }_0{\mathcal{F}}_{k-1}+{\beta }_2{\mathcal{G}}_{k-1}\right)\hfill \end{array}\end{array}

in which {\mathcal{F}}_k={({\mathcal{F}}_{1,k},{\mathcal{F}}_{0,k},{\mathcal{F}}_{-1,k})}^{\mathrm{T}}, {\mathcal{G}}_k=({\mathcal{G}}_{1,k}, {\mathcal{G}}_{0,k}, {\mathcal{G}}_{-1,k}){}^{\mathrm{T}}, where

\begin{array}{}\text{(19)}& \begin{array}{cc}{\mathcal{F}}_{\kappa ,k}=\hfill & {\displaystyle \sum _{j=0}^k}{\displaystyle \sum _{l=0}^{k-j}}{F}_{\kappa ,j}{\displaystyle \sum _{\alpha =-1}^1}{F}_{\alpha ,l}{F}_{\alpha ,k-j-l},\kappa =1,0,-1\hfill \\ {\mathcal{G}}_{1,k}=\hfill & {\displaystyle \sum _{j=0}^k}{\displaystyle \sum _{l=0}^{k-j}}{F}_{1,j}(F_{1,l}{F}_{1,k-j-l}+F_{0,l}{F}_{0,k-j-l}\hfill \\ \hfill & -F_{-1,l}{F}_{-1,k-j-l})+F{}_{-1,j}F{}_{0,l}F{}_{0,k-j-l}\hfill \\ {\mathcal{G}}_{0,k}=\hfill & {\displaystyle \sum _{j=0}^k}{\displaystyle \sum _{l=0}^{k-j}}{F}_{0,j}(F_{1,l}{F}_{1,k-j-l}+F_{-1,l}{F}_{-1,k-j-l}\hfill \\ \hfill & +\mathrm{\hspace{0.17em}2}{F}_{1,l}{F}_{-1,k-j-l})\hfill \\ {\mathcal{G}}_{-1,k}=\hfill & {\displaystyle \sum _{j=0}^k}{\displaystyle \sum _{l=0}^{k-j}}{F}_{-1,j}(F_{-1,l}{F}_{-1,k-j-l}+F_{0,l}{F}_{0,k-j-l}\hfill \\ \hfill & -F_{1,l}{F}_{1,k-j-l})+F{}_{1,j}F{}_{0,l}F{}_{0,k-j-l}\hfill \end{array}\end{array}

The system of Eq. (18) is a recursive one, which means we can get {\Phi }_1(x) once {\Phi }_0(x) is given; furthermore, we can get {\Phi }_2(x) when both {\Phi }_0(x) and {\Phi }_1(x) are derived, and so on. Unfortunately, only {\Phi }_0(x) can be obtained easily although all the equations are linearized. It is rather difficult to obtain the exact solutions for {\Phi }_k(x) when k>0.  Here, we use the Galerkin spectral method to find the approximate analytical solutions for {\psi }_k to Eq. (18). We note that a similar procedure is also employed in Refs. [41, 42], where this approach is named the spectral-homotopy analysis method.

For a given k, suppose that {\Phi }_j(x) and {\Lambda }_j(j=0,1,\mathrm{\cdots },k-1) have been derived and we now find {\Phi }_k(x) and {\Lambda }_k. We consider a basis of functions, \{{\varphi }_l(x)|l=0,1,\mathrm{\cdots },L\}, perhaps satisfying some set of boundary conditions over the interval (-\mathrm{\infty },+\mathrm{\infty }), which is complete enough for expanding the unknown function over that interval to sufficient accuracy. For simplicity we assume that these functions form an orthonormal basis for the space. It is a natural choice to take {\varphi }_l(x) as the eigenfunctions of the 1D harmonic oscillator

\begin{array}{}\text{(20)}& {\varphi }_l(x)={\mathrm{e}}^{-x^2/2}{H}_l(x)/\sqrt{2^{l}l!\sqrt{\pi }},l=0,1,\mathrm{\cdots },L\end{array}

where H_l(x) is the Hermite polynomial with order l. The basis functions \{{\varphi }_l|l=0,1,\mathrm{\cdots },L\} form an orthonormal set which satisfy

\begin{array}{}\text{(21)}& ⟨{\varphi }_l|{\varphi }_j⟩={\displaystyle {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}}{\varphi }_l(x){\varphi }_j(x)dx={\delta }_{lj}\end{array}

and

\begin{array}{}\text{(22)}& {\widehat{H}}_0{\varphi }_l(x)=(l+1/2){\varphi }_l(x)\end{array}

All the unknown functions {\Phi }_k(x)(k=0,1,\mathrm{\cdots }) now can be expressed as the linear combination of the basis functions

\begin{array}{}\text{(23)}& F_{\alpha ,k}(x)={\displaystyle \sum _{l=0}^L}{A}_{\alpha ,k,l}{\varphi }_l(x), \alpha =0,\pm 1;k=0,1,2,\mathrm{\cdots }\end{array}

where A_{\alpha ,k,l} are expansion coefficients, under which the boundary conditions F_{\alpha ,k}(x)=0 are satisfied automatically.

For the case of k=0, the first equation in Eq. (18) can be rewritten as

\begin{array}{}\text{(24)}& {\widehat{H}}_0{F}_{\alpha ,0}(x)={\lambda }_{\alpha ,0}{F}_{\alpha ,0}, \alpha =0,\pm 1\end{array}

from which we can choose F_{\alpha ,0}(x)=A_{\alpha ,0,0}{\varphi }_0(x) with A_{\alpha ,0,l}=0(l>0) and {\lambda }_{\alpha ,0}=\frac{1}{2} (\alpha =0,\pm 1).

For k=1, we have

\begin{array}{}\text{(25)}& \begin{array}{c}{A}_{\alpha ,1,l}(l+1/2)-{\lambda }_{\alpha ,0}{A}_{\alpha ,1,l}-{\lambda }_{\alpha ,1}{A}_{\alpha ,0,0}⟨{\varphi }_0|{\varphi }_l⟩\hfill \\ =p{C}_{\alpha ,0,0}⟨{\varphi }_0^3|{\varphi }_l⟩\hfill \\ \alpha =0,\pm 1;l=0,1,\mathrm{\cdots },L\hfill \end{array}\end{array}

from which one can get

\begin{array}{}\text{(26)}& {\lambda }_{\alpha ,1}=-{\displaystyle \frac{p{C}_{\alpha ,0,0}}{A_{\alpha ,0,0}}}⟨{\varphi }_0^3|{\varphi }_0⟩=-{\displaystyle \frac{p{C}_{\alpha ,0,0}}{A_{\alpha ,0,0}}}{\displaystyle \frac{1}{\sqrt{2\pi }}}\end{array}

and

\begin{array}{}\text{(27)}& \begin{array}{l}{A}_{\alpha ,1,l}=\frac{p{C}_{\alpha ,0,0}⟨{\varphi }_0^3\mid {\varphi }_l⟩}{l}\\ \begin{array}{cc}& \end{array} =\{\begin{array}{c}0,\begin{array}{cccccc}& & & & & \end{array}\begin{array}{ccc}\begin{array}{cc}& \end{array}& & l\text{ is odd}\end{array}\\ \frac{p{C}_{\alpha ,0,0}}{\sqrt{2\pi }}\frac{(-1{)}^{\frac{l}{2}}(l-1)!!}{l\sqrt{2^{l}l!}}\begin{array}{cc}& l\text{ is even}\end{array}\end{array}\\ \begin{array}{c}\begin{array}{cc}& \end{array} l=1,2,\cdots ,L\end{array}\end{array}\end{array}

where

\begin{array}{}\text{(28)}& \begin{array}{ccc}{C}_{1,0,0}\hfill & =\hfill & {\beta }_0{A}_{1,0,0}{\displaystyle \sum _{\alpha =-1}^1}{A}_{\alpha ,0,0}^2+{\beta }_2(A_{1,0,0}^2+A_{0,0,0}^2\hfill \\ \hfill & \hfill & -A_{-1,0,0}^2)A{}_{1,0,0}+\beta {}_{2}A{}_{-1,0,0}{A}_{0,0,0}^2\hfill \\ C_{0,0,0}\hfill & =\hfill & {\beta }_0{A}_{0,0,0}{\displaystyle \sum _{\alpha =-1}^1}{A}_{\alpha ,0,0}^2+{\beta }_2(A_{1,0,0}^2+A_{-1,0,0}^2\hfill \\ \hfill & \hfill & +2A_{1,0,0}{A}_{-1,0,0})A{}_{0,0,0}\hfill \\ C_{-1,0,0}\hfill & =\hfill & {\beta }_0{A}_{-1,0,0}{\displaystyle \sum _{\alpha =-1}^1}{A}_{\alpha ,0,0}^2+{\beta }_2(A_{-1,0,0}^2+A_{0,0,0}^2\hfill \\ \hfill & \hfill & -A_{1,0,0}^2)A{}_{-1,0,0}+\beta {}_{2}A{}_{1,0,0}{A}_{0,0,0}^2\hfill \end{array}\end{array}

The coefficient A_{\alpha ,1,0} can not be determined by now. Without loss of generality, we may adopt A_{\alpha ,1,0}=0 because they can be combined with A_{\alpha ,0,0}, which will be determined by the normalization condition (6), the magnetization (8) and the relation between {\mu }_1,{\mu }_0 and {\mu }_{-1}, i.e., Eq. (11), together. (Note that if the jth component of the wave functions is zero, then Eq. (11) is no more valid and it should be replaced by \int {|{\psi }_j|}^{2}dx=0.) For the same reason, we take A_{\alpha ,k,0}=0(k⩾1).

The similar way can be used to find {\Phi }_k(x) and {\Lambda }_k when k>1. As k or L increase, however, the computational quantity will increase rapidly. Fortunately, the procedure described above can be implemented easily with the aid of symbol calculation software, such as Maple or Mathematica. Finally, at the Kth order approximation we have the analytical solution

\begin{array}{}\text{(29)}& \begin{array}{c}\psi (x)\approx {\psi }^K(x)={\Phi }_0(x)+{\displaystyle \sum _{k=1}^K}{\Phi }_k(x)\hfill \\ \mu \approx {\mu }^K={\Lambda }_0+{\displaystyle \sum _{k=1}^K}{\Lambda }_k\hfill \end{array}\end{array}

which can be written in an explicit component form

\begin{array}{}\text{(30)}& \begin{array}{c}{\psi }_{\alpha }(x)\approx F_{\alpha ,0}(x)+{\displaystyle \sum _{k=1}^K}{F}_{\alpha ,k}(x)\hfill \\ {\mu }_{\alpha }\approx {\lambda }_{\alpha ,0}+{\displaystyle \sum _{k=1}^K}{\lambda }_{\alpha ,k}(x),\alpha =0,\pm 1\hfill \end{array}\end{array}

It is worth remarkable that the Gauss-Hermite quadrature, employed in Refs. [43, 44], is also used, which makes our method more efficient to implement by the computer algebraic software.

3 Results and discussion

To verify the accuracy of the HAM descried above and to explore the behavior of the spin-1 BEC system, we consider two systems that their stationary states can be obtained numerically. The two systems are composed by ⁸⁷Rb and ²³Na in the F=1 hyperfine state and have repulsive atomic interactions. The numerical results are obtained by the CNGF methods [20]. In practice, we take M=0.2 to perform the computation of HAM.

When the spin-1 BEC system is composed of ⁸⁷Rb atoms, then m\approx 1.445\times {10}^{-25} kg. The frequencies of the external magnetic potential can be adopted as {\omega }_x=2\pi \times 40 Hz, {\omega }_y={\omega }_z=2\pi \times 130 Hz [45]. Other parameters are taken to be \mathrm{\hslash }\approx 1.05457\times {10}^{-34} J\cdots, a_0=10a_{\mathrm{B}}\approx 0.52917 nm, a_2=2a_{\mathrm{B}}\approx 0.10583 nm, where a_{\mathrm{B}}\approx 0.052917 nm is the Bohr radius. Thus we have {\beta }_0\approx 9.420\times {10}^{-4}N and {\beta }_2\approx -5.383\times {10}^{-4}N, which indicates that {\beta }_0\approx 9.4 and {\beta }_2\approx -5.4 when N={10}^4. Under this case, it is obviously that the atom-atom interaction is repulsive and it increases as the total number of particles increase. On the other hand, the spin-exchange interaction is attractive, which can reduce the nonlinear interaction among atoms. If we use ²³Na as an example, then m\approx 3.8194\times {10}^{-26} kg. We take a_0=6a_B and a_2=10a_B, other parameters are adopted as before, then {\beta }_0\approx 8.994\times {10}^{-4}N and {\beta }_2\approx 1.384\times {10}^{-4}N, suggests that {\beta }_0\approx 9.0 and {\beta }_2\approx 1.4 when N={10}^4. Clearly, both the atom-atom interaction and the spin-exchange interaction are repulsive and the interaction increases as the total number of particles increase. The spin-exchange interaction can enhance or debase the nonlinear interaction among atoms.

In the frame of HAM, the solutions are expressed as a series of basis functions. Note that the solutions contains the adjustable parameter p, which provides us with a simply way to adjust and control the convergence of the solution series. In general, by means of the so-called p-curve, it is straightforward to choose an appropriate range for p which ensures the convergence of the series. In Refs. [31, 35], the authors proposed an optimal HAM to find the “best” p by minimizing the residual error. However, for our method described above, the computational quantity is terrible to perform such a calculation, especially for larger K and L. In this paper, we use two steps to find the values for p. Firstly, we find a rough value of p by observing the p-curve. Secondly, we adjust p to change the magnitude of the residual error. The value for p is accepted when the residual error is small enough. If the residual error can not be small enough, we have to increase the value of K to find the higher order approximation solutions. By this approach, the value of p is an “acceptable one” but not the “best one”. It should be remarkable that one of the difficulty in the HAM is to prove the convergence of the series. In Refs. [36–38], the author made some efforts to present a rigorous mathematical approach to show the convergence of the obtained series. It is very important to guarantee the convergence of homotopy-series. Once the series are convergent, the solutions obtained by the HAM has converged to the exact solutions. However, it is a pity that there do not exist any mathematical theorems which can guide us in details how to construct a good enough homotopy for any a given nonlinear equation. In theory, the convergence of homotopy-series is strongly dependent upon the initial approximation and the auxiliary linear operator as a whole [25].

In practice, we find that as the nonlinear interactions increase, the interval of convergence becomes smaller and smaller. Thus, our method described above fails to give a good result when the nonlinear interactions is large enough. This can be understood as following. When the nonlinear interactions become stronger, the interaction energy among atoms is much greater than the kinetic energy, which can be neglected under this case. Therefore, the auxiliary linear operator and the initial approximation we chosen can no longer be a good choice, which leads the present approach is no good for stronger interactions.

Figure 1 shows the numerical values of the energy E per particle and the chemical potential {\mu }_1, given by the CNGF method [20] with 2049 grids in the interval x,y\in [-6.4,6.4], and the analytical results, obtained by the HAM described above, for different values of {\beta }_0 and {\beta }_2. In Fig. 1(a) and (b), , while {\beta }_2>0 in Fig. 1(c) and (d). The values of p we choose can ensure that the series are convergent. For the case of , both the numerical results and the analytical results provide that {\mu }_0={\mu }_{-1}={\mu }_1. However, this is not necessarily the case if {\beta }_2>0. When both {\beta }_0 and {\beta }_2 are positive, the numerical results exhibit that the ground state solutions satisfy {\psi }_0(x)=0, which is confirmed by the analytical results. Under this case, {\mu }_0 is undefined and Eq. (11) is no more valid.

It can be seen from Fig. 1(a) and (b) that the analytical results are very close to the numerical results as {\beta }_0+{\beta }_2 is small. When {\beta }_0 ({\beta }_2) is fixed, the difference of the chemical potential between the analytical results and the numerical results increases as {\beta }_2 ({\beta }_0) increases. While it is not obvious that the error of energy E varies as the chemical potential does. In the cases of Fig. 1(c) and (d), we find that the chemical potential {\mu }_1 is no longer equal to {\mu }_{-1} for the ground state solutions. The values of {\mu }_1 and {\mu }_{-1} depend on the magnetization M. In Fig. 1(c), one can see that the analytical dimensionless energy is agreement very well with the numerical results. The difference between the chemical potential of the numerical results and the corresponding analytical results are also small (\sim 0.01) and the relative error is less than 0.1\%. For the case we considering in Fig. 1(c), the chemical potential {\mu }_1 ({\mu }_{-1}) approximately increases (decreases) linearly as {\beta }_2 increases. In Fig. 1(d), one can also see that the error increases as {\beta }_0 increases, just as shown in Fig. 1(b). The densities of the wave functions {\psi }_{\alpha }(x)(\alpha =0,\pm 1) given by our method and the corresponding numerical results, derived by the CNGF method, are shown in Fig. 2 for different values of {\beta }_0 and {\beta }_2. It also can be seen that the analytical results are in agreement very well with the numerical results when {\beta }_0+{\beta }_2 is not too large. We note that the density at the center of the trap is lower for bigger values of {\beta }_0+{\beta }_2, which is expected because of the corresponding higher inter-particle repulsion.

From Fig. 2, one can see clearly that the error, mainly induced by the truncation of the HAM and the spectral Galerkin method (the round-off error is small enough that can be neglected), increases as the value of {\beta }_0+{\beta }_2 increases. An effective way to decrease the error is increasing K when L is large enough, but the computational quantity will increase rapidly. To see this more transparent, Table 1 lists the difference between the numerical results and the analytical ones for different orders of the HAM approximations K, where the number of basis functions used in the spectral Galerkin method is fixed as L=8. In Table 1, we define ϵ=\underset{\genfrac{}{}{0pt}{}{j}{\alpha =0,\pm 1}}{\max }|{\psi }_{\alpha }^{\mathrm{N}}{(x_j)}^2-{\psi }_{\alpha }^{\mathrm{A}}{(x_j)}^2| and {ϵ}_r=\underset{\genfrac{}{}{0pt}{}{|{\psi }_{\alpha }^{\mathrm{N}}{(x_j)}^2|>0.01}{\alpha =0,\pm 1}}{\max }\frac{|{\psi }_{\alpha }^{\mathrm{N}}{(x_j)}^2-{\psi }_{\alpha }^{\mathrm{A}}{(x_j)}^2|}{|{\psi }_{\alpha }^{\mathrm{N}}{(x_j)}^2|} to measure the magnitude of density difference between the numerical results and analytical ones, the superscript N and A denotes numerical result and analytical result, respectively. ϵ and {ϵ}_r can be regarded as the maximum absolute error and relative error of the particle densities. It can be seen from Table 1 that the error decreases as the order of the HAM approximations increases. One can also see that the error increases as the nonlinear interactions {\beta }_0+{\beta }_2 increases, just as we can see from Figs. 1 and 2. When the atomic interactions is strong, i.e., {\beta }_0+{\beta }_2 is large enough, we have to take larger K and L to obtain an appropriated result, which will take a long time to perform such a calculation.

When {\beta }_0+{\beta }_2=0, the HAM provides a class of exact solutions to Eq. (10), which read

\begin{array}{}\text{(31)}& \begin{array}{c}{\psi }_1(x)={\displaystyle \frac{1+M}{2}}{\varphi }_0(x), {\psi }_0(x)=\sqrt{{\displaystyle \frac{1-M^2}{2}}}{\varphi }_0(x)\hfill \\ {\psi }_{-1}(x)={\displaystyle \frac{1-M}{2}}{\varphi }_0(x)\hfill \end{array}\end{array}

with {\mu }_1={\mu }_0={\mu }_{-1}=1/2. Solution (31) is the ground state solution with the Gaussian profiles for a set of linear Schrödinger equations that has a harmonic potential. One can note that the coefficients of atom-atom interaction and spin-exchange interaction, i.e., {\beta }_0 and {\beta }_2, are in absence of the expressions. The reason is that the atomic interaction are exactly balanced by the spin-exchange interaction when {\beta }_0+{\beta }_2=0. Heuristically, under this condition, we also obtain a set of excited state solutions to Eq. (10)

\begin{array}{}\text{(32)}& \begin{array}{c}{\psi }_1(x)={\displaystyle \frac{1+M}{2}}{\varphi }_{\nu }(x), {\psi }_0(x)=\sqrt{{\displaystyle \frac{1-M^2}{2}}}{\varphi }_{\nu }(x)\hfill \\ {\psi }_{-1}(x)={\displaystyle \frac{1-M}{2}}{\varphi }_{\nu }(x)\hfill \end{array}\end{array}

with {\mu }_1={\mu }_0={\mu }_{-1}=\nu +1/2(\nu =1,2,3,\mathrm{\cdots }), where ν is a quantum number. Solutions (32) cannot be obtained by the CNGF method because this approach is designed to seek the ground state solutions of the spin-1 system.

The HAM is usually used to find the approximate solutions for nonlinear problems. But now, we find that it may be applied to seek the exact solutions for some nonlinear evolution equations.

4 Conclusion

In summary, an analytical approach based on the homotopy analysis and Galerkin method is adopted to find the approximate analytical solutions of a set of Gross–Pitaevskii equations, which is used to describe the dynamical behavior, based on the mean-field theory, of spin-1 BEC in a harmonic potential. Comparisons between the analytical solutions we got and the numerical results, obtained by the CNGF method, have been made. The results indicate that our results are agreement very well with the numerical results if the nonlinear interaction is not too large. We also obtain a class of exact solutions to the GPEs with harmonic potential for a special case, which implies that the HAM may be applied to find a class of exact solutions for some nonlinear evolution equations. The method used in this paper can be generalized to solve many nonlinear mathematical physics problems in various fields.

References

[1] J. Stenger, S. Inouye, D. M. Stamper-Kurn, H. J. Miesner, A. P. Chikkatur, and W. Ketterle, Nature, 1998, 396(6709): 345

[2] D. M. Stamper-Kurn, M. R. Andrews, A. P. Chikkatur, S. Inouye, H. J. Miesner, J. Stenger, and W. Ketterle, Phys. Rev. Lett., 1998, 80(10): 2027

[3] M. S. Chang, C. D. Hamley, M. D. Barrett, J. A. Sauer, K. M. Fortier, W. Zhang, L. You, and M. S. Chapman, Phys. Rev. Lett., 2004, 92(14): 140403

[4] A. Sinatra, J. C. Dornstetter, and Y. Castin, Front. Phys., 2012, 7(1): 86

[5] X. W. Guan, Front. Phys., 2012, 7(1): 8

[6] T. L. Ho, Phys. Rev. Lett., 1998, 81(4): 742

[7] T. Ohmi and K.  Machida, J. Phys. Soc. Jpn., 1998, 67(6): 1822

[8] T. Isoshima, K. Machida, and T. Ohmi, Phys. Rev. A, 1999, 60(6): 4857

[9] P. G. Kevrekidis, G. Theocharis, D. J. Frantzeskakis, and B. A. Malomed, Phys. Rev. Lett., 2003, 90(23): 230401

[10] N. N. Klausen, J. L. Bohn, and C. H. Greene, Phys. Rev. A, 2001, 64(5): 053602

[11] E. G. M. van Kempen, S. J. J. M. F. Kokkelmans, D. J. Heinzen, and B. J. Verhaar, Phys. Rev. Lett., 2002, 88(9): 093201

[12] H. E. Nistazakis, D. J. Frantzeskakis, P. G. Kevrekidis, B. A. Malomed, and R. Carretero-Gonzlez, Phys. Rev. A, 2008, 77(3): 033612

[13] E. J. Mueller, Phys. Rev. A, 2004, 69(3): 033606

[14] J. Ieda, T. Miyakawa, and M. Wadati, Laser Phys., 2005, 16(4): 678

[15] B. J. Dabrowska-Wüster, E. A. Ostrovskaya, T. J. Alexander, and Y. S. Kivshar, Phys. Rev. A, 2007, 75(2): 023617

[16] P. Szańkowski, M. Trippenbach, E. Infeld, and G. Rowlands, Phys. Rev. A, 2011, 83: 013626

[17] R. Barnett, A. Turner, and E. Demler, Phys. Rev. Lett., 2006, 97(18): 180412

[18] H. Saito and M.  Ueda, Phys. Rev. A, 2005, 72(2): 023610

[19] W. Z. Bao and Y. Z.  Zhang, Methods Appl. Anal., 2010, 17: 49

[20] W. Z. Bao and F. Y.  Lim, SIAM J. Sci. Comput., 2008, 30(4): 1925

[21] D. S. Wang, S. W. Song, B. Xiong, and W. M. Liu, Phys. Rev. A, 2011, 84(5): 053607

[22] D. S. Wang, X. H. Hu, and W. M. Liu, Phys. Rev. A, 2010, 82(2): 023612

[23] D. S. Wang, X. H. Hu, J. P. Hu, and W. M. Liu, Phys. Rev. A, 2010, 81(2): 025604

[24] S. J. Liao, Beyond Perturbation: Introduction to Homotopy Analysis Method, Boca Raton: Chapman & Hall/CRC Press, 2003

[25] S. J. Liao, Homotopy Analysis Method in Nonlinear Differential Equations, Beijing: Springer & Higher Education Press, 2012

[26] S. J. Liao, PhD Thesis, Shanghai Jiao Tong University, 1992

[27] S. J. Liao, Commun. Nonlinear Sci. Numer. Simul., 2010, 15(8): 2003

[28] S. J. Liao, Commun. Nonlinear Sci. Numer. Simul., 2009, 14(4): 983

[29] Y. Cheng and S. J.  Liao, Commun. Nonlinear Sci. Numer. Simul., 2006, 47: 75

[30] S. J. Liao, Int. J. Non-linear Mech., 2004, 39(2): 271

[31] S. J. Liao, Commun. Nonlinear Sci. Numer. Simul., 2010, 15(8): 2003

[32] Y. L. Zhao, Z. L. Lin, Z. Liu, and S. J. Liao, Appl. Math. Comput., 2012, 218(17): 8363

[33] M. Turkyilmazoglu, Appl. Math. Lett., 2010, 23(10): 1226

[34] M. Turkyilmazoglu, Int. J. Mech. Sci., 2010, 52(12): 1735

[35] M. Turkyilmazoglu, Comput. Fluids, 2010, 39(5): 793

[36] M. Turkyilmazoglu, Int. J. Therm. Sci., 2011, 50(5): 831

[37] M. Turkyilmazoglu, Int. J. Nonlinear Sci. Numer. Simul., 2011, 12(1–8): 9

[38] M. Turkyilmazoglu, Commun. Nonlinear Sci. Numer. Simul., 2012, 17(11): 4097

[39] Y. R. Shi, X. J. Xu, Z. X. Wu, Y. H. Wang, H. J. Yang, W. S. Duan, and K. P. Lü, Acta Phys. Sin., 2006, 55: 1555 (in Chinese)

[40] Y. R. Shi and H. J.  Yang, Acta Phys. Sin., 2010, 59: 67 (in Chinese)

[41] S. S. Motsa, P. Sibanda, and S. Shateyi, Commun. Nonl. Sci. Numer. Simulat., 2010, 15(9): 2293

[42] P. Sibanda, S. S. Motsa, and Z. Makukula, Int. J. Numer. Method. H., 2012, 22(1): 4

[43] C. M. Dion and E.  Cancès, Phys. Rev. E, 2003, 67(4): 046706

[44] C. M. Dion and E.  Cancès, Comput. Phys. Commun., 2007, 177(10): 787

[45] S. B. Papp, J. M. Pino, and C. E. Wieman, Phys. Rev. Lett., 2008, 101(4): 040402

Rights & permissions

© Higher Education Press and Springer-Verlag Berlin Heidelberg 2013