1.   引　言

玻色-爱因斯坦凝聚体(Bose-Einstein condensate, BEC)由Bose和Einstein于1924年提出. 随后1995年科罗拉多大学(NIST-JILA) 的Cornell和Wieman等首次采用蒸发冷却技术实现中性铷原子BEC, 并于2001年获得诺贝尔物理学奖^[1]. 近年来, 大量的理论和实验工作致力于研究超冷玻色气体的动力学性质^[2-4]. 在极低温下, BEC的动力学行为可以用Gross-Pitaevskii方程(GPE)^[5-7] 来描述. 该方程是著名非线性薛定谔方程 (nonlinear Schrödinger equation, NLSE)^[8]的一种, 也是描述非线性色散介质演化的通用模型 ^[9-11]. 在低浓度BEC中, 原子间距离较大, 只需考虑两体相互作用. 然而, 对于浓度较高的情况, 三体相互作用变得有意义, 使GPE 具有五次非线性项 ^[12,13]. 随实验技术的不断发展, 三体相互作用的研究也得到了科研工作者的广泛关注 ^[14,15]. 由Feshbach共振^[16,17]调控的三体相互作用导致系统出现了与仅考虑两体相互作用不同的物理性质.

过去几十年里, 孤子在不同科研领域得到了大量研究. 在实验和理论研究中发现了亮孤子^[18]、暗孤子^[19]、灰孤子^[20]和涡旋孤子^[21,22]等. 最近, Makris等^[23]提出了存在于均匀介质界面的表面波, 这是一种沿不同介质界面传播的导波. 在固体物理^[24]、非线性光学^[25]和近场光学^[26,27]等领域内, 研究的一个重点问题就是非线性表面波^[28,29]. 非线性界面的存在极大地改变了非线性模型的特性, 产生了许多新奇的现象, 比如它使孤子结构变得复杂, 且改变孤子的存在区域和稳定条件. 已有许多工作涉及到周期势中表面带隙孤子的理论研究^[30] 和实验观察^[31,32]. 尽管表面带隙孤子已在不同领域中得到了深入研究, 但在均匀介质与半无限雅可比椭圆正弦势相互作用的界面上, 准一维BEC 中表面带隙孤子的存在性和稳定性仍然是一个开放而富有挑战性的问题.

本文对准一维BEC中表面带隙孤子及其稳定性进行了数值研究. 首先用布洛赫定理得到了系统的能带结构, 然后用牛顿共轭梯度法得到了表面孤子, 接着用线性稳定性分析和非线性动力学演化研究了孤子的稳定性, 结果表明三体相互作用会明显影响表面亮孤子的稳定性. 考虑到三体相互作用时, 系统中还会出现表面扭结孤子和气泡孤子. 表面扭结孤子既有稳定的也有不稳定的, 但表面气泡孤子均不稳定.

2.   理论模型

在平均场近似下, 具有两体和三体相互作用BEC的动力学行为可用GPE描述^[33]:

$$\begin{equation} \begin{split} \begin{aligned} {\rm{i}}\hbar\frac{\partial\varPsi}{\partial t}=-\frac{\hbar^2}{2m}{\nabla^2}\varPsi+V(\boldsymbol r)\varPsi+g_2|\varPsi|^2\varPsi+g_3|\varPsi|^4\varPsi, \end{aligned} \end{split} \end{equation}$$

(1)

式中, $\varPsi=\varPsi({\boldsymbol r}, t)$是波函数, $\hbar$是约化普朗克常数, m是原子质量, $V({\boldsymbol r})$表示外势. $g_2$和$g_3$分别表征两体和三体相互作用强度. $g_2=4\pi\hbar^2 a_{\rm{s}}/m$, 当$a_{\rm{s}} > 0$时, 原子间相互排斥;$a_{\rm{s}} < 0$时, 原子间相互吸引. g₃表示为^[34]

