重算温暗物质截断、模糊暗物质截断和孤子核心

先读温暗物质和模糊暗物质怎样抹掉小起伏和造出暗物质核心的三种方式。你需要 Python、NumPy 和 SciPy，也要能读懂常微分方程。除非某一步另有说明，下面都是红移为零时的线性理论。下一章怎样读银河系卫星星系给温暗物质和模糊暗物质的限制会把半模质量用到一份已发表的卫星星系普查上。

材料包里有什么

所有文件都放在本站的同一个文件夹里。冻结参考网格保存了冷暗物质功率谱，以及用玻尔兹曼程序 CLASS 3.4.2 算出的温暗物质转移函数。生成网格的脚本是唯一需要 CLASS 的文件。两个基准脚本要么读这个网格，要么不需要输入：线性截断基准和孤立晕基准。它们保存的输出是 linear-output.json 和 isolated-output.json。第三个脚本 compute_teaching_values.py 重算第2到第4章例题里的每个数，写出 teaching-values.json。

宇宙学参数取普朗克2018年的 TT,TE,EE+lowE+lensing+BAO 组合，即Planck Collaboration 2020, arXiv:1807.06209v4表2最后一列：$h=0.6766$，$\omega_{\mathrm b}=0.02242$，$\omega_{\mathrm c}=0.11933$，$A_{\mathrm s}=2.105\times10^{-9}$，$n_{\mathrm s}=0.9665$，$\tau=0.0561$。中微子按无质量处理，这是明说的简化。网格有141个波数，从 $0.05$ 到 $300\,h\,\mathrm{Mpc^{-1}}$ 按对数均匀分布。$k$ 的单位是 $h\,\mathrm{Mpc^{-1}}$，$P$ 的单位是 $(h^{-1}\,\mathrm{Mpc})^3$。温暗物质的三次计算把全部暗物质交给质量为2、3或6.5 keV 的热遗迹粒子，再调节它的温度，让它的密度等于冷暗物质的值。

保存的输出用的是 Python 3.9.6、NumPy 2.0.2 和 SciPy 1.13.1，我们只在这组版本上运行过。先确认你手里的网格和这里的一样：它的 SHA-256 摘要以 04e1c5d5 开头，以 3bd48e 结尾。

shasum -a 256 reference-grid.json
python3 -I benchmark_linear_cutoffs.py reference-grid.json linear-output.json
python3 -I benchmark_isolated_soliton.py isolated-output.json
python3 -I compute_teaching_values.py reference-grid.json linear-output.json teaching-values.json

python3 后面那个参数让脚本在隔离模式下运行，当前文件夹里混进来的同名文件就不会顶替 NumPy 被导入。重新生成网格本身需要 CLASS 的 Python 接口 classy。

python3 -I generate_frozen_reference_grid.py reference-grid.json

检验一：温暗物质半模和 CLASS 比

Viel 的拟合给出振幅比 $T(k)=[1+(\alpha k)^{2\nu}]^{-5/\nu}$，其中 $\nu=1.12$，$\alpha$ 由粒子质量决定。Viel et al. 2005, equations 6 and 7 半模有解析解 $k_{\mathrm{hm}}=(2^{\nu/5}-1)^{1/(2\nu)}/\alpha$。CLASS 那边，我们取网格里存的比值 $T=\sqrt{P_{\mathrm{WDM}}/P_{\mathrm{CDM}}}$，找它第一次降到 $1/2$ 以下的位置，在相邻格点之间对 $\log k$ 和 $\log T$ 做线性插值。半模质量接着由下式得到：

Nadler et al. 2021, equation 2 其中 $\Omega_{\mathrm m}=0.3096$，$\rho_{\mathrm{crit}}=2.775\times10^{11}\,h^2\,\mathrm M_\odot\,\mathrm{Mpc^{-3}}$。

波数单位是 $h\,\mathrm{Mpc^{-1}}$，质量单位是太阳质量。质量随 $k_{\mathrm{hm}}^{-3}$ 变化，所以在6.5 keV，波数差12%，质量就差45%。粒子越重，差距越大，而且 CLASS 的截断总比拟合落在更小的尺度上。Viel 等人选 $\nu=1.12$，是因为它对 $k<5\,h\,\mathrm{Mpc^{-1}}$ 的玻尔兹曼程序结果拟合得最好。这里的半模在20到86之间，在那个范围之外，本来就不该指望百分之一级别的吻合。

