EDP Sciences
Free Access
Issue
A&A
Volume 593, September 2016
Article Number A87
Number of page(s) 17
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201424930
Published online 27 September 2016

© ESO, 2016

1. Introduction

It is well established that giant molecular clouds within the interstellar medium (ISM) are the birthplace of stars and later planetary systems. However, understanding star formation remains one of the outstanding remaining challenges of modern astrophysics. In spite of significant progress in recent years (Crutcher 2010), fundamental questions about the basic physics of star formation remain unanswered.

The importance of magnetic fields in regulating the star formation process is still a topic of ongoing research. Self-gravitating clouds have long been believed to reduce support against collapse by magnetic support and triggering star formation (Mouschovias & Ciolek 1999) with time scales of collapse of up to 10 Myr, assuming ambipolar diffusion (Falgarone et al. 2008; Girart et al. 2009; Crutcher et al. 2010) indicating a significant influence of magnetic fields on the star formation process. Observations, however, suggest that molecular clouds are shorter lived and the cloud collapse is controlled by turbulence with negligible magnetic contribution (Elmegreen 2000; Klessen et al. 2004).

Observations of magnetic fields in molecular clouds are essential for understanding their role in the evolution from dense clouds to the formation of stars. One way to improve our knowledge about the involvement of magnetic fields in the star formation process is to trace its morphology with the help of dust grains. In addition to magnetic fields, dust grains also play a central role in the ISM – from heating and cooling processes, astro-chemistry, to the polarization effects of light (Das et al. 2010). Non-spherical rotating dust particles have the tendency to align with their shorter axis along the local magnetic field direction (Martin 1974; Voshchinnikov 2012). Therefore, mapping of light that is polarized by elongated and aligned dust grains in the near-infrared to sub-mm (Lucas et al. 2004; Gonçalves et al. 2005; Frau et al. 2011) offers one way to obtain information about the magnetic field morphology in regions associated with star formation. Additionally, at optical to near infrared wavelengths, light from embedded stars scattered as dust grains enables us to investigate the spatial structure of star-forming regions and can also constrain the dust properties (Lucas & Roche 1997).

However, to make comparisons between star formation theories with observational data, one requires an expedient understanding of the intrinsic radiative transfer (RT), the dust properties, and the physics of grain alignment. Unfortunately, the 3D dust RT problem is nonlinear, which makes it one of the hardest challenges in computational astrophysics. In idealized scenarios, the problem can be solved by using analytical radiative transfer techniques (e.g. Ivezic et al. 1997; van Zadelhoff et al. 2002; Steinacker et al. 2006). Environments with ongoing star formation, however, are regions where gas and dust density and magnetic fields form complex structures and their physical quantities cover several orders of magnitude. The Monte-Carlo (MC) method is a powerful approach to perform radiative transfer simulations independent of complexity. Here, nature is mimicked by sending photon packages along probabilistic pathways to obtain numerical solutions to RT problems. A broad variety of codes has been developed to cover problems of astrophysical importance such as the heating of dust, photoionization, and the calculation of polarization of scattered radiation (Juvela 1999; Wolf et al. 1999; Whitney & Wolff 2002; Niccolini et al. 2001; Ercolano et al. 2003; Wolf 2003; Steinacker & Henning 2003; Min et al. 2009; Whitney 2011; Baes et al. 2011; Dullemond 2012; Robitaille 2013; Harries 2014).

So far, vital questions concerning dust composition, alignment efficiency, and size distribution are unknown. With forthcoming observational equipment, e.g. SOPHIA/HAWC+ Dowell et al. (2013) or ALMA Brown et al. (2004), the need for accurate and expedient modeling of observational polarimetry data is increasing. Hence, we developed a new RT code from scratch and fully incorporated radiative transfer, dust heating, and polarization algorithms, as well as state-of-the-art dust grain alignment theories. Here, we implemented the classical imperfect Davis-Greenstein (IDG) alignment owing to paramagnetic relaxation (Aannestad & Greenberg 1983; Jones & Spitzer 1967), the radiative torque alignment (RAT) because of radiation-dust interaction (Hoang & Lazarian 2007a), and considered the effects of mechanical alignment (GOLD; Gold 1952; Lazarian 1994, 1995) as a result of gas streams.

This paper is structured as follows. First, we present the POLARIS code and its basic RT and dust heating algorithm in Sect. 2. Then, we briefly discuss the applied dust grain model and the implemented polarization and grain alignment features in Sect 2.2. The optimization methods that were considered in the POLARIS code are listed in Sect. 2.5. We describe the challenges of the scattering process on partially aligned dust grains in Sect 2.3. To demonstrate the functionality and predictability of POLARIS, we calculate synthetic intensity and polarization maps in Sect. 3. Here we considered both analytical models and complex magnetohydrodynamic (MHD) simulations of different stages of the star formation process. Finally, we summarize our results and give a short outlook to future projects in Sect. 4.

2. The code

POLARIS is a three-dimensional Monte-Carlo (MC) continuum radiative transfer (RT) code dedicated to post-process MHD simulations to determine the observability of the physical quantities. The code can make use of the full parameter set (density, temperature, velocity, and magnetic field distributions) delivered by MHD simulations to perform dust heating and complex polarized radiative transfer (RT) calculations. All RT simulations are performed on an adaptive octree grid, which is optimized to maintain the data structure of certain MHD codes, e.g. FLASH Fryxell et al. (2010). The grid allows a time efficient propagation of photon packages and provides a dynamical spatial resolution to resolve regions of high optical depth and simultaneously reduces the amount of required memory. The code is platform independent and completely written in C++ in a strict object orientated architecture making it easy to develop classes for future purposes such as photoionization, line transfer, or to implement further grid geometries.

The code works in three independent main operations modes:

  • 1.

    MC simulations of self-consistent heating of the dust grains byan arbitrary number of radiation sources.

  • 2.

    MC calculation of the mean energy density and anisotropy by the radiation field for radiative torque (RAT) alignment (see Sect. 2.2.4).

  • 3.

    Calculation of spectral energy distributions (SEDs) and synthetic multi-wavelength intensity and polarization maps for statistical analysis.

To access its full functionality, the code POLARIS parses script files. Predefined commands in a pseudo XML-style that is implemented in the code allows the user to create scripts with sequences of operation modes.

2.1. Radiative transfer

Monochromatic photon packages with a fixed energy ϵ0 from distinct sources are emitted into the model space. To cover the various radiation-dust interactions in the ISM, we provide classes of point (stars), diffuse (star field), thermal dust re-emission and external (interstellar radiation field and background objects) sources of radiation. Each source has its characteristic spatial emission, energy spectrum, and degree of linear and circular polarization. A higher amount of photon packages leads to an increase in calculation time but also results in a better S/N inherent in the MC process.

