Legendre配置譜方法求解Bose-Einstein凝聚態的基態解

劉文杰, 王漢權,2

(1.云南財經大學 統計與數學學院,昆明 650221;2.云南師范大學 數學學院,昆明 650504)

0 引 言

Bose-Einstein凝聚態是Bose氣體冷卻到接近絕對零度時的一種物態,是1920年前后Einstein在Bose分析光子行為的工作基礎上對有質量的粒子所作的預測．20世紀90年代以來,在3位物理學家(Chu(朱棣文)、Cohen、Phillips)的杰出工作下,激光冷卻與囚禁中性原子技術得到了極大發展,也為Bose-Einstein凝聚的實現提供了條件．1995年,第一批實現Bose-Einstein凝聚(BEC)的幾個研究小組分別來自美國科羅拉多大學實驗天體物理聯合研究所(JILA)、美國萊斯大學(Bradley小組)、麻省理工學院(MIT)(Wolfgang Ketterle小組),他們分別獨立宣告在實驗上觀察到了Bose-Einstein凝聚現象,在物理界引起了強烈反響,是Bose-Einstein凝聚研究歷史上的一個重要里程碑．此后,有關BEC的研究迅速發展,觀察到了一系列新的現象,如BEC的相干性、Josephson效應、蝸旋、超冷Fermi原子氣體．

為了研究Bose-Einstein凝聚態的基態解和動力學問題,2005年,Bao等[1]計算了旋轉Bose-Einstein凝聚態的基態、對稱及中心的漩渦態,并研究了其能量和化學勢的變化．2017年,馮悅[2]介紹了在多維勢阱下的一種時空自適應有限元方法求解Bose-Einstein凝聚態的基態解．同年,Liu等[3]提出了梯度法來求解旋轉的雙原子Bose-Einstein凝聚態基態解,并用大量的例子來證明其有效性．2018年,溫建蓉等[4]采用數值方法和Fermi近似來求解非線性Schr?dinger方程,并且研究了Bose子凝聚態的基態穩定性．2021年,Gaidamour等[5]提出了一個HPC譜解器(BEC2HPC),用來求解非線性的Schr?dinger方程和帶旋轉的BEC 基態解問題,該方法主要考慮基于快速Fourier變化的標準偽譜離散化方法,再采用預處理非線性共軛梯度方法來求解歸一化約束條件下能量泛函極小化問題的基態解．2020年,王智軍[6]介紹了正規梯度流和預條件共軛梯度法求旋轉Bose-Einstein凝聚體的基態．同年,Xu等[7]改進了時空自適應有限元方法求解Bose-Einstein凝聚態的基態解．2022年,Chen等[8]引入了兩種二階流作為約束非凸優化問題的能量最小化策略來求解帶旋轉的Bose-Einstein凝聚態的基態問題,并且討論了幾種數值離散方案．而針對非旋轉Bose-Einstein凝聚,1995年,Edwards等[9]提出了一種Runge-Kutta方法來解決方程的基態解．2003年,Bao和Du[10]利用歸一化的梯度流方法來計算Bose-Einstein凝聚態的基態解和第一激發態問題,并探討了離散系統的穩定性．2005年,Bao等[11-12]介紹了時間分裂譜方法,用來計算Bose-Einstein凝聚態的基態解,并驗證了該方法的有效性和準確性．2007年,舒級等[13]在二維空間中討論了一類擬線性Schr?dinger方程,并證明了該方程所對應初值問題的解在一定條件下爆破,同時利用變分方法,也得到了整體解存在的一個充分條件．2013年,Caliari等[14]利用MATLAB 的一個套件程序(GSGPES)來求解GPE系統的基態解問題．2011年,華冬英等[15]將具有徑向對稱的三維Bose-Einstein凝聚態問題簡化為一維的雪茄形問題,并提出了有限元虛數法求解基態解問題．2017年,Wu等[16]使用一種正則化Newton法來求解Bose-Einstein凝聚態的基態解．2018年,楊娜等[17]針對廣義帶導數的非線性Schr?dinger方程的精確解問題進行了研究分析,采用行波變換,將其化為常微分方程動力系統,并計算出該方程動力系統的首次積分,討論了系統在不同參數條件下的奇點與相圖,得到了對應的精確解．2019年,代猛等[18]研究了立方Schr?dinger方程的二階向后差分有限元方法(BDF2-FEM)的無條件最優誤差估計．2021年,曹蕊等[19]用一種迭代求解方法對所得非線性離散方程進行計算,與常規采用的線性化處理方法所得的數值結果進行了詳細比較和分析．結果表明,線性化求解法和迭代求解法這兩種算法均可用于求解基態解,計算所得能量均隨時間演化呈衰減趨勢．而求基態解的數值方法通常分為三類:第一類就是本文所用的方法,第二類是直接離散原來泛函極小值問題所對應的Euler-Lagrange方程,第三類是構造含時間的梯度流方法．三類方法各有優勢,第二類方法的優點是求解跟矩陣相關的非線性特征向量與特征值問題,第三類方法是求解含時間的偏微分方程,而本文這種方法的優點在于可以使用現有的最優化理論與方法[20]中介紹的constrained minimization方法來求解帶約束的優化問題．基于此,本文對Bose-Einstein凝聚態的基態解問題做了3個方面的研究．首先,對Bose-Einstein凝聚態的Gross-Pitaevskii方程(GPE)進行降維和無量綱化處理,將GPE問題轉換成能量泛函極小值問題．其次,嘗試通過Legendre配置譜方法[21]的離散方法對能量泛函極小值問題進行離散．最后,進行數值模擬實驗,并對實驗結果進行分析,得出結論．

