---
title: "亲手复算一段电离历史"
description: "从原始单元数组出发，检查光子收支，再用独立积分和傅里叶投影复算电离历史与空间起伏。"
canonical_url: "https://sandboxphysics.com/zh/stories/reionization/recompute-a-fixed-source-history/"
language: "zh-CN"
page_type: "story"
alternate_languages:
  en: "https://sandboxphysics.com/stories/reionization/recompute-a-fixed-source-history/"
last_modified: "2026-11-12"
---

[星光怎样改变星系间的氢](https://sandboxphysics.com/zh/stories/reionization/)

# 亲手复算一段电离历史

Sandbox Physics · 2026-11-12

复算先要认准同一份物理问题：哪些源在发光，光能逃出多少，氢气怎样复合，以及数组对应哪个时刻。文件能打开只是第一步。接着要把单元中的电离体积加起来，再从另一条数学路线计算这段历史。

以下步骤使用配套的[固定源数据与复算脚本](https://sandboxphysics.com/zh/stories/reionization/materials/fixed-source-analysis.zip)。先读[电离体积怎样变成亮温和功率](https://sandboxphysics.com/zh/stories/reionization/reading-a-neutral-hydrogen-map/)，再解压数据包。根目录的 records.json 列出五个预设的清单和物理输入；每份记录保存六十五个原生时刻。

进入解压后的 source 目录，把默认记录的路径记在两个环境变量中：

```
export AC48_LOCAL_MANIFEST="../records/reionization-default-a512-20261112/native-manifest.json"
export AC48_LOCAL_INPUTS="../records/reionization-default-a512-20261112/declared-inputs.json"
```

只读数组和显式投影的脚本使用 Python 标准库，不需要 NumPy 或 SciPy。后面的完整模型复算使用固定的 Python 3.12.12、NumPy 2.3.5 和 SciPy 1.17.1；数据包的 requirements.txt 记录了数值库版本。代码与原创合成数据的许可分别在包内列出，第三方软件和原论文图像没有捆绑在其中。

## 先看输入，再看一幅漂亮的图

记录把物理输入、时钟、单元数组和计算结果分开存放。先读源的位置与本征光子率，再读逃逸比例和共同复合参数。位置按周期盒边长的比例表示；盒子边长20 cMpc，没有乘或除小写哈勃参数。反应网络只计算氢，氢质量分数仍采用0.76。

同一组源的位置变了，总本征光子率和共同复合参数若保持原值，这个均匀模型的全局电离历史应当相同；电离区怎样分布仍会改变。逃逸比例变了，进入盒子的光子数也跟着变。这个预测可以在运行之前写下来，随后用两个记录检查。[全局电离体积模型的来源](https://arxiv.org/abs/astro-ph/9809058v1)

清单中有一个容易混淆的量：单元数组保存的是电离相占这个单元多少体积。记作 $f_j$。它的等体积平均值

$$
\overline f=\frac{1}{N^3}\sum_j f_j
$$

应当与同一时刻独立推进的电离体积分数 $Q$ 相符。这里 $N$ 是每条边的单元数。比较时，把这两项都对应到同一原生时刻。

## 读出原来的八字节数

每个原始数占八字节，使用小端序Float64。文件的排列是横向索引最快，然后是纵向，再到深度；写成三维数组时，轴顺序是深度、纵向、横向。这个约定决定某个源、探针和切片到底落在哪个单元。

压缩文件与解压后的数组各有一份SHA-256。它是文件字节的指纹。指纹吻合，说明读到了清单指定的那份数据；随后还要做物理和数学检查。压缩层只拆一次，保留原始双精度值再求平均，先不要把它们改成绘图用的低精度纹理。

从脚本所在目录，检查一个原生时刻：

```
python -S inspect_saved_record.py --manifest ../records/reionization-default-a512-20261112/native-manifest.json --epoch 24 --with-modes --output inspected-frame.json
```

--epoch 指原生索引，范围为零到六十四，索引本身没有时间单位。脚本报告原始文件校验、数组均值与同帧计数之差，以及显式傅里叶投影之差。它还保留源记录原来的资格说明。省略 --epoch 和 --with-modes，就会读取全部六十五帧，先检查每帧的字节与收支。

默认记录的第 24 帧给出 $Q=0.458014596$。同一物理时刻，低逃逸预设为 $Q=0.229007298$，高复合预设为 $Q=0.370305927$。分散源预设与默认记录的全局计数相同，空间排列可以不同。用这个对照检查读到的记录，也能把前面写下的预测接回实际数字。

## 四个计数怎样对应一幅图

用每个氢核来归一化，把留在电离相的数量、累计复合次数、进入盒子的光子数和饱和后剩余光子数分别记为 $q,r,e,u$。它们应满足

$$
e=q+r+u.
$$

这个模型直接推进四个计数，再检查它们之间的关系。若先把 $r$ 定义成 $e-q-u$，这条等式每次都会成立，就无法发现推进过程中多算或漏算的光子。

举一个单独的算术例子：进入盒子的光子有每氢核1.2个，累计复合0.3次，剩余0.1个，保留在电离相的就是每氢核0.8个。均匀密度假设下，它也对应0.8的电离体积分数。原始单元平均值还要分别和直接积分得到的 $q$ 比较；这一步把空间图与光子收支接起来。

## 用积分因子重新算一次

在达到完全电离之前，共同复合模型写成

$$
\frac{\mathrm dQ}{\mathrm dt}=s(t)-b(t)Q,\qquad
b(t)=\alpha_{\mathrm B}C\bar n_{\mathrm H,0}a^{-3}.
$$

$s$ 是每个氢核每秒得到的逃逸光子数，$b$ 的单位也是每秒。复合系数 $\alpha_{\mathrm B}$、共同参数 $C$ 和适当时刻的实际氢密度决定复合损失。

若初始电离体积为零，可以把某时刻新提供的电离贡献乘上此后仍然保留的比例，再把所有时刻相加：

$$
Q(t)=\int_{t_{\mathrm{on}}}^{t}
s(t')\exp\!\left[-\int_{t'}^{t}b(t'')\,\mathrm dt''\right]\mathrm dt'.
$$

这就是积分因子解。指数把共同复合损失积累起来。它与逐步推进候选微分方程使用不同的计算路线，适合作为独立参照。达到电离上限后，还需另外处理维持电离所消耗的光子和剩余光子；不能继续把这个无上限表达式当成体积分数。

独立参考类的输入仍是同一份声明的物理量。下面把参考计数先保存下来，再和原记录比较：

```
import json
import os
import warnings
import numpy as np
from scipy.integrate import IntegrationWarning
from independent.scalar_quadrature import QuadratureReference

with open(os.environ["AC48_LOCAL_INPUTS"]) as source:
    inputs = json.load(source)
with open(os.environ["AC48_LOCAL_MANIFEST"]) as source:
    record = json.load(source)
reference = QuadratureReference(
    inputs,
    record["settings"]["escapeFraction"],
    record["settings"]["commonRecombinationClumping"],
)
epochs = np.array([row["scaleFactor"] for row in record["epochs"]])
assert epochs[0] == reference.initial_scale
with warnings.catch_warnings():
    warnings.simplefilter("error", IntegrationWarning)
    counters, cap_scale_factor = reference.history(epochs)
with open("independent-counters.npz", "xb") as output:
    np.savez(output, scaleFactors=epochs, fourCounters=counters)
```

参考代码计算自己的背景密度、时间关系和积分，导入图与候选推进器分开。积分发出警告时应停下来检查，不能把一次未完成的积分当成通过。保留独立数值之后，再比较四列；归一化计数的原定允许误差是 $5\times10^{-8}$。这个数是检验门槛，具体误差要从真实运行报告读取。[复合系数的原始拟合与约定](https://arxiv.org/abs/astro-ph/9612232v1)

若要从输入重算整条候选历史，在固定版本环境中运行：

```
python recompute_history.py --manifest ../records/reionization-default-a512-20261112/native-manifest.json --physical-input ../records/reionization-default-a512-20261112/declared-inputs.json --output recomputed-history
```

这会重新推进四个计数和全部五百一十二步空间分配，保存原始返回场、各单元还能容纳的电离体积与过量电离量，再逐字比较六十五张原生场。每个数的八字节都相同，才算复现这份固定记录；参考计数则按上面的独立积分路线检查。脚本拒绝覆盖已有输出目录。

这份最终数据包实际解压到新目录后，用这条方法完整复算了五个预设：共三百二十五帧的原始数均一致，独立计数的最大绝对差为 $4.44\times10^{-15}$，低于原定 $5\times10^{-8}$ 门槛。包内证据是标明来源的摘要，完整原始收据另行保存在科学档案中。这个结果检查了固定输入的复算；共同复合、均匀密度与非因果分配的物理近似仍然保留。

时间也值得单独检查。尺度因子由今天归一化为一，经过的秒数却从源在红移十开启时算起。两者通过 $\mathrm dt=\mathrm da/[aH(a)]$ 联系。每年取 31557600 秒时，一段 $3.15576\times10^{15}\,\mathrm s$ 的时间是100 Myr。上面是单位换算练习；每个原生帧的时间仍以清单中的实际秒数为准。

## 守恒之后，再问形状

光子数吻合，仍有许多不同的空间分布可以满足它。增量分配模型会受网格与剩余电离量处理顺序影响。独立检查需要分别改变步长、顺序和空间分辨率，再把它们的误差列开。[光子守恒空间分配方法的原始研究](https://arxiv.org/abs/1807.00836v2)

单个源提供一个可以想清楚的几何参照。均匀气体中，独立电离体积给出等体积球半径

$$
R=L\left(\frac{3Q}{4\pi}\right)^{1/3}.
$$

只在半径小于半个盒子、没有周期绕回且其余均匀近似适用时，用无限空间球与原始单元比较。先把连续球在每个单元内所占体积算准，再看候选图像的误差。把参考积分分得更细，是检查参考本身；把物理网格分得更细，是检查空间计算。这两步回答不同的问题。

最后用同一幅原始场投影到正弦、余弦波。脚本的显式投影逐个单元相加，不调用FFT。系数对应清单里的三十二个带方向模式；四组平均会进一步丢掉方向和相位信息。保留原始场与复系数，可以查出只看四个功率数时漏掉的问题。

## 练习和解释

1. 一幅原始图的所有单元都没变，只把横向轴和深度轴互换。均值会发现错误吗？四组按对称方向平均的功率一定能发现吗？还可检查什么？
2. 候选程序用微分方程推进 $Q$，另一个程序调用候选程序得到每个时刻的导数，再用更小步长推进。它检验了什么？为什么还需要积分因子参考？
3. 把两个源移得很远，总发光率、逃逸比例和共同复合参数保持原值。预测全局计数与空间功率怎样变化。若改成真正逐点计算密度相关复合，哪一步推论需要重做？

第一题的均值保持原样，对称模式组功率也可能保持。可以放一个坐标已知的单元标记，或检查带方向标签的复系数与源、探针位置。数据排列属于物理对应关系的一部分。

第二题可以检查数值步长与推进精度，但两个程序仍共享同一份导数代码。独立积分因子计算还会检查这份导数怎样连接源、密度和时间。独立方法的作用在于多检查一条可能出错的联系。

第三题在当前共同预算模型下，全局四个计数相同，空间排列与功率可以改变。逐点复合若依赖局部密度和电离状态，源的位置也可能改变累计复合量，需要重新计算完整源与复合的收支。保存清楚的模型条件，才能知道一次输入变化应当改变哪些结果。

## 来源

1. [全局电离体积模型的来源](https://arxiv.org/abs/astro-ph/9809058v1)
2. [复合系数的原始拟合与约定](https://arxiv.org/abs/astro-ph/9612232v1)
3. [光子守恒空间分配方法的原始研究](https://arxiv.org/abs/1807.00836v2)

## 继续阅读

星光怎样改变星系间的氢 · 阅读路径

### 从光子、电离体积到无线电信号

1. [文章 **星系的光怎样改变宇宙中的氢**](https://sandboxphysics.com/zh/stories/reionization/starlight-and-hydrogen/)
2. [文章 **电离一半究竟是什么意思**](https://sandboxphysics.com/zh/stories/reionization/ionized-volume-and-local-fraction/)
3. [文章 **一片膨胀的气体需要多少光**](https://sandboxphysics.com/zh/stories/reionization/photons-in-an-expanding-box/)
4. [文章 **两个电离泡相遇以后，光去了哪里**](https://sandboxphysics.com/zh/stories/reionization/when-ionized-bubbles-overlap/)
5. [文章 **从中性氢的信号读出电离泡**](https://sandboxphysics.com/zh/stories/reionization/reading-a-neutral-hydrogen-map/)
6. [当前页面 **亲手复算一段电离历史**](https://sandboxphysics.com/zh/stories/reionization/recompute-a-fixed-source-history/)
7. [文章 **怎样读懂HERA给出的21厘米上限**](https://sandboxphysics.com/zh/stories/reionization/reading-hera-upper-limits/)