用冻结反应网格复算小网络，并传播主要反应率误差

什么是冻结反应网格

一个核合成网络的大部分计算量，花在只依赖温度的量上：热核反应率、中子质子之间的弱反应率，以及光子电子等离子体的能量密度。只要把弱反应率写成以 $1/\tau_n$ 为单位，这些量就都与重子密度、额外辐射和中子寿命无关。程序包在 173 个温度点 $\log_{10}T_9=-2.3,-2.275,\dots,2.0$ 上把它们算一次，存进 frozen-rate-grid.json。网络读入网格，在 $\log T_9$ 上插值，任意改变 $\eta_{10}$、$\Delta N_{\mathrm{eff}}$ 或 $\tau_n$ 重算一次都不到一秒。Lab 用的是逐字节相同的同一份网格。

程序包是原创代码，采用 MIT 许可，不含第三方代码；Serpico 的系数和 AME2020 的质量是转录的数字，在使用它们的文件里注明了出处。

下载链接：MANIFEST.json、input.json、expected.json、frozen-rate-grid.json、teaching-values.json 和 provenance.json。Python 文件的地址规律相同，每个文件名一个地址。

单位和约定

读任何数字之前，先把这些弄清楚。

• 温度：$T_9=T/10^9\,\mathrm K$，$k_BT=0.08617\,T_9\,\mathrm{MeV}$。网络从 $T_9=100$（$8.6\,\mathrm{MeV}$）算到 $T_9=0.01$（$0.86\,\mathrm{keV}$），步长在 $\ln T$ 上均匀。
• 反应率：$N_A\langle\sigma v\rangle$，单位 $\mathrm{cm^3\,mol^{-1}\,s^{-1}}$，以 $\log_{10}$ 存储。每个反应还存有 $\delta f/f$，即相对一倍标准差误差。
• 弱反应率：$\tau_n\Gamma_{n\to p}$ 和 $\tau_n\Gamma_{p\to n}$，无量纲。
• 丰度：n、p、d、t、${}^3\mathrm{He}$、${}^4\mathrm{He}$、${}^7\mathrm{Li}$ 和 ${}^7\mathrm{Be}$ 的 $Y_i=n_i/n_b$。重子数守恒要求 $\sum_iA_iY_i=1$，每次保存的计算都记录了这个和。
• 输出沿用 PRIMAT 的约定：$Y_p=4Y({}^4\mathrm{He})$；$\mathrm{D/H}=Y_d/Y_p$；${}^3\mathrm{He/H}$ 把氦三和氚一起算；${}^7\mathrm{Li/H}$ 把锂七和铍七一起算，因为氚和铍七之后会衰变。
• 输入：$\eta_{10}$ 是今天的重子光子比乘以 $10^{10}$；$\Delta N_{\mathrm{eff}}$ 是在三种中微子之外额外加入的辐射，按中微子种数计；$\tau_n$ 以秒为单位，默认 $878.4\,\mathrm s$，取自粒子物理评述 2025 年版 24.2 节。

第一步：下载并校验文件

把文件放进一个新建的空文件夹。如果用镜像站，替换主机名即可。

mkdir bbn-frozen-rate-network && cd bbn-frozen-rate-network
base=https://sandboxphysics.com/stories/primordial-nucleosynthesis/materials
curl -fsSO "$base/MANIFEST.json"
for name in $(python3 -c "import json; print(' '.join(f['name'] for f in json.load(open('MANIFEST.json'))['files']))"); do
  curl -fsSO "$base/$name"
done
python3 -c "import hashlib, json; bad=[f['name'] for f in json.load(open('MANIFEST.json'))['files'] if hashlib.sha256(open(f['name'],'rb').read()).hexdigest()!=f['sha256']]; print('digests OK' if not bad else bad)"

最后一行应当打印 digests OK。只要有一个对外提供的文件与 MANIFEST.json 中的大小和摘要不符，本站的构建就会停止；出现不符，就说明文件离开本站之后被改动过。

第二步：复算保存的结果

python3 -B replay.py actual-result.json

运行约两分钟，输出

wrote actual-result.json
PASS

