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 finite-difference 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]

$${\mathrm{i}}\hbar \frac{{\partial \psi (r,t)}}{{\partial t}} = - \frac{{{\hbar ^2}}}{{2{m^ * }}}{\nabla ^2}\psi (r,t) + V(r)\psi (r,t) ,$$

(1)

其中$\psi$是与空间r和时间t有关的波函数, m*是粒子的质量, $- \dfrac{{{\hbar ^2}}}{{2{m^ * }}}{\nabla ^2} + V$是哈密顿算符. 为避免在数值计算中使用复数, 将$\psi (r, t)$分解为实分量和虚分量, 即:

$$\psi (r,t) = {\psi _{\mathrm{R}}}(r,t) + {\mathrm{i}}{\psi _{\mathrm{I}}}(r,t) .$$

(2)

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

$$\begin{split}&\hbar \frac{\partial {\psi }_{{\mathrm{R}}}(r,t)}{\partial t}=-\frac{{\hbar }^{2}}{2m^*}\left[\frac{{\partial }^{2}{\psi }_{{\mathrm{I}}}(r,t)}{{\partial }^{2}x}+\frac{{\partial }^{2}{\psi }_{{\mathrm{I}}}(r,t)}{{\partial }^{2}y}+\frac{{\partial }^{2}{\psi }_{{\mathrm{I}}}(r,t)}{{\partial }^{2}z}\right]+V(r){\psi }_{{\mathrm{I}}}(r,t) ,\\ &\hbar \frac{\partial {\psi }_{{\mathrm{I}}}(r,t)}{\partial t}=\frac{{\hbar }^{2}}{2m^*}\left[\frac{{\partial }^{2}{\psi }_{{\mathrm{R}}}(r,t)}{{\partial }^{2}x}+\frac{{\partial }^{2}{\psi }_{{\mathrm{R}}}(r,t)}{{\partial }^{2}y}+\frac{{\partial }^{2}{\psi }_{{\mathrm{R}}}(r,t)}{{\partial }^{2}z}\right]-V(r){\psi }_{{\mathrm{R}}}(r,t) .\end{split}$$

(3)

根据文献[14], 采用3步3阶的显式辛算法, 即SFDTD(3, 4)方法数值求解含时薛定谔方程, 其实部和虚部波函数的迭代公式可离散成如下形式:

$$\begin{split} \psi _{{\mathrm{R}}}^{n + \frac{l}{m}}(i,j,k) =\;& \psi _{\text{R}}^{n + \frac{{(l - 1)}}{m}}(i,j,k) + \frac{{V(i,j,k)\Delta t}}{h} \times \psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i,j,k) \\ & - {\alpha _{{\mathrm{R}}x1}} \times \left[ {\psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i + 1,j,k) - 2\psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i,j,k) + \psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i - 1,j,k)} \right] \\ & - {\alpha _{{\mathrm{R}}x2}} \times \left[ {\psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i + 2,j,k) - 2\psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i,j,k) + \psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i - 2,j,k)} \right] \\ & - {\alpha _{{\mathrm{R}}y1}} \times \left[ {\psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i,j + 1,k) - 2\psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i,j,k) + \psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i,j - 1,k)} \right] \\ & - {\alpha _{{\mathrm{R}}y2}} \times \left[ {\psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i,j + 2,k) - 2\psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i,j,k) + \psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i,j - 2,k)} \right] \\ & - {\alpha _{{\mathrm{R}}z1}} \times \left[ {\psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i,j,k + 1) - 2\psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i,j,k) + \psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i,j,k - 1)} \right] \\ & - {\alpha _{{\mathrm{R}}z2}} \times \left[ {\psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i,j,k + 2) - 2\psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i,j,k) + \psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i,j,k - 2)} \right] \end{split} ,$$

(4)