本文的結構安排如下:第1節對Bose-Einstein凝聚態的能量泛函極小值模型進行了簡單的介紹; 第2節介紹了Legendre配置譜方法一維和二維的具體離散格式;第3節通過給出的數值例子對該問題進行數值模擬并進行了分析; 第4節對本文的工作做了一個簡單的總結．

1 Bose-Einstein凝聚態的能量泛函極小值模型簡介

自從稀Bose原子氣體中首次實驗的實現,BEC引起了原子、分子和光學(AMO)物理界和凝聚態物質界的極大興趣．在描述三維(3D)的GPE時[22-25],可以用非線性Schr?dinger方程(NLSE)或宏觀波函數ψ=ψ(x,t)來描述絕對零度和低溫狀態下的性質．然后通過使用降維和無量綱化[1,26]的手法對原問題進行適當降維并得到無量綱GP方程(對于非旋轉的BEC,即Ω=0時d=1,2,3):

(1)

其中β∈為無量綱化相互作用系數,V(x)為無量綱化實值外部捕獲勢．且波函數歸一化和能量泛函分別表示為

(2)

(3)

因此BEC的基態通常被定義為非凸極小化問題的最優值問題[27-29]:

(4)

其中球面約束S被定義為

(5)

我們能夠驗證問題(4)的變分形式是一個非線性特征值問題．

對于方程(3),令

(6)

存在實數μ,使得

(7)

根據變分法的基本引理

(8)

可得

(9)

即得到在限制條件下求解特征值問題:

(10)

(11)

這是一個正規化限制下的非線性特征值問題,任何特征值μ都可以用與之對應的φ(x)通過下式得到:

(12)

事實上,求得的特征函數是單位球面上能量泛函的臨界點．為了找到Bose-Einstein凝聚態的基態解,在單位球面上最小化能量泛函,即求φg∈S,使得滿足式(4)和(5)．

2 能量泛函極小值問題的離散化

2.1 一維情形

在齊次Dirichlet邊界條件下,對有界計算區域U上截斷的

(13)

和

(14)

進行離散化,用Lagrange插值多項式來逼近空間導數,用Legendre-Gauss-Lobatto積分公式求解定積分．接下來我們進行方程的離散．

