旋转的双原子–玻色爱因斯坦凝聚体基态数值模拟

Numerical Simulations on Ground States for Rotating Two-Component Bose-Einstein Condensates

摘要: 静态相耦合的Gross-Pitaevskii方程组的解描述了双原子的玻色爱因斯坦凝聚体在极低温度下的基态现象。我们提出一种十分有效的数值方法——梯度法来求解此基态解。我们提出的梯度法在数值上既保持总模量守恒又能使总能量递减；我们严格地证明我们提出的梯度法是一种获得能量函数在给定限制性条件下的最小值(也即基态解)的十分有效的方法。我们通过大量的例子来显示该方法的优点并且应用该方法去研究处于旋转状态下的双原子–玻色爱因斯坦凝聚体在极低温度下所呈现的复杂涡旋现象。

Abstract: The ground states of rotating two-component Bose-Einstein condensates (BEC) at extremely low temperature are solutions of time-independent coupled Gross-Pitaevskii equations. To compute the ground states, we propose an efficient numerical method—gradient flows with discrete nor-malization. In linear cases and under properly chosen initial data, we can prove that the gradient flows converge into the ground state which has the lowest energy. We show that the method is quite efficient and apply the method to study complicated vortex structure in the ground state solutions of rotating two-component BEC at extremely low temperature.

1. 引言

自从在双原子–玻色爱因斯坦凝聚体的实验中发现涡旋现象后 [1] - [6] ，人们就对旋转状态下的双原子–玻色爱因斯坦凝聚体所具有的一种新的现象——涡旋现象产生了浓厚的兴趣。根据平均场近似理论，静态且相耦合的Gross-Pitaevskii (GP)方程组被用来描述旋转中的双原子–玻色爱因斯坦凝聚体的基态现象 [7] [8] [9] [10] [11] 。事实上，许多关于无旋转的双原子–玻色爱因斯坦凝聚体的基态现象 [12] [13] [14] [15] ，如自然对称现象的破裂 [16] 、中心对称的涡旋 [17] [18] [19] [20] 以及稳定的skyrmions [21] 等揭示了许多实验现象，并解释了许多实验难以解决的问题。

最近我们在文 [22] 中研究了旋转中的双原子–玻色爱因斯坦凝聚体的基态现象，但是我们并没有考虑由外部力场引起的Josephson耦合力。在Josephson耦合力作用下的双原子–玻色爱因斯坦凝聚体的基态现象会不同于没有Josephson耦合力作用的双原子–玻色爱因斯坦凝聚体的基态现象。对于这一问题，就我们所知，只有较少人研究过。我们撰写本文的主要目的是研究在Josephson耦合力作用下的双原子–玻色爱因斯坦凝聚体的基态现象。该基态现象可以用静态且相耦合的GP方程组来描述。为此，我们首先提出一种十分有效的数值方法——梯度法来求解耦合的静态GP方程组。我们提出的梯度法在数值上既保持总模量守恒又能使总能量递减，从而保证用该方法求得能量函数在给定限制性条件下的最小值。我们严格地证明我们提出的梯度法是一种获得能量函数在给定限制性条件下的最小值的十分有效的方法。

本篇论文是按照如下方式组织的。在第二节中我们将介绍用来描述双原子–玻色爱因斯坦凝聚体的相耦合的GP方程组、如何把它们化为无量纲形式的方程组以及如何把三维方程组简化为二维方程组。在第三节中我们详细地讨论如何用梯度法求得相耦合的GP方程组的数值解，我们给出详细的算法并在数学上严格地证明我们提出的方法的有效性。在第四节中我们利用提出的数值方法研究旋转中的双原子–玻色爱因斯坦凝聚体在极低温度下所呈现的涡旋现象。在最后一节，也就是第五节，我们将探讨得到我们的结论以及我们提出的新方法可能用到的领域。

2. 动态相耦合的GP方程组

我们考虑处于旋转状态下的双原子–玻色爱因斯坦凝聚体(比如铷87原子及其同位素)。用来描述该双原子–玻色爱因斯坦凝聚体的数学模型便是以下相耦合的GP方程组 [23]

$i\hslash \frac{\partial {\psi }_1}{\partial t}=\left(-\frac{{\hslash }^2}{2m}{\nabla }^2+V_1\left(x\right)-\hslash \Omega L_z+g_{11}{\left|{\psi }_1\right|}^2+g_{12}{\left|{\psi }_2\right|}^2\right){\psi }_1-\hslash {\omega }_R{\psi }_2$ (2.1)

$i\hslash \frac{\partial {\psi }_2}{\partial t}=\left(-\frac{{\hslash }^2}{2m}{\nabla }^2+V_2\left(x\right)-\hslash \Omega L_z+g_{21}{\left|{\psi }_1\right|}^2+g_{22}{\left|{\psi }_2\right|}^2\right){\psi }_2-\hslash {\omega }_R{\psi }_1$ (2.2)

这里 ${\psi }_j={\psi }_j\left(x,t\right),j=1,2\left(x=\left(x,y,z\right)\in R^3\right)$ 是描述双原子–玻色爱因斯坦凝聚体的两个玻动函数。 $m$ 是⁸⁷Rb原子的质量，调和磁势井函数 $V_j\left(x\right)\left(j=1,2\right)$ 用

$V_j\left(x\right)=\frac{m}{2}\left({\omega }_{x,j}^2{\left(x-x_j\right)}^2+{\omega }_{y,j}^2{\left(y-y_j\right)}^2+{\omega }_{z,j}^2{\left(z-z_j\right)}^2\right)$

来表示。 ${\alpha }_j\left(\alpha =x,y,z\right)$ 是 $\alpha$ 轴在上磁势井的中心。角动量转动项为 $\Omega L_z=-i\Omega \left(x{\partial }_y-y{\partial }_x\right)$ 。双原子内部以及双原子之间的键力分别用 $g_{jj}=4\text{π}{\hslash }^2{a}_{jj}/m\left(j=1,2\right)$ 和 $g_{jl}=4\text{π}{\hslash }^2{a}_{jl}/m\left(j\ne l=1,2\right)$ 来表示。方程组(2.1) (2.2)中的最后一项表示由外部力场引起的Josephson耦合力，并且这里的 ${\omega }_R$ 是Rabi频率。该相耦合的GP方程组(2.1)~(2.2)保持系统的全部粒子数

