 SKIRT 9 advanced radiative transfer in dusty systems
Dust properties

Absorption and scattering by dust plays an important role in most SKIRT simulations. This concept note describes how SKIRT obtains and uses the optical and calorimetric dust properties that are appropriate for a particular astrophysical model.

# Physical quantities

SKIRT uses two related quantities for defining the amount of material present in a given volume: number density and mass density. The kind of physical entity being counted by the number density and the conversion from number density to mass density depend on the type of material, as indicated in the table below.

Material type Entity counted Mass conversion
Dust hydrogen atom dust mass per hydrogen atom
Electrons electron electron mass
Gas hydrogen atom gas mass per hydrogen atom

The following table lists some key quantities used in SKIRT simulations, with their corresponding SI units.

 Symbol Units Description $$\lambda$$ $$\text{m}$$ Wavelength $$\overrightarrow{r}$$ $$\text{m}$$ Position $$\Delta s$$ $$m$$ Distance along a path $$V$$ $$\text{m}^3$$ Volume $$n$$ $$\#\,\text{m}^{-3}$$ Number density (of entities) $$\mu$$ $$\text{kg}\,\#^{-1}$$ Mass per entity $$\varsigma$$ $$\text{m}^2\,\#^{-1}$$ Cross section per entity $$\mathcal{N}=n\Delta s$$ $$\#\,\text{m}^{-2}$$ Number column density $$N=nV$$ $$\#$$ Number (of entities) $$\rho=n\mu$$ $$\text{kg}\,\text{m}^{-3}$$ Mass density $$\Sigma=n\mu\Delta s$$ $$\text{kg}\,\text{m}^{-2}$$ Mass column density $$M=n\mu V$$ $$\text{kg}$$ Mass $$\kappa=\varsigma/\mu$$ $$\text{m}^2\,\text{kg}^{-1}$$ Mass coefficient $$k=n\varsigma$$ $$\text{m}^{-1}$$ Opacity $$\tau=n\varsigma\Delta s$$ $$1$$ Optical depth

# Class hierarchy

The set of material properties corresponding to a particular transfer medium is represented in SKIRT as an instance of a MaterialMix subclass. The following diagram shows a portion of the class inheritance tree (connections starting with a triangle denote inheritance). The MaterialMix class hierarchy allows fundamentally different material types (e.g. dust, electrons, and hydrogen-dominated gas) to be implemented as part of a single framework. The common public interface offers the material properties required for tracing photon packets through a material, in other words, for processing absorption and scattering. Material properties that may be needed to calculate secondary emission spectra (e.g. thermal emission from dust grains) are offered by additional interfaces that depend on the material type.

## MaterialMix class

The MaterialMix class defines an abstract common interface to obtain the following basic material properties:

• the mass per entity $$\mu$$
• the absorption cross section per entity $$\varsigma^{\text{abs}}_{\lambda}$$
• the scattering cross section per entity $$\varsigma^{\text{sca}}_{\lambda}$$
• the total extinction cross section per entity $$\varsigma^{\text{ext}}_{\lambda} = \varsigma^{\text{abs}}_{\lambda} + \varsigma^{\text{sca}}_{\lambda}$$
• the scattering albedo $$\varpi_\lambda = \varsigma_{\lambda}^{\text{sca}} / \varsigma_{\lambda}^{\text{ext}}$$
• scattering phase function properties; depending on the supported scattering mode this may be:
• the assymmetry parameter $$g$$ for the Henyey-Greenstein phase function
• a custom phase function $$\Phi_\lambda(\cos\theta)$$ that depends only on the cosine of the scattering angle $$\theta$$
• a custom phase function $$\Phi_\lambda(\theta,\phi)$$ that depends on both scattering angles $$\theta$$ and $$\phi$$, and on the polarization state of the incoming radiation
• the equilibrium temperature $$T_{\text{eq}}$$ for a given embedding radiation field, using average material properties and assuming local thermal equilibrium conditions

The SKIRT photon packet lifecycle allows a material mix to support one of three scattering modes:

