使用 discontinuous Galerkin 时域方法对 Casimir 力进行数值评估

《Physical Review A》：Numerical evaluation of Casimir forces using the discontinuous Galerkin time-domain method

　　卡西米尔力在纳米尺度上主导了电中性宏观物体之间的相互作用，因此对其定量描述对于现代技术和基础物理学具有重要意义。在这项工作中，我们提出了一种基于麦克斯韦应力张量形式主义的时域方案，用于计算真空分离的物体之间的卡西米尔力。特别是，我们使用基于有限元的不连续伽辽金时域方法提供了具体的实现方式。这种方法能够准确评估任意几何形状和一般材料模型下的有限温度卡西米尔-利夫希茨相互作用。该方法的核心是将电磁格林张量表示为系统对偶极激励的响应，从而将麦克斯韦应力张量转化为一系列由电偶极子和磁偶极子驱动的经典散射问题。我们通过参考平行半空间之间卡西米尔相互作用的计算结果，在零温度和有限温度下验证了该方法的有效性。此外，我们还证明了该方法适用于有限圆柱体与半空间之间的相互作用，对于这种相互作用，现有的半解析表达式并不适用，我们的结果与基于物理考虑的渐近预测以及其他数值方法都非常吻合。这些发现展示了该方法在研究与纳米器件设计和实验设置相关的微观和纳米尺度结构中的卡西米尔相互作用方面的潜力。

### 引言

由电磁场的量子涨落和热涨落引起的现象产生了超越经典电动力学范围的迷人效应。一个著名的例子是卡西米尔效应，即通常存在于电中性宏观物体之间的吸引力。这一效应最初由H. B. G. Casimir针对真空中的完美反射镜预测，并归因于量子电磁场的零点涨落[1]，后来Lifshitz将其扩展到了实际材料中，考虑了有限温度下的色散和耗散等特性[2]。卡西米尔力可以被视为整个系统量子涨落的宏观表现，包括材料内部的涨落。从这个意义上说，它包含了范德瓦尔斯相互作用，这一作用最初由范德瓦尔斯提出，并后来由London描述为间距接近的粒子之间的吸引力[3]。Casimir和Polder进一步发展了London的理论，包括了迟滞效应，描述了原子与宏观物体之间的相互作用[4]，由此产生的涨落效应通常被称为卡西米尔-波尔德相互作用。虽然在宏观尺度上可以忽略不计，但在亚微米距离下，卡西米尔力变得显著，并且随着距离的减小而迅速增加。除了理论上的相关性外，卡西米尔效应已通过实验观测得到了一致证实[6]。现代测量依赖于精密仪器，如扭摆[7]、原子力显微镜[8-10]和微机电扭振子[11-15]。实验灵敏度的持续提高使得检测有限温度修正成为可能[16-20]，但仍有一些相关争议尚未解决[13,21,22]。尽管已经报告了成功的平面-平面测量[23]，但实验装置经常使用靠近基底的曲面以最小化对齐误差。在这种配置中，通常使用近距力近似（PFA）[24]来估计相互作用，它将相互作用表面视作局部平行的，并积分相应的无限小平行板贡献。从技术角度来看，纳米结构设备的持续小型化使卡西米尔效应在纳米和微机电系统的设计中越来越重要，引起了实验和理论研究社区的兴趣[25,26]。一个后果是，由于其通常具有吸引力，卡西米尔效应可能导致粘附——即移动部件的不希望出现的粘连——从而破坏设备性能[11,27,28]。

在这种背景下，对卡西米尔相互作用的理论描述对于理解几何和材料依赖性效应、设计纳米结构设备以及优化实验装置至关重要。随着时间的推移，已经开发出基于散射方法的半解析方法来描述实际系统，包括复杂的几何形状、现实的材料模型和有限温度修正[29-35]。然而，这些方法仍然严重依赖于几何形状。作为替代方案，全数值技术提供了一个更灵活的框架，用于处理任意三维配置，而不会增加数学复杂性。常见的方法包括边界元方法（BEMs）[36-39]，以及频率域（FDFD）[40]和时域（FDTD）[41,42]中的有限差分技术[41,42]。这些方法考虑了温度效应[43]、广泛的材料模型和任意几何形状，揭示了有趣的现象，如非单调力行为[44,45]、圆柱体在穿孔板上方时的几何诱导排斥[39,46,47]，或流体中的材料诱导悬浮[48,49]。在这种背景下，对于具有任意几何形状和材料属性的完全三维系统进行数值精确的卡西米尔力计算，同时一致地纳入有限温度效应，在技术上仍然具有挑战性，因为存在多个相关的长度尺度。

在这项工作中，我们扩展了Kristensen等人在参考文献[50]中引入的数值框架，用于高精度卡西米尔-波尔德计算，将卡西米尔相互作用和有限温度修正纳入麦克斯韦应力张量形式主义中。我们的方法采用了不连续伽辽金时域（DGTD）方法[51]，这是一种基于有限元的时间域数值求解麦克斯韦方程的技术。作为代表性示例，我们考虑了图1中显示的几何形状。首先，我们研究了真空分离的金属半空间之间的卡西米尔相互作用，并将我们的结果与基于利夫希茨理论的半解析解进行了比较[2]。我们进一步考虑了由真空分隔的有限圆柱体与半空间之间的相互作用。尽管这种配置之前已经为具有小于分离距离的半径的完美导电无限薄盘半解析地处理过[52]，但有限圆柱体缺乏有效半解析处理所需的对称性。然而，最近有报道使用了BEM方法[38]和快速收敛的多重散射展开方法[53,54]来计算涉及色散和耗散材料的任意几何形状中的卡西米尔力。因此，我们通过确认与圆柱体半径相比的大和小的分离距离处的预期渐近行为，以及通过与其它已知数值方法的比较[55]来验证我们的结果。

图1. 使用Gmsh [5]生成的代表性四面体网格，用于（a）平面-平面和（b）圆柱-平面几何形状。其中，是最小表面到表面的分离距离，和分别表示圆柱体的半径和高度。

当前框架基于Rodriguez等人和McCauley等人在参考文献[41,42]中的早期方法，他们使用FDTD方法计算了人工耗散系统中展开源激发的场。这种方法重新制定了介电常数，有效地变形了复频率积分轮廓，从而减少了所需的计算时间。相比之下，我们的方法保留了原始的电磁问题，并在高斯求积方案的节点处放置了定制的宽带偶极源。这确保了在后处理步骤中准确高效地评估时间和空间积分。与依赖于空间和时间上规则交错网格的标准FDTD方法[56]相比，DGTD方法允许自适应网格划分，并通过局部多项式基函数和高效的时间步进方案实现了更高阶的收敛[51]。

