--- title: "用冻结反应网格复算小网络，并传播主要反应率误差" description: "下载冻结反应网格和Python脚本，重跑十二反应网络，对照已发表结果，再传播反应率误差。" canonical_url: "https://sandboxphysics.com/zh/stories/primordial-nucleosynthesis/recompute-a-frozen-rate-network/" language: "zh-CN" page_type: "story" alternate_languages: en: "https://sandboxphysics.com/stories/primordial-nucleosynthesis/recompute-a-frozen-rate-network/" last_modified: "2026-11-12" --- [宇宙最初几分钟与轻元素](https://sandboxphysics.com/zh/stories/primordial-nucleosynthesis/) # 用冻结反应网格复算小网络，并传播主要反应率误差 Sandbox Physics · 2026-11-12 读完这篇方法指南，你会在自己的电脑上重跑配套 Lab 背后的十二反应网络，用三份已发表的资料核对它，并重建它的核反应率误差预算。只需要 Python 3.9 或更新版本，不需要任何第三方包。请把三份资料放在手边：[Serpico 等人，arXiv:astro-ph/0408076v3](https://arxiv.org/abs/astro-ph/0408076v3) 附录 E 中的反应率拟合，[Pitrou 等人的 PRIMAT 论文，arXiv:1801.08023v3](https://arxiv.org/abs/1801.08023v3) 的表 V，以及 [Fields 等人，arXiv:1912.01132v1](https://arxiv.org/abs/1912.01132v1) 的表 II。 ## 什么是冻结反应网格 一个核合成网络的大部分计算量，花在只依赖温度的量上：热核反应率、中子质子之间的弱反应率，以及光子电子等离子体的能量密度。只要把弱反应率写成以 $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 的质量是转录的数字，在使用它们的文件里注明了出处。 **表 5.1** 文件清单，全部由本站提供。 | 文件 | 作用 | | --- | --- | | physical\_constants.py、nuclear\_data.py | 常数、AME2020 质量过剩和自旋、十二个反应 | | serpico\_rate\_fits.py、rate\_fit\_windows.py | 十一个 Serpico 拟合及其一倍标准差拟合，以及每个拟合可信的温度窗口 | | weak\_rates.py、plasma\_thermodynamics.py、numerical\_quadrature.py | 玻恩近似弱反应率；光子、电子和中微子热力学 | | build\_frozen\_rate\_grid.py | 由上面的模块生成网格 | | network\_solver.py | 只由网格驱动的网络 | | propagate\_rate\_uncertainties.py | 敏感度、一倍标准差预算、蒙特卡罗 | | replay.py、input.json、expected.json | 完整复算脚本、输入和保存的结果 | | teaching\_values.py、teaching-values.json | 课程各章引用的全部数值 | | README.md、provenance.json、LICENSE.txt、write\_manifest.py、MANIFEST.json | 命令和单位、来源、许可，以及每个文件的大小和 SHA-256 摘要 | 下载链接：[MANIFEST.json](https://sandboxphysics.com/zh/stories/primordial-nucleosynthesis/materials/MANIFEST.json)、[input.json](https://sandboxphysics.com/zh/stories/primordial-nucleosynthesis/materials/input.json)、[expected.json](https://sandboxphysics.com/zh/stories/primordial-nucleosynthesis/materials/expected.json)、[frozen-rate-grid.json](https://sandboxphysics.com/zh/stories/primordial-nucleosynthesis/materials/frozen-rate-grid.json)、[teaching-values.json](https://sandboxphysics.com/zh/stories/primordial-nucleosynthesis/materials/teaching-values.json) 和 [provenance.json](https://sandboxphysics.com/zh/stories/primordial-nucleosynthesis/materials/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 年版](https://pdg.lbl.gov/2025/reviews/rpp2025-rev-bbang-nucleosynthesis.pdf) 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$，反应率变为 $$ f(T)\to f(T)\,\bigl[1+\sigma(T)\bigr]^{z},\qquad \sigma(T)=\frac{\delta f}{f}(T). $$ 一倍标准差预算每次只对一个反应取 $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}}$。与观测比较时这些限制很重要。上面每一项检查都是同类相比，结论仍然成立。 ## 遇到下一个网络时的检查清单 - 它用哪套反应率、哪些弱反应率修正、怎样处理中微子脱耦？ - 它是否在舍入误差范围内守恒重子数，并记录了这项检查？ - 它是否与自己的来源表一致，又是否在相同输入下与一个独立计算一致？ - 是否用加倍的时间步重算过？ - 反应率误差是否逐个反应、带温度依赖地传播，并用固定种子的蒙特卡罗检验过？ - 在给出任何新物理限制之前，是否先说明了模型的绝对偏差？ ## 在Lab里试一试 改任何设置前，先核对Lab默认结果和保存的复算结果一致。 [打开轻元素反应网络Lab](https://sandboxphysics.com/zh/labs/big-bang-light-element-network/) 下一篇：[读懂原初核合成对额外辐射和锂问题的证据](https://sandboxphysics.com/zh/stories/primordial-nucleosynthesis/reading-lithium-and-extra-radiation-bounds/) ## 来源 1. [Serpico 等人，arXiv:astro-ph/0408076v3](https://arxiv.org/abs/astro-ph/0408076v3) 2. [Pitrou 等人的 PRIMAT 论文，arXiv:1801.08023v3](https://arxiv.org/abs/1801.08023v3) 3. [Fields 等人，arXiv:1912.01132v1](https://arxiv.org/abs/1912.01132v1) 4. [粒子物理评述 2025 年版](https://pdg.lbl.gov/2025/reviews/rpp2025-rev-bbang-nucleosynthesis.pdf) ## 继续阅读 宇宙最初几分钟与轻元素 · 阅读路径 ### 从第一秒走到轻元素 1. [文章 **温度和膨胀速度怎样决定今天的氦和氘**](https://sandboxphysics.com/zh/stories/primordial-nucleosynthesis/temperature-expansion-and-light-elements/) 2. [文章 **中子和质子什么时候停止互相转变**](https://sandboxphysics.com/zh/stories/primordial-nucleosynthesis/when-neutrons-and-protons-stop-converting/) 3. [文章 **氘为什么要等宇宙冷下来**](https://sandboxphysics.com/zh/stories/primordial-nucleosynthesis/why-deuterium-waits/) 4. [文章 **重子密度和额外辐射怎样移动轻元素丰度**](https://sandboxphysics.com/zh/stories/primordial-nucleosynthesis/baryon-density-radiation-and-abundance-sensitivity/) 5. [当前页面 **用冻结反应网格复算小网络，并传播主要反应率误差**](https://sandboxphysics.com/zh/stories/primordial-nucleosynthesis/recompute-a-frozen-rate-network/) 6. [文章 **读懂原初核合成对额外辐射和锂问题的证据**](https://sandboxphysics.com/zh/stories/primordial-nucleosynthesis/reading-lithium-and-extra-radiation-bounds/)