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凝聚态的基态解问题．