曹 蕊，華冬英，王 茜，張讀翠，李祥貴

(北京信息科技大學(xué) 理學(xué)院，北京 100192)

0 引言

玻色—愛因斯坦凝聚態(tài)是一種區(qū)別于人們所熟知的3種物質(zhì)形態(tài)即氣、液、固的新的物質(zhì)狀態(tài)。1924年，玻色提出新的光子的量子統(tǒng)計方法，愛因斯坦將其理論推廣到帶質(zhì)量的理想氣體，從理論上預(yù)言當(dāng)溫度在臨界溫度以下，玻色子(不存在相互作用的粒子)會在最低能量量子態(tài)上凝聚，達(dá)到可觀數(shù)量時便形成凝聚態(tài)，這就是著名的玻色—愛因斯坦凝聚態(tài)(BEC)[1-2]。這些處于BEC的粒子具有宏觀量子相關(guān)、隧穿和量子超流性等特殊的性質(zhì)。隨著激光冷卻等技術(shù)的發(fā)展，關(guān)于玻色—愛因斯坦凝聚態(tài)的實驗取得了突破性的進(jìn)展[3-4]，至此該研究領(lǐng)域飛速發(fā)展，數(shù)值模擬方法也在其相關(guān)的理論和實驗方面發(fā)揮重要作用，其中計算BEC的基態(tài)、第一激發(fā)態(tài)以及動力學(xué)性質(zhì)是BEC研究的基本問題。BEC現(xiàn)象的研究，不但在基礎(chǔ)研究領(lǐng)域有著重要意義，同時在納米技術(shù)、原子芯片技術(shù)、量子計算、原子干涉儀等領(lǐng)域有廣泛的應(yīng)用前景。近幾十年來，關(guān)于BEC基態(tài)的一系列分析和數(shù)值研究，國內(nèi)外學(xué)者們提出了很多不同的數(shù)值算法，其中廣泛采用的是一種虛時法和投影法[5-6]相結(jié)合的方法，該方法使求解基態(tài)問題轉(zhuǎn)為求Gross-Pitaevskii(簡稱G-P)方程最低能量問題，各種迭代格式可由不同的離散方法所得，比如向前/向后歐拉有限差分法，以及間斷伽遼金法、時間分裂偽譜法[7-12]等。

本文針對G-P方程比較和分析兩種數(shù)值計算方法求解玻色—愛因斯坦凝聚態(tài)的基態(tài)。對不同差分格式所得的完全離散格式分別采用完全非線性化迭代方法和線性化處理方法進(jìn)行計算，最后對不同離散格式下兩種算法的數(shù)值結(jié)果進(jìn)行比較和分析。

1 問題介紹

非線性薛定諤方程是量子力學(xué)中的經(jīng)典方程。本文將研究在溫度遠(yuǎn)小于臨界凝聚溫度時描述玻色—愛因斯坦凝聚態(tài)的G-P方程求解基態(tài)解的數(shù)值方法。從模型的角度出發(fā)，考慮以下d(d=1,2,3)維的無量綱G-P方程：

β|ψ(x,t)|2ψ(x,t)x∈d

(1)

初始條件為

ψ(x,0)=ψ0(x)x∈d

(2)

式中：i 為虛數(shù)單位；x為空間坐標(biāo)；t為時間；ψ為一個復(fù)值波函數(shù)；ψ0為給定的初值函數(shù)；Δ為拉普拉斯算子；V(x)≥0為一個實值電勢；β∈為一個任意的無量綱參數(shù)，它表示粒子間相互作用強(qiáng)弱。式(1)中波函數(shù)有兩個重要的不變量，即質(zhì)量 (歸一化即單位化后)

(3)

和粒子能量

E(ψ(·,0))t≥0

(4)

為找到式(1)的穩(wěn)態(tài)解，把穩(wěn)態(tài)解ψ(x,t)表示為

ψ(x,t)=φ(x)e-iμt

(5)

式中：μ為化學(xué)勢；φ為一個與時間無關(guān)的實函數(shù)。將式(5)代入式(1)中，可得：

β|φ(x)|2φ(x)x∈d

(6)

在歸一化條件下：

(7)

可以看出式(6)的本質(zhì)是帶約束條件式(7)的非線性特征值問題。特征值μ可由其對應(yīng)的特征函數(shù)φ(x)表示為

