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.plot(np.log10(sed.lam), sed.luminosity, label=f"{T} K")

plt.legend()
plt.show()

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

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.plot(np.log10(sed.lam), sed.luminosity, label=f"{T} K")

plt.legend()
plt.show()

[8, 1000] * um

unyt_array([   8, 1000], 'μm')

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=7)
    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.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 = "MW3.1"
grid_dir = "../../../tests/test_grid/"
grid = Grid(grid_name, grid_dir=grid_dir, read_lines=False)
print(grid.axes)

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Cell In[7], line 5
      3 grid_name = "MW3.1"
      4 grid_dir = "../../../tests/test_grid/"
----> 5 grid = Grid(grid_name, grid_dir=grid_dir, read_lines=False)
      6 print(grid.axes)

File ~/work/synthesizer/synthesizer/src/synthesizer/units.py:667, in accepts.<locals>.check_accepts.<locals>.wrapped(*args, **kwargs)
    664 for name, value in kwargs.items():
    665     kwargs[name] = _check_arg(units, name, value)
--> 667 return func(*args, **kwargs)

File ~/work/synthesizer/synthesizer/src/synthesizer/grid.py:167, in Grid.__init__(self, grid_name, grid_dir, spectra_to_read, read_lines, new_lam, lam_lims)
    164 self.parameters = {}
    166 # Get the axes of the grid from the HDF5 file
--> 167 self._get_axes()
    169 # Get the ionising luminosity (if available)
    170 self._get_ionising_luminosity()

File ~/work/synthesizer/synthesizer/src/synthesizer/grid.py:231, in Grid._get_axes(self)
    229 """Get the grid axes from the HDF5 file."""
    230 # Get basic info of the grid
--> 231 with h5py.File(self.grid_filename, "r") as hf:
    232     self.parameters = {k: v for k, v in hf.attrs.items()}
    234     # Get list of axes

File /opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/h5py/_hl/files.py:561, in File.__init__(self, name, mode, driver, libver, userblock_size, swmr, rdcc_nslots, rdcc_nbytes, rdcc_w0, track_order, fs_strategy, fs_persist, fs_threshold, fs_page_size, page_buf_size, min_meta_keep, min_raw_keep, locking, alignment_threshold, alignment_interval, meta_block_size, **kwds)
    552     fapl = make_fapl(driver, libver, rdcc_nslots, rdcc_nbytes, rdcc_w0,
    553                      locking, page_buf_size, min_meta_keep, min_raw_keep,
    554                      alignment_threshold=alignment_threshold,
    555                      alignment_interval=alignment_interval,
    556                      meta_block_size=meta_block_size,
    557                      **kwds)
    558     fcpl = make_fcpl(track_order=track_order, fs_strategy=fs_strategy,
    559                      fs_persist=fs_persist, fs_threshold=fs_threshold,
    560                      fs_page_size=fs_page_size)
--> 561     fid = make_fid(name, mode, userblock_size, fapl, fcpl, swmr=swmr)
    563 if isinstance(libver, tuple):
    564     self._libver = libver

File /opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/h5py/_hl/files.py:235, in make_fid(name, mode, userblock_size, fapl, fcpl, swmr)
    233     if swmr and swmr_support:
    234         flags |= h5f.ACC_SWMR_READ
--> 235     fid = h5f.open(name, flags, fapl=fapl)
    236 elif mode == 'r+':
    237     fid = h5f.open(name, h5f.ACC_RDWR, fapl=fapl)

File h5py/_objects.pyx:54, in h5py._objects.with_phil.wrapper()

File h5py/_objects.pyx:55, in h5py._objects.with_phil.wrapper()

File h5py/h5f.pyx:102, in h5py.h5f.open()

FileNotFoundError: [Errno 2] Unable to synchronously open file (unable to open file: name = '../../../tests/test_grid//MW3.1.hdf5', errno = 2, error message = 'No such file or directory', flags = 0, o_flags = 0)

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)
    plt.plot(
        np.log10(sed.lam),
        np.log10(sed.luminosity),
        label="{:.1e} Msun, <U>={:.2f}".format(mdust, model.u_avg),
    )

plt.xlabel("log10(lam/AA)")
plt.ylabel("log10(lnu/(erg/s))")
plt.ylim(40, 45)
plt.legend()
plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[8], line 3
      1 for mdust in [1e7, 1e8, 5e8, 1e9, 5e9]:
      2     model = IR_templates(
----> 3         grid, mdust=mdust * Msun, ldust=1e11 * Lsun, verbose=False
      4     )
      5     sed = model.get_spectra(lam)
      6     plt.plot(
      7         np.log10(sed.lam),
      8         np.log10(sed.luminosity),
      9         label="{:.1e} Msun, <U>={:.2f}".format(mdust, model.u_avg),
     10     )

NameError: name 'grid' is not defined