Chapter 3
The Grain Properties of Nearby Galaxies

In this chapter, we review how the models presented in Sect. 2.3 are used to derive the grain properties of nearby galaxies. The term “dust properties” is vague. In the literature, it often indistinctively encompasses the three following categories (Galliano, Galametz, & Jones, 2018).

The dust mixture constitution

is characterized by:

• the chemical composition of the bulk material and its stoichiometry (Sect. 1.1.4);
• the structure of the grains (crystalline, amorphous, porous, aggregated, etc.; Sect. 1.2.2.4);
• the presence of heterogeneous inclusions (Sect. 1.2.2.4);
• the presence of organic and/or icy mantles (Sect. 1.2.2.4);
• the shape of the grains (Sect. 1.2.2.5);
• their size distribution (Sect. 2.3.2);
• their abundance relative to the gas (Sect. 2.2.3).

The dust physical conditions

are the state a grain enters, when exposed to a particular environment:

• thermal excitation due to (i) radiative heating (equilibrium or stochastic; Sects. 1.2.4.2 – 1.2.4.3); and/or (ii) collisional heating in a hot plasma (Sect. 1.2.4.4);
• grain charging by exchange of electrons with the gas (Sect. 3.3.1.1);
• alignment of elongated grains with the magnetic field (Sect. 1.2.2.5);
• grain rotation (Sect. 1.2.2.5; Sect. 2.2.2.3).

The dust observables

arise from a grain mixture experiencing a particular set of physical conditions.

Emission

of partially polarized components (Sect. 2.2.2.2): (i) a thermal continuum (IR to mm; Sect. 2.2.2.1); (ii) molecular and solid-state features (MIR; Sect. 2.2.2.1); (iii) a possible microwave emission (cm; Sect. 2.2.2.3); (iv) a possible luminescence (visible; Sect. 2.2.2.3).

Absorption

of the light from a background source by (Sect. 2.2.1): (i) a continuum (X-ray to MIR; Sect. 2.2.1.1); (ii) atomic, molecular and solid-state features (X-rays, UV and MIR; Sect. 2.2.1.2), including DIBs (visible; Sect. 2.2.1.5) and ices (MIR; Sect. 2.2.1.2). The escaping light can be partially polarized as a result (Sect. 2.2.1.4).

Scattering

of the light from a bright source in our direction, and its polarization (X-rays to NIR; Sect. 2.2.1).

Depletion

patterns seen through the gas-phase elemental abundances (Sect. 2.2.3).

3.1 Spectral Energy Distribution Modeling

SED modeling is one of the main methods to empirically derive the dust properties of a region or a galaxy (cf. e.g. Galliano et al., 2018, for a review). The inherent complexity of astrophysical sources requires to account for the diversity of physical conditions within the studied region. The treatment of radiative transfer, even in an extremely approximated fashion, is thus necessary.

3.1.1 Radiative Transfer

Radiative transfer is the method solving the propagation of multiple rays of light, emitted by one or several sources, through a macroscopic heterogeneous medium. It accounts for the scattering, absorption and emission, at each point and along each direction, in the studied region.

3.1.1.1 Definition of the Main Radiative Transfer Quantities

Radiative transfer deals with several quantities that are often improperly defined or mixed together in the literature: intensity, flux, emissivity, brightness, etc. We have already seen some of these quantities in Sect. 1.2.4. We now define them and explicit their differences (cf. Chap. 1 of Rybicky & Lightman, 1979, for a complete review). In what follows, we assume stationary systems. The time variable, $t$, is used only to denote constant rates.

The moments of the specific intensity. The primary radiative transfer quantity is the specific intensity or brightness (cf. Fig. 3.1.a):

The specific intensity is the electromagnetic energy, $E$, per unit time, $t$, area, $A$, solid angle, $\Omega$, and frequency, $\nu$. It therefore quantifies the infinitesimal power carried by a monochromatic light ray. This quantity depends on the position in the region,$\overrightarrow{r}$, and on the direction of propagation, $\left(𝜃,\varphi \right)$. We adopt the spherical coordinate conventions used in physics (Fig. 3.1; cf. Appendix C.1):

• the polar angle, $0\le 𝜃<\pi$, is the angle with the $z$ axis;
• the azimuthal angle, $0\le \varphi <2\pi$, is the rotation angle in the $\left(x,y\right)$ plane;
• the solid angle element is $d\Omega \equiv d\cos <!--FUNCTION\; APPLICATION-->𝜃d\varphi =\sin <!--FUNCTION\; APPLICATION-->𝜃d𝜃d\varphi$, with:

  $${\iint }_{sphere}d\Omega ={\int \; <!--nolimits-->}_0^{2\pi }d\varphi {\int \; <!--nolimits-->}_{-1}^{1}d\cos <!--FUNCTION\; APPLICATION-->𝜃=4\pi .$$

  (3.2)

The first moments of the specific intensity, relative to its angular distribution, are physically meaningful.

The mean intensity

is the zeroth order moment of Eq. ($3.1$):

