对照萨哈、三能级原子与冻结的CLASS参考输出

参考对象是CLASS仓库的v3.3.4标签，提交号e85808324f51fc694d12e3ed7439552a3c3f9540。GitHub上的CLASS v3.3.4 我们从中读取四个文件。为这一页，我们没有编译或运行CLASS、HyRec或CosmoTherm。下文所有“参考值”，都是这个发行版自带文件里印好的数值。

文件、版本与单位

所有脚本只用Python 3标准库，运行环境是Python 3.9.6。材料放在本合集的材料文件夹里：

ac75_model.py：参考模型，实验室的TypeScript代码与它一一对应；
inputs.json：背景参数、求解器设置和五次能量释放；
run_recombination.py 及其保存的输出 recombination-output.json；
check_distortion_kernels.py 与 kernel-checks.json；
compare_class_files.py 与 class-comparison.json；
hyrec-reference-extract.json：HyRec-2文件的60行抽样，另有说明文件。

红移是无量纲的，脚本对 $s=\ln(1+z)$ 积分。密度以 $\mathrm{m^{-3}}$ 计，速率以 $\mathrm{s^{-1}}$ 计，能量在程序内部用焦耳、正文里用电子伏。频谱以 $\mathrm{kJy\,sr^{-1}}$ 计， $1\,\mathrm{kJy\,sr^{-1}}=10^{-23}\,\mathrm{W\,m^{-2}\,Hz^{-1}\,sr^{-1}}$ 。背景参数取 $h=0.6736$ 、 $T_0=2.7255\,\mathrm K$ 、 $\Omega_{\mathrm b}=0.04941$ 、 $\Omega_{\mathrm{cb}}=0.3124$ 、氦质量分数 $0.245$ 、 $N_{\mathrm{eff}}=3.046$ 。这些与CLASS里external/HyRec2020/input.dat一致，只是把那一种有质量的中微子当作无质量辐射处理。

mkdir ac75-check && cd ac75-check
base=https://sandboxphysics.com/stories/cosmic-transparency/materials
for name in ac75_model.py inputs.json run_recombination.py check_distortion_kernels.py compare_class_files.py hyrec-reference-extract.json; do curl -fsSLO "$base/$name"; done
python3 -I run_recombination.py inputs.json hyrec-reference-extract.json my-recombination-output.json
python3 -I check_distortion_kernels.py my-kernel-checks.json
git clone --depth 1 --branch v3.3.4 https://github.com/lesgourg/class_public.git
python3 -I compare_class_files.py class_public my-class-comparison.json

四个CLASS文件中只要有一个的SHA256与脚本里写的不同，对照脚本就拒绝运行。把你新生成的三个JSON文件和保存的版本比一比；换一台机器，最后几位数字可能不同。

第一步：三种复合历史在同一组红移上的值

脚本从红移1800处的萨哈平衡出发，向下求解三能级方程直到红移50，在 $s$ 上用 $5\times10^{-4}$ 的后向欧拉步长。后向欧拉用每一步终点处的速率，并用牛顿迭代解出它，这样能让早期的刚性阶段保持稳定。

**Table 5.1** recombination-output.json里的自由电子比例。
红移	萨哈平衡	三能级，因子1.14	三能级，因子1.00	HyRec-2文件	三能级与HyRec-2的相对差
$1500$	$0.947$	$0.9546$	$0.9546$	$0.9549$	$-0.0003$
$1300$	$0.210$	$0.5646$	$0.5646$	$0.5612$	$+0.0060$
$1100$	$0.00469$	$0.1424$	$0.1426$	$0.1447$	$-0.0160$
$900$	$1.6\times10^{-5}$	$0.01223$	$0.01249$	$0.01269$	$-0.0367$
$200$	$1.2\times10^{-53}$	$3.39\times10^{-4}$	$3.83\times10^{-4}$	$3.36\times10^{-4}$	$+0.0064$

看哪里：最后一列。红移1500及以下的采样点里，最大偏差是 $0.039$ ，出现在红移100。常见误读：因子1.14在红移1100几乎不起作用，却让红移200处的自由电子比例变了13%，因为它进入的是复合系数，而复合系数决定了自由电子最后冻结在什么水平。

接着看可见度概况。三能级历史的峰值在1081.0，在1092.0处达到光学深度1。对红移1600以下的HyRec-2文件积分，得到1080和1089.8，那里的氦已经无关紧要。普朗克2018列出 $z_\ast=1089.92\pm0.25$ 。Planck合作组，2018，表2 这里没有核对那篇论文里 $z_\ast$ 的定义。

第二步：数值检验与失败案例

步长减半后，光学深度达到1的红移移动 $0.078$ ，红移1100处的 $x$ 相对改变 $-5.2\times10^{-4}$ 。两者都远小于与HyRec-2之间的差距，所以这个差距来自物理，而不是步长。

有两个失败案例是特意保存的。用 $5\times10^{-3}$ 的显式欧拉步长，一开始就越过了萨哈平衡，到红移1500时得到nan。把三能级解的起点从红移1800改到1300，红移1100处得到 $x=0.105$ ，而正确值是 $0.142$ 。到红移1300时真实历史已经离开了萨哈平衡，起点本身就错了，这个误差一直带到红移1100。

第三步：完整的HyRec-2文件

对照脚本读取external/HyRec2020/output_xe.dat的全部8000行。在红移50到1600之间，最大相对差是 $+8.8\%$ ，出现在红移50。模型让气体温度一直等于光子温度，直到红移50；真实气体在黑暗时代后期冷得更快，会改变复合速率，所以这个假设是最可能的原因，但我们没有单独验证。在红移8000，文件给出 $x_{\mathrm e}=1.16341$ ，而氦完全电离时 $1+2f_{\mathrm{He}}=1.16343$ 。这说明文件计入了氦的电子，所以在红移1600以上比较时要考虑它们。

