Dust emission

Synthesizer can model thermal dust emission using an energy balance approach. This is done using dust emission models, which can then be handed to an EmissionModel to produce a dust emission spectrum, which can be combined with any other spectra in the model.

First we load some common modules, and define a wavelength array over which we will generate our dust emission SED

import matplotlib.pyplot as plt
import numpy as np
from unyt import Angstrom, K, Lsun, Msun, um

from synthesizer.emission_models.dust.emission import (
    Blackbody,
    Casey12,
    Greybody,
    IR_templates,
)

lam = 10 ** (np.arange(3.0, 8.0, 0.01)) * Angstrom

Parametric models

Below we show a number of simple models, including a blackbody, greybody (non-unity emissivity) and the model detailed in Casey (2012) (mid-IR power-law + far-IR greybody). By default a dust emission model provides a normalised spectrum.

for T in [10, 25, 50, 100, 200]:
    model = Blackbody(temperature=T * K)
    sed = model.get_spectra(lam)
    plt.loglog(sed.lam, sed.luminosity, label=f"T={T} K")

plt.title("Blackbody")

plt.xlabel(r"$ \lambda / \AA$")
plt.ylabel(r"$ L_{\nu} / L_{bol}$")

plt.xlim(1e4, 1e8)
plt.ylim(10**-3.0, 1.5)

plt.legend()
plt.show()

/opt/hostedtoolcache/Python/3.10.17/x64/lib/python3.10/site-packages/unyt/array.py:1832: RuntimeWarning: overflow encountered in exp
  out_arr = func(np.asarray(inp), out=out_func, **kwargs)
