© ESO, 2017

1. Introduction

Many astronomical objects contain or are shrouded by dust. Often, a non-negligible fraction of ultraviolet (UV) to near-infrared (NIR) radiation emitted by embedded sources is scattered or absorbed by dust grains before leaving the system. Scattered radiation is generally polarized. The polarization state of the light can be used to deduce information about the grains that would not be available using intensity measurements alone (Scicluna et al. 2015). There are indications that dust properties differ widely and systematically between galaxies (Fitzpatrick & Massa 1990; Gordon et al. 2003; Rémy-Ruyer et al. 2015; Dale et al. 2012; De Vis et al. 2016) and that they can vary significantly within a galaxy (Draine et al. 2014; Mattsson et al. 2014). Polarimetric studies can help in constraining these properties. Theoretical frameworks for modeling radiative transfer therefore usually include a section on polarization (see, e.g., Chandrasekhar 1950; Van De Hulst 1957).

Numerical simulations of dust radiative transfer most commonly use the Monte Carlo technique (see, e.g., Whitney 2011; Steinacker et al. 2013). Codes using this method track many individual photon packages as they propagate through the dusty medium, simulating emission, scattering, and absorption events based on random variables drawn from the appropriate probability distributions. While it is conceptually straightforward to track the polarization state of a photon package as part of this process, there are many details to be considered, and the implementation complexity depends on the assumptions and approximations one is willing to make. Moreover, the dust model used by the code must provide the extra properties necessary to calculate the changes to the polarization state for each interaction with a dust grain.

As a result, various authors have made different choices for implementing polarization in Monte Carlo radiative transfer (MCRT) codes. Most commonly, the MCRT codes consider only scattering by spherical dust grains (e.g., Bianchi et al. 1996; Pinte et al. 2006; Min et al. 2009; Robitaille 2011; Goosmann et al. 2014). Some codes include more advanced support for polarization by absorption and scattering off aligned spheroids (Whitney & Wolff 2002; Lucas 2003; Reissl et al. 2016) and/or for polarized dust emission (Reissl et al. 2016).

To verify the correctness of the various polarization implementations, authors sometimes compare the results between codes (e.g. Pinte et al. 2009). Because of the variations in assumptions and capabilities, however, such a comparison is tricky and the “correct” result is usually simply assumed to be the result obtained by a majority of the codes. Even when the basic assumptions about grain shape and alignment as well as the dust mixture are the same and the codes support the same polarization processes, comparing results is usually complicated.

In this paper we present a robust framework that is independent of a coordinate system for implementing polarization in a tree-dimensional (3D) MCRT code. The mathematical formulation and the numerical calculations in our method rely solely on reference frames determined by the physical processes under study (i.e., the propagation direction or the scattering plane) and not on those determined by the coordinate system (i.e., the z-axis). This approach avoids numerical instabilities for special cases (i.e., a photon package propagating in the direction of the z-axis or close to it) and enables a more streamlined implementation.

We have implemented this framework in SKIRT¹ (Baes et al. 2003, 2011; Camps & Baes 2015), a versatile multipurpose Monte Carlo dust radiative transfer code. It has been designed and optimized for systems with a complex 3D structure, as multiple components are configured separately and construct a more complex model for the dust and/or radiation sources (Baes & Camps 2015). The code is equipped with a range of efficient grid structures on which the dust can be spatially discretized, including octree, k-d tree, and Voronoi grids (Saftly et al. 2013, 2014; Camps et al. 2013). A powerful hybrid parallelization scheme has been developed that guarantees an optimal speed-up and load balancing (Verstocken et al., in prep.), and it opens up a wide range of possible polarization applications. In order to test the correct behavior and the accuracy of our implementation, we have developed a number of analytical test cases designed to validate polarization implementations in a structured manner. Furthermore, we carefully match our polarimetric conventions to the recommendations issued by the International Astronomical Union (IAU 1974).

In large-scale dust systems complex geometries arise and need to be handled by the codes. We first apply our method to some elementary models of dusty disk galaxies, enabling a qualitative comparison with the two-dimensional (2D) models of Bianchi et al. (1996). We also perform the polarization part of the Pinte et al. (2009) benchmark and compare with the published results.

We then implement spiral arms into dusty disk galaxy models and show that this produces a marked polarimetric signature tracing the positions of the arms. Our current implementation supports only scattering by spherical grains. Dichroic extinction and more complex grain shapes may also have a strong influence (Voshchinnikov 2012; Siebenmorgen et al. 2014; Draine & Fraisse 2009) and will be supported in future work.

In Sect. 2 we summarize the notation and conventions used in this paper to describe the polarization state of electromagnetic radiation, and we provide recipes for translation into other conventions. We then present our method and its implementation in Sect. 3, and the accompanying analytical test cases and their results in Sect. 4. The application of our method to benchmark tests is described in Sect. 5. The dusty spiral galaxy model is described and implemented in Sect. 6. We summarize and conclude in Sect. 7.

2. Polarization

Fig. 1 Illustration of the Stokes vector conventions recommended by IAU (1974) and used in this paper. The radiation beam travels along its propagation direction, k, out of the page. The electric field vector describes an ellipse over time. The linear polarization angle, γ, is given by the angle between the primary axis of the ellipse (green line) and the reference direction, d. The position angle of the electric field vector increases with time, the beam has right-handed circular polarization.

2.1. Stokes vector

The polarization state of electromagnetic radiation is commonly described by the Stokes vector, S (see, e.g., Van De Hulst 1957; Chandrasekhar 1950; Bohren & Huffman 1998; Mishchenko et al. 1999), (1)where I represents the intensity of the radiation, Q and U describe linear polarization, and V describes circular polarization. The degrees of total and linear polarization, P and P_L, can be written as a function of the Stokes parameters, The (linear) polarization angle, γ, can be written as (4)where arctan₂ denotes the inverse tangent function that preserves the quadrant. Combining Eqs. (3)and (4), we can also write The values of Q and U depend on the polarization angle γ, which describes the angle between the direction of linear polarization and a given reference direction, d, in the plane orthogonal to the propagation direction, k. The angle is measured counter-clockwise when looking at the source, as illustrated in Fig. 1. A linear polarization angle in the range 0 <γ<π/ 2 implies a positive U value. A radiation beam is said to have right-handed circular polarization (with V> 0) when the electric field vector position angle increases with time, and left-handed when it decreases.