In the dust heating mode we use a combined algorithm of continuous absorption (Lucy 1999) and immediate temperature correction (Bjorkman & Wood 2001). First, the mass specific emissivity (1)of the entire ensemble of dust grain sizes is pre-calculated and interpolated as a function of dust temperature for a user-defined temperature range using pre-tabulated dust grain properties (see Sect. 3). Here, Bλ(T) is the Planck function, Td is the temperature of the dust, and Cabs the cross-section of absorption. Afterwards, the photon package propagates a distance li in the grid between two points of radiation-dust interactions. At each cell wall, the absorption rate per time unit (2)is updated as presented in Lucy (1999), assuming constant temperatures and densities in each cell. We extended the presented algorithm by an additional therm ΔĖ to incorporate temperatures provided by MHD simulations into the RT simulations (see Sect. 2.4 for detail). Since we assume local thermodynamic equilibrium (amount of absorbed energy equals amount of re-emitted energy), the updated dust temperature is provided by solving (3)The optical depth τ is integrated over the passed distance l up to the next interaction (τ>τmax). The interaction distance is determined by τmax = −log 10(1−z) (see Wolf 2003). In this paper, z donates a randomly generated number with z ∈ [0;1 [. At each point of interaction the scattering-absorption probability is determined by the albedo Csca/Cext. In the case of scattering (Csca/Cext<z), a new direction (ϑ′,ϕ′) is sampled from the phase function of the dust grain. POLARIS uses optionally isotropic scattering, full Mie-scattering or, alternatively, the inverted Henyey-Greenstein phase function (4)with (5)as presented in Henyey & Greenstein (1941). For the latter, the probability distribution of scattering is described by one parameter, the anisotropy g = ⟨cos(ϑ′)⟩ ∈ [−1,1], g = 0 means isotropic scattering, where g = −1 backward scattering, and g = 1 forward scattering, respectively. The anisotropy parameter g and the full scattering matrix (see Sect. 2.3) are pre-calculated with DDSCAT (Draine & Flatau 2013, see Sects. 2.2.4 and 3) and is a function of grain size, grain orientation, and wavelength. Although alternative phase functions are currently discussed (see Sharma & Roy 2008), the HG phase function still provides a good approximation for the scattering probability. During all dust-heating calculations, the dust grain orientation is assumed to be random.

The dust temperature Td is the essential parameter for thermal radiation. Hence, the case of an absorption and immediate thermal re-emission event is treated by the instantaneous temperature correction technique of Bjorkman & Wood (2001). A photon package absorbed by dust will be immediately re-emitted, assuming a blackbody spectrum. The new wavelength λi of the thermal dust photon emitted in the ith cell is determined by (6)to ensure thermal equilibrium. This process is repeated until the photon package reaches the boundary of the grid.Once the correct dust temperature is achieved by the dust heating algorithm, intensity and polarization maps can be calculated.

The multi-wavelength MC-RT polarization simulation, including scattering, are calculated similarly to dust heating and can be performed in a monochromatic or multi-wavelength mode. In the case of polarization calculations, the dust temperature now remains constant and different dust-grain alignment mechanisms can be applied.

2.2. Dust modeling and alignment

Interstellar dust grains consist of silicate and carbon compounds and make up only about 1% of total mass of the ISM. However, the extinction of dust grains can lead to significant degrees of linear and circular polarization from the optical to millimeter regime because of scattering, dichroic extinction, and thermal re-emission.

POLARIS considers the dust as a mixture of different dust components. Each component consist of its own material, size distribution, sublimation temperature, and alignment behavior. When the dust temperature exceeds the sublimation temperature of a certain material, this component will be removed from the cell. For the mixing of different dust materials, the distinct cross-sections can simply be multiplied with their ratio of abundance and added up (see e.g. Das et al. 2010).

Non-spherical spinning dust grains align with their shorter axis, which is parallel to the local magnetic field direction. Thus, they preferentially block light perpendicular to the magnetic field direction, which results in a polarization parallel to the magnetic field while re-emitted light is polarized perpendicularly. The radiation emitted by dust will mostly be found in the mid-infrared, for a given dust temperature. Hence, the orientation of polarized light is related to the direction of the projected underlying magnetic field morphology, leading to a characteristic polarization pattern.

The implemented RT equations for extinction and polarization (Whitney & Wolff 2002) require the sum and difference of cross-sections (C = πa2 × Q) parallel Cext, | | and perpendicular Cext, ⊥ to the direction of the magnetic field (see Reissl et al. 2014, for detail). Here, a is the effective radius of the spherical grain with equivalent volume to that of a spheroid grain and Q is the efficiency of interaction.

We use the Stokes vector S = (I,Q,U,V)T to describe the change in intensity and degree of polarization where I is intensity, Q and U quantify the linear polarization, and V the circular polarization. Writing the radiative transfer equation in the Stokes vector formalism allows to calculate the effects of dichroic extinction, thermal re-emission (Martin 1974) with: (7)where Td and nd are the temperature and number density, respectively, of the dust. In this paper denotes the averaging over grain size distribution and X over orientation. So far, POLARIS is optimized for oblate dust grains but adjustments for arbitrary shapes can easily be applied. In Eq. (7), the required matrix of cross-sections depends on grain size, orientation, and wavelength. The matrix elements are determined by the values along the minor axis (Cext, | |) and major axis (Cext, ⊥), and are calculated for each dust grain size in each RT simulation with (8)(9)The angle ϑ is between incident light and magnetic field direction (see Fig. 1). For oblate dust grains and (12)The Rayleigh reduction factor R (Greenberg 1968) is introduced to take care of imperfect grain alignment and is defined as (13)Here, β is the alignment cone angle between the angular momentum J and magnetic field B and ζ is the internal alignment angle between the axis of the largest moment of inertia I|| and J, while G(x) = 1.5x−0.5. The Rayleigh reduction factor is defined as R ∈ [−0.5;1] , where positive values correspond to an alignment with the longer grain axis that is perpendicular to the magnetic field direction and vice versa for negative values. We note that the Rayleigh reduction factor can also be a function of wavelength and effective dust grain radius. Consequently, grain alignment and, subsequently, polarization is completely determined by first order moments weighted over a distribution function defined by . Unfortunately, alignment and internal alignment do not work independently. For exact solutions, a simultaneous integration of both distribution functions f(β) and f(ζ) over cone angle β and angle ζ, respectively, is required. Depending on the considered grain alignment mechanism, the distribution functions, in turn, are also functions of density, temperature, magnetic field strength, velocity, and direction of the incident light (see Sects. 2.2.2–2.2.4). Because of the enormous parameter space, the reduction in polarization cannot be pre-calculated without losing precision. The exact calculations that are required, however, for each photon package-dust interaction results in a further burden to 3D MC-RT simulations. Therefore, in POLARIS we approximate (14)with a correlation factor fc where fc = 0 stands for no correlation (Hoang et al. 2013, 2014). The correlation factor is of minor influence to the net polarization, meaning that this approach is well justified.

thumbnail Fig. 1

Geometrical configuration of a dust grain partially aligned with its angular momentum J to the magnetic field direction B. The precession of J around B defines the cone angle β, and the precession of maximum moment of inertia I around J defines the angle of internal alignment ζ. Here, the angle ϑ is defined as being between the direction of incident light k and B. The angles α and ϵ are between the direction of the supersonic velocity stream vg and B and between the anisotropy uλ in the radiation field and B, respectively. The vector is the projection of vg on the direction of the magnetic field direction B.

Open with DEXTER

Once, the Rayleigh reduction factor is known, the exact cross-sections for the RT simulation can be calculated with Eqs. (10) and (11) for a given size distribution n(a) between a minimal amin and maximal amax grain size with The same procedure is used for the cross-sections of scattering , absorption , and circular polarization to obtain the matrix required in Eq. (7).

2.2.1. Imperfect internal (II) alignment

The cause of internal alignment (see Purcell 1979) is the Barnett effect (Barnett 1917). Any change in angular momentum (e.g. gas-dust collisions) forces free electrons in the paramagnetic dust grain material to align their spin parallel to the axis of grain rotation. This induces a net magnetic moment and, subsequently, grain alignment. However, this net magnetic moment is not constant because of thermal fluctuations inside the dust grain. Subsequently, thermal fluctuations are transferred into a wobbling motion, leading to a precession of the largest moment of inertia I|| arround J. The distribution function for imperfect internal alignment is (17)(Lazarian 1996), where h = I||/I is the ratio of the momenta of inertia, Jeff is the rms value of the angular momentum, and N is a normalizing constant (see Appendix B for details).

2.2.2. Imperfect Davis-Greenstein (IDG) alignment

For the IDG, the alignment mechanism is paramagnetic relaxation (Davis & Greenstein 1951; Jones & Spitzer 1967; Purcell 1979; Spitzer & McGlynn 1979). In rotating non-spherical paramagnetic dust grains the electrons’ spin axis follows the external magnetic field direction. However, for a sufficiently high angular velocity, the electrons cannot fully respond. Spin-lattice interaction transfers the component of angular velocity perpendicular to the magnetic field into internal heat. In the ISM, paramagnetic relaxation and subsequently perfect alignment of the dust grains is opposed by random gas bombardment.

The ratio of the characteristic timescales of paramagnetic relaxation and gas-dust collisions provides a threshold for the effective grain radius a up to which dust grains align (18)where Tg is the gas temperature and ng is the number density of the gas. Dust grains with radii a above δ0 no longer significantly contribute to polarization. For the IDG, an analytical function for the distribution of the opening angle β (see Spitzer & McGlynn 1979) allows us to calculate (19)(Aannestad & Greenberg 1983), with

(20)and subsequently the Rayleigh reduction factor. To calculate the internal alignment, we use the approximation from Lazarian & Roberge (1997) for the angular momentum (21)where is the thermal angular momentum obtained from gas-dust collisions and s is the aspect ratio of the grains. This approach delivers reliable values in the limit of thermal equilibrium (Td/Tg = 1). However it slightly underestimates the internal alignment for s< 0.5 and Td/Tg< 0.2.

2.2.3. Gold (magneto mechanical) alignment

The alignment of dust grains by purely mechanical alignment was proposed by Gold (1952). This mechanism requires a predominant direction in the velocity field of the gas. The dust grains align with their longer axis parallel to the gas flow. In magneto-mechanical alignment (Lazarian 1994, 1995, 1997), the alignment behavior depends on both gas flow and magnetic field. Here, the reference frame is determined by the Barnett effect and, subsequently, the magnetic field direction, but not the gas stream. In this paper we also refer to magneto-mechanical alignment as GOLD alignment. The gas flow must be supersonic, otherwise the dust grain orientation remains randomized (see Gold 1952; Lazarian 1995).

A distribution function for the cone angle β was derived by Dolginov & Mitrofanov (1976). For spheroidal dust grains, the first order moment is completely defined by the anisotropy in the gas stream (22)where is the part of the relative drift velocity vrel projected on the magnetic field direction and g = h-1 is the grain non-sphericity (Lazarian 1994). In the case of oblate dust grains, g is negative and (23)Here, the Rayleigh reduction factor is independent of wavelength and grain size. The polarization can change sign for an angle α between magnetic field and a predominant gas stream lower than 54° (see Fig. 1). For other grain shapes we refer to Lazarian (1996).

For the internal alignment the angular momentum can be approximated (Lazarian 1997) by (24)with the molecular mass μmH of the gas. For the required conditions of supersonic velocity we locally calculate the Mach-number in each cell with (25)and demand M > 1. Otherwise, we consider the dust grains to be randomized.

Possible mechanical alignment of helical dust grains in sub-sonic environments is currently being discussed (e.g. Lazarian & Hoang 2007) but is not yet in POLARIS.

2.2.4. Radiative torque (RAT) alignment

It was noticed by Dolginov & Mitrofanov (1976) that radiation can both spin up and align dust grains. The alignment mechanism is also known as the Barnett effect.

Based on a numerical approach and the DDSCAT Code (Draine & Goodman 1993; Draine & Flatau 1994, 2000, 2013). Draine & Weingartner (1996, 1997), Weingartner & Draine (2003) demonstrate that RATs can significantly align a number of different grain shapes. The spin up of dust grains is also supported by laboratory experiments (see Abbas et al. 2004). Later, the analytical framework for RAT alignment was developed in Hoang & Lazarian (2007a,b), where they reproduced the numerical findings from Weingartner & Draine (2003).

With increasingly effective dust grain radius a, RATs become more efficient and the dust grains get suprathermally spun up. Analytical models by Hoang & Lazarian (2007a,b) show grain alignment at attractor points with high angular momentum (high-J) and low angular momentum (low-J). The high-J attractor point is stable and, with supra-thermal rotation (ωrad> 3 × ωth), dust grains can be considered to be perfectly aligned. Most dust grains drift to the unstable low-J attractor, where they are prone to thermal fluctuations and are subsequently partially disaligned. To determine the effective radius (aalg) at which the dust grains start to align, we follow the procedure described in Bethell et al. (2007) with (26)Here, QΓ(ϵ) is the alignment efficiency of the RAT, which depends on the angle ϵ between the predominate direction in the radiation field and the magnetic field direction. The alignment efficiency QΓ(ϵ) is also a function of grain size, wavelength, and grain composition, and can be pre-calculated. The shape-dependent parameter δ, as well as the characteristic timescales of gas bombardment tgas and radiation trad acting on the dust grains are introduced in detail in Draine & Weingartner (1996, 1997). Checking for supra-thermal conditions also requires information about the local anisotropy γλ and the mean energy density of the radiation field.

These quantities can be determined with POLARIS by a MC mode, where we store the energy density direction uλ in each cell. In this mode we consider all radiation sources including thermal dust re-emission. Following Lucy (1999) for the mean energy density per path length and wavelength we derive (27)This allows us to calculate the mean energy density (28)and the anisotropy parameter determined by (29)With RATs, gas bombardment can even increase grain alignment by lifting grains from the low-J to the high-J attractor point. At the current stage of its development no analytical function for the distribution of the cone angle β is available. While perfect alignment can be assumed at both high-J and low-J attractor points, internal aligned is just given for suprathermally rotating grains. At low-J the angular momentum can be assumed to be thermal (Hoang & Lazarian 2009) and the distribution of internal alignment (Eq. (17)) becomes a function of grain geometry alone because of h. Now, for the Rayleigh reduction factor, only the moment remains of relevance. When we consider fhigh−J to be the fraction of dust grain aligned at the high-J attractor point the Rayleigh reduction factor is(30)So far, fhigh−J cannot be exactly determined (Hoang & Lazarian 2009) and remains a free parameter.

Once, a dust grain starts to tumble in the presence of a magnetic field a Larmor precession acts on the grains because of the Barnett effect. For the RATs to work, the magnetic field strength has to be sufficiently high to resist disaligment by gas-dust interactions. This requires the characteristic timescale for Larmor precession (tL) to be lower than the gas dumping time (tgas): tL<tgas (Hughes et al. 2009; Hoang & Lazarian 2009), which gives us a limit of (31)for the magnetic field strength. This limitation is similar in nature to that of Eq. (18). Since the coupling with the magnetic field is also the Barnett effect, the same limitation for the field strength is applied to RT simulations with GOLD alignment.

2.3. Scattering on partially aligned dust grains

For each scattering event the coordinate system of the dust grain is determined by the magnetic field direction. The coordinate system of the Stokes vector must be transformed into the dust reference frame. For scattering from one direction (ϑ,ϕ) to a new direction (ϑ′,ϕ′), the Stokes vector S is modified via a 4 × 4 scattering matrix by . Non-spherical, partially aligned dust grains can require a full 16-element matrix and an asymmetry factor g, with both dependent on direction of inclination, scattering direction, wavelength, and grain size, which makes scattering a multi-dimensional problem that is fully supported by POLARIS.

The scattering matrix and asymmetry factor g can be pre-calculated with DDSCAT. For imperfectly aligned wobbling dust grains the incident (ϑ,ϕ) and scattering angles (ϑ′,ϕ′) become functions of alignment angles β and ζ, respectively. To obtain an exact scattering matrix for an average-sized dust grain one has to integrate the alignment distribution functions f(β) and f(ζ), respectively, as well as the grain size distribution n(a) simultaneously. Here, we face the same problem as for the calculation of cross-sections in Sect. 2.2. In the current stage of development, the exact distribution function for RAT alignment is still missing. Solving the exact integrals for GOLD or IDG alignment would push the timescales of 3D MC-RT simulations beyond accessible computational equipment. To handle the scattering on non-spherical imperfectly aligned dust grains, we apply the following approximation:

(32)The same approximation is used for the orientation and size-dependent asymmetry parameter g of the HG-Phase function (see Eq. (4)): (33)where and g(a,λ,ϑ,ϕ)⟩ are the entries of the scattering matrix and the asymmetry factor, respectively, for randomly aligned dust grains. With increasing efficiency in dust grain alignment the Rayleigh reduction factor reaches unity and the dust grains and the contribution of randomly aligned dust grains diminishes.

2.4. Radiative transfer with MHD data

The complexity of star formation can not be accurately mapped by simple analytical models. However, with the advantages of high-resolution MHD simulations, e.g. FLASH Fryxell et al. (2010), it is possible to provide complex 3D physical data as an input for MC-RT simulations. Not only can this data deliver density, magnetic field, and velocity distributions, but it can also provide the exact parameter of stars appearing in such simulations. This stars have to be considered as additional sources of radiation. For RT simulations with MHD data, we have to take the following into account: For the implemented alignment mechanisms to work, the characteristic timescales have to be taken into account. By post-processing MHD data, the shortest time steps in the MHD simulations must exceed the largest timescale of alignment mechanisms to ensure a steady-state of grain alignment (see Draine & Weingartner 1996, 1997; Weingartner & Draine 2003).

Finally, a dust temperature can also be provided by the actual MHD simulations. With temperatures resulting from MHD simulations, we have to deal with two effects of dust heating: RT-heating because of radiation sources and MHD-heating because of the compression and gas-dust interaction (Banerjee et al. 2006; Krumholz et al. 2007). Both temperatures are equally physically well-motivated, however, they are rarely fully implemented in MHD codes. To account for all effects of dust heating, POLARIS optionally performs an extended dust heating algorithm. Originally, the algorithm of Lucy (1999) assumes a cell’s energy content Ė to be empty at the beginning of the radiative dust heating process. Since we have already pre-calculated the emissivity of an ensemble of dust grains with a certain temperature, we can now perform the reverse calculation. Using Eq. (1), we determine the offset amount of energy ΔĖ that corresponds to the dust temperature delivered by the MHD data. We refer to this method as offset dust heating in the following.

2.5. Code optimization

An 3D MC-RT simulation is challenging when it comes to run time and the capabilities of available computational equipment. The POLARIS code is parallelized to run on shared memory machines using the OpenMP library with the propagation of several photon packages simultaneously through the grid. The run-time scales well with the number of used processors. However, complex environments provided by MHD simulations, the physically well-motivated dust models, and alignment mechanisms requires a huge parameter space to solve the RT problem, including polarization in 3D. For this reason POLARIS is equipped with several optimization algorithms to perform RT calculations in a reasonable time. To propagate photon packages in regions with extreme optical depth more efficiently, we implemented the flowing algorithms: The modified random walk (MRW), a technique that uses a diffusion approximation instead of an MC approach (Min et al. 2009; Robitaille 2010) for optically thick regions and the forced first scattering (e.g. Mattila 1970; Wood & Reynolds 1999) that increases the signal-to-noise ratio in optically thin dust environments (see Morris 1984). For the calculation of the RAT alignment, the MRW cannot be applied. While the mean energy density is correctly approached, the MRW leads to a loss of information about the anisotropy of the radiation field. However, since the effective grain radii are pre-calculated and therefore discrete values, a far lower amount of photon packages is required to converge against a distinct value of aalg in the RAT mode.

thumbnail Fig. 2

Dust temperature distribution in the mid-planes of a 3D MC-RT/MHD collapse simulation. Rows show the results of the same model, but of inclinations differing by 90°. In the left column, the dust was heated by 3D MC-RT simulation alone. The middle column shows the dust temperature distribution as a result of MHD shock heating. The dust temperature combined by offset dust heating, as described in Sect. 2.4, is shown in the right column. The length of 1400 AU is for scale. We applied an upper cut-off of 400 K (blue regions) to the mid-plane images for better illustration.

Open with DEXTER

Furthermore, for synthetic intensity and polarization maps, the S/N can be enhanced with the peel-off technique (Yusef-Zadeh et al. 1984) by sending a fraction of the radiation toward the observer at each scattering event. With longer wavelengths, scattering becomes increasingly insignificant and a ray-tracing algorithm can be applied. Here, we simply add up all contributions because of extinction and re-emission along the line of sight. Compared to a radiative transfer with scattering, the ray-tracing approach is more effective and allows us to achieve an excellent S/N (Robitaille 2010; Dullemond 2012).

Radiation sources emit photons in the full wavelength range of its characteristic SED. However, certain wavelengths can be neglected since they do not contribute significantly to dust heating and alignment compared to the total emitted energy of the source. A wavelength range selection algorithm reduces the number of wavelengths considered and subsequently the computational effort.

3. Applications

In this section, we show selected applications of the POLARIS code, demonstrating its features and predictive capability. In particular, we show RT simulations in astrophysical environments, representing different stages of star formation on multiple scales. We run 3D MC-RT simulations to calculate temperature distributions and synthetic polarization maps. A benchmark of POLARIS’ scattering calculations is provided in Appendix A.

The often used dust grain model of Mathis et al. (1977) suggests an upper cut-off in grain radius at amax = 250 nm. More recent models (e.g. Weingartner & Draine 2001), however, indicate an upper limit with μm-radii. Since RAT alignment is heavily dependent on the maximum value of grain size, we pre-calculated a database for oblate dust grains with a fixed aspect, 118 effective radii, and 100 wavelength logarithmically distributed between a ∈ [5 nm:2 μm] and λ ∈ [90 nm:2 mm], respectively. Here, we used the program DDSCAT (Draine & Flatau 2013) and approximated the oblate shape of the dust grains with 171 500 distinct dipoles and an error tolerance of 8.5 × 10-3. As materials, we considered astro-silicate and graphite, used the refraction indices of Weingartner & Draine (2000), and applied the 1/3−2/3 approximation to take care of the anisotropic graphite structure (Draine & Malhotra 1993). We weighted the efficiencies for each wavelength according to the power-law distribution to derive the average cross-sections over all effective radii (see Eq. (15)).

Once the photon package escapes the grid, the degree of linear polarization Pl is determined by the Stokes vector with (34)and its projected position angle χ(35)The aspect ratio of minor axis to major axis remains s = 0.5 in all cases. Higher aspect ratios would result in higher peak values and more contrast in the maps of linear and circular polarization. The dust heating is assumed to be independent of effective radius and aspect ratio, which is in good agreement with the findings of Hong & Greenberg (1980) and Voshchinnikov et al. (1999) within the considered range of parameters.

3.1. Offset dust heating

In this section, we post-processed an MHD collapse simulation with a total mass of 100 M (see Seifried et al. 2011, for a detailed description). We took a snapshot after a simulation time of 5000 yr and consider just the inner most region of the MHD simulation with an accretion disc of about 1400 AU. The dust temperature, as well as position, stellar radii, and surface temperatures of the stars were provided by the simulation. For dust material, we use a mixture of 62.5% astro-silicate and 37.5% graphite with an an upper cut-off for the size distribution at amax = 250 nm. The highest temperature in these simulations is below the sublimation temperature of the silicate and graphite grains, respectively. First, we ignore the MHD dust temperature and performed a 3D MC-RT simulation, with the stars as radiation sources alone. Then, we adjust the temperature with the method described in Sect. 2.4.

In Fig. 2, we show the resulting dust temperature distributions in the center x-y planes and x-z planes. For the RT simulations with no offset energy (Fig. 2 left panels), the dust temperature shows a diffuse distribution with the highest temperature near the stars and regions with lowest densities. The dust temperature distribution matches the density structure near the disk region, while the outflow regions remain hidden. Regions with high densities are shielded from stellar radiation and cannot be heated sufficiently. High density regions, however, gain temperature by the dynamical processes of shock heating in MHD simulations (Fig. 2 middle panels). The MHD dust temperature resembles the dust density in all regions while underestimating the temperature towards the center region of the accretion disk. By considering both RT and MHD leads to a shift in the profitability distribution of the wavelength of re-emission (see Eq. (6)) towards shorter wavelengths. In this regime of wavelengths, absorption of radiation is more effective because of higher cross-sections of absorption Cabs. This results in a net temperature that is higher than single temperatures alone, which would not have been the case by simply adding up the dust temperatures resulting from RT and MHD simulations.

3.2. Disk models

In this section, we trace the magnetic field morphology in circumstellar disks of different mass with POLARIS. The magnetic field plays an important role in the formation of stars, as well as in the formation and evolution of the circumstellar disk and, subsequently, in the formation of planets. Angular momentum of the gas component of the disk can, effectively, be transferred outwards by magnetic braking and interstellar outflows allow the fall of gas towards the center. Thus, to understand circumstellar disk and planet formation, we need to investigate the properties of the underlying magnetic field that is involved in the formation process.

By considering the dust grains to be partially aligned with the direction of the magnetic field, the field morphology may be investigated by polarization observations in the sub-mm regime of wavelength without the contribution of scattering in the polarization signal.

3.2.1. Model setup

The density distribution of our disk models are parametrized as presented in Shakura & Sunyaev (1973). The model assumes the density distribution to follows: (36)where (37)Here, r and z are cylindrical coordinates, h0 is the scale height at r0 and the quantities α and β are parameters that characterize the radial density profile and the disk flaring. We adjusted n0 for the entire disk mass to range between 10-3 [M] and 10-6 [M] for the different models. Theory and observations show a strong correlation with nBκ between magnetic field strength and density over several orders of magnitude (Seifried et al. 2011; Mestel 1966; Crutcher 1999). Here, we assume κ = 0.6 to keep the magnetic field strength between 10-8 T in the center and 10-17 T at the edges of the model space. The magnetic field is assumed to be dragged with the rotation of the disk because of ionized gas. Therefore, the field morphology can be approximated by a toroidal field (see Reissl et al. 2014).

Here, the considered dust consists of only astro-silicate as material with a cut-off in grain size at amin = 5 nm and amax = 2 μm, respectively with a power-law size distribution of dn(a) ∝ a-3.5da.

For the dust heating, we consider an offset dust temperature of 10 K within the boundaries of the disk and post-process the temperature with POLARIS and a central star as a radiation source for all disk models. The offset dust temperature is of minor influence to the following calculations but ensures a non-zero temperature in each cell in the denser regions of the model space, which enables us to reduce the number of photon packages. The model of Shakura & Sunyaev (1973) assumes the dust to be in thermal equilibrium with the gas (Tg = Td). In this case, Eq. (20) reaches unity and the IDG Rayleigh reduction factor (see Eq. (19)) is undefined. Subsequently, the IDG mechanism cannot be applied in the disk models. Hence, in this section, we only consider the effects of grain alignment according to RAT theory and assume the internal alignment to be perfect. As a second step, we apply the RAT mode of POLARIS to calculate the anisotropy and energy density of the radiation field in the disk, as described in Sect. 2.2.4, to determine the minimal dust grain radius aalg at which the dust grains start to align.

thumbnail Fig. 3

Distribution of grain sizes aalg, where dust particle start to align to the magnetic field direction B according to RAT theory. The plane shown is perpendicular to the mid-plane of the disks. The disk model D03 with 10-3 M is in the left panel and the model D05 with 10-5 M is in the right panel. Lower values of aalg result in higher degrees of linear polarization.

Open with DEXTER

Table 1

Physical quantities for the applied disk models.

thumbnail Fig. 4

Maps with pattern of linear polarization overlaid with normalized orientation vectors at an exemplary wavelength of λ = 515 μm. We applied an offset of 90° to the orientation vectors to match the projected toroidal magnetic field morphology. The model D03 with 10-3 M is in the right column and model D05 with 10-5 M is shown in the left column for inclination angles from i = 30° (top row) to i = 90° (bottom row) in steps of 15°.

Open with DEXTER

3.2.2. Results

Figure 3 shows the resulting distribution of the D03 and D05 disk models. The distribution of aligned radii is layered with smaller aligned dust grains, which are starting to align at the disk’s surface with larger grains towards the center plane. The model with the higher disk mass D03 on the left hand side of Fig. 3 shows a distinct distribution of imperfectly aligned dust grains. In contrast to model D03 for model D05 on the right hand side of Fig. 3, smaller dust grains are better aligned at the surface of the disk. However, with decreasing mass, the alignment of dust grains increases quickly towards the center plane and the outer regions. This is because of different dust and stellar contributions to the radiation field as well as a local dust temperature.

The optical depth decreases with lower disk mass. Therefore, the surface is more efficiently heated in the disks with higher mass. Whereas the grain alignment is suppressed near the sublimation radius for all models as a result of high temperatures in this region, the stellar radiation becomes more irrelevant in the outer regions of the disk. Consequently, the D03 model is less influenced by the central star. The radiation in the innermost regions in each disk originates from thermal photons emitted in the hot surface layers and not the star itself. Subsequently, the higher mass model D03 has aligned dust grains in much deeper layers as the lower mass model D6. The resulting layered distribution in dust grain alignment is consistent with the results of Hoang & Lazarian (2014).

Table 2

Physical parameter for the applied cloud models.

In Fig. 4, we show the effects of thermal dust re-emission, which is dependent on the total disk mass. We compare the resulting polarization maps for the disk models D03 and D05 as a function of inclination angle at an exemplary wavelength of λ = 515 μm. The pattern and orientation of linear polarization in Fig. 4 are in good agreement with the predictions by Cho & Lazarian (2005, 2007). Since all disk models are optically thin at that wavelength, the synthetic polarization maps appear to be symmetric.

Because lower values of the characteristic size of aligned grains aalg are associated with a higher degree of linear polarization (see Eq. (30)), the D03 disk model shows a more distinct pattern of linear polarization. A lower degree of linear polarization in low mass disks is consistent with the results shown in Fig. 3. This would also hold even if we considered imperfect internal alignment since the internal distribution function for RATs is independent of the local physical conditions of the disk models and just a function of grain geometry (see Sect. 30). The pattern in the normalized polarization vectors are identical for all models, which reveals the projected toroidal magnetic field geometry at all inclination angles since the influence of scattering can be neglected at λ = 515 μm. However, just for the D03 model, the degree of linear polarization (Pl> 0.5%) should be detectable in the entire disk, even for low inclination angles.

The patterns of linear polarization have sharp cut-offs with no polarization at the borders of the model. This is a result of the applied power-law correlation between density and magnetic field strength (nBκ). Hence, the magnetic field no longer satisfies the minimal required field strength (see Eq. (31)) in the thinner disk regions.

3.2.3. Discussion

Aligned dust grains as an explanation of polarized radiation in the far-infrared and submm from disks has already been indicated by the observations of Tamura et al. (1999). However, the extent to which the observation of polarized light reveals the underlying magnetic field morphology is still a matter for debate. Numerical calculations (Cho & Lazarian 2007; Hoang & Lazarian 2014) predict that polarized radiation should emerge from the circumstantial disk region because of RAT alignment. In contrast to predictions, observations of circumstellar systems (Hughes et al. 2009) show no detectable degree of linear polarization at all. The overestimation of RAT alignment, because of perfect internal alignment in earlier studies (Cho & Lazarian 2007), cannot account for this finding alone. An internal alignment with a conservative chosen parameter fhigh−J = 0.25 would still result in a detectable degree of linear polarization. However, fhigh−J can still not be analytically determined and might also vary depending on the predominant local conditions of each disk. Another explanation might be the shape and composition of the dust grains. We cannot expect a similar amount of non-spherical dust grains that can align to be equally present in all circumstantial disks since this may vary significantly from disk to disk.

The distribution of dust grain sizes is also not expected to be equal in the entire disk. Dust grain growth and dust settling are processes that redistribute dust grains of different sizes in the disk (see e.g. Tanaka et al. 2005). As a result, larger dust grains are more likely near the disk plane. Larger dust is also more likely to be aligned by RATs than smaller dust grains on the disk surface. Dust settling in disks might also help to explain the low amount of linear polarization and need to be considered in future studies with a more sophisticated dust model.

3.3. Cloud model

thumbnail Fig. 5

Schematic illustration of the cloud models C01 (left panel) and C02 (middle panel). Both models differ only by the position of a single star. The isocontour surface at log (nd) = 3.27 m-3 is in dark yellow. The direction of the common gas stream v (left panel) toward the center is shown in green lines and the analytical hourglass field morphology is shown in grey. Red lines and blue bars indicate the different angles between the radiation field and the magnetic field direction. The letters A–D indicate areas with expected characteristic features for RAT and GOLD alignment, respectively.

Open with DEXTER

thumbnail Fig. 6

Ratio of polarization cross-sections of imperfectly aligned dust grains Cx, to perfectly aligned dust grain Cp in the mid-plane of the model space for the C01 model (top row) and C02 model (bottom row). The x stands for IDG alignment (left column), RAT alignment (middle column), and GOLD alignment (right column). In the middle columns, dotted white lines indicate the regions where the alignment efficiency QΓ(ϵ) is at its minimum. The dotted red lines in the right column show the transition where GOLD alignment changes the sign of linear polarization. The areas AD correspond to that of Fig. 5.

Open with DEXTER

In this section, we examine the pattern of linear polarization emerging from different grain alignment theories separately. The main focus here is on the ambiguities associated with different grain alignment theories used to infer the underlying magnetic field morphology on a predefined cloud model.

Our cloud model consists of two spheres following a Bonnor-Ebert density profile (Ebert 1955; Bonnor 1956) with a distance of 0.375 ly from each other. The characteristic radii of the density distribution are 0.1 ly and 0.2 ly, with a number density for the gas of 1012 m-3 and 1013 m-3, respectively. Density and temperature are chosen to be consistent with theoretical models and observations (Launhardt & Henning 1997; Sipilä et al. 2011).

3.3.1. Model setup

We modeled an hourglass-shaped magnetic field geometry by using an analytical function (see Reissl et al. 2014), with its center in the denser sphere as expected by observations (e.g., Girart et al. 1999; Frau et al. 2011). The field strength was assumed to be constant with |B| = 1.5 × 10-8T (Falgarone et al. 2008; Beuther et al. 2010). To compare calculations performed by POLARIS with the expected predictions for the GOLD alignment mechanism, we also assume a supersonic velocity stream v with an in-fall direction to the common center of both Bonnor-Ebert spheres. In this section, we also consider the effects of the internal alignment with a correlation factor of fc = 0.6 and assume fhigh−J = 0.5. The initial dust temperature is adjusted with a single star as radiation source (see Table 2). We consider two different configurations concerning the position of the star. In model C01, the star is embedded in the center of the denser Bonnor-Ebert sphere, and in model C02, the star is 0.53 ly away from the common center of both spheres. A schematic illustration of the cloud models provided in Fig. 5.

For the dust model, we apply a range of effective radii of a ∈ [5 nm:250 nm] with a size distribution of nda-3.5 and consider a mixture of 62.5% silicate and 37.5% graphite (Mathis et al. 1977; Draine & Li 2001). However, just the silicate grains align with the direction of the magnetic field while the graphite remains randomized (Mathis et al. 1977). An initial sphere temperature of Td = 10 K was used for both cloud models, C01 and C02, and we corrected the temperature according to the method described in Sect. 2.4. Here, we assume dust and gas not to be in thermal equilibrium and assume a correaltion of Tg = 10 × Td. In a final step, we applied the RAT mode of POLARIS to calculate the distribution of the size cutoff aalg of dust grains to align. With this configuration the models allow the calculation of the characteristic linear polarization pattern for IDG, RAT, and GOLD alignment mechanism.

3.3.2. Results

Figure 6 shows the ratio of cross-sections of partially aligned dust grains to the cross-sections of perfectly aligned dust grain in the mid-plane of both models separately for all alignment mechanisms at an exemplary wavelength of λ = 723 μm. In both cloud models, the surrounding regions of the stars are also the regions of highest dust temperature. The IDG mechanism reacts sensitively to both values, the density, and the temperature, because of the threshold δ0 (see Eq. (18)). Therefore the efficiency of grain alignment is lowest towards the center of the clouds and near the star (see Fig. 6 left panels).

The expected alignment pattern according to RAT theory is contrary to that of the IDG alignment. The emission of heated dust and stellar radiation contribute to the local energy density uλ and the anisotropy parameter can reach a maximum value up to γλ = 0.9 in close proximity to the star. The characteristic dust grain radius aalg decreases and dust grains with radii above aalg become aligned. Since, the process of determining RAT alignment quantities requires two separate MC simulations (dust heating and calculation of the radiation field), the grain alignment of RATs has an inherent higher amount of noise. This becomes apparent in Fig. 6 (middle panel) in the inner cloud regions. Besides the amount of energy density and anisotropy, the grain alignment also highly depends on the RAT alignment efficiency QΓ(ϵ). This parameter, in turn, depends on the chosen dust grain geometry and has, in our dust model of oblate dust grains, its minimum at angles of ϵ = 0° and ϵ = 90° between radiation and magnetic field direction. This leads to areas of reduced grain alignment, indicated by white dotted lines in Fig. 6 (middle panel). In Fig. 6 (middle top panel), the low degree of grain alignment in the areas A and C is a direct result of the angle dependency of RAT theory because of QΓ(ϵ). However, the low grain alignment on the left half of the panel is a result of the shielding of this region by the dense clouds from stellar radiation and thus the overall energy density remains low. The expected low grain alignment in the area D cannot be found. This is because the overall energy density outweighs the angle dependency in this very thin region. In Fig. 6 (middle bottom panel), the low degree of grain alignment in the areas AD follows the angle between predominant direction of the radiation field and the magnetic field direction more clearly. This result is in excellent agreement with the expectations of the model setup shown in Fig. 5 and RAT theory.

Since we assume a constant supersonic velocity field, GOLD alignment alone barely depends on temperature and density. However, a dependency exists because of internal alignment (see Eq. (24)) and the Mach-limit (see Eq. (25)). As with RAT alignment, GOLD alignment is also limited by the local magnetic field strength (see Eq. (31)). Since the internal alignment mechanism is sensitive to high temperatures, this manifests itself as a sphere of completely randomized dust grains in the surrounding regions B (Fig. 6 top right panel) and A (Fig. 6 bottom right panel) around the star. By comparing Fig. 6 (left panels), Fig. 7 (left panels), and Eq. (23), we can see that the dominant parameter maintains the angle α between gas stream v and magnetic field direction B. The areas of minimal grain alignment appear along the transition where GOLD alignment changes the sign of the Rayleigh reduction factor, which is shown by red lines in Fig. 6 (middle panel).

The resulting maps of linear polarization for the cloud models C01 and C02 are presented in Fig. 7. The normalized orientation vectors of linear polarization have an offset angle of 90° to be parallel with the projected magnetic field direction. For the IDG mechanism (Fig. 7 top row), the center of depolarization moves with the position of the star because of the change in the dust temperature distribution, leading to a nearly rotationally symmetrical pattern in both cases.

The polarization pattern for RAT alignment in the middle column of Fig. 6 follows the regions of highest radiation. Therefore, a high polarization coincides with the position of the star. The degree of linear polarization is reduced in the C01 model along the regions of minimal grain alignment efficiency QΓ(ϵ). The same effect can be detected in the polarization map of the C02 model. However, here this effect is less evident leading to a butterfly-like pattern of linear polarization in the surrounding area of the star. Both, IDG and RAT alignment mechanisms resemble the underlying hourglass field morphology. The GOLD alignment in turn, can change the alignment angle along the line of sight and this, subsequently, changes the polarization direction.

thumbnail Fig. 7

Maps with pattern of linear polarization overlaid with normalized orientation vectors at an exemplary wavelength of λ = 723 μm. We applied an offset of 90° to the orientation vectors to match the projected hourglass magnetic field morphology. The cloud model (C01) is in the left column and the model (C02) is in the right column for dust grains aligned with the IDG mechanism (top row), RAT mechanism (middle row), and GOLD mechanism (bottom row).

Open with DEXTER

thumbnail Fig. 8

Degree of linear polarization of the C01 model (top row) and C02 model (bottom row) as a function of normalized intensity (P-I relation). The results with IDG alignment (right column), RAT alignment (middle column), and GOLD alignment (left column). Both models C01 and C02 (blue dots) were rotated by 90° around the y-axis of Fig. 5 (red dots).

Open with DEXTER

A correlation between intensity I and polarization Pl was recognized in Henning et al. (2001) following the power law P(I/Imax)α where α is a fit parameter. This finding is consistent in the literature of observations (see Matthews & Wilson 2002; Lai et al. 2003; Cho & Lazarian 2005). Since the PI relation is hyperbolic in nature, we expect to find a flat tail when I/Imax reaches unity, making it an ideal test for different grain alignment theories. Therefore, we calculated the increasing depolarization towards lower intensities for our cloud models C01 and C02, making use of the polarization maps of Fig. 7. Here, each data point presents a different line of sight. We also calculated the PI relations for our cloud models rotated by 90°. The PI relation of the IDG mechanism shown in Fig. 8 (left column) is consistent with observations. Since the anisotropy in radiation and grain alignment is directly correlated and stellar radiation is the dominant source of radiation in both cloud models, the PI relation in the RAT case is less definitive. In contrast to the IDG alignment, the PI relation of RAT alignment starts with an increase in polarization up to a point of I/Imax ≈ 0.1 and then remains almost constant at its maximum value. However, the RAT alignment mechanism matches the expected PI relation curvature in agreement with observations for larger fractions of I/Imax.

The GOLD alignment shows an exceptional behavior. In contrast to observational findings, the degree of linear polarization remains constant in the non-rotated case and increases even slightly for the rotated case. Data points show a broader scatter, compared with other alignment mechanisms.

The characteristic appearance of each grain alignment theory in the maps of linear polarization is often discussed to account for observational data (see Andersson 2015, for review). So far, RAT alignment is considered among the most promising theories of grain alignment. The observational fact that grains are better aligned in close proximity to a bright star is a strong indication of RAT alignment (Matsumura et al. 2011). A direct correlation between grain alignment and direction, as well as the strength of the radiation field, was also confirmed in Andersson & Potter (2010) and Andersson et al. (2011) next to the star HD 97300 in the Chamaeleon I star-forming region. These observational findings are in accordance with the results of the synthetic linear polarization maps that were calculated with RAT alignment, as presented in this section.

3.3.3. Discussion

The velocity-dependent flip in the orientation of polarization vectors is characteristic of GOLD alignment and may help to determine areas with supersonic gas streams by observations. From the observational point of view, a flip in polarization angles was reported along the direction of molecular outflows (e.g. Rao et al. 1998; Cortes et al. 2006). However, the rather unusual pattern of linear polarization, as presented in this section, can so far not be confirmed by observations. This may be a result of the limited set of parameters of the cloud models that were considered. The two main limiting factors are velocity and magnetic field strength (see Eqs. (23) and (31)), which were assumed to be constant for simplicity. Recent progress on the field of mechanical grain alignment also indicates that dust-grain alignment may also be efficient in subsonic environments (Lazarian & Hoang 2007).

The IDG alignment has still its place in accounting for the observed alignment of small dust grains Clayton et al. (1992), Hoang et al. (2014). The maximum values in linear polarization are in good agreement with observations (Gonçalves et al. 2005; Davidson et al. 2011) for all the applied alignment mechanisms. Taking the maximal degree of linear polarization as a measurement of the prevailing alignment mechanism, the polarization is clearly dominated by IDG in the outer parts followed by RAT near the star. Higher linear polarization degrees can easily be achieved for RATs by increasing the free parameter fhigh−J and the upper cut-off amax in the grain size distribution of the applied dust models. However, exact values of these parameter are currently speculative. In contrast to the IDG and RAT alignment, models considering GOLD alignment hold the lowest degree of linear polarization, which indicates that the observable net polarization of the cloud would lack any flip in polarization angles. This is also consistent with the behavior of the calculated PI relations. Both RAT and GOLD alignment show an unexpected distribution of linear polarization. This uncertainty and the deviations in calculated PI relations may be resolved with more sophisticated models from dedicated MHD simulations in forthcoming parameter studies.

4. Summary and outlook

We have presented POLARIS (POLArized RadIation Simulator), a newly developed 3D MC-RT radiative transfer code that was designed to perform self-consistent dust heating calculations and polarization simulations. The main novelty of the POLARIS code is the calculation of linear and circular polarization maps considering dichroic extinction, thermal re-emission, birefringence, and scattering on dust grains that are partially aligned with the direction of the magnetic field. In contrast to previous MC codes, POLARIS calculates the impact of imperfect grain alignment on light polarization by considering the major classes of state-of-the-art dust grain alignment theories (alignment owing to paramagnetic relaxation (IDG), alignment because of radiation-dust interaction (RAT), and mechanical alignment (GOLD) as a result of gas streams).

To calculate the dust temperature distribution, we use the combined techniques of continuous absorption and immediate re-emission. The code is openMP parallelized and RT simulations are performed in an adaptive octree cartesian grid, which allows for high performance and the consideration of MHD simulations as input data. Further advanced optimization techniques (ray-tracing, wavelength range selection, forced first scattering, peel-off technique, and modified random walk) enable us to perform complex 3D MC-RT simulations, including polarization on dust grains in a reasonable time. The code considers various sources of radiation to cover a broad variety of astrophysical problems.

The predictive capability of the POLARIS code was tested to run in specific complex environments resulting from MHD simulations, as well as in case studies with analytical density, temperature, velocity, and magnetic field distributions. Here, we applied a method for combining the dust temperatures. We tested the dust heating algorithm of POLARIS in a complex MHD collapse simulation. The resulting distribution of combined temperatures matches simultaneously the compressed dust distribution and the positions of the considered sources of radiation.

The capability to create polarization maps according to the RAT alignment mechanism was tested in circumstellar disk models that were consistent with theoretical predictions. Here, we modeled the disk with an analytical function and applied a toroidal magnetic field geometry. The disk was heated by a central star. Stellar radiation and thermal dust re-emission lead to an anisotropy in the radiation field. This resulted in characteristic layers of dust grain sizes in the disk with larger sizes towards the mid-plane. This is in agreement with RAT alignment theory and verified the predictive capability of of the POLARIS code.

Linear polarization greatly depends on the dominant alignment mechanism. Each mechanism has an effect, in a unique way, on the local physical parameter leading to a characteristic pattern of linear polarization. We tested this scenario in an analytical model of two Bonnor- Ebert spheres with their centers close to each other, a constant velocity field, and an hourglass-shaped magnetic field morphology. We heated the dust with a star as the radiation source to calculate the anisotropy in the radiation field for RAT alignment. Finally, we derived synthetic maps of linear polarization, which considered the IDG, GOLD, and RAT grain alignment mechanism. The resulting pattern of linear polarization resembled the physical parameter in accordance with applied grain alignment theories and showed the expected PI-relations known from observations.

In forthcoming papers, we will present the code extensions for state-of-the-art radiative line transfer, including Zeeman effect. Furthermore, we will apply POLARIS to interpret observational data by calculating the characteristic polarization pattern in MHD outflow simulations. Here, we follow questions about the potential of multi-wavelength polarization measurements to identify distinct sub-structures of a larger magnetic field morphology.

The code is continuously in development to keep up with scientific progress, especially in grain alignment theory. POLARIS is available only on request. However, sample files, scripts and documentation are publicly available from http://www1.astrophysik.uni-kiel.de/~polaris/

Acknowledgments

S.R. and S.W. thank the anonymous referee for providing useful feedback that improved this paper. We also wish to thank colleagues Daniel Seifried for providing MHD data, Florian Kirchschlager for assistance with DDSCAT, Jan Philipp Ruge and Gesa H.-M. Bertrang for fruitful discussions about disk models and grain alignment, Florian Ober and Peter Scicluna for their contributions to radiative transfer and dust heating algorithms. We also wish to thank the students Yong Kyung OH and Jan Thiel for serving as code testers. For this project the authors S.R. and S.W. acknowledge the support of the DFG: WO 857/11-1.

References

Appendix A: Benchmark

thumbnail Fig. A.1

Resulting dust temperature Td for the D03 disk model (see Table 1). The outer left panel shows the temperature distribution in a plane perpendicular to the mid-plane of the disk. The plots show the temperature along the edge of the disk, perpendicular to the mid-plane (left plot), along the mid-plane of the disk (middle plot), and the distribution along the disks surface (right plot). The red lines show the results from the MC3D code, the blue ones the results from the POLARIS code, and the gray lines show the errors between both codes.

Open with DEXTER

thumbnail Fig. A.2

Resulting distribution of the polarized flux Fp for the D03 disk model (see Table 1). The outer left panel shows the polarized flux distribution as a result of scattered stellar radiation under an inclination angle of i = 45° and a wavelength of λ = 730 nm. The plots show the polarized flux through the center of the disk along the horizontal direction (left plot), along the vertical direction (middle plot), and the bisecting line (right plot). The red lines show the results from the MC3D code, the blue ones the results from the POLARIS code, and the gray lines show the errors between both codes.

Open with DEXTER

thumbnail Fig. A.3

The same as Fig. A.1 for the D06 disk model.

Open with DEXTER

thumbnail Fig. A.4

The same as Fig. A.2 for the D06 disk model.

Open with DEXTER

To make predictions about the impact of imperfectly aligned dust grains on polarization measurements, the POLARIS code goes beyond previous approaches in this field. Its treatment concerning polarization as a result of dichroic extinction and thermal re-emission, combined with state-of-the-art dust grain alignment mechanisms is unique in that way. Therefore, essential features cannot be tested because of the absence of standardized benchmark tests. However, dust heating and polarization calculations of scattered light are implemented in many codes, but the POLARIS code is limited, so far, by its implemented octree grid geometry, making it inapplicable to benchmarks such as Pinte et al. (2009). Although, POLARIS is a completely new line of development, we used the established and well-tested MC code MC3D (see Wolf 2003) as a reference for light scattering and dust heating. Here, we confirm the accuracy of POLARIS by comparing the results of the disk models D03 – D06, presented in Sect. 3.2, with the results of MC3D. To provide identical test cases for both codes, we converted the standard MRN MC3D dust catalog with spherical dust grains in a POLARIS database file format. We processed the MC3D output of the density distribution to create an octree grid that is required by POLARIS. Since, the dust temperature in low density areas amounts to unrealistically high values, and MC3D does not control the sublimation temperature of the applied dust materials, we considered cells with a density of nd < 10-20 m-3 to be empty. For the MC calculations of the dust temperature distribution Td and the polarized flux Fp = P × Fλ, because of scattering, we used the parameter of Table 1. The deviation between the results of both codes was quantified with an error of e = 100 | x1x2 | /max(|x1 | , | x2|), where x stands for dust temperature and polarized flux, respectively.

The dust temperature distribution of the disk model C03 resulting from the MC3D code and POLARIS code is shown in Fig. A.1. In the outer regions, the temperature distributions match well, the difference in temperature is much higher at the inner edge of the disk. This difference is a result of the available grid geometries implemented in both codes. While MC3D performed the RT calculations on a spherical grid with smaller cells towards the center, POLARIS uses a Cartesian octree grid with a constant minimal cell size. The geometrical difference in both coordinate systems also results in a different shaped inner border of the disk. Because of this, POLARIS overestimates the temperature in the inner regions. For the same reason the dust temperature differs in the C06 disk model shown in Fig. A.3.

Finally, we compare the resulting polarized flux in Fig. A.2 for the C03 model and in Fig. A.4 for the C06 at a distance of 140 pc to the observer. While the overall trend along the selected lines matches quite well for both codes, the region of deviation remains near the inner radii in all disk models. The higher resolution in the center region of MC3D’s grid leads to a layer of higher dust density at the inner edge. Subsequently, more light is scattered towards the observer. In comparison to MC3D, the synthetic images of the polarized flux POLARIS have a higher S/N because of the implemented optimizations of the forced first scattering algorithm and the peel-of

technique (see Sect. 2.5). However, in all test cases, the results of POLARIS fit those of MC3D within the inherent limitations of applied grid geometries.

Appendix B: Internal alignment equation

Numerical integration is one of the challenges for performing MC-RT simulations in reasonable time. Here, we present the approximation we applied for the calculation of internal alignment. While the analytical solution for IDG, Gold, and RAT alignment are available (see Sect. 2.2), the internal alignment still requires numerical integration. However, the distribution function can be rewritten as (B.1)with α = J2/ 2I||kBTd and δ = h−1. The degree of alignment is determined by (B.2)and also has an exact solution. This enables us to calculate the internal alignment with x = α × δ(B.3)Here, erfi is the complex error function, which can easily be pre-calculated and interpolated. Obviously, when internal alignment is taken into account. This is consistent with physics because in the supra-thermal regime () and for disk-like dust grains (δ → 1), internal alignment becomes irrelevant (see Lazarian 1996; Lazarian & Roberge 1997). For x → 0, the first order moment remains positive and polarization vectors also cannot change their sign in IDG or RAT alignment as a result of imperfect internal alignment.

Appendix C: Combined Rayleigh reduction factor

As shown in this paper, several mechanism for grain alignment have been proposed over the decades. However, these mechanisms are not mutually exclusive and most are probably simultaneously at work inside the ISM. To investigate this scenario, POLARIS enabled us to combine several grain alignment mechanisms in a single 3D MC-RT simulation. Following Lazarian (1995), we use (C.1)for the combined Rayleigh reduction factor. Here, we need to emphasize that this approach remains a first approximation. It requires further research to investigate its validity and accuracy.

All Tables

Table 1

Physical quantities for the applied disk models.

Table 2

Physical parameter for the applied cloud models.

All Figures

thumbnail Fig. 1

Geometrical configuration of a dust grain partially aligned with its angular momentum J to the magnetic field direction B. The precession of J around B defines the cone angle β, and the precession of maximum moment of inertia I around J defines the angle of internal alignment ζ. Here, the angle ϑ is defined as being between the direction of incident light k and B. The angles α and ϵ are between the direction of the supersonic velocity stream vg and B and between the anisotropy uλ in the radiation field and B, respectively. The vector is the projection of vg on the direction of the magnetic field direction B.

Open with DEXTER
In the text
thumbnail Fig. 2

Dust temperature distribution in the mid-planes of a 3D MC-RT/MHD collapse simulation. Rows show the results of the same model, but of inclinations differing by 90°. In the left column, the dust was heated by 3D MC-RT simulation alone. The middle column shows the dust temperature distribution as a result of MHD shock heating. The dust temperature combined by offset dust heating, as described in Sect. 2.4, is shown in the right column. The length of 1400 AU is for scale. We applied an upper cut-off of 400 K (blue regions) to the mid-plane images for better illustration.

Open with DEXTER
In the text
thumbnail Fig. 3

Distribution of grain sizes aalg, where dust particle start to align to the magnetic field direction B according to RAT theory. The plane shown is perpendicular to the mid-plane of the disks. The disk model D03 with 10-3 M is in the left panel and the model D05 with 10-5 M is in the right panel. Lower values of aalg result in higher degrees of linear polarization.

Open with DEXTER
In the text
thumbnail Fig. 4

Maps with pattern of linear polarization overlaid with normalized orientation vectors at an exemplary wavelength of λ = 515 μm. We applied an offset of 90° to the orientation vectors to match the projected toroidal magnetic field morphology. The model D03 with 10-3 M is in the right column and model D05 with 10-5 M is shown in the left column for inclination angles from i = 30° (top row) to i = 90° (bottom row) in steps of 15°.

Open with DEXTER
In the text
thumbnail Fig. 5

Schematic illustration of the cloud models C01 (left panel) and C02 (middle panel). Both models differ only by the position of a single star. The isocontour surface at log (nd) = 3.27 m-3 is in dark yellow. The direction of the common gas stream v (left panel) toward the center is shown in green lines and the analytical hourglass field morphology is shown in grey. Red lines and blue bars indicate the different angles between the radiation field and the magnetic field direction. The letters A–D indicate areas with expected characteristic features for RAT and GOLD alignment, respectively.

Open with DEXTER
In the text
thumbnail Fig. 6

Ratio of polarization cross-sections of imperfectly aligned dust grains Cx, to perfectly aligned dust grain Cp in the mid-plane of the model space for the C01 model (top row) and C02 model (bottom row). The x stands for IDG alignment (left column), RAT alignment (middle column), and GOLD alignment (right column). In the middle columns, dotted white lines indicate the regions where the alignment efficiency QΓ(ϵ) is at its minimum. The dotted red lines in the right column show the transition where GOLD alignment changes the sign of linear polarization. The areas AD correspond to that of Fig. 5.

Open with DEXTER
In the text
thumbnail Fig. 7

Maps with pattern of linear polarization overlaid with normalized orientation vectors at an exemplary wavelength of λ = 723 μm. We applied an offset of 90° to the orientation vectors to match the projected hourglass magnetic field morphology. The cloud model (C01) is in the left column and the model (C02) is in the right column for dust grains aligned with the IDG mechanism (top row), RAT mechanism (middle row), and GOLD mechanism (bottom row).

Open with DEXTER
In the text
thumbnail Fig. 8

Degree of linear polarization of the C01 model (top row) and C02 model (bottom row) as a function of normalized intensity (P-I relation). The results with IDG alignment (right column), RAT alignment (middle column), and GOLD alignment (left column). Both models C01 and C02 (blue dots) were rotated by 90° around the y-axis of Fig. 5 (red dots).

Open with DEXTER
In the text
thumbnail Fig. A.1

Resulting dust temperature Td for the D03 disk model (see Table 1). The outer left panel shows the temperature distribution in a plane perpendicular to the mid-plane of the disk. The plots show the temperature along the edge of the disk, perpendicular to the mid-plane (left plot), along the mid-plane of the disk (middle plot), and the distribution along the disks surface (right plot). The red lines show the results from the MC3D code, the blue ones the results from the POLARIS code, and the gray lines show the errors between both codes.

Open with DEXTER
In the text
thumbnail Fig. A.2

Resulting distribution of the polarized flux Fp for the D03 disk model (see Table 1). The outer left panel shows the polarized flux distribution as a result of scattered stellar radiation under an inclination angle of i = 45° and a wavelength of λ = 730 nm. The plots show the polarized flux through the center of the disk along the horizontal direction (left plot), along the vertical direction (middle plot), and the bisecting line (right plot). The red lines show the results from the MC3D code, the blue ones the results from the POLARIS code, and the gray lines show the errors between both codes.

Open with DEXTER
In the text
thumbnail Fig. A.3

The same as Fig. A.1 for the D06 disk model.

Open with DEXTER
In the text
thumbnail Fig. A.4

The same as Fig. A.2 for the D06 disk model.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.