BOSA: log sSFR sweep at fixed log L_TIR

Sweep log10(sSFR / yr^-1) across the canonical 14-point BOSA grid (Boquien & Salim 2021) at fixed log10 L_TIR = 11 (typical LIRG luminosity). Higher sSFR → harder mid-IR colour and stronger PAH features; quiescent (low-sSFR) galaxies have a colder FIR peak.

Reference

Boquien M. & Salim S. 2021, A&A 653 A149, arXiv:2106.04595. Library: https://salims.pages.iu.edu/bosa/.

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/bosa_templates.h5"),
        Path("../data/bosa_templates.h5"),
        Path("../../data/bosa_templates.h5"),
    ):
        if p.exists():
            return str(p)
    return None


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

with h5py.File(_PATH, "r") as f:
    wave_aa = np.asarray(f["wavelength_aa"][:])
    log_ltir = np.asarray(f["log_ltir_grid"][:])
    log_ssfr = np.asarray(f["log_ssfr_grid"][:])
    spectra = np.asarray(f["spectra"][:])  # (n_ltir, n_ssfr, n_wave) normalised L_nu

wave_um = wave_aa * 1.0e-4
i_ltir = int(np.argmin(np.abs(log_ltir - 11.0)))

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(3.0, 1.0e3)
# BOSA templates are pre-normalised to ∫L_ν dν = 1 (CGS), so
# nu*L_nu peaks dimensionlessly at ~0.5.  Y-axis range covers the
# 3-4 decades of meaningful spectral structure (PAH features at
# 6-13 μm and the FIR peak at ~100 μm).
ax.set_ylim(1.0e-3, 2.0e0)

cmap = plt.get_cmap("viridis")
for k, lssfr in enumerate(log_ssfr):
    L_nu = spectra[i_ltir, k]
    ax.plot(
        wave_um, nu * L_nu,
        color=cmap(k / max(1, len(log_ssfr) - 1)),
        lw=1.3,
        label=rf"$\log_{{10}} \mathrm{{sSFR}} = {lssfr:+.1f}$",
    )

ax.legend(loc="lower left", frameon=False, fontsize=8, ncol=3)
ax.set_title(
    rf"BOSA (Boquien & Salim 2021) at $\log_{{10}} L_{{\rm TIR}}={log_ltir[i_ltir]:.1f}$",
    fontsize=11,
)

plt.show()