The reference direction, d, can be chosen arbitrarily as long as it is well defined and perpendicular to the propagation direction. However, when the polarization state changes as a result of an interaction (e.g., a scattering event), most recipes for properly adjusting the Stokes vector require the reference direction to have a specific orientation (e.g., lying in the scattering plane). Before applying the recipe, the existing reference direction must be rotated about the propagation direction to match this requirement. This is accomplished by multiplying the Stokes vector by a rotation matrix, R(ϕ), (6)A rotation about the direction of propagation by an angle ϕ, counter-clockwise when looking toward the source of the beam, is described by the matrix (7)To record the polarization state change for a scattering event, the Stokes vector is multiplied by the Müller matrix, M, corresponding to the event, assuming that the reference direction lies in the scattering plane (as well as in the plane orthogonal to the propagation direction). The Müller matrix components depend on the geometry of the scattering event and the physical properties of the scatterer, and they often depend on the wavelength. In general, the Müller matrix is (8)where θ is the angle between the propagation directions before and after the scattering event, and λ is the wavelength of the radiation. For clarity of presentation, we drop the dependencies from the notation for the individual Müller matrix components. Including the reference direction adjustments before and after the actual scattering event, ϕ and ϕ_new, the transformation of a Stokes vector for a scattering event can thus be written as(9)For scattering by spherical particles, the Müller matrix simplifies to (Krügel 2002) (10)again assuming that the reference direction lies in the scattering plane. The Müller matrices for a particular grain size and material can be calculated using Mie theory (see e.g., Voshchinnikov & Farafonov 1993; Bohren & Huffman 1998; Peña & Pal 2009).

For scattering by electrons, also called Thomson scattering, the Müller matrix is wavelength-independent and can be expressed analytically as a function of the scattering angle (Bohren & Huffman 1998), (11)

2.2. Conventions

In this paper we define the Stokes vector following the recommendations of the International Astronomical Union (IAU 1974), as presented in Sect. 2.1 and illustrated in Fig. 1. Historically, however, authors have used various conventions for the signs of the Stokes parameters U and V (Hamaker & Bregman 1996, see also a recent IAU announcement²). For example, the polarization angle γ is sometimes measured while looking toward the observer rather than toward the source, flipping the sign of both U and V. Reversing the definition of circular polarization handedness also flips the sign of V. Table 1 lists a number of references with the conventions adopted by the authors.

Assuming that the adopted conventions are properly documented, translating the values of the Stokes parameters from one convention into another is straightforward – by flipping the signs appropriately. When comparing or mixing formulas and recipes constructed for different conventions, changes in the signs of U and V affect the sign of the Müller matrix components (Eq. (8)) in the third row and column and fourth row and column, respectively. Mathematically this can be described by multiplying the Müller matrix by signature matrices, with σ and ς being + 1 or −1. In case of the rotation matrix, Eq. (7), this implies that the signs in front of the sine functions change based on the chosen convention.

3. Method

3.1. SKIRT code

SKIRT (Baes et al. 2003, 2011; Camps & Baes 2015) is a public multipurpose MCRT code for simulating the effect of dust on radiation in astrophysical systems. It offers full treatment of absorption and multiple anisotropic scattering by the dust, self-consistently computes the temperature distribution of the dust and the thermal dust reemission, and supports stochastic heating of dust grains (Camps et al. 2015). The code handles multiple dust mixtures and arbitrary 3D geometries for radiation sources and dust populations, including grid- or particle-based representations generated by hydrodynamical simulations (Camps et al. 2016).

SKIRT is predominantly used to study dusty galaxies (Baes et al. 2010; De Looze et al. 2012, 2014; De Geyter et al. 2014, 2015; Saftly et al. 2015; Mosenkov et al. 2016; Viaene et al. 2017), but it has also been applied to active galactic nuclei Stalevski et al. 2012, 2016), molecular clouds (Hendrix et al. 2015), and binary systems (Deschamps et al. 2015; Hendrix et al. 2016).

Employing the MCRT technique, SKIRT represents electromagnetic radiation as a sequence of discrete photon packages. Each photon package carries a specific amount of energy (luminosity) at a given wavelength. A SKIRT simulation follows the individual paths of many photon packages as they propagate through the dusty medium. The photon package life cycle is governed by various events determined stochastically by drawing random numbers from the appropriate probability distributions. A photon package is created at a random position based on the luminosities of the sources and is emitted in a random direction depending on the (an)isotropy of the selected source. Depending on the dust material properties and spatial distribution, the photon package then undergoes a number of scattering events at random locations (using forced scattering; see Cashwell & Everett 1959), and is attenuated by absorption along its path (using continuous absorption; see Lucy 1999; Niccolini et al. 2003).

To boost the efficiency of the simulation and reduce the noise levels at the simulated observers, SKIRT employs the common peel-off optimization technique (Yusef-Zadeh et al. 1984). Rather than waiting until a photon package happens to leave the system under study in the direction of one of the observers, a special photon package is peeled off in the direction of each observer at each emission and scattering event, including an appropriate luminosity bias for the probability that a photon package would indeed be emitted or scattered in that direction. Meanwhile, the original or main photon package continues its trajectory through the dust until its luminosity has become negligible and the package is discarded.

For the purposes of this paper, we assume that a newly emitted photon package represents unpolarized radiation, and its polarization state is not affected by the attenuation along its path through the dusty medium. We assume that the photon package is scattered by spherical dust grains (which does affect the polarization state). This leaves us with three types of events: scattering the main photon package into a new direction, peeling off a special photon package toward a given observer, and detecting a peel-off photon package at an observer. As a first step toward describing the procedures for each of these events, we discuss our approach for handling the Stokes vector reference direction.

3.2. Reference direction