$N:=N\left({\psi }_1,{\psi }_2\right)={\displaystyle \underset{j=1}{\overset{2}{\sum }}{\displaystyle \underset{R^3}{\int }{\left|{\psi }_j\left(x,t\right)\right|}^2\text{d}x}}=N_1+N_2$ (2.3)

不变。这里 $N_j={\displaystyle \underset{R^3}{\int }{\left|{\psi }_j\left(x,t\right)\right|}^2\text{d}x}$ 是双原子之一的粒子数， $N=N_1+N_2$ 是玻色爱因斯坦双原子总粒子数。该相耦合的GP方程组(2.1)~(2.2)也保持系统的能量

$\begin{array}{c}E\left({\psi }_1,{\psi }_2\right)={\displaystyle \underset{R^3}{\int }{\displaystyle \underset{j=1}{\overset{2}{\sum }}\left[\frac{{\hslash }^2}{2m}{\left|\nabla {\psi }_j\right|}^2+V_j{\left|{\psi }_j\right|}^2-\Omega {\bar{\psi }}_j{L}_z{\psi }_j+\frac{1}{2}\left(g_{j1}{\left|{\psi }_1\right|}^2+g_{j2}{\left|{\psi }_2\right|}^2\right){\left|{\psi }_j\right|}^2\text{d}x\right]}}\\ -{\displaystyle \underset{R^3}{\int }\hslash {\omega }_R\left({\bar{\psi }}_2{\psi }_1+{\bar{\psi }}_1{\psi }_2\right)\text{d}x}\end{array}$

不变。

2.1. 无量纲形式的GP方程组

GP方程组(2.1)~(2.2)不大方便直接用于数值模拟。如果我们引入：

$\tilde{t}={\omega }_{x,1}t$ , $\tilde{x}=\frac{x}{a_0}$ , ${\tilde{\psi }}_j\left(\tilde{x},\tilde{t}\right)=\frac{a_0^{3/2}}{\sqrt{N}}{\psi }_j\left(x,t\right)$ , $\tilde{\Omega }=\frac{\Omega }{{\omega }_{x,1}}$ , $a_0=\sqrt{\frac{\hslash }{m{\omega }_{x,1}}}$ ，

$\tilde{\lambda }=\frac{{\omega }_R}{{\omega }_{x,1}}$ 和 ${\tilde{\alpha }}_j=\frac{{\alpha }_j}{a_0}$ ( $\alpha =x,y,z$ )。 (2.4)

事实上这里我们取 $1/{\omega }_{x,1}$ 和 $a_0$ 分别作为时间及长度的无量纲单位。把参数(2.4)代入方程组(2.1)~(2.2)中，并且在方程组(2.1)~(2.2)的两边都乘以 $\frac{1}{m{\omega }_{x,1}^2{\left(N{a}_0\right)}^{1/2}}$ 。如果我们再去掉所有的~，则我们可以得到以下的三维的无量纲形式的GP方程组：

$i\frac{\partial }{\partial t}{\psi }_j\left(x,t\right)=\left(-\frac{1}{2}{\nabla }^2-\Omega L_z+V_{2,j}\left(x\right)+\underset{l=1}{\overset{2}{\sum }}{\beta }_{2,jl}{\left|{\psi }_l\right|}^2\right){\psi }_j\left(x,t\right)-\lambda {\psi }_k\left(x,t\right),j=1,2$ (2.5)

并且

$k=2$ ，当 $j=1$ ； $k=1$ 当 $j=2$ ，

$V_j\left(x\right)=\frac{1}{2}\left({\gamma }_{x,j}^2{\left(x-x_j\right)}^2+{\gamma }_{y,j}^2{\left(y-y_j\right)}^2+{\gamma }_{z,j}^2{\left(z-z_j\right)}^2\right)$ ，

${\gamma }_{x,j}=\frac{{\omega }_{x,j}}{{\omega }_{x,1}}$ ， ${\gamma }_{y,j}=\frac{{\omega }_{y,j}}{{\omega }_{x,1}}$ ， ${\gamma }_{z,j}=\frac{{\omega }_{z,j}}{{\omega }_{x,1}}$ (2.6)

${\beta }_{jl}=\frac{g_{jl}N}{a_0^3\hslash {\omega }_{x,1}}=\frac{4\text{π}{a}_{jl}N}{a_0},j,l=1,2$ 。

系统每个粒子的能量成为

$E_{\beta }\left({\psi }_1,{\psi }_2\right)={\displaystyle \underset{j=1}{\overset{2}{\sum }}{E}_{\beta ,j}\left({\psi }_1,{\psi }_2\right)}-\lambda {\displaystyle \underset{R^3}{\int }\left({\bar{\psi }}_1{\psi }_2+{\psi }_1{\bar{\psi }}_2\right)\text{d}x}$ (2.7)

这里

$\beta ={\max }_{1\le j,l\le 2}\left|{\beta }_{jl}\right|$

$\begin{array}{c}{E}_{\beta }\left({\psi }_1,{\psi }_2\right)={\displaystyle \underset{l=1}{\overset{2}{\sum }}{\displaystyle \underset{R^3}{\int }\left[\frac{1}{2}{\left|\nabla {\psi }_j\right|}^2+V_j{\left|{\psi }_j\right|}^2-\Omega {\bar{\psi }}_j{L}_z{\psi }_j+\frac{1}{2}{\displaystyle \underset{l=1}{\overset{2}{\sum }}{\beta }_{jl}{\left|{\psi }_l\right|}^2{\left|{\psi }_j\right|}^2}\right]\text{d}x}}\\ -\lambda {\displaystyle \underset{R^3}{\int }\left({\bar{\psi }}_1{\psi }_2+{\psi }_1{\bar{\psi }}_2\right)\text{d}x}\end{array}$