首先基于Legendre-Gauss-Lobatto節點,有如下的Lagrange插值多項式形式:

(15)

其中lj(s)為Lagrange插值基函數

(16)

且滿足Lagrange正交性．

(17)

其中

(18)

因此對原函數基于Legendre-Gauss-Lobatto積分公式進行如下的離散操作:

其中Φ=(φ(x1),φ(x2),…,φ(xN-1))．

同樣地,我們把約束條件也進行離散操作:

其中Φ=(φ(x1),φ(x2),…,φ(xN-1))．

于是得到一個普通優化問題:

(19)

2.2 二維情形

在齊次Dirichlet邊界條件下,對有界計算區域U上截斷的

和

進行離散化,用Lagrange插值多項式來逼近空間導數,用Legendre-Gauss-Lobatto積分公式求解定積分．接下來我們進行方程的離散．

首先基于Legendre-Gauss-Lobatto節點,有如下的Lagrange插值多項式形式:

(20)

其中li(s),ln(s)為Lagrange插值基函數

(21)

(22)

且滿足Lagrange正交性．

(23)

(24)

其中

(25)

(26)

因此對原函數基于Legendre-Gauss-Lobatto積分公式進行如下的離散操作:

其中Φ是一個(N-1)×(N-1)的矩陣．

同樣地,把約束條件也進行離散操作:

其中Φ是一個(N-1)×(N-1)的矩陣．

于是得到一個普通優化問題:

(27)

對于優化問題(17)和(25),我們對解的存在性進行了討論．首先針對約束優化問題考慮使用Lagrange函數:

L(Φ,λ)=f(Φ)-λ(g(Φ)-1),

(28)

其中λ為Lagrange乘子．

在Φ*,λ*處對Lagrange函數(26)求偏導:

?ΦL(Φ*,λ*)=?Φf(Φ*)-λ*?Φ(g(Φ*)-1)．

(29)

因此對于方程(27),如果對于任意解Φ*,λ*使得?ΦL(Φ*,λ*)≠0,則優化問題的解不存在．反之,如果對于任意解Φ*,λ*使得?ΦL(Φ*,λ*)=0,則此優化問題解存在,并且滿足最優性條件(KKT):

針對以上一維、二維的離散優化問題有解的情況,考慮對優化問題(17)和(25)使用現有的最優化理論與文獻[20]中介紹的內點法來求解．

3 數 值 計 算

前面我們已經對一維和二維的能量泛函方程進行了離散,下面對具體的計算例子進行數值計算和分析．

3.1 一、二維數值誤差分析

在文獻 [30-31]中已經證明了在近似空間中的誤差估計,在本小節中,我們首先考慮能量泛函極小值問題(4)的一種簡單情形,特別地,當β=0時,一維情形(d=1),能量泛函極小值問題(4)有真解,為φg(x)=(1/π1/4)e-x2/2．二維情形(d=2),能量泛函極小值問題(4)有真解,為φg(x)=(1/π1/2)e-(x2+y2)/2,并且考慮Legendre配置譜方法的近似和空間XN=span{li(x),i=1,2,…},其中li(x)為Lagrange基函數,假設φN為離散問題的最小化解,則

在XN空間下,有如下形式:

又由式(16)知,當N→∞時,minvN∈XN‖vN-φ‖H1→0,故當N→∞時,‖φN-φg‖L2→0,即是收斂的．

因此本文對離散誤差ε=‖φN-φg‖L2進行了數值模擬．表1和表2為一維和二維的誤差值與N的變化關系,從表中數據可以看出當N→∞時,空間分割較細,φN→φg,則‖φN-φg‖L2→0,因此是具有收斂性的．圖1為一維(-16,16)和二維(-8,8)×(-8,8)的誤差收斂性圖像,從圖中可以發現當N≥30時,誤差值很小,收斂率快,精確度很高．此外Legendre配置譜方法具有收斂性強,收斂速度快等特點．

表1 一維情況下,改變N的大小,誤差ε的變化情況Table 1 In the 1D case,changes of error ε with N

