advanced radiative transfer in dusty systems


SKIRT supports the effects of polarization when radiation is scattered by spherical dust grains. To enable polarization support in a SKIRT simulation:

In this text we present the methods used to implement polarized random walk scattering and peel-off in SKIRT. We orient the Stokes parameters with a coordinate system that flies along the photon package and use a vector-based approach to calculate angles and directions. Thus we can circumvent numerical exceptions that arise in other codes.

Formulas relevant for polarization

In this section we present our choices of coordinate systems and the basics of the Stokes formalism. We use the IAU convention for polarization and right-handed coordinate systems.

Orientation basics

In this document we use three reference frames. The first one is the one describing the grid in which the simulation is set up, the laboratory frame in cartesian coordinates,

\[{\Sigma_{L}} = (\mathbf{e}_{x},\mathbf{e}_{y},\mathbf{e}_{z})\ .\]

The second frame is the beam frame, which describes a right handed orthonormal coordinate system that flies with a beam of photons. The \(z\)-axis of this coordinate system is along the propagation direction \(\mathbf{k}\) of the beam. The \(x\)-axis is the vector \({\mathbf{d}_{s}}\) along which the Stokes parameters (see below) are defined. The \(y\)-axis of this system is \({\mathbf{n}_{s}}\), the normal to \(\mathbf{k}\) and \({\mathbf{d}_{s}}\). For practical reasons, we will use \({\mathbf{n}_{s}}\) for most of the calculations. Since the coordinate system is orthonormal, using \(\mathbf{k}\) and \({\mathbf{n}_{s}}\) also determines \({\mathbf{d}_{s}}\).

\[{\Sigma_{B}}=({\mathbf{d}_{s}}, {\mathbf{n}_{s}}, \mathbf{k})\]

The third reference frame is the observer frame, which is used by the observing instrument, with its \(x\)- and \(y\)-axis ( \(\mathbf{k}_\text{x} \perp \mathbf{k}_\text{y}\)) determined from the set up of the simulation and its \(z\)-axis is pointing from the simulation towards the observer, \(\mathbf{k}_\text{obs}\). Since the FullInstrument uses parallel projection, the direction to the observer is perpendicular to the \(x\)- and the \(y\)-axis, making the system orthonormal.

\[{\Sigma_{O}}=(\mathbf{k}_\text{x}, \mathbf{k}_\text{y},\mathbf{k}_\text{obs})\]

(Because of the orthonormality requirement of the observer frame as defined above, the PerspectiveInstrument cannot be extended to support polarization).

The Stokes parameters

Consider a beam of radiation traveling through space in the beam reference frame \({\Sigma_{B}}\). The polarization status of the radiation beam is characterized by the Stokes vector

\[ \mathbf{S} = \left(\begin{array}{c} I\\ Q\\ U\\ V \end{array}\right) \]

where \(I\) denotes the total intensity, \(Q\) and \(U\) describe the linear polarization state relative to \({\mathbf{d}_{s}}\) (the \(x\)-axis of the beam reference system), and \(V\) is the degree of circular polarization. The Stokes vector spans the space of unpolarized, partially polarized and fully polarized light (they do not form a preferred basis for this space, but are useful because they can be easily measured and calculated).

Useful quantities that can be derived from the Stokes parameters include the degree of linear polarization,

\[ P_{\text{L}} = \frac{\sqrt{Q^2+U^2}}{I}\leq 1, \]

and the direction of linear polarization \(\gamma\), often depicted as

\[ \gamma = \frac12\arctan\left(\frac{U}{Q}\right). \]

(In order to conserve the quadrant, the atan2 function should be used to calculate \(\gamma\)). The degree of total polarization

\[ P = \frac{\sqrt{Q^2 + U^2+V^2}}{I} \]

can be calculated from the Stokes parameters as well. These definitions imply that \(\gamma\) is the angle, measured clockwise while looking along \(\mathbf{k}\), between the Stokes vector direction \({\mathbf{d}_{s}}\) and the major axis of the polarization ellipse. For linearly polarized radiation ( \(V=0\)), the polarization ellipse is a line along \({\mathbf{d}_{s}}\) for \(Q>0\) and \(U=0\) and \(45^\circ\) from \({\mathbf{d}_{s}}\) for \(Q=0\) and \(U>0\).