Fig. 2 Geometry of a scattering event. The angle θ is between the incoming and outgoing propagation directions k and knew. The angle ϕ is given by the normals of the previous and the present scattering planes n and nnew.

As noted in Sect. 2, the Stokes vector describing the polarization state of a photon package is defined relative to a given reference direction, d, in the plane perpendicular to the propagation direction. We define a new direction, n, perpendicular to both the propagation direction, k, and the reference direction, d, which are perpendicular to each other as well, so that(14)assuming all three vectors are unit vectors. By definition, the scattering plane contains both the incoming and outgoing propagation directions k and k_new. Consider the situation before the event (also, see Fig. 2). If the reference direction d, which is always perpendicular to k, lies in the scattering plane as well, then n corresponds to the normal of the scattering plane. A similar situation applies after the scattering event. We store n rather than d with each photon package, and our procedures below are described in terms of n.

Some authors (e.g., Chandrasekhar 1950; Code & Whitney 1995; Gordon et al. 2001) choose to rotate the Stokes reference direction in between scattering events into the meridional plane of the coordinate system. Their procedure uses two rotation operations for each scattering event, one before and one after the event, and requires special care to avoid numerical instabilities when the propagation direction is close to the z-axis. The latter occurs because the meridional plane is then ill defined.

We leave the reference direction unchanged after a scattering event and instead perform a single rotation as part of the next scattering event. This method is also applied by Fischer et al. (1994) and Goosmann & Gaskell (2007) and illustrated in Fig. 2. The current scattering plane (red) includes the incoming and outgoing propagation directions k and k_new, defining the scattering angle θ. After the previous scattering event, the reference direction has been left in the previous scattering plane (blue), so the angle between the normals n and n_new to the previous and the current scattering planes determines the angle ϕ over which the Stokes vector must be rotated to end up in the current scattering plane. The transformation of the Stokes vector given in Eq. (9)can therefore be simplified to (15)Care must be taken to properly set a reference normal n for newly emitted photon packages that have not yet experienced a scattering event. Because we assume our sources to emit unpolarized radiation, we can pick any direction perpendicular to the propagation direction. We postpone the details and justification for this procedure to Sect. 3.7.

3.3. Scattering phase function

The probability that a photon package leaves a scattering event along a particular direction k_new for an incoming direction k is given by the phase function, Φ(k,k_new) ≡ Φ(θ,ϕ), where θ and ϕ represent the inclination and azimuthal angles of k_new relative to k, and where we omit the wavelength dependency from the notation. In the formulation of Sect. 2, we can say that the phase function is proportional to the ratio of the beam intensities, I and I_new, before and after the scattering event, (16)For spherical grains, combining Eqs. (7), (10), and (15)leads to (17)and therefore, (18)Using Eq. (4) and introducing a proportionality factor, N, we can write (19)The proportionality factor is determined by normalizing the phase function (a probability distribution) to unity. Integration over the unit sphere yields

3.4. Sampling the phase function

After scattering, a new direction of the photon package is determined by sampling random values for θ and ϕ from the phase function Φ(θ,ϕ). To accomplish this, we use the conditional probability technique. We reduce the phase function to the marginal distribution Φ(θ), (22)We sample a random θ value from this distribution through numerical inversion, that is to say, by solving the equation (23)for θ, where is a uniform deviate, that is, a random number between 0 and 1 with uniform distribution.

Once we have selected a random scattering angle θ, we sample a random azimuthal angle ϕ from the normalized conditional distribution, where the ratio S₁₂/S₁₁ depends on θ. This can again be done through numerical inversion, by solving the equation for ϕ, with being a new uniform deviate.

3.5. Updating the photon package

After randomly selecting both angles θ and ϕ, we can use Eq. (15)to adjust the main photon package’s Stokes vector. We can also calculate the outgoing propagation direction k_new and the new reference normal n_new from the incoming propagation direction k and the previous reference normal n (see Fig. 2). We use Euler’s finite rotation formula (Cheng & Gupta 1989) to rotate a vector v about an arbitrary axis a over a given angle β (clockwise while looking along a), (28)The last term of the right-hand side vanishes when the vector v is perpendicular to the rotation axis a.

In our configuration, the reference normal n rotates about the incoming propagation direction k over the azimuthal angle ϕ. Because n is perpendicular to k, we have (29)Furthermore, the propagation direction rotates in the current scattering plane, that is, about n_new, over the scattering angle θ. Again, k is perpendicular to n_new, so that (30)

3.6. Peel-off photon package

As described in Sect. 3.1, a common MCRT optimization is to send a peel-off photon package toward every observer from each scattering site. The peel-off photon package must carry the polarization state after the peel-off scattering event, and its luminosity must be weighted by the probability that a scattering event would indeed send the outgoing photon package toward the observer. To obtain this information, we need to calculate the scattering angles θ_obs and ϕ_obs, given the outgoing direction of the peel-off scattering event, or in other words, the direction toward the observer, k_obs. This is effectively the scattering problem in reverse, in which random angles were chosen based on their probability, and the new propagation direction was calculated from these angles.

Finally, when the peel-off photon package reaches the observer, its Stokes reference direction must be rotated so that it lines up with the direction of north in the observer frame, k_N, according to the IAU (1974) conventions. The scattering angle θ_obs is easily found through the scalar product of the incoming and outgoing directions, (31)Because 0 ≤ θ_obs ≤ π, the cosine unambiguously determines the angle.

To derive the azimuthal angle ϕ_obs, we recall (Fig. 2) that it is the angle between the normal to the previous scattering plane, n, and the normal to the peel-off scattering plane, n_obs. The latter can be obtained through the cross product of the incoming and outgoing directions, (32)We need both cosine and sine to unambiguously determine ϕ_obs in its 2π range. We easily have (33)Because k is perpendicular to both n and n_obs, the following relation holds, (34)or, after projecting both sides of the equation on k, (35)This derivation of ϕ_obs breaks down for a photon package that happens to be lined up with the direction toward the observer before the peel-off event. Indeed, in this special case of perfect forward or backward peel-off scattering, Eq. (32)is undefined. However, because the scattering plane does not change, it is justified to simply set ϕ_obs = 0 instead.