本文的组织结构如下：在第二节中，我们概述了在麦克斯韦应力张量形式主义内评估卡西米尔相互作用的理论框架。第三节描述了数值实现及其与时域经典散射问题的关系。第四节，我们介绍了两个案例研究：（i）平行半空间之间的卡西米尔力，说明了时间积分方案；以及（ii）有限金属圆柱体与半空间之间的相互作用，介绍了处理缺乏平移对称性系统所需的表面积分方案。最后，第五节总结了结果并概述了潜在的未来发展方向。

### 麦克斯韦应力张量方法

我们使用麦克斯韦应力张量形式主义来计算热平衡状态下系统中的涨落诱导力。具体来说，我们考虑了一组任意形状、不重叠的物体，它们嵌入在温度为T的真空中。任何给定物体上受到的净力是由电磁场的量子涨落和热涨落引起的，由在系统的热态下评估的麦克斯韦应力张量算子的期望值的表面积分给出（1）。这里，表示包围物体的任何封闭表面，在该表面上评估力，是位于位置处的向外指的单位法向量。在真空中，应力张量的期望值用同一位置和时间的波动宏观电场和磁场算子的量子相关函数表示，如下[57]：

（2）

其中表示克罗内克delta，和分别是真空的介电常数和磁导率。在热平衡状态下，场相关函数可以使用涨落耗散定理（FDT）来评估（3a）（3b）。尽管场算子在时间上波动，但热状态的稳态意味着它们的等时相关性是时间独立的；因此，等式（1）中的应力张量算子和卡西米尔力的期望值不依赖于时间，为了简便，在本文的其余部分省略了它们的显式时间依赖性。这些表达式将量子场和热场涨落关联到频率T处的电磁格林张量的虚部，其中区分了电场和磁场贡献。

格林张量编码了系统的几何形状和材料属性，并依赖于观察点T，场在此处被评估，以及源位置S，这在系统中引起了散射。对于嵌入在真空中的线性和局部物体，它们满足以下非均匀向量亥姆霍兹方程（4a）（4b），受到出射辐射边界条件的限制（60）。这里，和分别是分段均匀的介电常数和磁导率函数，表示单位二叉量。总的格林张量分解为，其中是真空贡献，在各向同性和均匀真空中解决了方程（4）。尽管在重合极限下会发散，但其对方程（1）的贡献消失，因为它变得独立于，并且在任何封闭表面上积分结果为零。另一方面，散射的格林张量表示在材料不均匀性（如几何边界和界面反射）存在的情况下，对位于S处的点状源的响应。由于不贡献于净卡西米尔力，我们在后续部分仅关注，并为了简便省略了下标。

### 数值实现

方程（3）中的FDT将方程（2）中的麦克斯韦应力张量的期望值与格林张量的虚部联系起来。这种关系将量子电动力学效应转化为由方程（4）控制的经典电动力学问题。然而，如上所述，这些方程的半解析解仅存在于高度对称的几何形状集合中。因此，挑战在于开发能够准确评估任意几何形状的散射格林张量的高效数值方法。一种实用策略是将二叉格林张量的每一列解释为系统对点状偶极源的场响应。具体来说，电格林张量的分量与沿轴排列且在处位于的无限小振荡电偶极子在位置处产生的电场的分量成正比[61]。相应的源电流密度由表示，其中是频域幅度。类似地，磁格林张量通过引入磁偶极子得到。在重合点，格林张量分量由场得出（5），其中和是由沿轴排列的偶极子引起的散射电场和磁场的分量[50]。

通过将方程（3）和（5）代入方程（2），可以用人散射电磁场来表示方程（1）中出现的麦克斯韦应力张量的期望值（6），其中几何依赖的散射场贡献定义为（7）。这里，求和遍及空间分量。在方程（6）中，我们引入了温度依赖函数（9），其中是Heaviside阶跃函数。在极限下，由于作为时间域电流响应在时间的一阶非零导数的阶数，会发散，其傅里叶变换为[50]。尽管其行为会发散，但方程（6）中的虚部算子“”与积分可交换，导致有限的结果。这可以通过将其解释为缓和分布来证明，因为它仅在对平滑、衰减的场函数进行积分时出现。

### 从频率域到时域