2.2. 二维的方程组

方程组(2.5)也可近似地变为二维形式的方程组。在磁势井是圆盘形式时，我们可假设 ${\psi }_j\left(x,y,z,t\right)\to {\psi }_j\left(x,y,t\right){\psi }_{ho}\left(z\right)$ ，这里的 ${\psi }_{ho}\left(z\right)=\frac{{\gamma }_{z,j}^{1/4}}{{\text{π}}^{1/4}}{\text{e}}^{-\frac{{\gamma }_{z,j}{\left(z-z_j\right)}^2}{2}}$ 。将之代入方程组(2.5)，我们可把偶合的GP方程组(2.5)转化成一个二维的无量纲形式的GP方程组。

$i\frac{\partial }{\partial t}{\psi }_j\left(x,t\right)=\left(-\frac{1}{2}{\nabla }^2-\Omega L_z+V_{2,j}\left(x\right)+{\displaystyle \underset{l=1}{\overset{2}{\sum }}{\beta }_{2,jl}{\left|{\psi }_l\right|}^2}\right){\psi }_j\left(x,t\right)-\lambda {\psi }_k\left(x,t\right)$

这里 $k=2$ ，当 $j=1$ ； $k=1$ ，当 $j=2$ 。 $V_{2,j}\left(x\right)=\frac{m_j}{2m_1}\left({\gamma }_{x,j}^2{\left(x-x_j\right)}^2+{\gamma }_{y,j}^2{\left(y-y_j\right)}^2\right)$ 。

总之，我们将讨论以下无量纲形式且相耦合的GP方程组，

$i\frac{\partial }{\partial t}{\psi }_j\left(x,t\right)=\left(-\frac{1}{2}{\nabla }^2-\Omega L_z+V_{d,j}\left(x\right)+{\displaystyle \underset{i=1}{\overset{2}{\sum }}{\beta }_{d,jl}{\left|{\phi }_l\right|}^2}\right){\psi }_j\left(x,t\right)-\lambda {\psi }_k$ (2.9)

这里 $d=2$ 或者3， ${\beta }_{3,ji}={\beta }_{jl}$ ， $V_{3,j}\left(x\right)=V_j\left(x\right)$ 。

2.3. GP方程组(2.9)的一些性质

为方便起见，本节我们将讨论GP方程组(2.9)的一些解析性质。

当 $\lambda \ne 0$ 时，(2.9)式中两个重要的守恒量。一个是总模量

$N\equiv N\left({\psi }_1,{\psi }_2\right)={\displaystyle \underset{j=1}{\overset{2}{\sum }}{{\displaystyle {\int }_{R^d}\left|{\psi }_j\left(x,t\right)\right|}}^2}\text{d}x=1,t\ge 0,j=1,2$ (2.10)

守恒；另一个是总能量

$E\equiv E_{{\beta }_d}\left({\phi }_1,{\phi }_2\right)={\displaystyle \underset{j=1}{\overset{2}{\sum }}{E}_{{\beta }_d,j}}\left({\phi }_1,{\phi }_2\right)-\lambda {\displaystyle \int \left({\bar{\phi }}_1{\phi }_2+{\bar{\phi }}_2{\phi }_1\right)\text{d}x}$ (2.11)

守恒，这里

${\beta }_d=\underset{1\le j,l\le 2}{\max }\left|{\beta }_{d,jl}\right|$ ，

$E_{{\beta }_d,j}\left({\phi }_1,{\phi }_2\right)={\displaystyle {\int }_{R^d}\left[\frac{1}{2}{\left|\nabla {\phi }_j\right|}^2+V_{d,j}\left(x\right){\left|{\phi }_j\right|}^2+\frac{1}{2}{\displaystyle \underset{l=1}{\overset{2}{\sum }}{\beta }_{d,jl}{\left|{\phi }_l\right|}^2{\left|{\phi }_j\right|}^2-\Omega {\bar{\phi }}_j{L}_z{\phi }_j}\right]\text{d}x}$ 。

3. 静态相耦合的GP方程组

为了寻求(2.9)式的一种静态解(它们通常被用来描述旋转中的双原子–玻色爱因斯坦凝聚体的平衡态现象)，为此，我们引入

${\psi }_j\left(x,t\right)={\text{e}}^{-i\mu t}{\varphi }_j\left(x\right),j=1,2$ (3.1)

这里 $\mu$ 是这两种粒子之间的化学势能， ${\varphi }_j$ 是与时间无关的函数。函数 ${\varphi }_j\left(x\right)\left(j=1,2\right)$ 被用来描述旋转中的双原子–玻色爱因斯坦凝聚体的平衡态现象。

将(3.1)带入(2.9)式得到以下只与 ${\varphi }_j\left(x\right)$ 有关且与时间无关的GP方程组

$\mu {\varphi }_1\left(x\right)=\left(-\frac{1}{2}{\nabla }^2+V_{d,1}\left(x\right)+{\displaystyle \underset{l=1`}{\overset{2}{\sum }}{\beta }_{d,1l}{\left|{\varphi }_l\right|}^2}-\Omega L_z\right){\varphi }_1\left(x\right)-\lambda {\varphi }_2\left(x\right)$ (3.2)

$\mu {\varphi }_2\left(x\right)=\left(-\frac{1}{2}{\nabla }^2+V_{d,2}\left(x\right)+{\displaystyle \underset{l=1`}{\overset{2}{\sum }}{\beta }_{d,2l}{\left|{\varphi }_l\right|}^2}-\Omega L_z\right){\varphi }_2\left(x\right)-\lambda {\varphi }_1\left(x\right)$ (3.3)

当对这两部分进行标准化时，我们需要总模量

${\displaystyle \underset{j=1}{\overset{2}{\sum }}{\displaystyle {\int }_{R^d}{\left|{\varphi }_j\left(x\right)\right|}^2\text{d}x}}=1$ (3.4)

我们知道化学势能