We insert the calculated θ_obs and ϕ_obs values into Eq. (15)to adjust the peel-off photon package’s Stokes vector, and we also update the reference normal to n_obs. When the photon package’s polarization state is recorded at the observer, its Stokes vector reference direction must be parallel to the north direction, k_N, in the observer frame. This is equivalent to orienting the reference normal along the east direction, k_E = k_obs × k_N. The angle, α_obs, between n_obs and k_E can be determined using a similar reasoning as for ϕ_obs, so that (36)and (37)The final adjustment to the Stokes vector is thus a rotation (see Eqs. (6) and (7)) with the matrix R(α_obs). The Stokes vector is indifferent to rotations by π. Using k_W = −k_E yields the same result.

3.7. Reference direction for new photon packages

We now return to the issue of selecting a reference direction, or more precisely, a reference normal, for newly emitted photon packages. We stated at the end of Sect. 3.2 that we can pick any direction perpendicular to the propagation direction, because we assume our sources to emit unpolarized radiation. Indeed, it is easily seen from Eq. (25)that the probability distribution for the azimuthal angle ϕ becomes uniform for unpolarized incoming radiation, that is, with P_L,in = 0. Consequently, our choice of reference normal in the plane perpendicular to the propagation direction will be completely randomized after the application of the scattering transformation (Eq. (15)).

We determine the reference normal, n = (n_x,n_y,n_z), perpendicular to the propagation direction, k = (k_x,k_y,k_z), using When k is very closely aligned with the Z-axis, the root in these equations vanishes, and we instead select n = (1,0,0) as the reference direction.

4. Validation

Fig. 3 Top: geometry used in the analytical test cases. The z-axis is toward the reader. For test case 1, a central point source illuminates thin slabs of electrons that are slanted by 45° toward the observer. The central source is replaced by a small blob of electrons for the other test cases. The blob is illuminated by a collimated beam from within the xy plane for test case 2 and from below the xy plane for test cases 3 and 4. Bottom: intensity map of test case 1 as seen from the observer.

4.1. Test setup

In order to confirm the validity of our method and its implementation in SKIRT, we develop four test cases for which the results can be calculated analytically. The analytical results are obtained solely using the formalisms of Sect. 2, so that taken together, the test cases verify most aspects of the procedures presented in Sect. 3.

The overall setup for the test cases is illustrated in Fig. 3. A central source illuminates two physically and optically thin slabs of material, which scatter part of the radiation to a distant observer. The slabs are arranged on the sides of a square rotated by 45° relative to the line of sight and are spatially resolved by the observer’s instrument. To simplify the calculations we only consider radiation detected close to the xy plane, essentially reducing the geometry to two dimensions. Because of the low scattering probability in the slabs, the path with the least number of scattering events will greatly dominate the polarization state at each instrument position. This allows us to deduce the dominating scattering angle, θ, corresponding to each instrument position along the x-axis. Referring to Fig. 3, geometrical reasoning leads to (41)with the plus sign for x> 0 and the minus sign for x< 0. In combination with analytical scattering properties for the slab material, it then becomes possible to derive a closed-form expression for the components of the Stokes vector at each position.

Because the observer is considered to be at “infinite” distance, we can use parallel projection, and the distance from the slabs to the observer does not vary with x. However, we do need to take into account the variations in the path length ℓ from the central source to the slabs because it affects the amount of radiation arriving at the slabs as a function of x. Geometrical reasoning in Fig. 3 again leads to (42)The scattering properties of the slabs and the makeup of the central source vary between the test cases. For test cases 1 through 3, the slabs contain electrons, with scattering matrix M_Th given by Eq. (11). We study the observed intensity, I, the degree of linear polarization, P_L, and the linear polarization angle, γ, of these test cases in Sect. 4.2. Scattering by electrons never causes circular polarization. Therefore in Sect. 4.3 we introduce slabs that contain synthetic particles with custom-designed scattering properties (test case 4). This allows us to study the observed circular polarization.

Test case 1 has a central point source. For the remaining test cases, the central source is replaced by a small blob of electrons illuminated by a collimated beam positioned at varying angles, so that the center becomes the site of first scattering.

A numerical implementation of the test cases will always discretize certain aspects of the theoretical test setup. For our implementation in SKIRT, we made the following choices. The spatial domain of the setup is divided into 601 × 601 × 61 cuboidal grid cells lined up with the coordinate axes (we use odd numbers to ensure that the origin lies in the center of a cell). The cells overlapping the slabs (and where applicable, the central blob) contain electrons, the other cells represent empty space. The edges of the cells and slabs are not aligned. Each detector pixel provides averages over multiple cells, and the slabs are optically thin (τ = 10^-3 along their depth). The length of the slabs is 1.141, their depth is 0.006, and their height is 2 × 10^-5. The detector in the observer plane has a resolution of 201 pixels along the x-axis. We use 10¹⁰ photon packages for each test case to minimize the stochastic noise characteristic of Monte Carlo codes. SKIRT uses precomputed tables of various quantities with a resolution of 1° in θ and ϕ, to help perform the numerical inversions in Eqs. (23)and (27).

4.2. Linear polarization

Fig. 4 Intensity (top row), linear polarization degree (middle row), and polarization angle (bottom row) of the observed radiation in test cases 1 through 3 (left to right columns). The top section in each panel shows the analytical solution (black) and the model results (dashed orange). The bottom section in each panel shows the absolute differences (blue) and relative differences (shaded area) of the analytic solution and the model. The green lines are calculated using twice the resolution of the θ scattering angle.

Test case 1 is designed to test the peel-off procedure described in Sect. 3.6. The slabs contain electrons, and the central point source emits unpolarized photon packages. When a photon package’s direction is toward one of the slabs, the forced interaction algorithm of SKIRT initiates a scattering event at the slab and a peel-off photon package is sent toward the observer. Because the slabs are optically thin, the luminosity of the scattered original photon package is small. It is subsequently deleted and a new photon package is launched.

