高階SF-SFDTD 方法在含時薛定諤方程求解中的應用研究*

謝國大 潘攀 任信鋼 馮乃星 方明 李迎松 黃志祥?

1) (安徽大學電子信息工程學院,合肥 230601)

2) (安徽大學,智能計算與信號處理教育部重點實驗室,合肥 230601)

1 引言

計算電磁學在當今各種電子設備中電磁效應(輻射、散射或傳播)的建模和仿真中占據重要地位[1,2].其應用包括但不限于通信、生物成像和診斷、電磁兼容性或干擾、信號完整性分析等領域[3,4].近年來,由于現代電子技術的飛速發展,集成電路相關元器件的特征尺寸已經縮小到納米級別,這大大提高了單位體積內物質存儲和信息處理的效率,為智能化、小型化電子設備的發展提供了強有力的支持[5,6].又因為納米器件本身尺寸極小,結構中粒子的量子效應明顯,這導致基于經典電磁理論的建模方案難以描述微觀粒子的量子特性.另外,采用實驗方法直接測量器件的各種性能又十分困難,且成本較高.因此,如何準確、高效地研究納米器件的相關特性,已成為當今納米器件建模和設計的關鍵挑戰.

對于分析大部分納米器件的特性,確定器件結構電子的本征值和本征態是一個重要的切入點.例如,在研究彈道輸運時,需要考慮傳導通道的橫向本征態,這一點與非平衡格林函數密切相關[7].此外,微觀電子傳輸中本征態的激發或不同本征態之間的相互作用,會產生各種有趣的現象,如Fano共振[8]、共振隧穿效應[9]等.因此,準確而有效地計算納米結構中的本征態和特征頻率對理解微納結構系統中各種微觀現象的物理機制和優化器件性能起著重要作用.而薛定諤方程可以描述物質粒子的運動和特性,通過數值求解薛定諤方程,可以確定納米器件的本征值與本征態,這就使得薛定諤方程的數值解變得越來越重要[7].在時間和空間上具有二階數值精度的時域有限差分(finite-difference time-domain,FDTD(2,2))方法被證明是求解含時薛定諤方程的最有效和最簡單的方法之一[10,11].這是因為FDTD(2,2)方法無需對薛定諤方程作進一步的推導,只需要直接對控制方程中的時間和空間偏導項進行中心差分離散,最后導出數值迭代公式.然而,FDTD(2,2)方法在時間和空間上只具有二階精度,在長時間的數值迭代中產生較大的誤差累計,導致仿真結果不正確或與理論結果差距較大.另外,FDTD(2,2)方法中時間步長的取值范圍受時間穩定性條件的影響,即最大時間步長受計算區域中網格尺寸大小的限制,這極大地限制了采用FDTD(2,2)方法求解薛定諤方程的計算效率.

為了克服傳統FDTD(2,2)方法數值色散大的問題,一些高階FDTD(2,4)方法被應用于求解含時薛定諤方程[12,13],數值計算誤差有明顯改善.然而高階FDTD(2,4)方法對控制方程的離散是一種非辛離散形式,在長時間的數值仿真中難以保證系統的能量守恒特性.由于薛定諤方程的時間演化,本質上是歸一化波函數的多次旋轉.文獻[14]提出了采用高階辛FDTD(3,4) (symplectic finitedifference time-domain,SFDTD(3,4))方法求解含時薛定諤方程,所推導出來的辛算子具有對稱屬性,確保離散方程在時間演化過程中保持能量守恒,沒有幅度誤差.然而,高階辛FDTD(3,4)方法也是一種顯式迭代算法,雖然具有比傳統FDTD方法更寬松的Courant-Friedrichs-Lewy (CFL)穩定性條件[15],但計算效率有待進一步的提高.為了突破CFL 穩定性條件的限制,一些隱式無條件穩定的FDTD 方法被開發出來并用于薛定諤方程的數值求解,但隱式算法往往需要對控制方程進行特殊的離散處理及繁瑣的公式推導和數值計算: 包括對時間偏導數的隱式離散、時間步長分裂處理及矩陣求逆計算等.另外,采用較大的時間步長會導致隱式FDTD 方法的數值色差較大,在實際應用中,需要根據實際的計算需求選擇合適的時間步長來平衡計算精度和計算效率.

