摘要:具有三体相互作用的玻色-爱因斯坦凝聚体(Bose-Einstein Condensate, BEC)束缚于雅可比椭圆周期势中, 在平均场近似下可用3—5次Gross-Pitaevskii方程(GPE)描述. 首先利用多重尺度法对该系统进行了理论分析, 将GPE化为一定态非线性薛定谔方程(Nonlinear Schr?dinger Equation, NLSE), 并给出了一类带隙孤子的解析表达式. 然后采用牛顿共轭梯度法数值得到了该系统中存在的两类带隙孤子, 发现孤子的振幅随着三体相互作用的增强而减小, 这与多重尺度法分析所得结论一致. 最后用时间劈裂傅里叶谱方法对GPE进行长时间动力学演化以考察孤子的稳定性, 发现系统中既存在稳定的带隙孤子, 也存在不稳定的带隙孤子, 且外势的模数会对孤子的结构和稳定性产生明显影响.
关键词: 玻色-爱因斯坦凝聚/
带隙孤子/
三体相互作用/
稳定性

English Abstract

全文HTML

--> --> -->

玻色-爱因斯坦凝聚(Bose-Einstein condensate, BEC)现象是在宏观尺度上能观察到的最显著的多体量子现象之一, 最早是由Bose和Einstein在1924年提出, 即理想的玻色子在非常低的温度下, 大部分粒子会突然跌落到最低能级上, 处于这种新状态的物质被称为BEC^[1,2]. 随着在$ ^{87}{\rm Rb} $和$ ^{23}{\rm Na} $等一系列碱金属原子气体实验中实现了BEC, 对于具有弱相互作用的原子气体中的BEC物质波孤子现象引起了研究者的注意^[3]. 孤子也是自然界中普遍存在的一种非线性现象, 广泛存在于水波、粒子物理、等离子体、分子生物学及纤维等各种非线性介质中^[4,5]. 作为一种非线性波, 孤子因其独特的传播性质及潜在的应用价值, 成为非线性科学研究领域的重要研究课题之一. 随着BEC和简并费米气体的实验实现, 研究表明超冷原子气体中也存在物质波孤子现象. 实验中已经相继发现物质波亮孤子、暗孤子及涡旋孤子等非线性现象^[6–9].
一般在低浓度BEC中, 原子间相互作用距离尺度远小于原子间距离, 这时只需考虑两体相互作用, s-波散射很重要^[10,11]. 但浓度较高时, 例如BEC在原子芯片和原子波导表面的发展将会涉及到强压缩和密度的提高, 则需考虑三体相互作用的影响^[12–16]. 近年来在铯原子超冷气体实验中已发现三体相互作用现象^[17]. 三体作用对于玻色凝聚气体的稳定性有着重要作用. 例如, 它可以明显改变凝聚体的稳定区域, 即使强度很小也可使得凝聚原子数量增加. 另一方面, 三体作用除了使凝聚体密度分布发生变化外, 也会改变集体振荡激发的光谱, 因为这时系统的可压缩性受到三体相互作用对基态能量贡献的影响^[18]. 考虑三体相互作用对于研究BEC在光晶格中呼吸子的性质也有重要意义^[19]. 在平均场近似下, 超低温下稀薄BEC的动力学行为可用Gross-Pitaevskii方程(GPE)描述, 它可从海森堡方程导出^[20,21]. GPE是一非线性方程, 一般情况下很难得到其精确解析解, 不过在特殊情形下却能由许多方法给出它的精确孤子解, 例如Darboux变换法^[22]、 Backlund变换^[23], Hirota直接法^[24]、Painleve展开方法^[25,26]和变分法^[27]等.
具有周期调制的系统中会出现一系列新的效应, 特别是线性能带谱的带隙中会存在一种结构丰富的孤子, 称为带隙孤子^[28,29]. 该类孤子可存在于不同类型的非线性系统中, 包括低维光子晶体、光子层状结构^[30,31]和光晶格势中的BEC^[32–34]等. 本文考虑具有三体相互作用的准一维BEC束缚于雅可比椭圆周期势中(实验中该势可由两束不同频率的激光叠加来近似), 在平均场近似下系统的动力学可由3—5次GPE描述. 利用多重尺度法对系统的带隙孤子进行了理论分析, 将3—5次GPE化为一定态非线性薛定谔方程(nonlinear Schr?dinger equation, NLSE), 并给出了一类带隙孤子的解析表达式. 然后采用牛顿共轭梯度法数值得到了该系统中存在的带隙孤子, 包括两类基本带隙孤子(on-site孤子和off-site孤子)与亚基本带隙孤子(sub-fundamental gap soliton). 最后, 用时间劈裂傅里叶谱方法对GPE进行长时间动力学演化以考察孤子的稳定性, 发现on-site孤子始终稳定, 同相偶极孤子和异相偶极孤子既有稳定的也有不稳定的. 外势的模数对孤子的结构和稳定性会产生明显影响.