$\mu =E_{{\beta }_d}\left({\varphi }_1,{\varphi }_2\right)+\frac{1}{2}{\displaystyle {\int }_{R^d}{\displaystyle \underset{j,l=1}{\overset{2}{\sum }}{\beta }_{d,jl}}}{\left|{\varphi }_j\left(x\right)\right|}^2{\left|{\varphi }_l\left(x\right)\right|}^2\text{d}x$ (3.5)

并且，双原子–玻色爱因斯坦凝聚体的能量可以表示为：

$E_{{\beta }_d}\left({\varphi }_1,{\varphi }_2\right)={\displaystyle \underset{j=1}{\overset{2}{\sum }}{E}_{{\beta }_d,j}}\left({\varphi }_1,{\varphi }_2\right)-\lambda {\displaystyle \int \left({\bar{\varphi }}_1{\varphi }_2+{\bar{\varphi }}_2{\varphi }_1\right)\text{d}x}$ (3.6)

函数 ${\varphi }_j\left(x\right)\left(j=1,2\right)$ 所满足的静态相耦合的GP方程组(3.2)~(3.4)实际是一个带有约束条件的非线性特征函数与特征值问题。静态相耦合的GP方程组(3.2)~(3.4)有许多解，

人们最关心的是下面的基态解。

3.1. 基态解的定义

旋转中的双原子玻色爱因斯坦凝聚体的基态解 ${\Phi }_g\left(x\right)=\left({\varphi }_{g,1}\left(x\right),{\varphi }_{g,2}\left(x\right)\right)$ 实际上是能量函数 $E_{{\beta }_d}\left({\varphi }_1,{\varphi }_2\right)$ 在限制条件(3.4)下的取得最小值时对应的函数，也就是

(A)找到 $\left({\mu }_g,{\Phi }_g=\left({\varphi }_{g,1},{\varphi }_{g,2}\right)\in U\right)$ 使得

$E_{{\beta }_d}\left({\Phi }_g\right)=\underset{\Phi \in U}{\min }{E}_{{\beta }_d}\left({\varphi }_1,{\varphi }_2\right)$ ，

${\mu }_g=E_{{\beta }_d}\left({\varphi }_{g,1},{\varphi }_{g,2}\right)+\frac{1}{2}{\displaystyle {\int }_{R^d}{\displaystyle \underset{j,l=1}{\overset{2}{\sum }}{\beta }_{d,jl}}}{\left|{\varphi }_{g,j}\left(x\right)\right|}^2{\left|{\varphi }_{g,l}\left(x\right)\right|}^2\text{d}x$ ，

这里的 $U$ 定义为：

$U=\left\{\Phi =\left({\varphi }_1,{\varphi }_2\right)|E_{{\beta }_d}\left({\varphi }_1,{\varphi }_2\right)<\infty ,{\displaystyle {\int }_{R^d}{\left|{\varphi }_j\left(x\right)\right|}^2\text{d}x}=1,j=1,2\right\}$ 。

我们可以证明 ${\Phi }_g\left(x\right)=\left({\varphi }_{g,1}\left(x\right),{\varphi }_{g,2}\left(x\right)\right)$ 满足：

$\mu {\varphi }_{g,1}\left(x\right)=\left(-\frac{1}{2}{\nabla }^2+V_{d,1}\left(x\right)+{\displaystyle \underset{l=1}{\overset{2}{\sum }}{\beta }_{d,1l}{\left|{\varphi }_{g,l}\right|}^2}-\Omega L_z\right){\varphi }_{g,1}\left(x\right)-\lambda {\varphi }_{g,2}(\; x\; )$

$\mu {\varphi }_{g,2}\left(x\right)=\left(-\frac{1}{2}{\nabla }^2+V_{d,2}\left(x\right)+{\displaystyle \underset{l=1}{\overset{2}{\sum }}{\beta }_{d,2l}{\left|{\varphi }_{g,l}\right|}^2}-\Omega L_z\right){\varphi }_{g,2}\left(x\right)-\lambda {\varphi }_{g,1}(\; x\; )$

且 ${\displaystyle \underset{j=1}{\overset{2}{\sum }}{\displaystyle {\int }_{R^d}{\left|{\varphi }_{g,j}\left(x\right)\right|}^2}}\text{d}x=1$ 。

很明显，静态相耦合的GP方程组(3.2)~(3.4)有许多解，基态解只是其中一类我们最关心的解。在物理学上，通常把GP方程组(3.2)~(3.4)除基态解以外的其它解称为激发态解，这是因为这些解所对应的能量比基态解所对应的能量 $E_{{\beta }_d}\left({\varphi }_{g,1},{\varphi }_{g,2}\right)$ 高。

3.2. 求基态解的数值方法

我们首先构建连续的梯度流

$\frac{\partial {\varphi }_1}{\partial t}=\frac{1}{2}{\nabla }^2{\varphi }_1-V_{d,1}{\varphi }_1-{\displaystyle \underset{l=1}{\overset{2}{\sum }}{\beta }_{d,1l}{\left|{\varphi }_l\right|}^2}{\varphi }_1+\Omega L_z{\varphi }_1+\lambda {\varphi }_2+{\mu }_{\Phi }{\varphi }_1$ (3.7)

$\frac{\partial {\varphi }_2}{\partial t}=\frac{1}{2}{\nabla }^2{\varphi }_2-V_{d,2}{\varphi }_2-{\displaystyle \underset{l=1}{\overset{2}{\sum }}{\beta }_{d,2l}{\left|{\varphi }_l\right|}^2}{\varphi }_2+\Omega L_z{\varphi }_2+\lambda {\varphi }_1+{\mu }_{\Phi }{\varphi }_2$ (3.8)

${\varphi }_j={\varphi }_j\left(x,t\right)=0,x\in \Gamma ,j=1,2$ (3.9)

${\varphi }_j\left(x,0\right)={\varphi }_{j,0}\left(x\right),x\in {\Omega }_x,{\displaystyle \underset{j=1}{\overset{2}{\sum }}{{\displaystyle {\int }_{{\Omega }_x}\left|{\varphi }_{j,0}\left(x\right)\right|}}^2\text{d}x}=1$ (3.10)