$${\alpha _{{\mathrm{R}}x1}} = \frac{4}{3}{c_l} \times {S_x},\;\;{\text{ }}{\alpha _{{\mathrm{R}}y1}} = \frac{4}{3}{c_l} \times {S_y},\;\;{\text{ }}{\alpha _{{\mathrm{R}}z1}} = \frac{4}{3}{c_l} \times {S_z} ,$$

(5)

$${\alpha _{{\mathrm{R}}x2}} = - \frac{1}{{12}}{c_l} \times {S_x},\;\;{\text{ }}{\alpha _{{\mathrm{R}}y2}} = - \frac{1}{{12}}{c_l} \times {S_y},\;\;{\text{ }}{\alpha _{{\mathrm{R}}z2}} = - \frac{1}{{12}}{c_l} \times {S_z} ,$$

(6)

$${S_x} = \frac{{\hbar {{\Delta }}t}}{{2{m^*}\Delta {x^2}}},\;\;{\text{ }}{S_y} = \frac{{\hbar {{\Delta }}t}}{{2{m^*}\Delta {y^2}}},\;\;{\text{ }}{S_z} = \frac{{\hbar {{\Delta }}t}}{{2{m^*}\Delta {z^2}}} ,$$

(7)

$$\begin{split} \psi _{{\mathrm{I}}}^{n + \frac{l}{m}}(i,j,k) =\;& \psi _{{\mathrm{I}}}^{n + \frac{{(l - 1)}}{m}}(i,j,k) - \frac{{V(i,j,k)\Delta t}}{h} \times \psi _{\text{R}}^{n + \frac{l}{m}}(i,j,k) \\ &+ {\alpha _{{\text{I}}x1}} \times \left[ {\psi _{\text{R}}^{n + \frac{l}{m}}(i + 1,j,k) - 2\psi _{\text{R}}^{n + \frac{l}{m}}(i,j,k) + \psi _{\text{R}}^{n + \frac{l}{m}}(i - 1,j,k)} \right] \\ & + {\alpha _{{\text{I}}x2}} \times \left[ {\psi _{\text{R}}^{n + \frac{l}{m}}(i + 2,j,k) - 2\psi _{\text{R}}^{n + \frac{l}{m}}(i,j,k) + \psi _{\text{R}}^{n + \frac{l}{m}}(i - 2,j,k)} \right] \\ & + {\alpha _{{\text{I}}y1}} \times \left[ {\psi _{\text{R}}^{n + \frac{l}{m}}(i,j + 1,k) - 2\psi _{\text{R}}^{n + \frac{l}{m}}(i,j,k) + \psi _{\text{R}}^{n + \frac{l}{m}}(i,j - 1,k)} \right] \\ & + {\alpha _{{\text{I}}y2}} \times \left[ {\psi _{\text{R}}^{n + \frac{l}{m}}(i,j + 2,k) - 2\psi _{\text{R}}^{n + \frac{l}{m}}(i,j,k) + \psi _{\text{R}}^{n + \frac{l}{m}}(i,j - 2,k)} \right] \\ & + {\alpha _{{\text{I}}z1}} \times \left[ {\psi _{\text{R}}^{n + \frac{l}{m}}(i,j,k + 1) - 2\psi _{\text{R}}^{n + \frac{l}{m}}(i,j,k) + \psi _{\text{R}}^{n + \frac{l}{m}}(i,j,k - 1)} \right] \\ & + {\alpha _{{\text{I}}z2}} \times \left[ {\psi _{\text{R}}^{n + \frac{l}{m}}(i,j,k + 2) - 2\psi _{\text{R}}^{n + \frac{l}{m}}(i,j,k) + \psi _{\text{R}}^{n + \frac{l}{m}}(i,j,k - 2)} \right] \end{split} ,$$

(8)

$${\alpha _{{\text{I}}x1}} = \frac{4}{3}{d_l} \times {S_x},\;\;{\text{ }}{\alpha _{{\text{I}}y1}} = \frac{4}{3}{d_l} \times {S_y},\;\;{\text{ }}{\alpha _{{\text{I}}z1}} = \frac{4}{3}{d_l} \times {S_z} ,$$