• HenyeyGreenstein: the material mix supplies the assymmetry parameter $$g$$ for the Henyey-Greenstein phase function. For a value of $$g=0$$, isotropic scattering is implemented directly rather than subsituting zero into the Henyey-Greenstein phase function.
• MaterialPhaseFunction: the material mix implements a custom phase function that depends only on the cosine of the scattering angle, $$\cos \theta$$. This mode does not support polarized radiation.
• SphericalPolarization: the material mix supports polarization through scattering by spherical particles. In this mode, the phase function depends on the polarization state of the incoming radiation, and the polarization state of the outgoing radiation must be updated appropriately. The phase function depends on both scattering angles $$\theta$$ and $$\varphi$$.

## DustMix class

The DustMix class implements the common interface defined by the MaterialMix base class for obtaining optical material properties. Depending on the scattering mode (selected by each subclass), the DustMix setup machinery requests the relevant optical dust properties from the subclass, and caches this information for later use. During the simulation, these properties can then be served up without further access to the subclass.

The properties served through this common interface correspond to the properties of a single grain population that is representative for the complete dust mix. In the context of tracking photon paths through a dusty medium, using integrated absorption and scattering cross sections and Mueller matrix coefficients is mathematically exact. In other words, for scattering modes MaterialPhaseFunction and SphericalPolarization, the representative grain approach does not involve an approximation. However, the calculation of a representative scattering asymmetry parameter $$g$$ for use with the HenyeyGreenstein scattering mode does involve a non-exact averaging procedure. Because the Henyey-Greenstein scattering phase function is non-physical to begin with, using a single approximate $$g$$ value for the complete dust mix is usually considered to be acceptable.

The DustMix class in turn has three abstract subclasses, each implementing a different approach to obtaining the relevant dust properties:

• TabulatedDustMix implements a dust mix described by tabulated properties for a single representative grain using the Henyey-Greenstein scattering mode. The optical dust properties are read from a user-supplied input file or can be listed directly inside the configuration file. This class and its subclasses are provided mostly for testing and benchmarking purposes.
• SingleGrainDustMix implements a dust mix described by a single representative grain, with or without support for polarization by scattering. The optical dust properties are retrieved from stored table resources provided with SKIRT. Subclasses must merely provide the names of the relevant resource files.
• MultiGrainDustMix implements a dust mix described by one or more grain populations, each with their own grain composition and size distribution, and with or without support for polarization by scattering. The class offers facilities to its subclasses to add dust grain populations to the dust mix. Based on this information it calculates the representative optical properties required for the common interface.

The information offered by the SingleGrainDustMix and TabulatedDustMix classes is insufficient to accurately calculate dust emission spectra for the dust mixture. This is so because the emission spectrum is a nonlinear function of the grain size (among many other things), and thus a single grain cannot accurately represent a population with a (potentialy large) range of grain sizes. The MultiGrainDustMix class therefore implements additional functionality to allow calculating the emission spectrum for multiple grain populations with given optical and calorimetric properties and with a given size distribution.

## MultiGrainDustMix class

As noted in the previous section, the MultiGrainDustMix class implements a dust mix described by one or more grain populations, each with their own grain composition and size distribution. The class obtains the following information obtained for each of the dust grain populations added by the subclass: the absorption efficiencies $$Q^{\text{abs}}(\lambda,a)$$, the scattering efficiencies $$Q^{\text{sca}}(\lambda,a)$$, the scattering phase function asymmetry parameter $$g(\lambda,a)$$, the Mueller matrix coefficients $$S^\text{xx}(\lambda,a,\theta)$$, the bulk density $$\rho_{\text{bulk}}$$ of the grain material, and the properly normalized grain size distribution per hydrogen atom $$\Omega(a)=(\frac{\text{d}n_\text{D}}{\text{d}a})/n_\text{H}$$ in the range $$[a_\text{min},a_\text{max}]$$.

### Calculating representative grain properties

The basic representative grain properties expected by the common MaterialMix interface are calculated by integrating over the grain size distribution $$\Omega(a)$$ using a builtin logarithmic grid and accumulating over all grain populations $$c$$, using the formulas listed below.

