{'fraction_non_alphanumeric': 0.031677173709506276, 'fraction_numerical': 0.016088259159852694, 'mean_word_length': 4.777587450414713, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Radiative transfer with POLARIS: I. Analysis of magnetic fields through synthetic dust continuum polarization measurements

Aims: We present POLARIS (POLArized RadIation Simulator), a newly developed three-dimensional Monte-Carlo radiative transfer code. POLARIS was designed to calculate dust temperature, polarization maps, and spectral energy distributions. It is optimized to handle data that results from sophisticated magneto-hydrodynamic simulations. The main purpose of the code is to prepare and analyze multi-wavelength continuum polarization measurements in the context of magnetic field studies in the interstellar medium. An exemplary application is the investigation of the role of magnetic fields in star formation processes. Methods: We combine currently discussed state-of-the-art grain alignment theories with existing dust heating and polarization algorithms. We test the POLARIS code on multiple scales in complex astrophysical systems that are associated with different stages of star formation. POLARIS uses the full spectrum of dust polarization mechanisms to trace the underlying magnetic field morphology. Results: Resulting temperature distributions are consistent with the density and position of radiation sources resulting from magneto-hydrodynamic (MHD) - collapse simulations. The calculated layers of aligned dust grains in the considered cirumstellar disk models are in excellent agreement with theoretical predictions. Finally, we compute unique patterns in synthetic multi-wavelength polarization maps that are dependent on applied dust-model and grain-alignment theory in analytical cloud models.


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. Selfgravitating 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, astrochemistry, 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 nearinfrared 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 Article number, page 1 of 17 arXiv:1604.05305v2 [astro-ph.IM] 10 May 2016 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 1994Lazarian , 1995 as a result of gas streams. This paper is structured as follows: First, we present the PO-LARIS 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 PO-LARIS 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 PO-LARIS, we calculate synthetic intensity and polarization maps in Sect. 3. Here we considered both analytical models and complex 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.

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 by an 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.

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 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, T d is the temperature of the dust, and C abs,λ the cross-section of absorption. Afterwards, the photon package propagates a distance l i in the grid between two points of radiation-dust interactions. At each cell wall, the absorption rate per time uniṫ 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 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 C sca /C ext . In the case of scattering (C sca /C ext <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 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 2.3) are pre-calculated with DDSCAT (Draine & Flatau (2013), see Sect. 2.2.4 and Sect. 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 T d 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 i − th cell is determined by 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.

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 crosssections 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 crosssections (C = πa 2 × Q) parallel C ext,|| and perpendicular C ext,⊥ 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: where T d and n d are the temperature and number density, respectively, of the dust. In this paper X denotes the averaging over grain size distribution and X over orientation. So far, PO-LARIS 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 (C ext,|| ) and major axis (C ext,⊥ ), and are calculated for each dust grain size in each RT simulation with The angle ϑ is between incident light and magnetic field direction (see Fig. 1). For oblate dust grains C ext = 0.5(C ext,x + C ext,y ), and C ext = 2C ext,|| + C ext,⊥ /3.
The Rayleigh reduction factor R (Greenberg 1968) is introduced to take care of imperfect grain alignment and is defined as 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 cos 2 (x) = f (x) sin(x)dx. 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 Sect. 2.2.2, 2.2.3 and 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 with a correlation factor f c where f c = 0 stands for no correlation (Hoang et al. 2013. The correlation factor is of minor influence to the net polarization, meaning that this approach is well justified. Once, the Rayleigh reduction factor is known, the exact crosssections for the RT simulation can be calculated with Eqs. 10 and 11 for a given size distribution n(a) between a minimal a min and maximal a max grain size with The same procedure is used for the cross-sections of scattering C sca , absorption C abs , and circular polarization C circ to obtain the matrix required in Eq. 7.

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 (Lazarian 1996), where h = I || /I ⊥ is the ratio of the momenta of inertia, J e f f is the rms value of the angular momentum, and N is a normalizing constant (see App. B for details).

Imperfect Davis-Greenstein (IDG) alignment
For the IDG, the alignment mechanism is paramagnetic relaxation (Davis & Greenstein 1951;Jones & Spitzer 1967; Purcell 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 v g and B and between the anisotropy u λ in the radiation field and B, respectively. The vector v rel is the projection of v g on the direction of the magnetic field direction B. 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 where T g is the gas temperature and n g 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 (Aannestad & Greenberg 1983), with and subsequently the Rayleigh reduction factor. To calculate the internal alignment, we use the approximation from Lazarian & S. Reissl, S. Wolf and R. Brauer: Radiative transfer with POLARIS: I. Analysis of magnetic fields through synthetic dust continuum polarization measurements Roberge (1997) for the angular momentum where J 2 th ≈ 2I || k B T g 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 (T d /T g = 1). However it slightly underestimates the internal alignment for s < 0.5 and T d /T g < 0.2.

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(Lazarian , 1995(Lazarian , 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 cos 2 (β) is completely defined by the anisotropy in the gas stream where v rel is the part of the relative drift velocity v rel projected on the magnetic field direction and g = h −1 is the grain nonsphericity (Lazarian 1994). In the case of oblate dust grains, g is negative and 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 with the molecular mass µm H of the gas. For the required conditions of supersonic velocity we locally calculate the Machnumber in each cell with and demand M > 1. Otherwise, we consider the dust grains to be randomized. Possible mechanical alignment of helical dust grains in subsonic environments is currently being discussed (e.g. Lazarian & Hoang 2007) but is not yet in POLARIS.

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 (a alg ) at which the dust grains start to align, we follow the procedure described in Bethell et al. (2007) with 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 shapedependent parameter δ, as well as the characteristic timescales of gas bombardment t gas and radiation t rad 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 u λ 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 This allows us to calculate the mean energy density and the anisotropy parameter determined by 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 G low−J cos 2 (ζ) remains of relevance. When we consider f high−J to be the fraction of dust grain aligned at the high-J attractor point the Rayleigh reduction factor is So far, f high−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 (t L ) to be lower than the gas dumping time (t gas ): t L < t gas (Hughes et al. 2009;Hoang & Lazarian 2009), which gives us a limit of 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.

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 4x4 scattering matrixM by S =MS . 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 matrixM and asymmetry factor g can be precalculated 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 Sec. 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 nonspherical imperfectly aligned dust grains, we apply the following approximation:M The same approximation is used for the orientation and sizedependent asymmetry parameter g of the HG-Phase function (see Eq. 4): where M i j (a, λ, ϑ , ϕ ) 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.

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, 1997Weingartner & 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: RTheating 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.

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 reason-  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 a alg in the RAT mode. 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.

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 PO-LARIS' 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 a max = 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 171500 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 P l is determined by the Stokes vector with and its projected position angle χ 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.

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 a max = 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 C abs,λ . 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.

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.

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: where Here, r and z are cylindrical coordinates, h 0 is the scale height at r 0 and the quantities α and β are parameters that characterize the radial density profile and the disk flaring. We adjusted n 0 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 n ∝ B −κ between magnetic field strength and We consider different disk masses and refer to the models with a total disk mass of 10 −3 M as D03, with 10 −4 M as D04, with 10 −5 M as D05, and with 10 −6 M as D06, respectively.
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 a min = 5 nm and a max = 2 µm, respectively with a power-law size distribution of dn(a) ∝ a −3.5 da.
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 (T g = T d ). 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 a alg at which the dust grains start to align.

Results
Fig . 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 . 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) and Cho & Lazarian (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 a alg 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 (P l > 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 (n ∝ B −κ ). Hence, the magnetic field no longer satisfies the minimal required field strength (see Eq. 31) in the thinner disk regions.

Discussion
Aligned dust grains as an explanation of polarized radiation in the FIR 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); ) 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 f high−J = 0.25 would still result in a detectable degree of linear polarization. However, f high−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.

Cloud model
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 10 12 m −3 and 10 13 m −3 , respectively. Density and temperature are chosen to be consistent with theoretical models and observations (Launhardt & Henning 1997;Sipilä et al. 2011).

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 −8 T (Falgarone et al. 2008;Beuther et al. 2010). To compare calculations performed by PO-LARIS 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 f c = 0.6 and assume f high−J = 0.5. The initial dust temperature is adjusted with a single star as radiation source (see Tab. 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 n d ∝ a −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 T d = 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 T g = 10 × T d . In a final step, we applied the RAT mode of POLARIS to calculate the distribution of the size cutoff a alg 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.   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 a alg decreases and dust grains with radii above a alg 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 u remains low. The expected low grain alignment in the area D cannot be found. This is because the overall energy density u outweighs the angle dependency in this very thin region. In Fig.6 (middle bottom panel), the low degree of grain alignment in the areas A − D 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 Article number, page 11 of 17 A&A proofs: manuscript no. paper of sight and this, subsequently, changes the polarization direction. A correlation between intensity I and polarization P l was recognized in Henning et al. (2001) following the power law P ∝ (I/I max ) −α 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 P − I relation is hyperbolic in nature, we expect to find a flat tail when I/I max 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 P − I relations for our cloud models rotated by 90 • . The P − I 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 P − I relation in the RAT case is less definitive. In contrast to the IDG alignment, the P − I relation of RAT alignment starts with an increase in polarization up to a point of I/I max ≈ 0.1 and then remains almost constant at its maximum value. However, the RAT alignment mechanism matches the expected P − I relation curvature in agreement with observations for larger fractions of I/I max . The GOLD alignment shows an exceptional behavior. In contrast to observational findings, the degree of linear polariza-tion 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.

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 Eq. 23 and Eq. 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); . 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 f high−J and the upper cut-off a max 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 P − I relations. Both RAT and GOLD alignment show an unexpected distribution of linear polarization. This uncertainty and the devi-ations in calculated P − I relations may be resolved with more sophisticated models from dedicated MHD simulations in forthcoming parameter studies.

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  ized 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). 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/.  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), along the mid-plane of the disk (middle), and the distribution along the disks surface (right). 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.  Table 1