玻色-爱因斯坦凝聚中的涡旋晶格形成模拟| COMSOL 博客

link管理

链接快照平台

输入网页链接，自动生成快照
标签化管理网页链接

“薛定谔方程”物理场接口

COMSOL Multiphysics 的半导体模块具有一个 薛定谔方程 物理场接口。在一个最简单的应用示例中，它描述了非相对论粒子在势能图景影响下的动力学（参考文献1 ）。它可以使用包络函数近似模拟量子约束固态系统，例如量子阱、量子线和量子点（参考文献2 ）。在本博客文章中，我们将演示如何将薛定谔方程转换成 Gross–Pitaevskii 方程（参考文献3和4 ），以及如何将其用于模拟玻色凝聚系统。

涡旋晶格形成实验

Madison、Chevy、Bretin 和 Dalibard 于 2001 年发表的论文（参考文献5 ）中展示了一系列引人注意的图像（论文中的图3），它生动地证明了在旋转激光场的作用下，玻色-爱因斯坦凝聚原子云中涡旋晶格的成核和形成。在同一图中，他们还绘制了椭圆率 |\tilde\alpha| 的时间演化图（请参阅下面的公式（7）），显示了当进入涡旋晶格的低能态之前，该系统经历了一段明显的动力学不稳定期。先是初始振荡，随后 |\tilde\alpha| 坍塌至接近零。

图3 摘自 Madison 等人的实验论文。（参考文献5 ）。

主要势阱由横向和纵向势阱频率 \omega_t 和 \omega_z 表征：

V_{trap}=m \omega_t^2(x^2+y^2)/2+m\omega_z^2 z^2/2

其中， m 是原子质量。

势阱频率之比 \lambda\equiv\omega_t/\omega_z 为 9.2，如参考文献1 中图 3 的图像说明文字所示。

旋转激光场中的光势可以近似为

V_{laser}=m\omega_t^2(\epsilon_X X^2+\epsilon_Y Y^2)/2

其中， X 和 Y 轴以一定频率 \Omega 围绕 z 轴旋转。

然后，总势 V_{trap}+V_{laser} 通过势阱频率 \omega_{X,Y}^2\equiv\omega_t^2(1+\epsilon_{X,Y}) 和 \omega_z 来表征。

可以调节激光场以产生不同的纵横比 \epsilon ，并将其定义为

\epsilon\equiv(\omega_X^2-\omega_Y^2)/(\omega_X^2+\omega_Y^2)=(\epsilon_X-\epsilon_Y)/(2+\epsilon_X+\epsilon_Y)

特征频率参数 \bar\omega 被定义为

\bar\omega \equiv \sqrt{(\omega_X^2+\omega_Y^2)/2}=\omega_t\sqrt{(2+\epsilon_X+\epsilon_Y)/2}

以作为旋转频率的参考 \Omega ： \Omega\equiv\tilde\Omega \bar\omega 。

例如，参考文献1 中的图3 是当 \tilde\Omega =0.7 时生成的。此外， \bar\omega 也用于量化椭圆率 |\tilde\alpha| ，我们将在下面详细介绍。但是，实验是在纵横比 \epsilon 随时间演化而上下起伏的情况下进行的。因此，如果定义（4）用于 \epsilon_X 和 \epsilon_Y 的任意组合，则在时间演化过程中的 \bar\omega 值也会上下变化。

为了维持 \bar\omega 的值为常数，以作为旋转频率 \Omega 和椭圆率 |\tilde\alpha| 的可靠参考，我们在模型中计算 \bar\omega 时使用 \epsilon_X 和 \epsilon_Y 的常数值分别为 0.03 和 0.09。参考值 0.03 和 0.09 是基于参考文献2 所引用的同一小组的早期实验论文。

最后，椭圆率 |\tilde\alpha| 定义如下。被捕获凝聚体的密度分布可以通过 Thomas–Fermi 分布很好地近似：

\rho_{TF}=max\left[0 , \frac{\mu_{TF}}{g}\left(1-\frac{X^2}{R_X^2}-\frac{Y^2}{R_Y^2}-\frac{z^2}{R_z^2}\right)\right]

其中， \mu_{TF} 是化学势（在 Thomas–Fermi 近似内）， g=4\pi\hbar^2a/m 是由散射长度 a 设置的耦合常数（在|F=2，mF=2 >基态中，对于 ⁸⁷ Rb， a =5.5nm ）。

