计算带Lee-Huang-Yang项的三维Gross-Pitaevskii方程基态的能量耗散投影梯度流方法

An energy-decreasing projected gradient flow method for computing ground states of the 3D Gross-Pitaevskii equation with Lee-Huang-Yang correction

摘要

本文提出了一种数值方法，用于计算带Lee-Huang-Yang（LHY）修正项的三维Gross-Pitaevskii方程（GPE）的基态。在三维情况下，GPE方程所对应的能量泛函中的LHY项为|\psi {|}^{\mathrm{5}}，被用于描述玻色-爱因斯坦凝聚中的量子涨落。为了高效稳健地计算带LHY项的GPE所对应的能量泛函的基态，本文首先提出了一种带拉格朗日乘子项的连续投影梯度方法，对这种方法在空间方向上使用傅里叶伪谱法进行离散，在时间方向上采用向后欧拉格式进行离散，从而得到一种线性隐式的向后欧拉格式（PGF-BESP）。该方法的核心思想是：通过引入一个正交条件来隐式地离散拉格朗日乘子，以保证质量守恒（归一化约束），并直接从正交条件给出拉格朗日乘子的显式表达式。这可以保留带拉格朗日乘子项的连续投影梯度流质量守恒和能量耗散的两个重要性质。接着，本文对格式进行了严格的理论分析，在时间步长满足合适的条件下证明了离散能量是单调非增的。然后，本文进行数值实验，以验证理论结果。结果表明，质量的确是守恒的，并且能量是单调衰减的。最后，本文使用该方法计算了不同势能项和不同相互作用参数下的基态构型，发现随着相互作用参数的增大基态密度分布的整体尺寸也随之增加。综上，理论和数值结果都表明PGF-BESP方法可以为计算带LHY修正的GPE基态提供一种稳定、高效且易于实现的数值工具。

Abstract

This paper presents a numerical method for computing the ground states of the three‑dimensional Gross-Pitaevskii equation （GPE） that includes the Lee-Huang-Yang （LHY） correction.The LHY term in the Gross-Pitaevskii energy functional is of the form |\psi {|}^{\mathrm{5}} and accounts for quantum fluctuations in Bose-Einstein condensates （BECs）.First， the continuous projected gradient flow is discretized with a Lagrange multiplier by using a Fourier pseudo-spectral method in space and a backward Euler scheme in time， obtaining a linearly implicit backward Euler scheme （PGF-BESP）.The key idea is to implicitly determine the Lagrange multiplier through an orthogonal condition， thereby guaranteeing mass conservation （the normalization constraint） and providing an explicit expression for the Lagrange multiplier directly from the orthogonal condition.This preserves the two fundamental properties of the continuous normalized gradient flow with a Lagrange multiplier： Mass conservation and energy dissipation.Next， we carry out a rigorous theoretical analysis of the scheme and prove that the discrete energy is monotonically non-increasing under a suitable condition on the time step size.We then perform numerical experiments to verify the theoretical results， and the results show that mass is indeed conserved and the energy decays monotonically.Finally， we use the proposed method to compute ground state configurations for different potential terms and different interaction parameters.We find that as the interaction parameter increases， the overall size of the ground state density distribution also increases， and under the optical lattice potential， multi-layer structures are observed.In summary， both theoretical and numerical results show that the proposed PGF-BESP method provides a stable， efficient and easy to implement numerical tool for ground state computations of GPE with LHY correction.

Graphical abstract

关键词

玻色爱因斯坦凝聚态 / 基态 / Lee-Huang-Yang修正 / 投影梯度流 / 傅里叶谱方法

Key words

Bose-Einstein condensate / ground state / Lee-Huang-Yang correction / projected gradient flow / Fourier spectral method（2020 MSC 65M60）

Gross-Pitaevskii方程（Gross-Pitaevskii Equation， GPE）是描述玻色-爱因斯坦凝聚态（Bose-Einstein Condensate， BEC）^［1］宏观量子行为的基础理论框架^［2-3］。借助平均场近似，它将复杂的多体问题简化为单个粒子波函数所满足的非线性演化方程。Lee-Huang-Yang（LHY）修正源于Bogoliubov理论，描述了量子涨落对系统基态能量的贡献，在超越平均场近似的BEC研究中起着关键作用。在不同空间维度下，LHY修正的数学形式及其物理效应存在本质差异：在二维情形下，LHY修正项呈对数形式；在一维和三维情形下，LHY修正项呈幂函数形式^［4-5］。

在温度T低于临界温度T_c时，带有LHY修正的三维BEC可由复值波函数\psi (\mathit{x},\mathrm{t})精确描述，其演化服从Gross-Pitaevskii方程，无量纲形式为^［6-8］：

\mathrm{i}{\partial }_t\psi =H_{\psi }\psi ,H_{\psi }:=-\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }+V(\mathit{x})+\beta \left|\psi {|}^{\mathrm{2}}+\delta \right|\psi {|}^{\mathrm{3}},\mathit{x}\in {\mathbb{R}}^{\mathrm{3}},t>\mathrm{0}（1）

其中t为时间，\mathit{x}={(x,y,z)}^{\mathrm{T}}为空间坐标，V(\mathit{x})为有界实值外势，且满足\underset{|\mathit{x}|\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}V(\mathit{x})=\mathrm{\infty }，通常取为调和势：

V_{\mathrm{h}\mathrm{a}\mathrm{r}}(\mathit{x})=\frac{\mathrm{1}}{\mathrm{2}}({\gamma }_x^{\mathrm{2}}{x}^{\mathrm{2}}+{\gamma }_y^{\mathrm{2}}{y}^{\mathrm{2}}+{\gamma }_z^{\mathrm{2}}{z}^{\mathrm{2}})（2）

其中{\gamma }_x>\mathrm{0},{\gamma }_{\mathrm{y}}>\mathrm{0}和{\gamma }_z>\mathrm{0}为无量纲常数，分别正比于x、y、z方向的囚禁频率，\beta和\delta表示非线性相互作用强度。

GPE的两个重要守恒量是波函数的模方（归一化）：

\begin{array}{l}N(\psi (\cdot ,t))={‖\psi (\cdot ,t)‖}^{\mathrm{2}}:=\\ {\int }_{{\mathit{ℝ}}^{\mathrm{3}}}|\psi (\cdot ,t){|}^{\mathrm{2}}\mathrm{d}\mathit{x}\equiv \mathrm{1},t\ge \mathrm{0}\end{array}（3）

以及能量：