(9)

$${\alpha _{{\text{I}}x2}} = - \frac{1}{{12}}{d_l} \times {S_x},\;\;{\text{ }}{\alpha _{{\text{I}}y2}} = - \frac{1}{{12}}{d_l} \times {S_y},\;\;{\text{ }}{\alpha _{{\text{I}}z2}} = - \frac{1}{{12}}{d_l} \times {S_z} ,$$

(10)

$${S_x} = \frac{{\hbar {{\Delta }}t}}{{2{m^*}\Delta {x^2}}},\;\;{\text{ }}{S_y} = \frac{{\hbar {{\Delta }}t}}{{2{m^*}\Delta {y^2}}},\;\;{\text{ }}{S_z} = \frac{{\hbar {{\Delta }}t}}{{2{m^*}\Delta {z^2}}} ,$$

(11)

其中c_l, d_l 为辛传播因子^[14]; n表示时间步, m代表单个时间迭代中子迭代步数, l表示在m个子迭代步数中的第l步计算, i, j, k表示坐标.

2.2   SFDTD(3, 4)方法稳定性条件的扩展

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

$$\frac{\partial }{{\partial t}}\left( {\begin{array}{*{20}{c}} {{\psi _{\text{R}}}} \\ {{\psi _{\text{I}}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} 0&{ - \dfrac{\hbar }{{2{m^ * }}}\left( {{\eta _x} + {\eta _y} + {\eta _z}} \right)} \\ {\dfrac{\hbar }{{2{m^ * }}}\left( {{\eta _x} + {\eta _y} + {\eta _z}} \right)}&0 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{\psi _{\text{R}}}} \\ {{\psi _{\text{I}}}} \end{array}} \right) ,$$

(12)

其中,

$$\begin{split} \eta _x^{} + \eta _y^{} + \eta _z^{}=\;& - 4 \times \bigg\{ {{{\left( {\frac{1}{{\Delta x}}} \right)}^2}\left[ {\frac{4}{3}{{\sin }^2}\left( {\frac{{{k_x}{{\Delta }}x}}{2}} \right) - \frac{1}{{12}}{{\sin }^2}\left( {{k_x}{{\Delta }}x} \right)} \right]} + {\left( {\frac{1}{{\Delta y}}} \right)^2}\bigg[ \frac{4}{3}{{\sin }^2}\left( {\frac{{{k_y}{{\Delta }}y}}{2}} \right) \\ & - \frac{1}{{12}}{{\sin }^2}\left( {{k_y}{{\Delta }}y} \right) \bigg] +{\left( {\frac{1}{{\Delta z}}} \right)}^2\left[ {\frac{4}{3}{{\sin }^2}\left( {\frac{{{k_z}{{\Delta }}z}}{2}} \right) - \frac{1}{{12}}{{\sin }^2}\left( {{k_z}{{\Delta }}z} \right)} \right] \bigg\}. \end{split}$$

(13)

这里k_x, k_y, k_z表示相应方向的波数.

对(12)式进行3阶时间辛积分近似处理, 则其系统时间演化矩阵的离散形式S为

$$\begin{split} & {\boldsymbol{S}}=\left[\begin{array}{cc}{l}_{11} & {l}_{12}\\ {l}_{21}& {l}_{22}\end{array}\right] \\ =\;& \underset{l=1}{\overset{m}{{{\displaystyle \prod }}^{\text{}}}}\left( \begin{array}{cc}1& 0\\ \dfrac{\hbar }{2{m}^{\ast }}\left({\eta }_{x}+{\eta }_{y}+{\eta }_{z}\right){d}_{l}\Delta t & 1\end{array} \right)\\ &\times \left( \begin{array}{cc}1& -\dfrac{\hbar }{2{m}^{\ast }}\left({\eta }_{x}+{\eta }_{y}+{\eta }_{z}\right){c}_{l}\Delta t\\ 0& 1\end{array} \right). \end{split}$$