The relation between the Stokes parameters and the parameters of the polarization ellipse can be presented by considering a 3D space, and looking at the vector with 3D Cartesian coordinates \((Q,U,V)\),

\[ \begin{split} Q &= I P\cos2\gamma \cos2\chi \\ U &= I P\sin2\gamma \cos2\chi \\ V &= I P\sin2\chi \end{split} \]

with the shape parameter \(\chi\) of the polarization ellipse. In other words, \(IP\), \(2\gamma\) and \(2\chi\) are the spherical coordinates of the 3D vector of Cartesian coordinates \((Q,U,V)\). The factor 2 in front of \(\gamma\) accounts for the fact that a polarization ellipse remains invariant for rotation by 180 \(^\circ\), whereas the factor 2 in front of \(\chi\) indicates that the polarization ellipse is invariant if we swap the lengths of the axes and simultaneously apply a rotation of 90 \(^\circ\).

Combining the latter equations with the formula for the degree of linear polarization we can write

\[ \begin{split} Q &= I P_{\text{L}} \cos2\gamma \\ U &= I P_{\text{L}} \sin2\gamma \ . \end{split} \]

Rotating the reference vector of the Stokes parameters

The Stokes parameters are defined with respect to the choice of \({\mathbf{d}_{s}}\) and \({\mathbf{n}_{s}}\) in the plane perpendicular to the propagation direction. If we rotate the photon coordinate system around its \(z\)-axis, the Stokes parameters will also alter. The total intensity \(I\) and the parameter \(V\) are invariant for such a rotation, but \(Q\) and \(U\) are not (logical, as they are defined with respect to \({\mathbf{d}_{s}}\)). If we rotate the axes \({\mathbf{d}_{s}}\) and \({\mathbf{n}_{s}}\) clockwise about \(\mathbf{k}\) over an angle \(\phi\), while looking towards the propagation destination, and we denote the Stokes vector in the new coordinate systems as \(\mathbf{S}'\), then

\[ \mathbf{S}' = {\bf{R}}(\phi)\,\mathbf{S} \]


\[ \left(\begin{array}{c} I\\ Q\\ U\\ V \end{array}\right) = \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & \cos2\phi & \sin2\phi & 0 \\ 0 & -\sin2\phi & \cos2\phi & 0 \\ 0 & 0 & 0 & 1 \end{array}\right) \left(\begin{array}{c} I\\ Q\\ U\\ V \end{array}\right) \]

Note the appearance of the factor 2 which is a logical consequence of the choice of the Stokes parameters as a basis for the space of polarised radiation.

Scattering in the Stokes formalism

When a radiation beam scatters, it changes direction and polarization status. Assume that the propagation direction before the scattering is denoted as \(\mathbf{k}\) and the one after the scattering as \(\mathbf{k}'\), with \(\theta\) the actual scattering angle,

\[ \cos\theta = \mathbf{k} \cdot \mathbf{k}'. \]

The change of the Stokes vector during a scattering event is described by Mueller matrices. These are defined with the Stokes vector direction \({\mathbf{d}_{s}}\) in the scattering plane (and thus \({\mathbf{n}_{s}}\) being the normal of the scattering plane). In general these matrices are

\[ {\bf{M}}(\theta) = \frac{1}{k^{2}r^{2}} \left(\begin{array}{cccc} S_{11} & S_{12} & S_{13} & S_{14} \\ S_{21} & S_{22} & S_{23} & S_{24} \\ S_{31} & S_{32} & S_{33} & S_{34} \\ S_{41} & S_{42} & S_{43} & S_{44} \end{array}\right) \]

and each of the elements of the Mueller matrix depends on the scattering angle \(\theta\). For MC simulations the factor \(\frac{1}{k^{2}r^{2}}\) is dropped, as the path of one photon package is followed (Bianchi et al. 1996).

Since there are multiple restrictions on Stokes vectors (like \(I^{2} \geq Q^{2}+U^{2}+V^{2}\)), there are at most seven independent variables (Kruegel 2002). For example for scattering at a sphere the formula above reduces to

\[ {\bf{M}}(\theta) = \left(\begin{array}{cccc} S_{11} & S_{12} & 0 & 0 \\ S_{12} & S_{11} & 0 & 0 \\ 0 & 0 & S_{33} & S_{34} \\ 0 & 0 & -S_{34} & S_{33} \end{array}\right)\ . \]

The Stokes vector after scattering \(\mathbf{S}'\) is calculated by multiplying the Stokes vector with the Mueller Matrix,

\[ \mathbf{S}' = \bf{M}\,\mathbf{S} \]

Implementation of the Scattering formalism

In this section we expound the steps that are necessary for a correct implementation of the scattering formalism in a 3D-MC radiative transfer code. First we explain our choice coordinate system. After this we describe how to obtain a new propagation direction when scattering. For this we need to calculate the phase function \(\Phi\), which we use to sample the angles \((\theta,\varphi)\). The state of the Stokes vector after scattering is derived by this pair of angles. The new propagation and reference directions are derived from the angles as well, with an explanation of the special case of peel-off photons towards the observer.

Orientation of the Stokes parameters

There are two different approaches/conventions to implement orientating Stokes parameters in between scatterings.

Orientation w.r.t. the meridional plane

Consider the scattering geometry as used in Code and Whitney (1995),


and the one used in Bianchi et al (1996),


Note that the angle \(i_{2}\) in the first illustration is dubiously oriented, it should probably be on the opposite side of the meridian.

The first option is to follow the original framework of Chandrasekhar (1960), as is also done by, for example, Code and Whitney (1995) and Gordon et al. (2001). They use a coordinate system where the orientation of the polarization ellipse is always stored with the direction along the meridional plane as the reference direction \({\mathbf{n}_{s}}\). In this case, the simulation of a scattering event consists of three consecutive operations. The first one is a rotation of the reference system so that the Stokes vector refers to the direction parallel to the scattering plane. The second step is the actual scattering by applying the Mueller matrix. The last step is to rotate the Stokes vector again so it refers to the direction along the meridional plane. This requires another rotation. In total we hence obtain

\[ \mathbf{S}' = \bf{R}(i_2-\pi) \, \bf{M}(\theta) \, \bf{R}(i_1)\,\mathbf{S} \]

(There is a sign difference to the book, as our rotation direction convention is different)

The same framework with a different notation is that of Bianchi et al. (1996). The reference \({\mathbf{n}_{s}}\) is first rotated by an angle \(\phi\), then the photon scatters and afterwards the reference frame is rotated back by \(-\gamma\). Comparing the first two pictures of this section, one can make the connection \(\phi = i_1\) and \(\gamma = \pi - i_2\). Applying this connection, the two formulas are equivalent.

Orientation w.r.t. the previous scattering plane

The scattering geometry as used in SKIRT:


The Chandrasekhar (1960) framework explained above has two drawbacks. On the one hand, there is a problem if the scattering occurs in the direction of the poles, as one cannot define the meridian in this case. A second problem is one of redundancy: for every scattering event, there are three operations that need to be performed: a rotation by \(i_1\) to go to the scattering plane, an actual scattering in the scattering plane, and a rotation by \(i_2-\pi\) to go back to the meridional frame. If several scattering events are chained, this means that we have two rotations in a row between each actual scattering: a first one by an angle \(i_2^{(1)}-\pi\) at the end of scattering event #1, and a second one by an angle \(i_1^{(2)}\) at the start of scattering event #2. Rather than executing these two rotations, it makes more sense to apply a single rotation over the angle

\[ \varphi \equiv i_2^{(1)}-\pi+i_1^{(2)}. \]

The angle \(\varphi\) is also the angle between the scattering planes corresponding to the two scattering events. So a scattering event is characterized as

\[ \mathbf{S}' = \bf{M}(\theta) \, \bf{R}(\varphi) \, \mathbf{S} \]

This essentially comes down to a different convention: rather than always using Stokes parameters that are defined with respect to the meridional plane, we just use Stokes vectors that are defined with respect to the previous scattering plane. This is the convention applied by Fischer et al. (1994) and Goosmann & Gaskell (2007), and we also apply it in SKIRT.

One special case is the first scattering, as there is not yet a previous scattering plane in this case. But this does not really matter: since the radiation is unpolarized before the first interaction, \(Q_{\text{in}} = U_{\text{in}} = V_{\text{in}} = 0\) and the rotation operator \({\bf{R}}(\varphi)\) is irrelevant.

Scattering in a random walk

Assume we want to simulate a scattering event in a Monte Carlo random walk step at a given location. We know the propagation direction of the photon package before the scattering event \(\mathbf{k}\), the reference direction \({\mathbf{n}_{s}}\) and the Stokes vector \(\mathbf{S}\). We want to generate a random new direction of the photon package, \(\mathbf{k}'\), with a new reference direction \({\mathbf{n}_{s}}'\) and the new Stokes vector \(\mathbf{S}'\). The information which new direction has which probability to occur is represented by the phase function \(\Phi\).

Contrary to unpolarized radiation, where the phase function only depends on the scattering angle \(\theta\), it is more complicated in the case of polarized radiation, as it depends on both \(\theta\), the azimuth \(\varphi\) and on the polarization status of the incoming photon package. To sample \(\theta\) and \(\varphi\), we follow the steps:

Calculating the phase function

The phase function is proportional to the intensity of the beam after the scattering event, \(I_{\text{out}}\) divided by the intensity \(I_{\text{in}}\) before the event.

\[ \Phi(\theta,\varphi) \propto \frac{I_{\text{out}}}{I_{\text{in}}} \]

The intensity after the scattering event can be calculated using multiplying the Stokes vector with the Mueller matrix for all tuples \((\theta,\varphi)\). Assuming the simplified Mueller matrix for scattering at a sphere, the first Stokes parameter after scattering is

\[ I_{\text{out}} = I_{\text{in}} S_{11} + S_{12} \left( Q_{\text{in}}\cos2\varphi+U_{\text{in}}\sin2\varphi \right) \]


\[ \Phi(\theta,\varphi) \propto S_{11} + S_{12}\left( \frac{Q_{\text{in}}}{I_{\text{in}}}\cos2\varphi+\frac{U_{\text{in}}}{I_{\text{in}}}\sin2\varphi \right) \ \ , \]

Using the equation for the degree of linear Polarization, we find

\[ \Phi(\theta,\varphi) \propto S_{11} + P_{\text{L,in}}\,S_{12}\cos2(\varphi - \gamma_{\text{in}}) \]

After normalization of the right side, both sides are equal.

Normalizing the phase function

There are different normalization conventions. We choose the Integral of the phase function over over the unit sphere to be normalized to \(4\pi\) sr. Another convention would be to normalize it to one. Depending on the implementation in the code, the normalization constant \(N\) needs to be adjusted. To normalize the phase function to \(4 \pi\), we divide \(4 \pi\) by the integral over the unit sphere,

\[ \begin{split} N&=\frac{4 \pi}{ \int_0^{2\pi} \int_0^{\pi} \left(S_{11}(\theta) + P_{\text{L,in}}S_{12}(\theta)\cos2(\varphi - \gamma_{\text{in}})\right)\sin\theta\, \text{d}\theta\, \text{d}\varphi }\\ &=\frac{2}{\int_0^\pi S_{11}(\theta)\sin\theta\, \text{d}\theta} \ . \end{split} \]

We can write the properly normalized phase function as

\[ \Phi(\theta,\varphi) = N\,S_{11}(\theta) \left[ 1 + P_{\text{L,in}}\,\frac{S_{12}(\theta)}{S_{11}(\theta)}\cos2(\varphi - \gamma_{\text{in}}) \right] \]

Sampling the scattering angle

In order to sample a random angle from the distribution function, we use the conditional probability technique. For this we reduce the phase function to the marginal distribution \(\Phi(\theta)\)

\[ \Phi(\theta) = \int_0^{2\pi} \Phi(\theta,\varphi)\,\text{d}\varphi = 2\pi\ N\,S_{11}(\theta) \]

Sampling a random \(\theta\) from this distribution can be done using a numerical inversion, i.e. by solving the equation

\[ {\cal{X}} =\frac{2\pi\,N \int_0^\theta S_{11}(\theta')\sin\theta'\,\text{d}\theta'}{4\pi} \]

for \(\theta\), where \({\cal{X}}\) is a uniform deviate, a random number between \(0\) and \(1\) with a uniform distribution function.

Sampling the azimuth angle

Once we know the scattering angle \(\theta\) with which we will scatter, we determine the random azimuth angle \(\varphi\) from the normalized conditional distribution

\[ \Phi_\theta(\varphi) = \frac{\Phi(\theta,\varphi)}{\int_0^{2\pi} \Phi(\theta,\varphi')\,\text{d}\varphi'} = \frac{1}{2\pi}\left[1+ P_{\text{L,in}}\,\frac{S_{12}(\theta)}{S_{11}(\theta)}\cos 2(\varphi - \gamma_{\text{in}})\right] \]

This can again be done using numerical inversion, so by solving the equation

\[ \begin{split} {\cal{X}} &= \frac{\int_{0}^{\varphi}\Phi_{\theta}(\varphi')\,\text{d}\varphi'}{2\pi} \\ &= \frac{1}{2\pi} \left[ \varphi + P_{\text{L,in}}\,\frac{S_{12}(\theta)}{S_{11}(\theta)} \sin\varphi \cos(\varphi - 2\gamma_{\text{in}})\right] \end{split} \]

for \(\varphi\) with \({\cal{X}}\) a new uniform deviate. Once this random value for \(\varphi\) has been obtained, we have the tuple \((\theta,\varphi)\) with which we will scatter. For a photon package with no linear polarization ( \(P_{\text{L,in}} = 0\)) all \(\varphi\) have the same probability.

New propagation direction

Now that we know the scattering angle tuple \((\theta,\varphi)\), we need to calculate the new direction of the photon package. For this we use the propagation direction \(\mathbf{k}\), the reference direction \({\mathbf{n}_{s}}\) and Rodrigues' rotation formula. It describes how any vector \(\mathbf{v}\) can be rotated around any axis \(\mathbf{a}\) by any angle \(\beta\) (clockwise while looking along \(\mathbf{a}\)), with the formula

\[ \mathbf{v}' = \mathbf{v} \cos\beta+ (\mathbf{a} \times \mathbf{v})\sin\beta+\mathbf{a}(\mathbf{a} \cdot \mathbf{v})(1-\cos\beta) \]

This formula is used to first calculate the new reference direction and then to calculate the new propagation direction.

Possible calculation of first reference direction

In case there is no reference direction \({\mathbf{n}_{s}}\) (usually because it is the first scattering event), it must be generated. There are \(2 \pi\) directions that are perpendicular to the propagation direction \(\mathbf{k}\). One could assume the choice of the reference direction could influence the randomness of the random walk. Luckily, if there is no reference direction \({\mathbf{n}_{s}}\), there is also no Stokes vector direction \({\mathbf{d}_{s}}\) (otherwise \({\mathbf{n}_{s}} = \mathbf{k} \times {\mathbf{d}_{s}}\)) and the linear polarization degree \(P_{\text{L,in}}\) is zero ( \(Q\) and \(U\) are only meaningful with a direction \({\mathbf{d}_{s}}\)). Applying this to the sampling of the azimuth results in a uniform probability distribution for the azimuth \(\varphi\). Our choice of the reference direction is therefore completely randomized when the azimuth is applied to generate the new direction.

We calculate a reference \({\mathbf{n}_{s}} = ({\mathbf{n}_{s}}_{x},{\mathbf{n}_{s}}_{y},{\mathbf{n}_{s}}_{z})\) that is perpendicular to the current direction of flight \(\mathbf{k} = (\mathbf{k}_{x},\mathbf{k}_{y},\mathbf{k}_{z})\),

\[ \begin{split} {\mathbf{n}_{s}}_{x} &= -\mathbf{k}_{x}\mathbf{k}_{z}/\sqrt[]{1-\mathbf{k}_z^{2}}\\ {\mathbf{n}_{s}}_{y} &= -\mathbf{k}_{y}\mathbf{k}_{z}/\sqrt[]{1-\mathbf{k}_z^{2}}\\ {\mathbf{n}_{s}}_{z} &= \sqrt[]{1-\mathbf{k}_z^{2}} \end{split} \]

these formulas are adapted from Bianchi et al. (1996) for an angle of \(90^\circ\) towards the direction of flight. In case \(\mathbf{k}\) is towards \(\pm \mathbf{e}_{z}\), the root degenerates and the reference direction \({\mathbf{n}_{s}} = \mathbf{e}_{x}\) is chosen instead.

Calculating the new reference direction


Depiction of the vectors and planes used here.

To calculate the new reference direction \({\mathbf{n}_{s}}'\), we use the physical meaning of the azimuth \(\varphi\). It describes the angle between the scattering planes, the angle around which the Stokes parameters are rotated so they are referring to an axis inside the scattering plane. As the reference direction \({\mathbf{n}_{s}}\) is perpendicular to the Stokes vector direction and the propagation direction \(\mathbf{k}\), it is perpendicular to the past scattering plane. For it to be perpendicular to the current scattering plane, it needs to be rotated by \(\varphi\) with the propagation direction as the rotation axis. The resulting vector is the new reference direction \({\mathbf{n}_{s}}'\). This rotation is calculated by using Rodrigues' Rotation formula,

\[ {\mathbf{n}_{s}}' = {\mathbf{n}_{s}} \cos\varphi+ (\mathbf{k} \times {\mathbf{n}_{s}})\sin\varphi+\mathbf{k}(\mathbf{k}\cdot{\mathbf{n}_{s}})(1-\cos\varphi) \]

As \({\mathbf{n}_{s}}\) is perpendicular to \(\mathbf{k}\), the third term is zero,

\[ {\mathbf{n}_{s}}' = {\mathbf{n}_{s}} \cos\varphi+ (\mathbf{k} \times {\mathbf{n}_{s}})\sin\varphi \]

and \({\mathbf{n}_{s}}'\) is the normal of the current scattering plane.

Calculating the new propagation direction

The second step is to rotate the current propagation direction \(\mathbf{k}\) around the normal vector of the current scattering plane \({\mathbf{n}_{s}}'\) by \(\theta\). The resulting vector is the new propagation direction \(\mathbf{k}'\). This is implemented the same way, using Rodrigues' Rotation formula and again the third term is zero,

\[ \mathbf{k}' = \mathbf{k} \cos\theta+ ({\mathbf{n}_{s}} \times \mathbf{k})\sin\theta \ . \]

Thus we obtain the new propagation direction \(\mathbf{k}'\) and the new reference direction \({\mathbf{n}_{s}}'\).

Scattering peel-off

Depending on the code there are photons peeled off towards the observer(s) whenever the photon package scatters (and in SKIRT this is the case). These photons need to convey the correct Stokes parameters to the observer. In contrast to random walk scattering the future direction is known, it is \(\mathbf{k}_\text{obs}\), while the scattering angles \((\theta_\text{obs},\varphi_\text{obs})\) are unknown. After calculating the direction towards the observer, the Stokes reference direction must be rotated so it refers to the \(y\)-axis of the observer frame, \(\mathbf{k}_\text{y}\).

Determining the direction angles towards the observer

The scattering angle \(\theta_\text{obs}\) is easily determined through the cosine of the peel-off scattering angle,

\[ \cos\theta_\text{obs} = \mathbf{k} \cdot \mathbf{k}_\text{obs} \]

as the scattering angle runs from \(0\) to \(\pi\), this relation is bijective and the \(\arccos\) can be used.

To compute \(\varphi_\text{obs}\), recall that it is the angle between the previous scattering plane and the current scattering plane, or equivalently, the angle between the normal vectors to these scattering planes. The normal of the last scattering plane is \({\mathbf{n}_{s}}\) (which might need to be generated, if it is the first scattering, using Calculating the new reference direction ). The normal of the current scattering plane \(\mathbf{n}\) can be calculated using the current propagation direction \(\mathbf{k}\) and the direction towards the observer \(\mathbf{k}_\text{obs}\),

\[ \mathbf{n} = \frac{\mathbf{k} \times \mathbf{k}_\text{obs}} {||\mathbf{k} \times \mathbf{k}_\text{obs}||} \]

Thus we have

\[ \cos\varphi_\text{obs} = {\mathbf{n}_{s}} \cdot \mathbf{n} \]

But since \(\varphi_\text{obs}\) can run from \(0\) to \(2 \pi\), the relation is not bijective. We need the sine of \(\varphi_\text{obs}\) as well. Since \(\mathbf{k}\) is perpendicular to both \({\mathbf{n}_{s}}\) and \(\mathbf{n}\), the relation

\[ \sin\varphi_\text{obs} \,\mathbf{k}= {\mathbf{n}_{s}} \times \mathbf{n} \]

holds. Or, after projecting both sides of the equation on \(\mathbf{k}\),

\[ \sin\varphi_\text{obs} = ({\mathbf{n}_{s}} \times \mathbf{n}) \cdot \mathbf{k} \]

so that \(\varphi_\text{obs}\) is fully defined. Once the angles are determined, the Stokes parameters are rotated by \(\varphi_\text{obs}\) (which updates \({\mathbf{n}_{s}}\) to \({\mathbf{n}_{s}}'\) as well). The Mueller matrix is applied with the scattering angle \(\theta_\text{obs}\) (and \(\mathbf{k}\) updated to \(\mathbf{k}'=\mathbf{k}_\text{obs}\)).

Calculating the normal of the current scattering plane for \(\mathbf{n}\) will cause an exception if the peel off is completely forwards or backwards, but in this case setting \(\varphi_\text{obs} = 0\) is reasonable, as the scattering plane does not need to be updated when scattering forwards or backwards.

Determining the orientation angle at the observer

Once the photon is oriented towards the observer, the Stokes vector reference direction must be oriented so the Stokes parameters refer to the \(y\)-axis of the observer frame, \(\mathbf{k}_\text{y}\), as their reference axis. This is equivalent to the reference direction after peeling off, \({\mathbf{n}_{s}}'\), being oriented towards \(\mathbf{k}_\text{x}\). The angle between \({\mathbf{n}_{s}}'\) and \(\mathbf{k}_\text{x}\) can be determined equivalently to the determination of \(\varphi_\text{obs}\), using

\[ \cos\alpha_\text{obs} = {\mathbf{n}_{s}}' \cdot \mathbf{k}_\text{x} \]


\[ \sin\alpha_\text{obs} = ({\mathbf{n}_{s}}' \times \mathbf{k}_\text{x}) \cdot \mathbf{k}' \]

with which \(\alpha_\text{obs}\) is fully defined.

Testing the simulation


A simple geometry to test the polarization method. A single point like source shines on thin slabs of electrons. In the tests for random walk and for scattering plane rotation the star is replaced with a single cell of electrons that is illuminated from different angles and is the first position of scattering.

In this section we show test cases to ensure the correct functionality of the implementation. They are designed to gradually use more of the routines so the incorrect implementation can be identified.

A very simple case of polarization arises from Thomson scattering of light at electrons. It is wavelength independent and the light cannot be absorbed. Only the direction and polarization state are changed. The simplified Mueller matrix for spheres is further reduced to

\[ {\bf{M}_{Th}}(\theta) = \left(\begin{array}{cccc} \frac{\cos^{2}\theta + 1}{2} & \frac{\cos^{2}\theta - 1}{2} & 0 & 0 \\ \frac{\cos^{2}\theta - 1}{2} & \frac{\cos^{2}\theta + 1}{2} & 0 & 0 \\ 0 & 0 & \cos\theta & 0 \\ 0 & 0 & 0 & \cos\theta \end{array}\right)\ . \]

Together with simple 2D-geometries like in the figure above, the resulting polarization state can be calculated analytically and compared to the output generated by SKIRT.

Testing the conservation of flux


Flux levels for different wavelengths and different optical depths of a ball of electrons around a central star. The flux is normalized to the original flux of the star.

If the implementation of the polarization routines is incorrect, photons might get lost or their intensity get changed. This can be tested by putting a star in an optically thick ball of electrons and checking whether the flux differs from a star without electrons around it. The figure above depicts the dependency of the total flux on the optical thickness \(\tau\). To better visualize it, the flux is normalized to the level of a simulation with no electrons around the star. As can be seen, the flux remains within 1 % of the original flux up to an optical thickness \(\tau=30\). As in this simulation there are only 50x50x50 grid cells in the whole simulation, the radius has around 25 cells for the optical thickness. Thus it seems that the optical thickness of a cell should not exceed unity, otherwise the flux starts to increase (to 105% at \(\tau_{\text{cell}}=4\)). Indeed, replications with a finer grid reveal that the flux is conserved better when using more cells.

Testing polarized peel off


Polarized peel off test. Comparison of the analytical result (blue) and the simulation with SKIRT (red).

Unpolarized light from a point-like source at \(P_{\star}\) shines onto thin slabs of electrons and is peeled off towards the observer (see Figure). As the slabs are very thin (optically and physically), the light is scattered only once and the observer can deduce the scattering angle \(\theta\) from the position \(x\) of the pixel. Easily compared are the degree of linear polarization \(P_{\text{L}}\), the direction of the polarization \(\gamma\) and the relative curve of the intensity \(I_\text{rel}\). For the calculation of \(I_\text{rel}\), the path length \(l\) to the slab must be taken into account as well. (Not from the slab to the observer, as the observer is at an `infinite' distance and thus all distances towards the observer are equal.)

For a position \(x \in [-x_{max}, x_{max}]\) at the observer, the angle \(\theta\) the photon was peeled off with at the slab can be deduced by geometrical reasoning to be

\[ \begin{split} \theta &= \frac{\pi}{2} \pm \arctan(x_{max}/|x| - 1) \end{split} \]

with the plus sign right of the origin and the minus sign left of it. The resulting curves are depicted in the figure above. SKIRT shows excellent agreement with the analytical solution. The intensity has some noise, but this seems to be noise rather than an error.

Testing the random walk routine


Test of the functionality of the random walk routine (see Testing the random walk routine ). The source illuminating the central electrons (see Figure) is pointing towards \(45^\circ\) south of the \(x\)-axis (to the bottom right).

The simple geometry can be used to check the functionality of the random walk routine as well. For this the central star is replaced with an electron cell in the center of the simulation. It is illuminated by a source that is placed in the plane of the slabs. To prevent the source from illuminating the slabs directly, the source's photons are all pointed towards the center. This way the light does a random walk scattering at the central object and the light is already polarized when it hits the electron slabs. As all components are in the same plane ( \(z \equiv 0\)), the scattering plane is always the \(xy\)-plane and thus this test does not rely on the functionality of the routines rotating the scattering plane. (This is not completely correct, as the orientation of the scattering plane will flip from \(+z\) to \(-z\) or vice versa for about half the pixels and it needs the implementation to determine the first reference direction \({\mathbf{n}_{s}}\). Nevertheless, this only relies on a fraction of the routines and only on the most basic cases. A flip by \(180^\circ\) does not change the Stokes parameters.) The two angles \((\theta_1,\theta_2)\) a photon takes to the observer are

\[ \begin{split} \theta_1 &=\theta_2 \pm (-45^\circ)\\ \theta_2 &= 90^{\circ} \pm \arctan(x_{max}/|x| - 1) \end{split} \]

with the plus signs right of the origin and the minus signs left of it. As can be seen in the figure above, the results of the analytical solution and SKIRT are equal (save some noise).

Testing the scattering plane rotation


Test of the scattering plane rotation routine. The source illuminating the central electrons (see Figure) is pointing towards \(15^\circ\) from the zenith (+ \(z\)-axis) and \(60^\circ\) from the \(y\)-axis (towards the top right).

By moving the source from the previous test out of the \(z=0\)-plane of the slabs, we can test the scattering plane rotation routines. The normal of the first scattering plane is tilted versus the \(z\)-axis, while the normal of the second is aligned with the \(z\)-axis. This makes a scattering plane rotation necessary.

As in this case obtaining the angles from geometrical reasoning is nearly impossible, the `expected' curves in the figure above are actually calculated using the three steps a photon takes to get to the observer (to the central source, to the slabs, to the observer) as vectors. The angles are determined from the vectors and the Stokes parameters are manipulated accordingly. This is similar to how SKIRT calculates some of the angles.