三體相互作用下準一維玻色-愛因斯坦凝聚體中的帶隙孤子及其穩定性*

唐娜 楊雪瀅 宋琳 張娟3) 李曉霖 周志坤 石玉仁?

1) (西北師范大學物理與電子工程學院, 蘭州 730070)

2) (甘肅省原子分子物理與功能材料重點實驗室, 蘭州 730070)

3) (蘭州工業學院基礎學科部, 蘭州 730050)

具有三體相互作用的玻色-愛因斯坦凝聚體(Bose-Einstein Condensate, BEC)束縛于雅可比橢圓周期勢中, 在平均場近似下可用3—5次Gross-Pitaevskii方程(GPE)描述.首先利用多重尺度法對該系統進行了理論分析, 將 GPE 化為一定態非線性薛定諤方程 (Nonlinear Schr?dinger Equation, NLSE), 并給出了一類帶隙孤子的解析表達式.然后采用牛頓共軛梯度法數值得到了該系統中存在的兩類帶隙孤子, 發現孤子的振幅隨著三體相互作用的增強而減小, 這與多重尺度法分析所得結論一致.最后用時間劈裂傅里葉譜方法對GPE進行長時間動力學演化以考察孤子的穩定性, 發現系統中既存在穩定的帶隙孤子, 也存在不穩定的帶隙孤子,且外勢的模數會對孤子的結構和穩定性產生明顯影響.

專題：非線性物理

1 引 言

玻色-愛因斯坦凝聚(Bose-Einstein condensate,BEC)現象是在宏觀尺度上能觀察到的最顯著的多體量子現象之一, 最早是由Bose和Einstein在1924年提出, 即理想的玻色子在非常低的溫度下,大部分粒子會突然跌落到最低能級上, 處于這種新狀態的物質被稱為 BEC[1,2].隨著在87Rb 和23Na 等一系列堿金屬原子氣體實驗中實現了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 三體相互作用BEC中的帶隙孤子

2.1 理論模型

具有三體相互作用的BEC束縛于外勢中, 在平均場近似下系統的動力學行為可用如下GPE描述[35]

其中Ψ=Ψ(r,t) 是系統的序參量 (波函數),r=(x,y,z)是位矢,g=4π2as/m表征原子間相互作用強度,m是原子質量.as是s波散射長度(as>0表示原子間相互排斥,as<0 表示原子間相互吸引), 可由 Feshbach共振調節[36].h表征三體相互作 用 強 度 ,V(r) 表 示 外 勢 , 總 粒 子 數實驗中, BEC 通常束縛于一諧振子勢阱內, 其中ωx,ωy,ωz分別表示x,y,z方向的頻率.當ωx≈ωy且ωz?ωx時, BEC在z方向被“凍結”于基態, 這時系 統 可 用 準 二 維 GPE描 述; 當ωy,ωz?ωx時,BEC 在y,z方向均被“凍結”, 將在x方向呈雪茄狀分布, 系統的動力學行為可用準一維GPE描述[37].

考慮BEC被囚禁于光晶格勢下的準一維情形, 并做如下無量綱化處理:

其中n0為無量綱化時粒子數密度的特征量, 此時總粒子數為為方便起見, 這里已省略變量上面的“~”.

取 外 勢V(x)=V0sn2(x,q) , 其 中 s n(x,q) 為 雅可比橢圓正弦函數,q是其模數 ( 0q1 ).當q=0時,V(x)=V0sin2x, 文獻 [38]中考慮的外勢即為此情形, 所以該勢可看作是對相關文獻的推廣.當q=1 時,V(x)=V0tanh2x, 此時V(x) 并非周期勢, 故后面不考慮此情形.若考慮該勢的實驗實現, 可用下面公式[39]

從中可以看出, 該勢在實驗中可用兩束不同頻率的激光來近似實現.當q<0.9 時, 其近似程度可達99%.

尋找方程(2)下列形式的定態解

其中ψ(x) 為實函數,μ為化學勢.代入方程 (2)后得

若波函數振幅為無窮小, 則可忽略非線性項, 得

方程(7)為一廣義馬丟(Mathieu)方程, 它的有界解被稱為布洛赫模(Bloch Modes), 對應的化學勢μ構成了布洛赫帶 (Bloch Bands).一般情況下, 它的有界解可寫為