\begin{array}{l}E(\psi (\cdot ,t)):=\left(\frac{\mathrm{1}}{\mathrm{2}}|\nabla \psi {|}^{\mathrm{2}}+V(\mathit{x})|\psi {|}^{\mathrm{2}}+\right.\\ \left.\frac{\beta }{\mathrm{2}}|\psi {|}^{\mathrm{4}}+\frac{\mathrm{2}\delta }{\mathrm{5}}|\psi {|}^{\mathrm{5}}\right)\mathrm{d}\mathit{x}\equiv E(\psi (,\mathrm{0})),t\ge \mathrm{0}\end{array}（4）

基态定义为使GPE能量泛函达到最小的宏观波函数{\varphi }_g={\varphi }_g(\mathit{x})，即

（5）

与之对应的Euler-Lagrange方程为

\left\{\begin{array}{l}\mu \varphi =H_{\varphi }\varphi ,\mathit{x}\in {\mathbb{R}}^{\mathrm{3}},\\ N(\varphi )=\mathrm{1}\end{array}\right.（6）

特征值\mu可由下式计算：

\begin{array}{l}\mu (\varphi )=\left(\frac{\mathrm{1}}{\mathrm{2}}|\nabla \varphi {|}^{\mathrm{2}}+V(\mathit{x})|\varphi {|}^{\mathrm{2}}+\beta |\varphi {|}^{\mathrm{4}}+\right.\\ \left.\delta |\varphi {|}^{\mathrm{5}}\right)\mathrm{d}\mathit{x}=E(\varphi )+\frac{\beta }{\mathrm{2}}{\int }_{{\mathit{ℝ}}^{\mathrm{3}}}|\varphi {|}^{\mathrm{4}}\mathrm{d}\mathit{x}+\\ \frac{\mathrm{3}\delta }{\mathrm{5}}{\int }_{{\mathit{ℛ}}^{\mathrm{3}}}|\varphi {|}^{\mathrm{5}}\mathrm{d}\mathit{x}\mathrm{。}\end{array}

目前，已有多种数值方法可被用于计算BEC的基态，如Bao和Tang^［9］提出利用有限元法直接极小化带质量约束的能量泛函，进而得到基态解。Bao^［10-12］等提出了归一化梯度流，并构造了不同的格式。Tang^［13-14］提出了预处理黎曼共轭梯度方法。此外，研究者还提出了Sobolev投影梯度流方法^［15-17］、正则化牛顿法^［18］、Gauss-Seidel型方法^［19］、二阶流方法^［20］等方法来计算BEC的基态。但是，对于带LHY项的三维BEC的基态问题，可用的方法相对比较少。本文提出了一种基于投影梯度流框架的线性化向后欧拉时间离散格式（PGF-BESP），并证明离散之后的原能量在时间步长满足某些合理条件下是下降的。

本文的剩余部分安排如下。第1节提出PGF-BESP方法，并对该方法进行时间和空间上的离散， 然后证明其在时间步长满足某些合理条件下离散能量是下降的。第2节验证能量的下降性，并计算不同势能项以及不同参数下的基态。最后，第3节给出全文所得的结论。

1 方法及性质

1.1 投影梯度流

我们使用如下的连续投影梯度流方法（Projected Gradient Flow，PGF）来计算带LHY项的三维BEC^［10-12］：

\begin{array}{c}\left\{\begin{array}{l}{\partial }_t\varphi =\frac{\mathrm{1}}{\mathrm{2}}\Delta \varphi -V\varphi -\beta |\varphi {|}^{\mathrm{2}}\varphi +\delta |\varphi {|}^{\mathrm{3}}\varphi +{\lambda }_{\varphi }(t)\varphi ,\\ \varphi (\mathit{x},\mathrm{0})={\varphi }_{\mathrm{0}}(\mathit{x}),\hspace{1em}\Vert {\varphi }_{\mathrm{0}}{\Vert }^{\mathrm{2}}=\mathrm{1}(\mathrm{7})\end{array}\right.\end{array}

其中{\lambda }_{\varphi }(t)是拉格朗日乘子，通过下面方式来计算：

{\lambda }_{\varphi }(t)={\int }_{{\mathit{ℝ}}^{\mathrm{3}}}\overline{\varphi }{H}_{\varphi }\varphi \mathrm{d}x（8）

加入拉格朗日乘子{\lambda }_{\varphi }(t)是为了确保质量守恒。

在连续情形下，容易证明PGF格式具有能量耗散和质量守恒的性质。但是，在离散情况下，大部分离散格式都不能够保证能量递减。文献［21］提出了一种隐式处理拉格朗日乘子的方法。在没有LHY项的情况下，该方法将拉格朗日乘子{\lambda }_{\varphi }(t)视为一个新的“未知量”而非已知量，随后对该未知量进行隐式离散，使其充当数值拉格朗日乘子，从而在不使用投影的情况下强制质量守恒。由于无需投影，可能更容易发展出使原始总能量递减的格式。因此，该方法可以同时保证质量守恒和能量衰减。但是，该方法在隐式离散拉格朗日乘子项的时候需要计算一个一元二次方程组，因而还需要从一元二次方程的解中选取适当的根，大大增加了算法实现的复杂性。为此我们将提出一种正交性条件来隐式地离散{\lambda }_{\varphi }(t)，以避免求解一元二次方程。同样，我们所提出的格式可以同时保证能量下降和质量守恒，详见1.3节。

1.2 空间离散及其性质

本节中我们对式（7）和（8）进行空间离散。由于外势的束缚作用，波函数\varphi当|\mathit{x}|\to \mathrm{\infty }时指数衰减，因而可将计算区域截断为有界区域\mathrm{𝒟}:=[L_x,R_x]\times [L_y,R_y]\times [L_z,R_z]，并施加周期边界条件。我们采用傅里叶伪谱方法离散\varphi^［22-24］，网格点取为：

{\mathrm{𝒟}}_N=\left\{{\mathit{x}}_{jkl}:=\left(x_j,y_k,z_l\right)=\left(L_x+j{h}_x,L_y+k{h}_y,L_z+l{h}_z\right),\hspace{0.25em}\left(j,k,l\right)\in O_N\right\},

其中网格步长h_{\nu }=(R_{\nu }-L_{\nu })/N_{\nu }，N_{\nu }为偶数，N:=N_x{N}_y{N}_z，O_N=\{(j,k,l)\in {\mathbb{N}}^{\mathrm{3}}\left|\mathrm{0}\le j\le N_x-\mathrm{1},\hspace{0.25em}\mathrm{0}\le k\le N_y-\mathrm{1},\hspace{0.25em}\mathrm{0}\le l\le N_z-\mathrm{1}\}\right.。此外，傅里叶基函数为：

\begin{array}{l}{W}_{p_x{p}_y{p}_z}\left(\mathit{x}\right)=\prod _{\nu =x,y,z}{\mathrm{e}}^{\mathrm{i}{\mu }_{p_{\nu }}^{\nu }\left(\nu -L_{\nu }\right)},\\ p_{\nu }=-\frac{N_{\nu }}{\mathrm{2}},\dots ,\frac{N_{\nu }}{\mathrm{2}}-\mathrm{1},{\mu }_{p_{\nu }}^{\nu }=\frac{\mathrm{2}\mathrm{\pi }{p}_{\nu }}{R_{\nu }-L_{\nu }},\end{array}

傅里叶伪谱逼近及其导数在{\mathit{x}}_{jkl}处定义为：

\begin{array}{l}\varphi ({\mathit{x}}_{jkl})\triangleq \tilde{\varphi }\left({\mathit{x}}_{jkl}\right):=\frac{\mathrm{1}}{N}\sum _{p_x,p_y,p_z}{\left(\tilde{\varphi }\right)}_{p_x{p}_y{p}_z}·\\ W_{p_x{p}_y{p}_z}\left({\mathit{x}}_{jkl}\right),\\ {\partial }_{\nu }\varphi ({\mathit{x}}_{jkl})\triangleq {(⟦{\partial }_{\nu }^{}⟧\tilde{\varphi })}_{jkl}:=\frac{\mathrm{1}}{N}\sum _{p_x,p_y,p_z}(\mathrm{i}{\mu }_{p_{\nu }}^{\nu }){\left(\tilde{\varphi }\right)}_{p_x{p}_y{p}_z}·\\ W_{p_x{p}_y{p}_z}({\mathit{x}}_{jkl}),\nu =x,y,z,\\ {\partial }_{\nu }^{\mathrm{2}}\varphi ({\mathit{x}}_{jkl})\triangleq {(⟦{\partial }_{\nu }^{\mathrm{2}}⟧\tilde{\varphi })}_{jkl}:=\\ -\frac{\mathrm{1}}{N}\sum _{p_x,p_y,p_z}({\mu }_{p_{\nu }}^{\nu }{)}^{\mathrm{2}}{\left(\tilde{\varphi }\right)}_{p_x{p}_y{p}_z}{W}_{p_x{p}_y{p}_z}\left({\mathit{x}}_{jkl}\right),\end{array}

其中的Fourier系数为：

{\widehat{(\tilde{\varphi })}}_{p_x{p}_y{p}_z}=\sum _{j,k,l}\tilde{\varphi }({\mathit{x}}_{jkl})\hspace{0.33em}{\overline{W}}_{p_x{p}_y{p}_z}\left({\mathit{x}}_{jkl}\right)\mathrm{。}

为了方便，下文依然使用\varphi表示\varphi ({\mathit{x}}_{jkl})构成的向量。则离散能量泛函为：

E\left(\varphi \right):=〈-\frac{\mathrm{1}}{\mathrm{2}}⟦\mathrm{\Delta }⟧\varphi +V\odot \varphi +\frac{\beta }{\mathrm{2}}\left(|\varphi {|}^{\mathrm{2}}\right)\odot \varphi +\frac{\mathrm{2}\delta }{\mathrm{5}}\left(|\varphi {|}^{\mathrm{3}}\right)\odot \varphi ,\hspace{0.33em}\varphi 〉,

离散特征值为：

{\lambda }_{\varphi }:=E(\varphi )+\frac{\beta }{\mathrm{2}}\Vert |\varphi {|}^{\mathrm{2}}{\Vert }^{\mathrm{2}}+\frac{\mathrm{3}\delta }{\mathrm{5}}〈|\varphi {|}^{\mathrm{2}},|\varphi {|}^{\mathrm{3}}〉,

这里\odot表示Hadamard乘积，|\varphi {|}^{\mathrm{2}}表示\varphi \odot \varphi，|\varphi {|}^{\mathrm{3}}表示\varphi \odot \varphi \odot \varphi。为简化记号，我们省略\odot。此外，〈\cdot ,\cdot 〉和‖\cdot ‖为欧式内积与欧式范数：

\begin{array}{l}〈U,V〉:=\mathrm{R}\mathrm{e}\left(h_x{h}_y{h}_z\sum _{j,k,l}{U}_{jkl}{\overline{V}}_{jkl}\right),\\ ‖U‖=\sqrt[]{〈U,U〉},\end{array}

其中Re表示取实部。离散的拉普拉斯算子⟦\mathrm{\Delta }⟧定义为⟦\mathrm{\Delta }⟧=⟦{\partial }_x^{\mathrm{2}}⟧+⟦{\partial }_y^{\mathrm{2}}⟧+⟦{\partial }_{\mathrm{z}}^{\mathrm{2}}⟧。为了之后的叙述，我们定义离散L^p和H^{\mathrm{1}}范数为：

\begin{array}{c}\Vert \varphi {\Vert }_{L^p}:={\left(h_x{h}_y{h}_z{\displaystyle \sum _{j=\mathrm{0}}^{N_x-\mathrm{1}}}{\displaystyle \sum _{k=\mathrm{0}}^{N_y-\mathrm{1}}}{\displaystyle \sum _{l=\mathrm{0}}^{N_z-\mathrm{1}}}|{\varphi }_{jkl}{|}^p\right)}^{\mathrm{1}/p}\end{array},

\Vert \varphi {\Vert }_{H^{\mathrm{1}}}:=\sqrt[]{{‖\varphi ‖}^{\mathrm{2}}+{‖⟦{\partial }_x⟧\varphi ‖}^{\mathrm{2}}+{‖⟦{\partial }_y⟧\varphi ‖}^{\mathrm{2}}+{‖⟦{\partial }_z⟧\varphi ‖}^{\mathrm{2}}}。

为方便，此文以下内容将舍掉\mathrm{符号}⟦⟧。

引理1.1 给定\varphi \in {\mathit{ℂ}}^N。对任意u,v\in {\mathrm{C}}^N，有

1） 算子\mathrm{\Delta }是自伴的，即〈u,\mathrm{\Delta }v〉=〈\mathrm{\Delta }u,v〉。因此离散哈密顿量 {\mathrm{\mathscr{H}}}_{\mathrm{\varphi }}=-\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }+V+\mathrm{\beta }|\mathrm{\varphi }{|}^{\mathrm{2}}+\mathrm{\delta }|\mathrm{\varphi }{|}^{\mathrm{3}}也是自伴随的，并且对任意自伴随算子\mathit{𝒢}有：

\begin{array}{c}\begin{array}{c}\begin{array}{l}\mathrm{2}〈\mathit{𝒢}u,u-v〉=\\ 〈\mathit{𝒢}u,u〉-〈\mathit{𝒢}v,v〉+〈\mathit{𝒢}\left(u-v\right),u-v〉,\end{array}\end{array}\end{array}

\begin{array}{l}\mathrm{2}〈\mathit{𝒢}v,u-v〉=〈\mathit{𝒢}u,u〉-〈\mathit{𝒢}v,v〉-\\ 〈\mathit{𝒢}(u-v),u-v〉\end{array}（9）

2） Gagliardo-Nirenberg不等式在p=\mathrm{4,6}时成立（参见文献［25］）：

\Vert v{\Vert }_{L^p}\le C\Vert v{\Vert }_{L^{\mathrm{2}}}^{\mathrm{3}/p-\mathrm{1}/\mathrm{2}}\Vert v{\Vert }_{H^{\mathrm{1}}}^{\mathrm{3}/\mathrm{2}-\mathrm{3}/p}（10）

证明 1） 分部积分可得：

〈u,{\partial }_{x}v〉=-〈{\partial }_{x}u,v〉,〈u,{\partial }_{y}v〉=-〈{\partial }_{y}u,v〉,〈u,{\partial }_{z}v〉=-〈{\partial }_{z}u,v〉\mathrm{。}

因

〈u,\mathrm{\Delta }v〉=-〈{\partial }_{x}u,{\partial }_{x}v〉-〈{\partial }_{y}u,{\partial }_{y}v〉-〈{\partial }_{z}u,{\partial }_{z}v〉=〈\mathrm{\Delta }u,v〉,

离散的哈密顿量{\mathrm{\mathscr{H}}}_{\varphi }是自伴随的。对任意的自伴算子\mathrm{𝒢}，令：

G_u:=\mathrm{2}〈Gu,u-v〉,G_v:=\mathrm{2}〈Gv,u-v〉。

我们可以得到

G_u+G_v=\mathrm{2}〈\mathit{𝒢}u,u-v〉+\mathrm{2}〈\mathit{𝒢}v,u-v〉=\mathrm{2}〈\mathit{𝒢}u,u〉-\mathrm{2}〈\mathit{𝒢}u,v〉+\mathrm{2}〈\mathit{𝒢}v,u〉-\mathrm{2}〈\mathit{𝒢}v,v〉=\mathrm{2}〈\mathit{𝒢}u,u〉-\mathrm{2}〈\mathit{𝒢}u,v〉+\mathrm{2}〈\mathit{𝒢}u,v〉-\mathrm{2}〈\mathit{𝒢}v,v〉=\mathrm{2}〈\mathit{𝒢}u,u〉-\mathrm{2}〈\mathit{𝒢}v,v〉,

\begin{array}{c}\begin{array}{l}{G}_u-G_v=\mathrm{2}〈\mathit{𝒢}u,u-v〉-\mathrm{2}〈\mathit{𝒢}v,u-v〉=\\ \mathrm{2}〈\mathit{𝒢}\left(u-v\right),u-v〉\mathrm{。}\end{array}\end{array}

因此，

\begin{array}{l}{G}_u=\frac{\left[\left(G_u+G_v\right)+\left(G_u-G_v\right)\right]}{\mathrm{2}}=\\ 〈\mathit{𝒢}u,u〉-〈\mathit{𝒢}v,v〉+〈\mathit{𝒢}\left(u-v\right),\left(u-v\right)〉,\end{array}

\begin{array}{l}{G}_v=\frac{\left[\left(G_u+G_v\right)-\left(G_u-G_v\right)\right]}{\mathrm{2}}=\\ 〈\mathit{𝒢}u,u〉-〈\mathit{𝒢}v,v〉-〈\mathit{𝒢}(u-v),(u-v)〉\mathrm{。}\end{array}

2） 的证明与文献［25］中的相应证明类似，为了文章的简便性此处省略。证毕。

1.3 时间离散及其性质

在本节中我们对已经空间离散完之后的（7）和（8）式进行时间离散。考虑如下线性化向后欧拉时间离散格式（PGF-BESP），该格式可显式求解：

求({\varphi }^{n+\mathrm{1}},{\lambda }^{n+\mathrm{1}})，使得：

\begin{array}{c}\left\{\begin{array}{l}\begin{array}{l}\frac{{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n}{\tau }=\frac{\mathrm{1}}{\mathrm{2}}\Delta {\tilde{\varphi }}^{n+\mathrm{1}}-\\ \left(V\left(\mathit{x}\right)+\beta \left|{\varphi }^n{|}^{\mathrm{2}}+\delta \right|{\varphi }^n{|}^{\mathrm{3}}\right){\tilde{\varphi }}^{n+\mathrm{1}}+{\lambda }^{n+\mathrm{1}}{\varphi }^n,\end{array}\\ \begin{array}{l}〈{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n,{\varphi }^n〉=\mathrm{0},\\ {\varphi }^{n+\mathrm{1}}=\frac{{\tilde{\varphi }}^{n+\mathrm{1}}}{‖{\tilde{\varphi }}^{n+\mathrm{1}}‖},n=\mathrm{0,1},....(\mathrm{11})\end{array}\end{array}\right.\end{array}

将（11）式的第一个方程重新整理得到：

\left(I+\tau {\mathrm{\mathscr{H}}}_{{\varphi }^n}\right){\tilde{\varphi }}^{n+\mathrm{1}}={\varphi }^n+\tau {\lambda }^{n+\mathrm{1}}{\varphi }^n，

其中，

{\mathrm{\mathscr{H}}}_{{\varphi }^n}=-\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }+V+\beta |{\varphi }^n{|}^{\mathrm{2}}+\delta |{\varphi }^n{|}^{\mathrm{3}}。

通过引入{\tilde{\varphi }}^{n+\mathrm{1}}={\varphi }_{\mathrm{1}}^{n+\mathrm{1}}+{\lambda }^{n+\mathrm{1}}{\varphi }_{\mathrm{2}}^{n+\mathrm{1}}，并令 {\mathrm{𝒜}}_n:=I+\tau H_{{\varphi }^n}，可得：

{\varphi }_{\mathrm{1}}^{n+\mathrm{1}}:={\mathrm{𝒜}}_n^{-\mathrm{1}}{\varphi }^n,{\varphi }_{\mathrm{2}}^{n+\mathrm{1}}=\tau {\varphi }_{\mathrm{1}}^{n+\mathrm{1}}。

因此，{\tilde{\varphi }}^{n+\mathrm{1}}=(\mathrm{1}+\tau {\lambda }^{n+\mathrm{1}}){\varphi }_{\mathrm{1}}^{n+\mathrm{1}}。由正交条件，可以得到{\lambda }^{n+\mathrm{1}}的显式表达式：

(\mathrm{1}+\tau {\lambda }^{n+\mathrm{1}})〈{\varphi }_{\mathrm{1}}^{n+\mathrm{1}},{\varphi }^n〉=\mathrm{1}\Rightarrow {\lambda }^{n+\mathrm{1}}=\frac{\mathrm{1}}{\tau }\left(\frac{\mathrm{1}}{〈{\varphi }_{\mathrm{1}}^{n+\mathrm{1}},{\varphi }^n〉}-\mathrm{1}\right)=\frac{〈{\varphi }^n-{\varphi }_{\mathrm{1}}^{n+\mathrm{1}},{\varphi }^n〉}{\tau 〈{\varphi }_{\mathrm{1}}^{n+\mathrm{1}},{\varphi }^n〉}（12）

计算可得\underset{\tau \to \mathrm{0}}{\mathrm{l}\mathrm{i}\mathrm{m}}{\lambda }^{n+\mathrm{1}}=〈{\varphi }^n,{\mathrm{\mathscr{H}}}_{{\varphi }^n}{\varphi }^n〉={\lambda }_{{\varphi }^n}。因此，离散的{\lambda }^{n+\mathrm{1}}和连续的{\lambda }_{{\varphi }^n}是相容的。

为了证明离散能量是下降的，我们首先需要对{\tilde{\varphi }}^{n+\mathrm{1}}做一个先验估计。我们有如下的引理：

引理1.2（先验估计） 设\{{\mathrm{\varphi }}^n{\}}_{n\ge \mathrm{0}}由格式（11）生成。当时，有：

{‖{\tilde{\mathrm{\varphi }}}^{n+\mathrm{1}}‖}_{H^{\mathrm{1}}}\le C\sqrt[]{\mathrm{1}+E^n}，

其中C是一个常数。

证明 由于H_{{\varphi }^n}正定，{\mu }_{\mathrm{m}\mathrm{i}\mathrm{n}}({\mathrm{𝒜}}_n)\ge \mathrm{1}，从而‖{\varphi }_{\mathrm{1}}^{n+\mathrm{1}}‖\le \mathrm{1}。展开〈{\varphi }_{\mathrm{1}}^{n+\mathrm{1}},{\varphi }^n〉得：

\begin{array}{l}〈{\varphi }_{\mathrm{1}}^{n+\mathrm{1}},{\varphi }^n〉=〈\left(I+\tau H_{{\varphi }^n}{)}^{-\mathrm{1}}{\varphi }^n,\right.\\ \left(I-\tau H_{{\varphi }^n}+\tau H_{{\varphi }^n}\right){\varphi }^n〉=\\ \mathrm{1}-\tau 〈\left(I+\tau H_{{\varphi }^n}{)}^{-\mathrm{1}}{\varphi }^n,H_{{\varphi }^n}{\varphi }^n\right]=\\ \mathrm{1}-\tau 〈H_{{\varphi }^n}{\varphi }^n,{\varphi }^n〉+{\tau }^{\mathrm{2}}〈H_{{\varphi }^n}{\varphi }^n,{\mathit{𝒜}}_n^{-\mathrm{1}}{H}_{{\varphi }^n}{\varphi }^n〉\ge \\ \mathrm{1}-\tau 〈H_{{\varphi }^n}{\varphi }^n,{\varphi }^n〉\mathrm{。}\end{array}

因为〈H_{{\varphi }^n}{\varphi }^n,{\varphi }^n〉\le \frac{\mathrm{5}}{\mathrm{2}}{E}^n，故〈{\varphi }_{\mathrm{1}}^{n+\mathrm{1}},{\varphi }^n〉\ge \mathrm{1}-\frac{\mathrm{5}}{\mathrm{2}}\tau E^n。取\tau \le \mathrm{1}/(\mathrm{5}{E}^n)可得〈{\varphi }_{\mathrm{1}}^{n+\mathrm{1}},{\varphi }^n〉\ge \mathrm{1}/\mathrm{2}，于是\mathrm{0}\le {\lambda }^{n+\mathrm{1}}\le \mathrm{1}/\tau。因此，

‖{\tilde{\varphi }}^{n+\mathrm{1}}‖=(\mathrm{1}+\tau {\lambda }^{n+\mathrm{1}})‖{\varphi }_{\mathrm{1}}^{n+\mathrm{1}}‖\le \mathrm{2}。

对于梯度项，利用H_{{\varphi }^n}\ge -\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }可得：

{‖\nabla {\varphi }_{\mathrm{1}}^{n+\mathrm{1}}‖}^{\mathrm{2}}\le \mathrm{2}〈H_{{\varphi }^n}{\varphi }_{\mathrm{1}}^{n+\mathrm{1}},{\varphi }_{\mathrm{1}}^{n+\mathrm{1}}〉\le \mathrm{2}〈H_{{\varphi }^n}{\varphi }^n,{\varphi }^n〉\le \mathrm{5}{E}^n\mathrm{。}

其中最后一步我们利用了非线性项可由能量控制的事实。最终，我们有：

{‖{\tilde{\varphi }}^{n+\mathrm{1}}‖}_{H^{\mathrm{1}}}=(\mathrm{1}+\tau {\lambda }^{n+\mathrm{1}}){‖{\varphi }_{\mathrm{1}}^{n+\mathrm{1}}‖}_{H^{\mathrm{1}}}\le \mathrm{2}\sqrt[]{\mathrm{1}+\mathrm{5}{E}^n}\mathrm{。}

证毕。

接下来，我们有如下的离散能量耗散性：

定理1.3 假设\{{\mathrm{\varphi }}^n{\}}_{n\ge \mathrm{0}}由格式（11）生成。则当

时能量非增，即E({\varphi }^{n+\mathrm{1}})\le E({\varphi }^n),\forall \hspace{0.25em}n\ge \mathrm{0}\mathrm{。}

证明 将（11）式的第一个方程与\mathrm{2}{\tilde{\varphi }}^{n+\mathrm{1}}-\mathrm{2}{\varphi }^n作离散内积得：

\mathrm{2}{\tau }^{-\mathrm{1}}{‖{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n‖}^{\mathrm{2}}=-\mathrm{2}〈H_{{\varphi }^n}{\tilde{\varphi }}^{n+\mathrm{1}},{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n〉\mathrm{。}

利用引理1.2的性质，我们可以得到：

\begin{array}{l}〈H_{{\varphi }^n}{\tilde{\varphi }}^{n+\mathrm{1}},{\tilde{\varphi }}^{n+\mathrm{1}}〉-〈H_{{\varphi }^n}{\varphi }^n,{\varphi }^n〉=\\ -\mathrm{2}{\tau }^{-\mathrm{1}}{‖{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n‖}^{\mathrm{2}}-\\ 〈H_{{\varphi }^n}\left({\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n\right),{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n〉\mathrm{。}\end{array}

由于V\ge \mathrm{0}，我们有：

\begin{array}{l}〈H_{{\varphi }^n}({\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n),{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n〉\ge \\ \frac{\mathrm{1}}{\mathrm{2}}{‖\nabla ({\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n)‖}^{\mathrm{2}}+\\ 〈\beta {\left|{\varphi }^n\right|}^{\mathrm{2}}+\delta {\left|{\varphi }^n\right|}^{\mathrm{3}},{\left|{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n\right|}^{\mathrm{2}}〉\end{array}

从而，

\begin{array}{l}〈H_{{\varphi }^n}{\tilde{\varphi }}^{n+\mathrm{1}},{\tilde{\varphi }}^{n+\mathrm{1}}〉-〈H_{{\varphi }^n}{\varphi }^n,{\varphi }^n〉\le \\ (-\mathrm{2}{\tau }^{-\mathrm{1}}+\mathrm{1}/\mathrm{2}){‖{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n‖}^{\mathrm{2}}-\\ \mathrm{1}/\mathrm{2}{‖{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n‖}_{H^{\mathrm{1}}}^{\mathrm{2}}-\\ 〈\beta {\left|{\varphi }^n\right|}^{\mathrm{2}}+\delta {\left|{\varphi }^n\right|}^{\mathrm{3}},{\left|{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n\right|}^{\mathrm{2}}〉\mathrm{。}\end{array}

现在我们来计算能量差。利用能量和Hamiltonian的表达式，有：

\begin{array}{l}E({\tilde{\varphi }}^{n+\mathrm{1}})-E({\varphi }^n)=\left(〈H_{{\varphi }^n}{\tilde{\varphi }}^{n+\mathrm{1}},{\tilde{\varphi }}^{n+\mathrm{1}}〉-\right.\\ \left.〈H_{{\varphi }^n}{\varphi }^n,{\varphi }^n〉\right)+\frac{\beta }{\mathrm{2}}\Vert |{\tilde{\varphi }}^{n+\mathrm{1}}{|}^{\mathrm{2}}-|{\varphi }^n{|}^{\mathrm{2}}{\Vert }^{\mathrm{2}}+\\ \frac{\delta }{\mathrm{5}}〈{\left|{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n\right|}^{\mathrm{2}},\mathrm{2}{\left|{\tilde{\varphi }}^{n+\mathrm{1}}\right|}^{\mathrm{3}}+\mathrm{4}{\left|{\tilde{\varphi }}^{n+\mathrm{1}}\right|}^{\mathrm{2}}\left|{\varphi }^n\right|+\\ \mathrm{6}\left|{\tilde{\varphi }}^{n+\mathrm{1}}\right|{\left|{\varphi }^n\right|}^{\mathrm{2}}+\mathrm{3}{\left|{\varphi }^n\right|}^{\mathrm{3}}〉\le \\ (-\mathrm{2}{\tau }^{-\mathrm{1}}+\mathrm{1}/\mathrm{2}){‖{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n‖}^{\mathrm{2}}-\mathrm{12}{‖{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n‖}_{H^{\mathrm{1}}}^{\mathrm{2}}+\beta {‖{\tilde{\varphi }}^{n+\mathrm{1}}|{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n|‖}^{\mathrm{2}}+\\ \frac{\delta }{\mathrm{5}}〈{\left|{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n\right|}^{\mathrm{2}},\mathrm{2}{\left|{\tilde{\varphi }}^{n+\mathrm{1}}\right|}^{\mathrm{3}}+\mathrm{4}{\left|{\tilde{\varphi }}^{n+\mathrm{1}}\right|}^{\mathrm{2}}\left|{\varphi }^n\right|+\\ \mathrm{6}\left|{\tilde{\varphi }}^{n+\mathrm{1}}\right|{\left|{\varphi }^n\right|}^{\mathrm{2}}-\mathrm{22}{\left|{\varphi }^n\right|}^{\mathrm{3}}〉\mathrm{。}\end{array}

下面估计非线性项。由Hölder不等式、Young不等式及引理1.1，有：

\begin{array}{l}\frac{\delta }{\mathrm{5}}〈{\left|{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n\right|}^{\mathrm{2}},\mathrm{2}{\left|{\tilde{\varphi }}^{n+\mathrm{1}}\right|}^{\mathrm{3}}+\mathrm{4}{\left|{\tilde{\varphi }}^{n+\mathrm{1}}\right|}^{\mathrm{2}}\left|{\varphi }^n\right|+\\ \mathrm{6}\left|{\tilde{\varphi }}^{n+\mathrm{1}}\right|{\left|{\varphi }^n\right|}^{\mathrm{2}}-\mathrm{2}{\left|{\varphi }^n\right|}^{\mathrm{3}}〉\le \\ \frac{\delta }{\mathrm{5}}\left(\frac{\mathrm{20}}{\mathrm{3}}{‖{\tilde{\varphi }}^{n+\mathrm{1}}‖}_{L^{\mathrm{6}}}^{\mathrm{3}}+\frac{\mathrm{10}}{\mathrm{3}}{‖{\varphi }^n‖}_{L^{\mathrm{6}}}^{\mathrm{3}}\right)\Vert {\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n{\Vert }_{L^{\mathrm{4}}}^{\mathrm{2}}\le C\delta \left(‖{\tilde{\varphi }}^{n+\mathrm{1}}‖{‖{\tilde{\varphi }}^{n+\mathrm{1}}‖}_{H^{\mathrm{1}}}^{\mathrm{3}}+{‖{\varphi }^n‖}_{H^{\mathrm{1}}}^{\mathrm{3}}\right)‖{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n‖\\ {‖{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n‖}_{H^{\mathrm{1}}}^{\mathrm{3}/\mathrm{2}}\le C\delta (\mathrm{1}+E^n)\\ ({ϵ}_{\mathrm{2}}^{\mathrm{4}}{‖{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n‖}^{\mathrm{2}}+{ϵ}^{-\mathrm{4}/\mathrm{3}}{‖{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n‖}_{H^{\mathrm{1}}}^{\mathrm{2}}\mathrm{。}\end{array}

取ϵ=\mathrm{4}C\delta (\mathrm{1}+E^n)。我们有：

\begin{array}{l}\frac{\delta }{\mathrm{5}}〈{\left|{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n\right|}^{\mathrm{2}},\mathrm{2}{\left|{\tilde{\varphi }}^{n+\mathrm{1}}\right|}^{\mathrm{3}}+\mathrm{4}{\left|{\tilde{\varphi }}^{n+\mathrm{1}}\right|}^{\mathrm{2}}\left|{\varphi }^n\right|+\\ \mathrm{6}\left|{\tilde{\varphi }}^{n+\mathrm{1}}\right|{\left|{\varphi }^n\right|}^{\mathrm{2}}+\mathrm{3}{\left|{\varphi }^n\right|}^{\mathrm{3}}〉\le \\ C{\delta }^{\mathrm{4}}(\mathrm{1}+E^n{)}^{\mathrm{6}}{‖{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n‖}^{\mathrm{2}}+\\ \frac{\mathrm{1}}{\mathrm{1}}{‖{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n‖}_{H^{\mathrm{1}}}^{\mathrm{2}}\mathrm{。}\end{array}

类似地，对四次方项有：

\begin{array}{l}\beta {‖{\tilde{\varphi }}^{n+\mathrm{1}}|{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n|‖}^{\mathrm{2}}\le C{\beta }^{\mathrm{4}}(\mathrm{1}+E^n{)}^{\mathrm{3}}\\ {‖{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n‖}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}}{‖{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n‖}_{H^{\mathrm{1}}}^{\mathrm{2}}\mathrm{。}\end{array}

因此，我们有

E({\tilde{\varphi }}^{n+\mathrm{1}})-E({\varphi }^n)\le \left(-\mathrm{2}{\tau }^{-\mathrm{1}}+C{\beta }^{\mathrm{4}}(\mathrm{1}+E^n{)}^{\mathrm{3}}+C{\delta }^{\mathrm{4}}(\mathrm{1}+E^n{)}^{\mathrm{6}}\right){‖{\tilde{\varphi }}^{n+\mathrm{1}}-{\varphi }^n‖}^{\mathrm{2}}\mathrm{。}

于是，若

\tau \le \mathrm{m}\mathrm{i}\mathrm{n}\left\{\frac{\mathrm{1}}{\mathrm{5}{E}^n},\hspace{0.33em}\frac{\mathrm{2}}{C{\beta }^{\mathrm{4}}(\mathrm{1}+E^n{)}^{\mathrm{3}}+C{\delta }^{\mathrm{4}}(\mathrm{1}+E^n{)}^{\mathrm{6}}}\right\},

则有E({\tilde{\varphi }}^{n+\mathrm{1}})\le E({\varphi }^n)。由于‖{\tilde{\varphi }}^{n+\mathrm{1}}‖\ge \mathrm{1}，故E({\varphi }^{n+\mathrm{1}})\le E({\tilde{\varphi }}^{n+\mathrm{1}})。因此，E({\varphi }^{n+\mathrm{1}})\le E\left({\varphi }^n\right)\mathrm{。}

最后，我们通过归纳法证明当

时能量耗散成立。当n=0时结论显然。假设对n\le k以及能量耗散成立。则\mathrm{0}\le E^i\le E^{\mathrm{0}}，从而{\tau }_{\mathrm{0}}\le {\tau }_i对i=\mathrm{0,1},\dots ,k成立。于是，当时，有\mathrm{0}\le E^{k+\mathrm{1}}\le E^k。因此，对n=k+\mathrm{1}以及能量依然耗散。于是PGF-BESP格式的能量耗散性得证，即

。

证毕。

2 数值模拟

在计算带LHY修正项的三维BEC的基态时，采用如下的算法终止准则：

\frac{\Vert {\varphi }^{n+\mathrm{1}}-{\varphi }^n{\Vert }_{\mathrm{\infty }}}{\tau }\le {ϵ}_r。

实际计算中，取{\mathrm{ϵ}}_r={\mathrm{10}}^{-\mathrm{7}}。算法在MATLAB （2022b）中实现，运行环境为Ubuntu GNU/Linux，CPU Intel（R）Xeon（R） Gold 6136R @ 3.00GHz，缓存25 MB。所有空间方向均采用均匀网格，h=h_x=h_y=h_z。本文考虑两类囚禁势：调和势（2）和光学晶格势：

V_d^{\mathrm{l}}(x)=V_{\mathrm{h}\mathrm{a}\mathrm{r}}(x)+{\sum }_{\nu }{\kappa }_{\nu }\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}(q_{\nu }\nu )（13）

其中v=x,y,z。计算区域取为{[-\mathrm{8,8}]}^{\mathrm{3}}，空间步长取为h=1/8。初值取为：

\varphi (\mathit{x})=\frac{\mathrm{1}}{\sqrt[\mathrm{4}]{{\pi }^{\mathrm{3}}}}{\mathrm{e}}^{-\frac{|\mathit{x}{|}^{\mathrm{2}}}{\mathrm{2}}}。

算例2.1 本算例旨在测试PGF-BESP格式的迭代步数及CPU时间。考虑调和势（2）下的带三维LHY修正项的BEC，取{\gamma }_x={\gamma }_y={\gamma }_z=\mathrm{1}，\beta =\mathrm{100},\delta =\mathrm{100}，结果参见表1和图1。可以看到，随时间步长的增大，迭代步数和CPU时间都减小，且格式在时间步长大于定理1中给出的时间步长时仍能稳定运行。

算例2.2 本算例旨在测试PGF-BESP格式的能量下降情况。考虑调和势（2）下的带三维LHY修正项的BEC，取{\gamma }_x={\gamma }_y={\gamma }_z=\mathrm{1}，\beta =\mathrm{100},\delta =\mathrm{100},\tau =\mathrm{0.01}，结果参见图2。可以看到，能量的确是随着时间下降的，与我们的分析相符。

算例2.3 考虑调和势（2）下带LHY修正的三维BEC，取 {\gamma }_x={\gamma }_y={\gamma }_z=\mathrm{1}，\beta =\mathrm{100}\mathrm{、}\mathrm{500},\delta =\mathrm{100}\mathrm{、}\mathrm{2000}，结果参见图3。可以看到，在调和势下基态呈现球形结构，且基态构型的整体尺寸随着β或者δ的增大而增大。

算例2.4 考虑三维光学晶格势（13）下带LHY修正的三维BEC，取{\gamma }_x={\gamma }_y={\gamma }_z=\mathrm{1},q_x=q_y=q_z=\mathrm{\pi }/\mathrm{4},{\kappa }_x={\kappa }_y={\kappa }_z=\mathrm{20},\beta =\mathrm{100}\mathrm{、}\mathrm{2000}，\delta =\mathrm{100}\mathrm{、}\mathrm{2000}，结果参见图4和图5。可以看到，在光学晶格势下，基态构型形成多层晶格结构，并且整体尺寸随着\beta或者\delta的增大而增大。

3 结论

本文针对带 LHY修正项的三维GPE的基态计算问题提出了一种基于投影梯度流框架的线性化向后欧拉时间离散格式PGF-BESP。该格式通过引入正交条件隐式处理拉格朗日乘子，既避免了已有方法中求解一元二次方程的繁琐步骤，同时又保证了离散质量守恒与能量递减性。理论分析表明，在适当的时间步长条件下，格式能够保证离散能量单调递减。数值实验验证了格式的能量耗散行为，并展示了不同相互作用参数下调和势与光学晶格势中基态密度分布的变化。本文的方法为带高阶非线性项的GPE基态计算提供了一种稳定、高效且易于实现的数值工具。

参考文献

[1] Bose S N.Plancks gesetz und lichtquantenhypothese ［J］.Zeitschrift für Physik， 1924， 26（1）： 178-181.

[2] Einstein A.Quantentheorie des einatomigen idealen gases ［J］.Sitz Ber Preuss Akad Wiss， 1924， 22： 261-267.

[3] Lieb E H， Seiringer R， Yngvason J.Bosons in a trap： A rigorous derivation of the Gross-Pitaevskii energy functional ［J］.Phys Rev A， 2000， 61（4）： 043602.

[4] Kadau H， Schmitt M， Wenzel M， et al.Observing the Rosensweig instability of a quantum ferrofluid ［J］.Nature， 2016， 530（7589）： 194-197.

[5] Lee T D， Huang K， Yang C N.Eigenvalues and eigenfunctions of a Bose system of hard spheres and its low-temperature properties ［J］.Phys Rev， 1957， 106（6）： 1135-1145.

[6] Petrov D S.Quantum mechanical stabilization of a collapsing Bose-Bose mixture ［J］.Phys Rev Lett， 2015， 115（15）： 155302.

[7] Semeghini G， Ferioli G， Masi L， et al.Self-bound quantum droplets of atomic mixtures in free space ［J］.Phys Rev Lett， 2018， 120（23）： 235301.

[8] Petrov D S， Astrakharchik G E.Ultradilute low-dimensional liquids ［J］.Phys Rev Lett， 2016， 117（10）： 100401.

[9] Bao W， Tang W.Ground state solution of Bose-Einstein condensate by directly minimizing the energy functional ［J］.J Comput Phys， 2003， 187（1）： 230-254.

[10] Bao W， Chern I， Lim F.Efficient and spectrally accurate numerical methods for computing ground and first excited states in Bose-Einstein condensates ［J］.J Comput Phys， 2006， 219（2）： 836-854.

[11] Bao W， Du Q， Zhang Y.Dynamics of rotating Bose-Einstein condensates and its efficient and accurate numerical computation ［J］.SIAM J Appl Math， 2006， 66（3）： 758-786.

[12] Bao W， Cai Y， Wang H.Efficient numerical methods for computing ground states of dipolar Bose-Einstein condensates ［J］.J Comput Phys， 2010， 229（20）： 7874-7892.

[13] Shu Q， Tang Q， Zhang S， et al.A preconditioned Riemannian conjugate gradient method for computing the ground states of arbitrary-angle rotating Bose-Einstein condensates ［J］.J Comput Phys， 2024， 512： 113130.

[14] Antoine X， Levitt A， Tang Q.Efficient spectral computation of the stationary states of rotating Bose-Einstein condensates by the preconditioned nonlinear conjugate gradient method ［J］.J Comput Phys， 2017， 343： 92-109.

[15] Chen Z， Lu J， Lu Y， et al.On the convergence of Sobolev gradient flow for the Gross-Pitaevskii eigenvalue problem ［J］.SIAM J Numer Anal， 2024， 62（2）： 667-691.

[16] Danaila I， Kazemi P.A new Sobolev gradient method for direct minimization of the Gross-Pitaevskii energy with rotation ［J］.SIAM J Sci Comput， 2010， 32（5）： 2447-2467.

[17] Henning P， Peterseim D.Sobolev gradient flow for the Gross-Pitaevskii eigenvalue problem： global convergence and computational efficiency ［J］.SIAM J Numer Anal， 2020， 58（3）： 1744-1772.

[18] Wu X， Wen Z， Bao W.A regularized Newton method for computing ground states of Bose-Einstein condensates ［J］.J Sci Comput， 2017， 73（1）： 303-329.

[19] Chang S， Lin W， Shieh S.Gauss-Seidel-type methods for energy states of a multi-component Bose-Einstein condensate ［J］.J Comput Phys， 2005， 202（1）： 367-390.

[20] Chen H， Dong G， Liu W， et al.Second-order flows for computing the ground states of rotating Bose-Einstein condensates ［J］.J Comput Phys， 2023， 475： 111872.

[21] Feng Z， Shu Q， Tang Q.Mass-preserving and energy-diminishing semi-implicit schemes for computing ground states of rotating Bose-Einstein condensates ［J］.J Sci Comput， 2026， 106（1）： 30.

[22] Antoine X， Bao W， Besse C.Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations ［J］.Comput Phys Commun， 2013， 184（12）： 2621-2633.

[23] Antoine X， Duboscq R.Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates ［J］.J Comput Phys， 2014， 258（1）： 509-523.

[24] Bao W， Shen J.A fourth-order time-splitting Laguerre-Hermite pseudospectral method for Bose-Einstein condensates ［J］.SIAM J Sci Comput， 2005， 26（6）： 2020-2028.

[25] Faou E， Jezequel T.Resonant time steps and instabilities in the numerical integration of Schrödinger equations ［J］.Differ Integral Equ， 2015， 28（3/4）： 221-238.

基金资助

国家重点研发计划(2024YFA1012803)