Paradigmatic examples for testing models of optical light polarization by spheroidal dust

TL;DR

This paper develops a comprehensive 3D Monte Carlo radiative transfer framework to model and test the polarization effects caused by spheroidal dust grains in astrophysical environments, including scattering, dichroic extinction, and birefringence.

Contribution

It introduces a new methodology for implementing polarization modeling of spheroidal dust in radiative transfer codes, validated by analytical benchmarks.

Findings

MCpol code achieves 0.1% accuracy in polarization simulations.

The framework effectively models polarization changes due to scattering, extinction, and birefringence.

Benchmark tests confirm the numerical precision of the implementation.

Abstract

We present a general framework on how the polarization of radiation due to scattering, dichroic extinction, and birefringence of aligned spheroidal dust grains can be implemented and tested in 3D Monte Carlo radiative transfer (MCRT) codes. We derive a methodology for solving the radiative transfer equation governing the changes of the Stokes parameters in dust-enshrouded objects. We utilize the M\"uller matrix, and the extinction, scattering, linear, and circular polarization cross sections of spheroidal grains as well as electrons. An established MCRT code is used and its capabilities are extended to include the Stokes formalism. We compute changes in the polarization state of the light by scattering, dichroic extinction, and birefringence on spheroidal grains. The dependency of the optical depth and the albedo on the polarization is treated. The implementation of scattering by spheroidal grains both for random walk steps as well as for directed scattering (peel-off) are described. The observable polarization of radiation of the objects is determined through an angle binning method for photon packages leaving the model space as well as through an inverse ray-tracing routine for the generation of images. We present paradigmatic examples for which we derive analytical solutions of the optical light polarization by spheroidal dust particles. These tests are suited for benchmark verification of MCpol and other such codes, and allow to quantify the numerical precision reached. We demonstrate that MCpol is in excellent agreement to within 0.1% of the Stokes parameters when compared to the analytical solutions.

Tables

Table 1. Table 1: Precision of MCpol for TC 4 by varying the number of sub-cubes and emitted photons. The minimum, maximum and the standard deviation σ 𝜎 \sigma of the model result minus the analytical solution in V / I 𝑉 𝐼 V/I is specified using Δ Δ \Delta as in Eq. ( 50 ) .

$N_{\mathrm{sub}}$	$N_{\gamma }$	$\min$	$\max$	$\sigma$
		$\left({10}^{-4}\right)$
1	$2.5\times {10}^9$	-35	28	18
1	$2.5\times {10}^{10}$	-26	23	18
3	$2.5\times {10}^9$	-63	64	20
3	$2.5\times {10}^{10}$	-25	22	9
11	$2.5\times {10}^{10}$	-16	13	5
11	$2.5\times {10}^{11}$	-9	10	3

Table 2. Table 2: Notation

Symbol	Description
$i,j,m$	numbering indices
$(x,y,z)$	axis of Cartesian system
$\alpha$	angle of incidence, Fig. 3
$\beta$	rotation angle, esp. between North direction and the plane of incidence, $\beta :=\sphericalangle ({𝒅}_{\text{E}},𝒏)$
$\gamma$	exit angle, Fig. 3
$(\theta$ , $\phi )$	pair of scattering angles
$({\theta }_{\text{obs}}$ , ${\phi }_{\text{obs}})$	pair of angles to observer
$N_{\theta }$, $N_{\phi }$	number of bins along scattering angles
$\zeta$	random number
$v_{ceil}$	“ceiling” value for rejection sampling (maximum probability)
$P(\theta ,\phi )$	probability density function
$s$, $s^{\prime }$	length of photon path through medium
$\Delta s_j$	length of incremental photon path in cube $j$
$\Delta$	${10}^4$ times the difference of $V/I$ of MCpol minus the analytical solution by P17
$𝒓$	position vector in the medium
$𝒐$	symmetry axis of grain along major axis, Fig. 3
$𝒌$	incoming photon direction, Fig. 3
${𝒌}^{\prime }$	outgoing photon direction
${𝒌}_{\text{obs}}$	photon direction towards observer
${𝒅}_{\text{N}}$	North direction in the particle frame (plane of incidence)
${𝒅}_{\text{N},init}$	initial North direction, Fig. 3
${𝒅}_{\text{N},sca}$	North direction in the plane of scattering before the scattering event, Fig. 3
${𝒅}_{\text{N}}^{\prime }$	North direction in particle frame (plane of departure), Fig. 3
${𝒅}_{\text{N},sca}^{\prime }$	North direction in the plane of scattering, outgoing after scattering process, Fig. 3
${𝒅}_{\text{E}}$	East direction
$𝒏$	normal to the plane of incidence
${𝒏}^{\prime }$	normal to the plane of departure
$I$	intensity
$I^{\prime }$	scattered intensity
$j$	anisotropic emissivity
$\Phi$	scattering phase function
$C_{\text{ext}}$	extinction cross section
$C_{\text{sca}}$	scattering cross section
$C_{\text{pol}}$	dichroism cross section
$C_{\text{cpol}}$	circular polarization cross section
$C_{\text{sca,pol}}$	polarization scattering cross section
$\Lambda$	albedo
$n$	density
$\tau$	optical depth
${\tau }_{\text{out}}$	optical depth from present position to outer cloud boundary
$𝑺$	Stokes vector
$(I,Q,U,V)$	components of Stokes vector $𝑺$
$𝑴$	$4\times 4$ Müller matrices
$𝑹$	rotation matrix
${𝑹}_{\text{obs}}$	rotation from cell to observer frame
$𝒁(𝒌,{𝒌}^{\prime })$	scattering matrix
$𝑲$	extinction matrix
${𝑬}_j$	dichroism and birefringence matrix

Equations

$$\begin{gathered} \frac{\text{d}I}{\text{d}s}(\mathbf{r},\mathbf{k})=j(\mathbf{r},\mathbf{k})-n(\mathbf{r})C_{\text{ext}}(\mathbf{r})I(\mathbf{r},\mathbf{k}) \\ +n(\mathbf{r})C_{\text{sca}}(\mathbf{r})\int \Phi (\mathbf{r},\mathbf{k},{\mathbf{k}}^{\prime })I(\mathbf{r},{\mathbf{k}}^{\prime })\text{d}{\Omega }^{\prime } \end{gathered}$$

$$I⟹\mathbf{S}=\left(\begin{array}{c}I\\ Q\\ U\\ V\end{array}\right).$$

$${\mathbf{S}}^{\prime }=\mathbf{M}\mathbf{S}=\left(\begin{array}{cccc}{M}_{11}& M_{12}& M_{13}& M_{14}\\ M_{21}& M_{22}& M_{23}& M_{24}\\ M_{31}& M_{32}& M_{33}& M_{34}\\ M_{41}& M_{42}& M_{43}& M_{44}\end{array}\right)\left(\begin{array}{c}I\\ Q\\ U\\ V\end{array}\right).$$

$${\mathbf{S}}^{\prime }=\mathbf{R}(\beta )\mathbf{S}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos 2\beta & \sin 2\beta & 0\\ 0& -\sin 2\beta & \cos 2\beta & 0\\ 0& 0& 0& 1\end{array}\right)\mathbf{S}.$$

$$\begin{gathered} n(\mathbf{r})C_{\text{sca}}(\mathbf{r})\int \Phi (\mathbf{r},\mathbf{k},{\mathbf{k}}^{\prime })I(\mathbf{r},{\mathbf{k}}^{\prime })\text{d}{\Omega }^{\prime } \\ ⟹n(\mathbf{r})\int \mathbf{Z}(\mathbf{r},\mathbf{k},{\mathbf{k}}^{\prime })\mathbf{S}(\mathbf{r},{\mathbf{k}}^{\prime })\text{d}{\Omega }^{\prime } \end{gathered}$$

$${\tilde{C}}_{\text{sca}}$$

$$=\frac{1}{I}\int \left(Z_{11}I+Z_{12}Q+Z_{13}U+Z_{14}V\right)\text{d}{\Omega }^{\prime }$$

$$=\int Z_{11}\text{d}{\Omega }^{\prime }+\frac{Q}{I}\int Z_{12}\text{d}{\Omega }^{\prime }$$

$$=C_{\text{sca}}+C_{\text{sca,pol}}\frac{Q}{I}.$$

$$-n(\mathbf{r})C_{\text{ext}}(\mathbf{r})I(\mathbf{r},\mathbf{k})⟹-n(\mathbf{r})\mathbf{K}(\mathbf{r},\mathbf{k})\mathbf{S}(\mathbf{r},\mathbf{k})$$

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

