面向物理系本科一年级 · 可执行教程
教程目标：理解Zhang et al. (2024) 如何用 eRASS:4 的 ~26,099 个 central galaxies 按 stellar mass 分 bin stacking，得到一个 20 kpc M31-mass template，最终变成 Figure 3 上的一个固定点。 Tutorial goal: Understand how Zhang et al. (2024) stack ~26,099 central galaxies from eRASS:4 by stellar mass bins, derive a 20 kpc M31-mass template, and turn it into one fixed point on Figure 3.
目标读者：物理系本科一年级，已学完普通物理（电磁学/光学），了解基本的原子物理概念（能级、跃迁），但不需要天文观测经验。 Target audience: First-year physics undergraduates who have completed general physics (electromagnetism/optics) and basic atomic physics, without requiring astronomical observing experience.
核心问题：eRASS:4 巡天观测了数万个外部星系的 hot CGM。Zhang+24 把 M31-mass bin (11.0 < log(M*/M☉) < 11.25) 的星系堆叠起来，用 projected beta law 拟合得到表面亮度 profile。在 R = 20 kpc 处评价这个 profile，经过距离抵消和 APEC 吸收转换，就得到 Figure 3 上的一个点。这个点不是 M31 的测量——它是 M31-mass 星系的 population stack。
本教程只需要 Python 标准科学栈。数据从 frozen ledger 直接读取——不需要独立的 CSV。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pathlib import Path
plt.style.use('seaborn-v0_8-whitegrid')
plt.rcParams['font.size'] = 12
plt.rcParams['figure.dpi'] = 150
plt.rcParams['savefig.dpi'] = 300
# NOTE: do NOT add a CJK font fallback here — see SKILL.md
# All matplotlib text labels stay English-only.
LEDGER = Path("assets/data/m31_cgmsum_conditional_prior_ledger.csv")
print("Environment ready!")Environment ready!温度在 \(10^6\)–\(10^7\) K 的 circumgalactic medium (CGM) 通过热轫致辐射和线发射发出软 X 射线（0.5–2.0 keV）。但单个外部星系的 CGM 信号极弱——在 eROSITA 的 4 次全天巡天（eRASS:4）中，单次曝光的信号几乎淹没在背景噪声中。
Zhang et al. (2024) 采用了 stacking 策略：
三维等离子体密度常用 beta 模型：
\[ n_e(r) = n_{e,0} \left[1 + \left(\frac{r}{r_c}\right)^2\right]^{-3\beta/2} \]
投影到天空平面上，X 射线表面亮度（假设均匀温度）近似为：
\[ S_X(R) = S_{X,0} \left[1 + \left(\frac{R}{r_c}\right)^2\right]^{-3\beta + 1/2} \]
Zhang+24 对 M31-mass bin 的拟合结果为：
| 参数 | 值 | 含义 |
|---|---|---|
| \(\log S_{X,0}\) | 37.1 | 中心表面亮度（erg s\(^{-1}\) kpc\(^{-2}\)） |
| \(r_c\) | 4 kpc | 核半径 |
| \(\beta\) | 0.37 | 斜率参数 |
graph TD
A["eRASS:4<br/>~26,099 central galaxies"] --> B["按 M* 分 bin<br/>M31-mass: 11.0 < log(M*/Msun) < 11.25"]
B --> C["逐星系提取<br/>X-ray surface brightness profile"]
C --> D["Stack: 对齐中心后叠加"]
D --> E["Projected beta law 拟合<br/>S_X(R) = S_X0 [1+(R/rc)^2]^(-3β+1/2)"]
E --> F["在 R=20 kpc 评价<br/>→ 1.725×10^36 erg s-1 kpc-2"]
F --> G["距离抵消 + 单位换算<br/>+ APEC 吸收转换"]
G --> H["Figure 3: Zhang+2024<br/>M31-mass stack = 0.828"]graph TD
A["eRASS:4<br/>~26,099 central galaxies"] --> B["按 M* 分 bin<br/>M31-mass: 11.0 < log(M*/Msun) < 11.25"]
B --> C["逐星系提取<br/>X-ray surface brightness profile"]
C --> D["Stack: 对齐中心后叠加"]
D --> E["Projected beta law 拟合<br/>S_X(R) = S_X0 [1+(R/rc)^2]^(-3β+1/2)"]
E --> F["在 R=20 kpc 评价<br/>→ 1.725×10^36 erg s-1 kpc-2"]
F --> G["距离抵消 + 单位换算<br/>+ APEC 吸收转换"]
G --> H["Figure 3: Zhang+2024<br/>M31-mass stack = 0.828"]
ledger = pd.read_csv(LEDGER)
zhang = ledger[ledger["prior_id"] == "zhang2024_m31_mass_nominal_20kpc"].iloc[0]
print(f"Loaded: {zhang['prior_id']}")
print(f"Reference: {zhang['reference']}")
print(f"DOI: {zhang['doi']}")
print(f"\nKey values from ledger:")
print(f" original_central = {zhang['original_central']:.6e} {zhang['original_units']}")
print(f" figure_central = {zhang['figure_central']:.15f}")
print(f" figure_low = {zhang['figure_low']:.15f}")
print(f" figure_high = {zhang['figure_high']:.15f}")
print(f" (low = high because this is a fixed nominal point, not a distribution)")
print(f"\nScope: {zhang['scope']}")
print(f"Caveat: {zhang['caveat']}")Loaded: zhang2024_m31_mass_nominal_20kpc
Reference: Zhang et al. (2024)
DOI: 10.1051/0004-6361/202449412
Key values from ledger:
original_central = 1.725299e+36 erg s-1 kpc-2, intrinsic 0.5-2 keV
figure_central = 0.828453687354712
figure_low = 0.828453687354712
figure_high = 0.828453687354712
(low = high because this is a fixed nominal point, not a distribution)
Scope: empirical M31-stellar-mass stack; characteristic 10-30 kpc value
Caveat: A nominal line, not a formal local prior; beta-parameter covariance is unavailable.Zhang+24 的 stacking 结果给出的是整个 profile，不是单个数字。Figure 3 需要一个特征值——在 R = 20 kpc 处评价（这是 M31 观测足迹的典型半径），得到一个名义点的 intrinsic 0.5–2.0 keV 表面亮度。
# Zhang+2024 M31-mass bin beta-profile parameters
log_SX0 = 37.1 # log10 of central normalization
S_X0 = 10**log_SX0 # erg s-1 kpc-2
r_c = 4.0 # core radius (kpc)
beta = 0.37 # beta slope
R_eval = 20.0 # evaluation radius (kpc)
# Projected beta law: S_X(R) = S_X0 * [1 + (R/r_c)^2]^(-3*beta + 1/2)
exponent = -3.0 * beta + 0.5
x = (R_eval / r_c)**2
S_X_20kpc = S_X0 * (1.0 + x)**exponent
print(f"S_X0 = {S_X0:.6e} erg s-1 kpc-2")
print(f"exponent = -3β + 1/2 = {exponent:.2f}")
print(f"(R/r_c)^2 = ({R_eval}/{r_c})^2 = {x:.2f}")
print(f"S_X(R=20 kpc) = {S_X_20kpc:.6e} erg s-1 kpc-2")
# Verify against the frozen ledger
LEDGER_ORIGINAL = zhang["original_central"]
print(f"\nFrozen ledger original_central = {LEDGER_ORIGINAL:.6e}")
# Note: log_SX0=37.1 is rounded from the paper; the exact S_X0
# that reproduces the ledger value is 1.258920e+37 (log ≈ 37.099998).
# The difference is ~5×10^-6 relative — negligible for this tutorial.
np.testing.assert_allclose(S_X_20kpc, LEDGER_ORIGINAL, rtol=1e-4)
print("✓ Matches frozen ledger to < 1e-4")
print(" (log S_X0=37.1 is a 2-significant-digit paper value)")S_X0 = 1.258925e+37 erg s-1 kpc-2
exponent = -3β + 1/2 = -0.61
(R/r_c)^2 = (20.0/4.0)^2 = 25.00
S_X(R=20 kpc) = 1.725306e+36 erg s-1 kpc-2
Frozen ledger original_central = 1.725299e+36
✓ Matches frozen ledger to < 1e-4
(log S_X0=37.1 is a 2-significant-digit paper value)R = np.logspace(-1, 2.2, 200)
S_R = S_X0 * (1.0 + (R/r_c)**2)**exponent
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.loglog(R, S_R, color="#2676b8", linewidth=2, label="Beta profile")
ax.axvline(R_eval, color="#d47a2c", linestyle="--", linewidth=1.5,
label=f"R = {R_eval:.0f} kpc")
ax.scatter([R_eval], [S_X_20kpc], color="#d47a2c", s=80, zorder=5,
label=f"S_X = {S_X_20kpc:.2e}")
# Mark the core radius
ax.axvline(r_c, color="#6f7782", linestyle=":", linewidth=1, alpha=0.7)
ax.annotate(f"$r_c$ = {r_c:.0f} kpc", xy=(r_c, 3e36), fontsize=10,
color="#6f7782")
ax.set_xlabel("Projected radius R (kpc)")
ax.set_ylabel("Surface brightness (erg s^-1 kpc^-2)")
ax.set_title("Zhang+2024 M31-mass stack: projected beta profile")
ax.legend(fontsize=9, loc="lower left")
plt.show()
观察：beta profile 在 \(R > r_c\) 处缓慢下降。R=20 kpc 处的值约为 \(1.73 \times 10^{36}\) erg s\(^{-1}\) kpc\(^{-2}\)——这是 intrinsic（吸收校正前）0.5–2.0 keV 的表面亮度。
对于 resolved 源（我们能分辨其空间结构的源），表面亮度不依赖距离。这是因为：
\[ S = \frac{L}{4\pi D^2}, \quad \Omega = \frac{A}{D^2} \]
所以表面亮度 \(S/\Omega = L/(4\pi A)\)——\(D^2\) 在分子和分母中抵消。Zhang+2024 的 stacking 给出的是 intrinsic luminosity per projected kpc²，需要转换为 flux per solid angle。
# Convert erg s-1 kpc-2 → erg s-1 cm-2 sr-1
# Start from the exact frozen ledger value (not the beta-evaluated approximation)
S_X_orig = zhang["original_central"] # 1.7252986024966184e+36
kpc_to_cm = 3.085677581491368e21 # 1 kpc in cm
# Step 2a: Divide by 4π to get per steradian
S_X_per_sr = S_X_orig / (4.0 * np.pi)
print(f"Step 2a: /4π → {S_X_per_sr:.6e} erg s-1 kpc-2 sr-1")
# Step 2b: Convert kpc^-2 → cm^-2
S_X_per_cm2_sr = S_X_per_sr / (kpc_to_cm**2)
print(f"Step 2b: /kpc_cm^2 → {S_X_per_cm2_sr:.6e} erg s-1 cm-2 sr-1")
# Why this works: S_X is luminosity per projected kpc^2 on the sky.
# Flux per solid angle = S_X / (4π) [erg s-1 kpc-2 sr-1]
# Then convert kpc^-2 → cm^-2 by dividing by (kpc_cm)^2.
# The D^2 cancels because S_X is already a surface brightness (per kpc^2),
# not a total luminosity.Step 2a: /4π → 1.372949e+35 erg s-1 kpc-2 sr-1
Step 2b: /kpc_cm^2 → 1.441960e-08 erg s-1 cm-2 sr-1Figure 3 使用 arcmin⁻² 作为立体角单位（因为 M31 的 CGM 在角分尺度上）。1 球面度（sr）= \((180/\pi \times 60)^2\) arcmin²：
sr_to_arcmin2 = (180.0 / np.pi * 60.0)**2
print(f"1 sr = {sr_to_arcmin2:.2f} arcmin^2")
S_X_per_arcmin2 = S_X_per_cm2_sr / sr_to_arcmin2
print(f"\nStep 3: /sr_to_arcmin2 → {S_X_per_arcmin2:.6e} erg s-1 cm-2 arcmin-2")1 sr = 11818102.86 arcmin^2
Step 3: /sr_to_arcmin2 → 1.220128e-15 erg s-1 cm-2 arcmin-2Figure 3 使用 \(10^{-15}\) erg cm\(^{-2}\) s\(^{-1}\) arcmin\(^{-2}\) 作为统一单位（“flux unit”）。
# Multiply by 1e15 to convert erg → 1e-15 erg
S_X_flux_unit_intrinsic = S_X_per_arcmin2 * 1e15
print(f"Step 4: ×1e15 → {S_X_flux_unit_intrinsic:.6f} flux unit (intrinsic)")Step 4: ×1e15 → 1.220128 flux unit (intrinsic)eRASS 探测的是 0.5–2.0 keV 波段的 X 射线。银河系中性氢（NH）会光致吸收部分光子，尤其软 X 射线（< 1 keV）更易被吸收。Zhang+2024 的 stacking 结果已经是 absorbed（经过银河系吸收后的）值——因为 stacking 用的是观测到的实际计数。
但 Figure 3 上的其他点（如 Locatelli+2024、Grayson+2025）用的是不同的吸收校正链。为了统一比较，所有点都用同一个 APEC 模板做吸收转换。
一个重要的近似：用单一的 v19 APEC 等离子体光谱代替 stack 中星系的多相位分布。真实的 stacking 信号来自不同温度、不同金属丰度的 CGM 混合体，但 Figure 3 的转换链假设可以用一个单一的 APEC 光谱（kT = 0.22 keV, Z = 0.3 solar, Lodders 丰度）来近似。
# v19 APEC absorbed/intrinsic ratio for 0.5-2.0 keV
# Computed with: kT=0.22 keV, Z=0.3 solar, Lodders abundances, NH=6.7×10^20 cm^-2
apec_ratio = 0.678988916
S_X_figure = S_X_flux_unit_intrinsic * apec_ratio
print(f"APEC absorbed/intrinsic ratio = {apec_ratio:.9f}")
print(f"Step 5: ×apec_ratio → {S_X_figure:.15f} flux unit (absorbed)")
# Verify against the frozen ledger
LEDGER_FIGURE = zhang["figure_central"]
print(f"\nFrozen ledger figure_central = {LEDGER_FIGURE:.15f}")
np.testing.assert_allclose(S_X_figure, LEDGER_FIGURE, rtol=1e-10)
print(f"✓ Matches frozen ledger to < 1e-10")APEC absorbed/intrinsic ratio = 0.678988916
Step 5: ×apec_ratio → 0.828453687357045 flux unit (absorbed)
Frozen ledger figure_central = 0.828453687354712
✓ Matches frozen ledger to < 1e-10graph LR
A["Beta profile<br/>S_X(R=20 kpc)<br/>1.725×10^36<br/>erg s-1 kpc-2"] --> B["÷ 4π<br/>→ per sr"]
B --> C["÷ kpc_cm^2<br/>→ cm-2"]
C --> D["÷ 11,818,102.86<br/>→ arcmin-2"]
D --> E["× 10^15<br/>→ flux unit"]
E --> F["× 0.678988916<br/>(v19 APEC)<br/>→ absorbed"]
F --> G["Figure 3<br/>Zhang+2024<br/>M31-mass stack<br/>= 0.828"]graph LR
A["Beta profile<br/>S_X(R=20 kpc)<br/>1.725×10^36<br/>erg s-1 kpc-2"] --> B["÷ 4π<br/>→ per sr"]
B --> C["÷ kpc_cm^2<br/>→ cm-2"]
C --> D["÷ 11,818,102.86<br/>→ arcmin-2"]
D --> E["× 10^15<br/>→ flux unit"]
E --> F["× 0.678988916<br/>(v19 APEC)<br/>→ absorbed"]
F --> G["Figure 3<br/>Zhang+2024<br/>M31-mass stack<br/>= 0.828"]
| 步骤 | 操作 | 值 | 单位 |
|---|---|---|---|
| 0 | Beta profile at R=20 kpc | \(1.725 \times 10^{36}\) | erg s\(^{-1}\) kpc\(^{-2}\) |
| 1 | ÷ 4π | \(1.373 \times 10^{35}\) | erg s\(^{-1}\) kpc\(^{-2}\) sr\(^{-1}\) |
| 2 | ÷ (3.086×10\(^{21}\))\(^2\) | \(1.442 \times 10^{-8}\) | erg s\(^{-1}\) cm\(^{-2}\) sr\(^{-1}\) |
| 3 | ÷ 11,818,102.86 | \(1.220 \times 10^{-15}\) | erg s\(^{-1}\) cm\(^{-2}\) arcmin\(^{-2}\) |
| 4 | × 10\(^{15}\) | 1.220 | flux unit (intrinsic) |
| 5 | × 0.678988916 | 0.828 | flux unit (absorbed) |
| 步骤 | 假设 | 来源 |
|---|---|---|
| Stacking 样本 | ~26,099 central galaxies from eRASS:4 | 原文 |
| M31-mass bin | \(11.0 < \log M_*/M_\odot < 11.25\), median \(z \approx 0.12\) | 原文 |
| 光谱模型 | Beta law 投影到天空 | 原文拟合 |
| 温度/金属丰度 | 单一 v19 APEC (kT=0.22 keV, Z=0.3) | 本项目选择 |
| 吸收柱密度 | NH = \(6.7 \times 10^{20}\) cm\(^{-2}\) (HI4PI) | 本项目选择 |
| 吸收截面库 | APEC 3.0.9 + Lodders 丰度 | 本项目 convention |
| 评价半径 | R = 20 kpc（名义值，不是物理测量） | 本项目选择 |
| 统计性质 | 固定 nominal 点（low = high = central） | 不是置信区间 |
| M31 关系 | 这是 population stack，不是 M31 的测量 | 原文 |
如果我们在 R = 10 kpc 而不是 20 kpc 评价 beta profile，Figure 3 上的值会变成多少？
At R = 10 kpc:
S_X(R) = 3.760052e+36 erg s-1 kpc-2
Figure value = 1.805501 flux unit
(At R = 20 kpc: 0.828454)
Ratio 10/20 kpc = 2.1794在 frozen ledger 中，figure_low == figure_high。为什么这个点没有误差棒？
Answer: This is a fixed nominal point, not a statistical estimate.
The beta-profile parameters (log S_X0, r_c, beta) have uncertainties,
but the covariance between them is not published in Zhang+2024.
Without the covariance matrix, we cannot propagate errors to R=20 kpc.
The Figure 3 point is therefore a nominal template, not a confidence interval.手动计算 \(1\ \mathrm{kpc} = 3.085677581491368 \times 10^{21}\ \mathrm{cm}\)。用这个值验证步骤2的转换。
#| label: ex3
#| echo: false
kpc_cm_check = 1e3 * 3.085677581491368e18 # 1 kpc = 1000 pc, 1 pc = 3.085677581491368e18 cm
print(f"1 kpc = {kpc_cm_check:.6e} cm")
print(f"Stored = {kpc_to_cm:.6e} cm")
assert abs(kpc_cm_check - kpc_to_cm) < 1e10
print("✓ Consistent")graph TD
A["eRASS:4<br/>~26,099 central galaxies"] --> B["M31-mass bin<br/>11.0 < log M*/Msun < 11.25"]
B --> C["Stack X-ray profiles<br/>+ fit beta law"]
C --> D["S_X(20 kpc)<br/>= 1.725×10^36 erg s-1 kpc-2"]
D --> E["Distance cancellation<br/>÷4π, ÷kpc_cm^2"]
E --> F["Solid angle conversion<br/>÷11,818,102.86 arcmin^2/sr"]
F --> G["Flux unit<br/>×10^15"]
G --> H["APEC absorption<br/>×0.678988916"]
H --> I["Figure 3<br/>Zhang+2024 M31-mass stack<br/>= 0.828"]graph TD
A["eRASS:4<br/>~26,099 central galaxies"] --> B["M31-mass bin<br/>11.0 < log M*/Msun < 11.25"]
B --> C["Stack X-ray profiles<br/>+ fit beta law"]
C --> D["S_X(20 kpc)<br/>= 1.725×10^36 erg s-1 kpc-2"]
D --> E["Distance cancellation<br/>÷4π, ÷kpc_cm^2"]
E --> F["Solid angle conversion<br/>÷11,818,102.86 arcmin^2/sr"]
F --> G["Flux unit<br/>×10^15"]
G --> H["APEC absorption<br/>×0.678988916"]
H --> I["Figure 3<br/>Zhang+2024 M31-mass stack<br/>= 0.828"]
一句话总结：Zhang+2024 把 eRASS:4 中 ~26,099 个 M31-mass 星系堆叠在一起，用 beta law 拟合 X 射线表面亮度，在 R = 20 kpc 评价。经过距离抵消、单位换算和 v19 APEC 吸收转换后，得到 0.828 flux unit——这是 M31-mass 星系（不是 M31 本身！）的 population stack 模板，作为一个固定点出现在 Figure 3 上。
教程结束 🎓 下一步：继续阅读 Grayson+2025 EAGLE/SIMBA tutorial，了解宇宙学模拟如何预测 M31-mass 星系的 CGM。