$g_3 = \dfrac{12\pi\hbar^2 a_{\rm{s}}^4}{m}\left[d_1 + d_2\mathrm{tan}\left(s_0\ln\dfrac{|a_{\rm{s}}|}{|a_0 |} + \dfrac{\pi}{2}\right)\right]$

其中, 常数$d_1, d_2, a_0$和$s_0$的值在文献[34, 35]中给出. $g_3$可以通过Feshbach共振技术在$(-\infty, +\infty)$内调控^[16,17,34]. 总粒子数为$N=\displaystyle \int|\varPsi|^2 {\rm{d}}^3\boldsymbol r$. 在实验中,BEC通常被约束在谐振子势$V(\boldsymbol r)=\dfrac{1}{2}m(\omega_x^2 x^2+ \omega_y^2 y^2+\omega_z^2 z^2)$中, $\omega_x$, $\omega_y$和$\omega_z$分别是$x, y$和z方向的角频率. 在饼状凝聚体中, $\omega_x\approx\omega_y$且$\omega_z\gg\omega_x$, 此时三维GPE可通过准二维近似变为二维GPE. 在雪茄状凝聚体中, $\omega_y, \omega_z\gg\omega_x$, 三维GPE 可通过准一维近似变为一维GPE. 本文仅考虑雪茄状凝聚体, 引入无量纲量

$$\begin{split} &\tilde t=\omega_xt,\;\tilde x=x/a_{\mathrm{ho}}, \tilde V = V/\hbar\omega_x, \;\tilde \varPsi=\frac{\varPsi}{\sqrt{n_0/a_{\mathrm{ho}}^3}},\\ &\tilde g_2=\frac{2a_{\rm{s}}n_0\sqrt{\omega_z\omega_y}}{a_{\mathrm{ho}}\omega_x},\;\tilde g_3=\frac{n_0\omega_z\omega_y}{3\pi^2a_{\mathrm{ho}}^3\hbar\omega_x^3}{g_3}, \\[-15pt]\end{split}$$

(2)

其中$a_{\mathrm{ho}}=\sqrt {\dfrac{\hbar }{{m{\omega _x}}}}$是谐振子特征长度, $n_0$是给定的粒子数密度. 此时总粒子数$N = {n_0}\displaystyle\int {{{| {\tilde \varPsi ({\tilde x, \tilde t})} |}^2}} {\rm{d}}\tilde x$. 利用文献[36]中所用方法并省略变量上的“$\sim$”号, 得以下准一维无量纲GPE

$${\rm{i}}\frac{\partial\varPsi}{\partial t}=-\frac{1}{2}\frac{\partial^2\varPsi}{\partial x^2}+V(x)\varPsi+g_2|\varPsi|^2\varPsi+g_3|\varPsi|^4\varPsi.$$

(3)

选择外势形式为^[30-37]