The absorption and scattering cross sections per hydrogen atom $$\varsigma_{\ell}^{\text{abs}}$$ and $$\varsigma_{\ell}^{\text{abs}}$$ for the $$\ell$$'th wavelength are calculated using

$\varsigma_{\ell}^{\text{abs}} = \sum_c \int_{a_{\text{min},c}}^{a_{\text{max},c}} \Omega_c(a)\, Q^{\text{abs}}_c(\lambda_\ell,a)\, \pi a^2\, {\text{d}}a$

and

$\varsigma_{\ell}^{\text{sca}} = \sum_c \int_{a_{\text{min},c}}^{a_{\text{max},c}} \Omega_c(a)\, Q^{\text{sca}}_c(\lambda_\ell,a)\, \pi a^2\, {\text{d}}a.$

The Mueller matrix coefficients provided by the grain population are assumed to be expressed as a cross section (in arbitrary units). They are thus integrated over the size distribution without again multiplying by the grain cross section, i.e. using

$S^\text{xx}_{\ell,\text{t}} = \sum_c \int_{a_{\text{min},c}}^{a_{\text{max},c}} \Omega_c(a)\, S^\text{xx}_c(\lambda_\ell,a,\theta_\text{t})\, {\text{d}}a$

The representative asymmetry parameter $$g_{\ell}$$ is averaged over the scattering cross section and thus calculated using

$g_{\ell} = \frac{1}{\varsigma_{\ell}^{\text{sca}}} \sum_c \int_{a_{\text{min},c}}^{a_{\text{max},c}} \Omega_c(a)\, g_c(\lambda_\ell,a)\, Q^{\text{sca}}_c(\lambda_\ell,a)\, \pi a^2\, {\text{d}}a.$

The dust mass per hydrogen atom $$\mu$$ is calculated by integrating the bulk density over the size distribution,

$\mu = \sum_c \int_{a_{\text{min},c}}^{a_{\text{max},c}} \Omega_c(a)\, \rho_{\text{bulk},c}\, \frac{4\pi}{3}\, a^3\, {\text{d}}a.$

### Calculating dust emission

The representative grain properties described above and offered by the public MaterialMix interface supported by this class are insufficient to accurately calculate dust emission spectra for the dust mixture. This is so because the emission spectrum is a nonlinear function of (among many other things) the grain size, and thus a single grain cannot accurately represent a population with a (potentialy large) range of grain sizes. Furthermore, smaller dust grains are often not in local thermal equilibrium, and instead are heated stochastically by individual photon absorption events. Modeling emission for these grains involves a temperature probability distribution rather than just an equilibrium temperature. The calculation needs calorimetric properties of the grain material in addition to optical properties.

It is numerically intractable to handle every possible grain size seperately. Instead, the MultiGrainDustMix class discretizes the grain size distribution for each type of grain material into a number of consecutive size bins (on a logarithmic scale), and calculates the optical and calorimetric properties of a representative grain for each of these bins. The number of bins for each type of grain material can be configured by the user. A larger number of bins improves the accuracy of the dust emission spectra. On the other hand, the calculation time scales roughly linearly with the number of bins.

The MultiGrainDustMix class uses one of two methods to calculate the emissivity of the dust mix:

• Assuming local thermal equilibrium for each representative grain (size bin): this method is fast but inaccurate because the equilibrium assumption is usually not justified. See the EquilibriumDustEmissionCalculator class for more information.
• Calculating a temperature probability distribution for each representative grain (size bin) to take into account stochastically heated grains: this second method is substantially more accurate but also much slower. See the StochasticDustEmissionCalculator class for more information.

## Turn-key dust mixes

SKIRT offers a range of built-in "turn-key" dust mixes representing sets of dust properties that have been published by various authors, including for example:

• MeanInterstellarDustMix: represents a population of identical dust grains approximating those of a mixture that is appropriate for the typical interstellar dust medium, based on publications by Draine, Li, and Weingartner. This dust mix is provided primarily for testing and tutorial purposes. Because it is a single-grain dust mix, it should not be used for calculating accurate dust emission spectra.
• DraineLiDustMix: represents a dust mixture of silicate, graphite, and PAH dust grains designed by Draine & Li 2007 such that the global dust properties accurately reproduce the extinction curve of the Milky Way.
• ZubkoDustMix: represents a realistic dust mixture of bare silicate, graphite, neutral PAH and ionized PAH dust grains designed by Zubko, Dwek & Arendt 2004 such that the global dust properties accurately reproduce the extinction, emission and abundance constraints on the Milky Way.
• ThemisDustMix: represents the THEMIS model for dust in the diffuse interstellar medium described by Jones et al. 2017 and the references therein. The model includes amorphous silicates with forsterite-normative composition and with enstatite-normative composition, and amorphous carbonaceous dust grains.

Refer to the respective classes (and other DustMix subclasses) for more information.

## Configurable dust mixes

The ConfigurableDustMix class represents a fully user-configurable dust mix described by one or more dust grain populations. Specifically, the class can be configured with a list of GrainPopulation instances, each of which represents a particular dust grain population with configurable grain composition, grain size distribution, and size bin discretization. This is illustrated in the following diagram (connections starting with a triangle denote inheritance; connections starting with a diamond denote aggregation). A GrainComposition object represents the optical and calorimetric properties of a particular grain material, according to a certain model. A GrainSizeDistribution object represents a particular grain size distribution function. For each grain population, a ConfigurableDustMix object holds a GrainPopulation object, which in turn holds a GrainComposition object and a GrainSizeDistribution object, and specifies the number of size bins used for discretizing the size distribution.

In addition, the user can define the amount of dust in each grain population by specifying one of the following three quantities:

• an absolute dust mass per hydrogen atom;
• a ratio of dust mass per hydrogen mass;
• a proportionality factor on the size distribution.

For example, one could mimic an MRNDustMix with the following ski file configuration:

<ConfigurableDustMix scatteringType="HenyeyGreenstein">
<populations type="GrainPopulation">
<GrainPopulation numSizes="10" normalizationType="FactorOnSizeDistribution"
factorOnSizeDistribution="7.762471166e-31">
<composition type="GrainComposition">
<DraineSilicateGrainComposition/>
</composition>
<sizeDistribution type="GrainSizeDistribution">
<PowerLawGrainSizeDistribution minSize="0.005 micron" maxSize="0.25 micron"
exponent="3.5"/>
</sizeDistribution>
</GrainPopulation>
<GrainPopulation numSizes="10" normalizationType="FactorOnSizeDistribution"
factorOnSizeDistribution="7.413102413e-31">
<composition type="GrainComposition">
<DraineGraphiteGrainComposition/>
</composition>
<sizeDistribution type="GrainSizeDistribution">
<PowerLawGrainSizeDistribution minSize="0.005 micron" maxSize="0.25 micron"
exponent="3.5"/>
</sizeDistribution>
</GrainPopulation>
</populations>
</ConfigurableDustMix>


# Normalizing dust media components

It is an important objective for the dust configuration in SKIRT to be understandable and transparent. Multiple dust populations should combine in a simple "linear" manner regardless of how they are included in the medium system. Based on a concrete example, this section demonstrates how this objective can be achieved.

To simplify the notation, we consider the dust properties at a single fixed wavelength so that we can omit the wavelength index. We also focus the discussion on the total extinction to avoid the upper index, and we leave the scattering asymmetry parameter aside.

A core issue is that when combining dust populations, the cross sections and masses per hydrogen atom can simply be added:

$\varsigma=\sum_c\varsigma_{c} \quad;\quad \mu=\sum_c\mu_{c},$

but this is not true for the mass coefficients:

$\kappa=\frac{\varsigma}{\mu}=\frac{\sum_c\varsigma_{c}}{\sum_c\mu_{c}} \neq \sum_c\frac{\varsigma_{c}}{\mu_c}=\sum_c\kappa_{c}.$

## Normalization schemes

