© ESO 2019

1. Introduction

Radiative transfer is a broad field in astronomy that aims to describe the interaction between radiation and matter. In an astronomical context, the Monte Carlo method is by far the most widely used method to solve radiative transfer problems. In the past decades, many different Monte Carlo codes have been developed to address different astrophysical radiative transfer problems, including photoionisation, absorption and scattering by cosmic dust, the origin of infrared emission, and resonant line scattering heating (e.g. Gordon et al. 2001; Ercolano et al. 2003; Robitaille 2011; Yajima et al. 2012; Whitney et al. 2013; Camps & Baes 2015; Reissl et al. 2016). General reviews on Monte Carlo transport can be found in Dupree & Fraley (2002) or Kalos & Whitlock (2008), for example, and dedicated reviews on radiative transfer in astrophysics include Whitney (2011) and Steinacker et al. (2013).

The essence of the Monte Carlo method is to represent the radiation field as the flow of a large, but finite, number of photon packages. The life cycle of each photon package is followed individually. Additionally, at every stage in this life cycle, the characteristics that determine the path of each photon package are determined in a probabilistic way by generating random numbers from the appropriate probability density function (PDF). At the end of the simulation, the radiation field, or more specifically the intensity of the radiation field, is recovered from a statistical analysis of the photon package paths.

An important ingredient of Monte Carlo radiative transfer is the knowledge of the appropriate PDF for a given characteristic, and the accurate and efficient sampling of random numbers from these PDFs. For some characteristics, the PDFs are simple and sampling random numbers from them is trivial. For example, most sources send out radiation isotropically, which implies that the generation of propagation directions after an emission event simply comes down to generating a random point on the unit sphere. The PDF that controls the random starting positions for the photon package is dictated by the 3D luminosity density of the sources, and specific methods have been developed to generate such random positions from a range of 3D density distributions (Baes & Camps 2015).

An aspect that is central to any Monte Carlo radiative transfer code is the random generation of the next interaction location. More specifically, if a photon package is emitted or scattered into a given direction, one needs to randomly generate a free path length s to the next interaction. To do so, we need to know the appropriate PDF p(s). In this context, the concept of optical depth τ plays a crucial role. In optical depth space, the PDF p(τ) is a simple exponential distribution (Cashwell & Everett 1959; Steinacker et al. 2013). This implies that the next interaction location can be found by randomly generating a random optical depth from an exponential distribution and by converting this optical depth to a physical path length, which is a critical point. We need to know the relation τ(s), or inversely s(τ), to properly calculate the next interaction location.

In radiative transfer problems where polarisation is not taken into account, τ(s) can immediately be calculated from the density and the optical properties along the path. Notably, it does not depend on the intensity of the radiation field. The situation becomes more complex when polarisation comes into play. In the case of spherical or randomly oriented particles, polarised Monte Carlo radiative transfer is only slightly more difficult to use than unpolarised Monte Carlo radiative transfer. The main added complexities are that a Stokes vector needs to be introduced to characterise the polarisation status of each photon package and that this Stokes vector can alter during a scattering event (e.g. Fischer et al. 1994; Code & Whitney 1995; Bianchi et al. 1996; Peest et al. 2017). The relation between optical depth and path length is the same as for unpolarised Monte Carlo radiative transfer.

The real complexity arises when the attenuating particles are non-spherical and aligned, as in the case of elongated dust grains in the interstellar medium. Such dust grains are aligned with respect to the magnetic fields through a variety of processes (Davis & Greenstein 1951; Jones & Spitzer 1967; Aannestad & Greenberg 1983; Lazarian 1994, 2007). An interesting feature is dichroism, which, in this context, means that radiation of different polarisation experience different amounts of extinction. The nature of dichroism complicates the relation between optical depth and physical path length.

In this paper, we investigate in detail whether an optical depth scale can still be used in a meaningful way in polarised Monte Carlo radiative transfer codes that take dichroic extinction into account. In Sect. 2 we discuss optical depth and the generation of random path lengths in standard unpolarised radiative transfer. In Sect. 3 we extend this discussion to polarised radiative transfer in the case of elongated aligned grains, in particular. Apart from the “standard” extinction optical depth, we introduce the total and dichroic optical depth scales, and we compare and analyse them. In Sect. 4 we discuss the implications of these findings for Monte Carlo radiative transfer codes that take dichroic extinction into account. In Sect. 5 we discuss these results and provide a summary.

2. Unpolarised radiative transfer

As discussed in the Introduction, one of the essential steps in the life cycle of a photon package in a Monte Carlo radiative transfer simulation is the calculation of the next interaction location, or equivalently, the physical path length s that is covered before the next interaction. In order to do this, we need to generate a random s from the appropriate probability distribution function p(s). The appropriate PDF p(s) can be determined by considering the variation of the specific intensity I(s) along the path. The probability that the photon package has not interacted along the path between 0 and s is equal to I(s)/I₀, while I₀ is the specific intensity of the photon package at the start of the path¹. Therefore, the density function P(s) can be written as follows:

$$\begin{array}{c}\hfill P(s)={{\displaystyle \int }}_0^{s}p(s\prime )\mathrm{d}s\prime =1-\frac{I(s)}{I_0}·\end{array}$$(1)

In general, we define the optical depth τ(s) as

$$\begin{array}{c}\hfill \tau (s)=-ln\left[\frac{I(s)}{I_0}\right]·\end{array}$$(2)

By combining these two equations, we find that

$$\begin{array}{c}\hfill P(s)=1-{\mathrm{e}}^{-\tau (s)},\end{array}$$(3)