鑒于傳統FDTD(2,2)方法及其優化算法在計算效率和計算精度兩方面所存在的問題,本文提出了一種穩定性條件可控的空間濾波(spatial filtering)-SFDTD(3,4)(SF-SFDTD(3,4))方法.在SFDTD(3,4)方法中采用不滿足穩定性條件的時間步長會導致空間頻域中高頻電磁分量幅值呈現指數倍增長并導致數值結果發散,因此,在數值迭代過程中采用空間濾波方法[16–18]能夠有效地濾除這些不穩定的高頻分量,保證計算結果的穩定性.又因為當時間步長的取值滿足CFL 穩定性條件時,高頻分量的幅值幾乎零,因此對計算結果的精確性影響較小.采用不滿足CFL 穩定性條件的時間步長導致高頻分量的幅值增加,但這一結果只會影響計算結果的穩定性,濾除這些不穩定的高頻分量幾乎不會對計算結果的精度造成影響.另外,SF-SFDTD(3,4)方法與傳統SFDTD(3,4)方法采用相同的迭代公式,只需要在每一次的數值迭代過程中加入濾波操作.因此,它保留了傳統SFDTD(3,4)方法的簡潔性和有效性,同時提高了計算效率.

論文的組織結構如下: 第2 節介紹了基于含時薛定諤方程的SF-SFDTD(3,4)方法的基本原理,包括擴展SFDTD(3,4)方法的CFL 條件限制,歸一化低通濾波器半徑的定義及濾波處理過程,第3節給出了SF-SFDTD(3,4)方法的數值色散分析,第4 節通過幾種典型的數值算例來驗證SF-SFDTD(3,4)方法的正確性和有效性.最后給出本文工作的總結與展望.

2 方 法

2.1 SFDTD 方法求解薛定諤方程的迭代公式

具有靜態電位V(r,t)=V(r)的含時薛定諤方程表達式為[9]

其中ψ是與空間r和時間t有關的波函數,m*是粒子的質量是哈密頓算符.為避免在數值計算中使用復數,將ψ(r,t) 分解為實分量和虛分量,即:

將(2)式代入(1)式中,得到以下方程組:

根據文獻[14],采用3 步3 階的顯式辛算法,即SFDTD(3,4)方法數值求解含時薛定諤方程,其實部和虛部波函數的迭代公式可離散成如下形式:

其中cl,dl為辛傳播因子[14];n表示時間步,m代表單個時間迭代中子迭代步數,l表示在m個子迭代步數中的第l步計算,i,j,k表示坐標.

2.2 SFDTD(3,4)方法穩定性條件的擴展

不失一般性,假定V(r)=0,等式(3)的空間四階差分形式可以寫成如下形式:

這里kx,ky,kz表示相應方向的波數.

對(12)式進行3 階時間辛積分近似處理,則其系統時間演化矩陣的離散形式S為

由于離散矩陣是行列式值為1 的辛矩陣,矩陣特征值λ 對應的特征方程可以寫為

(15)式的解為

穩定性條件要求|λ1,2|=1,矩陣S的跡滿足如下條件:

將(14)式依次相乘,得到

其中gl的值可以通過辛傳播因子cl和dl來計算,m是單個時間迭代中子迭代步數.

將其代入(17)式和(18)式,則關于時間步長Δt的不等式需滿足如下形式:

為確保數值穩定性,時間步長Δt應小于(21)式右邊函數的最小值,令函數Q中sin2(ku?u)=0(u=x,y,z),可得傳統SFDTD(3,4)算法的CFL 穩定性條件:

(22)式是確保SFDTD(3,4)方法中所有空間頻域場分量保持穩定,然而根據文獻[16,17],只有空間頻域中一部分的低頻分量對計算結果有影響,而剩余部分的高頻分量的幅值較小,同時對數值計算結果的正確性影響幾乎為零.然而,若采用較大的時間步長,會導致空間頻域中高頻分量的幅值呈現指數增加,并導致計算結果發散.基于此,將SFDTD(3,4)方法中空間頻域場分量限制在一個范圍,即在(21)式中,將波函數限制在一個以kmax為半徑的球形區域內:

然后將(23)式代入(20)式得

其中q=kmax?/2,則關于時間步長的不等式可寫成如下形式:

根據(25)式,時間步長Δt與函數Q中的θ 和φ有關,Q(θ,φ) 隨θ,φ的變化如圖1 所示,通過數值求解得到函數Q(θ,φ) 在θ=0.304π,φ=0.25π 取得最大值.將θ=0.304π,φ=0.25π 代入(20)式得到Qmax(0.304π,0.25π),然后根據(21)式得到關于時間步長的不等式為

圖1 kmax ? =0.1π 時,Q(θ,φ )隨θ,φ 的變化Fig.1.Variation of Q(θ,φ) with θ,φ at kmax ? =0.1π.

其中ΔtCFL是傳統SFDTD(3,4)方法求解薛定諤方程的最大離散時間步長;CE 表示傳統SFDTD(3,4)算法的CFL 條件被擴展的倍數.圖2 是關于CE 的數值變化曲線,可以發現CE 大于1.這表明通過濾除SFDTD(3,4)方法中空間頻域場分量中的高頻分量,可以使得傳統SFDTD(3,4)方法的CFL 得到進一步擴展.