As a result, the Stokes vector of the observed photon packages can be written as (43)reflecting from right to left a scattering transformation starting from an unpolarized state (and thus ϕ = 0), a rotation to align the reference direction with the observer frame, and the dependency on the path length from the source to the slab. We can now substitute Eqs. (7), (11), (41), and (42)into this equation. Because the arctangent of Eq. (41)is used as an argument for the cosine in Eq. (11), the final equations for the intensity and the linear polarization degree reduce to polynomials, The orientation of the polarization is perpendicular to the scattering plane, that is, north/south or along the y-axis in the observer frame.

In test case 2 we add a scattering event to verify part of the procedures for scattering the main photon package. To this end, we replace the central point source with a small blob of electrons at the same location, and illuminate this electron blob with a collimated beam positioned at (− 1,1,0) and oriented parallel to the slabs toward the bottom right (see Fig. 3). In this setup, a photon package emitted by the collimated source can reach the slabs only after being scattered by the electrons in the central blob. This effectively adds a forced first scattering event to all photon packages reaching the observer, with a scattering angle that can be deduced from the geometry. The peel-off scattering angle is still given by Eq. (41). The scattering angle for the first scattering event in the central electron blob is (46)again with the plus sign for x> 0 and the minus sign for x< 0. Because all components of the setup are in the same plane, the scattering plane is always the same (the xy-plane). The Stokes vector of the observed photon packages can thus be written as(47)The intensity and linear polarization degree for test case 2 are The orientation of the polarization is north/south.

In the third test case we move the collimated source below the xy-plane to , so that we can test the rotation of the Stokes vector reference direction between scattering events. The source is now placed at an inclination of 165° (relative to the z-axis) and an azimuthal angle of 30° (clockwise from the x-axis). It still points toward the central electron blob, that is, toward the top right and out of the page in Fig. 3. As a result, the normal of the first scattering plane is tilted, while the normal of the peel-off scattering plane remains aligned with the z-axis. With some trigonometry, we arrive at expressions for the angles involved in the first (main) and second (peel-off) scattering events, with the plus sign for x> 0 and the minus sign for x< 0. The expression provided in Eq. (51)for ϕ₁ is simplified and shifted by ±π for some x. The Stokes vector is invariant under rotations by π.

The Stokes vector of the observed photon packages can now be written as (54)To limit the complexity of presentation, we provide expressions for the Stokes parameters from which the linear polarization degree and angle can be calculated using Eqs. (3)and (4), Figure 4 compares the analytical solutions and SKIRT results for the observed intensity, linear polarization degree, and polarization angle for these three test cases. The intensity curves show a relative noise level of on average about 3%. The linear polarization degrees are identical to below 0.1% absolute for test cases 1 and 2 and below 1% absolute for test case 3. The polarization angles from north are correct to below 0.05° for test cases 1 and 2 and below 1° for test case 3. While the linear polarization degree and position angle can be determined from a simulation with relatively few photon packages, reducing the noise in the intensity requires significantly more photon packages. This is because the number of photon packages arriving at each pixel is subject to Poisson noise, whereas the path that each photon package takes to the same pixel is defined within tight boundaries. The intensity curve depends on the number of photons. The linear polarization degree and angle are independent of the number of photon packages. Their noise is due to the size of the pixels. Slightly different paths and scattering angles might contribute to the same pixel.

The polarization angle for test case 1 shows an intriguing spike near x = 0. As we can see in the linear polarization degree curve and from Eq. (45), the radiation arriving at x = 0 is unpolarized. This implies that both the Q and U components of the Stokes vector are zero and the polarization angle becomes undefined (see Eq. (4)). This in turn causes small numerical inaccuracies in the calculations for photon packages arriving very close to x = 0.

The relative differences for the other quantities are generally smaller than a fraction of a percent. In the results of test case 3 the residuals of the linear polarization degree and polarization angle contain spikes that are resolved by multiple pixels each and symmetric with respect to x = 0. The residual of the intensity curve shows a similar effect. It tends to be lower than expected in the outer regions and higher than expected in the inner region. These orderly deviations indicate a systematic difference rather than noise. In fact, these discrepancies are caused by resolution effects in the SKIRT implementation of our method (see Sect. 4.1 for a description of our discretization choices). For example, consider the interval 0.6 <x< 0.7 for test case 3, which is resolved by 10 pixels along the x-axis of the detector in the observer frame. The corresponding interval for scattering angle θ₁ (see Eq. (50)) is 75.0°<θ₁< 75.1°, that is, only a fraction of the 1° resolution in the SKIRT calculations related to θ. It is obvious that this lack of angular resolution relative to the output resolution will cause inaccuracies. We calculated the residuals in the polarization degree and angle from north using twice the θ resolution and show them in pale green in Fig. 4. They confirm that increasing the θ resolution in the calculations indeed reduces the discrepancies accordingly.

4.3. Circular polarization

To include circular polarization in test case 4, we use synthetic particles similar to electrons, but with a scattering matrix that mixes the U and V components of the Stokes vector: (59)We use the same geometry as for test case 3, replacing the electrons in the two slabs and in the central blob by particles described by Eq. (59). The angles are still described by Eqs. (50)through (53), Eq. (54)still holds, and the Stokes parameters of the observed photon packages become In Eq. (62) the plus sign is again for x > 0 and the minus sign for x < 0.

Fig. 5 Fraction of circular polarization, V/I, of the observed radiation in test case 4. The top panel shows the observed surface brightness overlaid with arrows indicating the handedness. The arrow length scales linearly with the magnitude. The middle and bottom panels compare the analytical solution with the model results, as in Fig. 4.

Figure 5 shows the (relative) circular polarization, V/I, for this test case, again comparing the analytical solutions with the SKIRT results. The relative differences between the analytical and simulated results for | x | < 0.7 are below one percent. The larger discrepancies for | x | > 0.7 are again due to the limited resolution in the SKIRT calculations related to θ. The pale green line again shows the residuals when calculating with 0.5° resolution and has a significantly smaller residual curve. Disregarding the outer part, the circular polarization is correct to 0.1% absolute.

5. Benchmark tests