and when we take the derivative of this cumulative density function,

$$\begin{array}{c}\hfill p(\tau )={\mathrm{e}}^{-\tau }.\end{array}$$(4)

This yields the well-know result that the PDF describing the next interaction location is an exponential distribution in optical depth space. Generating a random optical depth can easily be done by picking a uniform deviate 𝒳, setting τ = −ln(1 − 𝒳), and subsequently converting this random τ to a physical path length s.

The difficulty involved is that we need to know the solution I(s) to calculate the optical depth scale. This solution is found by solving the appropriate radiative transfer equation. For a photon package moving through an attenuating medium with density n(s) and extinction cross section C_ext(s), the radiative transfer equation reads

$$\begin{array}{c}\hfill \frac{\mathrm{d}I}{\mathrm{d}s}(s)=-n(s)C_{\mathrm{ext}}(s)I(s).\end{array}$$(5)

It is important to note that no additional emission along the path or scattering into the line-of-sight are included as this is not relevant for this particular photon package. The following equation is a simple first-order ordinary differential equation that is easy to solve:

$$\begin{array}{c}\hfill I(s)=I_0{\mathrm{e}}^{-{\tau }_{\mathrm{ext}}(s)}\end{array}$$(6)

with the extinction optical depth τ_ext(s) that is defined as

$$\begin{array}{c}\hfill {\tau }_{\mathrm{ext}}(s)={{\displaystyle \int }}_0^{s}n(s^{\prime })C_{\mathrm{ext}}(s^{\prime })\mathrm{d}{s}^{\prime }.\end{array}$$(7)

By comparing Eqs. (2) and (6), we see that τ(s) = τ_ext(s). In particular, the optical depth scale, from which random path lengths can be sampled, only depends on the density and optical properties of the material, and not on the properties of the photon package. In this simple case, the only property of the photon package that could matter is I₀.

While the strategy described above is conceptually very simple, some challenges need to be addressed. In particular, the conversion of extinction optical depth to physical path length is not usually a straightforward inversion. Except for some simple ideal cases, the conversion needs to be done numerically. In most Monte Carlo radiative transfer codes, the attenuating medium is subdivided into a large number of individual cells, each with a constant density and uniform properties. The codes are typically equipped with a routine to calculate paths through this tessellated medium. The routine returns an ordered list of all the cells m that the path crosses, as well as the length Δs_m of the path segments corresponding to the mth cell. Given this ordered list, we can calculate the running values for the path length s_m and the optical depth τ_ext, m at the exit point of each cell crossed by the path. The problem is thus reduced to finding the first cell in the array, for which τ_ext, m exceeds the randomly generated value of τ_ext, and subsequently to applying a linear interpolation to convert τ_ext to s².

These integrations through the dust grid often form the most time-consuming part of a radiative transfer simulation. In order to make these calculations as efficient as possible, especially in 3D geometries, advanced grid construction and grid traversal techniques are required (Niccolini & Alcolea 2006; Bianchi 2008; Lunttila & Juvela 2012; Camps et al. 2013; Saftly et al. 2013, 2014; Hubber et al. 2016).

3. Polarised radiative transfer

3.1. The Stokes formalism

The characterisation of the radiation field by the specific intensity and the corresponding radiative transfer Eq. (5) are no longer suitable when polarisation is considered. In order to take the polarisation state of radiation into account, one can use the Stokes formalism, which characterises the radiation field by means of the 4D Stokes vector

$$\begin{array}{c}\hfill \mathit{I}=\left(\begin{array}{c}I\\ Q\\ U\\ V\end{array}\right).\end{array}$$(8)

The first Stokes parameter, I, is still the specific intensity. The Stokes parameters Q and U describe the state of linear polarisation and V describes the state of circular polarisation of the radiation. The Stokes parameters are always defined with respect to a reference direction that is to be chosen freely from the plane perpendicular to the propagation direction. For a detailed description of the Stokes vector and its connection to the monochromatic transverse electromagnetic waves, we refer to Mishchenko et al. (2000).

When we consider the full Stokes vector, the simple radiative transfer Eq. (5) becomes

$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathit{I}}{\mathrm{d}s}(s)=-n(s)\mathbf{K}(s)\mathit{I}(s),\end{array}$$(9)

where K is now the extinction matrix, a 4 × 4 matrix that describes how the different Stokes components are affected when radiation passes through the medium.

3.2. Spherical grains

When dust grains are spherical, or non-spherical but arbitrarily oriented, the extinction matrix is a simple diagonal matrix, all components of the Stokes vector are affected in the same way, and

$$\begin{array}{c}\hfill \mathbf{K}=\left(\begin{array}{cccc}{C}_{\mathrm{ext}}& 0& 0& 0\\ 0& C_{\mathrm{ext}}& 0& 0\\ 0& 0& C_{\mathrm{ext}}& 0\\ 0& 0& 0& C_{\mathrm{ext}}\end{array}\right).\end{array}$$(10)

This implies that there is no mixture of the different Stokes components due to extinction, and the solution of the radiative transfer equation can be directly written as

$$\begin{array}{c}\hfill \mathit{I}(s)={\mathit{I}}_0{\mathrm{e}}^{-{\tau }_{\mathrm{ext}}(s)}.\end{array}$$(11)

In particular, the specific intensity I(s) still behaves according to Eq. (6), exactly as is the case for non-polarised radiative transfer. Also, in this case, we can immediately conclude that the PDF describing the net interaction location is an exponential distribution in extinction optical depth space.