通过沿实轴的频率域积分直接数值评估麦克斯韦应力张量受到高度振荡被积函数的阻碍[62]。这些振荡会降低方程（6）中宽带积分的准确性和效率，该积分考虑了所有波动频率模式的贡献。一个标准的替代方法是通过Wick旋转进行等高线变形，这涉及到沿着虚频轴评估Green张量。这种变换使被积函数平滑，将其转变为指数衰减的函数[63]。尽管频率域方法是Casimir型计算最早的数值方法之一[40,64,65]，但在许多离散频率上求解麦克斯韦方程会使得这种方法计算成本高昂。作为替代方案，时域方法可以从单次脉冲激励中提取完整的频谱响应[41,42]。这种效率促使我们选择通过时域求解麦克斯韦方程来评估Casimir相互作用。根据参考文献[41,42]，方程（6）中应力张量的期望值可以表示为时间卷积，(10)，其中表示从偶极激励开始经过的时间，下限积分是由因果关系决定的，这要求系统的响应在时间趋于零时消失，标记为激励开始的时间。在提出的框架内，平均麦克斯韦应力张量的评估简化为计算时间核函数和计算散射的电磁响应。前者可以通过解析方法处理，而后者则是通过我们自己实现的DGTD方法数值求解时变麦克斯韦方程得到的，参见图2。值得注意的是，这种数值方案是通用的，可以应用于其他时域电磁场求解器[66]。此外，散射响应与温度无关，因此只需要修改核函数就可以在后处理步骤中探索不同的温度，而无需额外的计算成本。图2. 计算工作流程的示意图。解析计算出的时域核函数和数值计算出的散射响应通过时间卷积结合在一起，然后在对封闭表面进行积分以得到Casimir力。1. 源电流密度和时间核函数核函数和散射响应的行为都与源电流密度的时间分布有关。在上述方法中，只要电流密度是因果的、连续的，并且产生的场是无散度的，这在形式上是任意的。因果性和物理可实现性要求电流密度为实数，并且在时间趋于零时消失。为了确保与麦克斯韦旋度方程和散度为零的演化一致，电流密度及其前两个导数必须在时间趋于零时消失[50,56]。在实际应用中，会选择特定的电流分布形式来加速方程（10）中时间积分的收敛。在参考文献[50]的基础上进行扩展，我们找到了一个指数衰减的多项式电流形式（11），其中整数确保电流及其前两个时间导数在时间趋于零时消失。振幅由决定，而因果性通过Heaviside阶跃函数来实现。参数控制衰减率，从而控制激励的频谱带宽，而离散参数（12）区分了电偶极（）和磁偶极（）激励，修改脉冲形状。实际上，选择（11）中的和可以确保在零温度和有限温度下，电偶极和磁偶极激励都能适当收敛，参见附录B。与参考文献[50]不同，在那里只有电场是相关的，由于存在长寿命的非振荡模式[67–69]，散射磁场的缓慢衰减会影响时域积分的收敛，因为它在长时间内会有不可忽视的贡献。核函数（9）取决于傅里叶变换后的电流密度，其表达式为（13）。通过将积分轮廓变形到复频平面，对方程（9）中的核函数进行傅里叶变换，可以得到（50）。由于因果性，核函数在上半平面是解析的，并满足Schwarz反射原理，这是任何实值时域函数傅里叶变换的要求。因此，当在上半平面闭合轮廓时，唯一对方程（9）有贡献的具有单极点的项是热因子（）。相应的极点位于Matsubara频率（）。应用残差定理得到（14），其中撇号表示该项乘以一个因子。对于方程（13）中的电流密度，方程（14）中的Matsubara求和可以表示为（15），其中和是对数函数的第二阶[70]。2. DGTD方法为了数值求解任意散射配置的麦克斯韦方程，我们采用了Hesthaven和Warburton引入的节点DGTD方法[71]。这种基于有限元的方法是为时域积分偏微分方程而设计的[72]。在这种方法中，计算域使用由有限元组成的非结构化网格离散化，这些有限元完全覆盖了模拟区域。在每个元素内部，通过用Lagrange多项式的一组局部基函数展开电磁场来局部求解麦克斯韦旋度方程，这些基函数的阶数为，展开系数与相应节点处的场振幅有关[51]。这种局部展开提高了内存效率和并行性[51,71,72]，从而减少了大规模三维模拟所需的计算资源。然后通过在元素接口引入所谓的上风数值通量来构建全局解，确保元素间的耦合和边界条件的满足[71]。通过空间离散化，麦克斯韦方程被转换为半离散形式，即一组一阶耦合的常微分方程，用于时间依赖的展开系数。这些方程使用显式的四阶低存储Runge-Kutta方案进行时间积分，共有14个阶段[73]，从而实现灵活的高阶时间步进技术。在数值时域方法中，由于频率依赖的介电常数不能像在频域那样通过点乘有效地纳入麦克斯韦方程；这样的实现需要在时间内进行计算成本高昂的卷积。相反，引入了一个辅助电流密度来表示材料的极化响应，该辅助电流密度满足从频域本构关系导出的辅助微分方程（ADE）。这个ADE与麦克斯韦方程在同一时间步进方案中同时求解[51,56]。在同一框架内，采用总场-散射场（TF-SF）技术来注入偶极激励产生的场分布，并将散射响应与总场分离。这种方法利用麦克斯韦方程的线性，通过将计算域划分为总场（TF）区域和散射场（SF）区域来分解电磁场，在TF区域中同时演化入射场和散射场，在SF区域中仅保留散射场。通过在TF-SF接口处的数值通量将入射场注入模拟中，通过添加或减去解析计算的源诱导场分布[51,56]。DGTD方法特别适用于这种形式，因为场值在元素接口的两侧是独立存储的。值得注意的是，当源放置在SF区域时，局部点源的奇异性不存在，这与之前介绍的散射Green张量形式一致。在建立了数值公式后，我们现在检验其计算效率和收敛特性。数值方法的性能通常由三个主要因素决定：（i）获得解的计算时间，（ii）相关的内存需求，以及（iii）所得解的准确性。对于显式的DGTD方案，计算复杂度随节点数（）和时间步数（）的平方增长，其中是由最小元素的半径决定的[51]，类似于FDTD方法中的Courant-Friedrichs-Lewy稳定性条件[56]。对于具有均匀网格细化的三维离散化，特征元素大小随（）增长，这意味着（）。因此，对于固定的物理模拟时间，这导致（），从而整个计算时间的增长为（）。对于局部细化的网格，如在Casimir力计算中通常使用的网格，DGTD方法可以支持高效的局部时间步进算法，允许在网格的较粗区域使用更大的时间步长[74]。这可以显著减少计算时间。相比之下，内存需求主要由每个元素内场未知数的存储和与相邻元素的固定数量局部通量相互作用决定，导致（）的渐近增长，与其他时域方法一致[56]。该方案的空间收敛性取决于特征边长（）和多项式阶数（）。增加（）可以在每个元素内添加更多的插值节点，从而提高空间分辨率。在这项工作中，我们使用直边四面体元素，这些元素特别适合于解析多尺度特征，并且比规则网格更自然地处理弯曲几何形状，减少了阶梯效应。当几何形状被网格很好地表示，即网格与几何形状一致时，该方法可以实现高阶收敛[51]。然而，对于具有未解决的曲率的非协调网格，收敛可能会受到限制[75]。对于平滑变化的边界，可以使用曲线元素来恢复高阶收敛[51,75]。IV. 应用在本节中，我们展示了所提方案的数值结果，重点研究了两个真空分离的金属物体之间的Casimir力。该程序中的一个关键步骤是在方程（1）中对封闭表面进行积分。与不需要数值应力张量积分的边界元素方法不同（39），当前方法通过离散化表面积分来评估力。原则上，这需要在对封闭表面的多个点评估应力张量，导致大量的时域模拟——每个评估点最多进行六次模拟，对应于三种正交的偶极取向和两种源类型。尽管这些模拟可以很容易地并行化，但高空间分辨率会显著增加所需的计算资源。虽然（）的几何形状在理论上可以是任意的，但某些选择在数值上可能比其他选择更高效。例如，如果一个物体被一个立方体包围，则积分的数值评估会大大简化。在无限扩展的立方体的情况下，唯一有贡献的表面是两个物体之间的面（），因为剩余表面上散射场的振幅随距离迅速衰减，在无穷远处它们的贡献消失。因此，只需在无限扩展的表面（）上进行积分，就可以完全捕捉到基本物理过程。进一步的简化可以依赖于特定系统的对称性。作为一个基准，我们首先考虑由真空间隙宽度（）分隔的金属半空间之间的Casimir相互作用。这种设置提供了对方程（10）中时间积分的洞察。此外，这种几何形状提供了计算上的 simplicity，并且允许基于Lifshitz理论（2）进行半解析处理，这为评估数值方法的准确性和收敛性提供了独立的参考计算。作为第二个例子，我们考虑一个位于半空间上方的有限大小的金属圆柱体，其轴垂直于表面。这种不太对称的配置突出了将应力张量与方程（1）中的力联系起来的表面积分方案。A. 平行半空间之间的Casimir力我们考虑两个由沿（）轴的距离（）分隔的相同金属半空间，参见图1(a)。金属的光学响应被建模为银的响应，由标准的Drude介电函数（16）描述，其等离子体频率为（），阻尼率为（）[76]。由于（）平面上的平移不变性，我们选择积分表面位于（）的中平面，与两个界面等距。这个表面在（）和（）方向上无限延伸，并且具有单位法向量（）。系统的对称性确保净Casimir力沿着（）轴方向，即（）。因此，应力张量的相关分量是（），需要在平面（）内的平面坐标上积分。然而，由于横向不变性，（）在（）上是空间均匀的。因此，单位面积的Casimir力（）可以在单个空间位置计算，即（17），从而将数值问题简化为单次时间积分。如上所述，系统的对称性允许进行半解析解，方程（17）可以根据Lifshitz公式（2,77）重新表示为（18），该公式描述了半空间上的压力，与平行于平面的电磁波矢量（）的分量有关，其模量为（），以及在虚Matsubara频率下评估的垂直于真空的波矢量分量（）。反射系数（），其中（），描述了半空间上场的散射。在当前背景下，这些结果由Fresnel公式[57,78]给出。如上所述，求和符号中的撇号表示该项具有权重。根据Lifshitz公式，这种相互作用可以通过其与两个特征长度尺度的相对大小来表征：等离子体波长λ_p，对于银来说大约为190 ?；以及热波长λ_T，大约对应于维恩位移定律中的特征长度尺度，在室温下约为280 ?。当λ_p远大于λ_T时，相互作用处于非延迟区域，在这种情况下可以忽略传播延迟。在这个极限下，方程(18)变为[78,79]：