2

-->

2.1.理论模型

具有三体相互作用的BEC束缚于外势中, 在平均场近似下系统的动力学行为可用如下GPE描述^[35]

$\begin{split}{\rm i}\hbar\frac{\partial\varPsi( r,t)}{\partial t} =\, & - \frac{\hbar^2}{2m}{\nabla^2}\varPsi( r,t)+V( r)\varPsi( r,t)\\ & +g|\varPsi( r,t) |^2\varPsi( r,t) \\ & +\eta |\varPsi( r,t)|^4\varPsi( r,t), \end{split} $

其中$ \varPsi = \varPsi( r, t) $是系统的序参量(波函数), $ r = $ $ (x, y, z) $是位矢, $ g = {4{\text{π}}\hbar^2 a_{\rm s}}/{m} $表征原子间相互作用强度, m是原子质量. $ a_{\rm s} $是s波散射长度($ a_{\rm s} > 0 $表示原子间相互排斥, $ a_{\rm s} < 0 $表示原子间相互吸引), 可由Feshbach共振调节^[36]. η表征三体相互作用强度, $ V( r) $表示外势, 总粒子数$ N \!=\! \displaystyle\int\!\mid\varPsi( r, t)\mid^2 {\rm d}^3 r $. 实验中, BEC通常束缚于一谐振子势阱$V( r) = $$ \dfrac{1}{2}m(\omega_x^2+\omega_y^2+\omega_z^2) $内, 其中$ \omega _x $, $ \omega _y $, $ \omega _z $分别表示x, y, z方向的频率. 当$ \omega _x \approx \omega _y $且$ \omega _z \gg \omega _x $时, BEC在z方向被“冻结”于基态, 这时系统可用准二维GPE描述; 当$ \omega _y, \omega _z\gg \omega _x $时, BEC在y, z方向均被“冻结”, 将在x方向呈雪茄状分布, 系统的动力学行为可用准一维GPE描述^[37].
考虑BEC被囚禁于光晶格势下的准一维情形, 并做如下无量纲化处理: $ \tilde x = x/a_0 $, $ \tilde t = \omega _xt $, $ \varPsi( r,\, t) = \tilde \varPsi(x,\, t)\left( {\dfrac{{{m^2}{\omega _y}{\omega _z}}}{{{{\text{π}} ^2}{\hbar ^2}}}} \right)^{1/4}{\rm e}^{-m(\omega_yy^2\,+\,\omega_zz^2)/2\hbar}\,,$ $ \tilde g = \dfrac{n_0\sqrt{\omega_z\omega_y}g}{2{\text{π}}{\hbar}a_0^{3}\omega_x} = \dfrac{2 a_sn_0\sqrt{\omega_z\omega_y}}{a_0\omega_x} $, $ \tilde\eta = \dfrac{n_0^{2}\omega_z\omega_y}{3{\text{π}}^2 a_0^3\hbar\omega_x^3}{\eta} $, $ a_0 = \sqrt{{\hbar}/({m\omega_x})}$, 便得

$ {\rm i}\frac{\partial\varPsi}{\partial t} = -\frac{1}{2}\frac{\partial^2\varPsi}{\partial x^2}+ V(x)\varPsi+g|\varPsi|^2\varPsi +\eta|\varPsi|^4\varPsi, $