Fig. 6 B-band surface brightness maps (color scale) overlaid with linear polarization maps (line segments) for inclinations of 20° (left), 75° (middle), and 90° (right column). The top row shows a model with only a stellar bulge, the middle row a model with only a stellar disk, and the bottom row a model with both stellar bulge and disk with a ratio of bulge to total luminosity of B/T = 0.5. The dust disk is the same in the three cases and has a central face-on V-band optical depth of 10.

Fig. 7 Polarization degree maps for two inclinations of a disk with high optical density (τ = 106). On the left are the linear polarization degree maps we calculated using SKIRT. The dust grain size is the same as the wavelength (1 μm), creating the intricate pattern of the polarization degree. On the right are cuts 1–6 through the maps along with results of various codes. Pinball (Watson & Henney 2001) in green, MCMax (Min et al. 2009) in black, MCFOST (Pinte et al. 2006) in blue, TORUS (Harries 2000) in dashed orange, and SKIRT in red. The data of the other codes were taken from the official website.

5.1. Disk galaxy

We compare results of our code to Bianchi et al. (1996) as a first test of our implication of dust scattering polarization (rather than just Thompson scattering). Bianchi and collaborators describe the polarization effects of scattering by spherical dust grains in monochromatic MCRT simulations of 2D galaxy models at the B band (0.44 μm) and I band (0.9 μm). The models include a stellar bulge, a stellar disk, and a dust disk. The stellar bulge is described by Jaffe (1983; scale radius 1.86 kpc, truncated at 14.8 kpc), and the stellar disk is a double-exponential disk (scale length 4 kpc, truncated at 24 kpc; scale height 0.35 kpc, truncated at 2.1 kpc). The relative strength of the stellar components and the stellar sources is varied. In all models, the dust is assumed to be homogeneously distributed in a double-exponential disk embedded in the stellar disk, with the same horizontal parameters as the stellar disk, but much thinner vertically (scale height 0.14 kpc, truncated at 0.84 kpc). The central face-on optical depth of the dust disk τ_V is a free parameter. The properties of the dust are taken from the MRN dust model (Mathis et al. 1977), which includes graphite grains and astronomical silicates with a size distribution n(a) ∝ a^-3.5, 0.005 μm < a < 0.25 μm. Polycyclic aromatic hydrocarbons (PAHs) and very small grains are not included.

Figure 6 displays our results, which correspond to the Bianchi et al. (1996) model shown in their Fig. 8. We use the same geometry, described above, and the same bulge-to-total ratio (B/T = 0.5), wavelength (λ = 0.44 μm), and optical depth (τ_V = 10). We use a different dust model. It includes a wider range of graphite and silicate grain sizes following the Zubko et al. (2004) size distribution, in addition to a small fraction of PAHs as described in Camps et al. (2015). The differences between the dust models do not affect the qualitative results at optical wavelengths, but do prevent a quantitative comparison.

Surface brightness maps (color scale) overlaid with a linear polarization maps (line segments) are shown in Fig. 6. The inclinations range from nearly face-on (20°, left) to edge-on (90°, right). The top and middle rows show models from which the stellar disk and the stellar bulge, respectively, was removed. The bottom row shows the B/T = 0.5 model, that is to say, the model including both stellar bulge and stellar disk. The dust disk is identical for the three models.

Overall, our results are compatible with those reported by Bianchi et al. (1996). The largest polarization degrees are observed near the major axis. For pure disks, the face-on view shows very little polarization. The polarization degree increases with inclination to a maximum near an inclination of about 80° and slightly decreases again when approaching 90°. The polarization degree averaged over all pixels is below 1% for all inclinations. For pure bulges, the maximum polarization degree is largely independent of the inclination. In the outer regions of the dust disk, the linear polarization degree is 22%. In the mixed B/T = 0.5 model, the polarization degrees are drastically reduced. We find values up to 1.6%, comparable to the results of Bianchi et al. (1996), who find 1 to 1.5%.

We also test the corresponding model for the I band (λ = 0.9 μm) with a color-adjusted bulge luminosity (B−I = 1 mag), again confirming overall agreement Bianchi et al. (1996).

Fig. 8 Spiral galaxy model described in Sect. 6 and Table 2 observed at a wavelength of 1 μm. Rows from top to bottom: face-on, inclined (53°), and edge-on surface brightness (color scale) overlaid with linear polarization degree and orientation (line segments); linear polarization degree of the edge-on view, averaged over the vertical axis. Columns from left to right: the reference model; the same model observed from a different azimuth angle (rotated clockwise by 120°); a model without the spiral arm perturbations in the stellar disks, with the rotated orientation.

5.2. Dusty disk around star

To test the performance of our code on a problem of a different scale, we compute polarization results of Pinte et al. (2009). In it, a thick dusty disk surrounds a central star. The star extends out to 2 AU and has a temperature of 4000 K. The dust consists of spherical grains with a radius of 1 μm, and the light has a wavelength of 1 μm as well. The dust density distribution ρ is cylindrical, (64)with a surface density Σ₀, scale radius R₀ = 100 AU, radial distance from the center R, vertical distance from the midplane z, and scale height h₀ = 10 AU. The disk is truncated at R_min = 0.1 AU and R_max = 400 AU. The surface density depends on the total dust mass m = 3 × 10^-5M_⊙ by (65)and the albedo of the dust is 0.6475 and the opacity 4752 cm²/g at 1 μm. We adopt the scattering matrix as provided by Pinte et al. (2009)³. The system is resolved at a distance of 140 pc by a detector with 251 × 251 pixels covering 900 × 900AU².

Figure 7 shows the linear polarization maps calculated by our code for inclinations of 69.5°and 87.1°. The flux difference from the borders to the center area is about 17 orders of magnitude. The maps are displayed in gray outside the truncation radius, where no flux was recorded. The intricate pattern of the polarization degree is a result of the uniform grain radius being equal to the wavelength. The phase function therefore contains resonances and steep gradients for small changes of the scattering angle. We compare our results to the results of the four polarization capable codes in Pinte et al. (2009) along six cuts through the maps. The first and third row show the polarization degree along the cuts, and the second and fourth row show the difference of the codes to the average of all results. The TORUS code (Harries 2000) is not included in calculating the averages because its signal-to-noise ratio is too low. In the central area the codes agree within 10% of absolute polarization, but as the true result is unknown, a quantitative analysis of the results of this benchmark is difficult. In general, the results of the SKIRT code are close to the average result, and SKIRT seems to agree particularly well with the Pinball code (Watson & Henney 2001). Pinball employs some of the same optimization techniques that SKIRT uses (e.g., forced interaction and peel-off, named “forced escape” in their paper).