这里

${\mu }_{\Phi }=E_{{\beta }_d}\left({\varphi }_1,{\varphi }_2\right)+\frac{1}{2}{\displaystyle {\int }_{R^d}{\displaystyle \underset{j,l=1}{\overset{2}{\sum }}{\beta }_{d,jl}}}{\left|{\varphi }_j\left(x\right)\right|}^2{\left|{\varphi }_l\left(x\right)\right|}^{2}dx$ (3.11)

并且 ${\Omega }_x$ 是 $R^d$ 中的一个有界区域。 $\Gamma =\partial {\Omega }_x$ 。

对于连续的梯度流(3.7)~(3.10)，我们有下面的定理，

定理3.1：在临界条件(3.9)和初始条件(3.10)下，连续的梯度流(3.7)~(3.8)满足

(1) ${\displaystyle \underset{j=1}{\overset{2}{\sum }}{\displaystyle {\int }_{{\Omega }_x}{\left|{\varphi }_j\left(x,t\right)\right|}^2\text{d}x}}={\displaystyle \underset{j=1}{\overset{2}{\sum }}{\displaystyle {\int }_{{\Omega }_x}{\left|{\varphi }_j\left(x,0\right)\right|}^2\text{d}x}}=1$ (3.12)

(2) $\frac{\partial }{\partial t}{E}_{{\beta }_d}\left({\varphi }_1\left(x,t\right),{\varphi }_2\left(x,t\right)\right)\le 0$ (3.13)

这个定理的证明虽然很长但证法很简单，这里我们将省略，不再陈述。

从定理3.1我们可以知道连续的梯度流(3.7)~(3.10)可保持总模量守恒和能量递减的性质。并且，当 $t\to \infty$ 时，我们可得到 ${\varphi }_1$ 和 ${\varphi }_2$ 接近于稳定状态时的解，这时的能量函数 $E_{{\beta }_d}\left({\varphi }_1\left(x,t\right),{\varphi }_2\left(x,t\right)\right)$ 在U空间上处于能量最低的状态。如果我们能为方程(3.9)引入一个合适的初始条件 $\left({\varphi }_{1,0}\left(x\right),{\varphi }_{2,0}\left(x\right)\right)$ ，就可以找到方程组(3.7)~(3.10)的稳态解，也就找到了问题(A)的基态解：

${\varphi }_j\left(x,t\right)\to {\varphi }_{g,j}\left(x\right),t\to \infty ,\left(j=1,2\right)$ 。

其次，从时间 $t_n=n\Delta t$ 到 $t_{n+1}=\left(n+1\right)\Delta t$ ，我们希望求得方程组(3.7)~(3.10)的数值解，为了避免每步计算化学势能，我们取而代之考虑下面的离散梯度流：

$\frac{\partial {\varphi }_1}{\partial t}=\left[\frac{1}{2}{\nabla }^2-V_{d,1}-{\displaystyle \underset{j=1}{\overset{2}{\sum }}{\beta }_{d,1l}{\left|{\varphi }_l\right|}^2+\Omega L_z}\right]{\varphi }_1\left(x,t\right)+\lambda {\varphi }_2,t_n\le t\le t_{n+1}$ (3.14)

$\frac{\partial {\varphi }_{21}}{\partial t}=\left[\frac{1}{2}{\nabla }^2-V_{d,2}-{\displaystyle \underset{j=1}{\overset{2}{\sum }}{\beta }_{d,2l}{\left|{\varphi }_l\right|}^2+\Omega L_z}\right]{\varphi }_2\left(x,t\right)+\lambda {\varphi }_1,t_n\le t\le t_{n+1}$ (3.15)

${\varphi }_j\left(x,t_{n+1}\right)={\varphi }_j\left(x,t_{n+1}^{+}\right)=\frac{{\varphi }_j\left(x,t_{n+1}^{-}\right)}{\sqrt{{\Vert {\varphi }_1\left(x,t_{n+1}^{-}\right)\Vert }^2+{\Vert {\varphi }_2\left(x,t_{n+1}^{-}\right)\Vert }^2}}$ (3.16)

${\varphi }_j\left(x,t_{n+1}\right)=0,x\in \Gamma ,j=1,2$ (3.17)

${\varphi }_j\left(x,0\right)={\varphi }_{j,0}\left(x\right),x\in {\Omega }_x$ (3.18)

从离散梯度流(3.14)~(3.16)，我们发现他们在能量减少时的某种条件下的性质。

我们令

${\tilde{\varphi }}_j(x,t)=\frac{{\varphi }_j\left(x,t\right)}{{\Vert {\varphi }_j\left(x,t\right)\Vert }^2},t_n\le t\le t_{n+1},n\ge 0,j=1,2$ (3.19)

代入到梯度流(3.14)~(3.16)，很容易建立以下基本事实：

引理3.1：假设 $V_{d,j}\ge 0,{\beta }_{d,ij}\ge 0\left(i,j=1,2\right)$ ，且 ${\Vert {\varphi }_{j,0}\Vert }^2=1$ ，那么

(i) 对于所有的 $t^{\prime }$ 和 $t$ 有 $t_n\le t^{\prime }\le t\le t_{n+1}$ ，我们有：

$E_{{\beta }_d}\left({\varphi }_1\left(x,t\right),{\varphi }_2\left(x,t\right)\right)\le E_{{\beta }_d}\left({\varphi }_1\left(x,t^{\prime }\right),{\varphi }_2\left(x,t^{\prime }\right)\right)$ 。

(ii) 如果 ${\beta }_{d,ij}=0\left(i,j=1,2\right)$ ，那么

$E_{{\beta }_d}\left({\tilde{\varphi }}_1\left(x,t\right),{\tilde{\varphi }}_2\left(x,t\right)\right)\le E_{{\beta }_d}\left({\tilde{\varphi }}_1\left(x,t_n\right),{\tilde{\varphi }}_2\left(x,t_n\right)\right)$ (3.20)