其中$ n_0 $为无量纲化时粒子数密度的特征量, 此时总粒子数为$ N = n_0\displaystyle\int|\varPsi(x, t)|^2 {\rm d}x $. 为方便起见, 这里已省略变量上面的“~”.
取外势$ V(x) = V_0{\rm sn}^2(x, q) $, 其中$ {\rm sn}(x, q) $为雅可比椭圆正弦函数, q是其模数($ 0\leqslant q\leqslant 1 $). 当$ q = 0 $时, $ V(x) = V_0\sin^2 x $, 文献[38]中考虑的外势即为此情形, 所以该势可看作是对相关文献的推广. 当$ q = 1 $时, $ V(x) = V_0\tanh^2 x $, 此时$ V(x) $并非周期势, 故后面不考虑此情形. 若考虑该势的实验实现, 可用下面公式^[39]

$ {\rm sn}(x,q) = \frac{2{\text{π}}}{qK(q)}\sum\limits_{n = 0}^\infty{\frac{\alpha^{n + 1/2}}{1 - \alpha^{2n + 1}}\sin \frac{(2n + 1){\text{π}} x}{2K(q)}}, $

其中$ \alpha = \exp\big\{-{\text{π}}[K(\sqrt{1-q^2})/K(q)]\big\} $, $ K(q) = $$\displaystyle\int_0^{\textstyle\frac{{\text{π}}}{2}} {\dfrac{{\rm d}\theta}{\sqrt{1-q^2\sin ^2\theta}}} $是第一类完全椭圆积分. $ V(x) $可近似为

$\begin{split}V(x)\approx\, &\frac{2V_0{\text{π}} ^2\alpha}{q^2K^2(q)(1-\alpha)^2}\Bigg[1 - \cos \frac{{\text{π}} x}{K(q)}\Bigg.\\& +\Bigg. \frac{2\alpha }{1 \!+\! \alpha \!+\! \alpha ^2}\left( \cos \frac{{\text{π}} x}{K(q)} \!-\! \cos \frac{2{\text{π}} x}{K(q)}\right)\Bigg], \end{split}$

从中可以看出, 该势在实验中可用两束不同频率的激光来近似实现. 当$ q < 0.9 $时, 其近似程度可达99%.
寻找方程(2)下列形式的定态解

$ \varPsi(x,t) = \psi(x){{\rm e}^{-{\rm i}\mu t}} , $

其中$ \psi(x) $为实函数, μ为化学势. 代入方程(2)后得

$ -\frac{1}{2}{\psi_{xx}}+V(x)\psi-\mu\psi+g\psi^3+\eta\psi^5 = 0 . $

若波函数振幅为无穷小, 则可忽略非线性项, 得

$ -\frac{1}{2}{\psi_{xx}}+V(x)\psi-\mu\psi = 0. $

方程(7)为一广义马丢(Mathieu)方程, 它的有界解被称为布洛赫模(Bloch Modes), 对应的化学势μ构成了布洛赫带(Bloch Bands). 一般情况下, 它的有界解可写为

$ \psi = p(x) = {{\rm e}^{{\rm i}kx}}\tilde{p}(x;\mu). $

将(8)式代入方程(7)后, 通过求解所得本征问题便可得色散关系$ \mu = \mu(k) $, 从而得到系统的能带结构. 这一点在文献[29]中有详细讨论, 此处不再赘述. 在Bands中, 线性波可以传播; 在Gaps中, 虽然线性波无法传播, 但可存在具有局域结构的非线性波, 即带隙孤子(gap soliton).
2 -->

2.2.带隙孤子的多重尺度法分析

若波函数振幅不是无穷小, 则方程(6)中非线性项不能忽略. 考虑化学势μ从Band的边界$ k = k_0 $, $ \mu_0 = \mu(k_0) $处进入带隙, 且$ k-k_0 $, $ \mu-\mu_0 $均为小量的情形, 可用多重尺度法对其进行理论分析. 引入多重尺度$ X_0 = x $, $ X_1 = {\varepsilon}x $, 将化学势与波函数展开为

$ \mu = \mu_0+\mu_2\varepsilon^2+\cdots , $