F_f(r, θ) = (1/4π?_0) * (Δz/3) * [ζ(1/2 + (z/λ_p^2) * (z/λ_T)^2)]

其中ζ(z)是黎曼ζ函数，对于z ≥ 0有ζ(z) = (1/2) * (ln(z)^(1/2))。注意，上述表达式不依赖于系统的温度，这表明在这个区域，量子涨落主导了热涨落。当λ_p与λ_T接近时，延迟效应变得显著，在零温度极限下，Lifshitz公式恢复为Casimir的原始结果[1]：

F_C(r, θ) = -(π/16) * (Δz^2 / (4π?_0 * λ_p^4 * λ_T^4))

在这个极限下，力与材料参数无关，这些参数仅通过条件λ_p ≤ λ_T起作用。当λ_p远大于λ_T时，由于有限温度引起的修正变得明显。在λ_p远大于λ_T的极限下，经典热涨落主导了相互作用，Lifshitz公式趋于[79]中的形式。

上述表达式和相关量是我们数值计算的基准。

1. 时间积分
由于系统的对称性，我们在以原点O为中心的圆柱形计算域内进行了计算，该域的半径为R，高度为H，边界条件采用Silver–Müller边界条件[51]。金属半空间被近似为半径为r、厚度为d的有限圆盘，构成TF区域。由于金属的耗散性质，选择d的值使得来自相对界面的贡献可以忽略，从而有效恢复了Lifshitz最初考虑的两个相对半空间的配置[2]。选择r的值是为了最小化从计算域边缘传来的虚假反射[50]。在金属圆盘内，结合Drude模型的ADE（参见第III A.2节）求解麦克斯韦方程组。SF区域对应于圆盘之间的真空间隙，TF-SF边界与金属-真空界面重合。由Gmsh[5]生成的四面体网格在评估位置附近局部细化，平均边长为e，范围在半径为R的圆柱区域内。网格向计算域的边界线性加密，最外层处的边长达到e_max。我们使用全局多项式阶数N，确保至少每个最小相关波长有两个插值节点[51,56]。为了有效地处理麦克斯韦方程组，使用无量纲单位也是有利的。这种方法涉及用参考场强E_0和长度尺度λ_0对所有场和长度进行缩放。对于这里考虑的应用，我们使用E_0 = 1 V/m和λ_0 = 100 ?。

激发源位于x = 0处，并通过TF-SF技术引入模拟方案中，参见第III A.2节。方位不变性意味着在x = 0处，我们只需要区分平行于平面的方向（记为θ_c）和垂直于平面的方向（记为θ_perp）。在图3(a)中，我们展示了平行取向的电偶极子（左列）和垂直取向的磁偶极子（右列）的电流密度随时间的行为。脉冲轮廓对应于方程(11)中的定义，对于电偶极子和磁偶极子我们都使用相同的形式[50]。脉冲以β的速率衰减，β是一个可调参数，通常设置为一个适当的值，该值控制着频谱带宽，如附录B中所讨论的。

图3. （左列）描述了两个相隔d的银半空间的时间积分方案的相关量，这些半空间由平行于平面的电偶极子激发。所有量都用参考场强E_0和长度尺度λ_0进行缩放。每个面板展示了方案的不同组成部分：(a) 电流密度的时间轮廓，如方程(11)中定义的，以E_0和d为单位表示；(b) 在x = 0处的散射电场，以E_0为单位表示；(c) 以E_0为单位表示的核函数，附录A中显示了零温度下的渐近行为；(d) 时间积分函数，以E_0为单位表示，附录中关注了长时间行为。（右列）与左列相同，但适用于垂直磁偶极子。所有量都用以下单位表示：(a) ...（此处省略了具体的单位信息）。

在图3(b)中，我们在左右两列中展示了在节点点x = 0直接记录的相应散射场的时间演化，没有进行空间插值。由于配置的对称性，当x = 0时，散射场E_s和E_p消失。这个性质也被下面定义的相关量所继承。注意，在图3(b)中，当x = 0时，场幅度保持为零，反映了激发到达界面并返回源位置的往返时间。当x ≠ 0时，出现瞬态反射，随后在初始脉冲充分衰减后进入衰减阶段。在没有源区域的图中，图3(b)左侧面板中的电场分量以接近表面等离子体共振频率f_s的频率展示指数衰减振荡，这对应于非延迟电格林张量的主导极点[77]。另一方面，图3(b)右侧面板中的磁场幅度呈指数衰减，没有振荡行为。

为了评估方程(10)的时间收敛性，我们定义了单位面积上局部Casimir力的分量贡献为[22]：

F_c(r, θ) = ∑[F_f(r, θ) * Δz / (4π?_0 * λ_p^4 * λ_T^4)]

对方程(22)的数值评估需要在足够大的x值处截断上限，这在计算效率和数值精度之间引入了权衡。在参考文献[41]中，选择x值使得场幅度在人工耗散介质中降低到指定的误差阈值以下，其中耗散是可调的。然而，将这一标准扩展到物理系统可能会在计算上要求很高，并可能导致较长的模拟时间，特别是对于腔体和谐振器[80]。为了在保持精度的同时减少计算时间，我们选择只捕获基于系统物理特性的相关场信息[80,81]。然后通过在x = 0到t_max的区间内将有限持续时间信号拟合到指数衰减振荡模式的线性叠加来外推长时间响应，其中t_max选择为初始脉冲充分衰减的时间，例如在图3(b)中。为此，我们使用基于滤波器对角化方法[82,83]的时间域谐波反演技术，该技术具有提供长时间响应场的半解析表达式的优势。模态参数是使用Harminv实现在一个有限频率窗口内提取的[84]。在图3(b)的左右两面板中，我们展示了重构场和模拟场之间的极好一致性。因此，通过从模拟数据中提取的谐波反演参数，从x = 0积分方程(22)到无穷大，获得了半解析的时间修正[23]。