/opt/hostedtoolcache/Python/3.10.17/x64/lib/python3.10/site-packages/unyt/array.py:1972: RuntimeWarning: overflow encountered in multiply
  out_arr = func(

We show here the Greybody spectrum for the optically-thin (solid lines) and optically-thick (dashed lines) case for different dust temperatures.

for ii, T in enumerate([25, 50, 100, 200]):
    color = plt.cm.tab10(ii)
    # Optically thin
    model = Greybody(temperature=T * K, emissivity=2.0, optically_thin=True)
    sed = model.get_spectra(lam)
    plt.loglog(sed.lam, sed.luminosity, label=f"T={T} K", color=color)

    # Optically thick below 100um
    model = Greybody(
        temperature=T * K, emissivity=2.0, lam_0=100 * um, optically_thin=False
    )
    sed = model.get_spectra(lam)
    plt.loglog(sed.lam, sed.luminosity, ls="dashed", color=color)

plt.title("Greybody")

plt.xlabel(r"$ \lambda / \AA$")
plt.ylabel(r"$ L_{\nu} / L_{bol}$")

plt.xlim(1e4, 1e8)
plt.ylim(10**-3.0, 1.5)

plt.legend()
plt.show()

for ii, T in enumerate([25, 50, 100, 200]):
    color = plt.cm.tab10(ii)
    model = Casey12(
        temperature=T * K, emissivity=2.0, alpha=2.0, lam_0=100 * um
    )
    sed = model.get_spectra(lam)
    plt.loglog(sed.lam, sed.luminosity, label=f"T={T} K", color=color)

plt.title("Casey12")

plt.xlabel(r"$ \lambda / \AA$")
plt.ylabel(r"$ L_{\nu} / L_{bol}$")

plt.xlim(1e4, 1e8)
plt.ylim(10**-3.0, 1.5)

plt.legend()
plt.show()

Adding CMB heating

At high redshift, the effect of CMB heating is important. To account for this any dust emission model can incorporate the impact of CMB heating (cmb_heating) at a given redshift (\(z\)). Below we compare the dust emission spectrum from the Casey+12 model with and without CMB heating:

for ii, T in enumerate([25, 50, 100, 200]):
    color = plt.cm.tab10(ii)
    model = Casey12(temperature=T * K, emissivity=2.0, alpha=1.4)
    model_cmb = Casey12(
        temperature=T * K,
        emissivity=2.0,
        alpha=1.4,
        cmb_heating=True,
        redshift=10,
    )
    sed = model.get_spectra(lam)
    sed_cmb = model_cmb.get_spectra(lam)
    L_ir_ratio = sed_cmb.measure_window_luminosity(
        window=[8, 1000] * um
    ) / sed.measure_window_luminosity(
        window=[8, 1000] * um
    )  # same as model_cmb.cmb_factor
    plt.scatter(
        T,
        L_ir_ratio,
        label=f"T={np.around(model_cmb.temperature_z, 2)}",
        color=color,
    )

plt.axhline(y=1)
plt.grid(ls="dotted")
plt.xlabel("Input dust temperature")
plt.ylabel(r"L$_{\rm IR}$ ratio")
plt.legend()
plt.show()

IR Templates

Draine & Li 2007 dust models

The more complex model present in Draine & Li 2007 is also available; this requires a grid to be loaded in.

The grids are parameterised using the Umin, qpah and alpha, which denote the star light intensity heating majority of the dust, PAH fraction and power-law index (alpha \(\ne\) 1), see Draine & Li 2007 for more details. There are other free parameters the user provides, such as the dust mass, dust luminosity, gamma (1-gamma is the fraction of starlight exposed to intensity Umin). All the dust spectra in the grid have Umax=1E7.

We also follow the approach in Magdis+2012, by assuming that there is a fixed relationship between dust luminosity and dust mass, by default Ldust/Mdust = 125. This makes the free parameter space smaller.

Begin by reading in the DL07 grids, which have been created by downloading the ASCII DL07 files and running

from synthesizer.utils import process_dl07_to_hdf5
process_dl07_to_hdf5()

However, you can download these preprocessed grids using the synthesizer-download command line tool, this will place the grid in your grids directory.

synthesizer-download --dust-grid

from synthesizer.grid import Grid

grid_name = "draine_li_dust_emission_grid_MW_3p1.hdf5"
grid = Grid(grid_name)
print(grid.axes)

['qpah', 'umin', 'alpha']

Below we show how to use the DL07 model directly, and plot the dust emission spectrum for a number of dust luminosities. You can also pass these models to an EmissionModel and they will be used automatically.

for mdust in [1e7, 1e8, 5e8, 1e9, 5e9]:
    model = IR_templates(
        grid, mdust=mdust * Msun, ldust=1e11 * Lsun, verbose=False
    )
    sed = model.get_spectra(lam)

    # Normalise the SED to the bolometric luminosity
    sed._lnu /= sed._bolometric_luminosity

    # And plot...
    plt.loglog(
        sed.lam,
        sed.luminosity,
        label="{:.1e} Msun, <U>={:.2f}".format(mdust, model.u_avg),
    )
    plt.title("Draine & Li 2007")

plt.xlabel(r"$ \lambda / \AA$")
plt.ylabel(r"$ L_{\nu} / L_{bol}$")
plt.ylim(10**-3, 1)
plt.legend()
plt.show()

/home/runner/work/synthesizer/synthesizer/src/synthesizer/emission_models/dust/emission.py:726: RuntimeWarning: No hydrogen gas mass provided, assuming adust-to-gas ratio of 0.01
  self.dl07()
/home/runner/work/synthesizer/synthesizer/src/synthesizer/emission_models/dust/emission.py:726: RuntimeWarning: Gamma, Umin or alpha for DL07 model not provided, using default values
  self.dl07()
/home/runner/work/synthesizer/synthesizer/src/synthesizer/emission_models/dust/emission.py:726: RuntimeWarning: Computing required values using Magdis+2012 stacking results
  self.dl07()
/home/runner/work/synthesizer/synthesizer/src/synthesizer/emission_models/dust/emission.py:726: RuntimeWarning: Gamma not provided, choosing default gamma value as 5%
  self.dl07()
/opt/hostedtoolcache/Python/3.10.17/x64/lib/python3.10/site-packages/spectres/spectral_resampling.py:104: RuntimeWarning: Spectres: new_wavs contains values outside the range in spec_wavs, new_fluxes and new_errs will be filled with the value set in the 'fill' keyword argument (by default 0).
  warnings.warn(