NBody 宇宙学模拟
Published:
基本方程
在宇宙学中，一个重要的假设是：宇宙是均匀且各向同性的，当然，实际的宇宙并没有我们所假设的那么绝对均匀，在较小的尺度上，物质的涨落是不可忽略的（毕竟如果是一个完全均匀的宇宙也不会形成各种各样的天体结构），当然，我们的宇宙整体而言还是相当均匀的，因此，我们使用微扰来分析物质涨落，这里我们只需要讨论一阶微扰即可得到较为精确的结果，即 \(\rho\left(\boldsymbol{x},t\right)=\bar{\rho}\left(t\right)\left(1+\delta\left(\boldsymbol{x},t\right)\right)\)，因此有微扰场：
\[\delta\left(\boldsymbol{x},t\right)=\frac{\rho\left(\boldsymbol{x},t\right)-\bar{\rho}\left(t\right)}{\bar{\rho}\left(t\right)}\]其中，\(\bar{\rho}\) 是宇宙的平均物质密度.
Friedmann 方程
背景宇宙的哈勃参数 \(H\) 和宇宙的标度因子 \(a\) 之间的关系由 Friedmann 方程给出：
\[H\left(a\right)=H_{0}\sqrt{\Omega_{m0}a^{-3}+\Omega_{\Lambda0}}=H_{0}E\left(a\right)\]其中，\(\Omega_{m0}\) 为今天的宇宙物质密度，\(\Omega_{\Lambda0}\) 为今天的宇宙暗能量密度. 另外，由于 \(\dot{a}=aH\left(a\right)\)，因此，有：
\[\dot{a}=aH_{0}\sqrt{\Omega_{m0}a^{-3}+\Omega_{\Lambda0}}\]同时，减速参数由下式给出：
\[q=\frac{\Omega_{m0}a^{-3}-2\Omega_{\Lambda0}}{2E^{2}\left(a\right)}\]傅里叶变换
我们考虑一个边长为 \(l\) 的共动体积 \(V=l\times l\times l\)，为了避免”边界效应”（由于我们模拟的体积是有限的，如果边界外没有物质，则靠近边界的粒子会受到一个额外的指向内部的引力，导致和现实有所不同），因此我们一般使用周期性边界条件，也就是 \(\delta\left(x\right)=\delta\left(x+l\right)\) （三个方向都满足），因此，我们可以对微扰场做傅里叶变换，有：
\[\hat{\delta}\left(\boldsymbol{k}\right)=\int_{\in V}\mathrm{d}^{3}x\ \exp\left(-\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{x}\right)\delta\left(\boldsymbol{x}\right)\]实空间的周期性使得傅里叶空间是离散的，也就是 \(\boldsymbol{k}\) 的取值为：
\[\boldsymbol{k}=2\pi\left(\frac{n_{1}}{l}\hat{\boldsymbol{i}}+\frac{n_{2}}{l}\hat{\boldsymbol{j}}+\frac{n_{3}}{l}\hat{\boldsymbol{k}}\right)\]其中：
\[n_{1},n_{2},n_{3}=0,\pm1,\pm2,\pm3,\dots\]功率谱
功率谱 \(P\left(k\right)\) 定义为：
\[\left\langle \hat{\delta}\left(\boldsymbol{k}\right)\cdot\hat{\delta}^{*}\left(\boldsymbol{k}'\right)\right\rangle =V\delta_{\boldsymbol{k},\boldsymbol{k}'}P\left(k\right)\]这里 \(\left\langle \dots\right\rangle\) 表示对里面的内容做平均，\(\delta_{\boldsymbol{k},\boldsymbol{k}'}\) 为 Kronecker 函数. 不难发现，功率谱 \(P\left(k\right)\) 描述的是，在特定尺度 \(k\) 上，密度涨落的振幅的平方的平均值. 换句话说也就是 \(P\left(k\right)\) 描述了宇宙在不同尺度上的”成团”的程度.
另外，上述定义中假定了功率谱是各向同性的，也就是其与方向无关，\(P\left(\boldsymbol{k}\right)=P\left(k\right)\).
粒子动力学
由于宇宙是在膨胀的，因此两点之间的距离会随着时间的增加而增加，即使两点相对背景网格都没有任何移动. 为了方便讨论，我们引入共动坐标 \(\boldsymbol{x}\)，若粒子相对背景网格没有移动，则其共动坐标不会发生变化. 与之相对的则是本动坐标 \(\boldsymbol{r}\)，两个坐标之间由标度因子 \(\boldsymbol{a}\) 联系：
\[\boldsymbol{r}=a\boldsymbol{x}\]并且，我们定义今天宇宙的标度因子为 \(a=1\).
对于一个粒子，其在引力场 \(\varphi\left(\boldsymbol{r}\right)\) 中的拉氏量为（记粒子质量 \(m=1\)）：
\[\begin{aligned} \mathcal{L} & =\frac{1}{2}\dot{\boldsymbol{r}}^{2}-\varphi\left(\boldsymbol{r}\right)\\ & =\frac{1}{2}\left(\dot{a}^{2}\dot{\boldsymbol{x}}^{2}+a^{2}\dot{\boldsymbol{x}}^{2}+2a\dot{a}\boldsymbol{x}\dot{\boldsymbol{x}}\right)-\varphi\left(\boldsymbol{x}\right)\\ & =\frac{1}{2}a^{2}\dot{\boldsymbol{x}}^{2}-\frac{1}{2}a\ddot{a}\boldsymbol{x}^{2}-\varphi\left(\boldsymbol{x}\right)+\frac{\mathrm{d}}{\mathrm{dt}}\left(\frac{1}{2}a\dot{a}\boldsymbol{x}^{2}\right) \end{aligned}\]上式中最后一项略去，因此有：
\[\mathcal{L}=\frac{1}{2}a^{2}\dot{\boldsymbol{x}}^{2}-\phi\left(\boldsymbol{x}\right)\]其中：
\[\phi\left(\boldsymbol{x}\right)=\varphi\left(x\right)+\frac{1}{2}a\ddot{a}\boldsymbol{x}^{2}\]对于一个均匀且各项同性的宇宙，有
\[\ddot{a}\boldsymbol{x}=-\frac{4}{3}\pi r^{3}\bar{\rho}\frac{G}{r^{2}}\hat{\boldsymbol{r}}\]即：
\[\ddot{a}=-\frac{4}{3}\pi a\bar{\rho}G\]由于 \(\varphi\) 满足泊松方程：
\[\nabla_{x}^{2}\varphi\left(\boldsymbol{x}\right)=4\pi G\rho a^{2}\]因此有：
\[\begin{aligned} \nabla_{x}^{2}\phi\left(\boldsymbol{x}\right) & =\nabla_{x}^{2}\left(\varphi\left(\boldsymbol{x}\right)+\frac{1}{2}a\ddot{a}\boldsymbol{x}^{2}\right)\\ & =\nabla_{x}^{2}\varphi\left(\boldsymbol{x}\right)+\frac{1}{2}a\ddot{a}\nabla_{x}^{2}\boldsymbol{x}^{2}\\ & =4\pi G\rho a^{2}+3a\ddot{a}\\ & =4\pi Ga^{2}\left(\frac{\rho}{\bar{\rho}}-1\right)\\ & =4\pi Ga^{2}\bar{\rho}\delta\left(\boldsymbol{x}\right) \end{aligned}\]之前的拉格朗日量可以写出对应的哈密顿量：
\[\mathcal{H}=\frac{\boldsymbol{v}^{2}}{2a^{2}}+\phi\left(\boldsymbol{x}\right)\]其对应的哈密顿方程给出：
\[\begin{aligned} \frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}t} & =\frac{\partial\mathcal{H}}{\partial\boldsymbol{v}}=\frac{\boldsymbol{v}}{a^{2}}\\ \frac{\mathrm{d}\boldsymbol{v}}{\mathrm{d}t} & =-\frac{\partial\mathcal{H}}{\partial\boldsymbol{x}}=-\frac{\partial\phi\left(\boldsymbol{x}\right)}{\partial\boldsymbol{x}}=-\nabla\phi\left(\boldsymbol{x}\right) \end{aligned}\]在模拟中，直接使用 \(t\) 进行计算不是很方便，我们一般使用标度因子 \(a\) 为变量，因此 \(\mathrm{d}t=\mathrm{d}a/\dot{a}=\mathrm{d}a/aH\left(a\right)\)，于是有：
\[\begin{aligned} \frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}a} & =\frac{\boldsymbol{v}}{a^{3}H\left(a\right)}\\ \frac{\mathrm{d}\boldsymbol{v}}{\mathrm{d}a} & =-\nabla\Phi\left(\boldsymbol{x}\right) \end{aligned}\]其中：
\[\Phi\left(\boldsymbol{x}\right)=\frac{\nabla\phi\left(\boldsymbol{x}\right)}{aH\left(a\right)}\]因此，有：
\[\begin{aligned} \nabla^{2}\Phi\left(\boldsymbol{x}\right) & =\frac{4\pi Ga\bar{\rho}}{H\left(a\right)}\delta\left(\boldsymbol{x},t\right)\\ & =\frac{3}{2}\frac{H_{0}^{2}\Omega_{m0}}{a^{2}H\left(a\right)}\delta\left(\boldsymbol{x},t\right) \end{aligned}\]Zel’dovich 近似
Zel’dovich 近似（简称为 ZA）描述了粒子最初的动力学演化过程. 首先我们定义一个位移场 \(\boldsymbol{S}\left(\boldsymbol{q},t\right)\)，假设一个位于拉格朗日坐标 \(\boldsymbol{q}\) 的粒子，其在时间 \(t\) 的欧拉坐标 \(\boldsymbol{x}\) 可以写为：
\[\boldsymbol{x}\left(\boldsymbol{q},t\right)=\boldsymbol{q}+\boldsymbol{S}\left(\boldsymbol{q},t\right)\]ZA 的核心假设是，这个位移场 \(\boldsymbol{S}\left(\boldsymbol{q},t\right)\) 可以写成一个不随时间改变的空间部分 \(\boldsymbol{s}\left(\boldsymbol{q}\right)\) 和一个只随时间改变的时间部分 \(D\left(t\right)\)（即线性增长因子）的乘积：
\[\boldsymbol{S}\left(\boldsymbol{q},t\right)=D\left(t\right)\cdot\boldsymbol{s}\left(\boldsymbol{q}\right)\]在宇宙的早期，我们可以使用线性近似，根据质量守恒，微扰场 \(\delta\) 与位移场的关系是：
\[\delta=-\nabla_{q}\cdot\boldsymbol{s}\left(\boldsymbol{q}\right)\]进一步假设初始流动是无旋的，则位移场 \(\boldsymbol{s}\) 可以写成一个标量势 \(\psi\) 的梯度:
\[\boldsymbol{s}\left(\boldsymbol{q}\right)=-\nabla_{q}\psi\left(\boldsymbol{q}\right)\]从而得到关于位移势的泊松方程：
\[\delta=-\nabla_{q}^{2}\psi\left(\boldsymbol{q}\right)\]求解上述泊松方程，从而得到 \(\boldsymbol{s}\) 并进一步得到粒子的坐标变化（这里我们将初始时刻的 \(D\left(t\right)=1\)，这总可以通过调整 \(\delta\) 得到，并近似认为初始的拉格朗日坐标就是我们均匀放置的粒子网格的欧拉坐标 \(\boldsymbol{x}\)）
\[\boldsymbol{x}\to\boldsymbol{x}-\nabla\nabla^{-2}\delta\left(\boldsymbol{x}\right)\]粒子的本动速度 \(\boldsymbol{v}_{p}\) 是对总位移 \(\boldsymbol{S}\) 对时间求导：
\[\boldsymbol{v}_{p}=\frac{\mathrm{d}\boldsymbol{S}}{\mathrm{d}t}=\dot{D}\left(t\right)\boldsymbol{s}\left(\boldsymbol{q}\right)=fH\boldsymbol{S}\left(\boldsymbol{q},t\right)\]其中定义增长律 \(f=\mathrm{d}\left(\ln D\right)/\mathrm{d}\left(\ln a\right)\)，因此 \(\dot{D}=fHD\). 本动速度转化成共动速度，得到：
\[\boldsymbol{v}=a\boldsymbol{v}_{p}=afH\boldsymbol{S}\left(\boldsymbol{q},t\right)\]因此：
\[\boldsymbol{v}\to-aHf\nabla\nabla^{-2}\delta\left(\boldsymbol{x}\right)\]初始条件生成
从理论力学中我们知道，在 \(M\) 维空间中，描述一个粒子的完整状态需要 \(2M\) 个坐标（\(M\) 个广义坐标和 \(M\) 个广义动量），对于 \(N\) 个粒子，则一共需要确定 \(2MN\) 个坐标. 因此，对于一个 NBody 宇宙学模拟，我们需要确定 \(2MN\) 的初始值. 一般而言，我们通过天文观测得到的 CMB 功率谱生成出早期的引力势，而位置坐标和动量坐标则通过 ZA 给出.
要通过功率谱生成初始的引力势，首先我们先生成一个高斯随机场，作为初始的引力势. 然后我们对其进行傅里叶变换，在傅里叶空间将场乘上 \(\sqrt{P\left(k\right)}\) 缩放模长，从而得到初始时刻的微扰场（的傅里叶表示）：\(\hat{\delta}\left(\boldsymbol{k}\right)\) . 为了保证微扰场 \(\delta\left(\boldsymbol{x}\right)\) 是实数，我们要求：\(\hat{\delta}\left(\boldsymbol{k}\right)=\hat{\delta}^{*}\left(\boldsymbol{k}\right)\).
在有了初始的微扰场后，我们通过 ZA 给出每个粒子的位置坐标和动量坐标：
\[\begin{aligned} \boldsymbol{x} & \to\boldsymbol{x}-\nabla\nabla^{-2}\delta\left(\boldsymbol{x}\right)\\ \boldsymbol{v} & \to-aHf\nabla\nabla^{-2}\delta\left(\boldsymbol{x}\right) \end{aligned}\]上述即给出每个粒子的初始位置和速度.
另外这里粒子的初始位置之所以不是通过直接对初始的微扰场 \(\delta\left(\boldsymbol{x}\right)\) 进行采样是因为后者在数学上是一个泊松采样过程，其会引入散粒噪声，这种噪声的功率谱是白噪声，即在所有尺度上都会引入这一噪声，使得我们生成的微扰场的功率谱和我们想要的不同.
模拟过程
蛙跳法
在前文中，我们给出了粒子的位置和速度随标度因子 \(a\) 的关系：
\[\begin{aligned} \frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}a} & =\frac{\boldsymbol{v}}{a^{3}H\left(a\right)}\to\frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}a}=f_{x}\left(\boldsymbol{v}\right)\\ \frac{\mathrm{d}\boldsymbol{v}}{\mathrm{d}a} & =-\nabla\Phi\left(\boldsymbol{x}\right)\to\frac{\mathrm{d}\boldsymbol{v}}{\mathrm{d}a}=f_{v}\left(\boldsymbol{x}\right) \end{aligned}\]我们的目标就是数值求解上述方程组. 在宇宙学模拟中我们使用蛙跳法: 蛙跳法的精髓在于其在一个交错的网格上进行计算：
位置 \(\boldsymbol{x}_{i}\) 在整数时间步进行计算：\(a_{0},a_{1},\dots\)
速度 \(\boldsymbol{v}_{i+1/2}\) 在半整数时间步上计算：\(a_{1/2},a_{3/2},\dots\)
我们想从当前位置 \(\boldsymbol{x}_{i}\)和中间点的速度 \(\boldsymbol{v}_{i+1/2}\) 得到下一个位置 \(\boldsymbol{x}_{i+1}\)，我们将位置 \(\boldsymbol{x}_{i+1}\)在中间点 \(a_{i+1/2}\) 进行展开：
\[\boldsymbol{x}_{i+1}=\boldsymbol{x}_{i+1/2}+\boldsymbol{x}_{t+1/2}^{'}\frac{\Delta a}{2}+\dots\]同理，将 \(\boldsymbol{x}_{i}\) 在中间点 \(a_{i+1/2}\) 进行展开：
\[\boldsymbol{x}_{i}=\boldsymbol{x}_{i+1/2}-\boldsymbol{x}_{t+1/2}^{'}\frac{\Delta a}{2}+\dots\]一式减二式，得到：
\[\boldsymbol{x}_{i+1}-\boldsymbol{x}_{i}=\boldsymbol{x}_{t+1/2}^{'}\Delta a+O\left(\left(\Delta a\right)^{3}\right)\]因此：
\[\boldsymbol{x}_{i+1}\to\boldsymbol{x}_{i}+f_{x}\left(\boldsymbol{v}_{i+1/2}\right)\Delta a\]同理，对于速度有：
\[\boldsymbol{v}_{i+1/2}\to\boldsymbol{v}_{i-1/2}+f_{v}\left(\boldsymbol{x}_{i}\right)\Delta a\]上述过程即构成了蛙跳法的循环：
Kick：\(\boldsymbol{v}_{i+1/2}\to\boldsymbol{v}_{i-1/2}+f_{v}\left(\boldsymbol{x}_{i}\right)\frac{\Delta a}{2}\)
Drift：\(\boldsymbol{x}_{i+1}\to\boldsymbol{x}_{i}+f_{x}\left(\boldsymbol{v}_{i+1/2}\right)\Delta a\)
Kick：\(\boldsymbol{v}_{i+1}\to\boldsymbol{v}_{i+1/2}+f_{v}\left(\boldsymbol{x}_{i+1}\right)\frac{\Delta a}{2}\)
这里之所以我们使用蛙跳法而不是更常用的精度更高的 RK4 等求解方法，一个很重要的原因是蛙跳法是一种辛算法，其在数值上保持了哈密顿系统的几何结构. 不难发现，蛙跳法的每一步都是相空间的一个剪切变换，其不改变相空间的体积，因此其保持了能量和动量等守恒量在长时间的守恒. 而更常见的 RK4 则没有这一点.
PM 法
在 NBody 宇宙学模拟中，我们需要求解数以百万计甚至百亿计的暗物质粒子在自身引力作用下的演化. 我们可以直接计算任意两个粒子之间的引力，但这需要 \(\mathcal{O}\left(N^{2}\right)\) 的计算量，这对于大规模模拟是不可接受的. 粒子-网格法（PM 法）则可以将长程引力的计算复杂度降低到 \(\mathcal{O}\left(N\log N\right)\)，极大地提升了模拟效率.
PM 方法的核心思想是：将粒子所代表的离散质量分布”绘制”到一个规则的计算网格上，然后在网格上利用快速傅里叶变换高效地求解引力势，再从引力势计算出引力场，最后将网格上的引力场插值回粒子所在的位置，从而求解粒子的受力.
在前文中，我们得到引力势满足的泊松方程：
\[\nabla^{2}\Phi=\frac{3}{2}\frac{H_{0}^{2}\Omega_{m0}}{a^{2}H\left(a\right)}\delta\left(\boldsymbol{x},t\right)\]同时，我们的速度满足：
\[\frac{\mathrm{d}\boldsymbol{v}}{\mathrm{d}a}=-\nabla\Phi\left(\boldsymbol{x}\right)\]我们的目标就是求解出 \(\Phi\) 并计算其梯度 \(\nabla\Phi\)，从而给出 \(\mathrm{d}\boldsymbol{v}/\mathrm{d}a\)，PM 法将这个求解过程分为以下四个主要步骤：
质量分配：将粒子的质量分配到网格点上，得到网格化的密度场.
泊松方程求解：利用 FFT 在傅里叶空间中求解泊松方程，得到网格点上的引力势.
引力计算：通过对引力势求梯度，计算出网格点上的引力场.
力插值：将网格点上的引力场插值到每个粒子的实际位置.
质量分配
PM 法需要通过傅里叶变换来加速，而在计算机计算中，我们只能对一个网格进行傅里叶变换，因此我们需要将位置随意的粒子分配到网格上. 我们假设一共有 \(N\) 个粒子，模拟的盒子是一个边长为 \(L_{tot}\) 的立方体，被划分为 \(M\times M\times M\) 个单元，每个单元的边长为：\(L=L_{tot}/M\).
令 \(\boldsymbol{x}_{p}\) 为第 p 个粒子的位置，其质量为 \(m_{p}\). 网格点的值为 \(\rho_{i,j,k}\) （\(i,j,k\) 为网格的索引）. 可以通过一个分配函数 \(W\left(\boldsymbol{x}\right)\) 来计算. 这个函数将粒子位置 \(\boldsymbol{x}_{p}\) 与网格点索引 \((i,j,k)\) 关联起来.
\[\rho_{i,j,k}=\frac{1}{L^{3}}\sum_{p}m_{p}W(\boldsymbol{x}_{i,j,k}-\boldsymbol{x}_{p})\]其中 \(\boldsymbol{x}_{i,j,k}\) 是网格点 \((i,j,k)\) 的坐标.
最简单的一种质量分配方法就是将粒子的质量分配到离他最近的网格处（最近邻网格法，NGP）. 但这会造成质量的”瞬移”，带来数值噪音. 一种更常用的方法是云中单元法（CIC 法）. 其将每个粒子看成一个同网格单元同样大小的均匀立方体”云”，然后根据该”云”与周围网格单元的重叠体积，将质量按比例分配给相邻的网格点.
对于一维情况，粒子在 \(x_{p}\) 处，其左右两个网格点为 \(x_{i}\) 和 \(x_{i+1}\)，则分配的权重为：
\[w_{i}=1-\frac{x_{p}-x_{i}}{L}\quad w_{i+1}=\frac{x_{p}-x_{i}}{L}\]在三维中，就是这三个方向权重的乘积.
泊松方程求解
获得网格上的密度场 \(\rho\) 后，将其转化为微扰场 \(\delta\)，下一步是求解方程
\[\nabla^{2}\Phi=\frac{3}{2}\frac{H_{0}^{2}\Omega_{m0}}{a^{2}H\left(a\right)}\delta\left(\boldsymbol{x},t\right)\]我们将上述方程进行傅里叶变换，在傅里叶空间，微分算子 \(\nabla\) 变成了代数乘子 \(\mathrm{i}\boldsymbol{k}\)。因此，\(\nabla^{2}\) 对应于 \(-\left|\boldsymbol{k}\right|^{2}\)，从而上述方程变为：
\[-\left|\boldsymbol{k}\right|^{2}\hat{\Phi}\left(\boldsymbol{k}\right)=\frac{3}{2}\frac{H_{0}^{2}\Omega_{m0}}{a^{2}H\left(a\right)}\hat{\delta}\left(\boldsymbol{k}\right)\]于是我们可以直接得到：
\[\hat{\Phi}\left(\boldsymbol{k}\right)=-\frac{3}{2}\frac{H_{0}^{2}\Omega_{m0}}{a^{2}H\left(a\right)}\frac{1}{\left|\boldsymbol{k}\right|^{2}}\hat{\delta}\left(\boldsymbol{k}\right)\]这是一个非常简单的代数方程. 注意，对于 \(\boldsymbol{k}=0\) 的模式（对应宇宙的平均密度），分母为零. 这个模式不产生引力，通常直接设 \(\hat{\Phi}(\boldsymbol{k}=0)=0\). 得到 \(\hat{\Phi}\left(\boldsymbol{k}\right)\) 后，我们再通过逆变换得到 \(\Phi\left(\boldsymbol{x}\right)\).
在实际计算中，我们使用快速傅里叶变换（FFT）算法来进行傅里叶变换，其计算复杂度为 \(\mathcal{O}\left(N\log N\right)\).
引力计算
得到了网格上的引力势 \(\Phi\left(\boldsymbol{x}\right)\) 后，\(\mathrm{d}\boldsymbol{v}/\mathrm{d}a\) 可以通过求其负梯度得到：
\[\frac{\mathrm{d}\boldsymbol{v}}{\mathrm{d}a}=-\nabla\Phi\left(\boldsymbol{x}\right)\]这个梯度运算同样可以在傅里叶空间进行. 因为 \(\nabla\leftrightarrow\mathrm{i}\boldsymbol{k}\)，所以我们可以在得到 \(\hat{\Phi}\left(\boldsymbol{k}\right)\) 后，直接乘以 \(-\mathrm{i}\boldsymbol{k}\)，然后进行一次逆傅里叶变换得到 \(\mathrm{d}\boldsymbol{v}/\mathrm{d}a\).
当然，我们也可以通过有限差分得到，例如，\(x\) 方向的引力（\(x\) 方向的速度变化率）可以表示为：
\[\left(\frac{\mathrm{d}v_{x}}{\mathrm{d}a}\right)\left(i,j,k\right)=-\frac{\Phi(i+1,j,k)-\Phi(i-1,j,k)}{2L}\]实际中可以使用更高精度的差分公式.
力插值
最后一步是将网格点上的引力插值到每个粒子所在的实际位置 \(\boldsymbol{x}_{p}\)，得到作用在每个粒子上的引力 \(\mathrm{d}\boldsymbol{v}_{p}/\mathrm{d}a\). 为了保证整个模拟系统的动量守恒，力插值必须使用与质量分配完全相同的插值函数来将网格的力”分配”回粒子. 前文中我们使用了 CIC 法进行质量分配，因此这里我们也使用 CIC 方案进行力插值.
设一个粒子 \(p\) 位于 \(\boldsymbol{x}_{p}\)，它被一个由8个网格点构成的立方单元所包围. 我们在之前计算出了这 \(8\) 个角点上的引力向量 \(\left(\mathrm{d}\boldsymbol{v}/\mathrm{d}a\right)_{i,j,k}\). 作用在粒子 \(p\) 上的引力 \(\left(\mathrm{d}\boldsymbol{v}/\mathrm{d}a\right)_{p}\) 就是这 \(8\) 个网格点上引力的加权平均：
\[\left(\frac{\mathrm{d}\boldsymbol{v}}{\mathrm{d}a}\right)_{p}=\sum_{l,m,n\in\{0,1\}}\left(\frac{\mathrm{d}\boldsymbol{v}}{\mathrm{d}a}\right)_{i+l,j+m,k+n}\cdot w_{i+l,j+m,k+n}\]这里的权重 \(w\) 与质量分配的权重完全相同. 例如，对于一维情况，粒子在 \(x_{p}\) 处，其左右两个网格点为 \(x_{i}\) 和 \(x_{i+1}\)，则作用在粒子上的力为：
\[\left(\frac{\mathrm{d}v}{\mathrm{d}a}\right)_{p}=\left(\frac{\mathrm{d}v}{\mathrm{d}a}\right)_{i}\cdot\left(1-\frac{x_{p}-x_{i}}{L}\right)+\left(\frac{\mathrm{d}v}{\mathrm{d}a}\right)_{i+1}\cdot\left(\frac{x_{p}-x_{i}}{L}\right)\]在三维中，这同样是三个方向权重的乘积.
得到 \(\mathrm{d}\boldsymbol{v}/\mathrm{d}a\) 后，即可通过蛙跳法进行粒子位置和动量的更新，并以此循环.