將(8)式代入方程(7)后, 通過求解所得本征問題便可得色散關系μ=μ(k) , 從而得到系統的能帶結構.這一點在文獻[29]中有詳細討論, 此處不再贅述.在 Bands中, 線性波可以傳播; 在 Gaps 中, 雖然線性波無法傳播, 但可存在具有局域結構的非線性波, 即帶隙孤子 (gap soliton).

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

若波函數振幅不是無窮小, 則方程(6)中非線性項不能忽略.考慮化學勢μ從 Band的邊界k=k0,μ0=μ(k0) 處進入帶隙, 且k?k0,μ?μ0均為小量的情形, 可用多重尺度法對其進行理論分析.引入多重尺度X0=x,X1=εx, 將化學勢與波函數展開為

其中ε=k?k0是一小量.將此展開式代入方程(6), 并假設η=ε?2=O(ε?2) .按e的同次冪項合并后得

方程(11)與(7)形式相同, 故具有下列形式的解

其中B(X1) 表示慢變包絡 (調制波),p(X0) 為快變載波, 顯然該解滿足Fredholm條件[40], 這里函數f1(x)與f2(x) 的內積定義為

“*”表示復共軛.方程 (12)的解可寫為

則有L0H=p′(X0) , 此時Fredholm條件也自動成立.

將ψ0和ψ1代入方程 (13) 得

再次應用Fredholm條件, 得

在X1=0 處B(X1) 取得極值, 故調制波的振幅由此可知在其他參數給定時, 隨著非線性相互作用g,η的增大, 帶隙孤子的振幅會單調遞減, 后面的數值結果也證實了該結論.另外, 為保證b為實數, 則μ2與D須異號,即 s gn(μ2)= ?sgn(D) , 這里 s gn(x) 為符號函數.該式表示此時的μ必位于帶隙內, 即孤子確為帶隙孤子.為保證解不存在奇性, 須d為正實數且δ+γ>0, 則有舉例來說, 如果要研究半無界帶隙內的孤子, 則此時D>0,μ2<0 , 若取g,η<0則兩個條件必同時滿足(充分非必要); 但若研究第一帶隙內接近第一Band的孤子, 則μ2>0,D<0, 取g,η>0 則兩個條件必同時滿足(充分非必要).以上分析對于數值研究帶隙孤子有一定的指導意義.

2.3 牛頓共軛梯度法尋找帶隙孤子

牛頓共軛梯度法(Newton-Conjugate-Gradient,NCG)[41]是一種高效的數值方法, 可用來求解非線性演化方程的孤立波解.該方法的主要思路是用牛頓迭代法結合共軛梯度法求解所得線性方程, 其收斂速度比共軛梯度法和牛頓法等其他現有的迭代法要快, 而且容易編程實現.文獻[42]也提出了牛頓法與共軛梯度法的組合方法, 并證明該方法的全局收斂性.NCG法是解決無約束最小優化問題的方法之一[43].文獻[44]中運用該方法討論了路徑約束動力學演化過程的優化問題.NCG法在一些物理模型上得到了廣泛應用, 例如具有周期外勢或無周期外勢的二維非線性薛定諤方程, Kortewegde Vries (KdV)方程和五階 Kadomtsev Petviashvili(KP)方程的求解等[41].它可用來尋找各種物理系統(如非線性光學, BEC和水波等)中的孤立波,而且既能尋找基態也能求解激發態, 可以作為處理這類問題的首選方法.

下面采用NCG尋找方程(6)的帶隙孤子.計算時, 需要將無窮區間x∈(?∞,+∞) 截斷為有限區間.通常可取外勢的足夠多個周期作為計算的空間范圍, 然后對該區間進行離散化后便可應用NCG求解.計算表明, 所得結果對迭代初值有一定的依賴性(但并不十分敏感).若迭代初值選擇不當, 則迭代過程會發散或者收斂于平凡解ψ=0 .計算時, 取以下多個高斯波包的疊加

作為迭代初始條件, 這里Aj,xj,Wj分別表示開始迭代時第j個波包的振幅、中心位置和寬度,N表示總波包數.通過選擇合適的Aj,xj,Wj,N, 便可得到所需帶隙孤子.若欲尋找結構更為復雜的帶隙孤子, 可嘗試選用其他形式的迭代初值.

圖1顯示了q=0.1 和0.99時不同參數條件下的單峰帶隙孤子, 迭代時取N=1,A1= 0.6,x1=0,W1=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|Ψ|為孤子振幅.圖2(a)和圖2(b)分別顯示了在吸引相互作用下半無界帶隙內和排斥相互作用下第一帶隙內單峰孤子振幅隨相互作用強度的變化.可以看出, 隨著兩體相互作用強度|g|和三體相互作用強度|h|的增大, 孤子的振幅確在單調遞減.

