劉文杰, 王漢權(quán),2

(1.云南財(cái)經(jīng)大學(xué) 統(tǒng)計(jì)與數(shù)學(xué)學(xué)院,昆明 650221;2.云南師范大學(xué) 數(shù)學(xué)學(xué)院,昆明 650504)

0 引 言

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

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

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

1 Bose-Einstein凝聚態(tài)的能量泛函極小值模型簡介

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

(1)

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

(2)

(3)

因此BEC的基態(tài)通常被定義為非凸極小化問題的最優(yōu)值問題[27-29]:

(4)

其中球面約束S被定義為

(5)

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

對于方程(3),令

(6)

存在實(shí)數(shù)μ,使得

(7)

根據(jù)變分法的基本引理

(8)

可得

(9)

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

(10)

(11)

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

(12)

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

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

2.1 一維情形

在齊次Dirichlet邊界條件下,對有界計(jì)算區(qū)域U上截?cái)嗟?/p>

(13)

和

(14)

進(jìn)行離散化,用Lagrange插值多項(xiàng)式來逼近空間導(dǎo)數(shù),用Legendre-Gauss-Lobatto積分公式求解定積分．接下來我們進(jìn)行方程的離散．

首先基于Legendre-Gauss-Lobatto節(jié)點(diǎn),有如下的Lagrange插值多項(xiàng)式形式:

(15)

其中l(wèi)j(s)為Lagrange插值基函數(shù)

(16)

且滿足Lagrange正交性．

(17)

其中

(18)

因此對原函數(shù)基于Legendre-Gauss-Lobatto積分公式進(jìn)行如下的離散操作:

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

同樣地,我們把約束條件也進(jìn)行離散操作:

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

于是得到一個(gè)普通優(yōu)化問題:

(19)

2.2 二維情形

在齊次Dirichlet邊界條件下,對有界計(jì)算區(qū)域U上截?cái)嗟?/p>

和

進(jìn)行離散化,用Lagrange插值多項(xiàng)式來逼近空間導(dǎo)數(shù),用Legendre-Gauss-Lobatto積分公式求解定積分．接下來我們進(jìn)行方程的離散．

首先基于Legendre-Gauss-Lobatto節(jié)點(diǎn),有如下的Lagrange插值多項(xiàng)式形式:

(20)

其中l(wèi)i(s),ln(s)為Lagrange插值基函數(shù)

(21)

(22)

且滿足Lagrange正交性．

(23)

(24)

其中

(25)

(26)

因此對原函數(shù)基于Legendre-Gauss-Lobatto積分公式進(jìn)行如下的離散操作:

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

同樣地,把約束條件也進(jìn)行離散操作:

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

于是得到一個(gè)普通優(yōu)化問題:

(27)

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

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

(28)

其中λ為Lagrange乘子．

在Φ*,λ*處對Lagrange函數(shù)(26)求偏導(dǎo):

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

(29)

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

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

3 數(shù) 值 計(jì) 算

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

3.1 一、二維數(shù)值誤差分析

在文獻(xiàn) [30-31]中已經(jīng)證明了在近似空間中的誤差估計(jì),在本小節(jié)中,我們首先考慮能量泛函極小值問題(4)的一種簡單情形,特別地,當(dāng)β=0時(shí),一維情形(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,…},其中l(wèi)i(x)為Lagrange基函數(shù),假設(shè)φN為離散問題的最小化解,則

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

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

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

表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 一維情形

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

圖2 β=5,10,100,1 000時(shí),一維Bose-Einstein凝聚態(tài)的基態(tài)解φ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,并且我們?nèi)=(-8,8)×(-8,8),然后進(jìn)行數(shù)值模擬．表4為當(dāng)N=30時(shí),Fourier譜方法和Legendre配置譜方法離散后得到的函數(shù)值與β的變化關(guān)系．在表中我們對兩種離散格式進(jìn)行了比較,其中特殊情況β=0時(shí)采取的初始條件為φ0(x,y)=(1/π1/2)e-(x2+y2)/2．通過比較可以發(fā)現(xiàn),Legendre配置譜方法離散后得到的函數(shù)值與Fourier譜方法離散后得到的函數(shù)值誤差也很小．圖3為在N=30的情況下β=5,10,100,1 000時(shí),二維Bose-Einstein凝聚態(tài)的基態(tài)解φ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時(shí)的φg(x,y)Fig.3 The φg(x,y) graphs for β=5,10,100,1 000

在上述數(shù)值算例中我們對Legendre配置譜方法求解的數(shù)值實(shí)驗(yàn)結(jié)果與Fourier譜方法求解的數(shù)值實(shí)驗(yàn)結(jié)果進(jìn)行了比較,發(fā)現(xiàn)Legendre配置譜方法相對而言計(jì)算精度更高,而此方法隨著精度更高的同時(shí)出現(xiàn)了計(jì)算效率不高的情況,但是計(jì)算精度的優(yōu)點(diǎn)遠(yuǎn)遠(yuǎn)超過計(jì)算效率慢的缺點(diǎn),因此針對非旋轉(zhuǎn)的Bose-Einstein凝聚態(tài)的基態(tài)解問題可以使用Legendre配置譜方法來求解．

4 結(jié) 論

本文嘗試?yán)肔egendre配置譜方法求解Bose-Einstein凝聚態(tài)的基態(tài)解．首先系統(tǒng)地介紹了Bose-Einstein凝聚的相關(guān)歷史背景與Bose-Einstein的物理模型,然后通過研究模型的推導(dǎo)過程,將GPE基態(tài)解問題轉(zhuǎn)換成能量泛函極小值問題,對其泛函使用Legendre配置譜方法進(jìn)行離散化,最后進(jìn)行數(shù)值模擬實(shí)驗(yàn)．通過數(shù)值模擬發(fā)現(xiàn)隨著數(shù)值解趨于穩(wěn)定,能量也會(huì)趨于穩(wěn)定變化且變化緩慢．在本文實(shí)驗(yàn)中,我們將Legendre配置譜方法求解的實(shí)驗(yàn)結(jié)果與Fourier譜方法求解的實(shí)驗(yàn)結(jié)果進(jìn)行了比較分析,最后通過實(shí)驗(yàn)分析結(jié)果得出了一個(gè)結(jié)論:針對非旋轉(zhuǎn)的Bose-Einstein凝聚態(tài)的基態(tài)解問題可以使用Legendre配置譜方法來求解,且誤差較?。?/p>

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