方程(23)右侧的第一个积分使用标准求积技术进行评估，其中使用模拟的散射场；第二个积分是解析计算的，其中ω_p是模式频率，φ_p是相位差，ζ_p是模式的衰减常数。直接数值积分与半解析长时间修正的结合加快了收敛速度，与仅基于方程(23)的第一项的方法相比。尽管谐波反演提供了长时间动力学的半解析描述，但它仍然是一种引入近似误差的数值技术。因此，方程(23)的准确性仍然取决于x的选择；因此，积分函数的足够快速衰减对于最小化误差至关重要。除了场耗散外，收敛性还受到方程(15)中定义的核函数k(x, t)的行为的支配，如图3(c)的左右两面板所示。如图3(d)所示，核函数在短时间内放大场贡献，而在较长时间内抑制这些贡献。值得注意的是，当t → ∞时，由于因果关系，散射场的发散是物理上规则的，因为散射场在t → ∞时消失。

将方程(23)的数值结果插入到方程(10)中的F_c分量中，得到沿x方向的单位面积Casimir力。图4显示了作为分离距离d函数的力的数值结果，归一化到Lifshitz公式(19)的渐近值。我们的数值方法（符号）与方程(18)中的半解析结果（实线）在广泛的分离距离和温度范围内非常吻合。当延迟可以忽略时（即x → ∞），数值数据趋近于单位值，对温度不敏感，并再现了方程(19)预测的预期渐近缩放。在较大的分离距离下，观察到两种不同的趋势。首先，当光速的有限性变得相关时[78,85]，力趋向于Casimir的原始表达式(20)，接近缩放关系F_c(r, θ) ∝ (δ/λ_p^4 * λ_T^4)。在这个区域，热涨落的贡献仍然足够小，尽管有可见的偏差。其次，当x → ∞时，热效应变得显著，并且根据方程(21)，数据逐渐恢复了预期的缩放关系。方程(21)中前因子的温度依赖性在图4中使用的归一化方案中清晰可见。

图4. 平行半空间之间的Casimir力作为它们分离距离d的函数。力归一化到非延迟的渐近值F_c(∞)，参见方程(19)。符号表示数值结果，而实线表示基于Lifshitz公式(18)的半解析参考计算。在短分离距离下，数值结果恢复了方程(19)预测的行为，并且与温度无关；而在大分离距离下，它们再现了方程(20)和(21)在零温度和有限温度下预期的依赖于d和T的依赖性。插图显示了使用Gmsh[5]生成的平面-平面几何的代表性网格。

2. 收敛性研究
为了量化数值方案的准确性和收敛性质，我们评估了单位面积上Casimir力的相对误差，分别为每个场分量计算[24]：

ΔF_c(r, θ) = ∑[F_f(r, θ) * Δz / (4π?_0 * λ_p^4 * λ_T^4)]

这里的F_f(r, θ)是从数值时间积分方案中获得的结果，如第IV A.1节所述，而ΔF_c表示通过复频域积分相应Green张量分量获得的参考值，参见附录C。

该方法的准确性受到不同误差源的限制，最显著的是任何基于网格的方法固有的有限空间离散化。为了评估其影响，我们进行了关于网格细化的收敛性研究，通过减小每个四面体元素的平均边长来统一减少网格细化，以及通过增加局部基函数的全局多项式阶数N来提高多项式细化[86]。

研究针对图4插图中的相同配置进行，并且结果是在22次名义上相同的模拟上平均得到的，这些模拟是通过在平面内移动x获得的不同的四面体分布[50]。这种平均减少了对于特定网格放置的敏感性，并提供了更稳健的空间离散化误差估计。图5展示了不同多项式阶数下的平均相对误差。图5(a)和5(b)中的电场贡献展示了预期的收敛行为[4]，与共形网格的理论预测一致，参见第III A.2节。相比之下，图5(c)和5(d)中磁场贡献的相对误差在超过某个多项式阶数和网格分辨率后趋于平稳，这表明过渡到了由时间截断误差而非空间离散化主导的区域，如下所述。在有限温度下，这种行为有所不同，例如在x → ∞时，时间核的更快衰减减少了截断误差，详见附录A。

图5. 在x = 0时，单位面积上分量Casimir力的相对误差，平均来自22次名义上相同的模拟，总持续时间为t_max。误差作为不同网格分辨率的多项式阶数的函数显示。网格大小定义为计算域中心区域的平均边长。结果针对零温度（蓝色）和室温（红色）分别展示。电磁场分量的单独贡献分别显示在面板(a)和(b)中：面板(a)和(b)分别对应于平行（δ）和垂直（δ_perp）电场分量，而面板(c)和(d)显示相应的磁场分量。

为了量化x对方程(23)中时间积分的影响，图6显示了在固定网格尺寸下，相对于x的相对误差作为分离距离d的函数，在零温度和有限温度下。在等式中（23），直接积分法与通过添加基于谐波反转的长时间贡献得到的结果进行了比较。对于电场分量，图6(a)和6(b)显示，随着的增加，误差会迅速趋近于由空间离散化决定的残差误差。由于在较大值时净偶极矩消失，因此对长时间贡献的依赖性较弱，导致等式（23）中的第二个积分在其渐近值周围振荡。相比之下，磁场分量，图6(c)和6(d)对时间截断更加敏感，使用谐波反转方法可以显著减少误差。进一步增加时间窗口可以通过更好地捕捉磁场分量的缓慢时间衰减特性来提高谐波反转的准确性[84]。随着的增加，直接积分和谐波反转的结果趋于相同的残差误差。中的局部最小值对应于误差的符号变化，并没有特殊意义。图6. 单位面积上各个分量的卡西米尔力的相对误差作为模拟时间的函数。误差分别针对、多项式阶数和网格尺寸进行展示。实线代表对时间域响应的直接积分，而虚线包括使用通过谐波反转获得的模态参数，从到进行的解析积分长时间贡献，根据等式（23）。结果显示了零温度（蓝色）和室温（红色）的情况。阴影区域表示22次名义上相同的模拟中相对误差的标准差，这些模拟的计算域足够大，以防止在模拟窗口内产生虚假的回射。图(a)和(b)分别展示平行（）和垂直（）电场分量，而图(c)和(d)展示相应的磁场分量。如参考文献[50]中所述，我们发现准确性可能还受到边界反射的限制。为了确认图5和图6中的残差误差不是由这种虚假回射主导的，我们使用更大的域进行控制模拟，在这些域中，从外边界到的距离超过了，从而防止边界反射在积分窗口内到达。图6中阴影区域表示的22次模拟的相对误差的标准差确认了残差误差主要是由空间离散化而非边界反射引起的。

