NBody 宇宙学模拟
Published:
（此文档未完成，仍在编写中）
Zeldovich
PM 法
在宇宙学 N-Body 模拟中，我们需要求解数以百万计甚至百亿计的暗物质粒子在自身引力作用下的演化. 直接计算任意两个粒子之间的引力（Direct Summation, DS）需要 O(N²) 的计算量，这对于大规模模拟是不可接受的. 粒子-网格（Particle-Mesh, PM）方法是一种经典且高效的算法，它将长程引力的计算复杂度降低到 O(N log N)，极大地提升了模拟效率.
PM 方法的核心思想是：将粒子所代表的离散质量分布”绘制”到一个规则的计算网格上，然后在网格上利用快速傅里叶变换（Fast Fourier Transform, FFT）高效地求解引力势，再从引力势计算出引力场，最后将网格上的引力场插值回粒子所在的位置，从而驱动粒子的运动.
PM 方法的数学框架
整个宇宙的物质密度可以表示为 \(\rho(\boldsymbol{x},t)\)，它产生的引力势 \(\Phi(\boldsymbol{x},t)\) 由下面的泊松方程（Poisson’s Equation）描述：
\[\nabla^{2}\Phi(\boldsymbol{x},t)=4\pi G\rho(\boldsymbol{x},t)\]在宇宙学中，通常使用共动坐标（comoving coordinates），此时泊松方程会包含宇宙膨胀因子 \(a(t)\)，形式变为：
\[\nabla^{2}\Phi=4\pi Ga^{2}\bar{\rho}\delta\]其中 \(\bar{\rho}\) 是宇宙的平均密度，\(\delta=(\rho-\bar{\rho})/\bar{\rho}\) 是密度扰动（density contrast）. 我们的目标就是求解这个方程，得到引力势 \(\Phi\)，然后通过 \(\boldsymbol{g}=-\nabla\Phi\) 计算引力场，最终更新粒子的速度和位置.
PM 方法将这个求解过程分为以下四个主要步骤：
质量分配（Mass Assignment）：将粒子的质量分配到网格点上，得到网格化的密度场.
泊松方程求解（Solving Poisson’s Equation）：利用 FFT 在傅里叶空间中求解泊松方程，得到网格点上的引力势.
引力计算（Force Calculation）：通过对引力势求梯度，计算出网格点上的引力场.
力插值（Force Interpolation）：将网格点上的引力场插值到每个粒子的实际位置.
质量分配 (Mass Assignment)
第一步是将代表连续质量分布的 N 个离散粒子”云”的质量分配到一个三维的规则网格上. 假设我们的模拟盒子是一个边长为 L 的立方体，被划分为 \(M\times M\times M\) 个网格单元，每个网格单元的边长为 \(H=L/M\).
令 \(\boldsymbol{x}_{p}\) 为第 p 个粒子的位置，其质量为 \(m_{p}\). 网格点的密度 \(\rho_{i,j,k}\) 可以通过一个分配函数 \(W(\boldsymbol{x})\) 来计算. 这个函数将粒子位置 \(\boldsymbol{x}_{p}\) 与网格点索引 \((i,j,k)\) 关联起来.
\[\rho_{i,j,k}=\frac{1}{H^{3}}\sum_{p}m_{p}W(\boldsymbol{x}_{i,j,k}-\boldsymbol{x}_{p})\]其中 \(\boldsymbol{x}_{i,j,k}\) 是网格点 \((i,j,k)\) 的坐标. 最常用的分配方案为 CIC 也就是云中单元分配方案. 它将每个粒子看作一个与网格单元同样大小的均匀立方体”云”，然后根据该”云”与周围网格单元的重叠体积，将质量按比例分配给相邻的 8 个网格点（在三维情况下）. 这等价于对粒子位置进行线性插值.
对于一维情况，粒子在 \(x_{p}\) 处，其左右两个网格点为 \(x_{i}\) 和 \(x_{i+1}\)，则分配的权重为：
\[w_{i}=1-\frac{x_{p}-x_{i}}{H}\quad\text{and}\quad w_{i+1}=\frac{x_{p}-x_{i}}{H}\]在三维中，就是这三个方向权重的乘积.
泊松方程求解 (FFT)
获得网格上的密度场 \(\rho(\boldsymbol{x})\) 后（为了简洁，我们省略下标 \(i,j,k\)），下一步是求解泊松方程 \(\nabla^{2}\Phi=4\pi G\rho\). 直接在实空间求解这个偏微分方程很困难，但在傅里叶空间（或 k 空间）它会变成一个简单的代数方程.
首先，对密度场和引力势进行离散傅里叶变换（Discrete Fourier Transform, DFT）：
\[\hat{\rho}(\boldsymbol{k})=\sum_{\boldsymbol{x}}\rho(\boldsymbol{x})e^{-i\boldsymbol{k}\cdot\boldsymbol{x}}\] \[\hat{\Phi}(\boldsymbol{k})=\sum_{\boldsymbol{x}}\Phi(\boldsymbol{x})e^{-i\boldsymbol{k}\cdot\boldsymbol{x}}\]其中 \(\boldsymbol{k}\) 是傅里叶空间中的波矢. 在傅里叶空间，微分算子 \(\nabla\) 变成了代数乘子 \(i\boldsymbol{k}\). 因此，\(\nabla^{2}\) 对应于 \(-|\boldsymbol{k}|^{2}\). 泊松方程就变成了：
\[-|\boldsymbol{k}|^{2}\hat{\Phi}(\boldsymbol{k})=4\pi G\hat{\rho}(\boldsymbol{k})\]于是，我们可以直接解出引力势的傅里叶分量：
\[\hat{\Phi}(\boldsymbol{k})=-4\pi G\frac{\hat{\rho}(\boldsymbol{k})}{|\boldsymbol{k}|^{2}}\]这个代数方程的求解非常简单. 我们定义傅里叶空间中的格林函数（Green’s function） \(\hat{G}(\boldsymbol{k})\)：
\[\hat{G}(\boldsymbol{k})=-\frac{4\pi G}{|\boldsymbol{k}|^{2}}\]于是，\(\hat{\Phi}(\boldsymbol{k})=\hat{G}(\boldsymbol{k})\hat{\rho}(\boldsymbol{k})\).
注意：对于 \(\boldsymbol{k}=0\) 的模式（对应宇宙的平均密度），分母为零. 这个模式不产生引力，通常设 \(\hat{\Phi}(\boldsymbol{k}=0)=0\). 此外，由于网格的离散性，拉普拉斯算子 \(\nabla^{2}\) 和波矢 \(\boldsymbol{k}\) 实际上需要使用它们的有限差分近似形式，以获得更高的精度. 例如，波数分量 \(k_{x}\) 可以用 \(k_{x}=\frac{2}{H}\sin(\frac{\pi i}{M})\) 来代替简单的 \(i/H\).
实际计算中，我们使用快速傅里叶变换（FFT）算法来执行 DFT，其计算复杂度为 O(M log M)，其中 M 是网格点的总数（在 NM͂ 的情况下，复杂度为 O(N log N)）.
引力计算
得到了网格上的引力势 \(\Phi(\boldsymbol{x})\) 后，引力场 \(\boldsymbol{g}(\boldsymbol{x})\) 可以通过求其负梯度得到：
\[\boldsymbol{g}(\boldsymbol{x})=-\nabla\Phi(\boldsymbol{x})\]这个梯度运算同样可以在傅里叶空间或实空间进行. 因为 \(\nabla\leftrightarrow i\boldsymbol{k}\)，所以引力场的傅里叶分量为：
\[\hat{\boldsymbol{g}}(\boldsymbol{k})=-i\boldsymbol{k}\hat{\Phi}(\boldsymbol{k})\]我们可以在得到 \(\hat{\Phi}(\boldsymbol{k})\) 后，直接乘以 \(-i\boldsymbol{k}\)，然后进行一次逆 FFT 得到 \(\boldsymbol{g}(\boldsymbol{x})\).
同时，我们也可以通过有限差分得到，例如，x 方向的引力分量可以表示为：
\[g_{x}(i,j,k)=-\frac{\Phi(i+1,j,k)-\Phi(i-1,j,k)}{2H}\]当然，实际中可以使用更高精度的差分公式.
力插值 (Force Interpolation)
最后一步是将网格点上的引力场 \(\boldsymbol{g}(\boldsymbol{x}_{i,j,k})\) 插值到每个粒子所在的实际位置 \(\boldsymbol{x}_{p}\)，得到作用在每个粒子上的力 \(\boldsymbol{F}_{p}\).
\[\boldsymbol{F}_{p}=m_{p}\boldsymbol{g}(\boldsymbol{x}_{p})=m_{p}\sum_{i,j,k}\boldsymbol{g}(\boldsymbol{x}_{i,j,k})W(\boldsymbol{x}_{i,j,k}-\boldsymbol{x}_{p})\]最后，使用一个时间积分方案来更新粒子的速度和位置，然后进入下一个时间步，重复上述整个过程.
\[\boldsymbol{v}_{p}(t+\Delta t/2)=\boldsymbol{v}_{p}(t-\Delta t/2)+\boldsymbol{F}_{p}(t)\Delta t/m_{p}\] \[\boldsymbol{x}_{p}(t+\Delta t)=\boldsymbol{x}_{p}(t)+\boldsymbol{v}_{p}(t+\Delta t/2)\Delta t\]