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 et al. (2012). By default a dust emission model provides a normalised spectrum.

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

plt.title("Blackbody")

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

plt.ylim(10**-3.0, 1.5)

plt.legend()
plt.show()

/opt/hostedtoolcache/Python/3.10.16/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.16/x64/lib/python3.10/site-packages/unyt/array.py:1972: RuntimeWarning: overflow encountered in multiply
  out_arr = func(

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

plt.title("Greybody")

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

plt.ylim(10**-3.0, 1.5)

plt.legend()
plt.show()

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

plt.title("Casey12")

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

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 T in [10, 25, 50, 100, 1000]:
    model = Casey12(T * K, 1.6, 2.0)
    model_cmb = Casey12(T * K, 1.6, 2.0, cmb_heating=True, redshift=15)
    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"{np.around(model_cmb.temperature_z, 2)}"
    )

plt.axhline(y=1)
plt.grid(ls="dotted")
plt.xlabel("Input dust temperature")
plt.ylabel(r"L$_{\rm IR}$ ratio")
plt.yscale("log")
plt.xscale("log")
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.

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.

synthesizer-download -d /path/to/destination --dust-grid

from synthesizer.grid import Grid

grid_name = "draine_li_dust_emission_grid_MW_3p1.hdf5"
grid_dir = "../../../tests/test_grid/"
grid = Grid(grid_name, grid_dir=grid_dir, read_lines=False)
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:732: 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:732: 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:732: RuntimeWarning: Computing required values using Magdis+2012 stacking results
  self.dl07()
/home/runner/work/synthesizer/synthesizer/src/synthesizer/emission_models/dust/emission.py:732: RuntimeWarning: Gamma not provided, choosing default gamma value as 5%
  self.dl07()
/opt/hostedtoolcache/Python/3.10.16/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(