 |
The SKIRT project
advanced radiative transfer for astrophysics
|
|
The SKIRT project
advanced radiative transfer for astrophysics
|
This page describes how SKIRT implements the effects on the polarization of radiation when it is scattered by spherical entities including dust grains and electrons. To actually enable polarization support in a SKIRT simulation, see Configuring for polarized radiation.
We first present the mathematical background of the formalism used in SKIRT. We employ coordinate-system-independent vector representations for directions and angles. This avoids numerical issues that may arise with angles referring to a fixed coordinate system when the propagation direction is (nearly) parallel to one of the coordinate axes. We further use the IAU conventions for polarization and right-handed coordinate systems.
We consider three different 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 instrument's configuration and its \(z\)-axis is pointing from the model origin towards the observer, \(\mathbf{k}_\text{obs}\). Assuming parallel projection as for all DistantInstrument subclasses, 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})\]
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 (it does not form a preferred basis for this space, but it is useful because its components 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} \]
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 \(\varphi\), 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}}(\varphi)\,\mathbf{S} \]
or
\[ \left(\begin{array}{c} I\\ Q\\ U\\ V \end{array}\right) = \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & \cos2\varphi & \sin2\varphi & 0 \\ 0 & -\sin2\varphi & \cos2\varphi & 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 polarized radiation.
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 Monte Carlo simulations the factor \(\frac{1}{k^{2}r^{2}}\) is dropped, as the path of one photon package is followed.
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. 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. Including the reference direction adjustments before and after the actual scattering event this yields,
\[ \mathbf{S}' = {\bf{R}}(\varphi'')\,\bf{M}(\theta)\,\mathbf{S}\,{\bf{R}}(\varphi')\]
As described in the previous section, each scattering event requires three operations: a rotation by \(\varphi'\) to adjust the orientation to the scattering plane, the Mueller transformation in the scattering plane, and a rotation by \(\varphi''\) to adjust the orientation to some new reference frame. If there are multiple consecutive scattering events, this means that we have two rotations in a row between each actual scattering. Rather than executing these two rotations, it makes more sense to combine these into a single rotation over the angle between the scattering planes corresponding to the two scattering events. A scattering event is then characterized as
\[ \mathbf{S}' = \bf{M}(\theta) \, \bf{R}(\varphi) \, \mathbf{S}. \]
In other words, rather than using some fixed reference direction for the Stokes vector, we simply use the previous scattering plane.
A special case is the first scattering event after emission, 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.
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:
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 by 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) \]
and therefore
\[ \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}}). \]
In SKIRT, the integral of the phase function over over the unit sphere should be normalized to \(4\pi\). Thus, in our case, the normalization constant \(N\) is given by,
\[ \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] \]
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.
Once we know the scattering angle \(\theta\), 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.
Now that we know the scattering angle tuple \((\theta,\varphi)\), we need to calculate the new direction of the photon package given its previous propagation direction \(\mathbf{k}\) and reference direction \({\mathbf{n}_{s}}\). We use Rodrigues' rotation formula that allows rotating any vector \(\mathbf{v}\) around any axis \(\mathbf{a}\) by any angle \(\beta\) (clockwise while looking along \(\mathbf{a}\)),
\[ \mathbf{v}' = \mathbf{v} \cos\beta+ (\mathbf{a} \times \mathbf{v})\sin\beta+\mathbf{a}(\mathbf{a} \cdot \mathbf{v})(1-\cos\beta). \]
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.
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}}'\).
Whenever a photon package scatters, new photon packets are peeled off towards the observer(s). In contrast to random walk scattering the new direction is already known, it is \(\mathbf{k}_\text{obs}\), while the scattering angles \((\theta_\text{obs},\varphi_\text{obs})\) are unknown and must be calculated. Also, before detecting the photon packet at the observer, the Stokes reference direction must be rotated so it refers to the \(y\)-axis of the observer frame, \(\mathbf{k}_\text{y}\).
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}}\). 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.
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} \]
and
\[ \sin\alpha_\text{obs} = ({\mathbf{n}_{s}}' \times \mathbf{k}_\text{x}) \cdot \mathbf{k}' \]
with which \(\alpha_\text{obs}\) is fully defined.