Polarization in Monte Carlo radiative transfer and dust scattering polarization signatures of spiral galaxies
Christian Peest, Peter Camps, Marko Stalevski, Maarten Baes, Ralf Siebenmorgen
AAstronomy & Astrophysics manuscript no. main c (cid:13)
ESO 2018September 22, 2018
Polarization in Monte Carlo radiative transfer and dustscattering polarization signatures of spiral galaxies
C. Peest (cid:63) , , P. Camps , M. Stalevski , , , M. Baes , and R. Siebenmorgen Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281 S9, B-9000 Gent, Belgium European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching b. M¨unchen, Germany Universidad de Chile, Observatorio Astronomico Nacional Cerro Calan, Camino El Observatorio 1515, Santiago, Chile Astronomical Observatory, Volgina 7, 11060 Belgrade, SerbiaPreprint online version: September 22, 2018
ABSTRACT
Polarization is an important tool to further the understanding of interstellar dust and the sources behind it. In this paper we describeour implementation of polarization that is due to scattering of light by spherical grains and electrons in the dust Monte Carlo radiativetransfer code SKIRT. In contrast to the implementations of other Monte Carlo radiative transfer codes, ours uses co-moving referenceframes that rely solely on the scattering processes. It fully supports the peel-o ff mechanism that is crucial for the e ffi cient calculationof images in 3D Monte Carlo codes. We develop reproducible test cases that push the limits of our code. The results of our program arevalidated by comparison with analytically calculated solutions. Additionally, we compare results of our code to previously publishedresults. We apply our method to models of dusty spiral galaxies at near-infrared and optical wavelengths. We calculate polarizationdegree maps and show them to contain signatures that trace characteristics of the dust arms independent of the inclination or rotationof the galaxy. Key words.
Polarization – Radiative transfer – Methods: numerical – ISM: dust, extinction – Galaxies: spiral
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 scat-tered or absorbed by dust grains before leaving the system.Scattered radiation is generally polarized. The polarization stateof the light can be used to deduce information about the grainsthat would not be available using intensity measurements alone(Scicluna et al. 2015). There are indications that dust propertiesdi ff er widely and systematically between galaxies (Fitzpatrick& Massa 1990; Gordon et al. 2003; R´emy-Ruyer et al. 2015;Dale et al. 2012; De Vis et al. 2016) and that they can vary sig-nificantly within a galaxy (Draine et al. 2014; Mattsson et al.2014). Polarimetric studies can help in constraining these prop-erties. Theoretical frameworks for modeling radiative transfertherefore usually include a section on polarization (see, e.g.,Chandrasekhar 1950; Van De Hulst 1957).Numerical simulations of dust radiative transfer most com-monly use the Monte Carlo technique (see, e.g., Whitney 2011;Steinacker et al. 2013). Codes using this method track many in-dividual photon packages as they propagate through the dustymedium, simulating emission, scattering, and absorption eventsbased on random variables drawn from the appropriate proba-bility distributions. While it is conceptually straightforward totrack the polarization state of a photon package as part of thisprocess, there are many details to be considered, and the imple-mentation complexity depends on the assumptions and approxi-mations one is willing to make. Moreover, the dust model usedby the code must provide the extra properties necessary to cal- (cid:63) [email protected] culate the changes to the polarization state for each interactionwith a dust grain.As a result, various authors have made di ff erent choicesfor implementing polarization in Monte Carlo radiative trans-fer (MCRT) codes. Most commonly, the MCRT codes consideronly scattering by spherical dust grains (e.g., Bianchi et al. 1996;Pinte et al. 2006; Min et al. 2009; Robitaille 2011; Goosmannet al. 2014). Some codes include more advanced support forpolarization by absorption and scattering o ff aligned spheroids(Whitney & Wol ff / orfor polarized dust emission (Reissl et al. 2016).To verify the correctness of the various polarization im-plementations, authors sometimes compare the results betweencodes (e.g. Pinte et al. 2009). Because of the variations in as-sumptions and capabilities, however, such a comparison is trickyand the ‘correct’ result is usually simply assumed to be the re-sult obtained by a majority of the codes. Even when the basicassumptions about grain shape and alignment as well as the dustmixture are the same and the codes support the same polarizationprocesses, comparing results is usually complicated.In this paper we present a robust framework that is indepen-dent of a coordinate system for implementing polarization in atree-dimensional (3D) MCRT code. The mathematical formu-lation and the numerical calculations in our method rely solelyon reference frames determined by the physical processes un-der 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 specialcases (i.e., a photon package propagating in the direction of the z -axis or close to it) and enables a more streamlined implemen-tation. a r X i v : . [ a s t r o - ph . I M ] F e b . Peest et al.: Polarization in RT simulations We have implemented this framework in SKIRT (Baes et al.2003; Baes et al. 2011; Camps & Baes 2015), a versatile multi-purpose Monte Carlo dust radiative transfer code. It has been de-signed and optimized for systems with a complex 3D structure,as multiple components are configured separately and constructa more complex model for the dust and / or radiation sources(Baes & Camps 2015). The code is equipped with a range ofe ffi cient grid structures on which the dust can be spatially dis-cretized, including octree, k-d tree, and Voronoi grids (Saftlyet al. 2013, 2014; Camps et al. 2013). A powerful hybrid par-allelization scheme has been developed that guarantees an op-timal speed-up and load balancing (Verstocken et al. in prep.),and it opens up a wide range of possible polarization applica-tions. In order to test the correct behavior and the accuracy ofour implementation, we have developed a number of analyticaltest cases designed to validate polarization implementations ina structured manner. Furthermore, we carefully match our po-larimetric conventions to the recommendations issued by theInternational Astronomical Union (IAU 1974).In large-scale dust systems complex geometries arise andneed to be handled by the codes. We first apply our methodto some elementary models of dusty disk galaxies, enabling aqualitative comparison with the two-dimensional (2D) modelsof Bianchi et al. (1996). We also perform the polarization partof the Pinte et al. (2009) benchmark and compare with the pub-lished results.We then implement spiral arms into dusty disk galaxy mod-els and show that this produces a marked polarimetric signaturetracing the positions of the arms. Our current implementationsupports only scattering by spherical grains. Dichroic extinctionand more complex grain shapes may also have a strong influ-ence (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 usedin this paper to describe the polarization state of electromagneticradiation, and we provide recipes for translation into other con-ventions. We then present our method and its implementationin Sect. 3, and the accompanying analytical test cases and theirresults in Sect. 4. The application of our method to benchmarktests is described in Sect. 5. The dusty spiral galaxy model is de-scribed and implemented in Sect. 6. We summarize and concludein Sect. 7.
2. Polarization
The polarization state of electromagnetic radiation is commonlydescribed by the Stokes vector, S (see, e.g., Van De Hulst 1957;Chandrasekhar 1950; Bohren & Hu ff man 1998; Mishchenkoet al. 1999), S = IQUV , (1)where I represents the intensity of the radiation, Q and U de-scribe linear polarization, and V describes circular polarization.The degrees of total and linear polarization, P and P L , can bewritten as a function of the Stokes parameters, P = (cid:112) Q + U + V / I , (2) P L = (cid:112) Q + U / I . (3) http: // k d −U +U γ el. fieldvector −Q +Q Fig. 1.
Illustration of the Stokes vector conventions recom-mended by IAU (1974) and used in this paper. The radiationbeam travels along its propagation direction, k , out of the page.The electric field vector describes an ellipse over time. The linearpolarization angle, γ , is given by the angle between the primaryaxis of the ellipse (green line) and the reference direction, d . Theposition angle of the electric field vector increases with time, thebeam has right-handed circular polarization.The (linear) polarization angle, γ , can be written as γ =
12 arctan (cid:32) UQ (cid:33) , (4)where arctan denotes the inverse tangent function that preservesthe quadrant. Combining Eqs. (3) and (4), we can also write Q = IP L cos 2 γ, (5a) U = IP L sin 2 γ. (5b)The values of Q and U depend on the polarization angle γ ,which describes the angle between the direction of linear polar-ization and a given reference direction, d , in the plane orthogonalto 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 < γ < π/ U value. A radiation beam is said to have right-handedcircular polarization (with V >
0) when the electric field vec-tor position angle increases with time, and left-handed when itdecreases.The reference direction, d , can be chosen arbitrarily as longas it is well defined and perpendicular to the propagation direc-tion. However, when the polarization state changes as a result ofan interaction (e.g., a scattering event), most recipes for prop-erly adjusting the Stokes vector require the reference directionto have a specific orientation (e.g., lying in the scattering plane).Before applying the recipe, the existing reference direction mustbe rotated about the propagation direction to match this require-ment. This is accomplished by multiplying the Stokes vector bya rotation matrix, R ( ϕ ), S new = R ( ϕ ) S . (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 R ( ϕ ) = ϕ sin 2 ϕ − sin 2 ϕ cos 2 ϕ
00 0 0 1 . (7)
2. Peest et al.: Polarization in RT simulations
To record the polarization state change for a scattering event,the Stokes vector is multiplied by the M¨uller matrix, M , corre-sponding to the event, assuming that the reference direction liesin the scattering plane (as well as in the plane orthogonal to thepropagation direction). The M¨uller matrix components dependon the geometry of the scattering event and the physical proper-ties of the scatterer, and they often depend on the wavelength. Ingeneral, the M¨uller matrix is M ( θ, λ ) = S S S S S S S S S S S S S S S S , (8)where θ is the angle between the propagation directions beforeand after the scattering event, and λ is the wavelength of theradiation. For clarity of presentation, we drop the dependenciesfrom the notation for the individual M¨uller matrix components.Including the reference direction adjustments before and afterthe actual scattering event, ϕ and ϕ new , the transformation of aStokes vector for a scattering event can thus be written as S new = R ( ϕ new ) M ( θ, λ ) R ( ϕ ) S . (9)For scattering by spherical particles, the M¨uller matrix sim-plifies to (Kr¨ugel 2002) M Sph ( θ, λ ) = S S S S S S − S S , (10)again assuming that the reference direction lies in the scatteringplane. The M¨uller matrices for a particular grain size and mate-rial can be calculated using Mie theory (see e.g., Voshchinnikov& Farafonov 1993; Bohren & Hu ff man 1998; Pe˜na & Pal 2009).For scattering by electrons, also called Thomson scattering,the M¨uller matrix is wavelength-independent and can be ex-pressed analytically as a function of the scattering angle (Bohren& Hu ff man 1998), M Th ( θ ) = cos θ + θ − θ − θ + θ
00 0 0 2 cos θ . (11) In this paper we define the Stokes vector following the rec-ommendations of the International Astronomical Union (IAU1974), as presented in Sect. 2.1 and illustrated in Fig. 1.Historically, however, authors have used various conventionsfor the signs of the Stokes parameters U and V (Hamaker &Bregman 1996, see also a recent IAU announcement ). For ex-ample, the polarization angle γ is sometimes measured whilelooking toward the observer rather than toward the source, flip-ping the sign of both U and V . Reversing the definition of cir-cular polarization handedness also flips the sign of V . Table 1lists a number of references with the conventions adopted by theauthors.Assuming that the adopted conventions are properly docu-mented, translating the values of the Stokes parameters from oneconvention into another is straightforward – by flipping the signsappropriately. When comparing or mixing formulas and recipes http: // iau.org / static / archives / announcements / pdf / ann16004a.pdf Table 1.
Conventions adopted by various authors regarding thesign of the Stokes parameters U and V relative to IAU (1974)( + U , + V ). + U − U + V IAU (1974) Chandrasekhar (1950)Martin (1974) Van De Hulst (1957)Tsang et al. (1985) Hovenier & Van der Mee (1983)*Trippe (2014) Fischer et al. (1994)*Code & Whitney (1995)*Mishchenko et al. (1999)*Gordon et al. (2001)*Lucas (2003)G´orski et al. (2005) − V Shurcli ff (1962) Bohren & Hu ff man (1998)Bianchi et al. (1996)* Mishchenko et al. (2002)* Convention indicated by citing papers with this convention constructed for di ff erent conventions, changes in the signs of U and V a ff ect the sign of the M¨uller matrix components (Eq. 8)in the third row and column and fourth row and column, respec-tively. Mathematically this can be described by multiplying theM¨uller matrix by signature matrices, M + U , + V = σ
00 0 0 ς M σ U ,ς V σ
00 0 0 ς (12) = S S σ S ς S S S σ S ς S σ S σ S S σς S ς S ς S σς S S , (13)with σ and ς being + −
1. In case of the rotation matrix,Eq. (7), this implies that the signs in front of the sine functionschange based on the chosen convention.
3. Method
SKIRT (Baes et al. 2003; Baes et al. 2011; Camps & Baes 2015)is a public multipurpose MCRT code for simulating the e ff ectof dust on radiation in astrophysical systems. It o ff ers full treat-ment of absorption and multiple anisotropic scattering by thedust, self-consistently computes the temperature distribution ofthe dust and the thermal dust reemission, and supports stochas-tic heating of dust grains (Camps et al. 2015). The code handlesmultiple dust mixtures and arbitrary 3D geometries for radia-tion 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 (Baeset 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.2016), 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 electro-magnetic radiation as a sequence of discrete photon packages.Each photon package carries a specific amount of energy (lu-minosity) at a given wavelength. A SKIRT simulation followsthe individual paths of many photon packages as they propagate
3. Peest et al.: Polarization in RT simulations nkn new ϕ k new θ Fig. 2.
Geometry of a scattering event. The angle θ is between theincoming and outgoing propagation directions k and k new . Theangle ϕ is given by the normals of the previous and the presentscattering planes n and n new .through the dusty medium. The photon package life cycle is gov-erned by various events determined stochastically by drawingrandom numbers from the appropriate probability distributions.A photon package is created at a random position based on theluminosities of the sources and is emitted in a random directiondepending on the (an)isotropy of the selected source. Dependingon the dust material properties and spatial distribution, the pho-ton package then undergoes a number of scattering events at ran-dom locations (using forced scattering; see Cashwell & Everett1959), and is attenuated by absorption along its path (using con-tinuous absorption; see Lucy 1999; Niccolini et al. 2003).To boost the e ffi ciency of the simulation and reduce thenoise levels at the simulated observers, SKIRT employs the com-mon peel-o ff optimization technique (Yusef-Zadeh et al. 1984).Rather than waiting until a photon package happens to leave thesystem under study in the direction of one of the observers, aspecial photon package is peeled o ff in the direction of eachobserver at each emission and scattering event, including anappropriate luminosity bias for the probability that a photonpackage would indeed be emitted or scattered in that direction.Meanwhile, the original or main photon package continues itstrajectory through the dust until its luminosity has become neg-ligible and the package is discarded.For the purposes of this paper, we assume that a newly emit-ted photon package represents unpolarized radiation, and its po-larization state is not a ff ected by the attenuation along its paththrough the dusty medium. We assume that the photon packageis scattered by spherical dust grains (which does a ff ect the po-larization state). This leaves us with three types of events: scat-tering the main photon package into a new direction, peeling o ff a special photon package toward a given observer, and detectinga peel-o ff photon package at an observer. As a first step towarddescribing the procedures for each of these events, we discussour approach for handling the Stokes vector reference direction. As noted in Sect. 2, the Stokes vector describing the polarizationstate of a photon package is defined relative to a given referencedirection, d , in the plane perpendicular to the propagation di-rection. We define a new direction, n , perpendicular to both the propagation direction, k , and the reference direction, d , whichare perpendicular to each other as well, so that n = k × d and d = n × k , (14)assuming all three vectors are unit vectors. By definition, thescattering plane contains both the incoming and outgoing prop-agation directions k and k new . Consider the situation before theevent (also, see Fig. 2). If the reference direction d , which isalways perpendicular to k , lies in the scattering plane as well,then n corresponds to the normal of the scattering plane. A sim-ilar situation applies after the scattering event. We store n ratherthan d with each photon package, and our procedures below aredescribed in terms of n .Some authors (e.g., Chandrasekhar 1950; Code & Whitney1995; Gordon et al. 2001) choose to rotate the Stokes referencedirection in between scattering events into the meridional planeof the coordinate system. Their procedure uses two rotation op-erations for each scattering event, one before and one after theevent, and requires special care to avoid numerical instabilitieswhen the propagation direction is close to the z -axis. The latteroccurs because the meridional plane is then ill defined.We leave the reference direction unchanged after a scatteringevent and instead perform a single rotation as part of the nextscattering 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 out-going propagation directions k and k new , defining the scatteringangle θ . After the previous scattering event, the reference direc-tion has been left in the previous scattering plane (blue), so theangle between the normals n and n new to the previous and thecurrent scattering planes determines the angle ϕ over which theStokes vector must be rotated to end up in the current scatteringplane. The transformation of the Stokes vector given in Eq. (9)can therefore be simplified to S new = M ( θ, λ ) R ( ϕ ) S . (15)Care must be taken to properly set a reference normal n fornewly emitted photon packages that have not yet experienced ascattering event. Because we assume our sources to emit unpo-larized radiation, we can pick any direction perpendicular to thepropagation direction. We postpone the details and justificationfor this procedure to Sect. 3.7. The probability that a photon package leaves a scattering eventalong a particular direction k new for an incoming direction k isgiven by the phase function, Φ ( k , k new ) ≡ Φ ( θ, ϕ ), where θ and ϕ represent the inclination and azimuthal angles of k new relativeto k , and where we omit the wavelength dependency from thenotation. In the formulation of Sect. 2, we can say that the phasefunction is proportional to the ratio of the beam intensities, I and I new , before and after the scattering event, Φ ( θ, ϕ ) ∝ I new ( θ, ϕ ) I . (16)For spherical grains, combining Eqs. (7), (10), and (15) leads to I new ( θ, ϕ ) = IS + S ( Q cos 2 ϕ + U sin 2 ϕ ) (17)and therefore, Φ ( θ, ϕ ) ∝ S + S (cid:18) QI cos 2 ϕ + UI sin 2 ϕ (cid:19) . (18)
4. Peest et al.: Polarization in RT simulations
Using Equation (5) and introducing a proportionality factor, N ,we can write Φ ( θ, ϕ ) = N S (cid:32) + P L S S cos 2( ϕ − γ ) (cid:33) . (19)The proportionality factor is determined by normalizing thephase function (a probability distribution) to unity. Integrationover the unit sphere yields N = (cid:82) π (cid:82) π ( S + P L S cos 2( ϕ − γ )) sin θ d θ d ϕ (20) = π (cid:82) π S sin θ d θ . (21) After scattering, a new direction of the photon package is de-termined by sampling random values for θ and ϕ from thephase function Φ ( θ, ϕ ). To accomplish this, we use the condi-tional probability technique. We reduce the phase function to themarginal distribution Φ ( θ ), Φ ( θ ) = (cid:90) π Φ ( θ, ϕ ) d ϕ = π N S = S (cid:82) π S sin θ (cid:48) d θ (cid:48) . (22)We sample a random θ value from this distribution through nu-merical inversion, that is to say, by solving the equation X = (cid:82) θ S sin θ (cid:48) d θ (cid:48) (cid:82) π S sin θ (cid:48) d θ (cid:48) (23)for θ , where X is a uniform deviate, that is, a random numberbetween 0 and 1 with uniform distribution.Once we have selected a random scattering angle θ , we sam-ple a random azimuthal angle ϕ from the normalized conditionaldistribution, Φ θ ( ϕ ) = Φ ( θ, ϕ ) (cid:82) π Φ ( θ, ϕ (cid:48) ) d ϕ (cid:48) (24) = π (cid:32) + P L S S cos 2( ϕ − γ ) (cid:33) , (25)where the ratio S / S depends on θ . This can again be donethrough numerical inversion, by solving the equation X = (cid:90) ϕ Φ θ ( ϕ (cid:48) ) d ϕ (cid:48) (26) = π (cid:32) ϕ + P L S S sin ϕ cos( ϕ − γ ) (cid:33) (27)for ϕ , with X being a new uniform deviate. After randomly selecting both angles θ and ϕ , we can useEq. (15) to adjust the main photon package’s Stokes vector. Wecan also calculate the outgoing propagation direction k new andthe new reference normal n new from the incoming propagationdirection k and the previous reference normal n (see Fig. 2). We use Euler’s finite rotation formula (Cheng & Gupta 1989) to ro-tate a vector v about an arbitrary axis a over a given angle β (clockwise while looking along a ), v new = v cos β + ( a × v ) sin β + a ( a · v )(1 − cos β ) . (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 aboutthe incoming propagation direction k over the azimuthal angle ϕ . Because n is perpendicular to k , we have n new = n cos ϕ + ( k × n ) sin ϕ. (29)Furthermore, the propagation direction rotates in the currentscattering plane, that is, about n new , over the scattering angle θ .Again, k is perpendicular to n new , so that k new = k cos θ + ( n new × k ) sin θ. (30) As described in Sect. 3.1, a common MCRT optimization is tosend a peel-o ff photon package toward every observer from eachscattering site. The peel-o ff photon package must carry the polar-ization state after the peel-o ff scattering event, and its luminositymust be weighted by the probability that a scattering event wouldindeed send the outgoing photon package toward the observer.To obtain this information, we need to calculate the scatteringangles θ obs and ϕ obs , given the outgoing direction of the peel-o ff scattering event, or in other words, the direction toward the ob-server, k obs . This is e ff ectively the scattering problem in reverse,in which random angles were chosen based on their probabil-ity, and the new propagation direction was calculated from theseangles.Finally, when the peel-o ff photon package reaches the ob-server, its Stokes reference direction must be rotated so that itlines 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 incomingand outgoing directions,cos θ obs = k · k obs . (31)Because 0 ≤ θ obs ≤ π , the cosine unambiguously determines theangle.To derive the azimuthal angle ϕ obs , we recall (Fig. 2) that itis the angle between the normal to the previous scattering plane, n , and the normal to the peel-o ff scattering plane, n obs . The lattercan be obtained through the cross product of the incoming andoutgoing directions, n obs = k × k obs || k × k obs || . (32)We need both cosine and sine to unambiguously determine ϕ obs in its 2 π range. We easily havecos ϕ obs = n · n obs . (33)Because k is perpendicular to both n and n obs , the following re-lation holds, sin ϕ obs k = n × n obs , (34)or, after projecting both sides of the equation on k ,sin ϕ obs = ( n × n obs ) · k . (35)
5. Peest et al.: Polarization in RT simulations
This derivation of ϕ obs breaks down for a photon package thathappens to be lined up with the direction toward the observerbefore the peel-o ff event. Indeed, in this special case of perfectforward or backward peel-o ff scattering, Eq. (32) is undefined.However, because the scattering plane does not change, it is jus-tified to simply set ϕ obs = θ obs and ϕ obs values into Eq. (15) toadjust the peel-o ff photon packages Stokes vector, and we alsoupdate the reference normal to n obs . When the photon package’spolarization state is recorded at the observer, its Stokes vectorreference direction must be parallel to the north direction, k N , inthe observer frame. This is equivalent to orienting the referencenormal 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 reasoningas for ϕ obs , so that cos α obs = n obs · k E (36)and sin α obs = ( n obs × k E ) · k obs . (37)The final adjustment to the Stokes vector is thus a rotation (seeEqs. 6 and 7) with the matrix R ( α obs ). The Stokes vector is indif-ferent to rotations by π . Using k W = − k E yields the same result. We now return to the issue of selecting a reference direction,or more precisely, a reference normal, for newly emitted photonpackages. We stated at the end of Sect. 3.2 that we can pick anydirection perpendicular to the propagation direction, because weassume our sources to emit unpolarized radiation. Indeed, it iseasily seen from Eq. (25) that the probability distribution for theazimuthal angle ϕ becomes uniform for unpolarized incomingradiation, that is, with P L,in =
0. Consequently, our choice ofreference normal in the plane perpendicular to the propagationdirection will be completely randomized after the application ofthe scattering transformation (Eq. (15)).We determine the reference normal, n = ( n x , n y , n z ), perpen-dicular to the propagation direction, k = ( k x , k y , k z ), using n x = − k x k z (cid:30)(cid:113) − k z (38) n y = − k y k z (cid:30)(cid:113) − k z (39) n z = (cid:113) − k z . (40)When k is very closely aligned with the Z -axis, the root in theseequations vanishes, and we instead select n = (1 , ,
0) as thereference direction.
4. Validation
In order to confirm the validity of our method and its implemen-tation in SKIRT, we develop four test cases for which the resultscan be calculated analytically. The analytical results are obtainedsolely using the formalisms of Sect. 2, so that taken together,the test cases verify most aspects of the procedures presented inSect. 3.The overall setup for the test cases is illustrated in Fig. 3.A central source illuminates two physically and optically thinslabs of material, which scatter part of the radiation to a dis-tant observer. The slabs are arranged on the sides of a square θ y xz θ θ TC1source slab slab ℓ ℓ ℓ
TC2 source
TC3/4source (-1,0,0) (0,-1,0)(0,1,0) (1,0,0) 021 F l u x [ a r b . un i t ] Fig. 3.
Top : Geometry used in the analytical test cases. The z -axis is toward the reader. For test case 1, a central point sourceilluminates thin slabs of electrons that are slanted by 45 ◦ towardthe observer. The central source is replaced by a small blob ofelectrons for the other test cases. The blob is illuminated by acollimated beam from within the xy plane for test case 2 andfrom below the xy plane for test cases 3 and 4. Bottom : Intensitymap of test case 1 as seen from the observer.rotated by 45 ◦ relative to the line of sight and are spatially re-solved by the observer’s instrument. To simplify the calculationswe only consider radiation detected close to the xy plane, es-sentially reducing the geometry to two dimensions. Because ofthe low scattering probability in the slabs, the path with the leastnumber of scattering events will greatly dominate the polariza-tion state at each instrument position. This allows us to deducethe dominating scattering angle, θ , corresponding to each instru-ment position along the x -axis. Referring to Fig. 3, geometricalreasoning leads to θ = π ± arctan (cid:32) | x | − (cid:33) , (41)with the plus sign for x > x <
0. In com-bination with analytical scattering properties for the slab mate-rial, it then becomes possible to derive a closed-form expressionfor the components of the Stokes vector at each position.Because the observer is considered to be at ‘infinite’ dis-tance, we can use parallel projection, and the distance from theslabs to the observer does not vary with x . However, we do needto take into account the variations in the path length (cid:96) from thecentral source to the slabs because it a ff ects the amount of radia-tion arriving at the slabs as a function of x . Geometrical reason-ing in Fig. 3 again leads to (cid:96) = (cid:112) x + (1 − | x | ) . (42)The scattering properties of the slabs and the makeup ofthe central source vary between the test cases. For test cases
6. Peest et al.: Polarization in RT simulations I n t e n s i t y ( a r b . un i t s ) . . . . . . . L i n e a r p o l a r i z a t i o n d e g r ee ( % ) − . − . . . . Test case 1: x extent − . − . − . . . . . P o l a r i z a t i o n a n g l e f r o m n o r t h ( ◦ ) − . − . . . . Test case 2: x extent − . − . − . . . . .
03 6080100120140160180 − . . ± rel. − . . ± rel. − . . ± rel. − . . ± . rel. − . . . ± . rel. − . − . . ± rel. − . − . . . . Test case 3: x extent − ± . rel. Fig. 4.
Intensity (top row), linear polarization degree (middle row), and polarization angle (bottom row) of the observed radiationin test cases 1 through 3 (left to right columns). The top section in each panel shows the analytical solution (black) and the modelresults (dashed orange). The bottom section in each panel shows the absolute di ff erences (blue) and relative di ff erences (shadedarea) of the analytic solution and the model. The green lines are calculated using twice the resolution of the θ scattering angle.1 through 3, the slabs contain electrons, with scattering ma-trix M Th given by Eq. (11). We study the observed intensity, I , the degree of linear polarization, P L , and the linear polariza-tion angle, γ , of these test cases in Sect. 4.2. Scattering by elec-trons never causes circular polarization. Therefore in Sect. 4.3we introduce slabs that contain synthetic particles with custom-designed scattering properties (test case 4). This allows us tostudy the observed circular polarization.Test case 1 has a central point source. For the remaining testcases, the central source is replaced by a small blob of electronsilluminated 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 dis-cretize certain aspects of the theoretical test setup. For our imple-mentation in SKIRT, we made the following choices. The spatialdomain of the setup is divided into 601 × ×
61 cuboidal gridcells lined up with the coordinate axes (we use odd numbers toensure that the origin lies in the center of a cell). The cells over-lapping the slabs (and where applicable, the central blob) con-tain electrons, the other cells represent empty space. The edgesof the cells and slabs are not aligned. Each detector pixel pro-vides averages over multiple cells, and the slabs are opticallythin ( τ = − along their depth). The length of the slabs is1.141, their depth is 0.006, and their height is 2 × − . The detec- tor in the observer plane has a resolution of 201 pixels along the x -axis. We use 10 photon packages for each test case to min-imize the stochastic noise characteristic of Monte Carlo codes.SKIRT uses precomputed tables of various quantities with a res-olution of 1 ◦ in θ and ϕ , to help perform the numerical inversionsin Eqs. (23) and (27). Test case 1 is designed to test the peel-o ff procedure describedin Sect. 3.6. The slabs contain electrons, and the central pointsource emits unpolarized photon packages. When a photon pack-age’s direction is toward one of the slabs, the forced interactionalgorithm of SKIRT initiates a scattering event at the slab and apeel-o ff photon package is sent toward the observer. Because theslabs are optically thin, the luminosity of the scattered originalphoton package is small. It is subsequently deleted and a newphoton package is launched.As a result, the Stokes vector of the observed photon pack-ages can be written as S TC1 = (cid:96) − R (cid:18) π (cid:19) M Th ( θ ) (1 , , , T , (43)reflecting from right to left a scattering transformation startingfrom an unpolarized state (and thus ϕ =
7. Peest et al.: Polarization in RT simulations − z e x t e n t ( a r b . un i t s ) +30 . − C i rc u l a r p o l a r i z a t i o n V / I ( % ) − . − . . . . x extent (arb. units) − . − . . . ± rel. − − Flux (Jy)
Fig. 5.
Fraction of circular polarization, V / I , of the observed ra-diation in test case 4. The top panel shows the observed surfacebrightness overlaid with arrows indicating the handedness. Thearrow length scales linearly with the magnitude. The middle andbottom panels compare the analytical solution with the modelresults, as in Fig. 4.the reference direction with the observer frame, and the depen-dency on the path length from the source to the slab. We cannow substitute Eqs. (7), (11), (41), and (42) into this equation.Because the arctangent of Eq. (41) is used as an argument forthe cosine in Eq. (11), the final equations for the intensity andthe linear polarization degree reduce to polynomials, I TC1 = x − | x | + x − | x | + , (44) P TC1L = x x − | x | + . (45)The orientation of the polarization is perpendicular to the scatter-ing plane, that is, north / south or along the y -axis in the observerframe.In test case 2 we add a scattering event to verify part of theprocedures for scattering the main photon package. To this end,we replace the central point source with a small blob of electronsat the same location, and illuminate this electron blob with acollimated beam positioned at ( − , ,
0) and oriented parallel tothe slabs toward the bottom right (see Fig. 3). In this setup, aphoton package emitted by the collimated source can reach theslabs only after being scattered by the electrons in the centralblob. This e ff ectively adds a forced first scattering event to allphoton packages reaching the observer, with a scattering anglethat can be deduced from the geometry. The peel-o ff scatteringangle is still given by Eq. (41). The scattering angle for the firstscattering event in the central electron blob is θ = θ ± (cid:18) − π (cid:19) , (46)again with the plus sign for x > x < xy -plane). The Stokes vector of the observed photon packages can thus be written as S TC2 = (cid:96) − R (cid:18) π (cid:19) M Th ( θ ) M Th ( θ ) (1 , , , T . (47)The intensity and linear polarization degree for test case 2 are I TC2 = x − | x | + x − | x | + x − | x | + , (48) P TC2L = x − | x | + x − | x | + x − | x | + x − | x | + . (49)The orientation of the polarization is north / south.In the third test case we move the collimated source belowthe xy -plane to ( − √ , − , − − √ ◦ (relative to the z -axis) and an azimuthal angle of 30 ◦ (clockwisefrom 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 aresult, the normal of the first scattering plane is tilted, while thenormal of the peel-o ff scattering plane remains aligned with the z -axis. With some trigonometry, we arrive at expressions for theangles involved in the first (main) and second (peel-o ff ) scatter-ing events, θ = arccos ± − + √ + (4 − √ | x | (cid:112) − | x | + x , (50) ϕ = ± arctan + √ (cid:112) − | x | + x − √ − | x | , (51) θ = π/ ± arctan(1 / | x | − , (52) ϕ = π/ , (53)with the plus sign for x > x <
0. Theexpression provided in Eq. (51) for ϕ is simplified and shiftedby ± π for some x . The Stokes vector is invariant under rotationsby π .The Stokes vector of the observed photon packages can nowbe written as S TC3 = (cid:96) − R ( ϕ ) M Th ( θ ) R ( ϕ ) M Th ( θ ) (1 , , , T . (54)To limit the complexity of presentation, we provide expressionsfor the Stokes parameters from which the linear polarization de-gree and angle can be calculated using Eqs. (3) and (4), I TC3 = (cid:96) (cid:104) (62 − √ x − (150 − √ | x | + (156 − √ x − (78 − √ | x | + − √ (cid:105) , (55) Q TC3 = (cid:96) (cid:104) (2 − √ x + (22 + √ | x | − (28 + √ x + (14 + √ | x | − (2 + √ (cid:105) , (56) U TC3 = (cid:96) (cid:104) (1 + √ x − (1 + √ | x | + √ (cid:105) , (57) V TC3 = . (58)Figure 4 compares the analytical solutions and SKIRT re-sults for the observed intensity, linear polarization degree, andpolarization angle for these three test cases. The intensity curvesshow a relative noise level of on average about 3%. The linearpolarization degrees are identical to below 0.1% absolute for test
8. Peest et al.: Polarization in RT simulations cases 1 and 2 and below 1% absolute for test case 3. The polar-ization angles from north are correct to below 0 . ◦ for test cases1 and 2 and below 1 ◦ for test case 3. While the linear polariza-tion degree and position angle can be determined from a simula-tion with relatively few photon packages, reducing the noise inthe intensity requires significantly more photon packages. Thisis because the number of photon packages arriving at each pixelis subject to Poisson noise, whereas the path that each photonpackage takes to the same pixel is defined within tight bound-aries. The intensity curve depends on the number of photons.The linear polarization degree and angle are independent of thenumber of photon packages. Their noise is due to the size of thepixels. Slightly di ff erent paths and scattering angles might con-tribute to the same pixel.The polarization angle for test case 1 shows an intriguingspike near x =
0. As we can see in the linear polarization de-gree curve and from Eq. (45), the radiation arriving at x = Q and U components ofthe Stokes vector are zero and the polarization angle becomesundefined (see Eq. 4). This in turn causes small numerical in-accuracies in the calculations for photon packages arriving veryclose to x = ff erences for the other quantities are generallysmaller than a fraction of a percent. In the results of test case 3the residuals of the linear polarization degree and polarizationangle contain spikes that are resolved by multiple pixels eachand symmetric with respect to x =
0. The residual of the inten-sity curve shows a similar e ff ect. It tends to be lower than ex-pected in the outer regions and higher than expected in the innerregion. These orderly deviations indicate a systematic di ff erencerather than noise. In fact, these discrepancies are caused by res-olution e ff ects in the SKIRT implementation of our method (seeSect. 4.1 for a description of our discretization choices). For ex-ample, consider the interval 0 . < x < . x -axis of the detector in theobserver frame. The corresponding interval for scattering angle θ (see Eq. 50) is 75 . ◦ < θ < . ◦ , that is, only a fractionof the 1 ◦ resolution in the SKIRT calculations related to θ . It isobvious that this lack of angular resolution relative to the outputresolution will cause inaccuracies. We calculated the residualsin the polarization degree and angle from north using twice the θ resolution and show them in pale green in Fig. 4. They con-firm that increasing the θ resolution in the calculations indeedreduces the discrepancies accordingly. To include circular polarization in test case 4, we use syntheticparticles similar to electrons, but with a scattering matrix thatmixes the U and V components of the Stokes vector: M syn ( θ ) = cos θ + θ − θ − θ + θ − θ sin θ θ sin θ θ . (59)We use the same geometry as for test case 3, replacing the elec-trons in the two slabs and in the central blob by particles de-scribed 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 I TC4 = I TC3 , (60) Q TC4 = Q TC3 , (61) U TC4 = ± (cid:96) (cid:104) (1 + √ | x | − (2 + √ x + (1 + √ | x | − √ (cid:105) , (62) V TC4 = (cid:96) (cid:104) − (1 + √ | x | + (1 + √ x − √ | x | (cid:105) . (63)In Eq. (62) the plus sign is again for x > x < V / I , forthis test case, again comparing the analytical solutions with theSKIRT results. The relative di ff erences between the analyticaland simulated results for | x | < . | x | > . θ . The pale greenline again shows the residuals when calculating with 0 . ◦ resolu-tion and has a significantly smaller residual curve. Disregardingthe outer part, the circular polarization is correct to 0.1% abso-lute.
5. Benchmark tests
We compare results of our code to Bianchi et al. (1996) as a firsttest of our implication of dust scattering polarization (rather thanjust Thompson scattering). Bianchi and collaborators describethe polarization e ff ects of scattering by spherical dust grains inmonochromatic MCRT simulations of 2D galaxy models at theB band (0.44 µ m) and I band (0.9 µ m). The models include astellar bulge, a stellar disk, and a dust disk. The stellar bulgeis described by Ja ff e (1983) (scale radius 1.86 kpc, truncated at14.8 kpc), and the stellar disk is a double-exponential disk (scalelength 4 kpc, truncated at 24 kpc; scale height 0.35 kpc, trun-cated at 2.1 kpc). The relative strength of the stellar componentsand the stellar sources is varied. In all models, the dust is as-sumed to be homogeneously distributed in a double-exponentialdisk embedded in the stellar disk, with the same horizontal pa-rameters as the stellar disk, but much thinner vertically (scaleheight 0.14 kpc, truncated at 0.84 kpc). The central face-on opti-cal depth of the dust disk τ V is a free parameter. The propertiesof the dust are taken from the MRN dust model (Mathis et al.1977), which includes graphite grains and astronomical silicateswith a size distribution n ( a ) ∝ a − . , 0 . µ m < a < . µ m.Polycyclic aromatic hydrocarbons (PAHs) and very small grainsare not included.Figure 6 displays our results, which correspond to theBianchi et al. (1996) model shown in their Figure 8. We usethe same geometry, described above, and the same bulge-to-totalratio (B / T = λ = . µ m), and optical depth( τ V = ff erent dust model. It includes a widerrange of graphite and silicate grain sizes following the Zubkoet al. (2004) size distribution, in addition to a small fraction ofPAHs as described in Camps et al. (2015). The di ff erences be-tween the dust models do not a ff ect the qualitative results at op-tical wavelengths, but do prevent a quantitative comparison.Surface brightness maps (color scale) overlaid with a linearpolarization maps (line segments) are shown in Fig. 6. The in-clinations range from nearly face-on (20 ◦ , left) to edge-on (90 ◦ ,right). The top and middle rows show models from which the
9. Peest et al.: Polarization in RT simulations − − E x t e n t ( k p c ) ◦ ◦ ◦ − − E x t e n t ( k p c ) ◦ ◦ ◦ − −
12 0 12 24
Extent (kpc) − − E x t e n t ( k p c ) ◦ − −
12 0 12 24
Extent (kpc) ◦ − −
12 0 12 24
Extent (kpc) ◦ -6 -5 -4 -3 -2 -1 S u r f a c e b r i g h t n e ss ( a r b . un i t s ) Fig. 6. B -band surface brightness maps (color scale) overlaid with linear polarization maps (line segments) for inclinations of20 ◦ (left), 75 ◦ (middle), and 90 ◦ (right column). The top row shows a model with only a stellar bulge, the middle row a model withonly 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 = V -band optical depth of 10.stellar disk and the stellar bulge, respectively, was removed. Thebottom row shows the B / T = ◦ and slightly decreases again when approaching 90 ◦ . The polar-ization degree averaged over all pixels is below 1% for all in-clinations. For pure bulges, the maximum polarization degree islargely independent of the inclination. In the outer regions of thedust disk, the linear polarization degree is 22%. In the mixedB / T = I band ( λ = . µ m) with a color-adjusted bulge luminosity ( B − I = To test the performance of our code on a problem of a di ff erentscale, we compute polarization results of Pinte et al. (2009). Init, a thick dusty disk surrounds a central star. The star extendsout to 2 AU and has a temperature of 4000 K. The dust consistsof spherical grains with a radius of 1 µ m, and the light has awavelength of 1 µ m as well. The dust density distribution ρ iscylindrical, ρ ( R , z ) = Σ R (cid:32) RR (cid:33) − / exp − (cid:32) zh (cid:33) (cid:32) RR (cid:33) − / , (64)with a surface density Σ , scale radius R =
100 AU, radial dis-tance from the center R , vertical distance from the midplane z , and scale height h =
10 AU. The disk is truncated at R min = . R max =
400 AU. The surface density depends on thetotal dust mass m = × − M (cid:12) by Σ = m / (cid:34)
125 (2 π ) / h R (cid:16) ( R min / R ) / − ( R max / R ) / (cid:17)(cid:35) (65)and the albedo of the dust is 0.6475 and the opacity 4752 cm / gat 1 µ m. We adopt the scattering matrix as provided by Pinte
10. Peest et al.: Polarization in RT simulations i =69 . °
123 TORUSMCFOSTMCMaxPinballSKIRT i =87 . °
456 0% 25% 50% 75% 100%Polarization degree 0% 20% 40% 60% 80%100% P o l a r i z a t i o n d e g r ee D i ff e r e n c e P o l a r i z a t i o n d e g r ee D i ff e r e n c e Fig. 7.
Polarization degree maps for two inclinations of a disk with high optical density ( τ = ). On the left are the linearpolarization degree maps we calculated using SKIRT. The dust grain size is the same as the wavelength (1 µ m), creating the intricatepattern 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 dashedorange, and SKIRT in red. The data of the other codes were taken from the o ffi cial website.et al. (2009). The system is resolved at a distance of 140 pc bya detector with 251 ×
251 pixels covering 900 ×
900 AU .Figure 7 shows the linear polarization maps calculated byour code for inclinations of 69.5 ◦ and 87.1 ◦ . The flux di ff erencefrom the borders to the center area is about 17 orders of mag-nitude. The maps are displayed in gray outside the truncationradius, where no flux was recorded. The intricate pattern of thepolarization degree is a result of the uniform grain radius beingequal to the wavelength. The phase function therefore containsresonances and steep gradients for small changes of the scat-tering angle. We compare our results to the results of the fourpolarization capable codes in Pinte et al. (2009) along six cutsthrough the maps. The first and third row show the polarizationdegree along the cuts, and the second and fourth row show thedi ff erence of the codes to the average of all results. The TORUS http: // ipag.osug.fr / pintec / benchmark / code (Harries 2000) is not included in calculating the averagesbecause its signal-to-noise ratio is too low. In the central area thecodes agree within 10% of absolute polarization, but as the trueresult is unknown, a quantitative analysis of the results of thisbenchmark is di ffi cult. In general, the results of the SKIRT codeare close to the average result, and SKIRT seems to agree par-ticularly well with the Pinball code (Watson & Henney 2001).Pinball employs some of the same optimization techniques thatSKIRT uses (e.g., forced interaction and peel-o ff , named “forcedescape” in their paper).
6. Application: spiral galaxy models
We study the polarization properties of a 3D galaxy model in-cluding spiral arms to investigate how the spiral arm structure isimprinted in the polarization structure. We also study this in the
11. Peest et al.: Polarization in RT simulations − − F a c e - o n v i e w ( k p c )
2% 2%2% − − I n c li n e d v i e w ( k p c )
2% 2%2% − − E dg e - o n v i e w ( k p c )
2% 2% − − Dusty spiral: x extent (kpc) . . . . . . . . E dg e - o n p o l a r i z a t i o n d e g r ee ( % ) − − Rotated dusty spiral: x extent (kpc) − − Stellar Disk: x extent (kpc) 10 -1 S u r f a c e b r i g h t n e ss ( M J y / s r ) 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; thesame model observed from a di ff erent azimuth angle (rotated clockwise by 120 ◦ ); a model without the spiral arm perturbations inthe stellar disks, with the rotated orientation.edge-on view when the structure is not easily characterized fromthe intensity alone.We still assume a homogeneous distribution for the stellarsources and the dust. The model includes a stellar bulge anddisk with an older star population, a second stellar disk with ayounger star population, and a dust disk. The relevant parametersare listed in Table 2. The bulge is modeled by a spheroidal den-sity distribution obtained by flattening a spherical S´ersic profile(S´ersic 1963), as implemented by Baes & Camps (2015). Thedistributions of the two stellar and the dust disks are truncateddouble-exponential disks. The spiral arm structure is introducedby adding a perturbation to the overall density profile, as pre-sented by Baes & Camps (2015). We have made the spiral armsin the older stellar population ‘lead’ those in the younger stellar population and in the dust disk by varying the spiral arm phasezero-points. The emission of the stellar populations is modeledas blackbody spectra at the indicated temperatures. We use thedust model of Sect. 5.1 (Camps et al. 2015). The total dust massis given in Table 2, the central face-on B -band optical depth ofthe dust disk is τ V ≈ . µ m. The top row shows the modelface-on, the second row inclined (53 ◦ ), and the third row edge-on. The polarization degree is up to 1% around the central partof the models and over the whole map the average polariza-tion degree is similar to the average polarization degree from theB / T = / T =
12. Peest et al.: Polarization in RT simulations
Fig. 9.
Linear polarization degree, averaged over the verticalaxis, of the rotated model (column B of Fig. 8) observed at 1 µ m,for inclinations ranging from edge-on (i = ◦ ) to face-on (i = ◦ ). Fig. 10.
Averaged linear polarization degree for one side of theedge-on view of the rotated model, observed at optical and near-infrared wavelengths from 0.28 µ m to 2.84 µ m. The red verticalbands correspond to the bands in Fig. 9 and trace the tangentpoints of the dust spiral arms.entation of the polarization is circular around the central bulge,and for the inclined view the polarization degree left and rightof the center increases, while it decreases behind and in frontof it. In contrast to the azimuthally uniform model, there is aspiral structure in the polarization map. The linear polarizationdegree is higher in the arm regions and disappears in the inter-arm region. The maximum polarization degree is slightly inwardfrom the regions of the arms with the highest flux. The panelin the bottom row plots the linear polarization degree for the Table 2.
Parameters for our spiral galaxy model, including theblackbody temperature and the bolometric luminosity of the stel-lar components, and the total dust mass in the dust component.
S´ersic profile BulgeS´ersic index 3E ff ective radius 1.6 kpcFlattening parameter 0.6Temperature 3500 KLuminosity 3 × L (cid:12) Exponential disks Old stars Young stars DustScale length 4 kpc 4 kpc 4 kpcHorizontal cuto ff
20 kpc 20 kpc 20 kpcScale height 0.35 kpc 0.2 kpc 0.2 kpcVertical cuto ff ◦ ◦ ◦ Radius zero-point 4 kpc 4 kpc 4 kpcPhase zero-point 0 ◦ ◦ ◦ Perturbation weight 0.25 0.75 0.75Arm-interarm size ratio 5 5 5Temperature 3500 K 10000 KLuminosity / mass 4 × L (cid:12) × L (cid:12) × M (cid:12) edge-on view, averaged over the vertical axis. This average isobtained by summing each individual component of the Stokesvector and calculating the polarization degree from these totals.Regions with higher linear polarization (up to 2%) clearly tracethe spiral arms and are prominent at all inclinations, includingthe edge-on view. The maxima in the polarization signature ofthe edge-on view match the positions of the spiral arms alongthe line of sight.To verify this, the middle column of Fig. 8 shows the samemodel from a di ff erent azimuth angle. The peaks in the polar-ization 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 armsperturbations from the stellar disks in the model. The polariza-tion signature remains essentially unchanged; the maxima areslightly higher (by a factor of up to 1.2), but the structure isthe same. The signature could also be produced by the di ff erentphase zero-points of the old stars and the dust (see Table 2). Wecalculated results using the same phase zero-point for all compo-nents (not shown here). The outer maxima are lower (by a factorof 0.8), while the inner maxima are unchanged. This confirmsthat the polarization signature is created by the distribution ofthe dust and not by the distribution of the sources.In Fig. 9 we further study the e ff ect of inclination on the ob-served polarization signature for the reference model. In the cen-tral region ( r (cid:46) ◦ . Outside this region, however, theform of the curves is very similar for all inclinations. Althoughthe polarization degree generally decreases toward lower incli-nations, the peaks remain at the spiral arm tangent points, andthe ratio between the maximum and minimum polarization de-gree remains roughly stable at a factor of about 2.In Fig. 10 we compare the edge-on polarization signature ofour reference model for various optical and near-infrared wave-lengths. The polarization peaks remain prominent over the wave-length range 0 . µ m (cid:46) λ (cid:46) µ m. At shorter wavelengths thegeneral polarization degree is higher and the signature is re-
13. Peest et al.: Polarization in RT simulations versed, light from within the arms is slightly less polarized thanlight from the inter-arm regions. The increased interaction crosssection causes the inter-arm dust to become e ffi cient at scatter-ing the stellar radiation, boosting the polarization degree. Wefind that in the spiral arms the ratio of once scattered to multi-ple times scattered light is 2.5:1, while in the inter-arm regionit is 3.3:1. The polarization orientation after multiple scatteringsis less uniform, which lowers the polarization degree in the armregions.At longer wavelengths, the signature retains the same form,but the reduced scattering e ffi ciency of the dust causes the polar-ization degree to be very low, so that the peaks become hard todiscern.Our results imply that polarization measurements could beused, at least in principle, to study the spiral structure of edge-onspiral galaxies, where intensity measurements alone have limiteddiagnostic power. The contrast would be highest at around 1 µ mand with an expected polarization degree of up to 0.6%, this iswell within the capabilities of current polarization capable tele-scopes.We note that the polarization degree of the edge-on galaxyNGC 891 was mapped by Montgomery & Clemens (2014) at1.6 µ m. They found polarization degrees of below 1% that variedalong the disk profile. We expect the orientation of polarizationdue to scattering to be perpendicular to the disk. Montgomery& Clemens (2014) found the orientation to be rather parallel tothe disk and attributed most of it to dichroic extinction. Our codedoes not yet support dichroic extinction, so we cannot comparethe strength of these two e ff ects.
7. Conclusion
We presented a robust framework that is independent of a coordi-nate system for implementing polarization in a 3D MCRT code,focusing on scattering by spherical dust grains. The mathemat-ical formulation and the numerical calculations in our methodrely solely on the scattering planes determined by the physicalprocesses rather than by the coordinate system. This approachavoids numerical instabilities for special cases and enables amore streamlined implementation. We described four test caseswith well-defined geometries and material properties, yieldinganalytical solutions. These setups are designed to validate thecalculated results in a structured manner, and they can serve asbenchmarks for other implementations as well.We reconstructed a selection of the 2D models of Bianchiet al. (1996), confirming that our implementation reproducestheir results at least qualitatively. A quantitative comparison isnot possible because of di ff erences in the dust model. We thencalculated results for the polarization part of the Pinte et al.(2009) benchmark and obtained similar results as the other codesthat took part in it. As an application of our code we constructeda 3D spiral galaxy model including a stellar bulge and disk withan older star population, a second stellar disk with a youngerstar population, and a dust disk. The stellar and dust distribu-tions feature an analytical spiral arm perturbation. We showedthat scattering of light at the dust in the spiral arms producesa marked polarimetric signature. It traces the tangent positionsof the arms for wavelengths in the range 0 . µ m (cid:46) λ (cid:46) µ m,regardless of inclination.It is fair to note, however, that our current implementationis limited to scattering by spherical dust grains. We plan to addsupport for polarized emission and for scattering and absorp-tion by (partially) aligned spheroidal grains in future work. Withthe relevant physics included, we can also study the influence of changes to dust properties on the polarization signature, whichwe have not addressed in this paper. Acknowledgements.
C.P., P.C., and M.B. acknowledge the financial sup-port from CHARM (Contemporary physical challenges in Heliospheric andAstRophysical Models), a Phase-VII Interuniversity Attraction Pole programorganized by BELSPO, the BELgian federal Science Policy O ffi ce. M.S. ac-knowledges support by FONDECYT through grant no. 3140518 and by theMinistry of Education, Science and Technological Development of the Republicof Serbia through the projects Astrophysical Spectroscopy of ExtragalacticObjects (176001) and Gravitation and the Large Scale Structure of the Universe(176003). We thank Rene Goosmann, Francesco Tamborra, and Frederic Marinfor many useful discussions and providing insights on the operation of theSTOKES code. We wish to thank the anonymous referee for their constructiveinput to this paper. References
Baes, M. & Camps, P. 2015, Astronomy and Computing, 12, 33Baes, M., Davies, J. I., Dejonghe, H., et al. 2003, MNRAS, 343, 1081Baes, M., Fritz, J., Gadotti, D. A., et al. 2010, A&A, 518, L39Baes, M., Verstappen, J., De Looze, I., et al. 2011, The Astrophysical JournalSupplement Series, 196, 22Bianchi, S., Ferrara, A., & Giovanardi, C. 1996, ApJ, 465, 127Bohren, C. F. & Hu ff man, D. R. 1998, Absorption and scattering of light bysmall particles (New York (N.Y.) : Wiley-Interscience)Camps, P. & Baes, M. 2015, Astronomy and Computing, 9, 20Camps, P., Baes, M., & Saftly, W. 2013, A&A, 560, A35Camps, P., Misselt, K., Bianchi, S., et al. 2015, Astronomy & Astrophysics, 580,A87Camps, P., Trayford, J. W., Baes, M., et al. 2016, MNRAS, 462, 1057Cashwell, E. & Everett, C. 1959, A practical manual on the Monte Carlo methodfor random walk problems, International tracts in computer science and tech-nology and their application (Pergamon Press)Chandrasekhar, S. 1950, Radiative transfer. (Clarendon Press, Oxford)Cheng, H. & Gupta, K. 1989, Journal of Applied Mechanics, 56, 139Code, A. D. & Whitney, B. A. 1995, ApJ, 441, 400Dale, D. A., Aniano, G., Engelbracht, C. W., et al. 2012, The AstrophysicalJournal, 745, 95De Geyter, G., Baes, M., Camps, P., et al. 2014, MNRAS, 441, 869De Geyter, G., Baes, M., De Looze, I., et al. 2015, MNRAS, 451, 1728De Looze, I., Baes, M., Bendo, G. J., et al. 2012, MNRAS, 427, 2797De Looze, I., Fritz, J., Baes, M., et al. 2014, A&A, 571, A69De Vis, P., Dunne, L., Maddox, S., et al. 2016, Monthly Notices of the RoyalAstronomical Society, stw2501Deschamps, R., Braun, K., Jorissen, A., et al. 2015, A&A, 577, A55Draine, B., Aniano, G., Krause, O., et al. 2014, The Astrophysical Journal, 780,172Draine, B. T. & Fraisse, A. A. 2009, The Astrophysical Journal, 696, 1Fischer, O., Henning, T., & Yorke, H. W. 1994, A&A, 284, 187Fitzpatrick, E. L. & Massa, D. 1990, The Astrophysical Journal SupplementSeries, 72, 163Goosmann, R. W. & Gaskell, C. M. 2007, A&A, 465, 129Goosmann, R. W., Gaskell, C. M., & Marin, F. 2014, Advances in SpaceResearch, 54, 1341Gordon, K. D., Clayton, G. C., Misselt, K., Landolt, A. U., & Wol ff , M. J. 2003,The Astrophysical Journal, 594, 279Gordon, K. D., Misselt, K., Witt, A. N., & Clayton, G. C. 2001, TheAstrophysical Journal, 551, 269G´orski, K. M., Hivon, E., Banday, A. J., et al. 2005, The Astrophysical Journal,622, 759Hamaker, J. & Bregman, J. 1996, Astronomy and Astrophysics SupplementSeries, 117, 161Harries, T. J. 2000, Monthly Notices of the Royal Astronomical Society, 315,722Hendrix, T., Keppens, R., & Camps, P. 2015, A&A, 575, A110Hendrix, T., Keppens, R., van Marle, A. J., et al. 2016, MNRASHovenier, J. & Van der Mee, C. 1983, Astronomy and Astrophysics, 128, 1IAU. 1974, Transactions of the IAU, ed. G. Contopoulos & A. Jappel, Vol. 15B(Springer Netherlands), 166Ja ff e, W. 1983, Monthly Notices of the Royal Astronomical Society, 202, 995Kr¨ugel, E. 2002, The Physics of Interstellar Dust, Series in Astronomy andAstrophysics (CRC Press)Lucas, P. 2003, Journal of Quantitative Spectroscopy and Radiative Transfer, 79,921, electromagnetic and Light Scattering by Non-Spherical ParticlesLucy, L. B. 1999, A&A, 344, 282Martin, P. 1974, The Astrophysical Journal, 187, 461
14. Peest et al.: Polarization in RT simulations
Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, The Astrophysical Journal,217, 425Mattsson, L., Gomez, H. L., Andersen, A., et al. 2014, Monthly Notices of theRoyal Astronomical Society, 444, 797Min, M., Dullemond, C. P., Dominik, C., de Koter, A., & Hovenier, J. W. 2009,A&A, 497, 155Mishchenko, M. I., Hovenier, J. W., & Travis, L. D. 1999, Light scattering bynonspherical particles: theory, measurements, and applications (Academicpress)Mishchenko, M. I., Travis, L. D., & Lacis, A. A. 2002, Scattering, absorption,and emission of light by small particles (Cambridge university press)Montgomery, J. D. & Clemens, D. P. 2014, The Astrophysical Journal, 786, 41Mosenkov, A. V., Allaert, F., Baes, M., et al. 2016, A&A, 592, A71Niccolini, G., Woitke, P., & Lopez, B. 2003, A&A, 399, 703Pe˜na, O. & Pal, U. 2009, Computer Physics Communications, 180, 2348Pinte, C., Harries, T., Min, M., et al. 2009, Astronomy & Astrophysics, 498, 967Pinte, C., M´enard, F., Duchˆene, G., & Bastien, P. 2006, A&A, 459, 797Reissl, S., Brauer, R., & Wolf, S. 2016, arXiv preprint arXiv:1604.05305R´emy-Ruyer, A., Madden, S., Galliano, F., et al. 2015, Astronomy &Astrophysics, 582, A121Robitaille, T. P. 2011, A&A, 536, A79Saftly, W., Baes, M., & Camps, P. 2014, A&A, 561, A77Saftly, W., Baes, M., De Geyter, G., et al. 2015, A&A, 576, A31Saftly, W., Camps, P., Baes, M., et al. 2013, A&A, 554, A10Scicluna, P., Siebenmorgen, R., Wesson, R., et al. 2015, Astronomy &Astrophysics, 584, L10S´ersic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La PlataArgentina, 6, 41Shurcli ff , W. A. 1962, Harvard University Press, Cambridge, Mass, 165Siebenmorgen, R., Voshchinnikov, N., & Bagnulo, S. 2014, Astronomy &Astrophysics, 561, A82Stalevski, M., Fritz, J., Baes, M., Nakos, T., & Popovi´c, L. ˇC. 2012, MNRAS,420, 2756Stalevski, M., Ricci, C., Ueda, Y., et al. 2016, MNRAS, 458, 2288Steinacker, J., Baes, M., & Gordon, K. 2013, ArXiv e-printsTrippe, S. 2014, arXiv preprint arXiv:1401.1911Tsang, L., Kong, J. A., & Shin, R. T. 1985, Theory of microwave remote sensingVan De Hulst, H. C. 1957, Light scattering by small particles (CourierCorporation)Verstocken, S., Camps, P., & Baes, M. in prep.Viaene, S., Baes, M., Tamm, A., et al. 2016, arXiv preprint arXiv:1609.08643Voshchinnikov, N. 2012, Journal of Quantitative Spectroscopy and RadiativeTransfer, 113, 2334Voshchinnikov, N. & Farafonov, V. 1993, Astrophysics and Space Science, 204,19Watson, A. M. & Henney, W. J. 2001, arXiv preprint astro-ph / ff , M. J. 2002, The Astrophysical Journal, 574, 205Yusef-Zadeh, F., Morris, M., & White, R. L. 1984, ApJ, 278, 186Zubko, V., Dwek, E., & Arendt, R. G. 2004, ApJS, 152, 211, M. J. 2002, The Astrophysical Journal, 574, 205Yusef-Zadeh, F., Morris, M., & White, R. L. 1984, ApJ, 278, 186Zubko, V., Dwek, E., & Arendt, R. G. 2004, ApJS, 152, 211