Consider a SKIRT configuration containing one or more dust media components, each with a specific geometry and a particular dust mix. Each dust mix has a number of grain populations. Let us denote dust components (i.e. geometries) with the index $$h$$ and grain populations for each component with the index $$c$$. The various dust mixes in this configuration thus define the cross sections per hydrogen atom $$\varsigma_{h,c}$$ and the dust masses per hydrogen atom $$\mu_{h,c}$$ for each dust component and dust population. In addition, the geometry for each dust component defines the dust mass density distribution $$\rho_{h}(\vec{r})$$ for that dust component, discretized over the dust grid, with arbitrary scaling (generally the distribution is normalized to unity).

The mass density distribution $$\rho_{h}(\vec{r})$$ can then be normalized to a given total mass $$M_h$$ or to a given optical depth $$\tau_h$$ along a specific path S using

$M_h = \int_\text{V}\rho_{h}(\vec{r})\,\text{d}V$

$\tau_h = \int_\text{S}\kappa_{h}\,\rho_{h}(\vec{r})\,\text{d}s$

where the derived dust mix property $$\kappa_{h}$$ is calculated according to

$\kappa_{h}=\frac{\sum_c\varsigma_{h,c}}{\sum_c\mu_{h,c}}.$

## A specific example

Now consider a configuration consisting of two dust populations, named 1 and 2, with identical spatial distribution. There are two distinct ways to configure this model in SKIRT:

• configuration A: a single dust component using a dust mix with the two populations 1 and 2;
• configuration B: two dust components with identical spatial distribution, the first using a dust mix with population 1 and the second using a dust mix with population 2.

Given appropriate normalization of the respective dust components, we expect the results of configurations A and B to be identical.

In both cases, the simulation obtains the cross sections $$\varsigma_1, \varsigma_2$$ and dust masses per hydrogen atom $$\mu_1, \mu_2$$ for each population. In configuration A there is a single dust component with mass density distribution $$\rho_\text{A}$$ and with mass coefficient $$\kappa_\text{A}=(\varsigma_1+\varsigma_2)/(\mu_1+\mu_2)$$. In configuration B there are two dust components with dust mass density distribution $$\rho_\text{B1}$$ and $$\rho_\text{B2}$$ and with mass coefficients $$\kappa_\text{B1}=\varsigma_1/\mu_1$$ and $$\kappa_\text{B2}=\varsigma_2/\mu_2$$.

### Normalization on mass

Assume that we are given the total normalization dust mass $$M$$ for configuration A. We'd like to find the normalization masses $$M_1$$ and $$M_2$$ for each of the dust components in configuration B so that the total optical depth along an arbitrary path is identical in both configurations.

In configuration A the normalization equation reads

$M = \int_\text{V}\rho_\text{A}(\vec{r})\,\text{d}V$

and the total optical depth along an arbitrary path P is given by

$\tau_\text{A} = \frac{\varsigma_1+\varsigma_2}{\mu_1+\mu_2} \int_\text{P}\rho_\text{A}(\vec{r})\,\text{d}s$

In configuration B the normalization equations read

$M_1 = \int_\text{V}\rho_\text{B1}(\vec{r})\,\text{d}V$

$M_2 = \int_\text{V}\rho_\text{B2}(\vec{r})\,\text{d}V$

and the total optical depth along an arbitrary path is given by

$\tau_\text{B} = \frac{\varsigma_1}{\mu_1} \int_\text{P}\rho_\text{B1}(\vec{r})\,\text{d}s + \frac{\varsigma_2}{\mu_2} \int_\text{P}\rho_\text{B2}(\vec{r})\,\text{d}s$

Since all geometries are identical, we can write $$\rho_\text{B1}(\vec{r})=b_1\rho_\text{A}(\vec{r})$$ and $$\rho_\text{B2}(\vec{r})=b_2\rho_\text{A}(\vec{r})$$ where $$b_1,b_2$$ are constants that don't dependent on $$\vec{r}$$. Requiring $$\tau_\text{A}=\tau_\text{B}$$ and $$M=M_1+M_2$$ then leads to the system of equations

$\frac{\varsigma_1+\varsigma_2}{\mu_1+\mu_2} = \frac{\varsigma_1}{\mu_1}b_1 + \frac{\varsigma_2}{\mu_2}b_2 \quad;\quad 1=b_1+b_2$