(8)

式中E(φ)為式(4)中給出的與φ相關(guān)的能量函數(shù)。

2 數(shù)值方法

在實際計算中，由于G-P方程(1)的解在遠(yuǎn)場快速衰減，具有齊次Dirichlet邊界條件，通常選擇截斷域足夠大則截斷誤差可忽略不計。為敘述及計算方便，本文只考慮空間一維的情況，則上述式(1)連同初始條件式(2)中x的范圍可以取成Ω=[a,b]。以下討論同理可推廣到二維和三維的情況。目前虛時法是計算G-P方程基態(tài)解的常用方法之一，即用-t替換式(1)中的it，使其變?yōu)槟芰亢纳⒌姆匠?

(9)

為保證歸一化條件式(7)成立，實值波函數(shù)φ在離散的每個時間間隔結(jié)束時，通過投射至單位球面來對該解進(jìn)行歸一化處理，即將時間步長設(shè)為Δt>0，時間方向進(jìn)行剖分tn=nΔt(n=0,1,2,…),在每個時間間隔[tn，tn+1]結(jié)束時讓波函數(shù)對單位球面作投影，即

(10)

φ(x,t)=0x∈?Ω

φ(x,0)=φ0(x)x∈Ω

2.1 空間離散

為方便起見，對一維空間區(qū)域Ω=[a,b]進(jìn)行等距剖分，正整數(shù)J表示剖分總數(shù)，設(shè)置節(jié)點下標(biāo)的集合為：

在空間方向上使用經(jīng)典的二階中心差分對式(9)進(jìn)行離散，可得式(9)的半離散格式為

β|φj(t)|2φj(t)j∈Ωh

(11)

式中Vj∶=V(xj)，初始條件及邊界條件分別為

φ0(t)=φJ(rèn)(t)=0

2.2 時間離散

2.2.1 向后歐拉格式

在時間方向上使用向后歐拉格式進(jìn)行離散得式(9)的完全離散格式：

(12)

初始條件及邊界條件分別為

(13)

此數(shù)值格式離散的能量和化學(xué)勢可表示為：

(14)

則上述離散方程就轉(zhuǎn)為更簡單的線性方程組求解，可直接采用追趕法。這種處理方法簡潔有效，但終究還是回避了原來更精確的格式(12)求解,所以本文仍采用原離散格式(12)進(jìn)行計算。對式中非線性項的處理設(shè)計一種算法——迭代求解法，此方法在理論上更直接，也符合實際理論要求。具體操作如下：

取整數(shù)m≥0表示迭代次數(shù)，離散格式(12)的迭代方法為：

(15)

根據(jù)式(15)，假設(shè)在第n個時間層t=tn時的數(shù)值解已知，求解第n+1個時間層時，將未知變量移到等號左側(cè)，已知變量移到等號右側(cè),化簡得

(16)

2.2.2 Crank-Nicolson格式

(17)

該數(shù)值格式的能量和化學(xué)勢為：

(18)

式(18)可直接用追趕法進(jìn)行求解。

同樣地，式(18)的計算過程回避了式(17)中原本的非線性項，理論上式(17)可使用迭代求解法直接進(jìn)行求解，具體步驟為

(19)

假設(shè)tn層數(shù)值解已知，求解tn+1層,化簡式(19)整理得

(20)

令

則有

令

則AΦn+1,m+1=Kn

式中

3 數(shù)值計算比較與分析

本節(jié)通過數(shù)值算例對上述數(shù)值算法進(jìn)行比較與分析。

表1 向后歐拉格式——兩種算法在不同β下所需總時長t和基態(tài)能量E的對比表

圖1 β=60時，向后歐拉格式—能量隨時間的演化

由表1可見，在此步長取值下，若β取值過大，迭代求解法失效。以下是簡單分析：

由于β很大，式(9)右端前兩項的作用可忽略，僅考慮非線性項作用：

-βφ3(x,t)tn≤t≤tn+1

(21)

計算得

(22)

此時tn+1-tn=τ。當(dāng)β很大，而τ不變，則左端變大，而對于同一時間層上相鄰的φn,m+1和φn,m相差很小，就可能會出現(xiàn)不匹配的問題，從而出現(xiàn)迭代失效。所以當(dāng)β較大的情況下，τ的取值需要相應(yīng)更小，才可保證算法順利進(jìn)行。