3.3. Spheroidal grains

Complexity arises when the dust grains are non-spherical and partially aligned. In this case, the extinction matrix K is not a diagonal matrix, but a full 4 × 4 matrix with 16 nonzero cross sections of which only seven are independent (van de Hulst 1957; Hovenier & van der Mee 1996). Fortunately, in the case of spheroidal grains, K is significantly less complex. If we denote the orientation of the grain alignment at a distance s along the path as ψ(s), K can be expressed as (Martin 1974; Wolf et al. 2002)

$$\begin{array}{c}\hfill \mathbf{K}(s)=\mathcal{L}(-\psi (s)){\mathbf{K}}_{\mathrm{ref}}(s)\mathcal{L}(\psi (s))\end{array}$$(12)

with ℒ(ψ) as the Mueller rotation matrix,

$$\begin{array}{c}\hfill \mathcal{L}(\psi )=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& cos2\psi & sin2\psi & 0\\ 0& -sin2\psi & cos2\psi & 0\\ 0& 0& 0& 1\end{array}\right),\end{array}$$(13)

and K_ref as the extinction matrix in the reference frame of the grain,

$$\begin{array}{c}\hfill {\mathbf{K}}_{\mathrm{ref}}=\left(\begin{array}{cccc}{C}_{\mathrm{ext}}& C_{\mathrm{pol}}& 0& 0\\ C_{\mathrm{pol}}& C_{\mathrm{ext}}& 0& 0\\ 0& 0& C_{\mathrm{ext}}& C_{\mathrm{cpol}}\\ 0& 0& -C_{\mathrm{cpol}}& C_{\mathrm{ext}}\end{array}\right).\end{array}$$(14)

We note that K_ref has only the following three independent elements: the extinction cross section C_ext, polarisation cross section C_pol, and circular polarisation cross section C_cpol (Mishchenko et al. 2000; Whitney & Wolff 2002).

The non-diagonal character of the extinction matrix has an important effect on the radiation field; the different components of the Stokes vector are coupled and will get mixed along the path. In particular, radiation that is initially unpolarised can develop linear and even circular polarisation just by propagating through the medium (Serkowski 1962; Martin 1972, 1974).

This dichroism complicates the relation between path length and optical depth. To stress the fact that we deal with dichroic extinction, we use the term dichroic optical depth and use the notation τ_dic(s) when we consider the optical depth in a dichroic medium. In general we can write that

$$\begin{array}{c}\hfill {\tau }_{\mathrm{dic}}(s)=-ln\left[\frac{I(s)}{I_0}\right],\end{array}$$(15)

where I(s) should now be seen as the first component of the Stokes vector I(s), which is obtained by solving the vector radiative transfer Eq. (9) with K(s), as given by (12).

In the case of spherical grains, we show that τ(s) = τ_ext(s). However, this equivalence does not hold for the general case of spheroidal grains. Indeed, with an extinction matrix given by expression (12), the radiative transfer equation is a set of four coupled first-order ordinary differential equations with varying coefficients. The full solution for I(s), and hence the dichroic optical depth τ_dic(s), depends on all elements of the extinction matrix and all initial Stokes components. As the extinction optical depth (7) only depends on one single element (K₁₁ = C_ext) of the extinction matrix, and as it is independent of the initial polarisation state I₀, it is impossible that the dichroic optical depth is equivalent to the extinction optical depth.

One could consider another optical depth scale that does depend on the initial polarisation state of the photon package. An interesting starting point is the total extinction cross section,

$$\begin{array}{c}\hfill {\stackrel{\sim }{C}}_{\mathrm{ext}}=K_{11}+\frac{Q_0}{I_0}{K}_{12}+\frac{U_0}{I_0}{K}_{13}+\frac{V_0}{I_0}{K}_{14},\end{array}$$(16)

introduced by Mishchenko et al. (2000) and Whitney & Wolff (2002). It properly describes the attenuation of an incident beam with the arbitrary initial polarisation state I₀. Based on this total extinction cross section, one can define an optical depth scale, which we denote as the total optical depth, in a similar way as for the extinction optical depth (7),

$$\begin{array}{c}\hfill {\tau }_{\mathrm{tot}}(s)={{\displaystyle \int }}_0^{s}n(s^{\prime }){\stackrel{\sim }{C}}_{\mathrm{ext}}(s^{\prime })\mathrm{d}{s}^{\prime },\end{array}$$(17)

or explicitly,

$$\begin{array}{cc}\hfill {\tau }_{\mathrm{tot}}(s)& ={\tau }_{\mathrm{ext}}(s)+\frac{Q_0}{I_0}{{\displaystyle \int }}_0^{s}n(s^{\prime })C_{\mathrm{pol}}(s^{\prime })cos2\psi (s^{\prime })\mathrm{d}{s}^{\prime }\hfill \\ \hfill & +\frac{U_0}{I_0}{{\displaystyle \int }}_0^{s}n(s^{\prime })C_{\mathrm{pol}}(s^{\prime })sin2\psi (s^{\prime })\mathrm{d}{s}^{\prime }.\hfill \end{array}$$(18)

Despite the increased complexity, it is computationally not much more difficult to employ than the extinction optical depth. Indeed, for a photon package in a Monte Carlo simulation with a given position, propagation direction, and polarisation state, it is relatively easy to compute the total optical depth τ_tot, m at the entry point of each cell crossed by the path.

