THEMIS: U_min sweep at fixed q_HAC

Sweep U_min over the full published 37-point CIGALE grid (0.10 to 80.0) at the fiducial q_HAC = 0.17 and alpha = 2. Higher U warms the dust → FIR peak shifts blueward and the MIR small-grain emission grows relative to the FIR cold peak.

Reference

Jones, A.P., Köhler, M., Ysard, N., et al. 2017, A&A 602, A46 (arXiv:1703.00775). Templates via CIGALE.

from pathlib import Path

import h5py
import jax
import matplotlib.pyplot as plt
import numpy as np

jax.config.update("jax_enable_x64", True)


def _find_h5():
    for p in (
        Path("data/themis_templates.h5"),
        Path("../data/themis_templates.h5"),
        Path("../../data/themis_templates.h5"),
    ):
        if p.exists():
            return str(p)
    return None


_PATH = _find_h5()
if _PATH is None:
    raise FileNotFoundError(
        "THEMIS HDF5 not found. Build with "
        "`python scripts/build_themis_hdf5.py --clone`."
    )

with h5py.File(_PATH, "r") as f:
    wave_aa = np.asarray(f["wavelength_aa"][:])
    qhac_grid = np.asarray(f["qhac_grid"][:])
    umin_grid = np.asarray(f["umin_grid"][:])
    single_u = np.asarray(f["single_u"][:])

wave_um = wave_aa * 1.0e-4
i_qhac = int(np.argmin(np.abs(qhac_grid - 0.17)))
c_aa_per_s = 2.99792458e18
nu = c_aa_per_s / wave_aa

fig, ax = plt.subplots(figsize=(8.0, 5.5), constrained_layout=True)
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel(r"$\lambda\ [\mu\mathrm{m}]$", fontsize=12)
ax.set_ylabel(
    r"$\nu L_\nu\ [\mathrm{normalised}\ \int L_\nu d\nu = 1]$",
    fontsize=11,
)
ax.set_xlim(2.0, 1.0e3)
ax.set_ylim(1.0e-26, 1.0e-22)

cmap = plt.get_cmap("plasma")
# Subsample 12 representative points across the full 37-point grid.
idx_show = np.linspace(0, len(umin_grid) - 1, 12).astype(int)
for k, iu in enumerate(idx_show):
    L_nu = single_u[i_qhac, iu]
    ax.plot(
        wave_um, nu * L_nu,
        color=cmap(k / max(1, len(idx_show) - 1)),
        lw=1.3,
        label=rf"$U_{{\rm min}}={umin_grid[iu]:.2f}$",
    )

ax.legend(loc="upper left", frameon=False, fontsize=8, ncol=2)
ax.set_title(
    rf"THEMIS (Jones+2017) at $q_{{\rm HAC}}={qhac_grid[i_qhac]:.2f}$, $\alpha=2$",
    fontsize=11,
)

plt.show()