证明：

(i)

$\frac{\partial }{\partial t}{E}_{{\beta }_d}\left({\varphi }_1\left(x,t\right),{\varphi }_2\left(x,t\right)\right)=-2{\displaystyle \underset{R^d}{\int }{\displaystyle \underset{j=1,2}{\sum }{\left|\frac{\partial }{\partial t}{\varphi }_j\left(x,t\right)\right|}^2}}\text{d}x\le 0$

这就证明了(i)中的结果。

(ii) 这个证明可以仿照 [24] 中的相关证明。这里，具体证明省略

从引理3.1，我们可直接得到：

定理3.2：设 ${\displaystyle \underset{j=1}{\overset{2}{\sum }}{\Vert {\varphi }_j\Vert }^2}=1$ ，如果 ${\beta }_{d,ij}=0\left(i,j=1,2\right)$ ，梯度流(3.14)~(3.16)满足

(1) ${\displaystyle {\sum }_{j=1}^2{\Vert {\varphi }_j\Vert }^2}=1$ ；

(2)

$\begin{array}{l}{E}_{{\beta }_d}\left({\varphi }_1\left(x,t_{n+1}\right),{\varphi }_2\left(x,t_{n+1}\right)\right)\le E_{{\beta }_d}\left({\varphi }_1\left(x,t_n\right),{\varphi }_2\left(x,t_n\right)\right)\\ \le \cdots \le E_{{\beta }_d}\left({\varphi }_1\left(x,0\right),{\varphi }_2\left(x,0\right)\right)=E_{{\beta }_d}\left({\varphi }_{1,0},{\varphi }_{2,0}\right),n\ge 0\end{array}$ (3.21)

事实上，方程中的标准化步骤(3.16)等同于准确地求解下面普通的常微分方程组

$\frac{\partial {\varphi }_j}{\partial t}\left(x,t\right)={\mu }_{\Phi }\left(t,k\right){\varphi }_j\left(x,t\right),x\in R^d,t_n<t<t_{n+1},n\ge 0$ 。 (3.22)

${\varphi }_j\left(x,t_n^{+}\right)={\varphi }_j\left(x,t_{n+1}^{-}\right),x\in R^d,j=1,2$ (3.23)

这里 $k=\Delta t_n$ 并且

${\mu }_{\Phi }\left(t,k\right)\equiv {\mu }_{\Phi }\left(t_{n+1},\Delta t_n\right)=-\frac{1}{2\Delta t_n}\ln \left({\displaystyle {\sum }_{j=1}^2{\Vert {\varphi }_j\left(x,t_{n+1}^{-}\right)\Vert }^2}\right),t_n<t<t_{n+1}$ (3.24)

因此，梯度流(3.14)~(3.16)可以看作作为下面离散带有间断系数的梯度流的一种一阶分裂法

$\frac{\partial {\varphi }_1}{\partial t}=\frac{1}{2}{\nabla }^2{\varphi }_1-V_{d,1}\left(x\right){\varphi }_1-{\displaystyle \underset{l=1}{\overset{2}{\sum }}{\beta }_{d,1l}}{\left|{\varphi }_l\right|}^2{\varphi }_1+\Omega L_z{\varphi }_1+\lambda {\varphi }_2+{\mu }_{\Phi }\left(t,k\right){\varphi }_1$ (3.25)

$\frac{\partial {\varphi }_2}{\partial t}=\frac{1}{2}{\nabla }^2{\varphi }_2-V_{d,2}\left(x\right){\varphi }_1-{\displaystyle \underset{l=1}{\overset{2}{\sum }}{\beta }_{d,2l}}{\left|{\varphi }_l\right|}^2{\varphi }_2+\Omega L_z{\varphi }_2+\lambda {\varphi }_1+{\mu }_{\Phi }\left(t,k\right){\varphi }_2$ (3.26)

对于上面的梯度流(3.25)，让 $k\to 0$ 并且注意到 ${\varphi }_j\left(x,t_{n+1}^{-}\right)$ (方程(3.23)的右边)在 $t_{n+1}=t_n+\Delta t_n$ 时是方程(3.23)的求解，我们得到

$\begin{array}{c}\underset{k\to 0^{+}}{\lim }{\mu }_{\Phi }\left(t,k\right)=\underset{\Delta t_n\to 0^{+}}{\lim }\frac{1}{-2\Delta t_n}\ln \left({\displaystyle \underset{i=1}{\overset{2}{\sum }}{\Vert {\varphi }_j\left(x,t_{n+1}^{-}\right)\Vert }^2}\right)=\underset{\Delta t_n\to 0^{+}}{\lim }=\frac{{\frac{\text{d}}{\text{d}\tau }\left({\displaystyle \underset{j=1}{\overset{2}{\sum }}{\Vert {\varphi }_j\left(x,t+\tau \right)\Vert }^2}\right)|}_{\tau =\Delta t_n}}{-2{\displaystyle \underset{j=1}{\overset{2}{\sum }}{\Vert {\varphi }_j\left(x,t+\Delta t_n\right)\Vert }^2}}\\ =\underset{\Delta t_n\to 0^{+}}{\lim }=\frac{{\mu }_{\Phi }^j\left({\varphi }_1\left(x,t+\Delta t_n\right),{\varphi }_2\left(x,t+\Delta t_n\right)\right)}{{\Vert {\varphi }_j\left(x,t+\Delta t_n\right)\Vert }^2}=\frac{\mu \left({\varphi }_1,{\varphi }_2\right)}{{\displaystyle \underset{j=1}{\overset{2}{\sum }}{\Vert {\varphi }_j\left(x,t\right)\Vert }^2}}\end{array}$ (3.27)

这里

$\mu \left({\varphi }_1,{\varphi }_2\right)=E_{{\beta }_d}\left({\varphi }_1,{\varphi }_2\right)+\frac{1}{2}{\displaystyle {\int }_{R^d}{\displaystyle \underset{j,l=1}{\overset{2}{\sum }}{\beta }_{d,jl}{\left|{\varphi }_j\left(x\right)\right|}^2{\left|{\varphi }_l\left(x\right)\right|}^2}\text{d}x}$ (3.28)

