from __future__ import annotations from pathlib import Path import matplotlib.pyplot as plt from matplotlib.patches import Ellipse import numpy as np import pandas as pd import scienceplots # noqa: F401 from astropy.coordinates import SkyCoord import astropy.units as u from scipy.integrate import trapezoid HERE = Path(__file__).resolve().parent ROOT = HERE.parent OUTPUT_DIR = HERE PRIMARY_CSV = ( ROOT / "independent_cgm_results_v19_mos14to20keV_excluded_hi4pi_nh_per_pointing" / "independent_cgm_v19_all22_measurements.csv" ) PRIOR_LEDGER_CSV = HERE / "m31_cgmsum_conditional_prior_ledger.csv" PREDICTIONS_CSV = HERE / "m31_cgmsum_total_difference_predictions.csv" FLUX_UNIT = 1.0e-15 NORTH_FLUX = 1.1062308841694078 NORTH_ERROR = 0.05357907323160643 SOUTH_FLUX = 0.893850973242182 SOUTH_ERROR = 0.2032861967092256 ALL_FLUX = 1.0923775268463574 ALL_ERROR = 0.05186950369566604 # For an APEC plasma with kT=0.175306 keV, Z=0.3 Z_sun, phabs nH=0.056, # angr abundances, and vern cross-sections, recomputed with Sherpa 4.18/XSPEC: # F_abs(0.4-1.25 keV) / F_intrinsic(0.5-2.0 keV). PRIMARY_TO_INTRINSIC_05_2_RATIO = 0.7719131901134512 KPC_CM = 3.0856775814913673e21 ARCMIN2_SR = (np.pi / (180.0 * 60.0)) ** 2 LUMINOSITY_SURFACE_DENSITY_TO_PRIMARY = ( ARCMIN2_SR / (4.0 * np.pi * KPC_CM**2) / FLUX_UNIT * PRIMARY_TO_INTRINSIC_05_2_RATIO ) def expected_output_paths() -> list[Path]: stems = ( "m31_cgmsum_mw_m31_constraint_plane", "m31_cgmsum_total_difference_predictive_plane", ) return [ Path(OUTPUT_DIR) / f"{stem}_{column}.{suffix}" for stem in stems for column in ("1col", "2col") for suffix in ("png", "pdf") ] def measurement_total_difference() -> tuple[float, float, np.ndarray]: total = 0.5 * (NORTH_FLUX + SOUTH_FLUX) difference = NORTH_FLUX - SOUTH_FLUX covariance = np.array( [ [0.25 * (NORTH_ERROR**2 + SOUTH_ERROR**2), 0.5 * (NORTH_ERROR**2 - SOUTH_ERROR**2)], [0.5 * (NORTH_ERROR**2 - SOUTH_ERROR**2), NORTH_ERROR**2 + SOUTH_ERROR**2], ] ) return total, difference, covariance def _primary_fields() -> pd.DataFrame: table = pd.read_csv(PRIMARY_CSV, dtype={"obsid": str}) return table.loc[table["quality_primary"]].copy() def _weighted_side_mean(table: pd.DataFrame, values: np.ndarray, side: str) -> float: error_name = "absorbed_flux_soft_0p4_1p25_staterr_erg_cm-2_s-1_arcmin-2" mask = table["side"].eq(side).to_numpy() weights = 1.0 / table.loc[mask, error_name].to_numpy() ** 2 return float(np.sum(weights * values[mask]) / np.sum(weights)) def _contrast_slope(table: pd.DataFrame, values: np.ndarray) -> float: north = _weighted_side_mean(table, values, "North/NW") south = _weighted_side_mean(table, values, "South/SE") return 2.0 * (north - south) / (north + south) def _zhang_profile(radius_kpc: np.ndarray) -> np.ndarray: core_radius_kpc = 7.39 beta = 0.37 return (1.0 + (radius_kpc / core_radius_kpc) ** 2) ** (0.5 - 3.0 * beta) def _locatelli_emission_shape( l_rad: np.ndarray, b_rad: np.ndarray, *, beta: float, c_norm: float, n0: float, radial_scale_kpc: float, height_scale_kpc: float, ) -> np.ndarray: distance_kpc = np.linspace(0.001, 350.0, 35001) solar_radius_kpc = 8.2 output = [] for longitude, latitude in zip(l_rad, b_rad, strict=True): cylindrical_radius = np.sqrt( solar_radius_kpc**2 + (distance_kpc * np.cos(latitude)) ** 2 - 2.0 * solar_radius_kpc * distance_kpc * np.cos(latitude) * np.cos(longitude) ) height = distance_kpc * np.sin(latitude) spherical_radius = np.sqrt(cylindrical_radius**2 + height**2) spherical_density = ( c_norm * spherical_radius ** (-3.0 * beta) if c_norm > 0.0 else np.zeros_like(distance_kpc) ) disk_density = n0 * np.exp(-cylindrical_radius / radial_scale_kpc) * np.exp( -np.abs(height) / height_scale_kpc ) output.append( trapezoid(spherical_density**2 + disk_density**2, distance_kpc) ) return np.asarray(output) def _apec_primary_band_transmission(nh_1e22: np.ndarray, kT_keV: float) -> np.ndarray: nh_grid = np.linspace(0.04, 0.08, 9) transmission_grids = { 0.15: np.array([0.729846506, 0.701996621, 0.675279178, 0.649645958, 0.625050780, 0.601449548, 0.578800060, 0.557062121, 0.536197223]), 0.178: np.array([0.744631128, 0.718073718, 0.692547309, 0.668009357, 0.644419095, 0.621737548, 0.599927377, 0.578952992, 0.558780250]), 0.225: np.array([0.772413462, 0.748317246, 0.725066336, 0.702628298, 0.680972008, 0.660067629, 0.639886507, 0.620401284, 0.601585631]), } return np.interp(nh_1e22, nh_grid, transmission_grids[kT_keV]) def model_contrast_slopes() -> dict[str, float | tuple[float, float]]: table = _primary_fields() m31_slope = _contrast_slope(table, _zhang_profile(table["rproj_kpc"].to_numpy())) # Locatelli et al. (2024), Table 1. The variants bracket the spatial # systematics used in that paper. Only relative projected morphology is # needed here. The absorbed-band prediction also includes the HI4PI nH # variation across the footprint through APEC transmission lookup tables. variants = ( ("disk", 0.5, 0.0, 0.017, 8.0, 3.3, 0.15), ("beta03", 0.3, 0.0009, 0.017, 8.0, 3.3, 0.15), ("beta05", 0.5, 0.046, 0.032, 6.2, 1.1, 0.15), ("beta07", 0.7, 0.176, 0.016, 9.9, 1.6, 0.15), ("swcx", 0.5, 0.032, 0.063, 3.9, 0.9, 0.178), ("instr", 0.5, 0.043, 0.055, 4.6, 0.9, 0.15), ("high_emissivity", 0.5, 0.0136, 0.0094, 6.2, 1.1, 0.225), ) galactic = SkyCoord( table["ra_deg"].to_numpy() * u.deg, table["dec_deg"].to_numpy() * u.deg, ).galactic l_rad = galactic.l.rad b_rad = galactic.b.rad nh_1e22 = table["nh_hi4pi_1e22_cm-2"].to_numpy() mw_slopes: dict[str, float] = {} for name, beta, c_norm, n0, radial_scale, height_scale, kT_keV in variants: values = _locatelli_emission_shape( l_rad, b_rad, beta=beta, c_norm=c_norm, n0=n0, radial_scale_kpc=radial_scale, height_scale_kpc=height_scale, ) values *= _apec_primary_band_transmission(nh_1e22, kT_keV) mw_slopes[name] = _contrast_slope(table, values) family = (min(mw_slopes.values()), max(mw_slopes.values())) return { "Zhang 2024 symmetric M31": m31_slope, "Locatelli 2024 MW central": mw_slopes["swcx"], "Locatelli 2024 MW family": family, "Uniform component": 0.0, } def build_prior_ledger() -> pd.DataFrame: rows: list[dict[str, object]] = [] hs_conversion = PRIMARY_TO_INTRINSIC_05_2_RATIO / 3600.0 * 1.0e3 rows.append( { "prior_id": "henley_shelton2013_high_latitude_population", "component": "MW", "label": "Henley & Shelton 2013 high-|b| population", "reference": "Henley & Shelton (2013)", "doi": "10.1088/0004-637X/773/2/92", "original_quantity": "intrinsic 0.5-2 keV halo surface brightness", "original_central": 1.5, "original_low": 0.5, "original_high": 7.0, "original_units": "1e-12 erg cm-2 s-1 deg-2", "conversion_to_primary": hs_conversion, "primary_central": 1.5 * hs_conversion, "primary_low": 0.5 * hs_conversion, "primary_high": 7.0 * hs_conversion, "scope": "population context; not a local M31-direction measurement", "caveat": "M31 lies outside the paper's |b|>30 degree domain; APEC conversion is conditional.", "provenance": "published range and median; converted with the v19 CGMsum APEC template", } ) rows.append( { "prior_id": "locatelli2024_m31_footprint_model_family", "component": "MW", "label": "Locatelli et al. 2024 O VIII morphology", "reference": "Locatelli et al. (2024)", "doi": "10.1051/0004-6361/202347061", "original_quantity": "O VIII-fitted 3D density-model family projected to the primary footprint", "original_central": np.nan, "original_low": np.nan, "original_high": np.nan, "original_units": "model parameters in their Table 1", "conversion_to_primary": np.nan, "primary_central": 1.435214486512026, "primary_low": 0.7278438371779999, "primary_high": 1.6419010865000614, "scope": "M31-direction model interpolation over the 14-field primary estimator", "caveat": "Not a broadband local measurement; depends on O VIII model geometry, thermochemistry, and all-foreground-screen conversion.", "provenance": "seven published spatial/thermochemical variants projected field by field with HI4PI nH", } ) zhang_luminosity_density = 1.7252986024966184e36 zhang_primary = zhang_luminosity_density * LUMINOSITY_SURFACE_DENSITY_TO_PRIMARY rows.append( { "prior_id": "zhang2024_m31_mass_nominal_20kpc", "component": "M31", "label": "Zhang et al. 2024 M31-mass stack", "reference": "Zhang et al. (2024)", "doi": "10.1051/0004-6361/202449412", "original_quantity": "nominal beta-profile luminosity surface density at 20 kpc", "original_central": zhang_luminosity_density, "original_low": zhang_luminosity_density, "original_high": zhang_luminosity_density, "original_units": "erg s-1 kpc-2, intrinsic 0.5-2 keV", "conversion_to_primary": LUMINOSITY_SURFACE_DENSITY_TO_PRIMARY, "primary_central": zhang_primary, "primary_low": zhang_primary, "primary_high": zhang_primary, "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.", "provenance": "published beta profile evaluated at 20 kpc and converted with the v19 APEC template", } ) grayson = ( ("eagle_agndt9", "EAGLE-AGNdT9", 4.704476594777882e35, 4.184141227254975e35, 5.1535212647748536e35), ("eagle", "EAGLE", 9.38177618116973e35, 8.236144961922412e35, 1.041199465882404e36), ("eagle_noagn", "EAGLE-NoAGN", 1.9709847221830294e36, 1.6858151590706473e36, 2.1874197384055262e36), ("simba", "SIMBA", 3.405962162605925e36, 2.838273867820679e36, 3.982108803986934e36), ("simba_noagn", "SIMBA-NoAGN", 2.2218741538234737e37, 2.1647474653904115e37, 2.340689955138248e37), ) for identifier, label, central, low, high in grayson: rows.append( { "prior_id": f"grayson2025_{identifier}", "component": "M31", "label": f"Grayson et al. 2025 {label}", "reference": "Grayson et al. (2025)", "doi": "10.3847/1538-4357/ae100f", "original_quantity": "digitized 10-30 kpc bin in the M31-mass stellar bin", "original_central": central, "original_low": low, "original_high": high, "original_units": "erg s-1 kpc-2, intrinsic 0.5-2 keV", "conversion_to_primary": LUMINOSITY_SURFACE_DENSITY_TO_PRIMARY, "primary_central": central * LUMINOSITY_SURFACE_DENSITY_TO_PRIMARY, "primary_low": low * LUMINOSITY_SURFACE_DENSITY_TO_PRIMARY, "primary_high": high * LUMINOSITY_SURFACE_DENSITY_TO_PRIMARY, "scope": "simulation prediction matched to the 11 None: for suffix in ("png", "pdf"): path = Path(OUTPUT_DIR) / f"{stem}_{column}.{suffix}" fig.savefig(path, dpi=400 if suffix == "png" else None, bbox_inches="tight") def _constraint_plane(priors: pd.DataFrame, column: str) -> None: width = 3.35 if column == "1col" else 7.0 height = 4.6 if column == "1col" else 5.0 fig, ax = plt.subplots(figsize=(width, height)) x = np.linspace(0.0, 2.05, 1000) constraints = ( (NORTH_FLUX, NORTH_ERROR, "#2166ac", "N/NW"), (SOUTH_FLUX, SOUTH_ERROR, "#b2182b", "S/SE"), ) for total, error, color, label in constraints: lower = total - error - x upper = total + error - x ax.fill_between(x, np.maximum(lower, 0.0), np.maximum(upper, 0.0), color=color, alpha=0.18) valid = total - x >= 0.0 ax.plot(x[valid], total - x[valid], color=color, lw=1.5, label=rf"{label}: ${total:.3f}\pm{error:.3f}$") valid = ALL_FLUX - x >= 0.0 ax.plot(x[valid], ALL_FLUX - x[valid], color="black", lw=1.1, ls="--", label=rf"All: ${ALL_FLUX:.3f}\pm{ALL_ERROR:.3f}$") mw = priors.loc[priors["component"].eq("MW")] hs = mw.loc[mw["prior_id"].str.startswith("henley")].iloc[0] loc = mw.loc[mw["prior_id"].str.startswith("locatelli")].iloc[0] ax.axvspan(hs.primary_low, hs.primary_high, color="0.45", alpha=0.055, zorder=0) ax.axvline(hs.primary_central, color="0.38", lw=1.0, ls=":") ax.axvspan(loc.primary_low, loc.primary_high, facecolor="#e69f00", alpha=0.075, hatch="///", edgecolor="#e69f00", lw=0.0, zorder=0) ax.axvline(loc.primary_central, color="#e69f00", lw=1.1, ls="-.") ybar = 0.045 ax.plot([hs.primary_low, hs.primary_high], [ybar, ybar], color="0.35", lw=2.3, solid_capstyle="butt") ax.plot(hs.primary_central, ybar, marker="|", color="0.2", ms=8) ax.text(hs.primary_low, ybar + 0.035, "HS13 population", fontsize=6.0 if column == "1col" else 7.2, color="0.25") ax.plot([loc.primary_low, loc.primary_high], [2.0 * ybar, 2.0 * ybar], color="#d17c00", lw=2.3, solid_capstyle="butt") ax.plot(loc.primary_central, 2.0 * ybar, marker="|", color="#a45f00", ms=8) ax.text(loc.primary_low, 2.0 * ybar + 0.035, "Locatelli+24 model family", fontsize=6.0 if column == "1col" else 7.2, color="#9a5a00") colors = { "grayson2025_eagle_agndt9": "#008080", "grayson2025_eagle": "#000080", "grayson2025_eagle_noagn": "#56b4e9", "grayson2025_simba": "#8b0000", "zhang2024_m31_mass_nominal_20kpc": "black", } label_text = { "grayson2025_eagle_agndt9": "EAGLE-AGNdT9", "grayson2025_eagle": "EAGLE", "grayson2025_eagle_noagn": "EAGLE-NoAGN", "grayson2025_simba": "SIMBA", "zhang2024_m31_mass_nominal_20kpc": "Zhang+24 stack", } for _, row in priors.loc[priors["component"].eq("M31")].iterrows(): identifier = str(row.prior_id) if identifier == "grayson2025_simba_noagn": continue color = colors[identifier] low = float(row.primary_low) high = float(row.primary_high) central = float(row.primary_central) if high > low: ax.axhspan(low, high, color=color, alpha=0.075, zorder=0) linestyle = "--" if identifier.startswith("zhang") else "-" ax.axhline(central, color=color, lw=1.0, ls=linestyle) ax.text( 2.02, central, f"{label_text[identifier]} {central:.2f}", color=color, fontsize=5.8 if column == "1col" else 7.0, va="center", ha="right", bbox={"facecolor": "white", "edgecolor": "none", "alpha": 0.72, "pad": 0.5}, ) if 0.0 <= ALL_FLUX - central <= 2.05: ax.plot( ALL_FLUX - central, central, marker="o", ms=3.3, markerfacecolor="white", markeredgecolor=color, markeredgewidth=0.9, zorder=6, ) noagn = priors.loc[priors["prior_id"].eq("grayson2025_simba_noagn")].iloc[0] ax.annotate( f"SIMBA-NoAGN: {noagn.primary_central:.1f} (off scale)", xy=(1.78, 2.095), xytext=(1.10, 2.095), arrowprops={"arrowstyle": "-|>", "color": "#fa8072", "lw": 0.9}, color="#c55449", fontsize=6.0 if column == "1col" else 7.2, va="center", ) ax.set_xlim(0.0, 2.05) ax.set_ylim(0.0, 2.2) ax.set_xlabel(r"MW CGM contribution, $x$") ax.set_ylabel(r"M31 inner-halo contribution, $y$") ax.text(0.02, 0.98, "LG set to zero; non-negative components", transform=ax.transAxes, ha="left", va="top", fontsize=5.8 if column == "1col" else 7.2) ax.text(0.02, 0.945, "Absorbed 0.4--1.25 keV; units $10^{-15}$ erg cm$^{-2}$ s$^{-1}$ arcmin$^{-2}$", transform=ax.transAxes, ha="left", va="top", fontsize=4.8 if column == "1col" else 6.2) ax.legend(loc="upper left", bbox_to_anchor=(0.02, 0.90), fontsize=5.3 if column == "1col" else 6.8, frameon=True) ax.tick_params(direction="in", top=True, right=True) _save_figure(fig, "m31_cgmsum_mw_m31_constraint_plane", column) plt.close(fig) def _covariance_ellipse(mean: np.ndarray, covariance: np.ndarray, delta_chi2: float, **kwargs: object) -> Ellipse: eigenvalues, eigenvectors = np.linalg.eigh(covariance) order = np.argsort(eigenvalues)[::-1] eigenvalues = eigenvalues[order] eigenvectors = eigenvectors[:, order] angle = np.degrees(np.arctan2(eigenvectors[1, 0], eigenvectors[0, 0])) scale = np.sqrt(delta_chi2) return Ellipse( xy=mean, width=2.0 * scale * np.sqrt(eigenvalues[0]), height=2.0 * scale * np.sqrt(eigenvalues[1]), angle=angle, **kwargs, ) def _residual_significance(slope: float, total: float, difference: float, covariance: np.ndarray) -> float: residual = difference - slope * total variance = covariance[1, 1] + slope**2 * covariance[0, 0] - 2.0 * slope * covariance[0, 1] return residual / np.sqrt(variance) def _predictive_plane(column: str) -> pd.DataFrame: width = 3.35 if column == "1col" else 7.0 height = 3.35 if column == "1col" else 4.5 fig, ax = plt.subplots(figsize=(width, height)) total, difference, covariance = measurement_total_difference() slopes = model_contrast_slopes() m31_slope = float(slopes["Zhang 2024 symmetric M31"]) mw_slope = float(slopes["Locatelli 2024 MW central"]) mw_low, mw_high = slopes["Locatelli 2024 MW family"] mw_low = float(mw_low) mw_high = float(mw_high) xx = np.linspace(0.0, 1.65, 500) ax.fill_between(xx, m31_slope * xx, mw_high * xx, color="#999999", alpha=0.12, label="Non-negative MW/M31 mixtures") ax.plot(xx, m31_slope * xx, color="#009e73", lw=1.5, label=rf"Zhang+24 symmetric M31 ($D/T={m31_slope:+.3f}$)") ax.fill_between(xx, mw_low * xx, mw_high * xx, color="#e69f00", alpha=0.18) ax.plot(xx, mw_slope * xx, color="#d17c00", lw=1.5, label=rf"Locatelli+24 MW ($D/T={mw_slope:+.3f}$)") ax.axhline(0.0, color="0.35", lw=1.0, ls="--", label="Uniform component") mean = np.array([total, difference]) ax.add_patch(_covariance_ellipse(mean, covariance, 5.991, facecolor="#2166ac", edgecolor="none", alpha=0.10)) ax.add_patch(_covariance_ellipse(mean, covariance, 2.279, facecolor="#2166ac", edgecolor="#2166ac", lw=1.0, alpha=0.22)) ax.errorbar(total, difference, fmt="o", color="#08306b", ms=4.2, zorder=5, label="Observed North-South estimator") model_rows = [] for name, slope, color in ( ("Uniform", 0.0, "0.25"), ("Zhang 2024 symmetric M31", m31_slope, "#009e73"), ("Locatelli 2024 MW central", mw_slope, "#d17c00"), ): predicted = slope * total significance = _residual_significance(slope, total, difference, covariance) ax.plot(total, predicted, marker="s", ms=4.0, color=color, zorder=4) model_rows.append( { "model": name, "contrast_slope_D_over_T": slope, "observed_total": total, "predicted_difference_at_observed_total": predicted, "observed_difference": difference, "residual_sigma_with_TD_covariance": significance, } ) ax.text( 0.03, 0.97, "Absorbed 0.4--1.25 keV; units $10^{-15}$ erg cm$^{-2}$ s$^{-1}$ arcmin$^{-2}$\nObserved: " + rf"$T={total:.3f}$, $D={difference:.3f}$; " + rf"M31-only offset={model_rows[1]['residual_sigma_with_TD_covariance']:.2f}$\sigma$", transform=ax.transAxes, ha="left", va="top", fontsize=5.5 if column == "1col" else 7.2, ) ax.text(0.98, 0.96, "68/95% covariance", transform=ax.transAxes, ha="right", va="top", color="#2166ac", fontsize=5.1 if column == "1col" else 6.8) ax.set_xlim(0.0, 1.65) ax.set_ylim(-0.48, 0.70) ax.set_xlabel(r"Side-symmetrized total, $T=(F_N+F_S)/2$") ax.set_ylabel(r"North--South difference, $D=F_N-F_S$") handles, _ = ax.get_legend_handles_labels() labels = [ "Non-negative mixtures", rf"Symmetric M31 ($D/T={m31_slope:+.3f}$)", rf"MW model ($D/T={mw_slope:+.3f}$)", "Uniform component", ] ax.legend(handles[:4], labels, loc="lower right", fontsize=5.2 if column == "1col" else 6.8, frameon=True) ax.tick_params(direction="in", top=True, right=True) _save_figure(fig, "m31_cgmsum_total_difference_predictive_plane", column) plt.close(fig) return pd.DataFrame(model_rows) def main() -> None: plt.style.use(["science", "no-latex"]) priors = build_prior_ledger() priors.to_csv(PRIOR_LEDGER_CSV, index=False) for column in ("1col", "2col"): _constraint_plane(priors, column) predictions = _predictive_plane("1col") _predictive_plane("2col") predictions.to_csv(PREDICTIONS_CSV, index=False) print(f"Wrote {PRIOR_LEDGER_CSV}") print(f"Wrote {PREDICTIONS_CSV}") for path in expected_output_paths(): print(path) if __name__ == "__main__": main()