网格怎样跟着折叠的物质求引力

粒子在节点之间，也能分清质量

若只把质量塞给最近节点，粒子刚越过半格，节点里的质量就会突然跳变。CIC把一颗粒子的质量按距离分给相邻节点：在一维是两处，在三维是八个角点。粒子仍留在原位置，分配的是求场用的权重。

把相邻两节点之间的分数位置记为 $f$ ，它从零到一。左、右权重分别是 $1-f$ 和 $f$ 。例如质量为一个单位、位置在四分之一路程处，两个节点分别获得 $0.75$ 与 $0.25$ 个质量单位。

w_0=1-f,\qquad w_1=f,\qquad w_0+w_1=1

权重非负且和为一。三维权重由三个方向的这类权重相乘，因此八个角点的和仍是一。

粒子处在四分之一和五分之四格位置时的CIC质量权重 — 自绘分配示意。白点是一颗粒子，彩色圆为分到节点的质量；箭头表示求场权重，不表示粒子运动。分数位置分别为 $0.25$ 与 $0.8$ ，每组左右权重合计一。查看完整图像

对应阅读：Springel (2005), The cosmological simulation code GADGET-2, astro-ph/0505010v1 · Original finite-mesh calculations and independent small-wave Python replay

先保住质量，再扣掉平均背景

设每轴有 $N_g$ 个节点，盒内有 $N_p$ 个等质量粒子。节点 $g$ 分到的粒子权重记为 $w_{gp}$ 。用一个网格体积来归一化，得到节点密度与同一时刻均密度之比。

把所有节点的质量加起来，应当恰好回到原来的总质量。先检查这一步，再定义密度对比度 $\delta=\rho/\bar\rho-1$ 。若先把结果强行重新归一化，就会藏起周期面上漏分配的质量。

\frac{\rho_g}{\bar\rho}=\frac{N_g^3}{N_p}\sum_p w_{gp},\qquad\sum_gw_{gp}=1

物理网格体积和平均密度中都有相同的膨胀因子，取比以后抵消。盒内总质量不随尺度因子改变。

对应阅读：Original finite-mesh calculations and independent small-wave Python replay · Springel (2005), The cosmological simulation code GADGET-2, astro-ph/0505010v1

每一条波，都有对应的势

周期密度可以拆成整数波矢 $\boldsymbol n$ 的波。其无量纲波数是 $\boldsymbol k_y=2\pi\boldsymbol n$ ；有长度单位的共动波数则是 $\boldsymbol k=2\pi\boldsymbol n/L$ 。整盒恒定的分量是零模式，均匀背景已经决定整体膨胀。

坐标篇使用的缩放势 $U$ 满足 $\nabla_y^2U=3\delta/2$ 。对单个Fourier波，空间二阶导数给出负的波数平方，因此可以逐波解出势。实际沉积会平滑波形，CIC回收力又带来一次平滑；这里对这两次平均窗口作显式修正。

定义 $\operatorname{sinc}x=\sin x/x$ ，零点取连续极限一。窗口 $W_{\boldsymbol n}$ 和保留范围写在下式。FFT只是快速完成离散变换的方法；保留哪些波仍由网格和截断决定。

W_{\boldsymbol n}=\prod_{j=1}^3\operatorname{sinc}^2\!\left(\frac{\pi n_j}{N_g}\right),\qquad U_{\boldsymbol n}=-\frac{3\delta_{\boldsymbol n}}{2|\boldsymbol k_y|^2W_{\boldsymbol n}^2}

此式用于非零且

|\boldsymbol n|\leq N_g/3

的模式；其他模式及零模式的势设为零。这个低通网格力仍有离散误差。

对应阅读：Springel (2005), The cosmological simulation code GADGET-2, astro-ph/0505010v1 · Original finite-mesh calculations and independent small-wave Python replay

把势的坡度取回粒子位置

变换回网格以后，力沿势下降的方向。这里对每条轴使用四点中心差分，网格间距 $h=1/N_g$ 按无量纲坐标量；再用与沉积相同的CIC权重，把节点力取回当前粒子位置。

以 $N_g=16$ 、 $L=20\,\mathrm{Mpc}$ 为例，节点间距是 $1.25\,\mathrm{Mpc}$ 共动长度，Fourier截断还会进一步限制短波。粒子可以挤到很小的间隔，求力却仍是这个较粗的模型。密集小点的内部细节因此需要另外的分辨率研究。

(\partial_jU)_i\simeq\frac{8(U_{i+1}-U_{i-1})-(U_{i+2}-U_{i-2})}{12h}

求势采用Fourier波数，求梯度采用这一离散差分；两者共同定义实际算子。这个受限PM没有叠加GADGET-2的短程树力。

对应阅读：Springel (2005), The cosmological simulation code GADGET-2, astro-ph/0505010v1 · Original finite-mesh calculations and independent small-wave Python replay

先推一下，再移动，再推一下

受力先改变动量，动量再推进位置。对称的kick–drift–kick方法在一步中先用旧位置的力更新一部分动量，再移动位置、重求引力，最后完成另一部分动量更新。

从 $a_0$ 走到 $a_1$ ，中点取 $a_m=\sqrt{a_0a_1}$ ，与对数尺度因子的等间隔一致。由于背景是EdS，位置和动量方程中的尺度因子可以直接积分。把两个因子记作 $\mathcal D$ 与 $\mathcal K$ ，以免和密度增长因子混淆。

\mathcal D(a_0,a_1)=2(a_0^{-1/2}-a_1^{-1/2}),\qquad\mathcal K(a_0,a_1)=2(\sqrt{a_1}-\sqrt{a_0})

它们分别来自对尺度因子积分位置方程与动量方程。两者都无量纲。

\boldsymbol u_h=\boldsymbol u_0+\mathcal K(a_0,a_m)\boldsymbol F(\boldsymbol y_0),\qquad\boldsymbol y_1=\boldsymbol y_0+\mathcal D(a_0,a_1)\boldsymbol u_h

新位置给出新力，再取

\boldsymbol u_1=\boldsymbol u_h+\mathcal K(a_m,a_1)\boldsymbol F(\boldsymbol y_1)

完成一步。

对应阅读：Springel (2005), The cosmological simulation code GADGET-2, astro-ph/0505010v1 · Original finite-mesh calculations and independent small-wave Python replay

步长细化，和网格细化分别检查

先做一个因子算例。区间 $a_0=0.25$ 到 $a_1=0.36$ 给出 $\mathcal D=2(2-1/0.6)=2/3$ ， $\mathcal K=2(0.6-0.5)=0.2$ 。实际轨迹把整个区间分成许多小步，并在每一步重算密度和力。

保持网格不变，把步数加倍，检查同一算子的时间误差；保持连续初值不变，再细化粒子与网格，检查空间与采样误差。独立自适应积分也能检验前者。前面平面例子的首穿越和后面的线性增长，则给出解析或独立力参考。

质量应当守恒。膨胀背景下，动能加引力势能一般不会保持常数；一次对称推进本身也不是精确能量守恒的证明。要信任某个读数，仍需看它在相应细化和独立参考下怎样变化。

对应阅读：Original finite-mesh calculations and independent small-wave Python replay · Springel (2005), The cosmological simulation code GADGET-2, astro-ph/0505010v1

网格怎样跟着折叠的物质求引力

粒子在节点之间，也能分清质量

先保住质量，再扣掉平均背景

每一条波，都有对应的势

把势的坡度取回粒子位置

先推一下，再移动，再推一下

步长细化，和网格细化分别检查

练习与解释性答案

动手检验预测

把位移参考和更新引力放在同一时刻比较

来源与计算范围