圖3顯示了不同參數情形下圖1中on-site孤子三體和兩體相互作用能量的比值, 其中可以看出, 當h固定時, |g|越大, 比值越小, 表明二體相互作用越占優; 但當g固定時, 隨著|h|的增大,比值也在增大, 表明三體相互作用逐漸增強.圖3(a)中兩種能量已處于可比擬的范圍; 圖3(b)中h足夠大時, 三體相互作用能量已超過兩體相互作用能量.這種情況下, 三體相互作用更不可忽略.

圖1中所示的帶隙孤子是結構最為簡單的一類基本孤子.除此之外, 也存在另一類如圖4所示的雙峰孤子.這類孤子具有兩個波峰, 兩波峰間的中心位置在外勢的最高點處, 一般稱其為off-site soliton.圖4(a)和圖(b)所示兩波峰具有相同的相位, 稱為同相偶極孤子 (迭代時取N=2 ,A1=A2=0.6,x1= 0,x2= –2K(q),W1=W2=K(q));圖4(c)和圖4(d)所示孤子則被稱為異相偶極孤子 (迭代時取N=2 ,A1= 0.6,A2= –0.6,x1= 0,x2= –2K(q),W1=W2=K(q)).該類孤子的振幅也隨著|g|和|h|的增大而減小.從波形結構上來看,off-site孤子似乎可以看作是兩個單峰on-site孤子的組合; 但從動力學穩定性來看, on-site孤子始終穩定而同相偶極孤子始終不穩定, 故同相偶極孤子也被視為另一類基本孤子.異相偶極孤子的穩定性則較為復雜, 一般依賴于系統的參數.