在下一小节里，我们给出离散梯度流(3.14)~(3.16)的一种差分格式。

离散梯度流(3.14)~(3.16)的一种差分格式

我们这里仅仅只考虑二维情形，三维情形的差分格式是类似的。假设我们计算的区域为 $\left[a,b\right]\times \left[c,d\right]$ 。我们选择时间步长 $k>0$ 以及空间步长 $h_x>0$ 和 $h_y>0$ 。并且 $h_x=\left(b-a\right)/M$ 和 $h_y=\left(d-c\right)/N$ (这里M和N均是偶的正整数)。空间上以及时间上的网格点分别为 $x_p,y_q$ 以及 $t_n$ ：

$x_p:=a+p{h}_x,p=0,1,\cdots ,M,y_q=c+q{h}_y,q=0,1,\cdots ,N,t_n=nk,n=0,1,\cdots$

设 ${\varphi }_{j,pq}^n$ 表示函数 ${\varphi }_j\left(x,y,t\right)$ 在点 $\left(x_p,y_q,t_n\right)$ 处的近似值( $j=1,2$ )。利用有限差分法，我们可得到下面的一种离散梯度流(3.14)~(3.16)的差分格式：

$\begin{array}{c}\frac{{\tilde{\varphi }}_{1,pq}-{\varphi }_{1,pq}^n}{k}=\frac{1}{2h_x{}^2}\left[{\tilde{\varphi }}_{1,\left(p+1\right)q}-2{\tilde{\varphi }}_{1,pq}+{\tilde{\varphi }}_{1,\left(p-1\right)q}\right]+\frac{1}{2h_y^2}\left[{\tilde{\varphi }}_{1,p\left(q+1\right)}-2{\tilde{\varphi }}_{1,pq}+{\tilde{\varphi }}_{1,p\left(q-1\right)}\right]\\ -V_{2,1}\left(x_p,y_q\right){\tilde{\varphi }}_{1,pq}-\left({\beta }_{2,11}{\left|{\varphi }_{1,pq}^n\right|}^2+{\beta }_{2,12}{\left|{\varphi }_{2,pq}^n\right|}^2\right){\tilde{\varphi }}_{2,pq}\\ +i\Omega y_q\frac{{\tilde{\varphi }}_{1,\left(p+1\right)q}-{\tilde{\varphi }}_{1,\left(p-1\right)q}}{2h_x}-i\Omega x_p\frac{{\tilde{\varphi }}_{1,p\left(q+1\right)}-{\tilde{\varphi }}_{1,p\left(q-1\right)}}{2h_y}+\lambda {\varphi }_{2,pq}^n\end{array}$ (3.29)

$\begin{array}{c}\frac{{\tilde{\varphi }}_{2,pq}-{\varphi }_{2,pq}^n}{k}=\frac{1}{2h_x{}^2}\left[{\tilde{\varphi }}_{2,\left(p+1\right)q}-2{\tilde{\varphi }}_{2,pq}+{\tilde{\varphi }}_{2,\left(p-1\right)q}\right]+\frac{1}{2h_y^2}\left[{\tilde{\varphi }}_{2,p\left(q+1\right)}-2{\tilde{\varphi }}_{2,pq}+{\tilde{\varphi }}_{2,p\left(q-1\right)}\right]\\ -V_{2,2}\left(x_p,y_q\right){\tilde{\varphi }}_{2,pq}-\left({\beta }_{2,21}{\left|\lambda {\varphi }_{1,pq}^n\right|}^2+{\beta }_{2,22}{\left|\lambda {\varphi }_{2,pq}^n\right|}^2\right){\tilde{\varphi }}_{2,pq}\\ +i\Omega y_q\frac{{\tilde{\varphi }}_{2,\left(p+1\right)q}-{\tilde{\varphi }}_{2,\left(p-1\right)q}}{2h_x}-i\Omega x_p\frac{{\tilde{\varphi }}_{2,p\left(q+1\right)}-{\tilde{\varphi }}_{2,p\left(q-1\right)}}{2h_y}+\lambda {\varphi }_{1,pq}^n\end{array}$ (3.30)

上述方程组(3.29)~(3.30)对所有的 $p=1,\cdots ,M-1,q=1,\cdots ,N-1$ ，都成立。

并且

${\tilde{\varphi }}_{j,0q}={\tilde{\varphi }}_{j,Mq}={\tilde{\varphi }}_{j,p0}={\tilde{\varphi }}_{j,pN}=0,p=0,\cdots ,M,q=0,\cdots ,N$ ，

${\varphi }_{j,pq}^{n+1}=\frac{{\tilde{\varphi }}_{j,pq}}{\sqrt{{\Vert {\tilde{\varphi }}_1\Vert }^2+{\Vert {\tilde{\varphi }}_2\Vert }^2}},p=0,1,\cdots ,M;q=0,1,\cdots ,N,n=0,1,\cdots$ ，

${\varphi }_{j,pq}^0={\varphi }_{j,0}\left(x_p,y_q\right),p=0,1,\cdots ,M;q=0,\cdots ,N$ ，

这里的 ${\Vert {\tilde{\varphi }}_j\Vert }^2\left(j=1,2\right)$ 的定义为 ${\Vert {\tilde{\varphi }}_j\Vert }^2=h_x{h}_y{\displaystyle \underset{p=1}{\overset{M-1}{\sum }}{\displaystyle \underset{q=1}{\overset{N-1}{\sum }}{\left|{\varphi }_{j,pq}\right|}^2}}$ 。

4. 数值实验结果

在本节的数值计算中，我们应用差分格式(3.29)~(3.30)来计算二维情形的基态解。因为当前成功的双原子–玻色爱因斯坦凝聚体实验主要考虑⁸⁷Rb及它的同位素 [25] 和²³Na及它的同位素 [3] ，所以在我们的数值计算中，我们仅仅考虑