in the unknowns $$b_1$$ and $$b_2$$. As can be easily verified by substitution, the solution of this system is $$b_1 = \frac{\mu_1}{\mu_1+\mu_2}$$ and $$b_2 = \frac{\mu_2}{\mu_1+\mu_2}$$. From the mass normalization equations above we see that $$M_1/M=b_1$$ and $$M_2/M=b_2$$ so that

$M_1 = \frac{\mu_1}{\mu_1+\mu_2}M \quad;\quad M_2 = \frac{\mu_2}{\mu_1+\mu_2}M.$

In other words, the normalization mass must be distributed over the dust components proportional to the dust mass of each dust population. A rather intuitive result!

### Normalization on optical depth

Assume that we are given the optical depth $$\tau$$ along a specific path S for normalizing configuration A. We'd like to find the optical depths $$\tau_1$$ and $$\tau_2$$ (along the same path) for normalizing each of the dust components in configuration B so that the total optical depth along an arbitrary path is identical in both configurations.

In configuration A the normalization equation reads

$\tau = \frac{\varsigma_1+\varsigma_2}{\mu_1+\mu_2} \int_\text{S}\rho_\text{A}(\vec{r})\,\text{d}s$

and the total optical depth along an arbitrary path is given by

$\tau_\text{A} = \frac{\varsigma_1+\varsigma_2}{\mu_1+\mu_2} \int_\text{P}\rho_\text{A}(\vec{r})\,\text{d}s$

In configuration B the normalization equations read

$\tau_1 = \frac{\varsigma_1}{\mu_1} \int_\text{S}\rho_\text{B1}(\vec{r})\,\text{d}s$

$\tau_2 = \frac{\varsigma_2}{\mu_2} \int_\text{S}\rho_\text{B2}(\vec{r})\,\text{d}s$

and the total optical depth along an arbitrary path is given by

$\tau_\text{B} = \frac{\varsigma_1}{\mu_1} \int_\text{P}\rho_\text{B1}(\vec{r})\,\text{d}s + \frac{\varsigma_2}{\mu_2} \int_\text{P}\rho_\text{B2}(\vec{r})\,\text{d}s$

Since all geometries are identical, we can write $$\rho_\text{B1}(\vec{r})=b_1\rho_\text{A}(\vec{r})$$ and $$\rho_\text{B2}(\vec{r})=b_2\rho_\text{A}(\vec{r})$$ where $$b_1,b_2$$ are constants that don't dependent on $$\vec{r}$$. Requiring $$\tau_\text{A}=\tau_\text{B}$$ and $$\tau=\tau_1+\tau_2$$ both lead to the equation

$\frac{\varsigma_1+\varsigma_2}{\mu_1+\mu_2} = \frac{\varsigma_1}{\mu_1}b_1 + \frac{\varsigma_2}{\mu_2}b_2$

in the unknowns $$b_1$$ and $$b_2$$. As can be easily verified by substitution, a solution of this equation is $$b_1 = \frac{\mu_1}{\mu_1+\mu_2}$$ and $$b_2 = \frac{\mu_2}{\mu_1+\mu_2}$$. If we also require that the total dust mass in the system remains the same, then this is the only solution (as shown in the previous subsection). Otherwise there is a family of solutions and we can arbitrarily select this one.

From the optical depth normalization equations above we see that $$\tau_1=b_1\frac{\varsigma_1}{\mu_1}\left(\frac{\varsigma_1+\varsigma_2}{\mu_1+\mu_2}\right)^{-1}\tau$$ and $$\tau_2=b_2\frac{\varsigma_2}{\mu_2}\left(\frac{\varsigma_1+\varsigma_2}{\mu_1+\mu_2}\right)^{-1}\tau$$ so that after some basic algebra we obtain

$\tau_1 = \frac{\varsigma_1}{\varsigma_1+\varsigma_2}\tau \quad;\quad \tau_2 = \frac{\varsigma_2}{\varsigma_1+\varsigma_2}\tau.$

In other words, the optical depth must be distributed over the dust components proportional to the cross section of each dust population. A rather intuitive result!