$$\begin{gathered} \frac{\text{d}\mathbf{S}}{\text{d}s}(\mathbf{r},\mathbf{k})=j(\mathbf{r},\mathbf{k})-n(\mathbf{r})\mathbf{K}(\mathbf{r},\mathbf{k})\mathbf{S}(\mathbf{r},\mathbf{k}) \\ +n(\mathbf{r})\int \mathbf{Z}(\mathbf{r},\mathbf{k},{\mathbf{k}}^{\prime })\mathbf{S}(\mathbf{r},{\mathbf{k}}^{\prime })\text{d}{\Omega }^{\prime }. \end{gathered}$$

$${\mathbf{d}}_{\text{E}}=\mathbf{k}\times {\mathbf{d}}_{\text{N}}\mathrm{and}{\mathbf{d}}_{\text{N}}={\mathbf{d}}_{\text{E}}\times \mathbf{k}.$$

$$\mathbf{n}=\frac{\mathbf{k}\times \mathbf{o}}{|\mathbf{k}\times \mathbf{o}|}$$

$$\cos \beta$$

$$\sin \beta$$

$${\mathbf{S}}^{\prime }=\mathbf{R}(\beta )\mathbf{S}$$

$$\frac{\text{d}\mathbf{S}}{\text{d}s}(\mathbf{r},\mathbf{k})=-n(\mathbf{r})\mathbf{K}(\mathbf{r},\mathbf{k})\mathbf{S}(\mathbf{r},\mathbf{k}).$$

$$\mathbf{S}(s)={\text{e}}^{-C_{\text{ext}}ns}\left(\begin{array}{c}{I}_0\cosh (C_{\text{pol}}ns)-Q_0\sinh (C_{\text{pol}}ns)\\ Q_0\cosh (C_{\text{pol}}ns)-I_0\sinh (C_{\text{pol}}ns)\\ U_0\cos (C_{\text{cpol}}ns)-V_0\sin (C_{\text{cpol}}ns)\\ V_0\cos (C_{\text{cpol}}ns)+U_0\sin (C_{\text{cpol}}ns)\end{array}\right).$$

$$\Delta s=\frac{-\ln (\xi )}{n{C}_{\text{ext}}}.$$

$$\tau (\Delta s)=C_{\text{ext}}n\Delta s.$$

$${\tau }^{\mathrm{pol}}(s)$$

$$=C_{\text{ext}}ns-\ln \left(\cosh (C_{\text{pol}}ns)-\sinh (C_{\text{pol}}ns)\frac{Q_0}{I_0}\right)$$

$$\frac{\Delta s}{\Delta s^{\mathrm{pol}}}\approx \frac{\tau (\Delta s)}{{\tau }^{\mathrm{pol}}(\Delta s)}.$$

$$\Delta s^{\mathrm{pol}}\approx \Delta s\frac{{\tau }^{\mathrm{pol}}(\Delta s)}{-\ln (\xi )}.$$

$$\Lambda =\frac{{\tilde{C}}_{\text{sca}}}{{\tilde{C}}_{\text{ext}}}=\frac{C_{\text{sca}}+C_{\text{sca,pol}}\cdot Q/I}{C_{\text{ext}}+C_{\text{pol}}\cdot Q/I}$$

$$p(\theta ,\phi )=\frac{I^{\prime }(\theta ,\phi )\sin \theta }{\int I^{\prime }({\theta }^{\prime },{\phi }^{\prime })\sin {\theta }^{\prime }\text{d}{\theta }^{\prime }\text{d}{\phi }^{\prime }}$$

$$\begin{gathered} I^{\prime }(\theta ,\phi )=Z_{11}(\theta ,\phi )I+Z_{12}(\theta ,\phi )Q+Z_{13}(\theta ,\phi )U \\ +Z_{14}(\theta ,\phi )V. \end{gathered}$$

$$\phi$$

$$\theta$$

$$v_{ceil}$$

$$=\max \left[\left(Z_{11}I+Z_{12}Q+Z_{13}U+Z_{14}V\right)\sin \theta \right]$$

$$\zeta ={\xi }_3\cdot v_{ceil}$$

$$\zeta \le I^{\prime }(\theta ,\phi )\sin \theta$$

Full text

1 Introduction

Dust processes the radiation that astronomical objects emit. Especially ultraviolet to optical light is easily scattered or absorbed. The dust also emits radiation, mostly at infrared to sub-millimetre wavelengths. Nearly all astronomical objects are seen through dust shrouds and their spectral energy distributions (SED) are altered by this. Dust also polarizes the radiation that is scattered, extinct or emitted by grains while passing through the medium. Observations have shown the polarization of radiation due to dust, for example around active galactic nuclei (Miller & Goodrich 1990), around supernovae (Tran et al. 1997), around single stars (Forrest et al. 1975), in the galactic interstellar medium (Serkowski et al. 1975), or other galaxies (Montgomery & Clemens 2014).

Dust clouds have complex morphologies and varying density profiles, and they contain non-spherical dust grains of various sizes and compositions, which are partially aligned. An analytic description of the signatures of dust on the radiation field is only possible under the assumption of strong simplifications, which generally do not hold. A common numerical technique to account for all these different dependencies of the dust–photon interactions is the Monte Carlo radiative transfer (MCRT) approach. In this procedure, photon packages propagate probabilistic through a simulation of the volume under study. Effects like absorption, scattering, and emission of photons are explicitly carried out and the effects on the dust temperature are recorded. A vast body of theoretical work has been developed for MCRT, a starting point can be the review by Steinacker et al. (2013). A large number of radiative transfer codes based on the Monte Carlo technique are available, and their scopes, sophistication, and application domains are extremely different.

The polarization of radiation influences how it interacts with the dust and the interaction with the dust influences the polarization. Several MCRT codes, therefore, started treating polarization in simplified schemes. A first step is often taken by calculating the polarization which is due to scattering at spherical dust grains and electrons. MCRT codes that include such polarization mechanisms have been presented by Voshchinnikov & Karjukin (1994), Bianchi et al. (1996), Harries (2000, TORUS), Watson & Henney (2001, Pinball), Pinte et al. (2006, MCFOST), Min et al. (2009, MCMAX), Robitaille (2011, HYPERION), Goosmann et al. (2014, STOKES), and Kataoka et al. (2015, RADMC-3D), to name a few. We have also treated polarization due to scattering on dust as presented in Peest et al. (2017), hereafter called P17.

Very few MCRT codes implement additional and more sophisticated processes that can lead to the polarization of radiation by dust. Among them are codes developed by Whitney & Wolff (2002) and Lucas (2003), which treat scattering, dichroic extinction and birefringence due to perfectly aligned spheroidal grains. The POLARIS code (Reissl et al. 2016) calculates the polarized emission, scattering, dichroic extinction and birefringence due to (imperfectly) aligned oblate grains. The MCRT code by Bertrang & Wolf (2017) uses in the first step spherical grains for the dust heating and scattering processes and then it uses non-spherical grains aligned by radiative torques and magnetic fields for the dust emission phase. MoCafe (Lee et al. 2008; Seon 2018) uses spherical grains including polarization by scattering and an empirical formula to emulate dichroism based on optical depth and magnetic field alignment. Vandenbroucke et al. (2021) implemented polarized emission by partially aligned spheroidal grains in the SKIRT MCRT code (Camps & Baes 2020).

A code that implements or improves an established functionality can be verified by comparing its results with previous numerical computations. Such benchmark tests are available for dust radiative transfer (RT) models of unpolarized light by Ivezic et al. (1997), Pascucci et al. (2004), Pinte et al. (2009), and Gordon et al. (2017). They provide numerical proof for the scattering, extinction, and emission of radiation due to dust in various environments. For treatments of dust polarization, such tests are generally not accessible except one, which covers the case of scattering by spherical dust grains (Pinte et al. 2009). It compares polarization images of a flared dust disk around a central star. The dust is spherical and extends optically thin above and below the optically thick disk. The geometry of the test provides mostly the accuracy of the polarization after single scattering events. As discussed in P17 when various of the above MCRT codes are applied to more general cases they show significant differences without the correct results being known a priori. Flexible and easily reproducible tests are necessary to benchmark and verify the current and future implementations of polarization due to non-spherical dust.

In this paper, we present a simple and efficient implementation of the polarization of radiation due to scattering, dichroic and birefringent extinction by spheroidal particles. We provide analytical test cases against which we verify our implementation and estimate the numerical accuracy of our code, which we call hereafter MCpol. The goal of this paper is to provide test cases to other teams enabling them to verify and estimate the numerical precision of their codes. This should allow groups from many different research areas to explore polarization including estimated numerical uncertainty.

In Sec. 2 we present the basic equations governing dichroic extinction and scattering of radiation by spheroidal dust grains. We illustrate the functionality of MCpol, and detail how we implement these polarization mechanisms (Sec. 3). We describe the validation methods in Sec. 4, which are applied to confirm that the implementations work as desired. We discuss our results and present our conclusions in Sec. 6.

2 Radiative transfer with polarization