(14)

由于离散矩阵是行列式值为1的辛矩阵, 矩阵特征值λ对应的特征方程可以写为

$${\lambda ^2} - {\text{tr}}\left( {\boldsymbol{S}} \right)\lambda + 1 = 0 ,$$

(15)

(15)式的解为

$${\lambda _{1,2}} = \frac{{{\text{tr}}\left( {\boldsymbol{S}} \right) \pm {{\text{j}}_0}\sqrt {4 - {\text{tr}}{{\left( {\boldsymbol{S}} \right)}^2}} }}{2} .$$

(16)

稳定性条件要求|λ_{1, 2}| = 1, 矩阵S的迹满足如下条件:

$$\left| {{\text{tr}}({\boldsymbol{S}})} \right| \leqslant 2 .$$

(17)

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

$${\text{tr}}({\boldsymbol{S}}) = 2 + \sum\limits_{l = 1}^m {{{\left( { - 1} \right)}^l}{g_l}{{\left[ {{{\left( {\frac{\hbar }{{2{m^*}}}} \right)}^2}\Delta {t^2}{{\left( {{\eta _x} + {\eta _y} + {\eta _z}} \right)}^2}} \right]}^l}} ,$$

(18)

$$\begin{split} \;&{g_l} = \sum\limits_{1 \leqslant {i_1} \leqslant {j_1} < {i_2} \leqslant {j_2} < \cdots < {i_l} \leqslant {j_l} \leqslant m} {{c_{{i_1}}}{d_{{j_1}}}{c_{{i_2}}}{d_{{j_1}}} \cdots {c_{{i_l}}}{d_{{j_l}}}} \\ & ~+ \sum\limits_{1 \leqslant {i_1} < {j_1} \leqslant {i_2} < {j_2} \leqslant \cdots \leqslant {i_l} < {j_l} \leqslant m} {{d_{{i_1}}}{c_{{j_1}}}{d_{{i_2}}}{c_{{j_1}}} \cdots {d_{{i_l}}}{c_{{j_l}}}} , \end{split}$$

(19)

其中 g_l 的值可以通过辛传播因子 c_l 和 d_l 来计算, m是单个时间迭代中子迭代步数.

令Δx = Δy = Δz = $\varDelta$和

$$\begin{split} Q =\;& \left\{ {\left[ {\frac{4}{3}{{\sin }^2}\left( {\frac{{{k_x}{{\Delta }}x}}{2}} \right) - \frac{1}{{12}}{{\sin }^2}\left( {{k_x}{{\Delta }}x} \right)} \right]} \right. \\ & + \left[ {\frac{4}{3}{{\sin }^2}\left( {\frac{{{k_y}{{\Delta }}y}}{2}} \right) - \frac{1}{{12}}{{\sin }^2}\left( {{k_y}{{\Delta }}y} \right)} \right] \\ & \left.+ \left[ {\frac{4}{3}{{\sin }^2}\left( {\frac{{{k_z}{{\Delta }}z}}{2}} \right) - \frac{1}{{12}}{{\sin }^2}\left( {{k_z}{{\Delta }}z} \right)} \right]\right\} , \end{split}$$

(20)

将其代入(17)式和(18)式, 则关于时间步长Δt的不等式需满足如下形式:

$$\Delta t \leqslant 2.26 \times \frac{{{m^*}{\varDelta ^2}}}{\hbar } \times \frac{1}{Q} .$$

(21)

为确保数值稳定性, 时间步长Δt应小于(21)式右边函数的最小值, 令函数Q中${\sin ^2}\left( {\dfrac{{{k_u}\Delta u}}{2}} \right) = 1$, ${\sin ^2}\left( {{k_u}\Delta u} \right) = 0\left( {u = x, y, z} \right)$, 可得传统SFDTD(3, 4)算法的CFL稳定性条件:

$$\Delta t \leqslant {0}{.565} \times \frac{{{m^*}{\varDelta ^2}}}{\hbar } .$$

(22)

(22)式是确保SFDTD(3, 4)方法中所有空间频域场分量保持稳定, 然而根据文献[16,17], 只有空间频域中一部分的低频分量对计算结果有影响, 而剩余部分的高频分量的幅值较小, 同时对数值计算结果的正确性影响几乎为零. 然而, 若采用较大的时间步长, 会导致空间频域中高频分量的幅值呈现指数增加, 并导致计算结果发散. 基于此, 将SFDTD(3, 4)方法中空间频域场分量限制在一个范围, 即在(21)式中, 将波函数限制在一个以k_max为半径的球形区域内:

$$\begin{split}& \left( {{k_x},{k_y},{k_z}} \right) \\ =\;& \left( {{k_{{\text{max}}}}\sin \theta \cos \varphi ,{k_{{\text{max}}}}\sin \theta \sin \varphi ,{k_{{\text{max}}}}\cos \theta } \right) .\end{split}$$

(23)

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

$$\begin{split} &Q\left( {\theta ,\varphi } \right) \\ =\;& \left\{ {\left[ {\frac{4}{3}{{\sin }^2}\left( {q\sin \theta \cos \varphi } \right) - \frac{1}{{12}}{{\sin }^2}\left( {2q\sin \theta \cos \varphi } \right)} \right]} \right. \\ & +\left[ {\frac{4}{3}{{\sin }^2}\left( {q\sin \theta \sin \varphi } \right) - \frac{1}{{12}}{{\sin }^2}\left( {2q\sin \theta \sin \varphi } \right)} \right] \\ &+ \left. {\left[ {\frac{4}{3}{{\sin }^2}\left( {q\cos \theta } \right) - \frac{1}{{12}}{{\sin }^2}\left( {2q\cos \theta } \right)} \right]} \right\} ,\\[-1pt] \end{split}$$

(24)

其中q = k_max$\varDelta$/2, 则关于时间步长的不等式可写成如下形式:

$$\Delta t \leqslant 2.26 \times \frac{{{m^*}{\varDelta ^2}}}{\hbar } \times \frac{1}{{Q\left( {\theta ,\varphi } \right)}} .$$

(25)

根据(25)式, 时间步长Δt与函数Q中的θ和$\varphi$有关, $Q\left( {\theta , \varphi } \right)$随θ, $\varphi$的变化如图1所示, 通过数值求解得到函数$Q\left( {\theta , \varphi } \right)$在θ = 0.304π, $\varphi$ = 0.25π取得最大值. 将θ = 0.304π, $\varphi$ = 0.25π代入(20)式得到${Q_{\max }}\left( {0.304{\text{π }}, 0.25{\text{π }}} \right)$, 然后根据(21)式得到关于时间步长的不等式为

$$\begin{split}\Delta {t_{{\text{SF}}}} \leqslant\;& 2.26 \times \frac{{{m^*}{\varDelta ^2}}}{\hbar } \times \frac{1}{{{Q_{\max }}\left( {0.304{\text{π}},0.25{\text{π}}} \right)}}\\ =\;&\Delta {t_{{\text{CFL}}}} \times {\text{CE}} ,\end{split}$$

(26)

$${\text{CE}} = {4}/({{{Q_{\max }}\left( {0.304{\text{π }},0.25{\text{π }}} \right)}}) ,$$

(27)

其中Δt_CFL是传统SFDTD(3, 4)方法求解薛定谔方程的最大离散时间步长; CE表示传统SFDTD(3, 4)算法的CFL条件被扩展的倍数. 图2是关于CE的数值变化曲线, 可以发现CE大于1. 这表明通过滤除SFDTD(3, 4)方法中空间频域场分量中的高频分量, 可以使得传统SFDTD(3, 4)方法的CFL得到进一步扩展.