We can, however, use analogous reasoning as above to argue that the total optical depth cannot be equivalent to the dichroic optical depth. Indeed, with an extinction matrix given by expression (12), the radiative transfer equation is a set of four coupled first-order ordinary differential equations with varying coefficients. The full solution for I(s), and hence the dichroic optical depth τ_dic(s), depends on all elements of the extinction matrix and all initial Stokes components. The expression (18) does indeed depend on the initial linear polarisation state of the incoming radiation and the polarisation cross section, but it does not involve any dependency on the circular polarisation cross section C_cpol or the initial circular polarisation V₀. As a result, it is impossible that the total extinction is fully equivalent to the dichroic optical depth.

If neither the extinction optical depth nor the total optical depth are equivalent to the dichroic optical depth, one inevitably questions what the correct expression for the dichroic optical depth is. For a medium with non-uniform grain alignment, the following answer is somewhat disappointing: given the non-trivial coupling of all four components of the radiative transfer equation, it is not possible to write down a closed formula for I(s) nor the dichroic optical depth τ_dic(s). Obtaining the solution I(s) can be achieved numerically using standard vector ordinary differential equation (ODE) solution techniques. One might wonder whether τ_ext or τ_tot, for which we have a closed expression, are reasonable approximations for the dichroic optical depth. To answer this question, we performed numerical tests.

We adopted a hypothetical example where n(s) = ψ(s) = s, that is, the density of material increases linearly with increasing distance, and the grain alignment rotates around the path. Furthermore, we assumed that the optical properties do not vary along the path, and we set C_ext = 1 and C_pol = C_cpol = 0.1 (representative values of the C_pol/C_ext ratio go up to 0.3 and more). On the one hand, we calculated the dichroic optical depth τ(s) by numerically solving the vector radiative transfer equation using an explicit Runge–Kutta method with variable step size control. On the other hand, we used Eqs. (7) and (17) to calculate τ_ext(s) and τ_tot(s), and we calculated the relative differences

$$\begin{array}{c}\hfill {\delta }_{\mathrm{ext}}=\frac{{\tau }_{\mathrm{ext}}-{\tau }_{\mathrm{dic}}}{{\tau }_{\mathrm{dic}}}\mathrm{and}{\delta }_{\mathrm{tot}}=\frac{{\tau }_{\mathrm{tot}}-{\tau }_{\mathrm{dic}}}{{\tau }_{\mathrm{dic}}}·\end{array}$$(19)

We calculated these quantities for initial Stokes vectors, that are 100% linearly polarised, with the linear polarisation angle

$$\begin{array}{c}\hfill \theta =\frac{1}{2}arctan\frac{U_0}{Q_0}=\frac{1}{2}arccos\frac{Q_0}{I_0}\end{array}$$(20)

by gradually changing between 0° and 90° in steps of 5°.

In the upper left panel of Fig. 1 we show δ_ext as a function of the dichroic optical depth τ_dic along the path. The extinction optical depth scale has the advantage that it is simple and does not explicitly depend on the polarisation state of the photon package, such that the relation between s and τ_ext can in principle be precomputed. It is clear, however, that τ_ext is a poor approximation for the dichroic optical depth. Depending on the value of Q₀/I₀, it can either underestimate or overestimate τ_dic by up to 20 percent. More importantly, the differences between τ_dic and τ_ext can be large even at small optical depths.

Fig. 1. Top left: relative differences δext between extinction optical depth τext and dichroic optical depth τdic for hypothetical example with density and grain alignment increasing linearly with increasing path length. The different curves correspond to different input Stokes vectors. In each case the photon package is initially fully linearly polarised, but the Stokes vector orientation θ increases in steps of 5°. The red lines corresponds to θ = 0° (Q0/I0 = 1) and the purple line to θ = 90° (Q0/I0 = −1), respectively. Bottom left: similar to top left panel, but now showing relative difference δtot between total optical depth τtot and dichroic optical depth. Right: same as panels on left, but now for model with uniform grain alignment.

The relative difference δ_tot in the bottom left panel shows very different behaviour. The total optical depth always overestimates the dichroic optical depth. A second important difference is that δ_tot is usually smaller than the absolute value of δ_ext. This is particularly true in the optically thin limit where τ_tot approximates τ_dic very well. At large optical depths (τ_dic >  10), it turns out that the total optical depth is not always a reliable approximation for the dichroic optical depth. It is interesting to note that τ_tot can even become a poorer approximation for the dichroic optical depth than the simpler approximation τ_ext.

3.4. Spheroidal grains with uniform alignment

We can gain more insight into these results by considering the special case ψ(s) = 0 in which the grains are all uniformly aligned along the path³. In this special case, the rotation matrices in (12) are the identity matrices, and K = K_ref. With this relatively simple block-diagonal extinction matrix, the four different components of the Stokes vector are paired, instead of fully coupled as follows: I and Q in addition to U and V are linked, respectively. The full solution of the radiative transfer equation can be written out as

$$I(s)={\text{e}}^{-{\tau }_{\text{ext}}(s)}\left[I_0\cosh {\tau }_{\text{pol}}(s)-Q_0\sinh {\tau }_{\text{pol}}(s)\right],$$(21a)

$$Q(s)={\text{e}}^{-{\tau }_{\text{ext}}(s)}\left[Q_0\cosh {\tau }_{\text{pol}}(s)-I_0\sinh {\tau }_{\text{pol}}(s)\right],$$(21b)

$$U(s)={\text{e}}^{-{\tau }_{\text{ext}}(s)}\left[U_0\cos {\tau }_{\text{cpol}}(s)-V_0\sin {\tau }_{\text{cpol}}(s)\right],\text{\hspace{1em}and}$$(21c)