$m_1=m_2,V_{2,1}\left(x,y\right)=V_{2,2}\left(x,y\right)=\frac{1}{2}\left(x^2+y^2\right)$ 。

关于相互作用的参数，我们也仅仅主要考虑下面三种情形

情形I： $\left(\begin{array}{l}{\beta }_{2,11}{\beta }_{2,11}\\ {\beta }_{2,21}{\beta }_{2,22}\end{array}\right)={\beta }_2\times \left(\begin{array}{l}1.0\delta \\ \delta 1.0\end{array}\right)\left(\delta >0\right)$ ，

情形II： $\left(\begin{array}{l}{\beta }_{2,11}{\beta }_{2,11}\\ {\beta }_{2,21}{\beta }_{2,22}\end{array}\right)={\beta }_2\times \left(\begin{array}{l}1.031.0\\ 1.00.97\end{array}\right)$ ，

情形III： $\left(\begin{array}{l}{\beta }_{2,11}{\beta }_{2,11}\\ {\beta }_{2,21}{\beta }_{2,22}\end{array}\right)={\beta }_2\times \left(\begin{array}{l}1.00.94\\ 0.940.97\end{array}\right)$ 。

其中情形II与情形III分别对应相应实验中的参数。在情形I中，两种原子内部以及两种原子之间均有相同的相互作用力；在情形II中，两种原子之间有相同的相互作用力，但是两种原子内部没有相同的相互作用力；在情形III中，两种原子之间有相同的相互作用力，但是两种原子内部没有相同的相互作用力。

Figure 1. Case I, ${\beta }_2=1000$ , $\delta =0.7$ , $\lambda =0.2$ , $\Omega =0.7$

图1. 情形I， ${\beta }_2=1000$ ， $\delta =0.7$ ， $\lambda =0.2$ ， $\Omega =0.7$

在图1中，N与E分别代表总模量与能量，它们是总模量(2.10)与能量(2.11)的数值近似。从图1中的结果，我们可以看出：我们的数值方法差分格式(3.29)~(3.30)在数值上能保证总模量守恒，且能量递减。

从图2中的结果，我们可以看出：(1) 在双原子–玻色爱因斯坦凝聚体的基态中，当两原子都含有数量较多的涡旋时且 $\delta <1$ ，它们分别形成涡旋网格；(2) 在双原子–玻色爱因斯坦凝聚体的基态中，一种原子的涡旋与另一种原子的涡旋镶嵌在一起；(3) 在双原子–玻色爱因斯坦凝聚体的基态中，随着 $\Omega$ 的增大，两种原子所含的涡旋数量也增大。

从图3中的结果，我们可以看出：(1) 在双原子–玻色爱因斯坦凝聚体的基态中，当两原子都含有数量较多的涡旋时且 $\delta >1$ ，它们分别形成带状的涡旋；(2) 在双原子–玻色爱因斯坦凝聚体的基态中，一种原子的涡旋与另一种原子的涡旋镶嵌在一起；(3) 在双原子–玻色爱因斯坦凝聚体的基态中，随着 $\Omega$ 的增大，两种原子所含的涡旋数量也增大。

从图4中的结果，我们可以看出：(1) 在双原子–玻色爱因斯坦凝聚体的基态中，两原子都含有数量较多的涡旋，一种原子的涡旋与另一种原子的涡旋镶嵌在一起；(2) 在双原子–玻色爱因斯坦凝聚体的基态中，随着 $\Omega$ 的增大，两种原子所含的涡旋数量也增大。

从图5中的结果，我们可以也看出：(1) 在双原子–玻色爱因斯坦凝聚体的基态中，两原子都含有数量较多的涡旋，一种原子的涡旋与另一种原子的涡旋镶嵌在一起。

在图6中，我们有较大的Josephson耦合力作用。从图6中的结果，我们可以也看出：(1) 在双原子–玻色爱因斯坦凝聚体的基态中，两原子都含有数量较多的涡旋；(2) 在双原子–玻色爱因斯坦凝聚体的基态中，一种原子的涡旋不是另一种原子的涡旋镶嵌在一起，而是它们都含有相同数量及同样形态的涡旋；(3) 在双原子–玻色爱因斯坦凝聚体的基态中，随着 $\Omega$ 的增大，两种原子所含的涡旋数量也增大。

Figure 2. Case I, ${\beta }_2=1000$ , $\delta =0.7$ , $\lambda =0.2$

图2. 情形I， ${\beta }_2=1000$ ， $\delta =0.7$ ， $\lambda =0.2$

Figure 3. Case I, ${\beta }_2=1000$ , $\delta =1.1$ , $\lambda =0.2$

图3. 情形I， ${\beta }_2=1000$ ， $\delta =1.1$ ， $\lambda =0.2$

Figure 4. Case II, ${\beta }_2=1000$ , $\lambda =0.2$

图4. 情形II， ${\beta }_2=1000$ ， $\lambda =0.2$

Figure 5. Case III, ${\beta }_2=1000$ , $\lambda =0.2$

图5. 情形III， ${\beta }_2=1000$ ， $\lambda =0.2$

Figure 6. Case III, ${\beta }_2=1000$ , $\lambda =2$

图6. 情形III， ${\beta }_2=1000$ ， $\lambda =2$

5. 结论

我们提出了一种梯度法来求解静态相耦合的GP方程组，并且应用该方法去研究旋转中的双原子–玻色爱因斯坦凝聚体在极低温度下所呈现的基态现象。我们提出的梯度法在数值上既保持总模量守恒又能使总能量递减。我们在数学上我们严格地证明我们提出的梯度法是一种获得能量函数在给定限制性条件下的最小值的十分有效的方法。在基态现象中，如果在双原子–玻色爱因斯坦凝聚体添加了Josephson耦合力后，我们发现复杂的涡旋结构。将来我们要做的工作是将我们新提出的梯度法应用到三维情形，因为它更有利于模拟双原子–玻色爱因斯坦凝聚体在极低温度下所呈现的基态现象，特别是涡旋现象的结构。

基金项目

国家自然科学基金项目(91430103)的资助。