$$J_{\nu }\left(\nu ,\overrightarrow{r}\right)\equiv \frac{1}{4\pi }{\iint }_{sphere}{I}_{\nu }\left(\nu ,\overrightarrow{r},𝜃,\varphi \right)d\Omega .$$

(3.3)

It is the specific intensity averaged over all directions. It is often used to quantify the ISRF, $4\pi J_{\nu }$ accounting for rays coming from all directions. In the isotropic case, we have $I_{\nu }\left(\nu ,\overrightarrow{r},𝜃,\varphi \right)=J_{\nu }\left(\nu ,\overrightarrow{r}\right)$, $\forall <!--FUNCTION\; APPLICATION-->\left(𝜃,\varphi \right)$.

The net flux

(cf. Fig. 3.1.b) is the first order moment of Eq. ($3.1$):

$$F_{\nu }\left(\nu ,\overrightarrow{r}\right)\equiv {\iint }_{sphere}{I}_{\nu }\left(\nu ,\overrightarrow{r},𝜃,\varphi \right)\cos <!--FUNCTION\; APPLICATION-->𝜃d\Omega .$$

(3.4)

It represents the monochromatic power per unit area passing through a surface element perpendicular to $ẑ$. The $\cos <!--FUNCTION\; APPLICATION-->𝜃$ factor is there to account for the reduction of the density of rays that are not perpendicular to the surface. In the isotropic case, $F_{\nu }=J_{\nu }2\pi {\mathrm{\int \; <!--nolimits-->}<!--FUNCTION\; APPLICATION-->}_{-1}^1\cos <!--FUNCTION\; APPLICATION-->𝜃d\cos <!--FUNCTION\; APPLICATION-->𝜃=0$, because there is the same amount of flux passing through the area in both directions.

The radiation pressure

(cf. Fig. 3.1.c) is the second order moment of Eq. ($3.1$):

$$p_{\nu }\left(\nu ,\overrightarrow{r}\right)\equiv \frac{1}{c}{\iint }_{sphere}{I}_{\nu }\left(\nu ,\overrightarrow{r},𝜃,\varphi \right){\cos <!--FUNCTION\; APPLICATION-->}^2𝜃d\Omega$$

(3.5)

It is the momentum flux carried by the photons, as $p=h\nu ∕c$ is the momentum of a single photon. The ${\cos <!--FUNCTION\; APPLICATION-->}^2𝜃$ term comes from two sources: (i) one $\cos <!--FUNCTION\; APPLICATION-->𝜃$ comes from the reduced fraction of inclined rays, similar to the flux; (ii) the second $\cos <!--FUNCTION\; APPLICATION-->𝜃$ factor comes from the fact the pressure is the momentum vector component that is perpendicular to the surface.

Specific energy density. The radiative energy within a volume element at a given time is the specific energy density:

The second equality comes from the fact that $dV=cdtdA$, $dV$ being a volume element (cf. Fig. 3.1.d). The specific intensity is a power per unit area, whereas the energy density is an energy per unit volume. Both quantities are linked, combining Eq. ($3.1$) and Eq. ($3.6$):

If we integrate Eq. ($3.6$) over all directions, we get, using Eq. ($3.3$):

Emission coefficient and emissivity. The monochromatic emission coefficient is the power radiated in a given direction, per unit volume and frequency:

We have also seen, in Sect. 1.2.4, the emissivity:

where $\rho$ is the mass density of the ISM, and $dm=\rho dV$, its mass element. The factor $4\pi$, in Eq. ($3.10$), makes ${𝜖}_{\nu }$ the solid angle fraction of monochromatic power emitted in a given direction, per unit mass. In Sect. 1.2.4, we were considering the volume and mass of a grain, whereas, here, we are considering the volume and mass of the ISM. Combining Eq. ($3.9$) and Eq. ($3.10$), we get:

Extinction coefficient and opacity. The amount of specific intensity absorbed and scattered along an infinitesimal path length, $dl$, in the direction $\left(𝜃,\varphi \right)$, is the extinction coefficient, $\alpha \ge 0$, defined such that:

Similarly to the convention we have adopted for $\kappa$ in Sect. 1.2.4, we pose $\alpha \equiv {\alpha }_{ext}={\alpha }_{abs}+{\alpha }_{sca}$ to distinguish absorption and scattering.. Eq. ($3.12$) can be expressed with microscopic quantities, assuming the absorbers and scatterers have a cross-section $C_{ext}\left(\nu ,\overrightarrow{r}\right)$ and a density $n\left(\overrightarrow{r}\right)$:

If the composition of the ISM is homogeneous, then $C_{ext}\left(\nu ,\overrightarrow{r}\right)=C_{ext}\left(\nu \right)$. The opacity, that we have seen in Sect. 1.2.4, is related to $\alpha$ by:

Mean free path. The mean free path of a photon with frequency, $\nu$, at position $\overrightarrow{r}$ can be defined as:

It is the average length a photon will be able to travel before being absorbed or scattered. We can make the same remark as for the emissivity: in Sect. 1.2.4, we were considering the cross-section of dust particles per mass of grain, whereas Eq. ($3.14$) gives the cross-section of the whole ISM per mass of ISM. Table 3.1 gives typical values of $l_{mean}$ in a homogeneous medium, assuming the dust constitution of the THEMIS model (cf. Sect. 2.3.2.2).

3.1.1.2 The Radiative Transfer Equation

The radiative transfer equation accounts for the variation of the specific intensity under the effects of absorption, scattering and emission (cf. Steinacker et al., 2013, for a review in the case of a dusty medium). It is schematically represented in Fig. 3.2. This equation can be written:

Scattering in the sightline:

this term is the integral of ${\alpha }_{sca}{I}_{\nu }$ (i.e. the scattered intensity) over the phase function, $\Phi$ (Eq. $1.40$). This expression depends on the relative angle, ${𝜃}^{\prime }$, between the incident rays and the scattered direction $\left(𝜃,\varphi \right)$. This is the term that makes this equation an integro-differential equation, coupling all the directions together. This is why, numerical methods are required to solve Eq. ($3.16$). If we assume isotropic scattering (i.e. $⟨\cos <!--FUNCTION\; APPLICATION-->𝜃⟩=0$, corresponding to the Rayleigh regime; cf. Sect. 1.2.2.3), this term simplifies and becomes: ${\alpha }_{sca}\left(\nu ,\overrightarrow{r}\right)J_{\nu }\left(\nu ,\overrightarrow{r}\right)$.

Emission:

we have assumed that the dust and stellar emissions are both isotropic, (i.e. independent of $𝜃$ and $\varphi$). This is a reasonable assumption in the ISM. If we also assume that there is only one grain species, and that this species is at equilibrium with the radiation field, we can explicit: $j_{\nu }^{dust}\left(\nu \right)=\alpha \left(\nu \right)B_{\nu }\left(\nu ,T_{eq}\right)$. This simplification still contains the problem that to determine the equilibrium temperature, $T_{eq}$, we need to integrate ${\alpha }_{abs}{I}_{\nu }$ over all directions, and all frequencies.

Eq. ($3.16$) can be rewritten, using the optical depth, $\tau$, defined such that $d\tau \left(\nu \right)=\alpha \left(\nu ,l\right)dl$, or:

Eq. ($3.15$) implies that $\tau \left(l_{mean}\right)=1$. At a given wavelength, a medium is: (i) optically thin or transparent, if $\tau «1$; and (ii) optically thick or opaque, if $\tau »1$. Replacing $dl$ by $d\tau ∕\alpha$ in Eq. ($3.16$), we obtain:

where the source function, $S_{\nu }$, includes all the terms added to the specific intensity:

We will discuss exact numerical solutions to this equation in Sect. 3.1.1.4. For now, we discuss trivial solutions, when some processes are assumed negligible.

Propagation in vacuum. Let’s assume we have a star of radius, $R$, and surface temperature, $T_{\star }$, located at $\overrightarrow{r}=\overrightarrow{0}$. The flux at $r=R$ has to be integrated only over the hemisphere where a surface element of the star emits: $F_{\nu }\left(R\right)=2\pi {\mathrm{\int \; <!--nolimits-->}<!--FUNCTION\; APPLICATION-->}_0^1{B}_{\nu }\left(T_{\star }\right)\cos <!--FUNCTION\; APPLICATION-->𝜃d\cos <!--FUNCTION\; APPLICATION-->𝜃=\pi B_{\nu }\left(T_{\star }\right)$. If there is no ISM around the star, Eq. ($3.16$) simply becomes $d{I}_{\nu }∕dl=0$. The solution is thus $I_{\nu }=B_{\nu }\left(T_{\star }\right)$ along the directions coming from the star and 0 in all other directions. At an arbitrary distance, $r$, from the star, the solid angle it occupies is ${\Omega }_{\star }=R^2∕r^2$. The flux is thus $F_{\nu }\left(r\right)=\pi B_{\nu }\left(T_{\star }\right)R^2∕r^2$, which is the classic $1∕r^2$ dilution of the flux.

Emission only. Let’s assume we are observing, at submm wavelengths, a molecular cloud constituted of equilibrium grains at $T\simeq 10$ K, with opacity $\kappa$. At these wavelengths, the extinction is negligible. Eq. ($3.16$) is therefore simply $d{I}_{\nu }∕dl=\rho \kappa B_{\nu }\left(T\right)$ within the cloud. The solution is thus $I_{\nu }\left(l\right)=\kappa B_{\nu }\left(T\right){\mathrm{\int \; <!--nolimits-->}<!--FUNCTION\; APPLICATION-->}_{l_0}^l\rho \left(l\right)dl$, where $l_0$ is the position of the edge of the cloud along the direction of the sightline. If the density is constant, we get $I_{\nu }\left(l\right)=\rho \kappa B_{\nu }\left(T\right)\times \left(l-l_0\right)$, which can be simplified as $I_{\nu }\left(l\right)=\tau \left(l\right)B_{\nu }\left(T\right)$. This expression of the brightness is often used in radio-astronomy.