$$V(s)={\text{e}}^{-{\tau }_{\text{ext}}(s)}\left[V_0\cos {\tau }_{\text{cpol}}(s)+U_0\sin {\tau }_{\text{cpol}}(s)\right],$$(21d)

with

$$\begin{array}{cc}& {\tau }_{\mathrm{pol}}(s)={{\displaystyle \int }}_0^{s}n(s^{\prime })C_{\mathrm{pol}}(s^{\prime })\mathrm{d}{s}^{\prime }\mathrm{and}\hfill \end{array}$$(22)

$$\begin{array}{cc}& {\tau }_{\mathrm{cpol}}(s)={{\displaystyle \int }}_0^{s}n(s^{\prime })C_{\mathrm{cpol}}(s^{\prime })\mathrm{d}{s}^{\prime }.\hfill \end{array}$$(23)

By combining Eq. (2) with the solution (21a) for the specific intensity, we find an explicit expression for the dichroic optical depth:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{dic}}(s)={\tau }_{\mathrm{ext}}(s)-ln[cosh{\tau }_{\mathrm{pol}}(s)-\frac{Q_0}{I_0}sinh{\tau }_{\mathrm{pol}}(s)]·\end{array}$$(24)

In the panels on the right-hand side of Fig. 1 we show δ_ext and δ_tot as a function of τ_dic along the path for the case with ψ = 0. Again, τ_ext is a poor approximation for the dichroic optical depth even at small optical depths. Based on the explicit expression (23), we can calculate the extreme values for δ_ext, which correspond to Q₀/I₀ = ±1,

$$\begin{array}{c}\hfill {\delta }_{\mathrm{ext}}^{\pm }=\mp \frac{C_{\mathrm{pol}}}{C_{\mathrm{ext}}\pm C_{\mathrm{pol}}}·\end{array}$$(25)

For our example, these differences run up to −17% and 25%. On the other hand, the total optical depth is a better approximation, in general, for the dichroic optical depth, especially in the optically thin limit τ ≪ 1. This can be understood by considering the Taylor expansion for expression (23) for τ_pol ≪ 1,

$$\begin{array}{c}\hfill {\tau }_{\mathrm{dic}}=\underset{{\tau }_{\mathrm{tot}}}{\underbrace{{\tau }_{\mathrm{ext}}+\frac{Q_0}{I_0}{\tau }_{\mathrm{pol}}}}-\frac{1}{2}[1-{\left(\frac{Q_0}{I_0}\right)}^2]{\tau }_{\mathrm{pol}}^2+\mathcal{O}({\tau }_{\mathrm{pol}}^3).\end{array}$$(26)

This result shows that the extinction optical depth approximates the dichroic optical depth to first order in τ_pol, whereas the total optical depth approximates it to second order. For small optical depths, it is hence not surprising that the dichroic optical depth is a much better approximation. Moreover, as the coefficient of the second order term is always negative, the total optical depth systematically overestimates the dichroic optical depth. It is important to note, however, that the total optical depth is not guaranteed to be a reliable estimator for the dichroic optical depth as at large optical depths, τ_tot overestimates τ_dic significantly.

4. Application in Monte Carlo radiative transfer

4.1. Random generation of path lengths

In the previous Section, we demonstrate that neither the simple extinction optical depth nor the total optical depth are equivalent to the dichroic optical depth. In particular, we see that, for a given physical path length, relative differences of up to several tens of percent are possible between these optical depth scales. It is interesting to note that for the total optical depth scale, these large differences only occur at large optical depths. As discussed in the Introduction, the optical depth in Monte Carlo radiative transfer simulations is mainly important for the random generation of the next interaction location. In practice, a random optical depth is generated from an exponential PDF, and this optical depth is translated to a physical path length, which immediately sets the next interaction location.

In order to find out to which degree the choice of the optical depth scale affects the determination of the path length, we performed a number of simple simulations. We adopted a similar setup as in Sect. 3.3, with n(s) = ψ(s) = s, C_ext = 1, and C_pol = C_cpol = 0.1, and we used a maximum path length s_max = 10, corresponding to τ_ext, max = 50. We generated one million random optical depths from an exponential PDF, and we converted these values to physical path lengths according to each of the three different optical depth scales.

The top row of Fig. 2 shows the histograms of the corresponding path lengths for three different input linear polarisations (Q₀/I₀ = 0, 0.5 and 1, from left to right). For initially unpolarised photon packages, the three distributions are very similar. In fact, the distributions corresponding the extinction and total optical depth are identical as τ_ext(s) ≡ τ_tot(s) for initially unpolarised radiation. When the initial linear polarisation degree increases, the histograms for the dichroic and total optical depth scales remain very similar, but they gradually deviate more from the histogram corresponding to the extinction optical depth. We applied two-sample Kolmogorov–Smirnov tests to quantify this observation. These tests demonstrate that there is indeed no significant difference between the path length distribution corresponding to the total optical depth scale and the dichroic optical depth scale, whereas the path length distribution corresponding to the extinction optical depth scale is significantly different in the case of initially linearly polarised photon packages. This result is not surprising, as the exponential PDF strongly favours small random optical depths, and in the optically thin regime, the total optical depth approximates the dichroic optical depth very well (see Fig. 1).

Fig. 2. Histograms for distribution of randomly generated path lengths, corresponding to different optical depth scales discussed in this paper. The model is described in Sect. 4.1, and the different columns correspond to different levels of initial linear polarisation. Top panels: optical depths randomly generated from an exponential distribution, bottom panels: optical depths generated by using the composite biasing approach discussed in Baes et al. (2016).