圖2 擴展因子CE 隨kmax ? 的變化Fig.2.Value of CE varying with kmax ? .

當CE 的值確定后,根據(24)式和(27)式可以得到kmax?的值,然后定義低通濾波器半徑R,其計算公式為

圖3 為采用不同CE (CE1R2>R3>R4).根據圖2 和圖3 可知,CE 的值越大時,R越小.這表明若采用較大的時間步長,需要濾除更多的空間頻域高頻分量來保證數值算法的穩定性.

圖3 不同CE 的值所確定的濾波半徑Fig.3.Filter radius determined by different values of CE.

3 色散分析

根據文獻[14],SFDTD(3,4)方法求解薛定諤方程的數值色散方程為

另外,薛定諤方程的解析色散關系式為[19]

定義數值色散誤差 Errω:

令?=A/kmax,每個波長的網格數由2π/A給出.同時將(26)式,(29)式和(30)式代入(31)式得

圖4(a),(b)分別為FDTD 方法和SFDTD 方法在CE=1 情況下的數值色散誤差.這里需要說明的,SF-SFDTD 方法在CE=1 情況下的色散誤差與SFDTD 方法的色散誤差是相同的.這是因為SF-SFDTD 方法與SFDTD 方法采用相同的數中額外加入濾波操作過程.若CE=1,此時時間步長滿足時間穩定性條件,則不需要采用濾波操作,則SF-SFDTD 方法與SFDTD 方法具有相同的計算過程,兩種方法所具有的數值色散誤差也是相同的.

圖4 FDTD(2,2),SFDTD(3,4)和SF-SFDTD(3,4)方法的數值色散誤差 (a) FDTD(2,2)(CE=1);(b) SFDTD(3,4)(CE=1);(c) SF-SFDTD(3,4)(CE=5);(d) SF-SFDTD(3,4)(CE=10)Fig.4.Curves of numerical dispersion error for the FDTD(2,2),SFDTD(3,4) and SF-SFDTD(3,4) method: (a) FDTD(2,2)(CE=1);(b) SFDTD(3,4)(CE=1);(c) SF-SFDTD(3,4)(CE=5);(d) SF-SFDTD(3,4)(CE=10).

為了進一步比較三種方法的數值色散誤差,圖4(c),(d)給出了SF-SFDTD 方法在CE=5 和CE=10 情況下的數值色散誤差.可以發現,SFSFDTD 方法的數值色散誤差略大于SFDTD 方法的數值色散誤差,這是因為隨著時間步長的增加,時間域上的采樣分辨率下降,導致數值色散誤差增加.但需要說明的是,隨著時間步長的增加,SFSFDTD(3,4)方法的數值色散誤差增加不明顯,與傳統方法的數值色散誤差處于同一量級.值得注意的是: 即使在CE=10 的情況下,SF-SFDTD(3,4)方法相較于FDTD(2,2)方法仍具有較低的數值色散誤差.這表明在SF-SFDTD(3,4)方法中采用較大的時間步長求解含時薛定諤方程,仍能得到十分準確的數值結果.

4 數值計算

為了驗證本文所提SF-SFDTD(3,4)方法的正確性和有效性,通過數值求解二維和三維薛定諤方程,計算量子阱結構的本征態和特征頻率,并將計算結果與其他方法所得計算結果進行對比.

4.1 二維量子阱

二維量子阱模型的尺寸為Lx=Ly=2.9 mm,空間步長Δx=Δy=0.1 mm,采用鏡像原理構建邊界條件[20],邊界條件ψ(0,y)=ψ(Lx,y)=ψ(x,0)=ψ(x,Ly)=0,系統的勢能V(r)=0.

圖5 和圖6 分別給出了傳統SFDTD(3,4)方法(CE=1)和SF-SFDTD(3,4)方法(CE=2,4,5)計算所得二維量子阱結構的特征頻率和本征態ψ2,2.數值結果表明,SF-SFDTD(3,4)方法采用不同的時間步長(CE=2,4,5)計算所得結果與傳統SFDTD(3,4)方法(CE=1)計算所得結果一致.另外,圖5 還給出了FDTD(2,2)方法的計算結果,結果表明SF-SFDTD(3,4)方法相較于FDTD(2,2)方法具有更高的數值計算精度.對比以上不同方法的計算結果,進一步驗證了本文所提SFSFDTD(3,4)方法的正確性.

圖5 二維薛定諤方程的特征頻率Fig.5.Eigenfrequencies of the two-dimensional Schr?dinger equation.

圖6 二維薛定諤方程在不同時間步長下的本征態(ψ2,2) (a) CE=1;(b) CE=2;(c) CE=4;(d) CE=5Fig.6.Eigenstates of two-dimensional Schr?dinger equation at different time steps (ψ2,2): (a) CE=1;(b) CE=2;(c) CE=4;(d) CE=5.