replay.py 会由 Python 模块重建网格并与保存的网格比较，重跑 input.json 中的每个情形，重新计算基准比较、敏感度、误差预算和蒙特卡罗，然后以 $10^{-7}$ 的相对容差逐个与 expected.json 比较。保存的 gridRebuildMaxDifference 恰好是 $0$：保存的网格正是这些模块生成的。

程序包只在 macOS 上的 CPython 3.9.6 运行过。换一个平台，math.exp 的最后几位可能不同，刚性求解器可能把差别放大到超过 $10^{-7}$。如果看到 FAIL，而差别出现在第七、八位有效数字，就改用下面基准比较里的物理容差判断；如果差在第三位，那就是真出了问题。

input.json 有五个情形：标准宇宙（$\eta_{10}=6.12$）、重子稀少（$2$）、重子丰富（$15$）、额外辐射（$\Delta N_{\mathrm{eff}}=1.5$）和辐射不足（$-1$）。每个情形都用 2000 步和 Lab 的 1000 步各算一次。标准情形在 expected.json 中是：

{
  "heliumMassFraction": 0.2423956078641409,
  "deuteriumToHydrogen": 2.5482022079050184e-05,
  "helium3ToHydrogen": 1.018552513743762e-05,
  "lithium7ToHydrogen": 4.5131264294261034e-10,
  "baryonNumberSum": 0.999999999999992
}

用 Lab 的 1000 步，D/H 是 $2.5570\times10^{-5}$；用 4000 步是 $2.5456\times10^{-5}$。所以 Lab 的时间步误差在 D/H 上约 0.4%，其他输出更小。

第三步：逐一阅读基准比较

expected.json 的 benchmarks 部分把网络和三份已发表的计算作比较。每项比较检验的东西不同，按下面的顺序读。

Serpico 等人，表 3 和表 4

这些结果用的是与网络相同的反应率拟合，所以检验的是转录和网络本身，而不是物理。Serpico 以 $\omega_b$ 给出结果，取 $\tau_n=885.7\,\mathrm s$，并有 $\eta_{10}=273.49\,\omega_b/(1-0.007\,Y_p)$；input.json 记下了这两点。在 $\omega_b=0.022$、$0.023$ 和 $0.024$ 处，网络复现的 D/H 相差不到 0.4%，氦三不到 0.8%，锂七偏低 1.1% 到 2.5%。常见的误读是把这种一致当作网络准确的证明。它只说明拟合抄对了、积分也对了；里面的物理仍是 2004 年的物理。

PRIMAT，表 V 的 Born 一行

PRIMAT 列出了弱反应率的每项修正怎样移动丰度，起点是玻恩近似，也就是教程实现的那一级。在 PRIMAT 的输入 $\eta_{10}=6.0913$、$\tau_n=879.5\,\mathrm s$ 下，网络给出 $Y_p=0.24258$，PRIMAT 是 $0.24262$，相差 $-0.018\%$。这检验的是弱反应率、等离子体热力学和时钟。这里只比较氦。PRIMAT 用的是自己的核反应率，所以它的 D/H $2.423\times10^{-5}$ 和锂 $5.635\times10^{-10}$ 不是基于 2004 年拟合的网络应当对上的目标；网络给出的是 $2.569\times10^{-5}$ 和 $4.47\times10^{-10}$。同一张表还显示，教程略去的修正会把 $Y_p$ 提高到 $0.24709$。

Fields 等人，表 II

Fields 在 $\eta_{10}=6.129$、$\tau_n=879.4\,\mathrm s$ 处列出对数敏感度。对于与该表共有的输入，网络复现了全部 60 个数，差别都不超过 $0.036$；最大的一项是锂对 $\eta$ 的敏感度，$2.058$ 对 $2.094$。敏感度检验的是网络响应的形状，而不是绝对值，误差预算用的正是它们。

另有两项检查比较每多一种中微子时氦的斜率。网络给出 $0.01341$；随 CLASS 程序发布的、基于 PRIMAT 的表 sBBN_2025_primat.dat 给出 $0.01345$，评述给出 $0.013$。

第四步：传播反应率误差