第四步：不重复使用公式来检验畸变形状

check_distortion_kernels.py拿闭合式 $G$ 、 $Y$ 、 $M$ 与它们不包含的量做比较：

用辛普森法则算能量和光子数，与精确的普朗克积分相比： $G$ 携带4倍黑体能量和3倍光子数； $Y$ 是4和 $4\times10^{-9}$ ； $M$ 是 $0.71395$ 和 $4\times10^{-5}$ 。 $0.71395$ 的倒数是 $1.4007$ ，正是 $\mu$ 用的系数。
把有限差分形式的康帕涅茨算符作用在普朗克占据数上，与 $Y$ 相差不到峰值的 $1.7\times10^{-7}$ 。
取 $\mu=10^{-6}$ 的玻色–爱因斯坦谱，重新拟合温度使光子数不变，与 $M$ 相差不到峰值的 $6.6\times10^{-5}$ 。

零点分别落在： $Y$ 在 $217.51\,\mathrm{GHz}$ ， $M$ 在 $124.5\,\mathrm{GHz}$ 。

第五步：格林函数里的能量

CLASS里的Greens_data.dat存放了120个红移上单位能量释放的频谱响应。CLASS的说明文件写明它改编自Jens Chluba公开的CosmoTherm 1.0.3；方法见Chluba，2013，第2节。把每个表列响应在文件的频率网格上积分，再除以黑体能量，就得到这一列承载了多少比例的释放能量。

**Table 5.2** 表列格林函数承载的能量。
释放红移	承载能量与释放能量之比
$1.0\times10^{3}$	$1.0000$
$1.1\times10^{5}$	$1.0004$
$2.6\times10^{5}$	$0.9991$
$5.6\times10^{5}$	$0.986$
$1.8\times10^{6}$	$0.931$
$3.3\times10^{6}$	$0.838$
$5.0\times10^{6}$	$1.013$

红移 $2.6\times10^{5}$ 以下，能量守恒到0.1%；再往上，偏差可达16%。按说明文件，这个文件另有几行，单独存放最终的黑体温度移动和它带走的能量。这里的积分没有把这几行算进去。早期的能量大多落进这个温度移动，所以它很可能是原因，但我们没有把它加回去验证。只用这一列做能量检验时，请限于 $2.6\times10^{5}$ 以下。

第六步：拆分畸变的两种方式

FIRAS_branching_ratios.dat对每个释放红移，列出单位 $\Delta\rho/\rho_\gamma$ 对应的 $\mu$ 、 $y$ 和温度移动。对照脚本在FIRAS从68到639吉赫的43个频率上，重复external/distortions/generate_PCA_files.py的步骤：按 $Y$ 、 $M$ 、 $G$ 的顺序做格拉姆–施密特正交化，再把每个格林函数投影上去。这一流程就是Lucca等人写下的式3.10到3.15。Lucca等，2020，第3.2.2节重复的结果与表格一致：在 $z=5\times10^{4}$ ，表格给出单位能量的 $\mu$ 为 $1.501$ ，我们在最近的表列红移 $4.96\times10^{4}$ 处投影得到 $1.488$ 。

实验室用的是平滑的分配拟合式，它们在过渡时期与表格差别很大： $z=5\times10^{4}$ 时是 $0.743$ ，表格是 $1.501$ ； $z=2\times10^{4}$ 时是 $0.177$ ，表格是 $0.341$ 。在 $5\times10^{5}$ 以上，两者的 $\mu$ 相差不到3%；在 $2\times10^{5}$ 相差不到9%；在 $5\times10^{3}$ 以下，两者的 $y$ 相差不到0.3%。

常见误读：把其中一个当成“正确的” $\mu$ 。投影结果取决于所用的频率和基矢的先后顺序。Lucca等人在脚注13里指出，换一种顺序， $y$ 、 $\mu$ 和温度三部分的相对大小就会变，而且 $y$ 与 $\mu$ 的区分本身是人为的。两种方法必须一致的，是总的频谱和它的能量。正因如此，实验室把 $\mu$ 和 $y$ 显示为一个明确说明的近似下的振幅，把FIRAS上限画成参考形状。

这些对照证明了什么，没证明什么

它们说明：实验室使用的氢复合历史和可见度，与一份详细的已发表计算相差在百分之几以内；步长已经收敛；畸变形状及其能量系数是对的。它们没有检验氦、红移200以后的气体温度、莱曼 $\alpha$ 辐射转移，也没有检验注入能量怎样逐级变成热。CLASS没有运行，所以这里没有任何内容检验CLASS自己怎样处理一次给定的能量释放。

下次做代码对照时可以复用的清单：

固定参考版本，记录文件校验和。
逐个对齐背景参数，并列出没能对齐的。
既在固定红移上比较，也比较可见度峰值这类导出量。
步长减半，确认与参考值的差距大于这一改变。
至少保存一个失败案例，让读者看到错误设置是什么样子。
每种形状都要和一个它不包含的量对照：一个矩、一个极限，或一个独立的算符。
比较拆分参数之前，先确认双方的拆分定义相同。

HyRec-2抽样和三个CLASS文件都注明了作者来源后使用。CLASS v3.3.4的目录树里没有找到许可证文件，所以Sandbox Physics只再分发那份60行的抽样，原始文件请读者从发行版下载。