非局域非线性相互作用下玻色–爱因斯坦凝聚体涡旋孤子动力学演化模式仿真研究

Simulation Study of the Dynamic Evolution Model of Vortex Solitons of Bose-Einstein Condensates with Non-Local Nonlinear Interactions

摘要: 基于二维含非局域非线性相互作用的Gross-Pitaevskii方程(GPE)，我们研究了二维玻色爱因斯坦凝聚体中在非局域的非线性相互作用下体系涡旋孤子的演化问题，我们的工作考虑了非局域非线性相互作用，也就是在传统的三次非线性相互作用之外的高阶非线性作用对涡旋孤子的演化影响，通过数值仿真的方法我们研究了体系涡旋孤子随时间演化的图样以及演化的一些非局域非线性作用的效应，得到了涡旋孤子可能的周期性演化与单调衰减的演化模式，相应地，数值计算与仿真结果与前期实验理论计算得出的结果吻合，我们的工作可以用以指导二维玻色爱因斯坦凝聚体中考虑进非局域非线性作用影响下体系涡旋孤子的演化的实验观测。

Abstract: Based on the two-dimensional Gross-Pitaevskii Equation with nonlocal nonlinear interactions, we study the evolution of vortex solitons in two-dimensional Bose-Einstein condensates under non-local nonlinear interactions. Our work takes into account the nonlocal nonlinear interactions, in other words, the evolution of vortex solitons influenced by higher-order nonlinear actions other than the traditional cubic nonlinear interactions. Through numerical simulation, we studied the evolution pattern of vortex solitons of the system over time and some effects of the evolution caused by nonlocal nonlinear actions, and obtained the evolution model of periodic and monotonic decay evolution of vortex solitons. The numerical and simulation results are in agreement with the previous experimental theoretical results, and our work can be used to guide the experimental observation of the evolution of vortex solitons in two-dimensional Bose-Einstein condensates under the influence of local nonlinear action.

1. 引言

在过去相当长一段时间内玻色爱因斯坦凝聚体相关的动力学演化问题备受瞩目，在实验和理论上都有较多研究，相关理论研究中，GPE方程被证明是可靠的理论模型，用于描述玻色爱因斯坦凝聚体中的孤子以及涡旋的演化，孤子、涡旋的动力学演化本身在实验展开理论上都备受瞩目，其中比较重要的现象是二维玻色爱因斯坦凝聚体中的涡旋演化问题，其中相应前期工作对它基本的只含一阶非线性项的GPE方程理论以及实验研究有相应的报导，但是，实际的一些场景中相应的非线性GPE方程需要考虑高阶以及非局域性的非线性影响，这个可以很好的解释非局域非线性项以及高阶非线性项。解释涡旋演化中的一些奇特的现象，例如，它的寿命以及调制不稳定性[1]-[5]。

先前关于玻色爱因斯坦凝聚体中的涡旋演化大多未考虑仅非局域相互作用[2]。在本项工作中，基于二维含非局域非线性相互作用的GPE方程[6]-[8]，通过数值计算与数值模拟的方法计算研究其中在一定初始条件下涡旋孤子的演化问题，本文采用了含非局域非线性相互作用的二维GPE方程，通过数值模拟与仿真研究在初始涡旋度为1的单涡旋孤子初始条件下体系随时间的演化，本文的数值模拟与仿真表明，在一定条件下，涡旋孤子会出现展宽与半径会出现周期性振荡，这是涡旋孤子的呼吸振荡模式，与前期实验与理论报导中的涡旋孤子振动模式较吻合，而且在不同的参数设置下还会有单调衰减模式，说明了本项工作的数值模拟相应的研究结果比较好地反映前期工作中报导的一些涡旋孤子的演化现象[9]-[12]。

本文的工作安排如下，下一部分介绍问题的理论模型，包括二维含非局域非线性相互作用的GPE方程模型(二维GPE模型)，介绍方程各项参数的意义，接下来一部分通过介绍了对二维GPE方程描述的二维玻色爱因斯坦凝聚体系的涡旋孤子的演化用数值的方法进行仿真模拟研究得出了涡旋孤子的周期性振荡与单调衰减模式，我们的研究结果可用于指导二维含非局域非线性相互作用的玻色爱因斯坦凝聚体系的涡旋孤子的实验观测。