propagate_rate_uncertainties.py 把每个反应率误差当作对数正态分布：同一反应在不同温度下完全相关，不同反应之间相互独立。对一个标准正态数 $z$，反应率变为

一倍标准差预算每次只对一个反应取 $z=+1$ 和 $z=-1$，记录 $\ln X$ 变化的一半；中子寿命移动 $\pm0.5\,\mathrm s$。在 $\eta_{10}=6.12$ 处，平方和总误差是 D/H 0.69%、氦三 0.75%、锂七 1.76%。氦四只变 0.04%，几乎全部来自中子寿命。

蒙特卡罗用固定种子 20261112 抽取 200 组 $z$，所以能精确复现。它给出的相对标准差，D/H 0.72%、锂七 1.78%，与平方和结果一致。这种一致检验的是线性近似；只有 200 个样本，标准差本身带有约 5% 的统计误差。

要在别的重子密度下计算误差预算，可以直接调用模块：

from network_solver import FrozenGrid, load_grid
from propagate_rate_uncertainties import one_sigma_budget

grid = FrozenGrid(load_grid('frozen-rate-grid.json'))
budget = one_sigma_budget(grid, {'eta10': 5.0, 'neutron_lifetime': 878.4, 'steps': 2000})
print(budget['quadratureTotal'])

第五步：重新生成课程中的数值

python3 -B teaching_values.py my-teaching-values.json

运行约十五秒，应当复现 teaching-values.json：$T_9=7.71$ 处的冻结点，$T_9=10$ 到 $0.1$ 的时钟，247 秒时氦形成一半，平衡估算的氘温度，各组比较计算，以及由观测氘推出的重子密度 $\eta_{10}=6.180$。

值得亲手复现的失败情形

failureCases 部分保留了这个程序包必须避开的错误。每一条都可以拿去检验别的程序。

• 拟合超出适用范围。在整个网格上计算，每个 Serpico 拟合都会在某处变为非正值，依反应不同，在 173 个网格点中占 6 到 65 个。拟合只在它声明的温度范围内有效。rate_fit_windows.py 在窗口下方用库仑势垒形式延续每个拟合，在上方保持常数。把所有下边界移动 $0.75$ 倍或 $1.33$ 倍，输出变化都小于 $2\times10^{-4}$，结果不依赖这个选择。
• 一个印刷符号。若照字面使用，$d(p,\gamma){}^3\mathrm{He}$ 拟合的势垒因子会随温度降低而增大，在 $T_9=0.1$ 处给出 $1452$ 而不是 $5.59$（单位同上）。程序包采用衰减形式，并在 serpico_rate_fits.py 和 provenance.json 中记录了核查；${}^3\mathrm{He}(\alpha,\gamma){}^7\mathrm{Be}$ 误差拟合下支的另一个符号也作同样处理。
• 一阶时间步。在同样 2000 步下，把步长加倍外推的后向欧拉格式换成普通后向欧拉，与 4000 步参考相比，D/H 偏高 6.4%，锂七偏低 8.5%。只看氦是发现不了的。
• 步数太少。250 步时 D/H 是 $2.651\times10^{-5}$，偏高 4%。在报出第四位数字之前，总要用加倍的步数再算一次。

程序包没有包括什么

它的弱反应率只算到玻恩近似，中微子瞬间脱耦，所以与完整计算相比 $Y_p$ 低约 $0.005$。它全程使用 2004 年的核反应率，所以给定 $\eta$ 下的 D/H 不是当前最好的预言。它没有等离子体的量子电动力学修正，并把额外辐射当作从一开始就存在的常数 $\Delta N_{\mathrm{eff}}$。与观测比较时这些限制很重要。上面每一项检查都是同类相比，结论仍然成立。

遇到下一个网络时的检查清单

• 它用哪套反应率、哪些弱反应率修正、怎样处理中微子脱耦？
• 它是否在舍入误差范围内守恒重子数，并记录了这项检查？
• 它是否与自己的来源表一致，又是否在相同输入下与一个独立计算一致？
• 是否用加倍的时间步重算过？
• 反应率误差是否逐个反应、带温度依赖地传播，并用固定种子的蒙特卡罗检验过？
• 在给出任何新物理限制之前，是否先说明了模型的绝对偏差？