6. Application: spiral galaxy models

Fig. 9 Linear polarization degree, averaged over the vertical axis, of the rotated model (column B of Fig. 8) observed at 1 μm, for inclinations ranging from edge-on (i = 90°) to face-on (i = 0°).

Fig. 10 Averaged linear polarization degree for one side of the edge-on view of the rotated model, observed at optical and near-infrared wavelengths from 0.28 μm to 2.84 μm. The red vertical bands correspond to the bands in Fig. 9 and trace the tangent points of the dust spiral arms.

We study the polarization properties of a 3D galaxy model including spiral arms to investigate how the spiral arm structure is imprinted in the polarization structure. We also study this in the edge-on view when the structure is not easily characterized from the intensity alone.

We still assume a homogeneous distribution for the stellar sources and the dust. The model includes a stellar bulge and disk with an older star population, a second stellar disk with a younger star population, and a dust disk. The relevant parameters are listed in Table 2. The bulge is modeled by a spheroidal density distribution obtained by flattening a spherical Sérsic profile (Sérsic 1963), as implemented by Baes & Camps (2015). The distributions of the two stellar and the dust disks are truncated double-exponential disks. The spiral arm structure is introduced by adding a perturbation to the overall density profile, as presented by Baes & Camps (2015). We have made the spiral arms in the older stellar population “lead” those in the younger stellar population and in the dust disk by varying the spiral arm phase zero-points. The emission of the stellar populations is modeled as blackbody spectra at the indicated temperatures. We use the dust model of Sect. 5.1 (Camps et al. 2015). The total dust mass is given in Table 2, the central face-on B-band optical depth of the dust disk is τ_V ≈ 1.3.

Figure 8 shows surface brightness maps (color scale) overlaid with linear polarization maps (line segments). The leftmost column is for this model at 1 μm. The top row shows the model face-on, the second row inclined (53°), and the third row edge-on. The polarization degree is up to 1% around the central part of the models and over the whole map the average polarization degree is similar to the average polarization degree from the B/T = 0.5 model from before. As in the B/T = 0.5 model, the orientation of the polarization is circular around the central bulge, and for the inclined view the polarization degree left and right of the center increases, while it decreases behind and in front of it. In contrast to the azimuthally uniform model, there is a spiral structure in the polarization map. The linear polarization degree is higher in the arm regions and disappears in the interarm region. The maximum polarization degree is slightly inward from the regions of the arms with the highest flux. The panel in the bottom row plots the linear polarization degree for the edge-on view, averaged over the vertical axis. This average is obtained by summing each individual component of the Stokes vector and calculating the polarization degree from these totals. Regions with higher linear polarization (up to 2%) clearly trace the spiral arms and are prominent at all inclinations, including the edge-on view. The maxima in the polarization signature of the edge-on view match the positions of the spiral arms along the line of sight.

To verify this, the middle column of Fig. 8 shows the same model from a different azimuth angle. The peaks in the polarization signature align with the tangent points of the spiral arms, which are now farther out from the center of the galaxy.

In the rightmost column of Fig. 8 we remove the spiral arms perturbations from the stellar disks in the model. The polarization signature remains essentially unchanged; the maxima are slightly higher (by a factor of up to 1.2), but the structure is the same. The signature could also be produced by the different phase zero-points of the old stars and the dust (see Table 2). We calculated results using the same phase zero-point for all components (not shown here). The outer maxima are lower (by a factor of 0.8), while the inner maxima are unchanged. This confirms that the polarization signature is created by the distribution of the dust and not by the distribution of the sources.

In Fig. 9 we further study the effect of inclination on the observed polarization signature for the reference model. In the central region (r ≲ 2 kpc), the bulge emission masks most polarization at inclinations above 40°. Outside this region, however, the form of the curves is very similar for all inclinations. Although the polarization degree generally decreases toward lower inclinations, the peaks remain at the spiral arm tangent points, and the ratio between the maximum and minimum polarization degree remains roughly stable at a factor of about 2.

In Fig. 10 we compare the edge-on polarization signature of our reference model for various optical and NIR wavelengths. The polarization peaks remain prominent over the wavelength range 0.5 μm ≲ λ ≲ 2 μm. At shorter wavelengths the general polarization degree is higher and the signature is reversed, light from within the arms is slightly less polarized than light from the inter-arm regions. The increased interaction cross section causes the inter-arm dust to become efficient at scattering the stellar radiation, boosting the polarization degree. We find that in the spiral arms the ratio of once scattered to multiple times scattered light is 2.5:1, while in the inter-arm region it is 3.3:1. The polarization orientation after multiple scatterings is less uniform, which lowers the polarization degree in the arm regions.

At longer wavelengths, the signature retains the same form, but the reduced scattering efficiency of the dust causes the polarization degree to be very low, so that the peaks become hard to discern.

Our results imply that polarization measurements could be used, at least in principle, to study the spiral structure of edge-on spiral galaxies, where intensity measurements alone have limited diagnostic power. The contrast would be highest at around 1 μm and with an expected polarization degree of up to 0.6%, this is well within the capabilities of current polarization capable telescopes.

We note that the polarization degree of the edge-on galaxy NGC 891 was mapped by Montgomery & Clemens (2014) at 1.6 μm. They found polarization degrees of below 1% that varied along the disk profile. We expect the orientation of polarization due to scattering to be perpendicular to the disk. Montgomery & Clemens (2014) found the orientation to be rather parallel to the disk and attributed most of it to dichroic extinction. Our code does not yet support dichroic extinction, so we cannot compare the strength of these two effects.

7. Conclusion