Absorption only. Let’s assume we are observing a background star through a cold molecular cloud, in the MIR. At these wavelengths, the albedo is close to 0 (Sect. 1.2.2.3). If we make the assumption that the background star is much brighter than the thermal emission of the cloud, Eq. ($3.16$) simply becomes $d{I}_{\nu }∕dl=-\alpha I_{\nu }$. The solution is therefore $I_{\nu }\left(l\right)=I_{\nu }^{\star }\exp <!--FUNCTION\; APPLICATION-->\left[-\tau \left(l\right)\right]$ along the direction coming from the star and 0 in all other directions.

Emission and absorption. We can merge the two previous cases. It could correspond to a hot molecular cloud, observed at MIR wavelengths. Its thermal emission is absorbed by the cloud itself. The solution of Eq. ($3.18$), in this case, is:

If the cloud contains a single grain species at temperature $T$, the solution becomes:

If we look in a direction away from the background star, we get the classical self-absorption formula:

This solution is displayed in Fig. 3.3. Eq. ($3.22$) has the following two limit regimes.

Optically thin:

if $\tau «1$, $I_{\nu }\simeq \tau B_{\nu }\left(T\right)$, which is the “emission only” solution. The emission from the cloud is a grey body (Sect. 1.2.4).

Optically thick:

if $\tau »1$, $I_{\nu }\simeq B_{\nu }\left(T\right)$, which means the cloud is not anymore a grey body, but a perfect one.

Scattering and absorption with central illumination. Let’s assume that we have a homogeneous spherical cloud, of radius $R_{cl}$, and a central isotropically illuminating source, with monochromatic luminosity, $L_{\nu }^{\star }\left(\lambda \right)$. If $L_{\nu }^{esc}\left(\lambda \right)$ is the monochromatic luminosity escaping the cloud, the escape fraction can be defined as:

The second equality defines ${\tau }_{e\mathrm{ff}}$ as the effective optical depth of the medium.

Optically thick:

In the case of pure isotropic scattering ($g=⟨\cos <!--FUNCTION\; APPLICATION-->𝜃⟩\simeq 0$), it can be shown that the net displacement of a photon after $N_{sca}$ interactions is: $l_{e\mathrm{ff}}=\sqrt{N_{sca}}{l}_{mean}$ (cf. Chap. 1 of Rybicky & Lightman, 1979). The probability a photon will be absorbed, at the end of a free path, is $1-\tilde{\omega }$. The mean number of free paths can thus be estimated by $N_{sca}\left(1-\tilde{\omega }\right)=1$, or: $N_{sca}=1∕\left(1-\tilde{\omega }\right)$. The optical depth of the cloud is $\tau \simeq R_{cl}∕l_{mean}$ (Eq. $3.15$ and Eq. $3.17$). The same way, the effective optical depth can be written ${\tau }_{e\mathrm{ff}}\simeq R_{cl}∕l_{e\mathrm{ff}}$. We thus get:

$${\tau }_{e\mathrm{ff}}\simeq \sqrt{1-\tilde{\omega }}\tau .$$

(3.24)

Optically thin:

in this case, there is a low probability of interaction. Eq. ($3.23$) tells us that ${\tau }_{e\mathrm{ff}}$ accounts only for photons that have been absorbed. We therefore simply have (cf. Chap. 1 of Rybicky & Lightman, 1979):

$${\tau }_{e\mathrm{ff}}=\left(1-\tilde{\omega }\right)\tau .$$

(3.25)

This formula simply subtracts among the few photons that may have interacted with the grains those which have been scattered. Contrary to Eq. ($3.24$), it is valid for any value of $g$.

Városi & Dwek (1999, hereafter VD99) have proposed an empirical approximation to interpolate the two regimes of Eq. ($3.24$) and Eq. ($3.25$), resulting in the following escape fraction:

with:

They have benchmarked this approximation with a Monte-Carlo radiative transfer model (cf. Sect. 3.1.1.4). Both are in very good agreement in most of the astrophysically relevant parameter space. This solution is displayed in Fig. 3.4.a.

Scattering and absorption with uniform illumination. VD99 have also proposed an approximation in the case where the internal illumination of the cloud is uniform. They start from the escape probability of a homogeneous sphere with uniform illumination, without scattering, given by Osterbrock (1989):

A demonstration of this formula is given in Appendix C of VD99. VD99 find, with a recursive argument, that the following expression is in relatively good agreement with a Monte-Carlo radiative transfer model:

VD99 show that Eq. ($3.29$) perfectly agrees with the exact solution for a particular value of $g\left(\tau \right)$. The largest discrepancies, of about $20\%$, are found at high optical depth, for $g$ close to 0. Eq. ($3.29$) is demonstrated in Fig. 3.4.b. We can see that, at a given ${\tau }_V$, this geometry has an overall larger escape fraction than central illumination.

Scattering and absorption with external illumination. There are also expressions in the case of a homogeneous spherical cloud, externally illuminated by an isotropic radiation field. VD99 derive the following approximation for the absorbed fraction:

The absorbed fraction is indeed more relevant in the case of external illumination, as it would be difficult to observationally separate the escaping radiation from the cloud and the ambient ISRF. On the contrary the absorbed fraction is meaningful if one wants to evaluate the heating of the cloud. This approximation is demonstrated in Fig. 3.4.c.

The three formulae for spherical clouds, given in Eq. ($3.26$), Eq. ($3.29$) and Eq. ($3.30$), do not allow us to model the internal heating of the cloud. These expressions are indeed global escape and absorbed fractions, but they do not account for the gradient of illumination within the cloud that would lead to a gradient of heating rate.

3.1.1.3 Approximations for Clumpy Media

The ISM is a highly heterogeneous medium, with contrast densities of several orders of magnitude. A useful approximation is the clumpy medium, composed of: (i) a diffuse, uniform InterClump Medium (ICM), characterized by its density, $n_{ICM}$; and (ii) dense, spherical clumps, with density, $n_C$, radius, $r_C$, and volume filling factor, $f_V$.

Effective optical depth of a clumpy medium. Let’s assume that we are looking at a background star through a cloud, and that the albedo of the grains is negligible. Along a given sightline, we have seen in Sect. 3.1.1.2 that the brightness in the direction of the star is $I_{\nu }=I_{\nu }^{\star }\times \exp <!--FUNCTION\; APPLICATION-->\left(-\tau \right)$. If the clumpy structure of the cloud is unresolved, the brightness measured in the telescope beam can be written as the sum of $N$ sightlines, some passing through the ICM, others through the clumps:

This is the general expression. In the case of an homogeneous medium, Eq. ($3.31$) simplifies: $I_{\nu }=I_{\nu }^{\star }\exp <!--FUNCTION\; APPLICATION-->\left(-{\tau }_{hom}\right)$, where ${\tau }_{hom}$ is the optical depth of the homogeneous medium. In order to have the same dust mass and opacity as in the clumpy medium, we need to have:

In the homogeneous medium, $\rho \times L$ is indeed the mass surface density of a cloud of depth $L$. If there is a statistical distribution of clumps in the beam, it must be identical in the clumpy medium. Thus, the brightness of the homogeneous medium is:

Invoking the arithmetic-mean/geometric-mean inequality ¹ (e.g. page 456 of Cauchy, 1821), we conclude that: $I_{\nu }^{clumpy}\ge I_{\nu }^{hom}$. An important consequence of this result is that, from an observational point of view, we can miss a large mass of dust hidden in clumps. VD99 discuss this result in more detail.

The mega-grains approximation. Neufeld (1991) proposed a simple approach to explain the leakage of Ly-$\alpha$ photons by galaxies. He treated dusty gas clumps, in an empty ICM, as large grains with their own albedo and asymmetry parameter. This idea was then further developed by Hobson & Padman (1993), in cases where the ICM is non empty. They named it the mega-grains approximation, as clumps are treated as grains, although they have macroscopic sizes. They applied this approach to a clumpy infinite slab, externally illuminated on one side, and compared the results to a Monte-Carlo radiative transfer model. VD99 then refined some of the expressions of Hobson & Padman (1993) and applied them to the case of a spherical clumpy cloud, with the three types of illuminations we have discussed in Sect. 3.1.1.2: (i) central; (ii) uniform; and (iii)  external. VD99 systematically benchmarked their results with a Monte-Carlo radiative transfer code. We briefly review their results in the rest of this section.

Escape fractions for a clumpy medium. VD99 derived a series of expressions, based on Eq. ($3.26$), Eq. ($3.29$) and Eq. ($3.30$), but replacing grain properties by effective mega-grains properties. These analytical approximations are all summarized in Sect. 5 of VD99. We demonstrate these analytical expressions for the three types of illuminations in Fig. 3.5. Overall, they provide a good agreement with Monte-Carlo radiative transfer calculations. They are also very easy to compute. The weakest point concerns the treatment of the grain heating. The mega-grains formalism allows us to separate the absorbed fractions in the clumps and in the ICM. It thus provides different heating rates in the two phases. In Fig. 3.5, we can clearly see that the clump emission (red) is significantly colder than that of the ICM (magenta). It however does not allow us to estimate the gradient of radiation field within the ICM and within clumps. This is the most dramatic in the case of central illumination. In order to obtain a more realistic SED for this particular case, VD99 used an ad hoc prescription, assuming a power-law distribution of equilibrium grain temperatures controlled by several tuning parameters depending on the grain type.

3.1.1.4 Rigorous Solutions

The radiative transfer equation (Eq. $3.16$) can be solved numerically. There are two main classes of methods (Steinacker et al., 2013, for a review).

Monte-Carlo Radiative Transfer

(MCRT) consists in simulating the random walk of photons, from their sources (stars, AGNs or thermal emission from dust grains) to the outside of the region, through their multiple scatterings.

Ray-tracing

numerically solves the radiative transfer equation on a discretized grid along multiple sightlines. It is more difficult to implement than MCRT, and can be computationally more intensive. It however allows the user a better assessment of the numerical errors.