However, it is not always advisable to randomly generate optical depths from an exponential distribution. One particularly relevant case is the penetration of radiation through an optically thick medium, which is a notoriously difficult task for Monte Carlo radiative transfer (e.g. Min et al. 2009; Gordon et al. 2017; Camps & Baes 2018). For a standard exponential PDF, huge numbers of photon packages need to be generated before a single one might cross the barrier. Path length stretching (Levitt 1968; Spanier 1970; Dwivedi 1982) is a Monte Carlo optimisation technique that artificially stimulates the generation of larger path lengths. Baes et al. (2016) combine this technique with the concept of composite biasing, which results in an approach that has the advantages of path length stretching (an increased probability for large path lengths) and minimises the disadvantages (no large weight factors for the photon packages). The bottom row of Fig. 2 shows the same histograms as the top row, but now the random optical depths are sampled from a biased distribution (for details, see Baes et al. 2016, Sect. 4.3). As the total optical depth is no longer guaranteed to be a good approximation for the dichroic optical depth at high optical depths, it is no surprise that the distributions are now more different, especially in the high path length tail. Kolmogorov–Smirnov (KS) tests indicate that the probability that they are drawn from the same distribution is negligible.

4.2. Practical calculation of the dichroic optical depth

The results from the previous subsections show that it is not a good idea to use the extinction or total optical depth scale to generate random path lengths in a Monte Carlo radiative transfer simulation with non-spherical dust grains. To avoid systematic errors, the dichroic optical depth scale should be used to translate randomly generated optical depths to physical path lengths. We hence need an efficient calculation for the dichroic optical depth τ_dic(s) along the path. However, in Sect. 3.3 we argue that it is not possible to write down a general closed formula for τ_dic(s), which is similar to the “simple” formulae (7) and (18) for the extinction and total optical depth scales, respectively. In general, a solution for τ_dic(s) needs to be obtained numerically using standard vector ODE solution techniques.

Fortunately, in the context of Monte Carlo radiative transfer simulations, in which the density of the medium is usually discretised on a grid with a uniform density and uniform optical properties within each grid cell, there is an efficient routine to solve the radiative transfer equation without the need for vector ODE solution methods (see also Whitney & Wolff 2002; Lucas 2003). Indeed, we can progressively solve the radiative transfer equation in each individual cell along the path. We assume again that the path is split into individual cells that are small enough so the density, optical properties, and grain orientation can be considered uniform within each cell. We denote the Stokes vector at the entry point of cell m as I_m − 1. We first rotate the Stokes vector over an angle ψ_m to align it with the grain orientation within cell m,

$$\begin{array}{c}\hfill {\mathit{I}}_{m-1}^{\prime }=\mathcal{L}({\psi }_m){\mathit{I}}_{m-1}.\end{array}$$(27)

As the grain orientation with the cell is constant, we can directly apply the solution (21) for the radiative transfer to calculate the Stokes vector at the exit point s = s_m of the cell,

$$\begin{array}{cc}& I_m^{\prime }={\mathrm{e}}^{-\mathrm{\Delta }{\tau }_{\mathrm{ext},m}}[I_{m-1}^{\prime }cosh\mathrm{\Delta }{\tau }_{\mathrm{pol},m}-Q_{m-1}^{\prime }sinh\mathrm{\Delta }{\tau }_{\mathrm{pol},m}],\hfill \end{array}$$(28)

$$\begin{array}{cc}& Q_m^{\prime }={\mathrm{e}}^{-\mathrm{\Delta }{\tau }_{\mathrm{ext},m}}[Q_{m-1}^{\prime }cosh\mathrm{\Delta }{\tau }_{\mathrm{pol},m}-I_{m-1}^{\prime }sinh\mathrm{\Delta }{\tau }_{\mathrm{pol},m}],\hfill \end{array}$$(29)

$$\begin{array}{cc}& U_m^{\prime }={\mathrm{e}}^{-\mathrm{\Delta }{\tau }_{\mathrm{ext},m}}[U_{m-1}^{\prime }cos\mathrm{\Delta }{\tau }_{\mathrm{cpol},m}-V_{m-1}^{\prime }sin\mathrm{\Delta }{\tau }_{\mathrm{cpol},m}],\hfill \end{array}$$(30)

$$\begin{array}{cc}& V_m^{\prime }={\mathrm{e}}^{-\mathrm{\Delta }{\tau }_{\mathrm{ext},m}}[V_{m-1}^{\prime }cos\mathrm{\Delta }{\tau }_{\mathrm{cpol},m}+U_{m-1}^{\prime }sin\mathrm{\Delta }{\tau }_{\mathrm{cpol},m}],\hfill \end{array}$$(31)

and

$$\begin{array}{c}\hfill \mathrm{\Delta }{\tau }_{\mathrm{xxx},m}=n_m{C}_{\mathrm{xxx},m}\mathrm{\Delta }{s}_m\end{array}$$(32)

When we rotate the resulting Stokes vector back to the laboratory frame, we find the desired result I_m,

$$\begin{array}{c}\hfill {\mathit{I}}_m=\mathcal{L}(-{\psi }_m){\mathit{I}}_m^{\prime }.\end{array}$$(33)

This recipe can be repeated for all cells along the path.

We tested this strategy using the same example as already discussed. The comparison between the brute-force Runge–Kutta approach and the algorithm described above is shown in Fig. 3 for an initially unpolarised Stokes vector. The top panel shows the evolution of the individual Stokes components, the bottom panels shows the Stokes ratios, as well as the degree of linear and total polarisation