第二个比较对象是别人发表的数。Nadler 等人写的是 $M_{\mathrm{hm}}=5\times10^8\,\mathrm M_\odot\,(m_{\mathrm{WDM}}/3\,\mathrm{keV})^{-10/3}$。Nadler et al. 2021, equation 5 他们的模拟用 $h=0.7$、$\Omega_{\mathrm m}=0.286$。把这组参数和温暗物质密度 $0.24$ 代进 Viel 拟合，3 keV 时得到 $4.81\times10^8\,\mathrm M_\odot$，比他们的值低4%。温暗物质密度取 $0.24$，是我们替他们的设置做的假设。

检验二：模糊暗物质截断和直接解增长方程比

Hu、Barkana 和 Gruzinov 把模糊暗物质的振幅比拟合成 $T=\cos x^3/(1+x^8)$，其中 $x=1.61\,m_{22}^{1/18}\,k/k_{\mathrm{J,eq}}$，$k_{\mathrm{J,eq}}=9\,m_{22}^{1/2}\,\mathrm{Mpc^{-1}}$，并说功率在 $k_{1/2}\approx4.5\,m_{22}^{4/9}\,\mathrm{Mpc^{-1}}$ 附近降到一半。Hu, Barkana and Gruzinov 2000, equations 8 and 9 这里 $m_{22}$ 是以 $10^{-22}\,\mathrm{eV}$ 为单位的粒子质量，这些波数的单位是 $\mathrm{Mpc^{-1}}$，不带 $h$。

独立的检验是直接积分单个波的线性增长。把波动方程写成流体形式，粒子的波长就表现为一种压强。用 $\delta$ 表示密度反差，以 $\ln a$ 为时间：

$\Omega_{\mathrm{cl}}(a)$ 是当时参与成团的物质占临界密度的比例。脚本在复合，即 $a=1/1100$ 之前只让暗物质成团，之后让全部物质成团，背景里有辐射、物质和宇宙学常数。每个波从进入视界时开始，取 $\delta=1$、$d\delta/d\ln a=1$，再把量子项去掉，把同一个波重算一遍。两次终值之比就是振幅比，脚本找出它的平方穿过 $1/2$ 的位置。

增长方程截断功率的波数比拟合小11%到17%。这个计算是简化模型：没有重子压强，没有和光子的耦合，在复合时刻是突然切换的。这张表只说明拟合和这个模型差这么多，说明不了谁更接近完整的模糊暗物质玻尔兹曼计算。下一节会看到，结果还取决于起始斜率，这也是原因之一。

检验三：把孤子当作精确基态来解

模糊暗物质晕的核心，是这种波在自身引力下能量最低的状态。在 $\hbar/m=1$、$G=1$ 的单位里，它满足

其中 $\rho=\psi^2$，$u$ 是引力势减去能量。脚本固定 $\psi(0)=1$，用二分法调 $u(0)$，直到 $\psi$ 降到零而不穿过零。方程在 $r\to r/\lambda$、$\rho\to\lambda^4\rho$ 下形式不变，所以对每个基态，$\rho_0r_{\mathrm c}^4$ 都是同一个数，$r_{\mathrm c}$ 是密度降到一半的半径。

Schive 等人用 $\rho_0=1.9\,(m/10^{-23}\,\mathrm{eV})^{-2}(r_{\mathrm c}/\mathrm{kpc})^{-4}\,\mathrm M_\odot\,\mathrm{pc^{-3}}$ 和形状 $[1+0.091(r/r_{\mathrm c})^2]^{-8}$ 去拟合他们模拟出来的核心。Schive et al. 2014, PRL, equation 3 换成这里的单位，他们的前因子是 $0.2223$。

两种容差 $10^{-7}$ 和 $10^{-10}$ 给出的数前六位相同。要让形状恰好在 $r_{\mathrm c}$ 处减半，系数应是 $2^{1/8}-1=0.0905$，$0.091$ 就是它四舍五入的结果，所以中心部分吻合得很好。偏差往外越来越大：在 $r_{\mathrm c}$ 处是0.4%，在 $2r_{\mathrm c}$ 处是4.0%，最大的5.9%出现在 $2.7r_{\mathrm c}$。那篇快报里说“精确到2%”的那句话，没有讲是和什么比。我们的5.9%是和精确基态比的，没有拿拟合去和他们的模拟核心比。