We presented a robust framework that is independent of a coordinate system for implementing polarization in a 3D MCRT code, focusing on scattering by spherical dust grains. The mathematical formulation and the numerical calculations in our method rely solely on the scattering planes determined by the physical processes rather than by the coordinate system. This approach avoids numerical instabilities for special cases and enables a more streamlined implementation. We described four test cases with well-defined geometries and material properties, yielding analytical solutions. These setups are designed to validate the calculated results in a structured manner, and they can serve as benchmarks for other implementations as well.

We reconstructed a selection of the 2D models of Bianchi et al. (1996), confirming that our implementation reproduces their results at least qualitatively. A quantitative comparison is not possible because of differences in the dust model. We then calculated results for the polarization part of the Pinte et al. (2009) benchmark and obtained similar results as the other codes that took part in it. As an application of our code we constructed a 3D spiral galaxy model including a stellar bulge and disk with an older star population, a second stellar disk with a younger star population, and a dust disk. The stellar and dust distributions feature an analytical spiral arm perturbation. We showed that scattering of light at the dust in the spiral arms produces a marked polarimetric signature. It traces the tangent positions of the arms for wavelengths in the range 0.5 μm ≲ λ ≲ 2 μm, regardless of inclination.

It is fair to note, however, that our current implementation is limited to scattering by spherical dust grains. We plan to add support for polarized emission and for scattering and absorption by (partially) aligned spheroidal grains in future work. With the relevant physics included, we can also study the influence of changes to dust properties on the polarization signature, which we have not addressed in this paper.

Acknowledgments

C.P., P.C., and M.B. acknowledge the financial support from CHARM (Contemporary physical challenges in Heliospheric and AstRophysical Models), a Phase-VII Interuniversity Attraction Pole program organized by BELSPO, the BELgian federal Science Policy Office. M.S. acknowledges support by FONDECYT through grant No. 3140518 and by the Ministry of Education, Science and Technological Development of the Republic of Serbia through the projects Astrophysical Spectroscopy of Extragalactic Objects (176001) and Gravitation and the Large Scale Structure of the Universe (176003). We are grateful for the generous support from the COST Action MP1104, “Polarization as a tool to study the Solar System and beyond”. We thank Rene Goosmann, Francesco Tamborra, and Frederic Marin for many useful discussions and providing insights on the operation of the STOKES code. We wish to thank the anonymous referee for their constructive input to this paper.

References

All Tables

All Figures

Fig. 1 Illustration of the Stokes vector conventions recommended by IAU (1974) and used in this paper. The radiation beam travels along its propagation direction, k, out of the page. The electric field vector describes an ellipse over time. The linear polarization angle, γ, is given by the angle between the primary axis of the ellipse (green line) and the reference direction, d. The position angle of the electric field vector increases with time, the beam has right-handed circular polarization.

In the text

	Fig. 2 Geometry of a scattering event. The angle θ is between the incoming and outgoing propagation directions k and knew. The angle ϕ is given by the normals of the previous and the present scattering planes n and nnew.
In the text

Fig. 3 Top: geometry used in the analytical test cases. The z-axis is toward the reader. For test case 1, a central point source illuminates thin slabs of electrons that are slanted by 45° toward the observer. The central source is replaced by a small blob of electrons for the other test cases. The blob is illuminated by a collimated beam from within the xy plane for test case 2 and from below the xy plane for test cases 3 and 4. Bottom: intensity map of test case 1 as seen from the observer.

In the text

Fig. 4 Intensity (top row), linear polarization degree (middle row), and polarization angle (bottom row) of the observed radiation in test cases 1 through 3 (left to right columns). The top section in each panel shows the analytical solution (black) and the model results (dashed orange). The bottom section in each panel shows the absolute differences (blue) and relative differences (shaded area) of the analytic solution and the model. The green lines are calculated using twice the resolution of the θ scattering angle.

In the text

	Fig. 5 Fraction of circular polarization, V/I, of the observed radiation in test case 4. The top panel shows the observed surface brightness overlaid with arrows indicating the handedness. The arrow length scales linearly with the magnitude. The middle and bottom panels compare the analytical solution with the model results, as in Fig. 4.
In the text

Fig. 6 B-band surface brightness maps (color scale) overlaid with linear polarization maps (line segments) for inclinations of 20° (left), 75° (middle), and 90° (right column). The top row shows a model with only a stellar bulge, the middle row a model with only a stellar disk, and the bottom row a model with both stellar bulge and disk with a ratio of bulge to total luminosity of B/T = 0.5. The dust disk is the same in the three cases and has a central face-on V-band optical depth of 10.

In the text

Fig. 7 Polarization degree maps for two inclinations of a disk with high optical density (τ = 106). On the left are the linear polarization degree maps we calculated using SKIRT. The dust grain size is the same as the wavelength (1 μm), creating the intricate pattern of the polarization degree. On the right are cuts 1–6 through the maps along with results of various codes. Pinball (Watson & Henney 2001) in green, MCMax (Min et al. 2009) in black, MCFOST (Pinte et al. 2006) in blue, TORUS (Harries 2000) in dashed orange, and SKIRT in red. The data of the other codes were taken from the official website.

In the text

Fig. 8 Spiral galaxy model described in Sect. 6 and Table 2 observed at a wavelength of 1 μm. Rows from top to bottom: face-on, inclined (53°), and edge-on surface brightness (color scale) overlaid with linear polarization degree and orientation (line segments); linear polarization degree of the edge-on view, averaged over the vertical axis. Columns from left to right: the reference model; the same model observed from a different azimuth angle (rotated clockwise by 120°); a model without the spiral arm perturbations in the stellar disks, with the rotated orientation.

In the text

	Fig. 9 Linear polarization degree, averaged over the vertical axis, of the rotated model (column B of Fig. 8) observed at 1 μm, for inclinations ranging from edge-on (i = 90°) to face-on (i = 0°).
In the text

	Fig. 10 Averaged linear polarization degree for one side of the edge-on view of the rotated model, observed at optical and near-infrared wavelengths from 0.28 μm to 2.84 μm. The red vertical bands correspond to the bands in Fig. 9 and trace the tangent points of the dust spiral arms.
In the text