圖1 3?5 次 GPE 的單峰帶隙孤子 ( V0=4 ) (a), (b) q =0.1 ; (c), (d) q =0.99 .陰影部分表示外勢 V (x) 低處Fig.1.Profiles of singl-hump gap solitons of the cubic-quintic GPE ( V0=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 單峰帶隙孤子的振幅隨相互作用強度的變化( V0=4 )Fig.2.Amplitudes of single-hump gap solitons v.s.nonlinear interaction strength.

圖3 on-site孤子三體相互作用能量和兩體相互作用能量的比Fig.3.The ratio of three-body energy to two-body energy for on-site solitons with different interaction strength.

圖4 3?5 次 GPE 的同相偶極孤子和異相偶極孤子 ( V0=4 ) (a), (b) q =0.1; (c), (d) q =0.99 .陰影部分表示外勢 V (x) 低處Fig.4.Profiles of double-hump gap solitons of the cubic-quintic GPE ( V0=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) .

3 帶隙孤子的穩定性

孤子的穩定性無論在理論上還是在實驗上都是一個重要的問題, 下面將用非線性動力學演化的方法來研究前面所得帶隙孤子的穩定性.在初始時刻, 給波函數一小擾動.若經過足夠長時間演化后孤子的振幅和波形沒有發生明顯變化, 則可認為該孤子動力學穩定; 否則認為它動力學不穩定.數值計算時, 所用方法為時間劈裂傅里葉譜方法[45].該方法具有效率高、精度高、計算過程中粒子數守恒,而且程序容易實現等優點.

對初始波函數做下面類型的擾動

其中ψ(x) 是前面所得定態帶隙孤子,β?1 是擾動參數, 計算時取β=0.01 .

圖5為不同參數情形下不同類型帶隙孤子動力學演化的等值線圖(V0=4 ).可以看出, 不同類型孤子的穩定性有所不同.圖5(a)為μ=1 ,q=0.1,g= –1,η=?1 (對應圖1(a)中情形)時 onsite孤子的動力學演化.可以看出, 經過一段時間演化后波形和振幅都沒有發生明顯變化, 故它是動力學穩定的.數值研究中, 通過對各種參數情形進行計算, 結果均表明無論在第一帶隙還是半無界帶隙中, on-site 孤子始終穩定.圖5(b)為μ=0.5 ,q=0.1,g= –1,η=?1 (對應圖4(a)中同相偶極孤子)的動力學演化.很明顯孤子演變成振蕩狀態,即能量在兩個相鄰光晶格之間周期性地轉移.初始時刻孤子空間結構的對稱性隨時間演化時被破壞,所以這是一種振蕩型不穩定.大量計算表明, 這種同相偶極孤子總是動力學不穩定的.即使q較大時, 盡管兩波峰相距較遠, 它仍然表現出不穩定(但不穩定性變弱).所以這種類型的孤子不宜視為兩個 on-site 孤子的組合.圖5(c)為μ=2 ,q= 0.5,g= 1,η=1 時第一帶隙內異相偶極孤子 (帶隙孤子結構與圖4(c)和圖4(d)中類似)的動力學演化.可以看出, 此時的孤子仍具有振蕩不穩定性.圖5(d)為q=0.99 (其他參數與圖5(c)中相同)時異相偶極孤子的動力學演化.可以看出, 此時的帶隙孤子是動力學穩定的.這一結果表明外勢的模數會影響帶隙孤子的穩定性; 同時也意味著其他參數固定時, 存在臨界值qc.當q?qc時, 異相偶極孤子穩定; 而當q

圖5 不同類型帶隙孤子的動力學演化 (a) on-site 孤子; (b) 同相偶極孤子; (c), (d) 異相偶極孤子Fig.5.Contour plots of | Ψ(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,η=?2,V0= 4 不同μ值時的亞基本帶隙孤子, 它有著和前面兩種類型的孤子明顯不同的結構特征.這類孤子在空間上呈奇對稱分布, 它的中心位置位于外勢的最低點處, 波峰則介于外勢最高點與最低點之間.圖6(b)和圖6(d)為該類帶隙孤子的動力學演化, 可以看出,μ=2.9 時孤子不穩定而μ=3 時孤子穩定.該結果意味著存在臨界值μc,當μ<μc時孤子不穩定而μ>μc時孤子穩定.數值結 果 表 明 在 該 組 參 數 下,μc≈ 2.992 .這 里 取ψ(x)=A1e?x2/W1sin(k1x)作為迭代初始條件, 計算時取A1= 1,W1= 4K(q),k1= 1.

三體相互作用強度對于帶隙孤子的穩定性也有一定影響.圖7 顯示了V0= 2,μ= 1.25,q= 0.1,g= 1而不同h時第一帶隙中同相偶極孤子的動力學演化.圖7(a)中η=0 , 可以看出, 經過一段時間演化后帶隙孤子呈現出動力學不穩定; 圖7(b)中η=?0.2, 可以看出此時的帶隙孤子是動力學穩定的.表明三體相互作用強度對于帶隙孤子的穩定性確實有著一定影響, 它可以改變帶隙孤子的穩定性區域, 該結論與文獻[18]中結論一致.注意這里的同相偶極孤子與圖3中有所不同, 它屬于文獻[29]中提到的(1, 1)結構, 并非從Band中分岔出來(迭代時取N=2 ,A1=A2= 0.6,x1= 0,x2=–2K(q),W1=W2=K(q) ).

圖6 (a), (c) 第一帶隙中的亞基本帶隙孤子 (紅色).藍線表示外勢; (b), (d) 亞基本帶隙孤子的動力學演化Fig.6.(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 | Ψ(x,t)| for perturbed gap solitons.

圖7 第一帶隙中同相偶極孤子 (a) η =0 和 (b) η =?0.2 時的動力學演化Fig.7.Contour plots of | Ψ(x,t)| for in-phase dipole solitons lie in the first bandgap with (a) η =0 and (b) η =?0.2 .

4 結 論

本文研究了準一維情形下束縛于雅可比橢圓函數周期勢中具有三體相互作用BEC系統中的帶隙孤子及其穩定性.在平均場近似下, 系統的動力學行為可由3—5次GPE描述.首先利用多重尺度法對系統的帶隙孤子進行了理論分析, 將3—5次GPE化為一定態NLSE, 并給出了一類帶隙孤子的解析表達式.然后采用NCG法數值得到了該系統中存在的帶隙孤子, 包括兩類基本帶隙孤子(on-site孤子和off-site孤子)與亞基本帶隙孤子 (sub-fundamental gap soliton).數值結果與理論分析均表明, 三體相互作用強度的增大將會導致帶隙孤子的振幅減小.最后, 用時間劈裂傅里葉譜方法對GPE進行長時間動力學演化以考察孤子的穩定性, 發現on-site孤子始終穩定, 同相偶極孤子和異相偶極孤子既有穩定的也有不穩定的.三體相互作用強度對于帶隙孤子的穩定性也有一定影響.外勢的模數對孤子的結構和穩定性會產生明顯影響, 故實驗中可通過調整外勢的模數來改變帶隙孤子的穩定性.