$$\begin{equation} \begin{split} \begin{aligned} V(x)=\left\{ \begin{aligned} &{V_0}{\rm{s}}{{\rm{n}}^2}\left( {x,q} \right), &x \geqslant 0,\\ &0, &x < 0, \end{aligned} \right. \end{aligned} \end{split} \end{equation}$$

(4)

其中${\rm{sn}^2}(x, q)$表示模数为$q(0\leqslant q < 1)$的雅可比椭圆正弦函数. 该势可看作是三角函数的推广, 实验中仅用两束激光便可很好地近似它^[38]. 当$q < 0.9$时, 其近似程度可达99%. 图1显示了当$q=0, 0.8, 0.99$时$V(x)$的图像, 其中$V_0=-2$. 当q接近1(但小于1)时会出现平台状结构. 图1插图显示了$x\geqslant0$时外势周期$2 K(q)$随模数q的变化($K(q)$是第一类完全椭圆积分). 可以看到$V(x)$的周期随着q从0增加到几乎为1时单调增加, 特别是当q接近1时迅速增加. 例如, 当$q=0$时, $V(x)$的周期为π, 而当$q=0.99$时, 周期约为6.7132. 后者是前者的两倍多.

3.   带隙结构

寻找方程(3)形式为$\varPsi(x, t)=\psi(x){\rm{e}}^{-{\rm{i}}\mu t}$的定态解, 其中化学势μ和实函数$\psi(x)$满足非线性微分方程

$$\begin{equation} \frac{1}{2}{\psi_{xx}}+\mu\psi-V(x)\psi-g_2\psi^3-g_3\psi^5=0. \end{equation}$$

(5)

当$\psi(x)$为无穷小时, 方程(5) 中的非线性项$\psi^3$和$\psi^5$可忽略不计, 得到一线性方程

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

(6)

如果$V(x)$在$x\in{\mathbb{R}}$上是周期函数, 那么方程(6)为广义马丢(Mathieu)方程, 其有界解称为布洛赫模(Bloch modes), 相应的化学势μ构成布洛赫带(Bloch bands). 基于布洛赫定理, 方程(6) 的有界解可以写为 ^[39]

$$\begin{equation} \psi=p(x)={{\rm{e}}^{{\rm{i}}kx}}\tilde{p}(x;\mu), \end{equation}$$

(7)

其中$\tilde{p}(x;\mu)$与外势$V(x)$具有相同的周期, $\mu=\mu(k)$是一维色散关系. Bloch带和色散关系可以用傅里叶配置法(Fourier collocation method)^[40]得到. 图2(a)显示了$q=0$时带隙结构随$V_0(\neq0)$的变化, 阴影区域和空白区域分别代表能带(band)和能隙(gap). 图2(a)中垂直虚线分别表示图2(b)—(d)中用到的$V_0$值. 从图2(b)—(d)可以看出布洛赫带被许多间隙分开, 这些间隙称为带隙(band-gap)^[41]. 通过比较不同$V_0$下的第一能带, 发现色散曲线的曲率随$|V_0 |$的增大而减小. 例如, 当$V_0=-8$时, 能带曲线变得非常平坦, 从而${\rm{d}}^2\mu/{\rm{d}} k^2 \to 0$. 线性Bloch波是能带的本征模, 它们不能存在于带隙中, 但非线性波(特别是孤子)有可能在带隙中激发. 下面用NCG方法数值地寻找孤子, 并讨论其稳定性. 关于NCG方法的详细介绍可参看文献[40-42].

4.   表面孤子及其稳定性

4.1   表面亮孤子

考虑化学势从能带边缘进入带隙时线性布洛赫波不再存在, 但可激发出孤立波即表面孤子^[30]. 采用NCG方法对表面亮孤子进行数值求解, 计算时迭代初始值取为

$$\begin{equation*} \psi(x)=\sum\limits_{j=1}^{\mathcal{N}}A_j {\rm{e}}^{-(x-x_j)^2/W_j}, \end{equation*}$$

其中$\mathcal{N}$为波包总数, $A_j, W_j$和$x_j$分别表示第j个高斯波包的振幅、宽度和中心位置. 通过调节$\mathcal{N}$, $A_j$, $W_j$和$x_j$, 有望得到不同结构类型的孤子.

图 3显示了用NCG方法得到的不同三体相互作用强度$g_3$下的表面孤子, 其中阴影区域表示外势$V(x)$低处, 图 3(a), (c)是单峰表面亮孤子, 图 3(b), (d)为双峰表面亮孤子. 可以看到, 表面亮孤子既可存在于半无限带隙中也可存在于第一带隙中, 也可看出它的振幅随三体相互作用强度$|g_3 |$的增大而减小, 这里孤子的振幅为$A=\max|\psi|$. 数值计算时, 图 3(a), (c)中$\mathcal{N}=1,\; A_1=0.6, \;x_1= K(q), W_1= K(q)$. 图 3(b), (d)中$\mathcal{N}=2, \;A_1= A_2=0.6, x_1=K(q),\;x_2=3 K(q),\; W_1=W_2=K(q)$. 表面亮孤子只有在$\mu < 0$时才能存在, 这是因为均匀介质中尾波呈指数衰减形式 ^[30]. 这一点也可做如下说明. 考虑$x\to-\infty$时$\psi(x)\to0$, 则方程(5)中非线性项可以忽略, 其退化为$\psi_{xx}+2\mu\psi=0$. 当$\mu < 0$时, 该方程的解$\psi(x)\sim {\rm{e}}^{\sqrt{-2\mu}x}$. 而当$\mu > 0$时, 其解为$\psi(x)\sim \cos\sqrt{2\mu}x$. 后者表示均匀介质中传播的波数为$\sqrt{2\mu}$的线性单色波, 其波长为$2\pi/\sqrt{2\mu}$.

系统中总粒子数为$N=n_0 P$, $P=\displaystyle\int{|\psi|^2 {\rm{d}}x}$. 在非线性光学中, P通常被称为功率或能量. 表面亮孤子只有当P超过某一阈值时才存在, 且其并非从布洛赫带的边缘分叉而出^[43]. 图4(a)显示了P随化学势μ的变化, 图4(b)—(g)为图4(a)中对应标记点处表面亮孤子的波形. 可见, 在半无界带隙, 三个分支从$\mu\approx-1.2775$处分岔而出; 而在第一带隙, 它们从$\mu\approx-0.99876$处分岔而出. 三个分支均非第一能带边缘分叉而出, 这与文献[43]中结论一致. 图4(a)中, 实线表示稳定分支, 虚线为不稳定分支(孤子的稳定性通过线性稳定性分析得到, 将在4.2节讨论). 最上面和中间两支是完全不稳定分支, 这里所有表面亮孤子都不稳定. 低处分支化学势μ存在一个或多个临界值. 例如, 在图中所取参数情形下, 半无限带隙中存在临界值$\mu_1^{\rm{c}}$, 仅当$\mu\leqslant\mu_1^{\rm{c}}$时表面亮孤子才稳定; 第一带隙中存在临界值$\mu_2^{\rm{c}}$和$\mu_3^{\rm{c}}$, 表面亮孤子在$\mu_2^{\rm{c}}\leqslant \mu\leqslant\mu_3^{\rm{c}}$时才稳定(图4(a)中部分图没显示), 否则存在扰动增长率$\lesssim 10^{-2}$的弱不稳定. 临界值也依赖于其他参数, 如$V_0$, $g_2$, $g_3$等. 在图4(a)所用参数情况下, 数值结果给出$\mu_1^{\rm{c}}\approx-1.4938,\; \mu_2^{\rm{c}}\approx-0.7185,\; \mu_3^{\rm{c}}\approx-0.57$.

4.2   表面暗孤子

4.1节讨论的表面亮孤子, 其满足边界条件$\psi(\pm\infty)=0$. 除这种情况外, 三体相互作用的BEC系统中还可存在另一种类型的表面孤子, 当${x\to\pm\infty}$时波函数并不趋于零, 即$\psi(-\infty)\neq0$或$\psi(+\infty)\neq0$. 一般地, 若$\psi(-\infty)\neq\psi(+\infty)$, 则称其为扭结孤子(kink); 若$\psi(-\infty)=\psi(+\infty)\neq0$, 则称其为气泡孤子(bubble). 这两类孤子统称暗孤子. 前一种情形下, 在两种介质分界面处会出现扭结结构, 称该类孤子为表面扭结孤子. 考虑$x\to-\infty$时$\psi(x)\to$常数, 则方程(5)化为$\mu\psi=g_2\psi^3+g_3\psi^5$, 由此可得$\psi_1=0$或$\psi_2^2=(-g_2\pm\sqrt{g_2^2+4\mu g_3})/2 g_3$, 后式要求$\mu g_3\geqslant q-g_2^2/4$ 且结果为正. 这些结论对于数值寻找该类孤子有着指导意义. 值得注意的是, 如果不考虑三体相互作用($g_3=0$), 则后一种情况不存在, 表明此时BEC系统中不存在表面暗孤子. 数值结果也证实了该结论.

用NCG方法寻找表面扭结孤子时, 迭代初值为

$$\psi(x) = a_1 (\tanh x + 1 ) [1 + a_2\cos2\pi(x-x_0)/T_x ]$$

(8)

或

$$\begin{equation} \psi(x) = \left\{ \begin{aligned} &a_3\tanh x\left[1+a_4\cos2\pi(x-x_0)/T_x\right],&x > 0, \\ &a_5, &x\leqslant0, \end{aligned}\right. \end{equation}$$

(9)

其中$T_x$表示外势的周期. (8)式和(9)式分别适合于寻找$\psi(-\infty)=0$和$\psi(-\infty)\neq0$的扭结孤子, 可以通过调节参数$a_i(i=1, 2, \cdots, 5), x_0$得到满意的结果. 一般地, (9)式中可以取$a_3, a_5=\psi_2$或$-\psi_2$. 实际计算表明, 数值结果对迭代初值并不十分敏感. 图5给出了NCG方法得到的不同参数下这两种扭结孤子的波形. 图5(a)是在半无界带隙中的扭结孤子, 可以看出, 随着三体相互作用强度$|g_3 |$ 的增加, 孤子的振幅逐渐减小, 这一点与亮孤子类似. 图5(b)是在第一带隙中的扭结孤子, 在其他参数固定的情况下, 随着μ的减小, 孤子的振幅也随之减小. 迭代初值中参数取为$a_1=a_2=a_4=0.8, x_0=T_x/2$. 后续稳定性分析表明, 图5中表面扭结孤子均线性和非线性动力学稳定.

图6为不同非线性相互作用强度$g_2, g_3$下扭结孤子振幅A随化学势μ 的变化, 其中子图为标记点处孤子的波形. 该类扭结孤子的振幅随着化学势μ的减小而减小; 它们不仅存在于带隙中, 也存在于能带中, 且在μ大于某一临界值时才存在. 如当$g_2=-1, g_3=0.5$时(图6(a)), 它们存在于$\mu\geqslant -1.622$的范围内, 该临界值依赖于非线性相互作用$g_2$和$g_3$.如图6(b)所示, 扭结孤子存在于第一能带和第一带隙中, 而在半无界带隙中则不存在. 从图6也可看出, 随三体相互作用强度$|g_3 |$的增大, 表面扭结孤子的振幅会随之减小, 这一结论与图5一致. 图6(a)中标记的黑点为相变点, 表示μ从小到大增加时, 一旦越过该点, 扭结孤子将由负无穷远处为0的类型变为另一种类型. 该相变点也依赖于非线性作用强度$g_2$和$g_3$.

4.3   线性稳定性分析

在非线性理论中, 孤子的稳定性是一个很重要的问题, 下面对表面带隙孤子进行线性稳定性分析. 作Bogoliubov展开

$$\begin{split} \varPsi(x,t)=\;&\{\psi(x)+[\upsilon(x)+w(x)]{\rm{e}}^{\lambda t}\\ &+[\upsilon^*(x)-w^*(x)]{\rm{e}}^{\lambda^*t}\}{\rm{e}}^{-{\rm{i}}\mu t}, \end{split}$$

(10)

其中$|\upsilon|, |w|\ll1$是无穷小扰动, λ是本征值, “*”表示复共轭. 若所有本征值的实部均非正数, 则孤子是线性稳定的. 但若至少有一个本征值的实部大于0, 则给定的孤子在线性化意义下不稳定. 将(10)式代入(3)式, 在线性意义下得到以下本征问题

$$\begin{split} &\begin{bmatrix} 0 & L_1 \\ L_2 & 0 \end{bmatrix} \begin{bmatrix} \upsilon \\ w \end{bmatrix} =-{\rm{i}}\lambda \begin{bmatrix} \upsilon \\ w \end{bmatrix},\\ & L_1=\dfrac{1}{2}\partial_{xx}+\mu-V(x)-g_2\psi^2-g_3\psi^4,\\ &L_2= \dfrac{1}{2}\partial_{xx}+\mu-V(x)-3 g_2\psi^2-5 g_3\psi^4 .\end{split}$$

(11)

一般很难解析得到该本征问题的精确解, 但可用有限差分法或傅里叶配置法对它进行数值求解^[40].结果表明, 图4(a)最下面分支中孤子的稳定性依赖于物理参数, 即稳定和不稳定的表面亮孤子都可以存在. 三体相互作用强度$g_3$对孤子的稳定性有明显影响. 例如图7(a)中, 化学势$\mu=-1.5$位于半无限带隙, $g_2=-1, g_3=0$时表面孤子线性不稳定. 当$g_3=0.6$且其他参数不变时, 表面孤子的结构没发生明显变化, 但变为线性稳定(图7(b)). 又如图7(c)中, $\mu=-0.7$位于第一带隙, $g_3=0$时的孤子线性稳定; 而图7(d)中$g_3=1$时的孤子则线性不稳定.

为清楚地看出$g_3$对孤子稳定性的影响, 我们系统地做了大量数值计算. 图8 显示了不同化学势μ下, 表面带隙孤子扰动最大增长率$\lambda_{\rm{m}} = \max [{\mathop{\rm Re}\nolimits} (\lambda)]$随$g_3$的变化. 在半无界带隙中(图8(a)), 两体相互作用为吸引的情形下($g_2=-1$), 当$g_3$增大到某临界值$g_3^{\rm{c}}$时(其他参量不变), 孤子从不稳定变为稳定. 在第一带隙中(图8(b)), 两体作用为排斥的情形($g_2=1$), 当$g_3$增大到某临界值$g_3^{\rm{c}}$时, 孤子从稳定变为不稳定. 临界值$g_3^{\rm{c}}$依赖于化学势μ与两体作用强度$g_2$. 如当$\mu=-1.52$时, $g_3^{\rm{c}}\approx-2.5$, 但当$\mu= -1.55$时, $g_3^{\rm{c}}\approx -6.7$. 因此, 在实验中可以通过调节三体相互作用强度来改变表面带隙孤子的稳定性. 值得注意的是, 图8中$\lambda_{\rm{m}}$很小($\sim 10^{-2}$), 意味着扰动的增长相当缓慢, 孤子的不稳定性很弱.

4.4   非线性动力学演化

前面分析了表面带隙孤子的线性稳定性, 下面采用时间劈裂傅里叶谱方法 ^[44,45]对(3)式进行长时间动力学演化, 以研究孤子的非线性动力学稳定性. 一般经过足够长时间演化后, 如果受扰动孤子的振幅与波形均不产生明显变化, 则孤子是动力学稳定的. 否则, 孤子是动力学不稳定. 实际计算中, 初始时刻对孤子施加了振幅为0.01的随机扰动.

图9为不同参数情况下表面亮孤子非线性动力学演化的等值线图, 初始时刻波形分别见图7子图. 在半无限带隙中($g_2=-1$), $g_3=0$时孤子不稳定(图9(a)), 而$g_3=0.6$时孤子稳定(图9(b)), 这与图7 中线性稳定性分析所得结论一致. 在第一带隙中($g_2=1$), $g_3=0$时孤子稳定(图9(c)), 而$g_3=1$时孤子不稳定(图9(d)), 也与线性稳定性分析结果一致.

图10(a)为$g_2=-1, g_3=0.3, \mu=-1.5$时, 半无界带隙中表面扭结孤子的线性稳定性谱, 子图为对应的波函数. 线性稳定性分析表明该扭结孤子稳定. 图10(b) 为受扰动扭结的非线性动力学演化结果, 也表明该表面扭结孤子动力学稳定, 与线性稳定性分析的结论一致. 我们用该方法系统地考察了图6中得到的表面扭结孤子的稳定性, 发现它们都是线性稳定.

通过仔细调整迭代初值, 我们也找到了表面气泡孤子, 其波形如图11(a)中子图所示. NCG方法寻找时, 迭代初值设为下列形式

$$\begin{split} \psi(x)=\;& b_1 [1+b_2\cos2\pi (x-x_1)/T_x ]\\ & \times [1-b_3 {\rm{e}}^{-(x-x_2)^2}]. \end{split}$$

(12)

通过调整参数$b_1, b_2, b_3, x_1, x_2$可以得到满意的结果. 图11(a)中, 取为$b_1=1.5, \;b_2=b_3=1,\; x_1=T_x/2, \;x_2= 3 T_x/2$. 图11(a), (b)分别为其线性稳定性谱和非线性动力学演化的等值线图, 两者均表明表面气泡孤子是不稳定的. 大量数值结果表明, 所有表面气泡孤子均不稳定, 这与文献[37,46]中所得结论一致. 此外, 也找到了一类不稳定的表面扭结孤子, 它们存在于μ的某个范围内. 图11(c), (d)分别显示了$g_2=-1, g_3=0.3, \mu=-1.5$时, 所得表面扭结孤子的线性稳定性谱和非线性动力学演化图, 其中子图为其波形. 两者均表明该扭结孤子是不稳定的. 寻找该孤子时, 迭代初值为(8)式中取$a_1= 0.5,\; a_2=0.4,\; x_0=T_x/2$.

外势模数q对表面孤子的稳定性也有一定影响. 图12(a)显示了第一带隙中表面亮孤子扰动最大增长率$\lambda_{\rm{m}}$随q的变化, 孤子的波形如图4(e)所示. 从图12(a)可以看出, 随着q的增加, 孤子逐渐从不稳定变为稳定. 对于给定的μ值, 存在一依赖于μ的临界值$q_{\mu}^{\rm{c}}$. 当$q > q_{\mu}^{\rm{c}}$时, 孤子稳定; 而$q < q_{\mu}^{\rm{c}}$时, 孤子不稳定. 前面图5所示的扭结型孤子, 数值结果表明该类孤子始终稳定, 其稳定性不受q的影响. 对图11中所示的气泡孤子和扭结孤子, 数值结果表明它们始终不稳定, 但增加q可以减弱其不稳定性. 图12(b)显示了不同化学势下表面气泡孤子(波形见图11(a))的$\lambda_{\rm{m}}$随q的变化. 可以看出, 当q增大时, $\lambda_{\rm{m}}$逐渐减小, 但始终不为零.

前面主要考虑了$V_0 < 0$的情况. 当$V_0 > 0$时, 系统中也存在类似的表面孤子. 图13(a)显示了$V_0= 2$时, 半无限带隙中的表面亮孤子, NCG方法寻找时取$\mathcal{N}=1,\; A_1=1.6,\; x_1=2 K(q),\; W_1= K(q)/2$. 图13(b)为第一带隙中的表面扭结孤子, 迭代初值为(9)式中取$a_3=1.4,\; a_4=0.08,\; x_0=T_x/2$. 此时孤子的稳定性仍可用前面方法进行分析, 结论也与前面类似, 故从略.

5.   结　论

数值研究了三体相互作用下均匀介质界面和半无限雅克比椭圆正弦势下BEC中的表面带隙孤子及其稳定性. 在平均场近似下, 其动力学行为可用3次-5次GP方程来描述. 用NCG方法数值求解了各种类型表面孤子, 发现表面亮孤子仅当化学势小于0时才可在带隙内激发产生. 采用线性稳定性分析和非线性动力学演化研究了表面孤子的稳定性, 结果表明三体相互作用会明显影响表面亮孤子的稳定性, 故实验中可通过调节它来改变孤子的稳定性. 另外, 当考虑到三体相互作用时, BEC系统中会激发产生表面扭结孤子和气泡孤子, 它们既可存在于带隙中也可存在于能带中. 数值结果表明表面扭结孤子既有稳定的也有不稳定的, 但表面气泡孤子均不稳定.