$$\begin{array}{c}\hfill p_{\mathrm{L}}=\frac{\sqrt{Q^2+U^2}}{I},\mathrm{and}p=\frac{\sqrt{Q^2+U^2+V^2}}{I}.\end{array}$$(34)

Fig. 3. Evolution of initially unpolarised Stokes vector propagating through medium of aligned grains where grain alignment rotates around path. The dots are the result from a Runge–Kutta numerical integration of the vector radiative transfer equation and the solid lines correspond to the method outlined in Sect. 3.

The two methods clearly agree. As the initially unpolarised photon propagates along the path, it gradually develops linear polarisation as a result of dichroic extinction and also circular polarisation later on. From τ ≳ 8, the circular polarisation starts to dominate the linear polarisation, and at τ = 50, the photon is almost 100% polarised.

The determination of a random path length now follows the same strategy as discussed for the unpolarised Monte Carlo radiative transfer in Sect. 2. From the solution of the Stokes vector I_m, we calculated the dichroic optical depth τ_dic, m at the exit point of each cell, and we searched for the first cell for which τ_dic, m exceeds the randomly determined τ. One additional difference needs to be taken into account. In the case of unpolarised radiation transfer, the increase in optical depth within each cell is directly proportional to the increase in path length within that cell, for example

$$\begin{array}{c}\hfill {\tau }_{\mathrm{ext}}(s)={\tau }_{\mathrm{ext},m-1}+n_m{C}_{\mathrm{ext},m}(s-s_{m-1})s_{m-1}\le s\le s_m.\end{array}$$(35)

To find the exact path length s corresponding to a randomly generated τ, we can therefore use simple linear interpolation as follows:

$$\begin{array}{c}\hfill s-s_{m-1}=\frac{\tau -{\tau }_{\mathrm{ext},m-1}}{n_m{C}_{\mathrm{ext},m}}=\left(\frac{s_m-s_{m-1}}{{\tau }_{\mathrm{ext},m}-{\tau }_{\mathrm{ext},m-1}}\right)(\tau -{\tau }_{\mathrm{ext},m-1}).\end{array}$$(36)

In the case of dichroic attenuation, the increase in optical depth within a cell is no longer proportional to the increase in path length, as

$$\begin{array}{cc}\hfill {\tau }_{\mathrm{dic}}(s)& ={\tau }_{\mathrm{dic},m-1}+n_m{C}_{\mathrm{ext},m}(s-s_{m-1})\hfill \\ \hfill & -ln[cosh(n_m{C}_{\mathrm{pol},m}(s-s_{m-1}))\hfill \\ \hfill & -\frac{Q_{m-1}^{\prime }}{I_{m-1}^{\prime }}sinh(n_m{C}_{\mathrm{pol},m}(s-s_{m-1}))]\hfill \\ \hfill & s_{m-1}\le s\le s_m.\hfill \end{array}$$(37)

For a randomly generated τ, we should in principle use this equation to determine the correct value of s, which can be done using standard root-finding algorithms. In practice, however, we suggest using linear interpolation. Given the approximations due to the discretisation itself, the gain in accuracy by applying an exact root finding algorithm which is probably not worth the additional computational cost. Only in cases where the individual cells have a high optical depth could it possibly be useful to consider a more advanced, and numerically more costly, higher-order interpolation scheme.

One could argue that, in spite of the errors made, it would still be advantageous to use the extinction or total optical depth instead of the dichroic optical depth because the calculation of the dichroic optical depth is numerically more demanding. Indeed, calculating τ_ext(s) or τ_tot(s) involves just a single summation along the path, whereas the calculation of τ_dic(s) requires the propagation of the entire Stokes vector, including rotations and hyperbolic function evaluations at every grid cell. However, it should be recognised that this operation has to be executed anyway in the Monte Carlo loop since the calculation of the Stokes vector up to the next interaction point is required. This is due to the albedo and the scattering matrix that explicitly depend on the polarisation state of the radiation (Mishchenko et al. 2000; Wolf et al. 2002). Rather than being a numerically expensive extra, the dichroic optical depth calculation is free. We can conclude that there is no benefit at all in using τ_ext, τ_tot, or any other approximation, instead of the dichroic optical depth to calculate the next interaction location.

5. Summary and outlook

We have performed an analysis of the attenuation of radiation when it passes through a medium of aligned spheroidal grains, fully taking into account the effects of dichroism. The most important conclusions from this analysis are as follows.

Firstly, in a dichroic medium, the dichroic optical depth is no longer equivalent to the usual extinction optical depth τ_ext, that is, the integral of the product of number density and extinction cross section along the path. For representative values of the optical properties of dust grains, the relative difference between both optical depth scales can be several tens of percent, even at low optical depths.

A potential option to account for dichroic attenuation could be to replace the extinction cross section by the total extinction cross section. The corresponding total optical depth τ_tot approximates the dichroic optical depth to first order, but it is always overestimated. Relative differences between total and dichroic optical depth are small at low optical depths, but they can also run up to several tens of percent at high optical depths. An accurate calculation of the dichroic optical depth requires the full solution of the intensity profile along the path. In the general case of a dichroic medium, the radiative transfer equation becomes a set of four coupled first-order differential equations with varying coefficients, and a closed expression for the dichroic optical depth cannot be derived. However, the exact solution corresponding to a medium with uniform grain alignment can be used to find the full solution in an elegant way without any further numerical integration. There is no benefit in using τ_ext, τ_tot, or any other approximation instead of the dichroic optical depth to calculate the next interaction location in a Monte Carlo radiative transfer simulation.