$\begin{split} \psi(x) =\, & \varepsilon\psi_0(X_0,X_1)+\varepsilon^2\psi_1(X_0,X_1)\\ &+\varepsilon^3\psi_2(X_0,X_1)+\cdots,\end{split}$

其中$ \varepsilon = k-k_0 $是一小量. 将此展开式代入方程(6), 并假设$ \eta = \hat{\eta}\varepsilon^{-2} = O(\varepsilon^ {-2}) $. 按ε的同次幂项合并后得

$ \varepsilon^1: L_0\psi _0 = 0 , $

$ \varepsilon^2: L_0\psi _1 = \frac{\partial ^2\psi _0}{\partial X_0\partial X_1} , $

$ \varepsilon^3: L_0\psi _2 = \frac{1}{2}\frac{\partial ^2\psi _0}{\partial X_1^2} + \frac{\partial ^2\psi _1}{\partial X_0\partial X_1} + \mu_2\psi _0 - g\psi _0^3 - \hat{\eta}\psi _0^5 , $

其中$ L_0 = -\dfrac{1}{2}\dfrac{\partial ^2}{\partial X_0^2} + V(X_0)-\mu_0 $.
方程(11)与(7)形式相同, 故具有下列形式的解

$ \psi_0(X_0,X_1) = B(X_1)p(X_0), $

其中$ B(X_1) $表示慢变包络(调制波), $ p(X_0) $为快变载波, 显然该解满足Fredholm条件^[40], 这里函数$ f_1(x) $与$ f_2(x) $的内积定义为

$ \langle{f_1, f_2}\rangle = \displaystyle\int_0^{4 K(q)}f_1(x)f_2^*(x){\rm d}x, $

“*”表示复共轭. 方程(12)的解可写为

$ \psi_1(X_0,X_1) = \frac{{\rm d}B}{{\rm d}X_1}H(X_0), $

则有$ L_0 H = p'(X_0) $, 此时Fredholm条件也自动成立.
将$ \psi_0 $和$ \psi_1 $代入方程(13)得