3. 计算效率
在确定了数值误差的主要来源之后，我们现在评估所提出的DGTD框架的计算效率。为了将其放在适当的视角下，我们将其与用于任意形状几何体的卡西米尔力计算的成熟数值方法的规模进行了比较，例如在开源软件Scuff-em[39,55]中实现的BEM。使用一系列均匀细化的网格，我们检查了两种方法在前一节中考虑的平行半空间配置的计算时间、最大内存需求和相对误差。

图7比较了基于四面体元素的多项式阶数为和3的体积DGTD离散化与BEMs采用的表面-only三角形离散化。为了考虑平行半空间的对称性，考虑了两种BEM公式。通常，Scuff-em将任意几何体建模为具有封闭表面的紧凑对象。对于当前的几何体，这意味着将半空间近似为足够大的有限圆盘，这与第IV A 1节中使用的有限计算域一致。或者，可以通过周期性边界条件利用系统的平移不变性，得到无限扩展半空间配置的精确数值表示。虽然对于这种几何体，扩展对象的公式是更自然的选择，但包括紧凑对象的公式是为了说明通用BEM实现对于有限对象的计算规模。在两种BEM公式中，通过评估多个虚频下的电磁响应来进行卡西米尔力的计算以执行轮廓积分。相比之下，当前配置的对称性允许DGTD方法仅通过四个独立的时间域模拟来计算单位面积的力，对应于两种偶极方向（平行和垂直）和两种源类型（电和磁），参见第IV A 1节。这些模拟的独立性使它们易于并行化；因此，图7中报告的经过时间和内存需求对应于最苛刻的单个模拟[87]。

图7. 比较了多项式阶数为和3的DGTD方法与Scuff-em中实现的BEM对于紧凑和扩展对象的计算性能[55]。结果显示了由真空间隙分开的平行半空间的单位面积上的卡西米尔力。图(a)和(b)显示了单个时间或频率域模拟的运行时间和内存需求作为网格元素数量的函数。图(c)显示了相对误差作为平均边长函数的函数。两种数值方法的计算行为可以通过它们的渐近缩放属性来理解。如第III A 2节所讨论的，显式DGTD方案的计算成本按照缩放，其中表示网格元素的数量。相比之下，Scuff-em的成本主要由涉及密集交互矩阵的操作主导，当使用密集的直接求解器时，对于足够大的，其渐近缩放为[39]。这种行为在图7(a)中得到了体现，其中Scuff-em的计算时间在均匀网格细化下增加得比DGTD更快。图7(b)中显示的内存需求也观察到了类似的行为。对于DGTD方法，内存需求按照缩放，参见第III A 2节。相比之下，BEM需要存储所有表面元素对之间的交互，导致系统矩阵密集，内存缩放为。

B. 圆柱对称物体位于半空间上方
我们现在考虑一个有限的、圆柱对称的物体，它位于半空间上方，两者之间由宽度为的真空间隙分隔。坐标系的选择使得轴与系统的对称轴重合，并且垂直于占据的基底。由于围绕轴的圆柱对称性，卡西米尔力的侧向分量相互抵消，导致产生一个沿着方向的净力。

再次，我们选择一个底面位于的封闭表面，该表面位于物体和基底之间的中点，而其余表面延伸到无限远，由于散射场的衰减，其贡献可以忽略不计，参见图8。因此，等式（1）中的积分简化为在处的有效开放表面积分，单位法向量为。利用系统的轴对称性，我们使用圆柱坐标。方位角依赖性可以解析积分，将应力张量的表面积分简化为一个关于径向坐标的一维积分，可以使用高斯求积方案进行评估[89]。尽管高斯-拉盖尔求积原则上可以处理半无限区间，但实际上，它引入了远离物体的节点，这些节点对最终值的贡献可以忽略不计，同时将计算域扩展到不切实际的大小。将积分截断为有限区间并使用高斯-勒让德求积方案可以实现高效计算[90]。

图8. 半空间上方不同间距处的圆柱体的麦克斯韦应力张量的空间分布。值归一化为平面-平面麦克斯韦应力张量。我们利用有限元方法的自适应网格能力，评估了对应于高斯-勒让德求积阶数为[36,77]的非均匀径向位置的平滑变化的应力张量，该求积由最大截止半径界定，超过该半径后假设应力张量可以忽略不计。这种公式使得使用标准的高斯求积程序在适当的变量变换下能够进行高效的表面积分，同时只需要少量的评估点。然后，将结果力近似为（25），其中是离散化的径向坐标，和分别是相应的高斯-勒让德权重和节点。

作为示例，我们考虑一个半径为、高度为的有限金属圆柱体，它位于半空间上方。材料响应由方程（16）中给出的德鲁德形式描述。计算域如图1(b)所示，网格是根据第IV A 1节中概述的程序生成的。为了解决圆柱体的曲率问题，在靠近侧面的地方将网格尺寸设置为。在域内，通过16点高斯-勒让德求积在径向位置评估局部应力张量贡献。在每个求积分点放置一个点偶极子，并随时间演化相应的散射场。在偶极子位置记录场幅度，并使用第IV A 1节中描述的方案随时间积分。图8中显示的总麦克斯韦应力张量是通过将各个场贡献相加得到的，如方程（10）所给出的。我们注意到，在圆柱体中心附近，配置类似于前一节讨论的平移不变平面-平面几何体，因此应力张量几乎均匀。对于，会出现一个峰值，这可以归因于来自圆柱体侧面的散射场的增加贡献。对于，应力张量以的倒数比例指数级衰减，从而事后证明了的选择是合理的。这些贡献的空间积分得到了最终的卡西米尔力。鉴于高斯求积技术的快速收敛[89]，对于足够大的，我们预期主要的数值误差将来自网格的空间离散化。对于非共形网格，如本例中所示，收敛速率可能受到限制，参见第III A 2节。

DGTD和BEM方法具有互补的优势和局限性；因此，这里提出的方法不应被视为现有数值方法的替代品，而应视为一个额外的补充框架。BEM实现，如Scuff-em，已成为卡西米尔力计算的标准和高效工具，因为它们只需要离散化材料界面，从而导致未知数的数量大大减少。当相关的物理长度尺度与系统的特征几何特征相当时，这种优势尤为明显，这时可以用相对粗糙的表面网格准确表示系统。相比之下，DGTD方法需要对整个计算域进行体积离散化，从而导致给定几何体的自由度数量增加。然而，其在网格细化下的更有利的缩放使其在需要精细空间分辨率时具有吸引力，因为其计算时间和内存需求增加的速度比密集矩阵BEM慢。圆柱体的半径为r，高度为h。作用力被归一化到方程(19)中Casimir压力的非延迟极限值，并乘以圆柱体底面积。符号表示数值结果，虚线代表半解析渐近预测。在短距离分离时，数值结果能够恢复方程(19)预测的行为，并且与温度无关；而在大距离分离时，它们再现了从Casimir-Polder相互作用的渐近表达式方程(28)和(27)在零温和有限温度下预期的依赖性。插图展示了使用Gmsh [5]生成的圆柱体-平面几何结构的代表性网格。

