用保存的输入复算一条星系团视线

先下载三个文件：保存的输入、脚本和保存的输出。脚本只需要Python 3和标准库；保存的输出是用Python 3.9.6生成的。物理内容见把光子和轴子看成两个会混合的态里的两态混合，以及随机磁畴、偏振与恒星冷却里的随机磁畴模型。

第一步：读输入和单位

先打开输入文件。units一栏在任何数字出现之前就把单位定死：光子能量用keV，耦合用GeV的倒数，磁场用微高斯，电子密度用每立方厘米，轴子质量用eV，长度用kpc，角度用弧度。从0.6到9 keV的七个能量只列一次，所有情形共用。

一共三种情形。均匀情形是一个磁场方向，分成12块，每块8 kpc。缠绕情形就是Lab的默认视线，24块，角度固定但不规则。共振情形把轴子质量设成气体的等离子体频率。每种情形都把各块角度明确写出来。Lab的角度来自一个带种子的随机数发生器；把角度存下来，复算时就不必重现那个发生器。

第二步：运行脚本

python3 -I recompute_sightline.py sightline-inputs.json my-outputs.json

-I参数让Python不从当前文件夹导入模块。把你得到的文件和保存的输出对比，应当一直吻合到最后几位。Lab自己的传递矩阵用另一种语言、另一种写法计算，与这些输出的存活概率相差不到 $10^{-10}$ ，这一点由仓库里的测试检查。

第三步：一块磁畴做了什么

每一块先把态转到磁场坐标系，在平行光子和轴子那个二阶矩阵的本征基里演化，给垂直偏振的光子加上等离子体相位，再转回来。核心函数如下：

def evolve_cell(state, plasma, axion, mixing, length, angle):
    ax, ay, a = state
    c, s = math.cos(angle), math.sin(angle)
    parallel = c * ax + s * ay
    perpendicular = (-s * ax + c * ay) * cmath.exp(1j * plasma * length)
    theta = 0.5 * math.atan2(2 * mixing, plasma - axion)
    mean = 0.5 * (plasma + axion)
    half_split = 0.5 * math.hypot(plasma - axion, 2 * mixing)
    ct, st = math.cos(theta), math.sin(theta)
    upper = (ct * parallel + st * a) * cmath.exp(1j * (mean + half_split) * length)
    lower = (-st * parallel + ct * a) * cmath.exp(1j * (mean - half_split) * length)
    parallel, a = ct * upper - st * lower, st * upper + ct * lower
    return [c * parallel - s * perpendicular, s * parallel + c * perpendicular, a]

参数plasma、axion和mixing分别是以千秒差距倒数为单位的 $\Delta_{\rm pl}$ 、 $\Delta_a$ 和 $\Delta_{a\gamma}$ 。脚本从CODATA常数推出它们：一千秒差距是 $1.5637\times10^{26}\,\mathrm{eV^{-1}}$ ，一特斯拉是 $195.35\,\mathrm{eV^2}$ ，每立方厘米一个电子时等离子体频率是 $3.713\times10^{-11}$ eV。这三个值写在输出的derivedConstants一栏里，读者能看出某个文件是用哪些常数算的。

非偏振光的处理方式是：两种纯偏振分别跑一遍，再把它们的密度矩阵各按一半权重平均。线偏振度为

\Pi_L=\frac{\sqrt{(\rho_{xx}-\rho_{yy})^2+4(\mathrm{Re}\,\rho_{xy})^2}}{\rho_{xx}+\rho_{yy}}.

第四步：必须通过的检查

概率守恒。 光子存活率加轴子比例必须等于一。21行保存结果里最大的偏离是 $3.6\times10^{-15}$ ，只是舍入误差。

均匀情形对公式。 一段均匀磁场里，非偏振光的存活率是 $1-P_{\gamma\to a}/2$ ，其中用两态概率。所有均匀情形的行都与这个闭式吻合到约 $10^{-15}$ 。6.6 keV处两者都是0.5572；3 keV处都是0.9729，正好落在振荡的一个节点附近。

共振时的平坦。 轴子质量等于等离子体频率时，七个能量上的存活率都是0.7761。

蒙特卡罗平均。 缠绕情形保存的角度只是一种实现。为了检验系综规律，脚本用Python自己的随机数发生器、种子79，重新抽2000组24个均匀角度，对3 keV处的轴子比例取平均：