表1 對比了SFDTD(3,4)方法和SF-SFDTD(3,4)方法的數值計算效率.由表1 可以看出,當CE=5 時,SF-SFDTD(3,4)方法的仿真時間小于傳統SFDTD(3,4)方法的仿真時間,這表明在取較大的時間步長情況下,SF-SFDTD(3,4)方法在計算效率方面相較于SFDTD(3,4)方法更具有優勢.的6 倍時,SF-SFDTD(3,4)方法和SFDTD(3,4)方法的數值計算結果仍較為一致,同時可以發現FDTD(2,2)方法的計算結果偏離其他兩種方法的計算結果,再一次驗證了本文所提方法的正確性.

為了定量分析FDTD(2,2)方法、SFDTD(3,4)方法和SF-SFDTD(3,4)方法計算特征頻率結果的相對誤差,定義相對計算誤差:

4.2 三維量子阱

三維量子阱模型的尺寸為Lx=Ly=Lz=3 mm,空間步長Δx=Δy=Δz=0.1 mm,邊界條件ψ(0,y,z)=ψ(Lx,y,z)=ψ(x,0,z)=ψ(x,Ly,z)=ψ(x,y,0)=ψ(x,y,Lz)=0,其勢能V(r)=0.

圖7 和圖8 分別給出了FDTD(2,2)方法(CE =1),SFDTD(3,4)方 法(CE =1)和SFSFDTD(3,4)方法(CE=2,4,6)計算所得三維量子阱結構的特征頻率和本征態ψ2,2,2.由圖7 和圖8 可知,當SF-SFDTD(3,4)方法所采用的時間步長是傳統SFDTD(3,4)方法最大離散時間步長

圖7 三維薛定諤方程的特征頻率Fig.7.Eigenfrequencies of the three-dimensional Schr?dinger equation.

圖8 三維薛定諤方程在不同時間步下的本征態(ψ2,2,2) (a) CE=1;(b) CE=2;(c) CE=4;(d) CE=6Fig.8.Eigenstates of three -dimensional Schr?dinger equation at different time steps (ψ2,2,2): (a) CE=1;(b) CE=2;(c) CE=4;(d) CE=6.

其中G(f) 和Ganalysis(f) 分別表示三種方法計算所得特征頻率結果和解析結果.如圖9 所示,雖然SF-SFDTD(3,4)方法在CE=2,4,6 情況下,其相對計算誤差大于SFDTD(3,4)方法(CE=1)的相對計算誤差,但保持在較低的誤差水平.另外,即使CE=6,SF-SFDTD(3,4)方法的相對計算誤差仍小于FDTD(2,2)方法的相對計算誤差.

圖9 FDTD(2,2),SFDTD(3,4)和SF-SFDTD(3,4)方法計算所得前四個特征頻率點的相對計算誤差Fig.9.Relative calculation errors of the first four characteristic frequency points were calculated by FDTD(2,2),SFDTD(3,4) and SF-SFDTD(3,4) methods.

表2 列出了SFDTD(3,4)方法和SF-SFDTD(3,4)方法取不同時間步長下的仿真時間.由表2可以看出,當CE=6 時,SF-SFDTD(3,4)方法的仿真時間小于傳統SFDTD(3,4)方法的仿真時間,這表明在三維仿真環境下,SF-SFDTD(3,4)方法在取較大的時間步長情況下相較于SFDTD(3,4)方法更具有計算優勢.

表2 SFDTD(3,4)和SF-SFDTD(3,4)方法數值求解三維薛定諤方程的運行時間的運行時間Table 2. Execution time(s) for SFDTD(3,4) and SF-SFDTD(3,4) methods for solving three-dimensional Schr?dinger equation.

5 結論

本文提出了一種新的高穩定性SF-SFDTD(3,4)方法用于數值計算含時薛定方程.相較于隱式無條件穩定的時域數值計算方法,SF-SFDTD(3,4)方法不改變傳統SFDTD(3,4)方法的迭代公式,只需要在每一次的數值計算過程中加入空間濾波操作,濾除因采用不滿足CFL 穩定性條件的時間步長而產生的空間高頻分量.因此,該方法與傳統SFDTD(3,4)方法具有較高的兼容性,數值實現較為容易.另外,文中給出了SF-SFDTD(3,4)方法的數值穩定性和色散誤差分析,并通過數值算例驗證了該方法的正確性和有效性.總的來說,本文所做工作為含時薛定諤方程的求解提供了一種高效的數值求解器,也拓展了時域數值方法在量子計算領域的應用范圍.未來可以將本文所提SFSFDTD(3,4)方法用于數值求解復雜的麥克斯韋-薛定諤耦合方程,解決半經典理論框架下的量子-電磁問題.