2. 二维玻色爱因斯坦凝聚体与含局域非线性相互作用效应GPE模型

二维玻色爱因斯坦凝聚体在各向同性谐振子外势并在非局域非线性相互作用下的GPE理论模型为：

$ih\frac{\partial }{\partial t}\psi \left(x,y,t\right)=-\frac{{\hslash }^2}{2m}{\nabla }^2\psi \left(x,y,t\right)+\frac{1}{2}m{\omega }^2\left(x^2+y^2\right)\psi \left(x,y,t\right)+R\left(x,y,t\right)\psi \left(x,y,t\right)$ (1)

其中$R\left(x,y,t\right)$ 描述非局域相互作用的关联函数(通常p函数取值$\propto \exp \left(-\frac{{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2}{2{\sigma }^2}\right)$ )，m为组成玻

色爱因斯坦凝聚体原子的质量，$\omega$ 为谐振子陷俘势强度。

3. 二维GPE数值仿真求解

这主要基于前一部分中二维涡旋孤子的GPE方程的一些数值仿真，数值求解方程(1)，我们要求的是波函数$\psi$ 的数值解，它有三个变量：x、y、t。其中，非线性项含非局域项R

$R\left(x,y,t\right)={\displaystyle \int \text{d}x}{\displaystyle \int \text{d}y}{\psi }^2\left(x^{\prime },y^{\prime },t\right)p\left(x,y,x^{\prime },y^{\prime }\right){\psi }^2\left(x,y,t\right)$

$p\left(x,y,x^{\prime },y^{\prime }\right)$ 为相互作用函数，$\left(x^{\prime },y^{\prime }\right)$ 不在$\left(x,y\right)$ 位置处，它影响了$\left(x,y\right)$ 的相互作用程度。这个函数的R是非局域项，一般来说，不方便用一些目前通用的数值计算软件来求解，例如MATLAB或是mathematica。对此我们用数值求解，首先设计算法，试着对这个系统演化进行仿真，程序运算结果为已经给出的波函数，即不同时刻，不同位置的波函数值。波函数初始值为

${\psi }_0\left(x,y,t\right)=C_0\left(x+iy\right){\text{e}}^{-\frac{x^2+y^2}{{\sigma }^2}}$

波函数取的初值为一个涡旋态，通过一个涡旋操控可以产生一个涡旋，${\psi }_0$ 涡旋态的形式，其中$C_0$ 为常数，初始态方程为(1) (2)，以下为求解的主要思路。

首先为了编程处理，我们假设一个波函数，将求得的波函数的数据$\psi \left(x_i,y_j,t_k\right)$ ，放在一个阵列中，在离散的点$\left(x_i,y_j,t_k\right)$ 的位置上，通过拟合可以得到在波函数在任意的空间位置，以及时间区间的任意一个时间点上的波函数值。求解空间上的阵列，我们取一个正方形，其中任意一个离散的点为$\left(x_i,y_j\right)$ ，然后取一个时间区间t，$t\in \left(0,T_0\right)$ 。$\left(x_i,y_j\right)$ 可以看成离散点的点数，例如取值为600 × 600、800 × 800或者900 × 900等，同时t也可以取一些离散的点。在这种情况下，存在$\left(x_i,y_j\right)$ 这些离散的点为均匀分布的特点，包括时间点的取值也是均匀的。然后波函数的波形图，在不同的位置处，它的变化率是不同的。例如基态，在x趋于无穷大或在0附近时，它的变化比较平缓，在这中间的某个区域变化幅度就比较大。对于这种情况，我们关心的主要是二阶导数，而非一阶导数。在数值求解时，一般来说，取得点数越多，精度越高，但计算量太大。因此为了保证精度，算法的设计显得尤为重要。关于波函数的求解，首先波函数$\psi$ 是一个复函数，它有实部和虚部，即

$\psi \left(x,y,t\right)=u\left(x,y,t\right)+iv\left(x,y,t\right)$

实际上我们在求解这个方程时，因为计算机处理也只是对实函数的一些操作，所以我们要对u和v分别求解。在$t_p$ 时刻，$\psi$ 的值为$\left\{\psi \left(x_i,y_j,t_q\right)\right\}$ ，其中$q\in \left(p-1,p-2,\cdots ,p-s\right)$ ，s为任意数，常取$s=3$ ，即在某一离散时刻$t_p$ 的值有$t_{p-1}$ 、$t_{p-2}$ 和$t_{p-3}$ 。这种拟合为最小二乘法拟合。在$t_p$ 时刻，根据方程(1)，我们可以看到函数含有关于$\psi$ 函数空间坐标的二阶导数，所以我们可以给出$t_p$ 时刻波函数的实部和虚部，即$u\left(x_i,y_j,t_p\right)$ 、$v\left(x_i,y_j,t_p\right)$ 。对方程(1)，从$t_{p-1}$ 到$t_p$ 进行积分。积分后先看实部的话，

$\begin{array}{l}u\left(x_i,y_j,t_p\right)-u\left(x_i,y_j,t_{p-1}\right)\\ ={\displaystyle {\int }_{t_{p-1}}^{t_p}-\frac{{\hslash }^2}{2m}\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right)u\left(x_i,y_j,t_p\right)\text{d}t}+{\displaystyle {\int }_{t_{p-1}}^{t_p}\frac{1}{2}m{\omega }^2\left(x_i^2+y_j^2\right)\psi \left(x_i,y_j,t\right)\text{d}t}\\ +{\displaystyle {\int }_{t_{p-1}}^{t_p}R\left(x_i,y_j\right)\psi \left(x_i,y_j,t\right)\text{d}t}\end{array}$