当CE的值确定后, 根据(24)式和(27)式可以得到k_max$\varDelta$的值, 然后定义低通滤波器半径R, 其计算公式为

$$R = \frac{{{k_{{\text{max}}}}\varDelta }}{{2{\text{π }}}} .$$

(28)

图3为采用不同CE (CE₁ < CE₂ < CE₃ < CE₄)的值所确定的滤波半径R (R₁ > R₂ > R₃ > R₄). 根据图2和图3可知, CE的值越大时, R越小. 这表明若采用较大的时间步长, 需要滤除更多的空间频域高频分量来保证数值算法的稳定性.

3.   色散分析

根据文献[14], SFDTD(3, 4)方法求解薛定谔方程的数值色散方程为

$$\tilde \omega \Delta t = \arccos \left[ {{\text{tr}}{{({\boldsymbol{S}})} /2}} \right] .$$

(29)

另外, 薛定谔方程的解析色散关系式为^[19]

$$\omega = \frac{\hbar }{{2{m^*}}}\left( {k_x^2 + k_y^2 + k_z^2} \right) .$$

(30)

定义数值色散误差${\text{Er}}{{\text{r}}_\omega }$:

$${\text{Er}}{{\text{r}}_\omega } = \left| {\frac{{\omega - \tilde \omega }}{\omega }} \right| .$$

(31)

令$\varDelta$ = A/k_max, 每个波长的网格数由2π/A给出. 同时将(26)式, (29)式和(30)式代入(31)式得

$$\begin{split}{\text{Er}}{{\text{r}}_\omega } =\;& \Bigg| 1 - \frac{1}{{0.2825{A^2}{\text{CE}}}}\arccos \Bigg( 1 + 0.5\sum\limits_{l = 1}^m {{{\left( { - 1} \right)}^l}{g_l}} \\ &\times{{\left[ {{{\left( {1.13 \times {\text{CE}} \times Q} \right)}^2}} \right]}^l} \Bigg) \Bigg| ,\\[-1pt] \end{split}$$

(32)

其中

$$\begin{split} Q=&\left\{ {\left[ {\frac{4}{3}{{\sin }^2}\left( {\frac{{A\sin \theta \cos \varphi }}{2}} \right) - \frac{1}{{12}}{{\sin }^2}\left( {A\sin \theta \cos \varphi } \right)} \right]} \right. \\ &+ \left[ {\frac{4}{3}{{\sin }^2}\left( {\frac{{A\sin \theta \sin \varphi }}{2}} \right) - \frac{1}{{12}}{{\sin }^2}\left( {A\sin \theta \sin \varphi } \right)} \right] \\ & + \left. {\left[ {\frac{4}{3}{{\sin }^2}\left( {\frac{{A\cos \theta }}{2}} \right) - \frac{1}{{12}}{{\sin }^2}\left( {A\cos \theta } \right)} \right]} \right\} .\\[-1pt] \end{split}$$

(33)

图4(a), (b)分别为FDTD方法和SFDTD方法在CE = 1情况下的数值色散误差. 这里需要说明的, SF-SFDTD方法在CE = 1 情况下的色散误差与SFDTD方法的色散误差是相同的. 这是因为SF-SFDTD方法与SFDTD方法采用相同的数中额外加入滤波操作过程. 若CE = 1, 此时时间步长满足时间稳定性条件, 则不需要采用滤波操作, 则SF-SFDTD方法与SFDTD方法具有相同的计算过程, 两种方法所具有的数值色散误差也是相同的.