MCRT is by far the most popular method. In this section, we briefly review its principle and apply it to an example.

Setting the model. To solve the radiative transfer equation (Eq. $3.16$), we need to specify the following physical ingredients.

A 3D spatial grid of dust density

has to be defined. Different coordinate systems can be chosen. In case there are steep density gradients (e.g. clumps), one has to make sure the transitions are finely enough sampled.

A 3D distribution of primary emitters,

as well as their SED needs to be defined. In this section, we will consider only stars, but AGNs can be treated the same way. It is possible to account both for: (i) discrete emission, such as an individual star or a cluster; and (ii) diffuse emission, such as unresolved stellar distribution.

The dust properties

need to include, at least: the opacity, the albedo and the asymmetry parameter. It is possible to account for a mixture of dust grains, with different compositions and sizes. The size-distribution-integrated quantities of Eqs. ($2.19$-$2.22$) are used for the transfer of stellar photons. It is however necessary to treat the thermal emission of individual species, especially if they are stochastically heated.

The principle of Monte-Carlo radiative transfer. To compute a MCRT model, we need to draw a large number of photons (typically $\simeq 10^6$ per wavelength bin), and execute the following steps. The procedure is schematically represented on Fig. 3.6.

1.

At each wavelength, photons are randomly drawn from primary emitters, proportionally to their specific intensity. The emission angle, $\left(𝜃,\varphi \right)$, is randomly chosen, if the emitters are isotropic, which is the case for stars.

2.

The interaction probability of a photon with a dust grain is then drawn from $1-\exp <!--FUNCTION\; APPLICATION-->\left(-\tau \right)$, along the original direction of the photon. On interaction, one accounts for the probabilities of: (i) absorption (proportional to $1-\tilde{\omega }$); and (ii) scattering (proportional to $\tilde{\omega }$). We could randomly choose between absorption and scattering, and thus terminate the path of the photons $1-\tilde{\omega }$ of the times. However, forced scattering, the way it is shown in Fig. 3.6, allows us to track the absorption and scattering probabilities at each interaction, in a numerically more efficient way. The scattering angle is randomly drawn from the scattering phase function (Eq. $1.43$), setting the new direction of the photon.

3.

We iterate this process as long as needed, until the photon exits the nebula. The average number of interactions increases with the effective optical depth of the nebula. At each interaction, $\left(1-\tilde{\omega }\right)\times I_{\nu }$ is absorbed by the grain, and $\tilde{\omega }\times I_{\nu }$ is scattered. After $N$ interactions, the scattered intensity is ${\tilde{\omega }}^N\times I_{\nu }$ and the absorbed power is proportional to $\left(1-\tilde{\omega }\right)\times {\tilde{\omega }}^{N-1}$. For a typical $\tilde{\omega }\simeq 0.5$, after three scatterings, the intensity is reduced by a factor $\simeq 0.13$.

4.

Once all the photons at all wavelengths have been drawn, the thermal emission of all dust species within each cell can be computed. IR photons are then drawn and scattered through the nebula, the same way as stellar photons. The overall opacity and albedo are usually much lower in the IR. The computation of IR radiative transfer is thus usually much faster. In the most embedded regions, the IR radiation absorbed by dust grains can be significant. One therefore has to recompute the IR transfer a few times, until an energy balance is reached.

These steps constitute the most basic implementation of MCRT. Numerous optimizations have however been proposed in the last fifty years (e.g. Witt, 1977a,b,c; Witt & Oshel, 1977; Yusef-Zadeh et al., 1984; Whitney & Hartmann, 1992; Wood, 1997; Wood & Jones, 1997; Városi & Dwek, 1999; Baes & Dejonghe, 2001; Gordon et al., 2001; Misselt et al., 2001; Steinacker et al., 2002, 2006; Wood et al., 2008; Baes et al., 2011; Camps & Baes, 2015; Siebenmorgen et al., 2015; Natale et al., 2015; Juvela, 2019). Improvements and optimizations include: (i) massive parallelization, and the use of GPU; (ii) the production of synthetic photometric images; (iii) the treatment of polarization by scattering.

Numerical method to randomly draw photons. To simulate the random walk of a photon, there are two sets of random variables to draw: (i)  random interaction events; and (ii) scattering angles.

Drawing interaction events

along the path of a photon consists in drawing a path length from the following Probability Density Function (PDF):

$$p\left(l\right)=1-\exp <!--FUNCTION\; APPLICATION-->\left(-{\int \; <!--nolimits-->}_0^l\kappa \rho \left(l^{\prime }\right)l^{\prime }d{l}^{\prime }\right).$$

(3.34)

For simplicity, $l$ corresponds to the path length along the direction of the photon, starting at $l=0$. Drawing a variable from this distribution can be achieved the following way (e.g. Városi & Dwek, 1999). It uses the rejection method (cf. Appendix C.2.3.1). These steps are represented on Fig. 3.7.

1.

First draw a uniform random variable between 0 and 1, ${\Theta }_1$.