由圖1、圖2可以看出在該步長取值下，比起線性化求解法，迭代求解法誤差更大。進(jìn)一步詳細(xì)分析當(dāng)β=60時兩種算法在不同h和τ的取值下的計算結(jié)果，如表2和表3所示。

由表2、表3可以看出，使用向后歐拉格式時，線性化求解法對時間步長和空間步長的取值要求不高，迭代求解法更容易受時間步長大小的影響，需較小的時間步長才可保證迭代順利進(jìn)行。

表2 當(dāng)τ=0.001時向后歐拉格式——兩種算法在不同空間步長下對比表

表3 當(dāng)h=1/128時向后歐拉格式——兩種算法在不同時間步長下對比表

使用C-N格式計算時，取與向后歐拉格式中相同的步長作為精確解，同樣發(fā)現(xiàn)在固定的網(wǎng)格比(τ=0.01，h=1/32)下，若β取值過大則迭代求解法無法求解,需取更小的時間步長來保證該算法求解順利。當(dāng)β=60時，在固定網(wǎng)格下使用C-N格式，兩種算法計算結(jié)果如圖3、圖4所示。同樣出現(xiàn)迭代求解法誤差比起線性求解法更大的現(xiàn)象。為進(jìn)一步詳細(xì)分析，使用C-N格式時，用兩種算法在不同h和τ的取值下的計算結(jié)果如表4和表5所示。

表4 當(dāng)τ=0.001時C-N有限差分格式——兩種算法在不同空間步長下對比表

表5 當(dāng)h=1/128時C-N有限差分格式——兩種算法在不同時間步長下對比表

圖3 C-N格式-能量隨時間演化

圖4 C-N格式 - 基態(tài)解

綜合數(shù)據(jù)可以看出，在相同網(wǎng)格比和時間步長情況下，C-N格式比向后歐拉格式計算時間要短。當(dāng)時間步長變大時，向后歐拉格式能量耗散更穩(wěn)定，而C-N格式會出現(xiàn)不穩(wěn)定情況。以下是簡單分析：

根據(jù)泰勒展開式，可得向后歐拉格式和C-N格式的數(shù)值近似分別為

(23)

(24)

可以看出式(23)中右端第2項是數(shù)值耗散項，式(24)中右端第2項是數(shù)值色散項，數(shù)值耗散會使所求基態(tài)解變得平滑，數(shù)值色散會使所求基態(tài)解振蕩。在例1中，通過實際數(shù)值測試，使用向后差分格式時，τβ的值控制在0～0.07；使用C-N差分格式時，τβ的值控制在0～0.14，可保證在相同β下取不同的τ使用迭代求解法求得的基態(tài)解誤差在10-3內(nèi)。

結(jié)合上述可知，無論是向后歐拉格式還是C-N格式，線性化求解法和迭代求解法這兩種算法均可以用于求解基態(tài)解，兩者計算的能量都隨時間演化呈衰減趨勢。在同樣較小的時間步長下兩種算法的計算結(jié)果較為接近，使用迭代求解法進(jìn)行計算是符合理論要求的。但在β變大時，往往需要取比線性求解法更小的時間步長τ才能使迭代求解法順利進(jìn)行，且所求結(jié)果更為準(zhǔn)確，而使用線性化求解法計算效果相對穩(wěn)定，在使用向后歐拉格式的時候無步長取值要求，總體而言所需總時間步長更少。

4 結(jié)束語

本文利用虛時法加投影法求解G-P方程的基態(tài)解，分別對向后歐拉差分格式及C-N差分格式的兩種數(shù)值計算方法進(jìn)行分析和比較。提出的迭代求解法進(jìn)行數(shù)值計算更符合理論要求，并且結(jié)果更為精確，但此求解過程需要根據(jù)β的大小來取適當(dāng)?shù)臅r間步長τ的值。而線性化求解法穩(wěn)定性較高，求解所需時間更短，實際數(shù)值算例驗證了分析結(jié)果。該討論與分析結(jié)果只針對本文一維G-P方程，其他G-P方程的結(jié)果還需進(jìn)一步的討論與分析。