表2 二維情況下,改變N的大小,誤差ε的變化情況Table 2 In the 2D case,changes of error ε with N

圖1 一維(左圖)、二維(右圖)情況下,改變N的大小,誤差ε的變化情況Fig.1 Changes of error ε with N in the 1D case (left) and the 2D case (right)

3.2 一維情形

針對具有強相互作用的非旋轉BEC,即β?1,初始解通常選擇Thomas-Fermi近似:

圖2 β=5,10,100,1 000時,一維Bose-Einstein凝聚態的基態解φg(x)Fig.2 Ground state solution φg(x) of the 1D Bose-Einstein condensates for β=5,10,100,1 000

表3 一維情況下Legendre配置譜方法的能量值Eβ(φg)和Fourier譜方法的能量值與β之間的變化情況Table 3 Energy value Eβ(φg) of the Legendre collocation spectrum method and energy value of the Fourier spectrum method in the 1D case,changing with β

3.3 二維情形

對于二維方程,初始解同樣選擇Thomas-Fermi近似:

其中μ=(β/π)1/2,并且我們取U=(-8,8)×(-8,8),然后進行數值模擬．表4為當N=30時,Fourier譜方法和Legendre配置譜方法離散后得到的函數值與β的變化關系．在表中我們對兩種離散格式進行了比較,其中特殊情況β=0時采取的初始條件為φ0(x,y)=(1/π1/2)e-(x2+y2)/2．通過比較可以發現,Legendre配置譜方法離散后得到的函數值與Fourier譜方法離散后得到的函數值誤差也很小．圖3為在N=30的情況下β=5,10,100,1 000時,二維Bose-Einstein凝聚態的基態解φg(x,y)．

表4 二維情況下Legendre配置譜方法的能量值Eβ(φg)和Fourier譜方法的能量值與β之間的變化情況Table 4 Energy value Eβ(φg) of the Legendre collocation spectrum method and energy value of the Fourier spectrum method in the 2D case,changing with β

(a) β=5 (b) β=10

(c) β=100 (d) β=1 000圖3 β=5,10,100,1 000時的φg(x,y)Fig.3 The φg(x,y) graphs for β=5,10,100,1 000

在上述數值算例中我們對Legendre配置譜方法求解的數值實驗結果與Fourier譜方法求解的數值實驗結果進行了比較,發現Legendre配置譜方法相對而言計算精度更高,而此方法隨著精度更高的同時出現了計算效率不高的情況,但是計算精度的優點遠遠超過計算效率慢的缺點,因此針對非旋轉的Bose-Einstein凝聚態的基態解問題可以使用Legendre配置譜方法來求解．

4 結 論

本文嘗試利用Legendre配置譜方法求解Bose-Einstein凝聚態的基態解．首先系統地介紹了Bose-Einstein凝聚的相關歷史背景與Bose-Einstein的物理模型,然后通過研究模型的推導過程,將GPE基態解問題轉換成能量泛函極小值問題,對其泛函使用Legendre配置譜方法進行離散化,最后進行數值模擬實驗．通過數值模擬發現隨著數值解趨于穩定,能量也會趨于穩定變化且變化緩慢．在本文實驗中,我們將Legendre配置譜方法求解的實驗結果與Fourier譜方法求解的實驗結果進行了比較分析,最后通過實驗分析結果得出了一個結論:針對非旋轉的Bose-Einstein凝聚態的基態解問題可以使用Legendre配置譜方法來求解,且誤差較小．

本文在利用Legendre配置譜方法求解Bose-Einstein凝聚態的基態解問題時也遇到了一些困難,例如,數值計算效率慢等．因此,接下來的工作我們會繼續對Legendre配置譜方法求解Bose-Einstein凝聚態的基態解問題進行相關的研究,嘗試找出一種既能提高效率又能增加計算精度的方法．此外我們也會關注更高維的情況以及帶旋轉型Bose-Einstein凝聚態的基態解問題．