从数值处理来说的话，我们可以通过拟合曲线得到u的二阶导数。对于离散点，它涉及到u对x和y的

二阶导函数。在t和$\frac{{\partial }^2}{\partial x^2}u$ 的拟合曲线中，取三个点，分别是$t_{p-1}$ 、$t_{p-2}$ 和$t_{p-3}$ ，考虑最小二乘法，用u关于

x的二阶导函数来拟合，结果是关于x的二次函数。如果是u关于y的二阶导函数，它的拟合曲线就是一个关于y的二次函数。关于导函数$u\left(x_i,y_j,t_p\right)$ 从$t_{p-1}$ 到$t_p$ 的积分，通过数据拟合可以得到含有三个参量及非线性项的二次函数，且包括谐振子势。

实部中涉及到三个数学积分，这三个数学积分的处理方法都是类似的。我们已知的一些数据点是离散的点，但积分是一个连续的区间，在连续区间内任意一点的函数值，可以通过拟合曲线，得到二次函

数。实际上我们可以这个积分变成一个二次函数的积分。例如$\frac{{\partial }^2}{\partial x^2}u\left(x_i,y_j,t\right)$ ，前面三个点就是$\frac{{\partial }^2}{\partial x^2}u\left(x_i,y_j,t_{p-3}\right)$ 、$\frac{{\partial }^2}{\partial x^2}u\left(x_i,y_j,t_{p-2}\right)$ 、$\frac{{\partial }^2}{\partial x^2}u\left(x_i,y_j,t_{p-1}\right)$ 。这三个点可以拟合出一个曲线。$\frac{{\partial }^2}{\partial x^2}u\left(x_i,y_j,t\right)$ 可以写成$a_x{t}^2+b_{x}t+c_x$ ，即

$\frac{{\partial }^2}{\partial x^2}u\left(x_i,y_j,t\right)=a_x{t}^2+b_{x}t+c_x$

然后通过拟合可以得到$a_x$ 、$b_x$ 和$c_x$ 。所以，