根据上式（5）横向尺寸 \alpha （ R_X 和 R_Y ）的密度分布来定义参数 a ：

\alpha\equiv\Omega\frac{R_X^2-R_Y^2}{R_X^2+R_Y^2}

然后，将椭圆率参数 \tilde\alpha 定义为

\tilde\alpha\equiv\alpha/\bar\omega=\tilde\Omega\frac{R_X^2-R_Y^2}{R_X^2+R_Y^2}

此模型的建模理论基于 Tsubota 等人的方法（参考文献2 ）。即通过求解具有唯象阻尼的旋转框架中的 Gross–Pitaevskii 方程来模拟这一壮观的时间演化过程（参考文献3 ）：

(i-\gamma)\hbar\frac{\partial\psi}
{\partial t}
=\left[-\frac{\hbar^2}
{2m}\nabla^2+V_{trap}+V_{laser}+g|\psi|^2-\mu-\Omega L_z\right]\psi

其中， \gamma 是阻尼系数， \mu 是要随时调整以保持凝聚原子数的化学势参数，并且 -\Omega L_z=i\hbar\Omega(x\partial_y-y\partial_x) 是旋转框架中的离心项。

Choi 等人（参考文献3 ）提供了估算阻尼系数 \gamma 的公式：

\gamma\approx\frac{4m(a k T)^2}
{\pi\hbar^3}
e^{2\mu/k T}\frac{\mu}
{k T}K_1\left(\frac{\mu}{k T}
\right)/\bar\omega

其中， K_1 是一个修正的贝塞尔函数，并且我们已经假设凝聚体与其周围的热原子之间达到平衡。

Thomas-Fermi 近似可以很好地估算许多重要参数，例如阻尼系数 \gamma 。圆柱形对称势阱中的峰值密度和尺寸参数总结如下：

\begin{align*}
\rho_0\equiv\frac{\mu_{TF}}{g} &= \frac{15N}{8\pi R_r^2 Rz} \\
R_r &= \left(\frac{15g \omega_z N}{4\pi m \bar\omega^3}\right)^{1/5}\\
R_z &= \left(\frac{15g \bar\omega^2 N}{4\pi m \omega_z^4}\right)^{1/5}
\end{align*}

其中， N 是凝聚体的原子数。

模拟涡旋晶格的形成

Gross-Pitaevskii 方程（8）可以使用半导体模块中的 薛定谔方程 物理场接口直接实现。在构建模型时，我们将模型参数与参考文献1 中的实际实验条件进行匹配，以使模拟的时间演变与已发布的数据完全吻合。

根据 Bretin 的论文，搅拌激光场瞬间开启，纵横比 \epsilon 在 20ms 内斜升，然后在 300ms 内保持恒定。因此，在此模型中，我们使搅拌激光场始终处于开启状态，并使纵横比 \epsilon 随时间变化。为了使瞬态和稳态研究共享同一个公式，对于瞬态研究而言，我们定义了一个与内置时间参数具有相同名称t的时间参数。分别将斜升时间 20ms 和停止时间 300ms 定义为参数 tau 和 t_off 。然后定义一个阶跃函数来斜升和斜降纵横比 \epsilon ，假设相同的斜降时间周期为 20ms。公式被设置为

epst=0.032*step1((t+tau)/tau)*(1-step1((t-t_off)/tau))

使斜升在 -20ms 处开始并在时间 t=0 处结束。基于参考文献1 和 Bretin 的论文， \epsilon 的大小被设置为 0.032。

根据论文，我们并不清楚激光场的长半轴和短半轴是如何随纵横比的变化而变化的。但是，假设椭圆的面积在斜升斜降期间保持恒定似乎是合理的。因此，我们从搅拌激光场中获得了光势的 \epsilon_X 和 \epsilon_Y 参数的公式：

epsX=(epst+sqrt(0.03*0.09+epst^2-0.03*0.09*epst^2))/(1-epst) ，

epsY=(-epst+sqrt(0.03*0.09+epst^2-0.03*0.09*epst^2))/(1+epst) 。

准备好这些纵横比参数后，按照实验部分中的说明输入势阱参数。凝聚体中的原子数选择为 1.5e5，这与实验范围一致，并且最适合于实验中的涡旋数。为了估计阻尼系数 \gamma ，根据参考文献1 ，使用的温度为 100nK。