2.1 Stokes vector and Müller matrices

The RT equation describes the interaction of radiation with matter. For a medium that absorbs, scatters, and emits radiation, the basic RT equation reads

[TABLE]

with $I$ the specific intensity of the radiation field111In our notation we do not write obvious wavelength dependencies, e.g., for the intensity and the cross sections., $n$ the matter density, $j$ the anisotropic emissivity of radiation, $C_{\text{ext}}$ the extinction cross section, $C_{\text{sca}}$ the scattering cross section, and $\Phi$ the phase function. The different quantities are depending on the position in the medium ${\boldsymbol{r}}$, and the direction of the radiation ${\boldsymbol{k}}$.

This description is incomplete, as it does not consider the polarization of the radiation. We use the Stokes formalism instead, in which the radiation field is the 4D Stokes vector ${\boldsymbol{S}}$,

[TABLE]

The first element describes the intensity of the radiation, the second and third the linear polarization and the fourth the circular polarization (Stokes 1852). There are different conventions in the literature concerning the handedness (Hamaker & Bregman 1996; Peest et al. 2017), we use the convention favored by the IAU (Contopoulos & Jappel 1974), which for example is not applied by Planck Collaboration et al. (2015).

Changes to the Stokes vector are described by $4\times 4$ Müller matrices ${\boldsymbol{M}}$,

[TABLE]

Linear polarization refers to a particular direction, which can be any direction perpendicular to the propagation direction. In analogy to observations, we call our choice of the reference direction North, ${\boldsymbol{d}}_{\text{N}}$, which shall not be confused with magnetic North. In our definition, North is “up” when looking against the propagation direction towards the source. In the plane of the sky East is counted “left” from the North as is common in astronomy. The Stokes vector changes depending on the choice of North. When rotating the North direction by an angle $\beta$, the Stokes vector must be multiplied with a rotation matrix ${\boldsymbol{R}}(\beta)$

[TABLE]

The ideal choice of reference direction depends on the phenomenon.

2.2 Scattering by spheroidal grains

For radiation that is scattered on spherical grains, the North is usually chosen in the plane of scattering, defined by the propagation direction before scattering ${\boldsymbol{k}}$ and the propagation direction of the photons after scattering ${\boldsymbol{k}}^{\prime}$ (see e.g., Chandrasekhar 1960). This allows for more efficient approaches when treating scattering events as the 16 entries of the scattering matrix (see below) can be simplified in the case of spherical grains to just four independent elements. The method is widely applied (e.g., Goosmann & Gaskell 2007).

For spheroidal grains, the Müller matrices typically have their simplest form when the North is in the plane of incidence (Mishchenko et al. 2002), defined by the propagation direction before scattering and the grain symmetry axis (Fig. 3).

Scattering events by spheroidal grains are described by the $4\times 4$ Müller matrix denoted as ${\boldsymbol{Z}}$ and called the scattering matrix. For each incoming photon ${\boldsymbol{k}}$, ${\boldsymbol{Z}}({\boldsymbol{k}},{\boldsymbol{k}}^{\prime})$ provides the probability that it is scattered towards an outgoing direction ${\boldsymbol{k}}^{\prime}$. The matrix elements depend on the grain shape, size, porosity, and optical constants, as well as the wavelength, and direction of incoming and outgoing radiation. The calculation of ${\boldsymbol{Z}}$ can be simplified using symmetry arguments. in the case of spheroidal grains in ${\boldsymbol{Z}}$ only seven of the 16 elements are independent (Bohren & Huffman 1998; Abhyankar & Fymat 1969). For spheroids, the separation of variables method has been established (Asano & Yamamoto 1975; Voshchinnikov & Farafonov 1993). For rotationally symmetric particles, the scattering matrix can be calculated using the so-called T-matrix method (Mishchenko et al. 1996; Mishchenko 1991; Vandenbroucke et al. 2020). For general non-spherical particles, ${\boldsymbol{Z}}$ can be calculated by binning the grains into a grid of dipoles (Purcell & Pennypacker 1973; Draine 1988; Draine & Flatau 1994). Several codes are publicly available or can be requested from the above-listed authors. The different methods provide consistent results in computing the cross sections up to size parameter $x=2\pi\/a/\lambda\sim 10$ and encounter numerical problems at x\mathbin{\lower 3.0pt\hbox{\hbox to0.0pt{\raise 5.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax}}20 (Siebenmorgen & Peest 2019; Draine & Hensley 2021).

Because of the rotational symmetry of spheroids, the direction of the incoming radiation can be described by one angle, the angle of incidence. It is defined as the angle between the direction of propagation and the grain symmetry axis.

The scattering part of the RT equation (Eq. 1) changes to a tensor equation when polarization is considered

[TABLE]

Following Whitney & Wolff (2002) and Mishchenko et al. (2002) we can calculate the total scattering cross-section $\tilde{C}_{\text{sca}}$. It is defined as the integral over all the scattered intensity $I^{\prime}$ relative to the intensity $I$ incident onto a spheroidal grain,

[TABLE]

The integrals over $Z_{13}$ and $Z_{14}$ are zero, because the grains are rotationally symmetric and we integrate over the entire unit sphere (Van De Hulst 1957, p. 47–51). The “classical” scattering cross section for unpolarized radiation is $C_{\text{sca}}$, while $C_{\text{sca,pol}}$ is the polarization scattering cross-section, which complicates matters and becomes important for polarized impinging light. Hence, the probability that a photon is scattered by a spheroid depends on the polarization status of the photon.

2.3 Extinction by spheroidal grains

Extinction by spherical grains is very straightforward and fully described by a single quantity, the extinction cross-section $C_{\text{ext}}$. For spheroidal grains, things are complicated, as the cross-section also depends on the polarization status of the radiation. The relevant effects are called dichroism and birefringence. Dichroism means that the extinction differs for differently polarized radiation. Birefringence means that the travel speed of the photons through the medium depends on their polarization status.

The extinction term in the RT equation can be expanded to contain the dichroism and birefringence effects,

[TABLE]

with ${\boldsymbol{K}}$ the extinction matrix. For spheroids with North in the plane of incidence, the extinction matrix has a block diagonal shape (Martin 1974; Mishchenko 1991; Whitney & Wolff 2002; Lucas 2003; Krügel 2008; Voshchinnikov 2012),

[TABLE]

The “classical” extinction cross section is $C_{\text{ext}}$. The dichroism cross section $C_{\text{pol}}$ describes how much more a radiation wave is polarized parallel to the North direction is extincted than a radiation wave polarized parallel to the East direction. The birefringence cross-section $C_{\text{cpol}}$ represents how much more a North polarized wave is slowed down than an East polarized wave. All cross sections depend on the grain shape, size, porosity, and optical constants, the wavelength, and the angle of incidence.

2.4 RT equation

When polarization due to spheroids is considered, the basic RT equation (1) becomes the following partial integrodifferential vector equation

[TABLE]

The goal of polarized radiative transfer is to derive the Stokes vector ${\boldsymbol{S}}$ at any position ${\boldsymbol{r}}$ and in any direction ${\boldsymbol{k}}$ for a given density $n$, emissivity $j$ and optical properties and alignment of the dust grains.

3 Monte Carlo solution of polarized radiative transfer

3.1 MCpol

We implement polarization in the MCRT code developed by Krügel (2008). The efficiency of the code was significantly improved by Heymann & Siebenmorgen (2012) who vectorized it using graphical processing units (GPUs) and applied for optically thin cells the treatment by Lucy (1999) and for optically thick cells the method by Fleck & Canfield (1984). A ray-tracing routine allows the computation of images of areas of interest at any wavelength. It uses the scattering and absorption events treated in the simulation and enable the calculation of SEDs of subdomains of the model. The code has been applied to describe the effective extinction curve when photons are scattering in and out of the observing beam (Krügel 2009). It has been used to investigate the structure of disks around T Tauri and Herbig Ae stars (Siebenmorgen & Heymann 2012), to create a two-phase AGN SED library (Siebenmorgen & Heymann 2012), to study the effects of a clumpy ISM on the extinction curve (Scicluna & Siebenmorgen 2015), and to examine the the appearance of dusty filaments at different viewing angles (Chira et al. 2016). A time-dependent MCRT version was used to discuss the impact of a circumstellar dust halo on the photometry of supernovae Ia (Krügel 2015).

So far polarization has not been considered in MCpol. We present a MCRT dust polarization implementation for spheroidal grains that keeps the code backwards compatible. This means that the logical order of the processing steps remains as before and that all calculations required for the polarization are in a module separate from the main program. The most important change compared to the previous MCRT code is that we add the Stokes formalism. Each photon package in the original version of the code was characterized by its frequency, origin, and propagation direction; in the new version, the East direction ${\boldsymbol{d}}_{\text{E}}$ and the Stokes $Q$, $U$ and $V$ parameters are stored as well, as motivated in P17.

In the following subsections, we describe how the effect of scattering and extinction by spheroidal grains is implemented in MCpol. First the Stokes vector must be oriented so that North is in the plane of incidence for the aligned dust grains (Sec. 3.2). A photon package propagating through a cloud of spheroids continuously changes its Stokes vector (Sec. 3.3). The relationship between physical path length through a cloud of aligned spheroids and the optical depth experienced by radiation is discussed in Sec. 3.4. In case the radiation interacts with the dust, a part of the intensity is either scattered or absorbed following its albedo. The albedo of the dust depends on the polarization of the incoming radiation, which further complicates the problem (Sec. 3.5). If a scattering takes place, the propagation direction of the photon after the scattering is obtained via rejection sampling of the scattering Müller matrix (Sect. 3.6). The polarization of photons exiting the model space is recorded (Sec. 3.7). Additionally, we implement a routine to calculate the probability of directed scattering (peel-off, Sec. 3.8). Finally, a routine is presented that is capabable of treating dichroism for inverse ray-tracing, which can be used for the generation of polarization maps (Sec. 3.9).

3.2 Orienting the Stokes vector

The polarization of the radiation changes as it propagates through the model space. As described in Sec. 2, these changes are encoded in Müller matrices. When dealing with aligned spheroidal grains, the Müller matrices have their simplest form when the North direction of the radiation is in the plane of incidence, defined by the propagation direction of the radiation ${\boldsymbol{k}}$ and the symmetry axis of the grains ${\boldsymbol{o}}$.

The propagation direction ${\boldsymbol{k}}$, the North direction ${\boldsymbol{d}}_{N}$, and the East directions ${\boldsymbol{d}}_{\text{E}}$ are unit vectors and describe a right handed coordinate system (Fig. 3),

[TABLE]

During the life cycle of a photon package, the orientation of the plane of incidence alters frequently: it changes either when the orientation of the grain’s changes, e.g., when a photon package enters a new cell where the grains are oriented in a different way as in the previous cell or when the propagation direction changes e.g., after a scattering event or when a new photon gets emitted by the dust. When either of these situations occurs, we need to rotate the Stokes vector to ensure that the North direction is in the plane of incidence.

We compute the normal to the plane of incidence ${\boldsymbol{n}}$ from the propagation direction ${\boldsymbol{k}}$ and the grain symmetry axis ${\boldsymbol{o}}$,

[TABLE]

Let $\beta$ be the angle that describes the rotation so that the North direction is in the plane of incidence. Then $\beta$ is also the angle such that the East direction is perpendicular to the plane of incidence. Therefore $\beta$ is the angle between ${\boldsymbol{d}}_{\text{E}}$ and ${\boldsymbol{n}}$. We calculate $\beta$ from the two directions and the vector algebra relations,

[TABLE]

We then rotate the Stokes vector ${\boldsymbol{S}}$ to the plane of incidence using $\beta$,

[TABLE]

We would like to stress again that the plane of incidence and the angle of incidence depend on both ${\boldsymbol{k}}$ and ${\boldsymbol{o}}$. The steps above must be repeated whenever either of them changes, that is when the photon is scattered or enters a cell with a different grain orientation or when a new photon is emitted. Note that Eq. 11 breaks down and cannot be applied when the angle of incidence is 0∘ or 180∘, that is, when ${\boldsymbol{k}}$ and ${\boldsymbol{o}}$ are parallel or anti-parallel. In this case, the photon path is oriented along the spheroid’s symmetry axis, and there is no preferential direction for the North or East direction. In that case, any arbitrary direction perpendicular to ${\boldsymbol{k}}$ can be chosen as ${\boldsymbol{n}}$.

3.3 Propagating the photon package

The polarization of a photon changes when travelling through a dichroic or birefringent medium. To mathematically describe this, we consider a beam of light propagating through such a medium. Without dust emission or scattering, the polarized RT equation (9) simplifies to

[TABLE]

As we assume that the density, grain orientation, and the optical properties within a single dust cell are constant, this linear differential matrix equation can be solved analytically within each individual cell. Using Eq. (8) we obtain two coupled systems and their solution is (Lucas 2003; Whitney & Wolff 2002)

[TABLE]

The equation describes the change of the Stokes vector when the photon travels through a single cell in a dichroic or birefringent medium. The change of the Stokes vector ${\boldsymbol{S}}(s)$ needs to be considered in MCRT computation when the photon packets propagate by $s$ in such a medium. We note that, whenever the photon package enters a new cell where the grain orientation is different, we need to apply a re-orientation of the Stokes vector (Sec. 3.2) before we can continue the propagation.

3.4 Optical depth and path length

One necessary exercise in MC treatments is computing the optical depth along the flight path of the photon package through the cell for determining whether and where it will interact with the dust. In MCpol the grid cells are cubes with constant density $n$ and grains have a fixed alignment. The flight path of the photon packets go along a straight line from the entry point into a cube up to either the interaction point or in case of no interaction up to the exit point of that cube.

Photon packets of unpolarised light in a non-dichroic medium travel the distance from the entry point of a cube to the interaction point $\Delta s$ corresponding to an optical depth $\tau(\Delta s$). This optical depth can be computed from a random exponential distribution such that the interaction of the photon package with dust in a cell follows a uniform distributed random number $\xi$ (Krügel 2008; Witt 1977; Lucy 1999; Steinacker et al. 2013),

[TABLE]

The optical depth along $\Delta s$ is

[TABLE]

In comparison to the unpolarized light, the photon packets of polarised light in a dichroic medium travel the distance from the entry point of a cube to the interaction point $\Delta s^{\rm pol}$ and corresponding optical depth $\tau^{\rm pol}(\Delta s^{\rm pol}$). We add suffix $\rm{pol}$ for quantities that depend on polarization. Following Baes et al. (2019), a direct connection between the physical path length and the optical depth as given for unpolarised light is lost in media with dichroism and birefringence. Considering polarization in a dichroic medium the optical depth cannot be calculated in closed form as in Eq. (18). Diminishing of the polarized radiation needs to consider the change of the Stokes vector ${S}$ when the photon packet travels through that cell and follows Eq. (16). In that solution, there is besides the extinction, which is given by an exponential decay, these additional terms of the Stokes components noted between the brackets. We visualize in Fig. 1 that dichroism indeed complicates the relation between path length and optical depth. We choose two examples of the exact solution for polarized light of Eq. 16 adopting $q_{o}=Q_{o}/I_{o}=0.5$ and $q_{o}=-0.1$ and show the solution for unpolarized light given by $\exp\left(-\tau\right)$. For better visualization, an extremely dichroic cell has been adopted with $n=0.0015$, $C_{\text{ext}}=1000$, and $C_{\text{pol}}=300$ (a.u.), hence the extinction optical depth $\tau(1)=1.5$. Note the vast variations of the interaction points ${\Delta s}$ and ${\Delta s^{\rm{pol}}}$ for the same optical depth either $\tau$ or $\tau^{\rm{pol}}$ at same random number $\xi$. The interaction points are at $\Delta s^{\rm{pol}}=0.47$ for $q_{o}=0.5$, $\Delta s^{\rm{pol}}=0.57$ for $q_{o}=-0.1$ and $\Delta s=0.54$ for unpolarized light (Fig. 1).

To stress the fact that we deal with dichroic extinction, we use the term polarization optical depth $\tau^{\rm{pol}}$. We use the definition of an optical depth and insert Eq. (16)

[TABLE]

The right-hand side of Eq. 20 can, except for special cases, not be simplified. For non-spherical dust particles, there is another complication that the orientation of the Stokes vector to the grain orientation, hence the dust cross sections need to be considered. In MCpol the grain alignment in a cell is fixed, which simplifies computing the angle of incidence $\alpha$, which is the angle between the Stokes vector of the incoming photon and the grain orientation (Fig. 3). The Stokes vector of the incoming photon needs to be rotated to the plane of incidence and $C_{\text{ext}}(\alpha)$ and $C_{\text{pol}}(\alpha)$ are computed for that $\alpha$. For clarity, the angle and frequency dependence has been dropped in Eq. 20.

Whenever taking dichroism into account the polarization state changes along the flight path of the photon (Sect.3.3) and an analytical solution for $\Delta s^{\rm{pol}}$ does not exist. Baes et al. (2019) showed that the solution can be expressed by a Taylor expansion and that the extinction optical depth $\tau$, this is the optical depth ignoring polarization approximates the polarization optical depth $\tau^{\rm{pol}}$ to first order, hence

[TABLE]

The density $n$ of the cube, the angle of incidence $\alpha$, and $C_{\text{ext}}(\alpha)$ are known. As indicated by the horizontal line in Fig. 1, the same optical depth for unpolarised and polarised light are chosen by the same random number, ${-\ln(\xi)=\tau(\Delta s)=\tau^{\rm{pol}}(\Delta s^{\rm{pol}})}$, which leads to different path lengths $\Delta s\neq\Delta s^{\rm{pol}}$. The distance ${\Delta s}$ is derived inserting the quantities $n$, $C_{\text{ext}}$, and $\tau(\Delta s)=-\ln(\xi)$ into Eq. 17. This allows computing ${\tau^{\rm{pol}}(\Delta s)}$ using Eq. 20. Hence, we derive the path length of the polarised light

[TABLE]

In the magnified examples of Fig. 1, the approximation of $\Delta s^{\rm{pol}}$ (Eq. 22) is within $0.3$ % of the correct solution. For typical ISM dust $C_{\text{pol}}$ is about a factor 100 smaller (Siebenmorgen 2022) and the uncertainty in $\Delta s^{\rm{pol}}$ is below $10^{-4}$. Whenever possible, the grid in MCpol is set up so that the total extinction optical depth of a cube, or when needed sub-cube, is about $\tau\la 1$. Radiation penetrating through a highly optically thick medium is challenging for MCRT codes, e.g., Min et al. (2009); Siebenmorgen & Heymann (2012); Gordon et al. (2017); Camps & Baes (2015). MC applications that consider polarization treatments at path length $\tau^{\rm{pol}}\gg 1$ can progressively solve the radiative transfer equation along the ray in each cube as outlined by Baes et al. (2019). This goes in hand with extraordinary boost in computing time and is not considered.

3.5 Dust albedo of polarized light

When the radiation interacts with a dust grain, it is either scattered or absorbed. The ratio of scattered over extincted intensity defines the albedo $\Lambda$. We rotate the Stokes vector so that its component $U=0$, apply Eqs. (6) - (8), and derive the albedo of a spheroid including polarization,

[TABLE]

This equation can be seen as a more general form of the albedo. For unpolarized radiation or if the polarization scattering cross-section and the dichroism cross-section are zero, Eq. (23) reduces to its common form. We also note that Eq. (23) is independent of $C_{\text{cpol}}$. This is consistent with the physical meaning of $C_{\text{cpol}}$: It describes a differential phase shift, a delay, between the North and East polarized parts of the radiation. Such a delay does not reduce the intensity of the radiation.

3.6 Sampling the scattering Müller matrix

In the case of non-polarized MCRT, simulating a scattering event is relatively straightforward. It essentially comes down to generating a random scattering angle $\theta$ from the scattering phase function $\Phi(\theta)$. Based on the original propagation direction ${\boldsymbol{k}}$ of the photon package, this scattering angle can be converted to a new propagation direction ${\boldsymbol{k}}^{\prime}$, and this describes the scattering event. In polarized MCRT, computing a scattering event is more complicated: not only is the random generation of a new propagation direction more complex, but we also must update the polarization status and the reference direction of the photon package.

Rather than a single scattering angle, we must determine a new pair of angles $(\theta,\varphi)$. In this pair, $\theta$ is still the scattering angle, that is the angle between the original propagation direction ${\boldsymbol{k}}$ and the new propagation direction ${\boldsymbol{k}}^{\prime}$. The angle $\varphi$ is the azimuth of ${\boldsymbol{k}}^{\prime}$ in a coordinate system with ${\boldsymbol{k}}$ as polar direction and the North direction as the reference direction with azimuth zero.

The appropriate probability distribution from which the random couple $(\theta,\varphi)$ needs to generated is

[TABLE]

where

[TABLE]

This formula shows two characteristics. Firstly, the probability density function is a true bi-variate distribution that depends explicitly on $\theta$ and $\varphi$ in a nontrivial way. Secondly, it does not only depend on the characteristics of dust grains and the grain alignment, but also on the polarization state of the photon package. To generate a random couple $(\theta,\varphi)$ from the bi-variate probability distribution function (24) we using the rejection sampling method (Von Neumann 1951). This method is comparatively simple and is readily applicable to multivariate probability density functions (Devroye 2013; Baes & Camps 2015). A scattering angle pair ($\theta$, $\varphi$) is generated and accept if another third random number $\xi_{3}$ is lower than the probability $P(\theta,\varphi)$ of scattering in this particular direction (Fig. 2). A ceiling value $v_{ceil}$ gives the maximum probability for any of the angles and is used to scale the third random number. The actual implementation begins with calculating two scattering angles based on uniform deviates $\xi_{i}$,

[TABLE]

It is important to note here, that we sample uniform in $\theta$, which means that we sample more angles per surface area in the forward and backward directions. This is motivated by the fact that dust grains scatter preferably forward or backwards. The sampling density must be taken into account when we calculate the probability of scattering in the direction $\theta$, $\varphi$. We calculate a ceiling value $v_{ceil}$, which is representative for the highest probability of scattering towards a direction $\varphi$, $\theta$. Following Eq. (3) with the scattering matrix ${\boldsymbol{Z}}$ the outgoing intensity $I^{\prime}$ for an incoming photon with the Stokes vector ${\boldsymbol{S}}$ is

[TABLE]

with $Z_{ij}(\lambda,\alpha,\varphi,\theta)$ the elements of the scattering matrix. The factor $\sin{\theta}$ compensates for oversampling the forward and backward regions and reducing the ceiling value $v_{\text{ceil}}$. The ceiling value needs to be recomputed for each scattering event, when the Stokes parameters $Q$, $U$, and $V$, the angle of incidence $\alpha$, or the wavelength $\lambda$ change. The decision of whether the angle pair is accepted is based on the third random number,

[TABLE]

with $0\leq\xi_{3}\leq 1$. The angles $\varphi$ and $\theta$ are accepted, if

[TABLE]

A low ceiling value $v_{ceil}$, will increase the probability of accepting an angle pair. Our method of changing the sampling density is a simple approach to reduce the average number of draws until a pair is accepted. There can be more sophisticated methods, using the scattering behavior of the dust mixture used in the simulations.

After generating a random angle pair $(\theta,\varphi)$, we need to update the characteristics of the photon package. The obvious characteristic to update is the propagation direction. The new propagation direction after scattering is calculated following Eq. (30) of P17,

[TABLE]

The polarization state of the photon package also needs to be updated since the scattering event does affect the state. The new Stokes vector becomes

[TABLE]

where the factor between the brackets guarantees that the total intensity of the photon package is conserved during the scattering event. Finally, also the East direction is updated. The convention we use is that the new North direction is in the plane of departure after scattering. This plane is given by the direction after scattering ${\boldsymbol{k}}^{\prime}$ and the symmetry axis of the grain ${\boldsymbol{o}}$, similar to the plane of incidence. Therefore, the East direction after scattering is

[TABLE]

If the outgoing direction and the symmetry axis are parallel or anti-parallel, the photon package continues its path through the model space without the need to rotate the Stokes vector.

3.7 Detection of escaped photons

Eventually, photon packages will leave the model space. We record their exiting directions ${\boldsymbol{k}}$ and wavelengths. Photons that exit with the same wavelength and with similar angles to the simulation $z$-axis are binned together. For axisymmetric geometries, this method allows for the quick calculation of the spectral energy distribution of the model under different viewing angles. The Stokes vector of the photon packages is rotated upon leaving, such that North is in the plane given by the escape direction and the $z$-axis of the model space. This permits binning photon packages by adding their Stokes vectors and intensities.

3.8 Directed scattering (peel-off)

MCRT codes commonly allow viewing the model space from (arbitrary) directions ${\boldsymbol{k}}_{\text{obs}}$. This simulates an observation by a distant observer. The chance of a photon scattering directly towards the observer is low. We, therefore, employ the peel-off method (Yusef-Zadeh et al. 1984) in which the probability of scattering towards the observer is calculated explicitly. During the model run, we store for all scattering events the position, direction ${\boldsymbol{k}}$ and Stokes vector ${\boldsymbol{S}}$ of the photon packages before scattering. After the simulation finishes we use the scattering matrix ${\boldsymbol{Z}}$ to calculate the Stokes vector ${\boldsymbol{S}}_{\text{obs}}$ of the radiation that would have scattered towards ${\boldsymbol{k}}_{\text{obs}}$,

[TABLE]

where $C_{\text{sca}}^{-1}$ is used as a normalization factor (see below), and $\varphi_{\text{obs}}$ and $\theta_{\text{obs}}$ the angles by which the photon is scattered towards the observer. The scattering angle $\theta_{\text{obs}}$ is calculated from the direction before and after scattering,

[TABLE]

and the azimuth $\varphi_{\text{obs}}$ from the East direction ${\boldsymbol{d}}_{\text{E}}$ and the normal to the scattering plane calculated from ${\boldsymbol{k}}$ and ${\boldsymbol{k}}_{\text{obs}}$,

[TABLE]

The normalization $C_{\text{sca}}^{-1}$ of the scattering matrix is essential for calculating the peel-off probability using Eq. (35). The scattering matrix is stored internally as a multi-dimensional array, with $N_{\theta}$ and $N_{\varphi}$ elements along the respective angles $(\theta,\varphi)$. Following Eq. (6), the sum of the $Z_{11}$ elements over the unity sphere is $N_{\theta}\ N_{\varphi}\ C_{\text{sca}}$. The scattering probability has already been considered during the MC run and so one divides by $C_{\text{sca}}$.

The intensity of the radiation that is scattered towards the observer is then the first component of the Stokes vector resulting from Eq. (35). The intensity of the radiation that would reach the observer has to be reduced to account for (dichroic) extinction between the point of scattering and the observer. For applying this properly we developed an inverse ray-tracing routine discussed in the following section.

3.9 Inverse ray-tracing

An inverse ray tracer is developed for calculating spatially resolved maps of the model space. For this, it is necessary to know the optical depth of a scattering event towards the observer. We compute a Stokes map by sending rays in the direction $-{\boldsymbol{k}}_{\text{obs}}$ from the observer to an area of the model. Along the path of a ray, it encounters cells that are numbered by $i=1,...,m$. We calculate for each cell $i$ the amount of radiation that would be scattered and emitted towards ${\boldsymbol{k}}_{\text{obs}}$ and attenuate this according to the optical depth $\tau^{\text{out}}$ from the present position to the entry point of the ray in the model space. The map is created by varying the position at which the ray enters the model and under the assumption that the complete model space is sampled (Heymann & Siebenmorgen 2012).

Without considering polarization, the intensity of the radiation leaving the simulation is determined by the optical depth to the edge of the model space, $\tau^{\text{out}}$. By stepping through the system along $-{\boldsymbol{k}}_{\text{obs}}$ the optical depth increases. The ray crosses the cells $1,...,m$ and in each cell $i$, the optical depth to the edge increases by the product of the extinction cross-section of the dust $C_{\text{ext},i}$, its density $n_{i}$ and the path length $(\Delta s)_{i}$,

[TABLE]

The radiation that is scattered or emitted in cell $i$ towards the observer with an intensity $I_{i}$, will exit the model and reach the observer with the reduced intensity $I_{i}\,\textup{e}^{-\tau^{\text{out}}_{i}}$.

When we take dichroism into account, the optical depth through cell $i$ depends on the polarization of the radiation. In addition, the polarization of the radiation will change along ${\boldsymbol{k}}_{\text{obs}}$. Both effects are described by Eq. (16). One also has to consider the orientation of the grains, as they can be different for each cell. Consider the change of the Stokes vector ${\boldsymbol{S}}$ of the photon package on its way from cell $i$ along direction ${\boldsymbol{k}}_{\text{obs}}$ out of the model space. It is given by an initial rotation into the frame of the cell, ${\boldsymbol{R}}$, and then an alternating application of a dichroic extinction step towards the edge of that cell, followed by a rotation into the plane of incidence of the next cell, and so forth until it eventually leaves the model space. This can be written as a single equation by combining Eq. (16) and Eq. (4),

[TABLE]

For each cell $j$ with $0\leq j\leq i$, the rotation into the cell $j-1$ is ${\boldsymbol{R}}_{j}$, and ${\boldsymbol{R}}_{\text{obs}}$ is the rotation of North into the reference frame of the observer. The optical depth inside $j$ is $\tau_{\text{ext},j}=C_{\text{ext},j}\,n_{j}\,(\Delta s)_{j}$ and the dichroism and birefringence matrix is ${\boldsymbol{E}}_{j}$,

[TABLE]

where for cell $j$ the dichroism and birefringent optical depths are given by

[TABLE]

The ray-tracing equation, Eq. (40), assumes that the photon package is emitted or scattered at the edge between the cell $i$ and cell $i+1$. This is correct for optically thin cells ($\tau<0.1$). In the case of optically thicker cells, the emission is distributed along the path in cell $i$ and needs to be integrated within cell $i$.

The advantage of Eq. (40) is that it can be evaluated while stepping along the pencil beam. The product of the matrices of the previous cells is sufficient to calculate the Stokes vector of the next cell. For cell $i$, the matrix of the previous cell $i-1$ is multiplied by ${\boldsymbol{R}}_{i}\,\textup{e}^{-\tau_{\text{ext},i}}\,{\boldsymbol{E}}_{i}$. One can store the compound matrix that is updated when entering the next cell. The compounded exponential factor of Eq. (40) can be stored separately to keep the entries of the matrix closer to unity and to prevent numerical instabilities.

4 Validation

MCRT codes need to be carefully validated. For doing this, one can often use benchmark results from existing codes for comparison. However, as scattering, dichroism, and birefringence due to spheroidal dust grains are uncommon capabilities for MCRT codes, there are no such benchmarks to reproduce. Ideally, the different functionalities are tested individually. Such a procedure simplifies the identification of potential shortcomings or even mistakes and it enables those codes with different sets of functionalities may reuse the same tests.

As advocated by P17, it is particularly advantageous to compare the results of our implementations to analytical solutions. Analytical solutions are easy to reproduce and can be used to estimate the errors of the numerical treatment. Analytical solutions are also of interest to other teams aiming to verify their MCRT codes. In this Section, we develop analytical test cases to verify our numerical implementation of scattering, dichroism, and birefringence due to spheroidal dust. The MCRT treatment of scattering by spheroidal dust (Sec. 4) is confirmed using renewed versions of analytical test cases by P17 developed for spherical grains. Additional analytical test cases for estimating the numerical accuracy of MCRT codes that are treating dichroism and birefringence mechanisms are presented in Sec. 5.4.

In P17 we considered radiation scattering on spherical grains and developed several test cases to validate the numerical procedure for solving this radiative transfer problem of polarized light. The analytic solutions can only be expressed exactly because single scattering in a simplified geometry is considered (Fig. 4) and electrons with their uncomplicated Müller matrix is applied. Four test cases (TC) are distinguished with analytical solutions given by P17 (Eq. 44-63). Circular polarization is not considered in TC 1 - TC 3 ($V=0$) but is treated in TC 4.

TC 1: In the first scenario a central point source emits unpolarized radiation. For testing the peel-off scattering procedure we select two slabs of electrons. The photons scatter once at these optically and physically thin electron slabs, which lie in the $xy$ plane. A distant observer records the intensity, polarization degree and polarization angle of the radiation scattered off the slabs. TC 1 verifies the polarization by single scattering and the peel-off mechanism.

In the next test scenarios (TC 2 - TC 4), the central light source is replaced by a tiny cloud of electrons. They are illuminated by a collimated beam. The analytical solutions of P17 were derived for photons that encounter a first scattering event at the center and a second scattering event at the slabs (Fig. 4).

TC 2: In the second test case the direction of the collimated beam is in the same plane as the direction of the slabs. The second scattering is from the electron slab to the observer. In this configuration, scattering, the peel-off mechanism, as well as a random walk step of the photons are tested.

TC 3: In the third test case, the initial beam direction is at an angle to the plane of the slabs. The plane of scattering of the initial scattering is rotated to the plane of scattering off the slabs. The rotation leads to a variable polarization angle at the observer.

TC 4: In the final test case, the scattering properties of the particles are changed. The Müller matrix of the electrons are hypothetically changed so that radiation scattering on them can become circular polarized, hence $M_{34}\neq 0$ and $M_{43}\neq 0$ (Eq. 3), which are for electrons zero otherwise.

We expand the test scenarios of P17, to cover the case of scattering at spheroidal dust particles. For spheres and sphere-like particles, the scattering matrix is most simple when North is in the scattering plane before and after the scattering event. In that case, the scattering matrix ${\boldsymbol{Z}}$ reduces to a block diagonal shape, which depends for electrons only on the scattering angle $\theta$ and is independent of $\varphi$, hence

[TABLE]

The geometry and our notation of the scattering process on spheroidal grains is shown in Fig. 3 and summarized in Table 2. In contrast to spheres, for spheroidal particles the scattering matrix is usually given for North in the plane of incidence before the scattering event and North in the plane of departure after the scattering event.

In the test cases for spheroidal grains, we use spherical grains which are treated in the computations as if they were spheroids. Therefore, an orientation ${\boldsymbol{o}}$ is artificially assigned to these spherical grains. The orientation ${\boldsymbol{o}}$ can be chosen either fixed or even at random without altering the result of the scattering computation. In addition, the scattering matrix of the spheres is multiplied by two rotation matrices ${\boldsymbol{R}}$ to account for the different orientations of the North direction during the scattering process on spheroidal grains. These two rotations describe the change of the North direction from the plane of incidence to the plane of scattering by the angle $\varphi$ and the rotation from the plane of scattering to the plane of departure by the angle $\gamma$ (Fig. 3). The scattering matrix ${\boldsymbol{Z}}_{\text{sph}}$ of the so artificially assigned orientation of spherical particles, therefore,

[TABLE]

where $x=\cos 2\gamma$, $y=\sin 2\gamma$, $p=\cos 2\varphi$, and $q=\sin 2\varphi$. The angle $\gamma$ is calculated using the angle of incidence $\alpha$ and the scattering angles $\theta$ and $\varphi$ as described in Appendix A.

5 Results

5.1 Numerical set-up

We verified the numerical treatment of the analytical test cases TC 1 - TC 4 as implemented in MCpol, while P17 verified them using SKIRT. The MC treatments of both codes differ in many aspects. In SKIRT various acceleration methods such as forced scattering and forced absorption are included that are not implemented in MCpol. SKIRT (Camps & Baes 2015) uses a simplified Stokes formalism, which is derived assuming spherical particles; whereas MCpol treats spheroidal-shaped particles. SKIRT and MCpol use different vectorization technologies. The model grids in both codes are different. P17 used for the test cases a grid of $(601,601,60)$ cuboids each having a length of $(0.003,0.003,3\times 10^{-7})$ in arbitrary units and for each cube that includes electrons an optical depth of $\tau=10^{-3}$. MCpol uses a cartesian coordinate system with cubes for which any of these cubes can be divided into sub-cubes. We find for the test cases a good MCpol set-up using $(587,587,3)$ cubes with a side length of each cube of one. The cube in the center and cubes of the slabs are filled by electrons at an optical depth of $\tau=0.5$, other cubes are free of electrons and have $\tau=0$. Hence, the probability of multiple scattering events in cells filled with electrons is unlikely but not zero. Photons that scatter a second time either at the center or at the slabs are ignored so that one can compare the numerical results with the analytical solutions. Therefore, the choice of $\tau$ shall be done with some care. A too-low value of $\tau$ reduces the scattering probability of the photons while an increase of $\tau$ enhances the chance of multiple scattering and will lead to photons that must be ignored.

The number of photon packets launched is unless otherwise stated $N_{\gamma}=2.5\times 10^{10}$. Results when decreasing or increasing the number of emitted photon packets and applying different sampling of the model space are discussed in Sect. 5.3. We repeat test cases TC 1 - TC 3 of P17 using the scattering matrix ${\boldsymbol{Z}}_{\text{sph}}$ and pretending that sphere-like grains have an orientation and therefore, need to be treated as spheroids. We choose the grain orientation along $+\hat{e}_{z}$.

5.2 Testcases TC 1 - TC 4

We compare the numerical results of MCpol for the test cases TC 1 - TC 4 against the analytical solutions as provided in Eq. 44 - 63 by P17. In Fig. 5 and Fig. 6 the analytic solutions are shown as black lines and the result of MCpol with red lines along the instrumental position $x$ of the detector (Fig. 3). P17 derived the scattering angle $\theta$ for a given detector position, which we also label in Fig. 5 and Fig. 6, by

[TABLE]

The intensity profiles in arbitrary units are shown in the top panels of (Fig. 5). At the bottom of each sub-panel, the absolute deviation from the analytic solution is shown with blue lines and the grey-shaded areas mark relative deviations with a typical error as labelled. The left column of Fig. 5 shows the results of the first test case (TC 1). The result of MCpol jitters around the analytic solution. The middle graph shows the linear polarization degree computed using Eq. 3 in P17. MCpol reproduces the analytic solution of TC 1 to better than $\pm 0.1\%_{\rm{rel}}$. We attribute residual deviations to sampling errors in the initial direction of the photons from the source and the finite size of the model grid. The analytical curves are calculated for an infinitesimal “height” of the slabs. In the numerical treatment, the photon packages that arrive at a given detector pixel have taken paths with a small but finite height. In the outer parts ($\|x\|>0.3$), these paths lead to a lower average polarization degree. In the inner parts ($\|x\|<0.3$), they lead to a higher average polarization degree. In the bottom graph, the polarization angle computed using Eq. 4 in P17 is compared to the analytic result which is $0{{}^{\circ}}$. Here MCpol deviates from the analytical result in the very central region ($\|x\|<0.1$). Such inaccuracies cannot be avoided and are expected due to the set-up of the numerical grid, which is the finite height of the slabs, and the amplification of the noise for polarization degrees close to zero.

The middle column of Fig. 5 shows the result of the second test case (TC 2). The intensity curve includes the noise of several per cent and is higher than in TC 1. This is because the chance of scattering at the central cell is lower than unity. Some photons do not scatter so fewer photons propagate towards the slabs. This leads to a lower photon count at the detector. The polarization degree for TC 2 is correct to $0.1\%_{\rm{rel}}$ (Fig. 5). As in the case of TC 1, we attribute the remaining differences to the resolution of the grid. The analytic solution of the polarization angle is zero and is shown together with the numerical solution by MCpol at the bottom of Fig. 5.

The right column of Fig. 5 shows the result of the third test case (TC 3). The intensity curve is plotted in the first row. The model results follow the analytic solution. The noise in the numerical solution is comparable to the noise found in TC 2 and with the same explanations. The polarization degree derived by MCpol is for TC 3 correct to well below $0.2\%_{\rm{rel}}$. As we noted in P17, the area around $0.6\leq\|x\|\leq 0.7$ is very difficult to treat at high precision. In this area, the scattering angle from the initial direction to the slabs changes only by $0.1^{\circ}$, whereas the polarization degree by $\sim 40\%_{\rm{abs}}$. The scattering matrix is tabulated for integer angles and linearly interpolated for scattering angles in between. This comparatively sparse sampling is sufficient to calculate the correct result. In comparison to P17, we even attain a higher precision. We attribute this to the fact that the scattering matrix for spheroids not only contains the scattering but also, the rotations of the plane of scattering. In P17 we calculated the effect of the plane rotations separately. Variations of the polarization angle are shown in the third row of Fig. 5. MCpol applied to this test case follows the analytic result to the high precision of $0.1\%_{\rm{rel}}$.

Scattering at electrons does not lead to circular polarization, as $d(\theta)=0$ in Eq. (44). In P17 we introduced a hypothetical particle for which

[TABLE]

otherwise, these test particles behave as electrons. Such particles lead to circular polarization. They are applied in TC 4, in which the geometry of TC 3 is used. In Fig. 6 the results of TC 4 are displayed for MCpol launching the large number of $N_{\gamma}=2.5\times 10^{11}$ photon packets. The intensity, the reduced Stokes parameters $Q/I$ and $U/I$, the linear polarization degree, the polarization angle and the circular polarization expressed as $V/I$ is shown. MCpol (red) reproduces the analytic (black) solutions of Stokes parameters of TC 4 at high precision and for the circular polarization to better than $0.03\%_{\rm{rel}}$.

5.3 Numerical precision

The precision of MCpol is exemplified for computations of circular polarized light by varying the number of sub-cubes $N_{\rm{sub}}$ of cells along the $x,y$, and $z$ coordinates and the number of emitted photons $N_{\gamma}$. The difference between the numerical model ${V}/{I}_{{\rm{{\sc{MCpol}}}}}$ and the analytical solution ${V}/{I}_{{\rm{ana}}}$ (P17) is computed for the reduced Stokes parameter for circular polarized light of TC 4

[TABLE]

The minimum, maximum and standard deviation $\sigma$ of that difference $\Delta$ are given for the parameter variations in Table 1. For the minimum number of emitted photons ($N_{\gamma}=2.5\times 10^{9}$) an increase of $N_{\rm{sub}}$ from 1 to 3 does not improve the $1\sigma$ noise in the model and $N_{\gamma}$ shall be increased. Launching 10 times more photons the noise in the model without sub-cubes ($N_{\rm{sub}}=1$) remains whereas it improves by a factor of 2 for the $N_{\rm{sub}}=3$ and almost by a factor of 4 for the $N_{\rm{sub}}=11$ model. In models without sub-cubes the location of the first scattering event, this is the blob of electrons in the center of the model space, is not well sampled. This causes a systematic error in backward scattering ($\theta\sim 180^{\rm o}$) at $x\sim 0$. The effect is reduced by increasing the number of sub-cubes and for $N_{\rm{sub}}=11$ decreases to $\Delta\ (x\sim 0)<10^{-3}$ as shown in the bottom panel of Fig. 6. Even for our choice of optically thin cubes, the detailed fine sampling impacts the precision that can be reached, and which cannot be improved by increasing simply the photon statistics (Table 1). Obviously, with an even finer grid and more photons, the precision of MCpol could be improved further at the expense of a boost in computer power. In our server environment, the latter model took (already) $\sim 200$ h. For reaching similar precisions both codes show similar run times.

5.4 Dichroism and birefringence

In addition to the polarization due to scattering, we have to confirm the implementations of polarization due to dichroism and birefringence. For this a spherical distribution of dust with density $n$ and radius $r$ around a central source are considered. The test is performed using hypothetical dust particles which dichroically absorb and slows radiation, but which do not scatter. This is to remove any side effects from the scattering implementation. Such hypothetical dust particles are chosen to be analytically simple,

[TABLE]

The angle $\alpha$ is the angle of incidence (Fig. 3) and the sine and cosine make it such that the transition is smooth for sight lines around $\alpha=0$. Consider initially right-handed circular polarized radiation, ${\boldsymbol{S}}=(1,0,0,1)$, that is travelling through the dust cloud. Its direction and the grain orientation have a constant angle of incidence $\alpha$. Following Eq. (16), upon leaving the simulation area the Stokes vector of the photon package will be

[TABLE]

Note that in this scenario the Stokes vector depends only on $\alpha$. The analytic solutions for the reduced Stokes parameters are,

[TABLE]

under the different viewing angles $\alpha$ towards the z-axis. We use this closed form validates the implementation of dichroism and birefringence in MCpol.

In Fig. 7 we compare the results of MCpol with the analytic solutions. In the top panel, the observed Stokes parameters are plotted, using lines for the analytic solution and symbols for the numerical results achieved by MCpol. The lower panel gives the relative error between the numerical and analytical solutions. As we bin the exiting photons using the cosine of their exit angles, there are more data points for large viewing angles than for small viewing angles. Generally, the analytic solution is reproduced to better than $5\%_{\rm{rel}}$. There are a few aliasing effects at $18^{\circ}$, $33^{\circ}$, $42^{\circ}$ and $61^{\circ}$ $U/I$ shows a bit larger noise. We attribute this to sampling errors caused by the finite resolution of the grid.

5.5 Albedo test case

The final test confirms the correct implementation of the albedo. The albedo of spheroids and spheroid-like particles depends on the polarization of the radiation interacting with the grains. MCpol handles this using Eq. (23). The albedo connects the effects of scattering and extinction, and a test of the implementation needs to consider both effects simultaneously.

We set up a test in which a collimated beam of right-handed circular polarized radiation propagates along an elongated dust cloud of length $r$. The dust cloud dichroically extinct and scatters the radiation. The orientation ${\boldsymbol{o}}$ of the dust particles is perpendicular to the beam ($\alpha=90^{\circ}$, Fig. 3) and no radiation scatters into the beam. The geometry of the albedo test case is illustrated in Fig. 8. We use an artificial dust grain that scatters isotropic while preserving the Stokes parameters. This is obtained by using the 4D unity matrix as the scattering matrix,

[TABLE]

The cross-sections are similar to the previous test (Eq. 51), but this time with the scattering cross sections differing from zero:

[TABLE]

The Stokes vector along $0\leq s\leq r$ in radial direction $r$ of the dust cloud before scattering is

[TABLE]

As the beam is collimated, the inverse square law does not apply here. The scattering probability is given by Eq. (6),

[TABLE]

The Stokes components of the radiation arriving at the detector are given by multiplying $\tilde{C}_{\text{sca}}$ with ${\boldsymbol{S}}(s)$,

[TABLE]

The change of the Stokes components can be illustrated as follows: The component of the radiation that is polarized parallel to the orientation ${\boldsymbol{o}}$ of the grains (out of the plane of the paper in Fig. 8) is extinct stronger than the component perpendicular to ${\boldsymbol{o}}$. Along the path through the dust cloud, the remaining radiation is more polarized in the plane of the paper, hence $Q<0$. Following Eq. 58d the circular polarization degree $V/I$ decreases along $s$ (Fig. 8). As the incidence of the radiation is perpendicular to the alignment ($\alpha=\pi/2$), the phase between the two components stays constant ($C_{\text{cpol}}=0$) and the circular polarization does not turn into linear polarization ($U\equiv 0$).

In Fig. 9 we show the results of the albedo test case solved by MCpol along with the analytic solutions. The intensity reduces by $\sim 20$ % along the profile of the dust cloud. The circular polarization degree ($V/I$) of the radiation reduces from 1 to about 0.27. These effects follow the analytic expectation with high precision, the relative error is below 0.3%. In summary, MCpol reproduces the correct behavior of the analytic solution of the albedo test case and the numerical imprecision of the Stokes vector is low (\mathbin{\lower 3.0pt\hbox{\hbox to0.0pt{\raise 5.0pt\hbox{$\char 60\relax$}\hss}\mathchar 29208\relax}}1 %).

If the albedo would not favor reflecting radiation polarized perpendicular to the alignment, the reflected intensity would behave very differently. As an example, we show the analytic solution of the intensity for test particles having $C_{\text{sca,pol}}=0$. In that case, the total intensity of the beam decreases exponentially as often applied in astronomy.

6 Summary

We have implemented polarization of radiation by spheroidal dust considering scattering, dichroic extinction, and birefringence in the Monte Carlo radiative transfer code MCpol. We provide a detailed description of the methodology used, for which the Stokes formalism was selected. We developed paradigmatic examples for testing such numerical polarization radiative transfer models using spheroidal dust against analytical solutions. They are used to verify the correct functionality of MCpol and they may be applied by other teams to verify their codes. Further, the comparison of MCpol against the analytical solutions provides means for estimating the numerical precision. We find a typical deviation of MCpol from the analytical solution for the intensity, the linear polarization degree and angle of 0.1% and in the circular polarization of 0.03%. However, situations, e.g., at a small angle of incidence $\alpha\sim 0$ or for backward scattering ($\theta\sim 180^{\rm o}$), we notice larger deviations that are due to sampling errors caused by the finite size of the applied grid of the model volume.

The analytical test cases are suited to other teams for estimating the precision of their dust polarization treatments. So far, only a limited number of MCRT codes support polarization including more complicated particle geometries than spheres. This is even though extinction due to spheres cannot explain the well-established Serkowski curve (Serkowski et al. 1975) of the diffuse ISM. We hope that the presented examples with their attached analytical solutions will help expand the number and ease the verification of MCRT codes that treat dust polarization of more complex grain shapes.

In realistic applications, one needs to consider more appropriate dust geometries, e.g., replacing spherical with spheroidal grain shapes as a starting point. For a mix of such particles that are of different kind of carbon and silicate materials and which span a range of particle sizes and grain porosity, their optical cross sections and the scattering matrix needs to be computed, furthermore some kind of grain alignment need to be implemented. MCpol is capable of tracing the polarization due to spheroidal dust grains. The here presented and tested code is open to the community in a collaborative effort.

Acknowledgements.

C.P. and M.B. acknowledge the financial support from CHARM (Contemporary physical challenges in Heliospheric and Astrophysical Models), a Phase-VII Inter-university Attraction Pole program organized by BELSPO, the Belgian federal Science Policy Office.

Appendix A Calculation of the exit angle

We begin with the vectors of the symmetry axis of the grain ${\boldsymbol{o}}$, the propagation direction before scattering ${\bf k}$ and the normal to the plane of incidence ${\boldsymbol{n}}$. By rotating ${\boldsymbol{n}}$ around ${\boldsymbol{k}}$ by $\varphi$ we obtain the normal to the scattering plane ${\boldsymbol{n}}_{\text{scat}}$. Rotating ${\boldsymbol{k}}$ around ${\boldsymbol{n}}_{\text{scat}}$ by $\theta$ results in the propagation direction after scattering ${\boldsymbol{k}}^{\prime}$. The (as of yet unknown) exit angle $\gamma$ is used to rotate ${\boldsymbol{n}}_{\text{scat}}$ around ${\boldsymbol{k}}^{\prime}$ into the normal of the plane of departure ${\boldsymbol{n}}^{\prime}$. These rotations can be calculated with Euler’s finite rotation formula (Cheng & Gupta 1989) and simplified by taking into account that many pairs of these vectors are perpendicular. With the substitutions $c_{x}=\cos x$ and $s_{x}=\sin x$ we can write

[TABLE]

This normal must be perpendicular to the symmetry axis of the grain, ${\boldsymbol{o}}\cdot{\boldsymbol{n}}^{\prime}=0$. This leads to

[TABLE]

which can be solved for $\gamma$,

[TABLE]

The arc-tangent function is undefined for $0/0$. This happens if either $s_{\varphi}=0$ and $\alpha=\theta$, or if $s_{\alpha}=0$ and $s_{\theta}=0$. The first case means that the scattering plane is the plane of incidence (because $\varphi=0$), therefore the plane of scattering is the plane of departure as well, $\gamma=0$ (or $\gamma=\pi$). In the second case, the propagation direction before scattering is (anti-)parallel to the grain orientation and the scattering plane will be the plane of departure, again $\gamma=0$.