**表 5.1** 3 keV处蒙特卡罗平均与随机角度规律的比较。
量	数值
2000个实现的平均轴子比例	$0.2434$
这个平均值的标准误差	$0.0027$
24块的乘积规律	$0.2399$
差距，以标准误差计	$1.3$

差1.3个标准误差，和随机涨落相符。如果超过三左右，就说明传播计算或公式里有错。

第五步：四个失败情形

输出里还存了四个故意犯的错，每个都配一个数字说明损害有多大。

**表 5.2** 保存的失败情形及其说明。
错误	测得的后果	教训
共振处把等离子体项的符号弄反	存活率随能量的变化幅度从 $2\times10^{-15}$ 升到 $0.013$	共振要求 $\Delta_{\rm pl}$ 和 $\Delta_a$ 同号
只取12个能量而不是400个	插值最多差 $0.28$	单个实现的起伏很窄
在一块强磁畴里用16步前向欧拉	总概率变成 $1.055$	欧拉步不是幺正的
把一条视线当成平均	0.5到10 keV内最大差距 $0.22$	实现只能和实现比

实践中最常见的是欧拉那一条。任何不保持模长的积分器都会慢慢凭空造出或吞掉光子，误差随路径变长而增大。精确的单块传递矩阵没有这个问题。仓库测试里还有一个独立的四阶龙格–库塔积分器，步长2 pc，与它吻合到 $10^{-6}$ 。

第六步：复算一条从Lab导出的视线

Lab的导出按钮会写出一个格式为sandboxphysics.axion-photon-conversion.v1的JSON记录。它的input一栏把耦合、密度和质量存成对数，键名形如logCouplingPerGev。每个都取十的该次幂换回来，domains列表原样照抄，然后作为新情形加进你的输入文件。导出记录里还分开存了spectrum.simulatedSurvival和spectrum.idealReference。把复算的结果和模拟曲线比较；理想参照是另一样东西：均匀模式下是闭式公式，缠绕模式下是系综平均。

第七步：带着前提去读冻结的上限

Lab的恒星标尺上标了四条已发表的上限。它们照原文抄录，不在这里重算，所以每一条都带着作者的前提。

**表 5.3** Lab中显示的已发表上限及其依赖的前提。
上限	数值	适用范围	依赖
NGC 1275的Chandra光谱	$7.9\times10^{-13}\,\mathrm{GeV^{-1}}$ ，99.7%	$m_a<10^{-12}$ eV	英仙座磁场模型B
NGC 1275的Fermi-LAT光谱	$5\times10^{-12}\,\mathrm{GeV^{-1}}$ ，95%	0.5到5 neV	六年LAT数据和假定的英仙座磁场模型
CAST太阳轴子望远镜	$5.8\times10^{-11}\,\mathrm{GeV^{-1}}$ ，95%	$m_a<0.02$ eV	标准太阳模型给出的太阳轴子流量
水平分支星计数	$6.6\times10^{-11}\,\mathrm{GeV^{-1}}$ ，95%	到几keV	恒星模型和氦丰度

X射线的值是Reynolds等2020年论文的式(5)。伽马射线的值来自Ajello等2016年论文的摘要，我们没有通读全文。CAST和恒星的值分别是PDG轴子综述2025的式(89.20)和式(89.16)。Lab里恒星上限的质量边界保守地取1 keV。Ayala等2014年论文讨论几keV以下的粒子，但没有印出明确的截断。

置信度不同，95%和99.7%，所以这些数不能直接互换。Reynolds等人对磁场模型A给出 $6.3\times10^{-13}\,\mathrm{GeV^{-1}}$ 。光是换磁场模型，X射线上限就会移动百分之几十。

下一条视线的检查清单

每个输入都定好单位，并写进文件。
明确保存磁场角度，而不只存一个种子。
每一行都检查归一化。
用均匀情形核对闭式公式，用共振情形核对平坦。
把许多实现的平均和系综规律比较，报告偏离了几个标准误差。
能量取样要足够密，能分辨单个实现的起伏。
每条曲线都标明是模拟、理想参照还是观测；已发表的上限要连同置信度、质量范围和磁场模型一起保留。

到Lab里试试：从Lab导出一条视线，与保存的输出对照。