$\begin{array}{c}{\displaystyle {\int }_{t_{p-1}}^{t_p}\frac{{\partial }^2}{\partial x^2}u\left(x_i,y_j,t\right)\text{d}t}={\displaystyle {\int }_{t_{p-1}}^{t_p}\left(a_x{t}^2+b_{x}t+c_x\right)\text{d}t}\\ =\frac{a_x}{3}\left(t_p^3-t_{p-1}^3\right)+\frac{b_x}{2}\left(t_p^2-t_{p-1}^2\right)+c_x\left(t_p-t_{p-1}\right)\end{array}$

类似地，

$\frac{{\partial }^2}{\partial y^2}u\left(x_i,y_j,t\right)=a_y{t}^2+b_{y}t+c_y$

对于$\frac{1}{2}m{\omega }^2\left(x_i^2+y_j^2\right)\psi \left(x_i,y_j,t\right)$ ，它要拟合的数据点就不是导函数了，它是关于$\psi$ 函数在$\psi \left(x_i,y_j,t_{p-3}\right)$ 、$\psi \left(x_i,y_j,t_{p-2}\right)$ 、$\psi \left(x_i,y_j,t_{p-1}\right)$ 这三个数据点的值，它们对应的水平坐标分别是$t_{p-3}$ 、$t_{p-2}$ 和$t_{p-1}$ ，然后拟合得到

$\psi \left(x_i,y_j,t\right)=a_h{t}^2+b_{h}t+c_h$

因此，

$\begin{array}{l}{\displaystyle {\int }_{t_{p-1}}^{t_p}\frac{1}{2}m{\omega }^2\left(x_i^2+y_j^2\right)\psi \left(x_i,y_j,t\right)\text{d}t}={\displaystyle {\int }_{t_{p-1}}^{t_p}\frac{1}{2}m{\omega }^2\left(x_i^2+y_j^2\right)\left(a_h{t}^2+b_{h}t+c_h\right)\text{d}t}\\ =\frac{1}{2}m{\omega }^2\left(x_i^2+y_j^2\right)\left[\frac{a_h}{3}\left(t_p^3-t_{p-1}^3\right)+\frac{b_h}{2}\left(t_p^2-t_{p-1}^2\right)+c_h\left(t_p-t_{p-1}\right)\right]\end{array}$

关于第三个积分项，也是可以通过拟合已得的早些时间的数据点得到，通过函数拟合，然后可以积分出来。通过这个步骤，就可以把所有的数据点，也就是三个维度的自变量取任意值时的函数值求出来。需要注意的是，通过数值仿真求得的不同时刻的数据点是离散的。通过以上步骤就可以得到一个完整的关于体系演化的描述，即通过波函数$\psi$ 的数值仿真与模拟得到不同时刻的数值演化的情形。

其中关于$t=t_p$ 时波函数数值结果的更新，注意波函数的实部和虚部。u和v都是x、y和t的函数，起始条件为$t=t_1$ 、$t_2$ 时，

$\psi \left(x_i,y_j,t_{k=1、2}\right)={\psi }_0\left(x_i,y_j,0\right)$

由初始条件确定，$t=t_p$ 时刻的波函数可以通过对t的积分，算出来

$\psi \left(x_i,y_j,t_p\right)=u\left(x_i,y_j,t_p\right)+iv\left(x_i,y_j,t_p\right)$

后续可以进行迭代，例如$t=t_{p+1}$ 时刻的计算。我们需要计算在时间序列上，$t=t_p$ 时刻的相应的一些

量，包括$t=t_p$ 时刻的一些导函数，例如$\frac{{\partial }^{2}u}{\partial x^2}$ 、$\frac{{\partial }^{2}u}{\partial y^2}$ 、$\frac{{\partial }^{2}v}{\partial x^2}$ 和$\frac{{\partial }^{2}v}{\partial y^2}$ 。确定这四个量在$t=t_p$ 时刻的取值，也