当距离大于紧凑物体的半径时，PFA的准确性下降，最终失去有效性，数值结果也与半解析估计值偏离，参见图10。在这个范围内，相互作用越来越多地受到圆柱体多极散射贡献的控制[93]。对于圆柱体尺寸与其与表面距离相比较小的情况，可以忽略高阶多极，偶极响应主导了散射的电磁辐射。在这种情况下，系统可以有效地描述为一个极化物体与表面的相互作用，Casimir力逐渐过渡到Casimir-Polder极限[4]。与方程(18)中的Lifshitz表达式不同，Casimir-Polder相互作用直接与物体的极化率相关。虽然对于小圆柱体的极化率没有封闭形式的表达式[94]，但可以通过将其与等体积椭球体的极化率关联来提供一个可靠的近似[49]。在这种方法和特定情况下，一个体积为V的圆柱体可以大致用一个等体积的球体来近似。在这些条件下，频率依赖的极化率遵循Clausius-Mossotti关系[57,95]（26），其中ε_d是方程(16)中给出的介电函数，因此极化率在ε_d=0时具有共振特征。

图10显示了圆柱体与半空间之间的Casimir力随它们之间距离的变化。温度固定为T_0，同时考虑了不同的圆柱体半径r和相应的高度h。作用力被归一化到方程(20)中Casimir压力的延迟极限值，并乘以圆柱体底面积。与图9一样，符号表示DGTD计算，虚线代表渐近预测。此外，实线描述了基于PFA的半解析近似，而点划线对应于Scuff-em计算。在短距离分离时，数值结果恢复方程(19)预测的行为；而在大距离分离时，它们再现了方程(28)中的Casimir-Polder相互作用的依赖性，其中F(x, y)随x和y的变化而变化。

前面的考虑使我们能够估计大距离分离时力的行为。在有限温度下，即T≠0的情况下，如上所述，Casimir-Polder相互作用受到热波动的主导，产生热极限（27）。物体的有限横向范围导致一个幂律衰减，其衰减速度比相应极限下的Casimir力快一个阶。这在图9中的数值数据趋势中可以明显看到，与图4中观察到的情况相反。随着距离的增加，数据趋近于方程(27)中的渐近极限，该极限由系统的温度决定，这验证了数值结果和背后的物理推理。在零温度下，大距离分离时的相互作用具有不同的渐近行为。Casimir-Polder力的延迟极限由T=0定义。因此，对于r≠0和h≠0的情况，相互作用遵循（85,95,96）（28），这也表现出一个特征性的幂律衰减，其衰减速度比零温度下Casimir力快一个阶，参见方程(20)。我们在图10中观察到了这种行为，其中我们考虑了在不同r值下的圆柱体。为了突出与PFA预测（实线）的偏差，归一化是相对于方程(20)中的Casimir压力乘以圆柱体底面积来选择的。对于每个圆柱体尺寸，使用DGTD方法（符号）和在Scuff-em [39,55]中实现的BEM（点划线）获得的数值结果逐渐接近方程(28)描述的渐近行为。正如预期的那样，r的值越小，数值开始偏离PFA的距离就越短。渐近表达式对r的依赖性直接源于关系（28）。

最后，值得强调的是，交叉区域（r<
**结论**
在这项工作中，我们提出了一个基于时间域有限元的高精度框架，用于计算真空分隔物体之间的casimir力，该框架采用maxwell应力张量形式主义。数值方案采用不连续galerkin时间域方法来评估由通用材料模型描述的任意形状物体之间的有限温度casimir相互作用。通过将电磁green张量表达为系统对偶极激发的响应，应力张量被重新构建为一系列由电偶极源和磁偶极源在时间域驱动的经典散射问题。对于真空分隔的半空间之间的casimir相互作用的基准计算，在多个长度尺度上与lifshitz理论有极好的一致性。实现这一结果的相关步骤包括为电和磁源设计特定的时间轮廓以及开发基于谐波反演的程序来外推电磁场的长时行为。关于收敛性，我们的分析表明，准确性主要受时间截断和空间离散化的限制，而计算成本和内存需求分别随r和h的增加而增加，在均匀网格细化下。

我们进一步将这一框架应用于金属圆柱体与嵌入真空中的半空间之间的casimir相互作用，这种几何结构缺乏半解析表达式。数值结果在短距离和长距离极限下都恢复了预期的渐近行为，同时在中间范围内与成熟的数值方法也显示出良好的一致性，在这个范围内，几何效应、色散效应和热效应以非平凡的方式相互作用。我们方法的主要技术挑战之一是在封闭积分表面上高效评估maxwell应力张量通量。对于这项工作中的特定配置，我们使用了高斯求积方案来达到所需的精度水平，而所需的求积点数量相对较少。尽管目前的计算利用了圆柱对称性，但可以通过优化的二维高斯求积或积分方案[97–99]直接扩展到更一般的几何形状。尽管这样的扩展会随着求积点数量的增加而线性增加所需的模拟次数，但每次运行的独立性使得该方法非常适合在高性能计算集群上进行并行化。因此，只要有足够的计算资源，墙钟时间就由最具挑战性的单个dgtd模拟的成本决定，而不是由总的求积点数量决定。这突出了dgtd方法在研究非平凡casimir现象（包括几何诱导的排斥力和非单调距离依赖性[39,40,44–46]）时的灵活性，这些现象超出了常用近似方法（如邻近力近似和更严格的多极散射方法）的范围。尽管目前的实现仅限于真空分隔的物体，但底层框架实际上也适用于无色散介质辅助的配置，对于这些配置已经提出了替代的电磁应力张量公式[100,101]。对于色散介质和基于abraham的公式[102]，方程（9）中的核通过材料响应函数ε_d和μ_d获得了额外的频率依赖性。这种依赖性可以按照参考文献[41]中的方法进行数值处理，或者类似于参考文献[50]中对极化率的处理进行解析处理。此外，在存在色散背景介质的情况下，数值实现还需要在tf-sf边界重建入射偶极场[56]。

