---
title: "From Locatelli 2024 to Figure 3: The Complete X-Ray Measurement Conversion Pipeline for the Milky Way Circumgalactic Medium"
subtitle: "An Executable Tutorial for First-Year Physics Undergraduates"
author: "M31 CGM Research Group"
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
---
# From Locatelli 2024 to Figure 3: The Complete X-Ray Measurement Conversion Pipeline for the Milky Way Circumgalactic Medium
> **Tutorial Goal**: Walk you step by step through the conversion from a raw astronomy-paper measurement to the value used in our paper's Figure 3. The entire pipeline runs in Jupyter, with live code, visualizations, and detailed annotations.
> **Target Audience**: First-year physics undergraduates (have completed general physics and calculus; no astronomy background required)
> **Core Question**: Locatelli et al. measured the O VIII line intensity of the Milky Way halo with eROSITA, but what we need is the **0.5–2.0 keV broadband, Galactic-absorption-corrected surface brightness** in the **direction of M31 (Andromeda Galaxy)**. How many conversion steps separate the two? What physical assumption does each step encode?
---
## 0. Setup: Environment and Data
First, install the required Python packages (if running locally):
```bash
pip install numpy pandas matplotlib astropy scipy jupyter
```
On the cloud (e.g. Google Colab) you can run directly without installation.
```{python}
#| label: setup
#| echo: true
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pathlib import Path
# Set up plotting style
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 to load a CJK font if available
try:
plt.rcParams['font.sans-serif'] = ['DejaVu Sans', 'SimHei', 'Microsoft YaHei']
except:
pass
print("Environment ready!")
print(f"NumPy: {np.__version__}")
print(f"Pandas: {pd.__version__}")
```
---
## 1. Background: What Are We Measuring?
### 1.1 Scientific Goal: The Hot Gas Halo of M31
M31 (the Andromeda Galaxy) is our nearest large galactic neighbour. It should be surrounded by a hot gas halo (the circumgalactic medium, CGM) at a temperature of about $10^6$ K, which radiates predominantly soft X-rays (0.5–2.0 keV).
**But there is a big complication**: we observe M31 from inside the Milky Way, so every line of sight passes through our Galaxy's own hot gas halo (the MW CGM) and a cold absorbing layer. **The X-rays we see toward M31 = M31 signal + Galactic foreground + background.** To extract the M31 signal we must accurately model and subtract the Galactic foreground.
### 1.2 What Did Locatelli 2024 Do?
Locatelli et al. (2024, A&A) used **eROSITA** all-sky survey data from the first cycle (eRASS1) to measure the intensity distribution of the **O VIII emission line** (0.654 keV, 80 eV window) over the **western Galactic hemisphere** ($180^\circ < l < 360^\circ$).
> **Key point**: what they measured is **line intensity** (units: photons s$^{-1}$ cm$^{-2}$ sr$^{-1}$), not broadband continuum surface brightness; and it covers only the western hemisphere, while **M31 lies in the eastern hemisphere ($l \approx 121^\circ$)** — completely outside their fit domain!
### 1.3 Their Model: Spherical Halo + Exponential Disk
They fitted the O VIII intensity map with two geometric components:
1. **Spherical halo**: $n_h(r) = C \cdot r^{-3\beta}$, with $\beta=0.5$ (reference model)
2. **Exponential disk**: $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. **Emissivity**: $\epsilon \propto n_h^2 + n_d^2$ (**no cross term**; this is the author-designated reference model)
4. **Plasma parameters**: $kT = 0.15$ keV, $Z = 0.1 Z_\odot$
---
## 2. Load the Raw Data: Predictions for 14 M31 XMM Fields
Our project has already integrated Locatelli's reference geometry (Combined $\beta=0.5$) along the 14 M31 sightlines starting from the Solar position, then performed a chain of conversions. Let us first look at the final data table:
```{python}
#| label: load-data
#| echo: true
# Read the CSV (these are the 14-field predictions converted from the Locatelli model)
data_path = Path("m31_cgmsum_locatelli2024_reference_beta0p5_m31_footprint_predictions.csv")
df = pd.read_csv(data_path)
print(f"Data shape: {df.shape}")
print(f"Columns:\n{list(df.columns)}")
print()
print("First 5 rows of key columns:")
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. Step One: From 3D Density to Emission Measure (EM)
### 3.1 Physical Principle
The X-ray surface brightness is proportional to the **emission measure**:
$$EM = \int n_e n_H \, ds \approx \int n^2 \, ds$$
where $n$ is the electron density and $s$ is the line-of-sight distance. For a fully ionized plasma, $n_e \approx 1.2 n_H$.
Locatelli's model gives the densities $n_h(r)$ and $n_d(R,z)$; our code integrates along each M31 sightline starting from the Solar position ($R_\odot = 8.2$ kpc):
$$EM_i = \int_0^{s_\text{max}} [n_h^2(s) + n_d^2(s)] \, ds$$
### 3.2 Visualization: EM Distribution of the 14 Fields
```{python}
#| label: em-plot
#| echo: true
#| fig-cap: "Emission measure (EM) distribution of the 14 M31 XMM fields, coloured by North/South side."
#| fig-width: 8
#| fig-height: 4
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
# Left panel: EM vs Galactic latitude
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('Galactic latitude $b$ (deg)')
axes[0].set_ylabel('Emission measure EM (kpc cm$^{-6}$)')
axes[0].set_title('EM vs Galactic latitude')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# Right panel: EM vs Galactic longitude
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('Galactic longitude $l$ (deg)')
axes[1].set_ylabel('Emission measure EM (kpc cm$^{-6}$)')
axes[1].set_title('EM vs Galactic longitude')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"EM range: {df['reference_emission_measure_n2_kpc_cm-6'].min():.2e} - {df['reference_emission_measure_n2_kpc_cm-6'].max():.2e} kpc cm^-6")
```
---
## 4. Step Two: From EM to Intrinsic O VIII Line Intensity
### 4.1 Physical Principle: Line Emission from Thermal Plasma
For a CIE (collisional-ionization equilibrium) plasma at temperature $kT = 0.15$ keV and abundance $Z = 0.1 Z_\odot$, the **emissivity coefficient** of the O VIII line (0.654 keV) can be computed with the APEC model:
$$\text{O VIII intensity (L.U.)} = EM \times \Lambda_{\text{O VIII}}(kT, Z)$$
where L.U. = Line Unit = photons s$^{-1}$ cm$^{-2}$ sr$^{-1}$.
The Locatelli paper approximates eROSITA's energy response to O VIII with an 80 eV window. Our code convolves the APEC line with an **80 eV FWHM Gaussian** to mimic the same instrumental response.
### 4.2 Data Validation
```{python}
#| label: o8-intensity
#| echo: true
o8_col = 'reference_intrinsic_o8_0p614_0p694_intensity_lu'
print(f"Intrinsic O VIII intensity range: {df[o8_col].min():.3f} - {df[o8_col].max():.3f} L.U.")
print(f"Mean: {df[o8_col].mean():.3f} L.U.")
# Check the linear relationship between EM and O VIII
em = df['reference_emission_measure_n2_kpc_cm-6']
o8 = df[o8_col]
ratio = o8 / em
print(f"O VIII / EM ratio range: {ratio.min():.3f} - {ratio.max():.3f} L.U. / (kpc cm^-6)")
print(f"(should be a constant, since the same T/Z model)")
```
---
## 5. Step Three: O VIII Normalization Anchor
### 5.1 Why Anchor?
Locatelli's fit targets the **eROSITA O VIII map**. The model parameters ($C$, $n_0$, etc.) are pinned down by matching the **overall normalization** of that map. When we extrapolate the model toward M31, we **must keep this normalization fixed** — it is the only bridge that connects "the paper's measurement" to "our prediction".
### 5.2 Our Approach
1. Compute the model's intrinsic O VIII intensity at every field (Step 2)
2. Convolve the APEC line with an **80 eV FWHM Gaussian** to mimic the eROSITA energy response
3. This yields an **emissivity closure scale factor** (the `o8_response_matched_emissivity_closure_scale` column in the data table)
4. The factor is about **0.9499** and is identical for every field (same model, same response function)
```{python}
#| label: o8-anchor
#| echo: true
scale_col = 'o8_response_matched_emissivity_closure_scale'
print(f"O VIII response-matched scale factor: {df[scale_col].unique()}")
print(f"(identical for every field: same model + same instrumental response)")
# What this factor does: intrinsic O VIII × scale ≈ what eROSITA would have measured
df['o8_observed_matched'] = df[o8_col] * df[scale_col]
print(f"Matched O VIII range: {df['o8_observed_matched'].min():.3f} - {df['o8_observed_matched'].max():.3f} L.U.")
```
---
## 6. Step Four: The Line-to-Broadband Bridge (the Key Conversion!)
### 6.1 The Core Difficulty
- **The paper's measurement**: O VIII **line** intensity (0.614–0.694 keV, narrow window)
- **What Figure 3 needs**: **broadband** 0.5–2.0 keV surface brightness (after absorption)
Two physical conversions sit between them:
1. **Line → continuum**: the same plasma emits both lines and continuum, with a ratio that depends on temperature and abundance
2. **Abundance change**: the paper uses $Z=0.1$, while our Figure 3 uniformly uses **$Z=0.3$** (closer to typical halo-gas abundance)
### 6.2 Conversion Formula
$$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}$$
where $\Lambda(E)$ is the APEC emissivity spectrum. Note that the denominator uses $Z=0.1$ (the original model) and the numerator uses $Z=0.3$ (the target band).
### 6.3 Corresponding Columns in the Data Table
- `reference_response_matched_absorbed_0p5_2p0_figure_fluxunit`: response-matched only; **abundance unchanged**, **no absorption**
- `line_normalized_z0p3_absorbed_0p5_2p0_figure_fluxunit`: **abundance changed to Z=0.3, but no absorption**
```{python}
#| label: line-to-broadband
#| echo: true
#| fig-cap: "Effect of the line-to-broadband conversion: raising abundance from 0.1 to 0.3 substantially boosts the broadband flux."
#| 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. Before/after the conversion
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 reference line
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, no absorption (units: 10$^{-15}$ erg s$^{-1}$ cm$^{-2}$ arcmin$^{-2}$)')
axes[0].set_ylabel('Z=0.3, no absorption')
axes[0].set_title('Effect of the abundance change')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 2. Distribution of conversion factors
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'mean = {ratio_z.mean():.3f}')
axes[1].set_xlabel('Conversion factor (Z=0.3 / Z=0.1)')
axes[1].set_ylabel('Number of fields')
axes[1].set_title('Distribution of line-to-broadband factors')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
# 3. Relation to O VIII intensity
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 (no absorption)')
axes[2].set_title('O VIII intensity vs broadband flux')
plt.colorbar(axes[2].collections[0], ax=axes[2], label='Galactic latitude $b$ (deg)')
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"Z=0.1 no-absorption range: {df[col_no_z].min():.4f} - {df[col_no_z].max():.4f}")
print(f"Z=0.3 no-absorption range: {df[col_z03].min():.4f} - {df[col_z03].max():.4f}")
print(f"Mean conversion factor: {ratio_z.mean():.4f}")
```
---
## 7. Step Five: Galactic Absorption (HI4PI Full-Sky Column Density)
### 7.1 Physical Principle: Photoelectric Absorption
X-rays are absorbed as they pass through the cold neutral-hydrogen (HI) layer of the Milky Way; the cross section roughly scales as $E^{-3}$. The absorption factor is:
$$\text{Transmission} = \exp[-\sigma(E) \times N_H]$$
where $N_H$ is the hydrogen column density (cm$^{-2}$) and $\sigma(E)$ is the photoelectric absorption cross section (using the `tbabs` or `phabs` model).
We use the per-field $N_H$ from the **HI4PI all-sky survey** (the `nh_hi4pi_1e22_cm-2` column in the data table, units $10^{22}$ cm$^{-2}$).
### 7.2 Effect of Absorption
```{python}
#| label: absorption
#| echo: true
#| fig-cap: "Effect of the HI4PI absorption correction: high-column-density fields are attenuated more strongly."
#| 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. Before/after absorption
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, no absorption')
axes[0].set_ylabel('Z=0.3, with HI4PI absorption')
axes[0].set_title('Effect of the absorption correction')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 2. Transmission factor 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('Transmission (after / before absorption)')
axes[1].set_title('Absorption vs column density')
plt.colorbar(axes[1].collections[0], ax=axes[1], label='Galactic latitude $b$ (deg)')
axes[1].grid(True, alpha=0.3)
# 3. Distribution of final values
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('Final predicted value (10$^{-15}$ erg s$^{-1}$ cm$^{-2}$ arcmin$^{-2}$)')
axes[2].set_ylabel('Number of fields')
axes[2].set_title('Distribution of final predictions across 14 fields')
axes[2].legend()
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"Final prediction range: {df[col_abs].min():.6f} - {df[col_abs].max():.6f}")
print(f"Mean: {df[col_abs].mean():.6f}")
```
---
## 8. Step Six: Weighted Average to the All-field Value
### 8.1 Why a Weighted Average?
The 14 fields differ in exposure time, background level, and signal-to-noise. We use the **same inverse-variance weights as the measured CGMsum** so that the model-side prediction and the observation-side total are compared in an identical statistical frame.
```{python}
#| label: weighted-mean
#| echo: true
# Load the CGMsum weights (inverse-variance from the measured data)
cgmsum_path = Path("m31_cgmsum_v19_primary_measurements_public.csv")
cgmsum = pd.read_csv(cgmsum_path)
# Weights are stored in some column, or computed as inverse variance
# Demonstration: simple mean vs weighted mean
simple_mean = df[col_abs].mean()
# Assume weights correlate with exposure time (real application uses inverse variance)
weights = cgmsum['pha_exposure_s'] if 'pha_exposure_s' in cgmsum.columns else np.ones(len(df))
weights = weights[:len(df)] # Ensure length matches
weighted_mean = np.average(df[col_abs], weights=weights)
print(f"Simple mean: {simple_mean:.6f}")
print(f"Weighted mean (using exposure time as weights for demo): {weighted_mean:.6f}")
print(f"All-field value in the registry: 1.290382")
print()
print("Note: the actual weights are the CGMsum inverse-variance weights;")
print("exposure time is used here only for demonstration, so the result is close but not identical.")
```
---
## 9. Final Comparison: Model vs Observation
### 9.1 Summary of Key Numbers
| Quantity | Value | Meaning |
|------|------|------|
| Locatelli all-field prediction | **1.290382** | Model-side prediction, including all conversion assumptions |
| 14-field footprint span | **1.157655 – 1.341463** | Deterministic spatial variation, not a confidence interval |
| Measured CGMsum total | **0.964727** | Observation under the same weights |
| Observed − Model | **−0.325655** | **Negative value = model tension, not negative M31 emission** |
### 9.2 Visualization: Where It Sits in the Figure-3 Context
```{python}
#| label: figure3-context
#| echo: true
#| fig-cap: "Where the Locatelli prediction sits in the Figure-3 context: it exceeds the measured total, forming a tension boundary."
#| fig-width: 9
#| fig-height: 5
fig, ax = plt.subplots(figsize=(9, 5))
# Model prediction distribution
ax.hist(df[col_abs], bins=12, alpha=0.5, label='14-field predictions',
color='orange', edgecolor='k', density=True)
# Key reference lines
ax.axvline(1.290382, color='red', linewidth=3, label='All-field prediction (1.290382)')
ax.axvline(0.964727, color='blue', linewidth=3, label='Measured CGMsum (0.964727)')
ax.axvline(1.157655, color='red', linestyle='--', alpha=0.7, label='Footprint range')
ax.axvline(1.341463, color='red', linestyle='--', alpha=0.7)
# Comparison with Ueda and HaloSat (from the registry)
ax.axvline(0.311904, color='green', linewidth=2, linestyle=':', label='Ueda prediction (~0.312)')
ax.axvline(0.565386, color='cyan', linewidth=2, linestyle=':', label='HaloSat extrapolation (~0.565)')
ax.set_xlabel('Absorbed 0.5–2.0 keV surface brightness (10$^{-15}$ erg s$^{-1}$ cm$^{-2}$ arcmin$^{-2}$)')
ax.set_ylabel('Density (normalized)')
ax.set_title('Locatelli prediction vs measured total vs other MW models')
ax.legend(loc='upper right', fontsize=10)
ax.grid(True, alpha=0.3)
# Annotate the 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. Key Assumption Checklist (Must Remember!)
Every conversion step encodes an assumption. **The final value 1.290382 is not "the Locatelli-measured M31 foreground"**, but:
| Step | Assumption | Source |
|------|------|------|
| 1. Geometry extrapolation | Combined $\beta=0.5$ model extrapolated from the western hemisphere to the M31 direction | **Project augmentation**; the paper does not publish an M31 prediction |
| 2. Emissivity implementation | $n_h^2 + n_d^2$, **no cross term** | Paper's Eq. (7) omits the path-length $s$; project follows the upstream formula |
| 3. Temperature/abundance | $kT=0.15$ keV, $Z=0.1 \to 0.3$ | Paper uses $Z=0.1$; Figure 3 uniformly uses $Z=0.3$ |
| 4. Line-to-broadband | APEC CIE model, 80 eV Gaussian response matching | **Project augmentation**; the paper only publishes the O VIII line |
| 5. Absorption model | HI4PI full-screen $N_H$, `phabs`/`tbabs` | **Conservative assumption**; the paper treats absorption more carefully |
| 6. Weights | Same inverse-variance weights as CGMsum | **Project augmentation**, for comparability |
> **Safe conclusion**: Locatelli provides strong support for **local disk-like O VIII emission**, but this **named-reference conversion cannot serve as an error-free foreground prior toward M31**. Changing the geometry, line normalization, abundance, or absorber placement changes the tension; with no published parameter covariance, **the displayed range is not a posterior interval**.
---
## 11. Hands-On Exercises: Compute It Yourself!
### Exercise 1: Verify the Line-to-Broadband Factor
```{python}
#| label: exercise1
#| echo: true
# Hints: use the astropy APEC model or read a precomputed table
# Simplified: assume you know the integrated emissivity at kT=0.15 keV for Z=0.1 and Z=0.3
# Pseudocode:
# from astropy.modeling import models
# or use a precomputed XSPEC table
#
# ratio = flux(0.5-2.0 keV, Z=0.3) / flux(O VIII line, Z=0.1)
# should approximately equal the observed ~1.6-1.8 factor
print("Exercise: try using XSPEC or SPEX to compute at kT=0.15 keV")
print("the APEC emissivity ratio: broadband(0.5-2.0, Z=0.3) / O VIII line(0.614-0.694, Z=0.1)")
print("Expected: roughly ~1.7")
```
### Exercise 2: Compute the Absorption Factor
```{python}
#| label: exercise2
#| echo: true
# Given N_H = 5.8e20 cm^-2 (a typical field), compute the transmission at 0.5 and 2.0 keV
# Using the phabs model: sigma(E) ~ E^-3
#
# transmission = exp(-N_H * sigma(E))
#
# Hint: 0.5 keV is strongly absorbed, 2.0 keV weakly
# This explains why absorption changes the effective band shape
import numpy as np
# Simplified: Morrison & McCammon cross-section approximation
# sigma(E) ≈ 3e-22 * (E/1keV)^-3 cm^2 (rough)
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: transmission = {t:.3f} ({100*(1-t):.1f}% absorbed)")
print("\nNotice: low energies are strongly absorbed, reshaping the spectrum — that is why absorption correction must be done energy-by-energy.")
```
### Exercise 3: Change Assumptions and Watch the Tension Move
```{python}
#| label: exercise3
#| echo: true
# Suppose we change abundance from Z=0.3 to Z=0.5, or disk scale height from 1.1 kpc to 2.0 kpc
# Qualitative predictions:
#
# 1. Z up -> broadband flux up (more metal lines) -> tension up
# 2. Thicker disk -> more disk emission toward M31 -> tension up
# 3. beta from 0.5 to 0.3 -> flatter halo -> halo contribution toward M31 changes
# 4. Partial-screen instead of full-screen absorption -> less absorption -> tension up
#
print("Qualitative prediction exercise:")
print("1. Z=0.3 -> Z=0.5: broadband flux ______ (up/down), tension ______")
print("2. Disk scale height 1.1 -> 2.0 kpc: M31 disk emission ______, tension ______")
print("3. Full-screen -> partial-screen absorption: transmission ______, tension ______")
print()
print("Answers:")
print("1. Up, larger (metal lines boost the broadband)")
print("2. Up, larger (more disk gas along the line of sight)")
print("3. Up, larger (less absorption, more model flux observed)")
```
---
## 12. Summary: The Complete Paper-to-Figure-3 Chain
```{mermaid}
graph TD
A["Locatelli 2024 raw data<br/>eRASS1 O VIII 0.614-0.694 keV map<br/>western hemisphere 180<l<360°"] --> B["Fit geometric model<br/>spherical halo β=0.5 + exponential disk<br/>n_h²+n_d², no cross term"]
B --> C["Reference Combined β=0.5 params<br/>C=0.046, n₀=0.032, Rh=6.2, zh=1.1<br/>kT=0.15 keV, Z=0.1"]
C --> D["Step 1: integrate along 14 M31 sightlines from the Solar position<br/>yielding EM and intrinsic O VIII (4.14-4.56 L.U.)"]
D --> E["Step 2: 80 eV Gaussian matching the eROSITA response<br/>emissivity closure scale = 0.9499"]
E --> F["Step 3: Line-to-Broadband bridge<br/>fix O VIII normalization, Z=0.1 → Z=0.3 APEC<br/>convert to 0.5-2.0 keV"]
F --> G["Step 4: HI4PI full-screen absorption<br/>per-field N_H via phabs/tbabs"]
G --> H["Step 5: CGMsum same-weight average<br/>All-field = 1.290382"]
H --> I["Figure 3: orange vertical line<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. Appendix: Complete Reproducible Code
If you want a complete local reproduction (needs CIAO/Sherpa/XSPEC):
```python
# Full reproduction script (pseudocode)
# The actual project contains this in generate_reference_presentation_figures.py and similar scripts
import numpy as np
from astropy.modeling import models
# or use sherpa/astro_models
# 1. Define the density model
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. Integrate along the line of sight (needs coordinate transforms: celestial -> Galactic -> 3D)
# This part uses a dedicated geometry library in the project code
# 3. APEC emissivity calculation
# Use XSPEC: apec(kT=0.15, Z=0.1) for the O VIII line emissivity
# then apec(kT=0.15, Z=0.3) for the 0.5-2.0 keV broadband emissivity
# 4. Absorption
# Use XSPEC: phabs(NH) * apec(...)
# 5. Weighted average
# weights = 1 / variance_of_measured_CGMsum
# final = sum(w * model) / sum(w)
```
---
## References
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 (upstream source of the cross-term path-length $s$)
3. **HI4PI Collaboration 2016**, *HI4PI: A full-sky H I survey*, A&A 594, A116
4. **M31 CGM paper** (this project): `paper_apj_v19_cgmsum_draft/`
---
> **End of tutorial** 🎓
> If you can run every code cell, understand the physics of each step, and explain why the final tension is $-0.326$ rather than "a negative M31 halo", congratulations — you have mastered one of the most central skills in X-ray astrophysics: **foreground modelling and conversion**!
>
> Suggested next step: read the analogous conversions for Ueda 2022 and HaloSat 2020, and compare how different geometric models predict different values toward M31.
---
*This tutorial is generated from the M31 CGM research group's reproducible analysis pipeline. All data, code, and intermediate products are available in the project repository.*