Our results have implications for Monte Carlo radiative transfer coders that wish to incorporate the attenuation by elongated dust grains. If scattering polarisation by spherical grains already adds some complexity to Monte Carlo radiative transfer codes, dealing with non-spherical grains increases this complexity to a new level. Compared to spherical grains, the scattering process is significantly more complex. The scattering properties of spherical grains are fully described by just the albedo and the scattering phase function; for elongated grains, a full 4 × 4 scattering matrix comes into play (Mishchenko et al. 2000), and the random determination of a new propagation direction after a scattering event is not trivial (Wolf et al. 2002; Whitney & Wolff 2002; Lucas 2003). A related complexity concerns the amount of data that needs to be stored and accessed. Each of the elements of the extinction matrix and the scattering matrix is not only dependent on grain material, size, and wavelength, but also on shape and incidence angle. Moreover, each element of the scattering matrix needs to be discretised on the unit sphere. Finally, the process of dichroic extinction adds yet another level of complexity, and the results in this paper show that this also affects the random generation of the next interaction location.

So far, only a limited number of Monte Carlo radiative transfer codes have attempted to actually calculate dichroic attenuation by non-spherical aligned grains⁴. Wolf et al. (2002) were the first to include non-spherical aligned grains in their Monte Carlo radiative transfer calculations. They presented multiple light scattering calculations and demonstrated that the incorporation of elongated grains is important to explain the circular polarisation of light. They discussed the concepts of dichroism and birefringence, but they did not include these effects in the simulations they presented. Whitney & Wolff (2002) presented a Monte Carlo code that models the effects of scattering and dichroic absorption by aligned grains in circumstellar environments. While their code is presented for general geometries, they only discussed models with a uniform grain alignment. Simpson et al. (2013) applied this code to massive young stellar objects with more complex magnetic field configurations. Lucas (2003) presented a third Monte Carlo code with more or less the same characteristics as the code by Whitney & Wolff (2002), and also with the modelling of young stellar objects as the prime science objective. Applications of this code were presented by Lucas et al. (2004) and Chrysostomou et al. (2007). Peest et al. (in prep.) discuss the implementation of polarisation by elongated grains in the vectorised Monte Carlo code MC3D developed by (Krügel 2008; Heymann & Siebenmorgen 2012). Finally, a new Monte Carlo radiative transfer code, POLARIS, is presented by Reissl et al. (2016). It handles dichroic extinction and polarised emission, and it is optimised to handle data that results from sophisticated magneto-hydrodynamic simulations (Brauer et al. 2017; Reissl et al. 2017, 2018).

The work we present here fits into a broader effort to fully integrate the attenuation, polarisation, and thermal emission by elongated interstellar dust grains into the publicly available radiative transfer code SKIRT (Camps & Baes 2015). The advantage of implementing elongated grains in SKIRT is that the code can then use many of the useful ingredients that are already available, such as a suite of optimisation techniques (Baes et al. 2011, 2016), a library of input geometries for sources and sinks (Baes & Camps 2015), advanced spatial grids and grid traversal techniques (Camps et al. 2013; Saftly et al. 2013, 2014), the coupling to the output of grid-based and particle-based hydrodynamic codes (Saftly et al. 2015; Camps et al. 2016), and hybrid parallelisation techniques for shared and distributed memory machines (Verstocken et al. 2017).

Contrary to the currently available radiative transfer codes that incorporate elongated grains, SKIRT mainly focuses on galaxy-wide scales (e.g. De Looze et al. 2012, 2014; Viaene et al. 2017; Trayford et al. 2017). The magnetic and turbulent energy densities in nearby galaxies are found to be roughly in equipartition, and therefore magnetic fields are expected to be important for the evolution of galaxies (Boulares & Cox 1990; Beck et al. 1996). High-resolution cosmological zoom simulations have recently started to take magnetic fields into account (Pakmor et al. 2014, 2017; Grand et al. 2017). Radiative transfer codes that can fully incorporate dichroic extinction and emission by elongated aligned grains could be important tools to compare such simulations to observations.

References

All Figures

Fig. 1. Top left: relative differences δext between extinction optical depth τext and dichroic optical depth τdic for hypothetical example with density and grain alignment increasing linearly with increasing path length. The different curves correspond to different input Stokes vectors. In each case the photon package is initially fully linearly polarised, but the Stokes vector orientation θ increases in steps of 5°. The red lines corresponds to θ = 0° (Q0/I0 = 1) and the purple line to θ = 90° (Q0/I0 = −1), respectively. Bottom left: similar to top left panel, but now showing relative difference δtot between total optical depth τtot and dichroic optical depth. Right: same as panels on left, but now for model with uniform grain alignment.

In the text

Fig. 2. Histograms for distribution of randomly generated path lengths, corresponding to different optical depth scales discussed in this paper. The model is described in Sect. 4.1, and the different columns correspond to different levels of initial linear polarisation. Top panels: optical depths randomly generated from an exponential distribution, bottom panels: optical depths generated by using the composite biasing approach discussed in Baes et al. (2016).

In the text

	Fig. 3. Evolution of initially unpolarised Stokes vector propagating through medium of aligned grains where grain alignment rotates around path. The dots are the result from a Runge–Kutta numerical integration of the vector radiative transfer equation and the solid lines correspond to the method outlined in Sect. 3.
In the text