总之，dgtd方法为时间域中计算casimir型力提供了一个准确且灵活的替代方案。它能够适应任意几何形状、色散材料模型和有限温度，特别适合研究微观和纳米尺度系统中的波动诱导力。未来的发展可能包括实现曲线元素以保持高阶收敛性[75]，这对于实验设置中常见的球形几何形状尤其有益。在材料方面，该方法可以扩展到包含空间色散模型的纳米结构几何形状，这在casimir物理学中变得越来越重要[103,104]。dgtd方法可以通过ade公式[50,105]适应空间色散模型，为多尺度问题提供了重要优势[106]。相比之下，边界单元方法需要额外的公式来结合非局部体块响应[107]，而有限差分时域方法可能需要极其精细的结构化网格来同时解决亚纳米级的非局部长度尺度和较大的几何特征[108]。 **结论** 在这项工作中，我们提出了一个基于时间域有限元的高精度框架，用于计算真空分隔物体之间的casimir力，该框架采用maxwell应力张量形式主义。数值方案采用不连续galerkin时间域方法来评估由通用材料模型描述的任意形状物体之间的有限温度casimir相互作用。通过将电磁green张量表达为系统对偶极激发的响应，应力张量被重新构建为一系列由电偶极源和磁偶极源在时间域驱动的经典散射问题。对于真空分隔的半空间之间的casimir相互作用的基准计算，在多个长度尺度上与lifshitz理论有极好的一致性。实现这一结果的相关步骤包括为电和磁源设计特定的时间轮廓以及开发基于谐波反演的程序来外推电磁场的长时行为。关于收敛性，我们的分析表明，准确性主要受时间截断和空间离散化的限制，而计算成本和内存需求分别随r和h的增加而增加，在均匀网格细化下。 我们进一步将这一框架应用于金属圆柱体与嵌入真空中的半空间之间的casimir相互作用，这种几何结构缺乏半解析表达式。数值结果在短距离和长距离极限下都恢复了预期的渐近行为，同时在中间范围内与成熟的数值方法也显示出良好的一致性，在这个范围内，几何效应、色散效应和热效应以非平凡的方式相互作用。我们方法的主要技术挑战之一是在封闭积分表面上高效评估maxwell应力张量通量。对于这项工作中的特定配置，我们使用了高斯求积方案来达到所需的精度水平，而所需的求积点数量相对较少。尽管目前的计算利用了圆柱对称性，但可以通过优化的二维高斯求积或积分方案[97–99]直接扩展到更一般的几何形状。尽管这样的扩展会随着求积点数量的增加而线性增加所需的模拟次数，但每次运行的独立性使得该方法非常适合在高性能计算集群上进行并行化。因此，只要有足够的计算资源，墙钟时间就由最具挑战性的单个dgtd模拟的成本决定，而不是由总的求积点数量决定。这突出了dgtd方法在研究非平凡casimir现象（包括几何诱导的排斥力和非单调距离依赖性[39,40,44–46]）时的灵活性，这些现象超出了常用近似方法（如邻近力近似和更严格的多极散射方法）的范围。尽管目前的实现仅限于真空分隔的物体，但底层框架实际上也适用于无色散介质辅助的配置，对于这些配置已经提出了替代的电磁应力张量公式[100,101]。对于色散介质和基于abraham的公式[102]，方程（9）中的核通过材料响应函数ε_d和μ_d获得了额外的频率依赖性。这种依赖性可以按照参考文献[41]中的方法进行数值处理，或者类似于参考文献[50]中对极化率的处理进行解析处理。此外，在存在色散背景介质的情况下，数值实现还需要在tf-sf边界重建入射偶极场[56]。>
**结论**
在这项工作中，我们提出了一个基于时间域有限元的高精度框架，用于计算真空分隔物体之间的casimir力，该框架采用maxwell应力张量形式主义。数值方案采用不连续galerkin时间域方法来评估由通用材料模型描述的任意形状物体之间的有限温度casimir相互作用。通过将电磁green张量表达为系统对偶极激发的响应，应力张量被重新构建为一系列由电偶极源和磁偶极源在时间域驱动的经典散射问题。对于真空分隔的半空间之间的casimir相互作用的基准计算，在多个长度尺度上与lifshitz理论有极好的一致性。实现这一结果的相关步骤包括为电和磁源设计特定的时间轮廓以及开发基于谐波反演的程序来外推电磁场的长时行为。关于收敛性，我们的分析表明，准确性主要受时间截断和空间离散化的限制，而计算成本和内存需求分别随r和h的增加而增加，在均匀网格细化下。

我们进一步将这一框架应用于金属圆柱体与嵌入真空中的半空间之间的casimir相互作用，这种几何结构缺乏半解析表达式。数值结果在短距离和长距离极限下都恢复了预期的渐近行为，同时在中间范围内与成熟的数值方法也显示出良好的一致性，在这个范围内，几何效应、色散效应和热效应以非平凡的方式相互作用。我们方法的主要技术挑战之一是在封闭积分表面上高效评估maxwell应力张量通量。对于这项工作中的特定配置，我们使用了高斯求积方案来达到所需的精度水平，而所需的求积点数量相对较少。尽管目前的计算利用了圆柱对称性，但可以通过优化的二维高斯求积或积分方案[97–99]直接扩展到更一般的几何形状。尽管这样的扩展会随着求积点数量的增加而线性增加所需的模拟次数，但每次运行的独立性使得该方法非常适合在高性能计算集群上进行并行化。因此，只要有足够的计算资源，墙钟时间就由最具挑战性的单个dgtd模拟的成本决定，而不是由总的求积点数量决定。这突出了dgtd方法在研究非平凡casimir现象（包括几何诱导的排斥力和非单调距离依赖性[39,40,44–46]）时的灵活性，这些现象超出了常用近似方法（如邻近力近似和更严格的多极散射方法）的范围。尽管目前的实现仅限于真空分隔的物体，但底层框架实际上也适用于无色散介质辅助的配置，对于这些配置已经提出了替代的电磁应力张量公式[100,101]。对于色散介质和基于abraham的公式[102]，方程（9）中的核通过材料响应函数ε_d和μ_d获得了额外的频率依赖性。这种依赖性可以按照参考文献[41]中的方法进行数值处理，或者类似于参考文献[50]中对极化率的处理进行解析处理。此外，在存在色散背景介质的情况下，数值实现还需要在tf-sf边界重建入射偶极场[56]。

总之，dgtd方法为时间域中计算casimir型力提供了一个准确且灵活的替代方案。它能够适应任意几何形状、色散材料模型和有限温度，特别适合研究微观和纳米尺度系统中的波动诱导力。未来的发展可能包括实现曲线元素以保持高阶收敛性[75]，这对于实验设置中常见的球形几何形状尤其有益。在材料方面，该方法可以扩展到包含空间色散模型的纳米结构几何形状，这在casimir物理学中变得越来越重要[103,104]。dgtd方法可以通过ade公式[50,105]适应空间色散模型，为多尺度问题提供了重要优势[106]。相比之下，边界单元方法需要额外的公式来结合非局部体块响应[107]，而有限差分时域方法可能需要极其精细的结构化网格来同时解决亚纳米级的非局部长度尺度和较大的几何特征[108]。>