核心半径以内的质量乘以这个半径，也被同一个标度关系固定住。对 $m=10^{-23}\,\mathrm{eV}$，精确基态给出 $M_{\mathrm c}r_{\mathrm c}=5.41\times10^9\,\mathrm M_\odot\,\mathrm{kpc}$。把 Schive 公式从 $0$ 积到 $r_{\mathrm c}$，得到 $5.30\times10^9$，低2%，和前因子的差距差不多。

同一个脚本还检验自相互作用模型用到的等温球。数值解和级数 $x^2/6-x^4/120+x^6/1890$ 在 $x\le0.4$ 范围内相差不超过 $2.4\times10^{-8}$。在远处，密度乘以 $x^2/2$ 应当趋于1，在 $x=200$、$600$、$2000$ 处读数是 $1.016$、$1.049$、$0.997$。等温球本来就是这样一边摆动一边慢慢逼近的。密度在 $x=2.2746$ 处减半。

每项检验可能怎样出错

下面每种出错情形都存进了 linear-output.json，你可以确认自己能复现。

模糊暗物质拟合会变号。$\cos x^3$ 在 $x=(\pi/2)^{1/3}$ 以上变负，对 $m_{22}=1$ 就是 $k=6.50\,\mathrm{Mpc^{-1}}$。这个拟合只描述截断的形状，过了这一点只能用 $T^2$。在那里取 $\log T$ 的代码会报错停下。

定步长积分器会发散。在 $k=30\,\mathrm{Mpc^{-1}}$、$m_{22}=1$ 时，早期的量子项非常大，方程变得刚性。自适应求解器给出振幅比 $3.69\times10^{-4}$，40个固定步长的四阶龙格–库塔给出 $-1.86\times10^{26}$。结果大成这样或者变成负数，说明步长太长，而不是物理变了。

起始条件有影响。默认起点下，$m_{22}=1$ 的半功率点是 $3.80\,\mathrm{Mpc^{-1}}$。把初始斜率设为零，得到 $3.22$。把起点的尺度因子提前十倍，得到 $3.31$。光是起点就带来15%的差别，这又是一个不该把表5.2当成精确结果的理由。

温暗物质拟合在深尾部失效。在 CLASS 第一次降到 $T\approx0.017$ 的格点上，Viel 拟合已经是2 keV 的 $0.0027$、3 keV 的 $0.0019$、6.5 keV 的 $0.0011$，低了6到17倍。即使在 $T$ 大于 $0.25$ 的范围内，三种质量的拟合偏差也分别达到16%、23%和31%。尾部要紧时，直接读 CLASS 的结果。

自己动手时的检查单

1. 开始之前先确认网格摘要。
2. 运行三个脚本，把你的输出文件和保存的文件比较，目标是前四位有效数字一致。
3. 一次只改一个输入。比如在生成脚本里试4 keV 的粒子，或在增长基准里试 $m_{22}=2$，运行之前先预测变化方向。
4. 收紧容差再跑一次，报告哪几位数字没有动。
5. 写下每项比较没有检验到的东西。两个计算如果有同样的简化，彼此吻合也说明不了太多。

实验任务

1. 预测。把暗物质晕实验设成3 keV 的温暗物质，读半模质量。它是表5.1里的拟合值，还是 CLASS 值？
2. 改一个条件。把模糊暗物质粒子质量从 $m_{22}=1$ 调到 $m_{22}=4$，读孤子核心半径。晕质量不变时，它应缩小到四分之一，因为 $r_{\mathrm c}\propto m^{-1}$。
3. 解释。用表5.3说明，为什么实验用 Schive 公式画孤子，不会把核心的位置画得明显不对。

练习

1. 从解析解出发，证明在宇宙学参数固定时，温暗物质半模波数按 $m_{\mathrm{WDM}}^{1.11}$ 变化，半模质量按 $m_{\mathrm{WDM}}^{-3.33}$ 变化。
2. 用表5.2给增长方程的半功率点拟合一个幂律。它还保持 Hu 式9的指数 $4/9$ 吗？
3. 40步的求解器在 $k=30\,\mathrm{Mpc^{-1}}$ 失败。估计早期量子项相对引力项有多大，并解释为什么在 $\ln a$ 上取固定步长时，步子太长了。