为了进一步比较三种方法的数值色散误差, 图4(c), (d)给出了SF-SFDTD方法在CE = 5和CE = 10情况下的数值色散误差. 可以发现, SF-SFDTD方法的数值色散误差略大于SFDTD方法的数值色散误差, 这是因为随着时间步长的增加, 时间域上的采样分辨率下降, 导致数值色散误差增加. 但需要说明的是, 随着时间步长的增加, SF-SFDTD(3, 4)方法的数值色散误差增加不明显, 与传统方法的数值色散误差处于同一量级. 值得注意的是: 即使在CE = 10的情况下, SF-SFDTD(3, 4)方法相较于FDTD(2, 2)方法仍具有较低的数值色散误差. 这表明在SF-SFDTD(3, 4)方法中采用较大的时间步长求解含时薛定谔方程, 仍能得到十分准确的数值结果.

4.   数值计算

为了验证本文所提SF-SFDTD(3, 4)方法的正确性和有效性, 通过数值求解二维和三维薛定谔方程, 计算量子阱结构的本征态和特征频率, 并将计算结果与其他方法所得计算结果进行对比.

4.1   二维量子阱

二维量子阱模型的尺寸为L_x = L_y = 2.9 mm, 空间步长Δx = Δy = 0.1 mm, 采用镜像原理构建边界条件^[20], 边界条件ψ(0, y) = ψ(L_x, y) = ψ(x, 0) = ψ(x, L_y) = 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)方法具有更高的数值计算精度. 对比以上不同方法的计算结果, 进一步验证了本文所提SF-SFDTD(3, 4)方法的正确性.

表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)方法更具有优势.

4.2   三维量子阱

三维量子阱模型的尺寸为L_x = L_y = L_z = 3 mm, 空间步长Δx = Δy = Δz = 0.1 mm, 边界条件ψ(0, y, z) = ψ(L_x, y, z) = ψ(x, 0, z) = ψ(x, L_y, z) = ψ(x, y, 0) = ψ(x, y, L_z) = 0, 其势能V(r) = 0.

图7和图8分别给出了FDTD(2, 2)方法(CE = 1), SFDTD(3, 4)方法(CE = 1)和SF-SFDTD(3, 4)方法(CE = 2, 4, 6)计算所得三维量子阱结构的特征频率和本征态ψ_{2, 2, 2}. 由图7和图8可知, 当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)方法计算特征频率结果的相对误差, 定义相对计算误差:

$${G_{\text{r}}} = \frac{{\left| {G(f) - {G_{{\text{analysis}}}}(f)} \right|}}{{{G_{{\text{analysis}}}}(f)}} ,$$

(34)

其中$G(f)$和${G_{{\text{analysis}}}}(f)$分别表示三种方法计算所得特征频率结果和解析结果. 如图9所示, 虽然SF-SFDTD(3, 4)方法在CE = 2, 4, 6情况下, 其相对计算误差大于SFDTD(3, 4)方法(CE = 1)的相对计算误差, 但保持在较低的误差水平. 另外, 即使CE = 6, SF-SFDTD(3, 4)方法的相对计算误差仍小于FDTD(2, 2)方法的相对计算误差.

表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)方法更具有计算优势.

5.   结　论

本文提出了一种新的高稳定性SF-SFDTD(3, 4)方法用于数值计算含时薛定方程. 相较于隐式无条件稳定的时域数值计算方法, SF-SFDTD(3, 4)方法不改变传统SFDTD(3, 4)方法的迭代公式, 只需要在每一次的数值计算过程中加入空间滤波操作, 滤除因采用不满足CFL稳定性条件的时间步长而产生的空间高频分量. 因此, 该方法与传统SFDTD(3, 4)方法具有较高的兼容性, 数值实现较为容易. 另外, 文中给出了SF-SFDTD(3, 4)方法的数值稳定性和色散误差分析, 并通过数值算例验证了该方法的正确性和有效性. 总的来说, 本文所做工作为含时薛定谔方程的求解提供了一种高效的数值求解器, 也拓展了时域数值方法在量子计算领域的应用范围. 未来可以将本文所提SF-SFDTD(3, 4)方法用于数值求解复杂的麦克斯韦-薛定谔耦合方程, 解决半经典理论框架下的量子-电磁问题.