一种基于电荷的蒙特卡洛粒子合并方法，用于单元格内粒子模拟

《Computer Physics Communications》：A Monte Carlo particle merge method based on the charge for particle-in-cell simulations

　　摘要：本文提出了一种基于电荷的、极其简单且计算效率高的蒙特卡洛粒子合并方法，用于粒子在单元格（PIC）模拟中，该方法产生的数值扰动最小。该算法在有效控制宏观粒子数量的同时，能够在统计意义上保持关键物理量的守恒，并维持相空间结构。与传统方法不同，这种新方法不需要粒子的位置和速度信息，也避免了复杂的粒子间数据依赖关系，因此与经典PIC方法的粒子级并行性天然兼容。多项数值基准测试及实际应用案例表明，该方法在物理精度与计算效率之间取得了良好平衡。

引言：在经典的多粒子系统中，粒子的加速度由其位置和速度决定。因此，具有相似位置和速度信息的多个粒子可以被合并为一个宏观粒子。通过追踪多个宏观粒子的位置和速度，可以预测粒子位置-速度分布函数的演化。一个成功的典型例子就是将麦克斯韦方程与粒子运动方程相结合的粒子在单元格模拟方法，该方法已被广泛用于等离子体物理问题的模拟[[1], [2], [3], [4], [5], [6], [7]]。然而，在一些实际应用中（如多极器放电和气体击穿）[[8], [9], [10], [11], [12], [13], [14], [15], [16], [17], [18]]，系统中带电粒子的数量会随时间呈指数增长。只有通过在模拟过程中不断合并宏观粒子来控制其数量，才能将计算负荷降至可接受的水平。

一种显而易见的合并方法是把位置和速度相似的宏观粒子合并在一起。但这种方法会带来两个典型问题：其一，寻找位置和速度相似的粒子通常需要大量计算资源；其二，当模拟在速度-位置相空间中存在重叠且传播方向相似的多个粒子流时，宏观粒子的合并可能会导致原本相邻但差异较大的多个峰被错误地简化为一个个峰。为了解决这些问题，一些学者尝试了一些不直接合并位置和速度相似粒子的解决方案。

在参考文献[19]中，采用了一种统计方法，直接对速度-位置空间中的密度分布函数进行重采样，从而使用更少的宏观粒子。虽然这种方法能够很好地保持多种物理量以及局部意义上的速度分布函数的守恒，但在粒子合并过程中需要读取局部粒子群中所有粒子的位置、速度和权重信息，同时还依赖于粒子的速度和位置信息来重建局部的速度-位置相空间分布。因此，在处理当前粒子时，该算法不仅取决于当前粒子本身的信息，还依赖于局部粒子群中其他粒子的信息。另一种不同的尝试方法是参考文献[20]中提出的不可知保守降采样方法。该方法能够完美解决在相空间中重叠且传播方向相似的粒子流的合并问题；此外，该算法还能确保任意物理量的守恒。不过，要实现严格的守恒，需要从局部粒子群中提取信息，构建与需要守恒的物理量对应的约束向量，然后生成与这些约束向量正交的平衡方向向量，以此进行保守性权重调整。换言之，单个粒子的处理结果不再仅由其自身属性决定，还取决于同一粒子群内其他粒子的状态和约束作用。

因此，这两种方法在算法实现过程中都会引入复杂的粒子间数据依赖关系，从而降低了它们与经典PIC方法底层计算逻辑的兼容性。在经典PIC方法的粒子推动和粒子到网格加权过程中，处理单个宏观粒子时无需读取其他粒子的数据。此外，这种不可知保守降采样方法每轮合并操作只能减少一个宏观粒子。如果要减少多个宏观粒子，则需要进行多轮重复操作，这可能会涉及大量的计算工作。本文将继续基于蒙特卡洛框架研究不可知降采样合并算法，以解决上述两个典型问题。

粒子模拟方法本身是通过众多宏观粒子的运动状态统计来反映结果的。因此，在合并过程中并不一定需要严格保持物理量的守恒，只需在统计层面上确保物理量得到守恒即可。如果放弃严格守恒物理量的理念，就可以实现一种非常简单且计算效率高的蒙特卡洛宏观粒子合并方法（初步MC算法），该方法已在我们关于微波气体击穿的研究中得到成功应用[13]。这种算法在合并宏观粒子时无需知道它们的位置和速度信息，也不需要获取或等待其他宏观粒子的数据，它仅利用简单的蒙特卡洛方法来减少宏观粒子的数量，同时保持了经典PIC方法避免复杂粒子间依赖关系的特点。

本文通过研究统计波动，对该方法（初步MC算法）进行了分析和改进，从而得到了一种基于宏观粒子电荷信息、统计波动最小的改进型MC宏观粒子合并方法（改进MC算法）。这样一来，改进型MC算法对速度-位置空间中密度分布函数造成的数值扰动，在宏观粒子数量相同的情况下几乎与没有合并操作时的扰动相当，而且该算法在保持较低相对误差水平的同时，还能实现较低的计算成本。此外，通过使用CHIPIC软件进行PIC模拟可以发现，改进型MC合并算法不会影响进入Vaughan二次电子发射模型的初级电子的入射能量分布，也不会影响多极器放电和气体放电的宏观演化过程。

第2节介绍了初步MC算法和改进型MC算法，第3节则通过一些示例验证了改进型MC算法的有效性。

片段内容：

初步MC算法：用于粒子合并的蒙特卡洛算法的初步版本可见我们之前关于气体击穿的论文[13]。该算法如图1所示，其中N为粒子数量，Nmax为N的上限值，Pi为第i个粒子灭绝的概率。如果N>Nmax，则第i个粒子的电荷Qi会增加到Qi/(1-Pi)，并且该粒子有概率Pi= (NNmax)/N而灭绝。处理完所有粒子后……

一维速度分布：随后，我们对初步MC算法、改进MC算法、参考文献[20]中描述的不可知保守降采样算法以及不进行宏观粒子合并的情况进行了的一维数值模拟，并比较了这些模拟的结果。模拟设置如下：模拟时间为2000个时间步，每个时间步的长度为1×10?10秒。在每个时间步中，有500个电荷量为Q = 1 ×10?11C的宏观粒子。

结论：针对粒子在单元格模拟中由于宏观粒子数量呈指数增长而导致的计算负荷过大的问题，本文提出了一种简单的基于电荷的蒙特卡洛宏观粒子合并方法。该方法符合PIC模拟的统计特性，且与经典PIC方法的粒子级并行性兼容。该算法仅利用粒子自身的局部属性以及少量预先计算好的参数来评估合并概率……

代码获取地址：DOI: https://doi.org/10.5281/zenodo.21335035

CRediT作者贡献说明：
黄欣怡：写作——初稿撰写、可视化、验证、研究、正式分析。
王慧慧：写作——审稿与编辑、验证、监督、项目管理、方法论、研究、资金获取、正式分析、概念设计。
刘大刚：软件、资源支持。
刘拉群：验证工作。
李志杰：监督工作。

黄欣怡|王慧慧|刘大刚|刘拉群|李志杰
中国电子科技大学物理学院，中国成都市610054