由于我们使用二维模型近似三维凝聚体，因此采用 Thomas-Fermi 公式（10）计算合理的面外厚度。选择面外厚度的标准，以使二维模型中的峰值密度与三维中的 Thomas-Fermi 峰值密度相匹配，我们获得以下面外厚度 L 的公式：

L=\frac{N}{\rho_0 \frac{\pi}{2} R_r^2}

初始静止状态

首先，用与模型示例玻色-爱因斯坦凝聚的 Gross-Pitaevskii 方程相同的两个研究来求解凝聚体的稳态。这为随后的瞬态研究提供了初始条件。

对于公式（8）中的相互作用能项 g|\psi|^2 ，将内置变量 schr.Pr 除以面外厚度 L 以得到三维的粒子密度。（请注意， schr.Pr 的意义与从常规薛定谔方程意义上所说的通常的“概率密度”不同。）在这里，因变量 \psi 代表 Gross–Pitaevskii 凝聚的阶次参数，并且 schr.Pr （ \equiv |\psi|^2 ) 表示二维原子密度。）在稳态研究的波函数 N=\int|\psi|^2 的归一化全局方程中，总能量由 Thomas–Fermi 化学势 \mu_{TF} 标定，因此要求解的全局变量具有统一的阶次。

下图比较了 X 轴上的最终粒子密度分布与 Thomas-Fermi 近似的分布。正如预期的那样，它们完全一致。

X 轴上的粒子密度分布的计算稳态解（蓝色）和通过 Thomas–Fermi 近似求出的一个稳态解（绿色）。

旋转框架、耗散和归一化

获得稳态解后，可以使用了 COMSOL Multiphysics 5.6 版本中提供的两个特征为瞬态研究添加更多的物理场特征。一种是 旋转框架 特征，如下面的截图所示。

旋转框架特征的设置 窗口。

对于这样的二维模型，旋转轴固定在面外方向上。对于三维模型，用户可以选择任意方向上的轴。

另一个是唯象阻尼的耗散特征，这对于使系统弛豫到涡旋晶格的低能态至关重要。请参见下面的屏幕截图。

耗散 特征的 设置 窗口。

依据参考文献2 ，与稳态研究类似，使用全局方程通过调节化学势来维持凝聚体中的原子数，并将要求解的全局变量按与 Thomas–Fermi 化学势 μ _TF 的能量标度统一的阶次进行缩放。

动态不稳定性

在达到低能量的涡旋晶格状态前，该系统经历了一段时间的动态不稳定。这种物理过程的随机本质会导致每次运行的模拟时间历史发生重大变化。所得晶格中的涡旋数也可能变化。为了提高数值收敛性，对求解器设置进行了一些调整。由于初始条件是来自稳态研究的物理解，因此可以禁用一致初始化。它通常有助于将代数状态排除在误差控制之外。具有大量迭代次数的自动牛顿法有助于克服非线性严重时的不稳定时期。

下面的动画显示了作为时间函数计算出的粒子密度分布。在凝聚体振荡/旋转的初始阶段之后，旋涡开始在外围形成。随着涡旋随机移动，一个动态不稳定周期随之而来（在此时间段内，动画会变慢。）最终，系统稳定到涡旋晶格的低能量状态。请欣赏希下面的演示！

作为时间函数的计算粒子密度分布。

下图显示了该动画的一些快照。

作为时间函数的计算粒子密度分布。

由于实验装置中光学成像系统的实际限制，在原子仍被捕获的情况下，无法获得如上图所示的密度分布图。取而代之的是，在实验中，原子从势阱中释放出来，并允许云在 25 ms 内自由扩展到大约 300 um 大小。在自由扩展前后，云的纵横比也发生了巨大变化。最初的雪茄形状变成最终的薄饼形状，在扩展前后长尺寸和短尺寸发生了互换。在将模拟的势阱内密度分布与扩展后的原子云的已发布图像进行比较时，应牢记这一点。