$\begin{split} L_0\psi_2 =\, & \frac{{\rm d}^2B}{{\rm d}X_1^2}\left[H'(X_0)+\frac{1}{2}p(X_0)\right]-gB^3p^3(X_0) \\ & -\hat{\eta}B^5p^5(X_0)+\mu_2Bp(X_0),\\[-15pt] \end{split}$

再次应用Fredholm条件, 得

$ -D\frac{\partial^2B}{\partial X_1^2} - \mu_2 B +g\alpha _1B^3+\hat{\eta}\alpha {}_2B^5 = 0, $

其中$ \alpha_1 = \dfrac{\displaystyle\int_0^{4 K(q)}p^4(x){\rm d}x}{\displaystyle\int_0^{4 K(q)}p^2(x){\rm d}x} $, $ \alpha_2 = \dfrac{\displaystyle\int_0^{4 K(q)}p^6(x ){\rm d}x}{\displaystyle\int_0^{4 K(q)}p^2(x){\rm d}x} $均为正数, $ D = \dfrac{{\langle {p/2 + H', p} \rangle }}{{\langle {p, p} \rangle }} = \dfrac{1}{2}\dfrac{{{{\rm d}^2}\mu }}{{{\rm d}{k^2}}}\left| {_{k = {k_0}}} \right. $^[38]. 方程(15)为一定态NLSE, 它有局域型孤立波解

$ B(X_1) = \frac{\pm{\rm sech}(\beta X_1)}{\sqrt{\delta+\gamma{{\rm sech}^2}(\beta X_1)}}, $

其中$ \beta = \sqrt{\dfrac{-2\mu_2}{D}} $, $ \delta = \dfrac{\sqrt{3 g^2\alpha_1^2+16\hat{\eta}\alpha_2\mu_2}}{2\sqrt{3}|\mu_2 |} $, $\gamma = $ $ \dfrac{g\alpha_1- 2\delta\mu_2}{4\mu_2} $. 当$ \hat{\eta} = 0 $时, 方程(15)退化为3次NLSE; 相应地, 解(16)也退化为3次NLSE的孤立波解.
在$ X_1 = 0 $处$ B(X_1) $取得极值, 故调制波的振幅$ A\sim\dfrac{1}{\sqrt{\delta+\gamma}} = \sqrt{\dfrac{4\mu_2}{g\alpha_1+2\mu_2\delta}} $. 由此可知在其他参数给定时, 随着非线性相互作用$ g, \eta $的增大, 带隙孤子的振幅会单调递减, 后面的数值结果也证实了该结论. 另外, 为保证β为实数, 则$ \mu_2 $与D须异号, 即$ {\rm sgn}(\mu_2) = -{{\rm sgn}}(D) $, 这里$ {\rm sgn}(x) $为符号函数. 该式表示此时的μ必位于带隙内, 即孤子确为带隙孤子. 为保证解不存在奇性, 须δ为正实数且$ \delta+ $ $\gamma > 0 $, 则有$ \hat{\eta}\mu_2 > -\dfrac{3 g^2\alpha _1^2}{16\alpha _2} $. 举例来说, 如果要研究半无界带隙内的孤子, 则此时$ D > 0, \mu_2 < 0 $, 若取$ g, \eta < 0 $则两个条件必同时满足(充分非必要); 但若研究第一带隙内接近第一Band的孤子, 则$ \mu_2 > 0, D < 0 $, 取$ g, \eta > 0 $则两个条件必同时满足(充分非必要). 以上分析对于数值研究带隙孤子有一定的指导意义.
2 -->

2.3.牛顿共轭梯度法寻找带隙孤子

牛顿共轭梯度法(Newton-Conjugate-Gradient, NCG)^[41]是一种高效的数值方法, 可用来求解非线性演化方程的孤立波解. 该方法的主要思路是用牛顿迭代法结合共轭梯度法求解所得线性方程, 其收敛速度比共轭梯度法和牛顿法等其他现有的迭代法要快, 而且容易编程实现. 文献[42]也提出了牛顿法与共轭梯度法的组合方法, 并证明该方法的全局收敛性. NCG法是解决无约束最小优化问题的方法之一^[43]. 文献[44]中运用该方法讨论了路径约束动力学演化过程的优化问题. NCG法在一些物理模型上得到了广泛应用, 例如具有周期外势或无周期外势的二维非线性薛定谔方程, Kortewegde Vries (KdV)方程和五阶Kadomtsev Petviashvili (KP)方程的求解等^[41]. 它可用来寻找各种物理系统(如非线性光学, BEC和水波等)中的孤立波, 而且既能寻找基态也能求解激发态, 可以作为处理这类问题的首选方法.
下面采用NCG寻找方程(6)的带隙孤子. 计算时, 需要将无穷区间$ x\in(-\infty, +\infty) $截断为有限区间. 通常可取外势的足够多个周期作为计算的空间范围, 然后对该区间进行离散化后便可应用NCG求解. 计算表明, 所得结果对迭代初值有一定的依赖性(但并不十分敏感). 若迭代初值选择不当, 则迭代过程会发散或者收敛于平凡解$ \psi = 0 $. 计算时, 取以下多个高斯波包的叠加

$ \psi(x) = \sum\limits_{j = 1}^{{\cal N}}{A_j{\rm e}^{-(x-x_j)^2/W_j}} $

作为迭代初始条件, 这里$ A_j, x_j, W_j $分别表示开始迭代时第j个波包的振幅、中心位置和宽度, $ {\cal N} $表示总波包数. 通过选择合适的$ A_j, x_j, W_j, {\cal N} $, 便可得到所需带隙孤子. 若欲寻找结构更为复杂的带隙孤子, 可尝试选用其他形式的迭代初值.
图1显示了$ q = 0.1 $和0.99时不同参数条件下的单峰带隙孤子, 迭代时取$ {\cal N} = 1,$ A₁ = 0.6, $ x_1 = 0, $W₁ = K(q) . 该孤子的峰值位于外势的最低点处, 一般称其为on-site soliton. 图1(a)和图1(b)为$ q = 0.1 $时吸引相互作用下半无界带隙中的带隙孤子, 图1(c)和图1(d)为$ q = 0.99 $时排斥作用下第一带隙中的带隙孤子. 从图1可以看出, 当三体相互作用强度不变而两体相互作用变强时, 带隙孤子的振幅明显降低; 也可看出, 当两体相互作用强度不变而三体相互作用变强时, 带隙孤子的振幅也会变小. 这与前面理论分析结论一致. $ q = 0.1 $时孤子的结构较为复杂, 但$ q = 0.99 $时孤子为钟形, 故外势模数对孤子结构有一定影响. 为进一步研究孤子振幅随相互作用强度的变化, 定义$ A = \max{|\varPsi|} $为孤子振幅. 图2(a)和图2(b)分别显示了在吸引相互作用下半无界带隙内和排斥相互作用下第一带隙内单峰孤子振幅随相互作用强度的变化. 可以看出, 随着两体相互作用强度|g|和三体相互作用强度|η|的增大, 孤子的振幅确在单调递减.
图 1 3?5次GPE的单峰带隙孤子($V_0=4$)　(a), (b) $q=0.1$; (c), (d) $q=0.99$. 阴影部分表示外势$V(x)$低处
Figure1. Profiles of singl-hump gap solitons of the cubic-quintic GPE ($V_0=4$): (a), (b) $q$ = 0.1; (c), (d) $q$ = 0.99. Shaded regions represent lattice sites, i.e., regions of low potential values $V(x)$.

图 2 单峰带隙孤子的振幅随相互作用强度的变化($V_0=4$)
Figure2. Amplitudes of single-hump gap solitons v.s. nonlinear interaction strength.

图3显示了不同参数情形下图1中on-site孤子三体和两体相互作用能量的比值, 其中$ E_{{\rm three}} = \dfrac{1}{3}\eta{\displaystyle\int|\psi(x)|^6}{\rm d}x $, $ E_{{\rm two}} = \dfrac{1}{2}g{\displaystyle\int|\psi(x)|^4}{\rm d}x $. 可以看出, 当η固定时, |g|越大, 比值越小, 表明二体相互作用越占优; 但当g固定时, 随着|η|的增大, 比值也在增大, 表明三体相互作用逐渐增强. 图3(a)中两种能量已处于可比拟的范围; 图3(b)中η足够大时, 三体相互作用能量已超过两体相互作用能量. 这种情况下, 三体相互作用更不可忽略.
图 3 on-site孤子三体相互作用能量和两体相互作用能量的比
Figure3. The ratio of three-body energy to two-body energy for on-site solitons with different interaction strength.

图1中所示的带隙孤子是结构最为简单的一类基本孤子. 除此之外, 也存在另一类如图4所示的双峰孤子. 这类孤子具有两个波峰, 两波峰间的中心位置在外势的最高点处, 一般称其为off-site soliton. 图4(a)和图(b)所示两波峰具有相同的相位, 称为同相偶极孤子(迭代时取$ {\cal N} = 2$, A₁ = A₂ = 0.6, x₁ = 0, x₂ = –2K(q), W₁ = W₂ = K(q)); 图4(c)和图4(d)所示孤子则被称为异相偶极孤子(迭代时取$ {\cal N} = 2$, A₁ = 0.6, A₂ = –0.6, x₁ = 0, x₂ = –2K(q), W₁ = W₂ = K(q)). 该类孤子的振幅也随着|g|和|η|的增大而减小. 从波形结构上来看, off-site孤子似乎可以看作是两个单峰on-site孤子的组合; 但从动力学稳定性来看, on-site孤子始终稳定而同相偶极孤子始终不稳定, 故同相偶极孤子也被视为另一类基本孤子. 异相偶极孤子的稳定性则较为复杂, 一般依赖于系统的参数.
图 4 3?5次GPE的同相偶极孤子和异相偶极孤子($V_0=4$)　(a), (b) $q=0.1; $ (c), (d) $q=0.99$. 阴影部分表示外势$V(x)$低处
Figure4. Profiles of double-hump gap solitons of the cubic-quintic GPE ($V_0=4$) (a), (b) $q$ = 0.1; (c), (d) $q$ = 0.99. Shaded regions represent lattice sites, i.e., regions of low potential values $V(x)$.

孤子的稳定性无论在理论上还是在实验上都是一个重要的问题, 下面将用非线性动力学演化的方法来研究前面所得带隙孤子的稳定性. 在初始时刻, 给波函数一小扰动. 若经过足够长时间演化后孤子的振幅和波形没有发生明显变化, 则可认为该孤子动力学稳定; 否则认为它动力学不稳定. 数值计算时, 所用方法为时间劈裂傅里叶谱方法^[45]. 该方法具有效率高、精度高、计算过程中粒子数守恒, 而且程序容易实现等优点.
对初始波函数做下面类型的扰动

$ \varPsi(x,0) = \psi (x){{\rm e}^{{\rm i}\beta x}}, $

其中$ \psi(x) $是前面所得定态带隙孤子, $ \beta\ll1 $是扰动参数, 计算时取$ \beta = 0.01 $.
图5为不同参数情形下不同类型带隙孤子动力学演化的等值线图($ V_0 = 4 $). 可以看出, 不同类型孤子的稳定性有所不同. 图5(a)为$ \mu = 1$, q = 0.1, g = –1, $\eta = -1 $ (对应图1(a)中情形)时on-site孤子的动力学演化. 可以看出, 经过一段时间演化后波形和振幅都没有发生明显变化, 故它是动力学稳定的. 数值研究中, 通过对各种参数情形进行计算, 结果均表明无论在第一带隙还是半无界带隙中, on-site孤子始终稳定. 图5(b)为$ \mu = 0.5$, q = 0.1, g = –1, $\eta = -1 $ (对应图4(a)中同相偶极孤子)的动力学演化. 很明显孤子演变成振荡状态, 即能量在两个相邻光晶格之间周期性地转移. 初始时刻孤子空间结构的对称性随时间演化时被破坏, 所以这是一种振荡型不稳定. 大量计算表明, 这种同相偶极孤子总是动力学不稳定的. 即使q较大时, 尽管两波峰相距较远, 它仍然表现出不稳定(但不稳定性变弱). 所以这种类型的孤子不宜视为两个on-site孤子的组合. 图5(c)为$ \mu = 2 $, q = 0.5, g = 1, $\eta = 1 $时第一带隙内异相偶极孤子(带隙孤子结构与图4(c)和图4(d)中类似)的动力学演化. 可以看出, 此时的孤子仍具有振荡不稳定性. 图5(d)为$ q = 0.99 $(其他参数与图5(c)中相同)时异相偶极孤子的动力学演化. 可以看出, 此时的带隙孤子是动力学稳定的. 这一结果表明外势的模数会影响带隙孤子的稳定性; 同时也意味着其他参数固定时, 存在临界值$ q_{\rm c} $. 当$ q\geqslant q_{\rm c} $时, 异相偶极孤子稳定; 而当$ q < q_{\rm c} $时, 异相偶极孤子不稳定. 对于图5(c)和图5(d)参数情形, $ q_{\rm c}\approx0.968 $. 在g, η, $ V_0 $等参数不变的情况下, 临界值$ q_{\rm c} $随μ值而变. 例如对于$ g = -1 $, $ \eta = -1 $, $ V_0 = 4 $时的异相偶极孤子, 当$ \mu = 0.55 $时$ q_{\rm c}\approx0.956 $; 而当$ \mu = 0.65 $时, $ q_{\rm c}\approx0.754 $.
图 5 不同类型带隙孤子的动力学演化　(a) on-site孤子; (b) 同相偶极孤子; (c), (d) 异相偶极孤子
Figure5. Contour plots of $|{\varPsi(x, t)}|$ for perturbed gap solitons: (a) On-site soliton; (b) in-phase dipole soliton; (c), (d) out-phase dipole soliton.

值得说明的是, 前面仅讨论了两类基本结构的带隙孤子, 即on-site孤子与off-site孤子. 实际上, 方程(6)存在无穷多结构各异的带隙孤子, 如三峰、四峰等结构更为复杂的孤子. 很多结构复杂的带隙孤子, 可视为这两类基本带隙孤子的组合. 这类孤子的特点是它们的峰值总位于外势的最低点处. 数值结果表明, 方程(6)也存在另外一类带隙孤子, 它们的峰值并不位于外势的最低点处. 文献中将该类孤子称为亚基本带隙孤子(sub-fundamental gap soliton)^[46]. 图6(a)和图6(c)显示了q = 0.1, g = –1, $\eta = -2,$ V₀ = 4不同μ值时的亚基本带隙孤子, 它有着和前面两种类型的孤子明显不同的结构特征. 这类孤子在空间上呈奇对称分布, 它的中心位置位于外势的最低点处, 波峰则介于外势最高点与最低点之间. 图6(b)和图6(d)为该类带隙孤子的动力学演化, 可以看出, $ \mu = 2.9 $时孤子不稳定而$ \mu = 3 $时孤子稳定. 该结果意味着存在临界值$ \mu_{\rm c} $, 当$ \mu < \mu_{\rm c} $时孤子不稳定而$ \mu > \mu_{\rm c} $时孤子稳定. 数值结果表明在该组参数下, $ \mu_{\rm c}\approx2.992 $. 这里取$ \psi(x) = A_1 {\rm e}^{-x^2/W_1}\sin({k_1}x) $作为迭代初始条件, 计算时取 A₁ = 1, W₁ = 4K(q), k₁ = 1.
图 6 (a), (c) 第一带隙中的亚基本带隙孤子(红色). 蓝线表示外势; (b), (d) 亚基本带隙孤子的动力学演化
Figure6. (a), (c) Profiles of sub-fundamental gap solitons (red lines) lie in the first bandgap. The solid blue lines denote the external potential; (b), (d) contour plots of $|{\varPsi(x, t)}|$ for perturbed gap solitons.

三体相互作用强度对于带隙孤子的稳定性也有一定影响. 图7显示了V₀ = 2, μ = 1.25, q = 0.1, g = 1而不同η时第一带隙中同相偶极孤子的动力学演化. 图7(a)中$ \eta = 0 $, 可以看出, 经过一段时间演化后带隙孤子呈现出动力学不稳定; 图7(b)中$ \eta = -0.2 $, 可以看出此时的带隙孤子是动力学稳定的. 表明三体相互作用强度对于带隙孤子的稳定性确实有着一定影响, 它可以改变带隙孤子的稳定性区域, 该结论与文献[18]中结论一致. 注意这里的同相偶极孤子与图3中有所不同, 它属于文献[29]中提到的(1, 1)结构, 并非从Band中分岔出来(迭代时取$ {\cal N} = 2$, A₁ = A₂ = 0.6, x₁ = 0, x₂ = –2K(q), W₁ = W₂ = K(q) ).
图 7 第一带隙中同相偶极孤子　(a) $\eta=0$和(b) $\eta=-0.2$时的动力学演化
Figure7. Contour plots of $|{\varPsi(x, t)}|$ for in-phase dipole solitons lie in the first bandgap with (a) $\eta=0$ and (b) $\eta=-0.2$.

本文研究了准一维情形下束缚于雅可比椭圆函数周期势中具有三体相互作用BEC系统中的带隙孤子及其稳定性. 在平均场近似下, 系统的动力学行为可由3—5次GPE描述. 首先利用多重尺度法对系统的带隙孤子进行了理论分析, 将3—5次GPE化为一定态NLSE, 并给出了一类带隙孤子的解析表达式. 然后采用NCG法数值得到了该系统中存在的带隙孤子, 包括两类基本带隙孤子(on-site孤子和off-site孤子)与亚基本带隙孤子(sub-fundamental gap soliton). 数值结果与理论分析均表明, 三体相互作用强度的增大将会导致带隙孤子的振幅减小. 最后, 用时间劈裂傅里叶谱方法对GPE进行长时间动力学演化以考察孤子的稳定性, 发现on-site孤子始终稳定, 同相偶极孤子和异相偶极孤子既有稳定的也有不稳定的. 三体相互作用强度对于带隙孤子的稳定性也有一定影响. 外势的模数对孤子的结构和稳定性会产生明显影响, 故实验中可通过调整外势的模数来改变带隙孤子的稳定性.