是通过拟合。这里拟合可以分两种情形：一种是以$\left(x_i,y_j,t_p\right)$ 为中心，其中$2\le i\le N-2$ ，$2\le j\le N-2$ ，$i、j\in \left(0,N\right)$ 。然后插值进行拟合，在三个数据点$\left\{\left(x_{i-1},y_{j-1},t_p\right),\left(x_i,y_j,t_p\right),\left(x_{i+1},y_{j+1},t_p\right)\right\}$ 上拟合u和v，得到

$u\left(x,y,t_p\right)=a_{xx}^u{x}^2+b_{xx}^{u}x+c_{xy}^u+a_{yy}^u{y}^2+b_{yy}^{u}y$

其中，$x\in \left(x_{i-1},x_i\right)$ ，$y\in \left(y_{j-1},y_j\right)$ 。

这样我们可以得到

$\frac{{\partial }^{2}u}{\partial x^2}=2a_{xx}^u,\frac{{\partial }^{2}u}{\partial y^2}=2b_{yy}^u$

然后就可以进行下一步的迭代，计算$t=t_{p+1}$ 时刻的情况，重复整个流程，计算出$t=t_{p+1}$ 时刻相应的波函数，即对t进行从$t_p$ 到$t_{p+1}$ 的积分，进入下一个循环的流程。对于虚部的函数，也可以用同样的方法求出来，

$\frac{{\partial }^{2}v}{\partial x^2}=2a_{yy}^v,\frac{{\partial }^{2}v}{\partial y^2}=b_{yy}^v$

这就是数值计算的算法步骤。

数值计算结果如图1与图2所示。同样由于涡旋孤子在演化过程中处于圆对称状态，在演化过程中涡旋形状没有突变，相邻图之间极角维度方向没有明显变化。图1中的图1(b)和图1(c)一个周期内涡旋半径变大，到图1(d)又回到周期起始最小形状状态，为周期性振荡模式。图2仿真运行较长时间内显示出的单调衰减模式，相应的数值仿真结果可与相关理论计算结果图样吻合。

4. 结论

在本项工作中，我们研究了二维玻色爱因斯坦在凝聚体中，在非局域非线性相互作用下涡旋孤子演化的模式研究，基于二维含非局域非线性相互作用的GPE方程，我们用数值仿真与模拟的方法研究了在一定条件下体系涡旋孤子的演化模式，发现对于涡度为1的轴对称的主涡旋孤子在一定条件下可以存在周期性振荡模式与单调衰减模式，显示了非局域非线性相互作用的稳定调制效应。我们研究的结果与前期报道的实验或理论周期性演化模式主涡旋孤子演化模式比较吻合，我们的研究结果可用以指导二维玻色爱因斯坦凝聚体中在非局域非线性相互作用下涡旋孤子演化的实验研究。

Figure 1. Numerical simulation results show the vortex soliton patterns at four time points ((a) 0, (b) $T_0/4$ , (c) $T_0/2$ , (d) $3T_0/4$ ) within one cycle ($T_0=2\pi /\omega$ ) of the periodic oscillation mode (with both horizontal and vertical coordinates ${\sigma }_0/1.5$ as unit length)

图1. 数值仿真计算得到周期振荡模式一个周期($T_0=2\pi /\omega$ )内四个时间点((a) 0, (b) $T_0/4$ , (c) $T_0/2$ , (d) $3T_0/4$ )的涡旋孤子图样(横竖坐标均以${\sigma }_0/1.5$ 为单位长度)

Figure 2. Numerical simulation results show the vortex soliton patterns (with horizontal and vertical coordinates ${\sigma }_0/1.5$ as unit length) at four time points ((a) 0, (b) $T_0/2$ , (c) $T_0$ , (d) $3T_0/2$ ) in a monotonic decay mode ($2T_0$ time range)

图2. 数值仿真计算得到单调衰减模式($2T_0$ 时间范围)四个时间点((a) 0, (b) $T_0/2$ , (c) $T_0$ , (d) $3T_0/2$ )的涡旋孤子图样(横竖坐标均以${\sigma }_0/1.5$ 为单位长度)

基金项目

本项工作受国家自然科学基金11791240178与11547024资助。

NOTES

^*通讯作者。