Normalizing the mass and/or optical depth of the dust in a simulation with multiple dust medium components can be confusing. The (advanced) thought experiment described on this page attempts to clarify the issues involved. To simplify the notation, we consider the total extinction of the dust at a single fixed wavelength. For the meaning of the employed symbols, refer to the table in the introduction of Material mixes (dust, electrons, gas).

A fundamental 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}. \]

Now consider a model consisting of two populations of a different dust type, 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}(\vec{r})\) 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}(\vec{r})\) and \(\rho_\text{B2}(\vec{r})\) 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!