上面的动画和图形中显示的时间演化过程可以简化为单个椭圆率参数 |\tilde\alpha| ，该参数通过将粒子密度分布拟合为一个简单函数来提取椭圆的长轴 R_X 和短轴 R_Y ，然后应用于等式（7）。对于上图所示的模拟的势阱内密度分布，Thomas–Fermi 近似提供了一个很好的拟合函数。通过将其拟合到每个时间点的密度分布，我们可以计算椭圆率参数 |\tilde\alpha| 以作为时间的函数。下图显示了计算结果（蓝点）。初始振荡的时间尺度和的最终塌陷与参考文献1 中图3 显示的数据吻合良好。虽然大小略有不同，但是考虑上述自由扩展前后可能发生的形状变化，这是可以理解的。

椭圆率参数和每个原子的角动量（单位为 \hbar ）。

表征从振荡/旋转的全云到涡旋晶格过渡的另一个重要参数是角动量，上图中对角动量也进行了绘制（绿色曲线）。初始振荡和最终获得一定角动量（与旋涡数量成比例）的一般行为与 Tsubota 等人的模拟结果一致（参考文献2 中的图3 ）。但是，这里模拟结果的时间尺度与实验数据非常接近。

优化模块用于拟合粒子密度分布。为了检查拟合的质量，我们可以绘制拟合数据的等值线（模拟密度分布图）和拟合函数的等值线（Thomas-Fermi密度分布图）并进行比较，如下图所示。

拟合数据（模拟密度分布，灰度）和拟合函数（Thomas-Fermi 密度分布，彩色）。

如下面的截图所示，优化方案的设置从拟合参数的定义开始。

优化方案的拟合参数。

拟合参数分别是：拟合密度分布图的第一轴 RXfit 、第二轴 RYfit 、倾角 thetafit 和峰值密度 rho0fit 。还定义了解参考的索引参数 index 和椭圆率参数 alphafit 的拟合值，使用拟合参数和公式（7）计算。

基于 Thomas–Fermi 密度分布的拟合函数被定义为变量 fit_fn 。然后，将计算的数据与拟合函数之间的差值求平方并在模拟域中求平均值，以作为优化研究要最小化的目标。为了防止目标的大小变得太大而使优化器无法处理，我们通过 Thomas-Fermi 峰值密度（公式（10）中的 \rho_0 ）来缩放差异，以使产生的目标接近于1 的量级。该变量被定义为以下截图中的变量 q0 。

通过优化研究最小化拟合密度分布和目标。

然后在优化研究设置中使用目标 q0 和 4 个拟合参数。使用 Thomas-Fermi 近似中的值，为每个拟合参数提供适当的初始值、比例和边界。参见以下截图。

优化研究设置。

首先使用参数 index 设置参数化扫描，然后选取适合每个时间点的解。在虚设的稳态研究步骤中，设置未求解的变量值栏，以使用 index 参数选择在每个时间步中的瞬态解。详可参见下面截图

使用参数 index 进行参数化扫描。

用于栏设置的带有“未求解变量的值”的虚设稳态研究步骤，以使用参数选择每个时间步中的瞬态解。

本文，我们使用一个由被捕获的冷原子形成的玻色-爱因斯坦凝聚中的涡旋晶格形成的动力学过程模型演示了 COMSOL 软件中的“薛定谔方程”物理场接口。欢迎您在下面的留言区评论您如何使用此物理场接口处理其他有趣的现象！

获取教程模型

1. L. I. Schiff, 《Quantum Mechanics》, McGraw-Hill, ed. 3, 1968.

2. P. Harrison, 《Quantum Wells, Wires and Dots》, Wiley, ed. 3, 2009.

3. E.P. Gross, “Structure of a quantized vortex in boson systems”, Il Nuovo Cimento, vol. 20, no. 3, pp. 454–457, 1961.

4. L.P. Pitaevskii, “Vortex lines in an imperfect Bose gas”, Sov. Phys. JETP, vol. 13, no. 2, pp 451–454, 1961.

5. K. W. Madison, F. Chevy, V. Bretin, and J. Dalibard, “Stationary States of a Rotating Bose-Einstein Condensate: Routes to Vortex Nucleation”, Phys. Rev. Lett., 86, 4443, 2001.

6. M. Tsubota, K. Kasamatsu, and M. Ueda,《“Vortex lattice formation in a rotating Bose-Einstein condensate”,, Phys. Rev., A 65, 023603 (2002).

7. S. Choi, S. A. Morgan, and K. Burnett, “Phenomenological damping in trapped atomic Bose-Einstein condensates”, Phys. Rev., A 57, 4057, 1998.