2.

Set $l_0=-\ln <!--FUNCTION\; APPLICATION-->\left(1-{\Theta }_1\right)∕\kappa {\rho }_{max}$, where ${\rho }_{max}$ is the maximum of the density along the path.

3.

Draw a second uniform random variable between 0 and 1, ${\Theta }_2$.

4.

If ${\Theta }_2\le \rho \left(l_0\right)∕{\rho }_{max}$, then the interaction is accepted, and a new scattering angle can be drawn. On the contrary, if ${\Theta }_2>\rho \left(l_0\right)∕{\rho }_{max}$, the interaction is rejected. One then needs to go back to the first step, starting from $l_0$, this time.

Scattering angles

are drawn from the scattering phase function (Eq. $1.43$; Fig. 3.8.a). If we place a spherical coordinate reference frame at the position of interaction with its $z$-axis aligned with the incoming direction of the photon, then the two spherical angle, $\left(𝜃,\varphi \right)$, that will determine the new direction of the photon are drawn the following way.

1.

The new polar angle, $𝜃$, is drawn by inverting the Cumulative Distribution Function (CDF) (cf. Appendix C.2.3.2). In case we use the Henyey & Greenstein (1941) phase function, the inverse of the CDF can be analytically derived (cf. Fig. 3.8.b):

$$F^{-1}\left({\Theta }_3\right)=\frac{1}{2g}\left[1+g^2-{\left(\frac{1-g^2}{1-g+2g{\Theta }_3}\right)}^2\right],$$

(3.35)

where ${\Theta }_3$ is a uniform random variable between 0 and 1.

2.

The new azimuthal angle, $\varphi$, is simply drawn from a uniform distribution between 0 and $2\pi$, because there is a symmetry of revolution around the direction of the photon.

Example of a clumpy medium. For this manuscript, we have developed a MCRT model, following the procedure previously described. We have applied it to a spherical clumpy cloud, similar to those discussed in Sect. 3.1.1.3. The radius is $R_s=1$ pc, with an ICM density, $n_{ICM}=1000$ cm${}^{-3}$. The $r_C=0.05$ pc clumps have a density, $n_C=10^5$ cm${}^{-3}$, with a filling factor, $f_V=20\%$. We assume THEMIS grain constitution, and assume all grains are at thermal equilibrium. This cloud is centrally illuminated by a $T_{\star }=1.5\times 10^4$ K star, with $L_{\star }=3.3\times 10^4{L}_{\odot }$.

The projected map

of the dust mass surface density is shown in Fig. 3.9.a. The average absorbed fraction, $⟨P_{abs}⟩$ is shown in Fig. 3.9.b. The latter quantity is the projected average of the absorbed fraction of photons passing through the cell. This fraction is highly concentrated in the central region, as most of the power is absorbed by clumps close to the star. Those clumps are hot, whereas clumps on the outskirt are colder, because they are essentially heated by fewer photons that have been multiply scattered. Fig. 3.10 represents a few photon paths at three different wavelengths.

The total SED

of this cloud is shown in Fig. 3.11. We notice that a fraction of the dust emission is self-absorbed. This is mainly the 10 and 18 $\mu m$ silicate emission of the hottest central clumps being absorbed by the outer ones. We have also compared our model to the mega-grains approximation (Sect. 3.1.1.3; in dashed blue). We see this approximation is quite good in the UV-to-NIR range. It is however discrepant in the IR. This is because there is a very strong gradient of heating conditions, seen in Fig. 3.9.b, whereas the mega-grains approximation accounts only for the total absorbed fraction in the clumps and the ICM ².

Despite its intensive numerical requirements, MCRT is thus the most flexible and accurate way to compute the SED of an interstellar cloud. It is however important to make sure that the spatial resolution of the density grid is fine enough to resolve the shortest mean free path of photons. Otherwise, we would be smearing out a potential sub-grid temperature gradient. For the present model, we have 0.01 pc cells for $l_{mean}\left(U\right)\simeq 0.005$ pc.

3.1.2 Approximate Treatments of the Mixing of Physical Conditions

Ideally, each time we study a Galactic region or a galaxy, we should solve the radiative transfer equation (Eq. $3.16$). This is however, most of the time, impossible, because of the lack of constraints on the actual 3D structure of the region. Even if we have a collection of high-angular resolution multiwavelength images of our object, the matter and stellar distributions along the sightline have to be inferred. This inference is possible when the large-scale geometry of the object is quite regular, for instance: a protostellar disk and its jet, or a galactic disk and its bulge, etc. We will discuss the MCRT modeling of disk galaxies in Sect. 3.1.3.2. Otherwise, most often, we need to adopt empirical approaches that allows us to constrain the dust properties, despite our uncertainty of the spatial structure of the region.

3.1.2.1 The Historical Model: the MBB

The MBB (cf. Sect. 1.2.4.1) is historically the most widely-used dust SED model. It is controlled by the three following parameters (Eq. $1.70$ and Eq. $1.72$).

1.

