---
title: "从 Locatelli 2024 到 Figure 3：银河系晕气体的 X 射线测量转换全过程"
subtitle: "面向物理系一年级本科生的可执行教程"
author: "M31 CGM 研究组"
date: "2026-07-17"
format:
html:
theme:
light: [cosmo, custom.scss]
dark: [cosmo, custom-dark.scss]
toc: true
toc-depth: 3
code-fold: true
code-tools: true
fig-cap-location: bottom
fig-width: 7
fig-height: 4
fig-align: center
embed-resources: true
css: styles.css
execute:
echo: true
eval: true
warning: false
message: false
jupyter: python3
---
# 从 Locatelli 2024 到 Figure 3：银河系晕气体的 X 射线测量转换全过程
> **教程目标**：带你一步步理解如何把一篇天文学论文的原始测量结果，转换成我们论文 Figure 3 中使用的数值。全程可在 Jupyter 中运行，包含实时代码、可视化图表和详细注释。
> **目标读者**：物理系一年级本科生（已修普通物理、微积分，无天文背景也可跟上）
> **核心问题**：Locatelli 等人用 eROSITA 测到了银河系晕的 O VIII 谱线强度，但我们要的是 **0.5–2.0 keV 宽波段、经过银河系吸收后的表面亮度**，且是在 **M31（仙女座星系）方向**上。这中间隔了多少步？每一步都在做什么物理假设？
---
## 0. 准备工作：环境与数据
首先安装必要的 Python 包（如果在本地运行）：
```bash
pip install numpy pandas matplotlib astropy scipy jupyter
```
在云端（如 Google Colab）直接运行即可，无需安装。
```{python}
#| label: setup
#| echo: true
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
plt.rcParams['axes.unicode_minus'] = False
# 尝试加载中文字体（如果可用）
try:
plt.rcParams['font.sans-serif'] = ['DejaVu Sans', 'SimHei', 'Microsoft YaHei']
except:
pass
print("环境就绪！")
print(f"NumPy: {np.__version__}")
print(f"Pandas: {pd.__version__}")
```
---
## 1. 背景知识：我们在测什么？
### 1.1 科学目标：M31 的热气体晕
M31（仙女座星系）是我们最近的大型邻居星系。它周围应该有一个热气体晕（CGM, Circumgalactic Medium），温度约 $10^6$ K，主要辐射软 X 射线（0.5–2.0 keV）。
**但有一个大麻烦**：我们在银河系内部观测 M31，视线穿过银河系自己的热气体晕（MW CGM）和冷气体吸收层。**我们在 M31 方向看到的 X 射线 = M31 信号 + 银河系前景 + 背景**。要提取 M31 信号，必须精确建模并减去银河系前景。
### 1.2 Locatelli 2024 做了什么？
Locatelli 等人 (2024, A&A) 使用 **eROSITA** 全天巡天第一轮数据（eRASS1），测量了 **O VIII 谱线**（0.654 keV，80 eV 窗口）在 **银河系西半球**（$180^\circ < l < 360^\circ$）的强度分布。
> **关键点**：他们测的是 **谱线强度**（单位：photons s$^{-1}$ cm$^{-2}$ sr$^{-1}$），不是宽波段连续谱表面亮度；且只覆盖西半球，**M31 在东半球（$l \approx 121^\circ$）**，完全在他们的拟合域之外！
### 1.3 他们的模型：球形晕 + 指数盘
他们用两个几何组件拟合 O VIII 强度图：
1. **球形晕**：$n_h(r) = C \cdot r^{-3\beta}$，$\beta=0.5$（参考模型）
2. **指数盘**：$n_d(R,z) = n_0 \exp(-R/R_h) \exp(-|z|/z_h)$
- $n_0 = 0.032$ cm$^{-3}$, $R_h = 6.2$ kpc, $z_h = 1.1$ kpc
3. **发射率**：$\epsilon \propto n_h^2 + n_d^2$（**无交叉项**，这是作者指定的 reference 模型）
4. **等离子体参数**：$kT = 0.15$ keV, $Z = 0.1 Z_\odot$
---
## 2. 加载原始数据：14 个 M31 XMM 视场的预测值
我们的项目已经把 Locatelli 的 reference 几何模型（Combined $\beta=0.5$）从太阳位置沿 14 条 M31 视线积分，并做了一系列转换。让我们先看看最终数据表：
```{python}
#| label: load-data
#| echo: true
# 读取 CSV 数据（这是从 Locatelli 模型转换而来的 14 个 M31 视场预测）
data_path = Path("m31_cgmsum_locatelli2024_reference_beta0p5_m31_footprint_predictions.csv")
df = pd.read_csv(data_path)
print(f"数据形状: {df.shape}")
print(f"列名:\n{list(df.columns)}")
print()
print("前 5 行关键列:")
key_cols = ['obsid', 'side', 'galactic_l_deg', 'galactic_b_deg',
'reference_emission_measure_n2_kpc_cm-6',
'reference_intrinsic_o8_0p614_0p694_intensity_lu',
'line_normalized_z0p3_absorbed_0p5_2p0_figure_fluxunit']
print(df[key_cols].head().to_string())
```
---
## 3. 第一步：从 3D 密度到发射量（Emission Measure, EM）
### 3.1 物理原理
X 射线表面亮度与 **发射量** 成正比：
$$EM = \int n_e n_H \, ds \approx \int n^2 \, ds$$
其中 $n$ 是电子密度，$s$ 是视线距离。对于全电离等离子体，$n_e \approx 1.2 n_H$。
Locatelli 的模型给出密度 $n_h(r)$ 和 $n_d(R,z)$，我们的代码从太阳位置（$R_\odot = 8.2$ kpc）出发，沿每条 M31 视线积分：
$$EM_i = \int_0^{s_\text{max}} [n_h^2(s) + n_d^2(s)] \, ds$$
### 3.2 可视化：14 个视场的 EM 分布
```{python}
#| label: em-plot
#| echo: true
#| fig-cap: "14 个 M31 XMM 视场的发射量（EM）分布。颜色按 North/South 分组。"
#| fig-width: 8
#| fig-height: 4
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
# 左图：EM vs 银纬
colors = {'North/NW': '#2196F3', 'South/SE': '#F44336'}
for side in df['side'].unique():
subset = df[df['side'] == side]
axes[0].scatter(subset['galactic_b_deg'],
subset['reference_emission_measure_n2_kpc_cm-6'],
label=side, color=colors[side], s=80, alpha=0.8, edgecolor='k')
axes[0].set_xlabel('银纬 $b$ (deg)')
axes[0].set_ylabel('发射量 EM (kpc cm$^{-6}$)')
axes[0].set_title('EM vs 银纬')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 右图：EM vs 银经
for side in df['side'].unique():
subset = df[df['side'] == side]
axes[1].scatter(subset['galactic_l_deg'],
subset['reference_emission_measure_n2_kpc_cm-6'],
label=side, color=colors[side], s=80, alpha=0.8, edgecolor='k')
axes[1].set_xlabel('银经 $l$ (deg)')
axes[1].set_ylabel('发射量 EM (kpc cm$^{-6}$)')
axes[1].set_title('EM vs 银经')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"EM 范围: {df['reference_emission_measure_n2_kpc_cm-6'].min():.2e} - {df['reference_emission_measure_n2_kpc_cm-6'].max():.2e} kpc cm^-6")
```
---
## 4. 第二步：从 EM 到 O VIII 内在谱线强度
### 4.1 物理原理：热等离子体的谱线辐射
对于温度 $kT = 0.15$ keV、丰度 $Z = 0.1 Z_\odot$ 的 CIE（碰撞电离平衡）等离子体，O VIII 谱线（0.654 keV）的**发射率系数** 可以用 APEC 模型计算：
$$\text{O VIII intensity (L.U.)} = EM \times \Lambda_{\text{O VIII}}(kT, Z)$$
其中 L.U. = Line Unit = photons s$^{-1}$ cm$^{-2}$ sr$^{-1}$。
Locatelli 论文用 80 eV 窗口近似 eROSITA 对 O VIII 的能量响应。我们的代码用 **80 eV FWHM 高斯函数** 卷积 APEC 谱线，模拟相同的仪器响应。
### 4.2 数据验证
```{python}
#| label: o8-intensity
#| echo: true
o8_col = 'reference_intrinsic_o8_0p614_0p694_intensity_lu'
print(f"O VIII 内在强度范围: {df[o8_col].min():.3f} - {df[o8_col].max():.3f} L.U.")
print(f"均值: {df[o8_col].mean():.3f} L.U.")
# 检查 EM 与 O VIII 的线性关系
em = df['reference_emission_measure_n2_kpc_cm-6']
o8 = df[o8_col]
ratio = o8 / em
print(f"O VIII / EM 比值范围: {ratio.min():.3f} - {ratio.max():.3f} L.U. / (kpc cm^-6)")
print(f"(应为常数，因为同一温度/丰度模型)")
```
---
## 5. 第三步：O VIII 归一化锚定
### 5.1 为什么要锚定？
Locatelli 的拟合是针对 **eROSITA 观测到的 O VIII 图** 的。模型参数（$C, n_0$ 等）是通过匹配观测图的**整体归一化**确定的。当我们把模型外推到 M31 方向时，**必须保持这个归一化不变**——这是连接“原论文测量”和“我们的预测”的唯一桥梁。
### 5.2 我们的做法
1. 计算模型在每个视场的 intrinsic O VIII 强度（步骤 2）
2. 用 **80 eV FWHM 高斯** 卷积 APEC 谱线，模拟 eROSITA 的能量响应
3. 这给出一个 **emissivity closure scale factor**（数据表中的 `o8_response_matched_emissivity_closure_scale`）
4. 该因子约为 **0.9499**，对所有视场相同（因为同一模型、同一响应函数）
```{python}
#| label: o8-anchor
#| echo: true
scale_col = 'o8_response_matched_emissivity_closure_scale'
print(f"O VIII 响应匹配尺度因子: {df[scale_col].unique()}")
print(f"(所有视场相同，因为同一模型+同一仪器响应)")
# 这个因子的作用：intrinsic O VIII × scale = 模拟 eROSITA 测得的 O VIII
df['o8_observed_matched'] = df[o8_col] * df[scale_col]
print(f"匹配后 O VIII 范围: {df['o8_observed_matched'].min():.3f} - {df['o8_observed_matched'].max():.3f} L.U.")
```
---
## 6. 第四步：Line-to-Broadband 桥梁（关键转换！）
### 6.1 核心难题
- **原论文测量**：O VIII **谱线**强度（0.614–0.694 keV，窄窗）
- **Figure 3 需要**：**宽波段** 0.5–2.0 keV 表面亮度（吸收后）
这中间隔了两个物理转换：
1. **谱线 → 连续谱**：同一等离子体既发谱线也发连续谱，比例取决于温度和丰度
2. **丰度变化**：原论文用 $Z=0.1$，我们的 Figure 3 统一用 **$Z=0.3$**（更接近典型晕气体丰度）
### 6.2 转换公式
$$F_{\text{0.5-2.0 keV}} = F_{\text{O VIII}} \times \frac{\int_{0.5}^{2.0} \Lambda(E; kT, Z=0.3) dE}{\int_{\text{O VIII line}} \Lambda(E; kT, Z=0.1) dE}$$
其中 $\Lambda(E)$ 是 APEC 发射率谱。注意分母用 $Z=0.1$（原模型），分子用 $Z=0.3$（目标波段）。
### 6.3 数据表中的对应列
- `reference_response_matched_absorbed_0p5_2p0_figure_fluxunit`：只做了响应匹配，**未改丰度**，**未吸收**
- `line_normalized_z0p3_absorbed_0p5_2p0_figure_fluxunit`：**改了丰度到 Z=0.3，但未吸收**
```{python}
#| label: line-to-broadband
#| echo: true
#| fig-cap: "Line-to-broadband 转换的效果：丰度从 0.1 到 0.3 显著增加宽波段通量。"
#| fig-width: 10
#| fig-height: 4
col_no_z = 'reference_response_matched_absorbed_0p5_2p0_figure_fluxunit'
col_z03 = 'line_normalized_z0p3_absorbed_0p5_2p0_figure_fluxunit'
fig, axes = plt.subplots(1, 3, figsize=(14, 4))
# 1. 转换前后对比
for side in df['side'].unique():
subset = df[df['side'] == side]
axes[0].scatter(subset[col_no_z], subset[col_z03],
label=side, color=colors[side], s=80, alpha=0.8, edgecolor='k')
# 1:1 线
mx = max(df[col_no_z].max(), df[col_z03].max())
mn = min(df[col_no_z].min(), df[col_z03].min())
axes[0].plot([mn, mx], [mn, mx], 'k--', alpha=0.5, label='1:1')
axes[0].set_xlabel('Z=0.1, 无吸收 (单位: 10$^{-15}$ erg s$^{-1}$ cm$^{-2}$ arcmin$^{-2}$)')
axes[0].set_ylabel('Z=0.3, 无吸收')
axes[0].set_title('丰度变化的影响')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 2. 转换因子分布
ratio_z = df[col_z03] / df[col_no_z]
axes[1].hist(ratio_z, bins=10, edgecolor='k', alpha=0.7, color='steelblue')
axes[1].axvline(ratio_z.mean(), color='red', linestyle='--',
label=f'均值 = {ratio_z.mean():.3f}')
axes[1].set_xlabel('转换因子 (Z=0.3 / Z=0.1)')
axes[1].set_ylabel('视场数')
axes[1].set_title('Line-to-broadband 因子分布')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
# 3. 与 O VIII 的关系
axes[2].scatter(df[o8_col], df[col_z03], c=df['galactic_b_deg'],
cmap='RdBu_r', s=80, edgecolor='k')
axes[2].set_xlabel('Intrinsic O VIII (L.U.)')
axes[2].set_ylabel('Z=0.3 Broadband (无吸收)')
axes[2].set_title('O VIII 强度 vs 宽波段通量')
plt.colorbar(axes[2].collections[0], ax=axes[2], label='银纬 $b$ (deg)')
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"Z=0.1 无吸收范围: {df[col_no_z].min():.4f} - {df[col_no_z].max():.4f}")
print(f"Z=0.3 无吸收范围: {df[col_z03].min():.4f} - {df[col_z03].max():.4f}")
print(f"平均转换因子: {ratio_z.mean():.4f}")
```
---
## 7. 第五步：银河系吸收（HI4PI 全天柱密度）
### 7.1 物理原理：光电吸收
X 射线穿过银河系冷中性氢（HI）层时会被吸收，截面随能量大致按 $E^{-3}$ 下降。吸收因子：
$$\text{Transmission} = \exp[-\sigma(E) \times N_H]$$
其中 $N_H$ 是氢原子柱密度（cm$^{-2}$），$\sigma(E)$ 是光电吸收截面（用 `tbabs` 或 `phabs` 模型）。
我们用 **HI4PI 全天调查** 提供的每个视场的 $N_H$（数据表中的 `nh_hi4pi_1e22_cm-2`，单位 $10^{22}$ cm$^{-2}$）。
### 7.2 吸收的效果
```{python}
#| label: absorption
#| echo: true
#| fig-cap: "HI4PI 吸收校正的效果：高柱密度视场衰减更强。"
#| fig-width: 10
#| fig-height: 4
col_abs = 'line_normalized_z0p3_absorbed_0p5_2p0_figure_fluxunit'
fig, axes = plt.subplots(1, 3, figsize=(14, 4))
# 1. 吸收前后对比
for side in df['side'].unique():
subset = df[df['side'] == side]
axes[0].scatter(subset[col_z03], subset[col_abs],
label=side, color=colors[side], s=80, alpha=0.8, edgecolor='k')
mx = max(df[col_z03].max(), df[col_abs].max())
mn = min(df[col_z03].min(), df[col_abs].min())
axes[0].plot([mn, mx], [mn, mx], 'k--', alpha=0.5)
axes[0].set_xlabel('Z=0.3, 无吸收')
axes[0].set_ylabel('Z=0.3, 含 HI4PI 吸收')
axes[0].set_title('吸收校正效果')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 2. 吸收因子 vs N_H
transmission = df[col_abs] / df[col_z03]
axes[1].scatter(df['nh_hi4pi_1e22_cm-2'], transmission,
c=df['galactic_b_deg'], cmap='RdBu_r', s=80, edgecolor='k')
axes[1].set_xlabel('$N_H$ ($10^{22}$ cm$^{-2}$)')
axes[1].set_ylabel('透射率 (吸收后 / 吸收前)')
axes[1].set_title('吸收随柱密度变化')
plt.colorbar(axes[1].collections[0], ax=axes[1], label='银纬 $b$ (deg)')
axes[1].grid(True, alpha=0.3)
# 3. 最终值分布
for side in df['side'].unique():
subset = df[df['side'] == side]
axes[2].hist(subset[col_abs], bins=8, alpha=0.6, label=side,
color=colors[side], edgecolor='k')
axes[2].set_xlabel('最终预测值 (10$^{-15}$ erg s$^{-1}$ cm$^{-2}$ arcmin$^{-2}$)')
axes[2].set_ylabel('视场数')
axes[2].set_title('14 个视场最终预测分布')
axes[2].legend()
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"最终预测范围: {df[col_abs].min():.6f} - {df[col_abs].max():.6f}")
print(f"均值: {df[col_abs].mean():.6f}")
```
---
## 8. 第六步：加权平均得到 All-field 值
### 8.1 为什么要加权平均？
14 个视场的曝光时间、背景水平、信噪比不同。我们用 **与实测 CGMsum 相同的逆方差权重** 来加权平均，确保模型侧预测与观测侧在同一统计框架下可比。
```{python}
#| label: weighted-mean
#| echo: true
# 读取 CGMsum 权重（来自实测数据的 inverse-variance weights）
cgmsum_path = Path("m31_cgmsum_v19_primary_measurements_public.csv")
cgmsum = pd.read_csv(cgmsum_path)
# 权重通常存储在某列，或者我们计算 inverse variance
# 这里演示：简单平均 vs 加权平均
simple_mean = df[col_abs].mean()
# 假设权重与曝光时间相关（实际应用中用 inverse variance）
weights = cgmsum['pha_exposure_s'] if 'pha_exposure_s' in cgmsum.columns else np.ones(len(df))
weights = weights[:len(df)] # 确保长度匹配
weighted_mean = np.average(df[col_abs], weights=weights)
print(f"简单平均: {simple_mean:.6f}")
print(f"加权平均 (用曝光时间作权重演示): {weighted_mean:.6f}")
print(f"注册表中的 all-field 值: 1.290382")
print()
print("注意：实际权重是 CGMsum 的 inverse-variance weights，")
print("这里用曝光时间仅作演示，结果接近但不完全相同。")
```
---
## 9. 最终结果对比：模型 vs 观测
### 9.1 关键数字汇总
| 量值 | 数值 | 含义 |
|------|------|------|
| Locatelli all-field 预测 | **1.290382** | 模型侧预测，含所有转换假设 |
| 14 场 footprint 范围 | **1.157655 – 1.341463** | 确定性空间变化，非置信区间 |
| 实测 CGMsum 总量 | **0.964727** | 同一权重下的观测值 |
| Observed − Model | **−0.325655** | **负值 = 模型张力，非负 M31 发射** |
### 9.2 可视化：Figure 3 语境下的位置
```{python}
#| label: figure3-context
#| echo: true
#| fig-cap: "Locatelli 预测在 Figure 3 语境下的位置：超出实测总量，形成 tension boundary。"
#| fig-width: 9
#| fig-height: 5
fig, ax = plt.subplots(figsize=(9, 5))
# 模型预测分布
ax.hist(df[col_abs], bins=12, alpha=0.5, label='14 个视场预测',
color='orange', edgecolor='k', density=True)
# 关键标记
ax.axvline(1.290382, color='red', linewidth=3, label='All-field 预测 (1.290382)')
ax.axvline(0.964727, color='blue', linewidth=3, label='实测 CGMsum (0.964727)')
ax.axvline(1.157655, color='red', linestyle='--', alpha=0.7, label='Footprint 范围')
ax.axvline(1.341463, color='red', linestyle='--', alpha=0.7)
# Ueda 和 HaloSat 对比（来自注册表）
ax.axvline(0.311904, color='green', linewidth=2, linestyle=':', label='Ueda 预测 (~0.312)')
ax.axvline(0.565386, color='cyan', linewidth=2, linestyle=':', label='HaloSat 外推 (~0.565)')
ax.set_xlabel('吸收后 0.5–2.0 keV 表面亮度 (10$^{-15}$ erg s$^{-1}$ cm$^{-2}$ arcmin$^{-2}$)')
ax.set_ylabel('密度 (归一化)')
ax.set_title('Locatelli 预测 vs 实测总量 vs 其他 MW 模型')
ax.legend(loc='upper right', fontsize=10)
ax.grid(True, alpha=0.3)
# 标注 tension
ax.annotate('Tension: Model − Obs = +0.326',
xy=(1.13, 1.5), xytext=(1.13, 2.5),
arrowprops=dict(arrowstyle='->', color='red'),
fontsize=12, color='red', fontweight='bold')
plt.tight_layout()
plt.show()
```
---
## 10. 关键假设清单（必须记住！）
每一步转换都引入了假设，**最终数值 1.290382 不是 "Locatelli 测得的 M31 前景"**，而是：
| 步骤 | 假设 | 来源 |
|------|------|------|
| 1. 几何外推 | Combined $\beta=0.5$ 模型从西半球外推到 M31 方向 | **项目增强**，原论文无 M31 预测 |
| 2. 发射实现 | $n_h^2 + n_d^2$，**无交叉项** | 原论文 Eq. (7) 漏印 path-length $s$，项目按上游公式修正 |
| 3. 温度/丰度 | $kT=0.15$ keV, $Z=0.1 \to 0.3$ | 原论文 $Z=0.1$，Figure 3 统一用 $Z=0.3$ |
| 4. Line-to-broadband | APEC CIE 模型，80 eV Gaussian 响应匹配 | **项目增强**，原论文只给 O VIII 谱线 |
| 5. 吸收模型 | HI4PI 全屏 $N_H$，`phabs`/`tbabs` | **保守假设**，原论文用更复杂的吸收处理 |
| 6. 权重 | 与 CGMsum 相同的 inverse-variance weights | **项目增强**，确保可比性 |
> **安全结论**：Locatelli 强有力支持**本地盘状 O VIII 发射**，但这个 **named reference 转换不能作为 M31 方向无误差的前景先验**。改变几何、谱线归一化、丰度或吸收位置，tension 大小都会变；没有发表的参数协方差，**现有范围不是后验区间**。
---
## 11. 动手练习：自己算一遍！
### 练习 1：验证 Line-to-Broadband 因子
```{python}
#| label: exercise1
#| echo: true
# 提示：用 astropy 的 APEC 模型或读取预计算表
# 这里给出简化版：假设你已知 kT=0.15 keV, Z=0.1 和 Z=0.3 的发射率积分
# 伪代码：
# from astropy.modeling import models
# 或使用预计算的 XSPEC 表
#
# ratio = flux(0.5-2.0 keV, Z=0.3) / flux(O VIII line, Z=0.1)
# 应该约等于观测到的 ~1.6-1.8 倍
print("练习：尝试用 XSPEC 或 SPEX 计算 kT=0.15 keV 下的")
print("APEC 发射率积分比：宽波段(0.5-2.0, Z=0.3) / O VIII 线(0.614-0.694, Z=0.1)")
print("预期结果：~1.7 左右")
```
### 练习 2：吸收因子计算
```{python}
#| label: exercise2
#| echo: true
# 给定 N_H = 5.8e20 cm^-2 (典型视场)，计算 0.5 keV 和 2.0 keV 的透射率
# 用 phabs 模型：sigma(E) ~ E^-3
#
# transmission = exp(-N_H * sigma(E))
#
# 提示：0.5 keV 吸收强，2.0 keV 吸收弱
# 这解释了为什么吸收会改变有效波段形状
import numpy as np
# 简化计算：用 Morrison & McCammon 截面近似
# sigma(E) ≈ 3e-22 * (E/1keV)^-3 cm^2 (粗略)
E = np.array([0.5, 1.0, 2.0]) # keV
sigma = 3e-22 * (E / 1.0)**-3 # cm^2
NH = 5.8e20 # cm^-2
trans = np.exp(-NH * sigma)
for e, t in zip(E, trans):
print(f"E = {e} keV: 透射率 = {t:.3f} ({100*(1-t):.1f}% 被吸收)")
print("\n可见：低能端吸收强，改变谱形，这是吸收校正必须逐能段做的原因。")
```
### 练习 3：改变假设看 tension 如何变
```{python}
#| label: exercise3
#| echo: true
# 假设我们把丰度从 Z=0.3 改为 Z=0.5，或者把盘 scale height 从 1.1 kpc 改为 2.0 kpc
# 定性预测：
#
# 1. Z 增加 → 宽波段通量增加（金属线增多）→ tension 变大
# 2. 盘更厚 → M31 方向盘发射更多 → tension 变大
# 3. beta 从 0.5 改为 0.3 → 晕更扁平 → M31 方向晕贡献变化
# 4. 用 partial-screen 吸收而非 full-screen → 吸收减少 → tension 变大
#
print("定性预测练习：")
print("1. Z=0.3 -> Z=0.5: 宽波段通量 ______ (增加/减少)，tension ______")
print("2. 盘 scale height 1.1 -> 2.0 kpc: M31 方向盘发射 ______，tension ______")
print("3. Full-screen -> Partial-screen 吸收: 透射率 ______，tension ______")
print()
print("答案：")
print("1. 增加，变大（金属线增强宽波段）")
print("2. 增加，变大（更多盘气体在视线内）")
print("3. 增加，变大（吸收变少，观测到更多模型通量）")
```
---
## 12. 总结：从论文到 Figure 3 的完整链路
```{mermaid}
graph TD
A["Locatelli 2024 原始数据<br/>eRASS1 O VIII 0.614-0.694 keV 图<br/>西半球 180<l<360°"] --> B["拟合几何模型<br/>球形晕 β=0.5 + 指数盘<br/>n_h²+n_d² 无交叉项"]
B --> C["Reference Combined β=0.5 参数<br/>C=0.046, n₀=0.032, Rh=6.2, zh=1.1<br/>kT=0.15 keV, Z=0.1"]
C --> D["步骤 1: 从太阳位置积分 14 条 M31 视线<br/>得到 EM 和 intrinsic O VIII (4.14-4.56 L.U.)"]
D --> E["步骤 2: 80 eV Gaussian 匹配 eROSITA 响应<br/>得到 emissivity closure scale = 0.9499"]
E --> F["步骤 3: Line-to-Broadband 桥梁<br/>固定 O VIII 归一化，Z=0.1 → Z=0.3 APEC<br/>转换到 0.5-2.0 keV"]
F --> G["步骤 4: HI4PI 全屏吸收<br/>每视场 N_H 做 phabs/tbabs"]
G --> H["步骤 5: CGMsum 同权重加权平均<br/>All-field = 1.290382"]
H --> I["Figure 3: 橙色竖线<br/>Footprint 1.158-1.341<br/>Tension boundary at -0.326"]
style A fill:#e3f2fd
style I fill:#ffebee
style H fill:#fff3e0
```
---
## 13. 附录：完整可复现代码
如果你想在本地完整复现（需要 CIAO/Sherpa/XSPEC 环境）：
```python
# 完整复现脚本伪代码
# 实际项目中在 generate_reference_presentation_figures.py 等脚本中
import numpy as np
from astropy.modeling import models
# 或用 sherpa/astro_models
# 1. 定义密度模型
def n_halo(r, C=0.046, beta=0.5):
return C * r**(-3*beta)
def n_disk(R, z, n0=0.032, Rh=6.2, zh=1.1):
return n0 * np.exp(-R/Rh) * np.exp(-np.abs(z)/zh)
# 2. 沿视线积分 (需要坐标变换：天球 -> 银河 -> 3D)
# 这部分在项目代码中用专门的几何库实现
# 3. APEC 发射率计算
# 用 XSPEC: apec(kT=0.15, Z=0.1) 计算 O VIII 线发射率
# 再用 apec(kT=0.15, Z=0.3) 计算 0.5-2.0 keV 宽波段发射率
# 4. 吸收
# 用 XSPEC: phabs(NH) * apec(...)
# 5. 权重平均
# weights = 1 / variance_of_measured_CGMsum
# final = sum(w * model) / sum(w)
```
---
## 参考资料
1. **Locatelli et al. 2024**, *The warm-hot circumgalactic medium of the Milky Way as seen by eROSITA*, A&A 681, A78
2. **Miller & Bregman 2013**, *The Milky Way's Hot Gas Halo*, ApJ 770, 118 (交叉项 path-length s 的上游来源)
3. **HI4PI Collaboration 2016**, *HI4PI: A full-sky H I survey*, A&A 594, A116
4. **M31 CGM 论文** (本项目): `paper_apj_v19_cgmsum_draft/`
---
> **教程结束** 🎓
> 如果你能跟着跑通所有代码、理解每一步的物理意义，并能解释为什么最终 tension 是 −0.326 而不是“负的 M31 晕”，恭喜你——你已经掌握了 X 射线天体物理中最核心的**前景建模与转换**技能！
>
> 下一步建议：阅读 Ueda 2022 和 HaloSat 2020 的类似转换过程，对比不同几何模型在 M31 方向的预测差异。
---
*本教程基于 M31 CGM 研究组的可复现分析管线生成。所有数据、代码和中间产物均可在项目仓库中找到。*
