面向物理系本科一年级 · 可执行教程
教程目标：理解Henley & Shelton (2013) 的86条高银纬视线如何经过吸收转换，变成Figure 3上的一个population median点。 Tutorial goal: Understand how 86 high-latitude sightlines from Henley & Shelton (2013) are absorption-converted into one population median 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.
核心问题：银河系的高银纬热气体在不同方向上亮度不同——论文发表了intrinsic（吸收校正后）的0.5–2.0 keV表面亮度。Figure 3需要的是absorbed（经过银河系中性氢吸收后的）值。这86条视线如何从intrinsic变absorbed？86个值又如何浓缩成一个median点？
本教程只需要 Python 标准科学栈。所有数据文件已随教程提供。
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.
DATA = Path("assets")
print("Environment ready!")Environment ready!温度在 \(10^6\)–\(10^7\) K 的稀薄气体（密度约 \(10^{-3}\)–\(10^{-4}\) cm\(^{-3}\)）主要通过热轫致辐射（thermal bremsstrahlung）和线发射（line emission）发出软X射线（约0.1–2 keV）。Henley & Shelton (2013) 用XMM-Newton卫星观测了110条高银纬（\(|b| > 30^\circ\)）视线上的 diffuse X-ray background，并用Raymond-Smith等离子体模型拟合每条视线上的银河系热气体（Galactic halo）成分。
X射线光子穿过银河系中性气体时，会被原子（主要是氢、氦和少量重元素）光致电离吸收。能量越低的光子越容易被吸收，因此软X射线（< 1 keV）受到的吸收最显著。
吸收强度用等效中性氢柱密度 \(N_\mathrm{H}\)（单位 cm\(^{-2}\)）衡量。Henley & Shelton 使用 Leiden/Argentine/Bonn (LAB) 全天21 cm巡天数据获取每条视线的 \(N_\mathrm{H}\)。
论文发表的是 intrinsic（吸收校正后）的0.5–2.0 keV表面亮度——即”假设没有银河系吸收时，该视线上的热气体应该有多亮”。而Figure 3需要的是 absorbed（经过吸收后实际到达我们的）值，这样才能与M31方向的实测值比较。
对每条视线 \(i\)：
\[ S_{\mathrm{abs},i} = S_{\mathrm{intr},i} \times T_{0.5-2.0}(kT_i, N_{\mathrm{H},i}) \]
其中 \(T_{0.5-2.0}(kT, N_\mathrm{H})\) 是该视线温度 \(kT\) 和柱密度 \(N_\mathrm{H}\) 下，Raymond-Smith等离子体的0.5–2.0 keV absorbed-to-intrinsic flux ratio。
graph TD
A["86条高纬视线<br/>86 high-latitude sightlines"] --> B["每条视线有自己的<br/>kT 和 N_H"]
B --> C["Raymond-Smith × phabs<br/>计算 T_0.5-2.0(kT, N_H)"]
C --> D["S_abs = S_intr × T"]
D --> E["deg⁻² → arcmin⁻²"]
E --> F["取 median + IQR<br/>= Figure 3 上的一个点"]graph TD
A["86条高纬视线<br/>86 high-latitude sightlines"] --> B["每条视线有自己的<br/>kT 和 N_H"]
B --> C["Raymond-Smith × phabs<br/>计算 T_0.5-2.0(kT, N_H)"]
C --> D["S_abs = S_intr × T"]
D --> E["deg⁻² → arcmin⁻²"]
E --> F["取 median + IQR<br/>= Figure 3 上的一个点"]
henley = pd.read_csv(DATA / "m31_cgmsum_henley_shelton2013_absorbed_0p5_2p0.csv")
print(f"Loaded {len(henley)} sightlines")
print(f"Columns: {list(henley.columns)}")
henley.head(8)Loaded 86 sightlines
Columns: ['seq', 'obsid', 'glon_deg', 'glat_deg', 'nhi_1e20_cm-2', 'temperature_1e6_K', 'temperature_fixed_flag', 'surface_brightness_intrinsic_0p5_2p0_1e12', 'surface_brightness_limit', 'absorbed_to_intrinsic_0p5_2p0_ratio', 'absorbed_0p5_2p0_fluxunit']| seq | obsid | glon_deg | glat_deg | nhi_1e20_cm-2 | temperature_1e6_K | temperature_fixed_flag | surface_brightness_intrinsic_0p5_2p0_1e12 | surface_brightness_limit | absorbed_to_intrinsic_0p5_2p0_ratio | absorbed_0p5_2p0_fluxunit | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 302580501.0 | 5.676000 | -77.682999 | 1.18 | 1.77 | NaN | 3.58 | NaN | 0.914851 | 0.909769 |
| 1 | 2 | 550460801.0 | 6.479000 | -55.914001 | 1.39 | 2.42 | NaN | 2.14 | NaN | 0.914395 | 0.543557 |
| 2 | 3 | 304080501.0 | 12.896000 | -78.051002 | 1.44 | 2.10 | NaN | 4.21 | NaN | 0.903952 | 1.057121 |
| 3 | 4 | 556210301.0 | 54.428001 | 42.050999 | 1.70 | 3.42 | NaN | 1.00 | NaN | 0.915718 | 0.254366 |
| 4 | 5 | 505880101.0 | 55.563999 | 42.457001 | 1.73 | 2.11 | NaN | 1.84 | NaN | 0.886193 | 0.452943 |
| 5 | 6 | 505010501.0 | 56.664001 | 43.826000 | 1.43 | 1.76 | NaN | 1.70 | NaN | 0.897609 | 0.423871 |
| 6 | 7 | 504100401.0 | 62.380001 | 65.692001 | 1.06 | 2.18 | NaN | 2.37 | NaN | 0.929956 | 0.612221 |
| 7 | 8 | 556212601.0 | 63.987999 | 42.877998 | 1.07 | 2.15 | NaN | 1.36 | NaN | 0.928713 | 0.350847 |
数据包含86条视线（原文110条中有87条halo detections，排除1条无唯一LAB \(N_\mathrm{H}\)的Seq=103）。每条视线有：
| 列名 | 含义 |
|---|---|
seq |
原文序号 |
glon_deg, glat_deg |
银道坐标 |
nhi_1e20_cm-2 |
LAB \(N_\mathrm{H}\)（\(10^{20}\) cm\(^{-2}\)） |
temperature_1e6_K |
拟合的halo温度（\(10^6\) K） |
surface_brightness_intrinsic_0p5_2p0_1e12 |
intrinsic 0.5–2.0 keV亮度（\(10^{-12}\) erg cm\(^{-2}\) s\(^{-1}\) deg\(^{-2}\)） |
absorbed_to_intrinsic_0p5_2p0_ratio |
已算好的 \(T_{0.5-2.0}\) 转换比 |
absorbed_0p5_2p0_fluxunit |
absorbed值（Figure 3单位：\(10^{-15}\) erg cm\(^{-2}\) s\(^{-1}\) arcmin\(^{-2}\)） |
每条视线有自己的温度和吸收柱密度，因此转换比 \(T\) 各不相同。温度更高的等离子体在0.5–2.0 keV波段有更大比例的硬光子，受吸收影响较小；\(N_\mathrm{H}\) 更大的视线则吸收更强。
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4.5), constrained_layout=True)
# Left: ratio vs temperature, colored by N_H
sc = ax1.scatter(
henley["temperature_1e6_K"],
henley["absorbed_to_intrinsic_0p5_2p0_ratio"],
c=henley["nhi_1e20_cm-2"],
cmap="viridis",
s=40,
edgecolor="none",
alpha=0.85,
)
cbar = fig.colorbar(sc, ax=ax1)
cbar.set_label("N_H (10^20 cm^-2)")
ax1.set_xlabel("Temperature (10^6 K)")
ax1.set_ylabel("Absorbed / Intrinsic ratio")
ax1.set_title("Ratio vs Temperature")
# Right: intrinsic vs absorbed, 1:1 line
intr = henley["surface_brightness_intrinsic_0p5_2p0_1e12"] / 3.6 # deg^-2 → arcmin^-2, 1e-12 → fluxunit
absd = henley["absorbed_0p5_2p0_fluxunit"]
limit = max(intr.max(), absd.max()) * 1.05
ax2.scatter(intr, absd, c=henley["nhi_1e20_cm-2"], cmap="viridis", s=40, edgecolor="none", alpha=0.85)
ax2.plot([0, limit], [0, limit], color="#6f7782", linestyle="--", linewidth=0.8, label="1:1 (no absorption)")
ax2.set_xlabel("Intrinsic brightness (flux unit)")
ax2.set_ylabel("Absorbed brightness (flux unit)")
ax2.set_title("Intrinsic vs Absorbed")
ax2.legend(fontsize=9)
plt.suptitle("Henley & Shelton 2013: 86 sightline conversion", fontsize=13, fontweight="bold")
plt.show()
观察： - 几乎所有点都在1:1线下方——吸收使观测到的亮度低于intrinsic值。 - 温度较高（>2 MK）的视线转换比更接近1，因为它们的发射谱更硬。 - \(N_\mathrm{H}\) 较大的视线（黄色）转换比更低。
Henley & Shelton的86条视线不指向M31方向（M31在 \(b \approx -21.6^\circ\)，不在样本的 \(|b|>30^\circ\) 域内）。这些视线分布在高银纬天空的不同方向，银河系热气体在不同方向的亮度有天然散布。取median（中位数）比取mean更能抵抗少数极亮或极暗视线的拉扯，因此更适合作为”population context”的代表值。
absd = henley["absorbed_0p5_2p0_fluxunit"].to_numpy()
quantiles = np.quantile(absd, [0.16, 0.25, 0.5, 0.75, 0.84])
print(f"Median (50th): {quantiles[2]:.6f}")
print(f"IQR (25-75): [{quantiles[1]:.6f}, {quantiles[3]:.6f}]")
print(f"p16-p84: [{quantiles[0]:.6f}, {quantiles[4]:.6f}]")
# Verify against the frozen Figure 3 ledger
LEDGER_MEDIAN = 0.38461005704507567
LEDGER_IQR_LO = 0.27354775115143726
LEDGER_IQR_HI = 0.5490818824107443
np.testing.assert_allclose(quantiles[2], LEDGER_MEDIAN, rtol=0, atol=1e-6)
np.testing.assert_allclose(quantiles[1], LEDGER_IQR_LO, rtol=0, atol=1e-6)
np.testing.assert_allclose(quantiles[3], LEDGER_IQR_HI, rtol=0, atol=1e-6)
print("✓ Matches frozen Figure 3 ledger to < 1e-6")Median (50th): 0.384610
IQR (25-75): [0.273548, 0.549082]
p16-p84: [0.246097, 0.690362]
✓ Matches frozen Figure 3 ledger to < 1e-6fig, ax = plt.subplots(figsize=(8, 4.5))
ax.hist(absd, bins=18, color="#2676b8", alpha=0.7, edgecolor="white", linewidth=0.5)
ax.axvline(quantiles[2], color="#d47a2c", linewidth=2, label=f"Median = {quantiles[2]:.3f}")
ax.axvspan(quantiles[1], quantiles[3], color="#d47a2c", alpha=0.12, label=f"IQR [{quantiles[1]:.3f}, {quantiles[3]:.3f}]")
ax.axvline(0.964727, color="#b74842", linestyle="--", linewidth=1.5, label="Observed CGMsum = 0.965 (M31)")
ax.set_xlabel("Absorbed 0.5-2.0 keV brightness (10^-15 erg s^-1 cm^-2 arcmin^-2)")
ax.set_ylabel("Number of sightlines")
ax.set_title("H&S13: 86 sightline absorbed brightness distribution")
ax.legend(fontsize=9)
plt.show()
关键观察： - 86条视线的absorbed亮度散布很宽——从约0.1到约1.1 flux unit。 - M31方向实测值（0.965）落在分布的上尾部——这说明M31方向的银河系前景可能比高银纬典型值更亮，但H&S13样本本身不覆盖M31方向，所以不能把它当作M31方向的”测量值”。
Figure 3使用统一的单位：\(10^{-15}\) erg cm\(^{-2}\) s\(^{-1}\) arcmin\(^{-2}\)。
原始数据经历了两次单位变换：
最终： \[ S_{\mathrm{fig}} = S_{\mathrm{abs}} \times \frac{1000}{3600} = \frac{S_{\mathrm{abs}}}{3.6} \]
# Demonstrate the unit conversion chain
intr_1e12 = henley["surface_brightness_intrinsic_0p5_2p0_1e12"].iloc[0]
ratio = henley["absorbed_to_intrinsic_0p5_2p0_ratio"].iloc[0]
abs_fluxunit = henley["absorbed_0p5_2p0_fluxunit"].iloc[0]
# Manual reconstruction
abs_manual = intr_1e12 * ratio / 3.6 # /3.6 = *1000/3600
print(f"Sightline 1:")
print(f" Intrinsic: {intr_1e12:.4f} × 10^-12 erg s^-1 cm^-2 deg^-2")
print(f" Ratio: {ratio:.6f}")
print(f" Absorbed (manual): {abs_manual:.6f}")
print(f" Absorbed (CSV): {abs_fluxunit:.6f}")
np.testing.assert_allclose(abs_manual, abs_fluxunit, rtol=1e-10)
print(" ✓ Manual reconstruction matches CSV")Sightline 1:
Intrinsic: 3.5800 × 10^-12 erg s^-1 cm^-2 deg^-2
Ratio: 0.914851
Absorbed (manual): 0.909769
Absorbed (CSV): 0.909769
✓ Manual reconstruction matches CSV| 步骤 | 假设 | 来源 |
|---|---|---|
| 光谱模型 | Raymond-Smith等离子体 + phabs吸收 | 原文spectral fitting convention |
| 元素丰度 | angr (Anders & Grevesse 1989) | 原文使用 |
| 截面库 | bcmc (Balucinska-Church & McCammon 1992) | 原文使用 |
| \(N_\mathrm{H}\) | LAB 21 cm巡天 | 原文Table |
| 样本排除 | 排除Seq=103（无唯一LAB \(N_\mathrm{H}\)） | 本项目决定 |
| 统计量 | median + IQR，不是mean ± std | 本项目选择（抗离群） |
| M31关系 | 86条视线不指向M31；这是population context | 原文样本域 \(|b|>30^\circ\) |
选择第50条视线，用它的intrinsic brightness、温度和\(N_\mathrm{H}\)，手动验证absorbed值。提示：absorbed_to_intrinsic_0p5_2p0_ratio列已算好转换比。
Seq 68: manual=0.550924, CSV=0.550924
✓ Match计算86条视线的温度与absorbed-to-intrinsic ratio之间的Pearson相关系数。这个相关性是正还是负？为什么？
Pearson r = 0.0279
Positive — hotter plasmas have harder spectra, less affected by absorption.graph TD
A["VizieR J/ApJ/773/92<br/>110 sightlines"] --> B["筛选：有halo detection<br/>+ 有唯一 LAB N_H → 86条"]
B --> C["逐视线：Raymond-Smith(kT)<br/>× phabs(N_H) → T_0.5-2.0"]
C --> D["S_abs = S_intr × T"]
D --> E["deg⁻² → arcmin⁻² (/3600)<br/>10⁻¹² → 10⁻¹⁵ (×1000)"]
E --> F["86个值的 median + IQR"]
F --> G["Figure 3: H&S13<br/>population median = 0.385"]graph TD
A["VizieR J/ApJ/773/92<br/>110 sightlines"] --> B["筛选：有halo detection<br/>+ 有唯一 LAB N_H → 86条"]
B --> C["逐视线：Raymond-Smith(kT)<br/>× phabs(N_H) → T_0.5-2.0"]
C --> D["S_abs = S_intr × T"]
D --> E["deg⁻² → arcmin⁻² (/3600)<br/>10⁻¹² → 10⁻¹⁵ (×1000)"]
E --> F["86个值的 median + IQR"]
F --> G["Figure 3: H&S13<br/>population median = 0.385"]
一句话总结：Henley & Shelton (2013) 给出了86条高银纬视线上银河系热气体的intrinsic X射线亮度。经过视线-specific的吸收转换和单位换算后，取median得到0.385——这是银河系热气体在高银纬天空的典型亮度，不是M31方向的直接测量。
教程结束 🎓 下一步：继续阅读 Ponti+2023 eFEDS tutorial，了解远天区固定SWCX情景的band reconstruction。