The dust mass, $M_{dust}$, is a scaling parameter.

2.

The equilibrium temperature, $T$, controls the emission peak wavelength.

3.

The emissivity index, $\beta$ (Eq. $1.72$), controls the long-wavelength slope of the SED. This is demonstrated in Fig. 3.12.a.

Its physical assumptions are simplistic: (i) the IR emission is optically thin; (ii) the dust is made of a single species of grains at thermal equilibrium with the ISRF; and (iii) the opacity is a power-law. By inferring both $T$ and $\beta$, this model is designed to constrain both the dust excitation and its optical properties. This model was popularized by Hildebrand (1983), in the IRAS days (cf. Sect. 2.1.2). At the time, it was well adapted, being a simple, but still physical model, with only three parameters to fit four broadbands (the four IRAS bands at 12, 25, 60 and 100 $\mu m$). It has however several limitations that are often disregarded in the literature.

The mixing of physical conditions

is not accounted for by the MBB. It means that the SED fit of a complex region, containing a gradient of temperatures, is biased (e.g. Juvela & Ysard, 2012; Galliano et al., 2018). This is demonstrated in Fig. 3.12.b. The grey-filled curve is the SED coming from a power-law distribution of temperatures ($T=15$ to 60 K; the color curves). The intrinsic emissivity index of the grains making up this region is $\beta =1.79$. Yet, fitting such a SED with a single MBB leads to a compromise temperature ($T=45\pm 8$ K) and a systematically lower emissivity index, $\beta =1.12\pm 0.30$. This is because a gradient of temperature tends to flatten the FIR slope, similarly to a lower $\beta$ value.

$\Rightarrow$ The emissivity index derived from a single MBB fit is always lower than its true, intrinsic value.

Stochastic heating

is formally equivalent to a gradient of temperatures (cf. Sect. 1.2.4.3). Stochastically heated grains dominate the emission at MIR wavelengths (cf. Fig. 2.27). A single MBB is thus biased at short wavelengths by small grains, the same way it is biased at long wavelengths by cold dust. Solving this issue by fitting a linear combination of two or three MBBs can palliate this problem, but is usually not sufficient. In addition, fitting an out-of-equilibrium emission with equilibrium grains is physically incorrect. It renders the interpretation of the parameters of the hot MBB unreliable.

Realistic opacities are more complex

than a power-law. Laboratory data show that emissivity indices of actual materials are wavelength-dependent quantities, $\beta \left(\lambda \right)$ (cf. e.g. Fig. 2.22). Additionally, the somehow arbitrary choice of the two other parameters in Eq. ($1.72$) – the reference wavelength, ${\lambda }_0$, and the scaling of the opacity at this wavelength, ${\kappa }_0$ – have dramatic consequences on the derived dust mass. It is thus important to limit the potential range of variations of these parameters. The Kramers-Kronig relations (Sect. 1.2.1.6) impose that $\beta \ge 1$, and most compounds studied in the laboratory have $\beta <2.5$ (cf. e.g. the Jena database). We recommend calibrating ${\kappa }_0$ on laboratory data or on well-constrained dust models (such as Eq. $2.26$). Contrary to what Hildebrand (1983) recommended, it is probably better to choose a reference wavelength, ${\lambda }_0$, around the peak of the FIR SED (between 100 and 300 $\mu m$), rather than in the submm. This way, variations of $\beta$ will impact only moderately the derived dust mass.

$\Rightarrow$ A MBB fit infers parameters whose physical meaning is difficult to assess.

The noise-induced degeneracy

between $T$ and $\beta$ is well-documented (e.g. Shetty et al., 2009; Kelly et al., 2012; Galliano, 2018). A false negative correlation arises when a series of MBB fits are performed with standard least-squares, maximum likelihood, or non-hierarchical Bayesian methods. This degeneracy prevents to explore the potential correlation of these parameters that laboratory data and solid-state models suggest (e.g. Mennella et al., 1998; Meny et al., 2007). Luckily, hierarchical Bayesian methods are efficient at removing false, noise-induced correlations (e.g. Kelly et al., 2012; Galliano, 2018). Chap. 5 is entirely devoted to fitting techniques.

3.1.2.2 A Phenomenological, Composite Approach

Dust models provide useful frameworks to model SEDs (Sect. 2.3). Without the possibility to compute the radiative transfer, we are however facing the problem of the mixing of physical conditions. A prescription, proposed by Dale et al. (2001), has proven to be a powerful solution to this issue. It consists in assuming that the dust mass is distributed in regions with different starlight intensities, $U$, following a power-law:

The idea is that the shape of the observed SED is used to constrain this distribution of ISRFs, assuming a dust mixture constitution. By lack of a better term, we call this approach the composite model. The free parameters are:

• the total dust mass, $M_{dust}$, acting as a scaling parameter;
• the power-law index, $\alpha$;
• the minimum starlight intensity, $U_{-}$;
• the width of the starlight intensity distribution, $\Delta U$.

It thus provides a way to account for the potential complexity of the region without having to model the radiative transfer. The model SED is then simply: