Issue 
A&A
Volume 645, January 2021



Article Number  A143  
Number of page(s)  18  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202039133  
Published online  27 January 2021 
The scattering order problem in Monte Carlo radiative transfer
University of Kiel, Institute of Theoretical Physics and Astrophysics, Leibnizstrasse 15, 24118 Kiel, Germany
email: akrieger@astrophysik.unikiel.de
Received:
10
August
2020
Accepted:
4
November
2020
Radiative transfer simulation is an important tool that allows us to generate synthetic images of various astrophysical objects. In the case of complex threedimensional geometries, a Monte Carlobased method that simulates photon packages as they move through and interact with their environment is often used. Previous studies have shown, in the regime of high optical depths, that the required number of simulated photon packages strongly rises and estimated fluxes may be severely underestimated. In this paper we identify two problems that arise for Monte Carlo radiative transfer simulations that hinder a proper determination of flux: first, a mismatch between the probability and weight of the path of a photon package and second, the necessity of simulating a wide range of high scattering orders. Furthermore, we argue that the peeloff method partly solves these problems, and we additionally propose an extended peeloff method. Our proposed method improves several shortcomings of its basic variant and relies on the utilization of precalculated sphere spectra. We then combine both peeloff methods with the Split method and the Stretch method and numerically evaluate their capabilities as opposed to the pure Split & Stretch method in an infinite planeparallel slab setup. We find that the peeloff method greatly enhances the performance of these simulations; in particular, at a transverse optical depth of τ_{max} = 75 our method achieved a significantly lower error than previous methods while simultaneously saving > 95% computation time. Finally, we discuss the inclusion of polarization and Miescattering in the extended peeloff method, and argue that it may be necessary to equip future Monte Carlo radiative transfer simulations with additional advanced pathfinding techniques.
Key words: methods: numerical / radiative transfer / scattering / dust, extinction / opacity / polarization
© ESO 2021
1. Introduction
Radiative transfer (RT) simulations based on the Monte Carlo (MC) method play a leading role in the calculation of the radiative field and corresponding observable properties of a wide variety of astrophysical objects. They enable simulations of complex systems, such as protoplanetary or circumbinary disks, Bok globules, filaments, or planetary atmospheres, which are subjects of studies of various Monte Carlo radiative transfer (MCRT) codes, for instance, RADMC3D (Dullemond et al. 2012), POLARIS (Reissl et al. 2016), and SKIRT 9 (Camps & Baes 2020). The core principle of this MC approach is the simulation of photon packages that follow individual probabilistically determined paths through the model space in which they interact with matter, like dust or gas. This approach can then be used to calculate temperature distributions and flux maps of the simulated object.
However, the computational demand for these simulations rises significantly with the optical depth of the system, and temperature and flux estimates may, despite this, show high levels of noise (e.g., Baes et al. 2016; Gordon et al. 2017; Camps & Baes 2018; Krieger & Wolf 2020). Several methods that aimed to solve these problems at high optical depths have been proposed, for instance, the partial diffusion approximation (Min et al. 2009), the modified random walk (Min et al. 2009; Robitaille 2010), biasing techniques (Baes et al. 2016), and precalculated sphere spectra (Krieger & Wolf 2020). However, despite the success of these methods in their respective areas of application, the overall problem of high optical depth has still not been solved.
In general, these methods assist at least one of two types of MCRT simulations: the calculation of temperature distributions and/or flux maps. The aim of this paper is to study the latter in greater detail. Gordon et al. (2017) conducted a study, in which the outcome of different RT codes for a simple threedimensional setup were compared. The setup is composed of a homogeneous slab of dust that is illuminated by a single blackbody source. Already at a transverse optical depth of τ_{max} = 30 the accounted transmitted flux differences between these RT codes reached ∼100%. This was further investigated by Camps & Baes (2018), who called this the failure of MCRT at medium to high optical depths. In order to grasp the core problem that leads to this failure, they further simplified the setup to an infinite planeparallel slab that is illuminated on one side and study the intensity of radiation that is transmitted through the slab. A comparison of the basic MC approach and two biasing techniques, the Split method and the Stretch method, showed the best results when both biasing techniques are used simultaneously. However, despite the relative success of the Split & Stretch (S&S) method, the equivalent serial computation time that was required to properly simulate a transverse barrier of τ_{max} = 75 reached 52 days. The authors concluded that a proper simulation at high optical depths requires an enormous amount of simulated photon packages, and they suggested equipping MCRT simulations with statistical methods to measure the convergence of their results.
In this study we further investigate the core problem of MCRT simulations at high optical depths and demonstrate how the peeloff method greatly enhances their performance. In Sect. 2 we briefly describe the setup for our simulations, which is adopted from Camps & Baes (2018). We then examine and identify the main difficulty MCRT simulations face in the domain of high optical depths. Based on these results, we give reasons for the great success of the S&S method but also show its shortcomings. In Sect. 3 we argue that the peeloff method bridges one of the core problems: the need to calculate a broad range of scattering orders for a proper flux estimation. Additionally, we introduce the concept of the extended peeloff method, which may greatly boost the computation speed of a MCRT simulation, that is based on the basic peeloff method. Moreover, we describe how emission properties of spheres, which are crucial for this method, as well as their projections onto a detector can be precalculated in order to speed up MCRT simulations. We then argue why using precalculated emission spectra is preferred to using spectra based on the modified random walk (MRW; Min et al. 2009; Robitaille 2010). In Sect. 4 we apply the basic peeloff method as well as the extended peeloff method to the infinite slab setup, and compare their results to those of the S&S method. In this context we discuss the advantages and potential disadvantages of using the peeloff method. In particular, we look into its capabilities of performing proper high optical depth simulations with the use of Miescattering and a full treatment of the polarization state of photon packages. Finally, in Sect. 5 we briefly summarize our results.
2. Setup and problem identification
In this section the infinite planeparallel homogeneous slab setup is described, and a MC based probabilistic numerical implementation and a nonprobabilistic implementation are summarized. In Sect. 2.3 two core problems of MCRT simulations at high optical depths are demonstrated on the basis of the slab experiment. Subsequently, the success of the S&S method at optical depths τ_{max} < 75 and its failure at higher optical depths is explained, highlighting the need for more advanced MCRT methods, which will be discussed in detail later in Sect. 3.
2.1. Setup: Infinite planeparallel slab
To study the performance and capabilities of MCRT simulations at high optical depths, we simulate an infinite planeparallel homogeneous slab, which is illuminated on the top side by an isotropic monochromatic source of radiation. Our goal is to study this setup using various biasing techniques, which are described in this section, and compare them to the methods described by Camps & Baes (2018). To ensure comparability, we adopt the setup, i.e., a homogeneous slab with an albedo of A = 0.5 under the assumption of isotropic scattering, and simulate the transmission through the slab. For the implementation we make use of optical depth units, which are the natural units of this problem, and simulate slabs of different transverse optical depths τ_{max}. Since this setup was described by Camps & Baes (2018), we narrow our description down to the key features.
During the simulation N photon packages are created that penetrate the top layer of the slab and travel through it, until they leave one side of the slab. The initial direction of any photon package is chosen randomly according to a probability density function (PDF) given by
where μ_{0} = cos(θ_{0}) is the cosine of the initial penetration angle θ_{0}; in particular, μ_{0} = 0 and μ_{0} = 1 correspond to a direction parallel to the slab’s surface layer and perpendicular to it, respectively. A sketch of the setup can be seen in Fig. 1. Upon entry, a photon package travels a distance τ that is chosen randomly according to the PDF
Fig. 1.
Sketch of the infinite planeparallel slab setup. Radiation enters the slab of transverse optical depth τ = τ_{max} through the top layer. The observer detects radiation that is transmitted through the slab and is located below the bottom layer of the slab. The initial direction of a photon package is described by the penetration angle θ_{0}. If the photon package scatters, its direction is chosen randomly according to its corresponding PDF and the new direction is described by the angle θ. 
At the resulting new position the photon package has a probability of absorption (1 − A) and a probability of scattering (A). In the case of a scattering event, a new direction μ is assigned to the photon package. The PDF for isotropic scattering is given by
The steps of translation, defined by Eq. (2), and redirection, defined by Eq. (3), of scattered photon packages are repeated until the translation distance τ of the photon package is large enough to leave the slab through one of its boundary layers. Photon packages that penetrate the bottom layer of the slab contribute to the transmission intensity. In order to estimate the contribution of a photon package to the transmitted intensity, it carries a weight w, which is modified during the MC procedure and eventually adds up at the detector to the overall transmitted intensity. During its creation, a photon package is neutral, meaning that its weight is initialized with w = 1. For a significant boost in performance we use the S&S method (Camps & Baes 2018). When the Split method is used, the weight of the photon package is reduced for every event of interaction by the factor w^{split} = A. A photon package with reduced weight then continues to propagate through the slab. This biasing technique prohibits an early removal of the photon package due to absorption. When the Stretch method is used, the traversed path length is chosen according to the PDF for compositebiasing (Baes et al. 2016), which is given by
where T is the distance to the surface layer of the slab in the current direction of the photon package. In order to compensate for this modification of the path length, the weight of the photon packages needs to be readjusted by the factor
The detector is defined such that transmitted radiation is sampled with respect to its direction of propagation (given by the variable μ) when leaving the slab, resulting in a sampled directional transmission spectrum. The detector is composed of M equalsized bins for the directions μ that cover the interval [0,1]. The jth detector bin of the detector measures the intensity I_{j} given by
where is the contribution (i.e., weight) of the nth photon package and μ_{j} the central value of the jth detector bin.
2.2. Nonprobabilistic numerical implementation
The conceptually simple problem of determining the transmission of radiation through an infinite planeparallel slab (see Sect. 2.1) has been studied extensively in the past. This setup is particularly interesting as it is challenging for the MC method, but at the same time solvable with a nonprobabilistic method, making it a suitable testbed for MCRT simulations. Of the different approaches to solve this problem, the spherical harmonics method has shown to be particularly efficient (Baes & Dejonghe 2001). For comparability, we adopt the method for solving this problem described by Camps & Baes (2018), which is based on work done by Roberge (1983) and its optimization proposed by di Bartolomeo et al. (1995). This method is based on a solution of the integral
where τ is the transverse optical depth and Ψ(μ, μ′) the normalized redistribution function for incoming (μ′) and outgoing (μ) radiation. Unless specifically stated otherwise, we assume an albedo of A = 0.5 and isotropic scattering (i.e., Ψ(μ, μ′) = 1) throughout this paper. When assuming no additional internal sources, this differential equation is determined by its boundary conditions. In particular, the method relies on the series expansion of I(τ, μ) in Legendre polynomials, which is truncated at some odd order L. The solution of I(τ, μ) is then valid for (L + 1)/2 directions μ, which correspond to the zeros of the Legendre polynomial of order L + 1. For comparability, we chose L = 81. A more detailed description of this method can be found in the works of Roberge (1983), di Bartolomeo et al. (1995), Baes & Dejonghe (2001), and Camps & Baes (2018). In later sections of this paper this method is referred to as the nonprobabilistic solution to the infinite slab experiment.
2.3. MCRT simulations at high optical depths
High optical depths have been proven to cause major difficulties for MCRT simulations (e.g., Baes et al. 2016; Gordon et al. 2017; Camps & Baes 2018; Krieger & Wolf 2020). Camps & Baes (2018) attempted to investigate the cause of this problem and find that at high transverse optical depths results of MCRT simulations strongly depend on the number of simulated photon packages. Additionally, the need for a theoretical understanding of this issue was highlighted and the use of statistical tests for convergence was recommended. The aims of this section aims are the deepening of our understanding of the core problem of MCRT simulations at high optical depths, and thus laying the foundation for a potential solution to it. To this end, we first give a brief description of flux determination in MCRT simulations before particularly analyzing the slab experiment.
Photon packages are emitted from one or more sources and subsequently follow their individual paths through the model space before eventually being detected. In order to determine flux accurately, in an ideal simulation, every possible path needs to be taken into account. Since this is not feasible in a realistic simulation, the number of simulated photon packages is limited. The complexity of potential paths increases with the number of interactions it undergoes. In other words, the volume of paths increases with the number of potential interactions. However, the contribution of a path to the measured intensity is highly dependent on its particular properties. In the MCRT simulation described in Sect. 2.1 the contribution of a path is represented by the weight of the corresponding photon package. In order to determine a correct solution it is particularly important to simulate paths with high weights. To identify such paths it is mandatory to identify factors that reduce the weight of a path, i.e., that pose some form of resistance on the path of the photon package. In the case of the infinite planeparallel slab, for instance, it is obvious that the paths of least resistance that contribute to detector bins of low μ values are those that quickly reach the bottom layer of the slab and then scatter in the direction of the corresponding bin. This implies that the travel direction of least resistance does not necessarily coincide with the direction of the detector. In a complex setup it is not clear at the outset which are the paths of least resistance. As a testbed, the slab experiment allows us to make general statements regarding the interplay of albedo A, scattering order n, and the traveled path length Δl. In particular, we study the contribution of different n to the total intensity. Solutions for the transmitted radiation through an illuminated slab have already been derived in the past (e.g., Adamson 1975). However, we derive such solutions with a novel approach (to our knowledge) for the onedimensional slab in Appendix A as a new perspective may have the potential to lead to new insights.
Detected photon packages traveled through the onedimensional slab, scattered in general several times, and eventually left the slab toward the detector. The probability p^{(n)} that a photon package undergoes n scattering events before its detection is given by
where I^{(n)} and I^{tot} are the nth order contribution and the total intensity, which can be calculated using Eqs. (A.26) and (A.30), respectively.
Figure 2 shows p^{(n)} for the case where A = 1 (i.e., for an absorptionfree medium) for five different transverse optical depths τ_{max} ∈ {5,10,20,40,80}. As described in Appendix A.5, we used analytical solutions of the first few orders of Eq. (A.32) as a measure for convergence, and truncated the sum when an accuracy of 0.1% was reached. In the case of τ_{max} = 80, however, we truncated the sum at the order k = 10^{9}, before reaching its convergence criterion for n ≤ 4. However, the displayed results show only solutions for significantly higher scattering orders, where the sum converges faster (see Appendix A.5 for details). Figure 2 shows a clear increase in the scattering order of highest probability n^{max} for increasing optical depths. The mean scattering order of detected photon packages, displayed as dashed vertical lines in the plot, was calculated using Eq. (A.34), and shifts likewise with increasing τ_{max} toward higher scattering orders and increasingly away from its corresponding n^{max} value. This shift is caused by the increasing contribution of higher scattering orders for increased optical depth (i.e., for n > n^{max} the slope of p^{(n)} increases with increasing τ_{max}). Furthermore, we find that low scattering orders n ≤ n^{max} are suppressed since they require a rapid displacement of the photon package toward the observer; in other words, they exhibit a distinctly directed overall movement. However, high scattering orders n ≥ n^{max} are also suppressed because these paths are confined to the volume of the slab (i.e., they experience an overall highly undirected displacement). This problem can be viewed as a pathfinding problem. Moreover, we find that for higher transverse optical depth τ_{max}, the range of contributing scattering orders increases strongly, which we refer to as the scattering order problem. In the case where A < 1, the contribution of the nth scattering order is linked to the contribution displayed in Fig. 2 by multiplying p^{(n)} with A^{n} and rescaling. This results in a shift of n^{max} toward smaller n values, meaning that the lower the albedo, the more quickly the photon package has to traverse the slab in order to have a significant contribution to its transmission spectrum. In the case where A = 0.5, for instance, the mean number of interactions a detected photon undergoes for τ_{max} = 80 is about N^{mean} ≈ 30, while for the same slab and A = 1 the mean number of interactions is ∼1119.
Fig. 2.
Detection probability for transmitted photon packages of scattering order n through an onedimensional slab. Plotted are the results for five different transverse optical depths τ_{max} ∈ {5,10,20,40,80} under the assumptions A = 1 and using isotropic scattering. Vertical dashed lines mark the mean scattering order of the corresponding optical depths. For details, see Appendix A.5. 
On the basis of these considerations it is possible to specify the advantages and drawbacks of the compositebiasing technique. The main effect of this technique is to increase the average distance a photon package travels between two consecutive scattering events. This affects the path in two ways: first, by reducing the overall distance the photon package travels before being detected and second, by reducing the scattering order of the photon package. This can be seen by considering the path of a photon package that originates in the center of a dense homogeneous sphere, which travels and scatters multiple times inside it before leaving the sphere through its rim. The number of interactions N_{int} it undergoes before leaving the sphere increases quadratically with the optical thickness τ_{sph} of the sphere: (e.g., Krieger & Wolf 2020). Increasing the traveled distance by some constant factor f_{x} is equivalent to reducing the optical depth of the sphere by that factor (i.e., N_{int} ∼ (τ_{sph}/f_{x})^{2}), thus it is effectively reducing the scattering order quadratically. The traveled distance of the photon package inside the dense sphere is on average given by ΔT = N_{int}Δτ, where Δτ is the average traveled distance between two consecutive events of interaction and satisfies Δτ = f_{x}. Therefore, the traveled distance satisfies meaning that it is reduced by the factor f_{x}. As a result the compositebiasing technique generates paths of low scattering orders that reach the detector in a relatively short distance. Nonetheless, it is clear that a too high value of τ_{max} and a high value of A may lead to a situation in which the generated paths do not result in the desired range of scattering orders. Furthermore, it seems inevitable that very high τ_{max} and A values require the simulation of a broad range of observed scattering orders, as seen in Fig. 2. In order to sample photon package paths according to their importance to the total observed intensity it is thus mandatory to generate high numbers of photon packages; in addition, they should exhibit a high number of scattering events. This, however, adversely affects MCRT simulations that are equipped with the S&S method.
3. Methods
3.1. Basic peeloff approach
The S&S method is expected to encounter issues at high optical depth and high albedo. In particular, Sect. 2.3 demonstrates the necessity for covering a wide range of scattering orders, which is not feasible with the S&S method alone. There is, however, a promising method for solving this scattering order problem: the peeloff method. This method is already used by many MCRT codes, for example, Mol3D (Ober et al. 2015) and POLARIS (Reissl et al. 2016). It simply generates an additional photon package, the peeloff photon package, at every point of interaction of the original photon package. The peeloff photon package is then sent immediately (i.e., interactionfree) to the detector and carries a weight that is calculated according to the probability for the original photon package to be scattered in this direction. This weight is also reduced due to the extinction optical depth between the point of interaction and the observer. The weight of the peeloff photon package which is detected by the jth detector bin is given by
where w_{org} is the weight carried by the original photon package and τ_{tv} is the transverse optical depth between the point of interaction and the bottom boundary of the slab. Since every scattering event of a path is represented by an individual photon package, every simulated photon path increments the measured flux at the detector for a range of scattering orders. In particular, the range of simulated scattering orders ranges from zero to the total number of scattering events the photon package undergoes before leaving the slab. In later sections of this paper this peeloff method is referred to as the basic peeloff method. In Sect. 4 we present results for a comparison between the S&S method and the Split & Peeloff & Stretch (SPS) method.
3.2. Extended peeloff
According to the findings in Sect. 2.3 it is clear that the basic peeloff method has its computational limits since the required number of simulated scattering events rapidly grows with τ_{max} (i.e., with the volume of the slab), as can be seen in Fig. 2. Therefore, it is expected that the SPS method fails to properly simulate arbitrarily high transverse optical depths and a more advanced method is needed. In order to compensate for the increased complexity of photon paths at higher τ_{max}, we propose a novel method, which we call the extended peeloff method. Later we combine this method with the S&S method and compare this Split & extended Peeloff & Stretch (SePS) method against the SPS method. Instead of generating a single peeloff photon package at a point of interaction, this method creates several peeloff photon packages that originate from an extended region. When the method is used, a homogeneous sphere needs to be defined with the point of interaction defining its center. Subsequently, the photon package is moved immediately to the rim of the sphere and the emission properties of the sphere are used to emit this original photon package and to generate peeloff photon packages at different points on the surface of the sphere. These features eventually enhance the performance of a MCRT simulation and improve the quality of its outcome (see Sect. 4). For this method it is crucial to precalculate the emission properties of spheres of certain sizes since they are used for the photon package path determination and for the flux calculation. This method aims to solve a few of the problems mentioned in Sect. 2.3. First, emission spectra of spheres have to be precalculated only once and can then be used repeatedly during one simulation and even across different simulations. Therefore, if sufficient time is spent on the calculation of the spectra, all subsequent MCRT simulations will benefit. Second, photon packages move more quickly through the model space since they are repeatedly emitted from extended spheres, which saves computation time and allows the number of simulated photon packages to be increased. Third, since sphere spectra are composed of photon packages of different scattering orders, any peeloff photon package is also composed of contributions from multiple scattering orders. In particular, if the original photon package scatters for the nth time, the generated peeloff photon package contains information of all orders n_{po} ∈ {n,n+1,n+2,…}, unlike the basic peeloff, which generates a pure nth order photon package. This is also the case for the original photon package that is emitted from the surface of the sphere. Moreover, Fig. 2 shows an increasingly negligible contribution of low scattering orders with increasing optical depth. Therefore, simulations that are based on the basic peeloff method are required to simulate a high number of weakly contributing scattering orders before reaching the range of highly contributing scattering orders. Simulations that use the extended peeloff method, on the contrary, contain contributions of arbitrarily high scattering orders already from the first application of the method. As a result, the extended peeloff method results in an increased amount of computational time spent on actually relevant interaction events.
Before describing the steps we took for the precalculation of spheres and before showing their results, we briefly comment on the MRW and relate it to the extended peeloff method. The MRW is similar to the extended peeloff method insofar as it aims to solve the problem of high optical depths by using predetermined solutions of a sphere’s spectrum and by quickly moving photon packages from the center to the rim of that sphere. In theory, its analytically obtained emission spectra could be used for the extended peeloff method and, consequently, the numerical MC based precalculations of these spectra could be avoided. However, there are two reasons why we decided against it. First, the MRW is limited in the type of scattering it is based on. Originally, it was developed for isotropic scattering (Fleck 1984) and then extended to anisotropic scattering (Ishimaru 1978) where the scattering probability only depends on the angle between the incident and outgoing particle direction. A further extension to Miescattering (Mie 1908), however, is not possible, making a proper treatment of polarization impossible. Secondly, the emission properties of photon package off the surface of the sphere are not defined, rendering a precise extended peeloff method impossible too. In the following section, we describe how emission spectra of spheres can be precalculated and in Sect. 4.2 we discuss the inclusion of polarization to the extended peeloff method.
3.2.1. Precalculation: Emission spectra
We precalculate spectra of spheres with optical depth radii . Transverse optical depths below are well represented by basic peeloff. The size of the largest precalculated sphere depends on the largest optical depth that can be encountered during a MCRT simulation, which is τ_{max} = 300 in this paper. We note, on the one hand, that a higher number of precalculated spectra generally results in a more optimized simulation of photon package paths. On the other hand, it increases the time spent on precalculations. While a study of the optimal number and spacing of the values of τ^{sph} would lead to further improvement of the extended peeloff method, it is not part of this study. However, while the chosen number of precalculated spectra is rather low and thus unproblematic, the values chosen for τ^{sph} are very useful, as will become clear in the course of this paper. In order to precalculate emission spectra that are based on isotropic scattering, we make use of the symmetry to reduce the amount of data that has to be stored. In Sect. 4 we discuss a more general case, as well as the inclusion of polarization in the extended peeloff method.
The precalculation setup is always composed of a sphere of certain radius τ^{sph} and a source of luminosity 1 in its center. The particular unit of this luminosity is irrelevant, since we are solely interested in the relative luminosity that leaves the sphere. The luminosity is carried by a number of photon packages (we chose 10^{8}) that are generated by the source, potentially scatter multiple times at various optical depth radii τ_{r}, and then eventually leave the sphere across its surface at optical depth radius τ^{sph}. The smallest sphere (k = 1) is precalculated using the SPS method. Larger spheres (k > 1), on the other hand, use the SePS method by utilizing spectra of already precalculated smaller spheres. The luminosity that is emitted through the surface is stored in 100 linearly sampled μ_{s} bins, where μ_{s} ∈ [0,1] is the cosine of the penetration angle θ_{s} (i.e., μ_{s} = cos(θ_{s})). Here, μ_{s} = 1 and μ_{s} = 0 are parallel and perpendicular, respectively, to the corresponding outwardly directed normal vector of the sphere.
In order to apply the SePS method already during the precalculation of a larger sphere, a precalculated smaller sphere has to be embedded inside it (see Fig. 3). It is important to note that the SePS method in general can only be used if a photon package interacts and if the point of interaction and the closest border of the simulated homogeneous region are at a sufficiently large distance Δτ_{min}, in particular . During the precalculation this condition becomes . If the distance is not sufficient, we instead use the SPS method.
Fig. 3.
Sketch of a radiating sphere at radial position τ_{r} embedded in a larger sphere with radius τ^{sph}. Radiation originating on the surface of the embedded sphere penetrates the surface of the larger sphere about the penetration angle θ_{s}, and thus contributes to its emission spectrum. For details, see Sect. 3.2.1. 
Subsequent to a scattering event, there is a probability that the photon package will leave the sphere across its surface independent of the particular scattering direction. Hence, the peeloff luminosity can be calculated in every direction and assigned to the corresponding μ_{s}bin of the larger sphere in order to calculate the contribution to its emission spectrum. To do this we define the polar angle θ_{ν} in the reference frame of photon package, such that θ_{ν} = 0 and θ_{ν} = π correspond to the direction of the smallest distance and of the largest distance to the surface of the larger sphere, respectively. In Fig. 3θ_{ν} = 0 and θ_{ν} = π respectively correspond to the north and the south pole of the larger sphere. Additionally, we define ν = cos(θ_{ν}) as the cosine of this angle. Next we describe how the peeloff emission is calculated for the case of the SPS method and the SePS method.
The SPS method. When the SPS method is applied during the precalculation the peeloff emission of the photon package in the direction ν contributes to the larger sphere’s μ_{s}bin of direction
This implies that a scattering event at radius τ_{r} results in a nonzero contribution for all . Thus, only scattering events close to the surface of the sphere contribute to low μ_{s}bins. The calculation of the full contribution to a particular μ_{s}bin can be obtained by solving the integral
where is the distance to the surface in direction ν. The portion of w_{org}, which crosses the surface of the sphere in a specific μ_{s}bin is obtained by choosing the νlimits in the integral that correspond to the μ_{s}bin by using Eq. (10). An integral of Eq. (11) is given by
where Ei(x) is the exponential integral defined as
The SePS method. If the distance between a scattering position and the surface of the sphere exceeds , a smaller precalculated sphere can be placed inside. This scenario is depicted in Fig. 3. Hereafter the spectrum of the small sphere can be used to determine its transmission through the surface of the larger sphere it is embedded in.
The transmission spectrum of such an embedded sphere strongly depends on τ_{r}. However, its calculation involves the solution of a computationally expensive multidimensional integral since radiation is emitted from every point of the surface of the embedded spheres and in every direction of its corresponding hemisphere. As displayed in Fig. 3, the transmitted luminosity must also be sampled in the larger spheres μ_{s}bins. Therefore, before precalculating the spectrum of a larger sphere, we first precalculate the contribution of all embedded spheres to its emission spectrum for a set of selected τ_{r} values. In this way we can quickly update the emission spectrum of a larger sphere during its precalculation by accessing and interpolating these previously precalculated transmission spectra.
To ensure a proper sampling of precalculated τ_{r} values, we make use of the fact that any embedded sphere is always chosen to be as large as the geometrical limitations allow. Therefore, the deeper the scattering event of the original photon package occurs, the larger the embedded sphere is chosen. As a consequence, the interval for radial positions τ_{r} ∈ [0,τ^{sph}] can be divided into disjoint intervals, each corresponding to the size of an embedded sphere. In particular, larger values of τ_{r} correspond to smaller embedded spheres as the distance to the surface of the larger sphere is smaller. In the interval (τ^{sph}−,τ^{sph}], however, no sphere can be embedded, and instead the transmission spectrum is given by the basic peeloff by Eq. (12). As pointed out earlier, transmission spectra of embedded spheres only have to be evaluated within their corresponding τ_{r}intervals. First we determine the contribution to the emission spectrum for an embedded sphere at the infimum and the supremum of its corresponding τ_{r}interval. Considering that different μ_{s}bins may be accessible depending on the position of the embedded sphere, we then additionally calculate the contribution at those points where the number of accessible bins changes. In other words, the closer the embedded sphere is placed to the surface of the larger sphere, the more μ_{s}bins receive flux from the radiating embedded sphere. Using this particular sampling of evaluated transmission spectra ensures that during their interpolation μ_{s}bins receive a contribution only if they actually have a nonzero contribution. The number of evaluated radial positions for different sphere sizes τ^{sph} therefore differs, and in particular the number of evaluated positions inside spheres of size , , , and are 97, 161,199, and 221, respectively.
In Fig. 4, we show as an example for τ^{sph} = the results for the total contribution δL_{tot} to the emission spectrum. For this calculation we used 10^{8} photon packages. We note that the transmission spectrum is sampled with respect to the spheres’ μ_{s}bins; however, δL_{tot} is the sum of all contributions. The blue curve shows the total luminosity that is transmitted due to basic peeloff for a neutral photon package that has interacted at the radial position τ_{r}. In general, δL_{tot} increases with increasing τ_{r} and reaches its maximum of δL_{tot} ≳ 0.5 at τ_{r} = . Compared to the basic peeloff, the extended peeloff, which is shown as the orange curve, generally reaches higher values of δL_{tot}. This is a consequence of the fact that it takes into account higher scattering orders than the basic peeloff. Moreover, the curve shows discontinuities which result from the use of a discrete set of precalculated spheres. Deeper inside (i.e., for lower values of τ_{r}) larger precalculated spheres can be embedded. The first vertical gray line represents the highest value of τ_{r} at which spheres of size can be embedded. The second vertical gray line marks the highest radius for the embedment of a sphere of size . Past this radius none of the precalculated spheres fit inside and the basic peeloff method is used.
Fig. 4.
Total transmitted peeloff contribution δL_{tot} of radiating embedded spheres or points, inside a larger sphere of size = 10^{1.5}. The orange and the blue curve show the contribution when using the extended and the basic peeloff method, respectively. The three gray vertical lines show the transition radii τ_{r}, where a differently sized sphere can be embedded. In particular, the left gray line gives the transition from using a sphere of size to a sphere of size . The middle line separates the τ^{sph} region from the region where the basic peeloff is used. Close to the rim of the sphere (i.e., to the right of the last gray line for radii τ_{r} ≲ ) none of the precalculated spheres fit and the basic peeloff method is applied. 
The calculation of the transmitted spectrum and the corresponding total transmitted luminosity δL_{tot} is performed by the use of MC integration. That means, noninteracting photon packages are sent from the surface of the embedded sphere in a probabilistically determined direction and cross the surface of the larger sphere with a certain penetration angle. The weight these photon packages carry after traversing the required optical depth to reach the surface, is then added to the corresponding μ_{s}bin. In this case, the launching position on the surface of the embedded sphere is chosen isotropically and the direction of emission is chosen randomly according to the emission properties of the embedded sphere. The resulting contribution to the emission spectrum is eventually normalized, such that it corresponds to an originally neutral photon package. It is important to note that these results are only used for the precalculation itself, for which they are specifically performed. During a realistic MCRT simulation, only the emission spectra of the precalculated spheres will be used.
Precalculated sphere spectra. Figure 5 finally shows the obtained emission spectra of five precalculated spheres of optical depth radius for k ∈ {1,2,3,4,5}. The highest intensity of any sphere is emitted in the radial direction (μ_{s} = 1) and decreases toward lower μ_{s}. For this calculation we used the SePS method whenever possible. It is crucial to use a method like this to precalculate the spectra, since the use of the basic peeloff strongly underestimates the resulting emission spectra, which can be seen in Fig. 6. It shows the ratio of the emitted intensity for different sphere sizes when using the SePS method to the SPS method during the precalculation. Small sphere sizes (i.e., small values of k) do not show a strong deviation between the two methods. However, the relative differences increase for increasing sphere sizes and the intensity is clearly underestimated for the full range of μ_{s} values. The plot also shows in brackets the total estimated luminosity a sphere loses when utilizing the SPS method rather than the SePS method, which for k = 4 (i.e., = 100) already reaches ∼68% and for k = 5 (i.e., ≈ 310) even > 99%.
Fig. 5.
Precalculated emission spectra for spheres of optical depth radius for k ∈ {1,2,3,4,5}. The internal source of luminosity is neutral and isotropically radiating. It is placed in the center of the sphere and propagates through the sphere until it leaves it through its rim. The intensity I_{SePS} leaves the surface of the sphere and depends on the cosine of the penetration angle μ_{s} = cos(θ_{s}). The intensity is highest in the direction radially outward and overall decreases toward a smaller μ_{s}. 
Fig. 6.
Ratio of emission spectra I_{SePS} to I_{SPS}, which use the SePS method and the SPS method, respectively, during the precalculation. The SPS method reaches increasingly lower intensities for increasing sphere sizes (i.e., for higher values of k). Shown in brackets is the total loss of emitted luminosity when using the SPS method compared to the SePS method. For details, see Sect. 3.2.1. 
3.2.2. Precalculation: Projection onto detector
One crucial step in a realistic MCRT simulation is the determination of flux. To this end, a detector that measures incoming flux needs to be defined as an object of certain geometrical size and orientation. Most MCRT codes, in particular the previously mentioned codes, usually implement the detector as a geometrically confined plane object at a great distance with an orientation that can be characterized by a constant normal vector. The infinite slab setup introduced in Sect. 2.1, however, has a curved detector with nonconstant normal vector. In this regard, the detector of this specific problem is comparably complex.
First we describe a method for propagating extended peeloff emission from an embedded sphere to the detector. To this end, the surface of the sphere is divided into patches whose radiation is projected onto the detector. Subsequently, this method is applied to the detector of the infinite slab setup. The application to a flat detector at great distance is discussed toward the end of this chapter and is merely a simplified variant of the presented method.
Whenever the original photon package interacts and its distance to any boundary of the slab is sufficient, the point of interaction defines the center of a precalculated sphere of size for some value k that is chosen to be as large as possible. This sphere is then embedded in the slab, and every point of its surface radiates in its corresponding hemisphere. During the application of the peeloff method, radiation whose direction satisfies μ > 0 in the frame of the detector stays interaction free and is thus eventually observed by the detector. Since the sphere in embedded in a region of nonzero density, different points on its surface may result in different detected peeloff intensities, and the particular emission properties may also have angular dependencies. In order to account for these features, we divide the surface of the sphere into N_{p} patches. In order to properly represent the emission of a sphere by patches, they are required to be limited by a maximum deviation of the normal vectors and to be limited by a maximum extent in optical depth units. The first criterion is met if any two points p_{1} and p_{2} in a patch and their respective normal vectors n_{1} and n_{2} fulfill n_{1} · n_{2} ≤ C_{n} for some chosen constant C_{n} > 0. The second criterion can be satisfied by imposing the condition p_{1} − p_{2} ≤ C_{p} for another chosen constant C_{p} > 0. For increasingly large spheres, the first directional criterion is met easily as the curvature of the surface decreases with its radius. Therefore, the spatial extent of the patches becomes the dominant restriction. However, due to the azimuthal symmetry of the infinite slab setup, we apply a simpler criterion and choose patches such that any two surface points and with equal azimuthal position within a patch fulfill
where we choose . The condition Eq. (14) results in circular patches at the poles and ringlike patches between these poles, which are symmetric with respect to the transverse direction of the slab. These resulting patches have a constant extent in their polar angles. In the case of the sphere size , , , , and the number of patches is N_{p} = 12, 20, 36, 64, and 112, respectively.
It is possible to precalculate the luminosity that is emitted from any patch p in the direction of the jth detector bin. Every patch then emits this luminosity from its center, which is defined by the average μ value of the patch, and sends it toward the detector bins. It is convenient to use luminosity as a measure for the extended peeloff method rather than intensities since the weight of a photon package is associated with its carried luminosity. In a MCRT simulation this peeloff luminosity is additionally attenuated by the material on the path between the patch and the detector. The luminosity that the jth detector bin measures due to the peeloff emission of the patch p, is given by
where τ_{p} is the transverse optical depth between the launching point on the patch and the boundary of the slab. Figure 7 shows the example of a sphere of size and its total patch luminosity per patch size, which is a measure for the emitted flux of the sphere that is observed by the detector. It shows a strong increase in emitted flux toward higher values of μ_{p}, which is the μ value that describes the location on the sphere. In other words, patches that face the detector contribute more strongly to the detected intensity. However, this plot also shows that all patches contribute and how the μ_{p} range of different patches varies as a consequence of the definition in Eq. (14). We note that the area below the curve equals the total observed luminosity .
Fig. 7.
Emitted luminosity per patch size dL^{p, det}/dμ_{p} as a function of the position μ_{p} = cos(θ_{p}) on the surface of the sphere. The area beneath the curve corresponds to the total emitted luminosity of the sphere. Values of μ_{p} increase when nearing a faceon view of the detector; resulting in a rising extended peeloff contribution. 
In order to precalculate the projected luminosity of any patch onto the detector we use a MC method in which we generate 10^{7} photon packages on the surface of a precalculated sphere. The launching position is chosen randomly according to an isotropic distribution. Then the direction of a photon package is determined randomly according to the spectrum of the sphere, which has an angular dependence as can be seen in Fig. 5. Next, the μ value of that photon package is determined in the frame of the detector. If it is directed toward a detector bin, the luminosity that is carried by the photon package is added to the projected luminosity of the sphere. Due to the symmetry, half of these photon packages would not contribute to the projected intensity, since their traverse direction is opposite to that of the detector. To increase the performace and get a better estimate for the projected spectra even of patches that barely radiate toward the detector, we used a biasing technique and only chose photon package directions that are pointing toward the detector. To compensate for this the weight of these photon packages is reduced accordingly. Contrary to the precalculated transmission spectra described in Sect. 3.2.1 that were solely used during the precalculation, the projected peeloff luminosities that we precalculate and that are described in this section are repeatedly used during the actual MCRT simulation as we describe later in Sect. 4.
It should be noted that this elaborate method for calculating the projected luminosity is a consequence of the curvature of the detector. In the previously mentioned case of a flat detector at a large distance that can be described with a unique normal vector, the projected luminosity can be estimated with the flux in exactly that unique direction of the detector. Therefore, this case does not require a MC integration method.
In addition to the described extended peeloff method, we test a simplified version of the extended peeloff method, later referred to as the simple SePS method. In this variant, the whole sphere is represented by a single peeloff photon package, rather than N_{p} peeloff photon packages, and the luminosity leaving the sphere toward the jth detector bin (L_{j}) is obtained by assuming a great distance to the detector. Hence, L_{j} can be written as
where is the total luminosity emitted by the sphere. The peeloff photon package is then emitted from the point that is geometrically closest to and in the direction of the detector (i.e., the point satisfying μ_{p} = μ_{i}).
In Sect. 4 we present results for a comparison between the S&S method, the SPS method, the SePS method, and the simple SePS method.
4. Results and discussion
In this section the previously described SPS method (Sect. 3.1) as well as the SePS and the simple SePS method (Sect. 3.2) are applied to the infinite planeparallel slab problem (see Sect. 2). Furthermore, we present a brief discussion of the role of the albedo and of the extension of the SePS method with regard to anisotropic scattering and polarization. For clarity, a short overview with descriptions of the different applied MCRT methods is provided in Table B.1.
4.1. Role of peeloff in MCRT simulations
In Fig. 8 the results of simulations of the infinite slab with τ_{max} = 75 are shown. The dotted curve is the result of a nonprobabilistic numerical calculation (see Sect. 2.2) and is used for comparison with the results of different MC based methods. A total of 10^{9} photon packages are used in all MCRT simulations. The plot clearly shows that the three simulations based on the peeloff method result in overall smoother results for the transmitted flux as opposed to the S&S method, which results in a noisier transmission spectrum. This is a consequence of the fact that during every scattering event of an original photon package a peeloff based simulation sends individual peeloff photon packages toward every detector bin, whereas the S&S method solely increments the detected flux of a single bin per original photon package. The SPS method shows an overall better performance than the S&S method in the sense that its overall detected luminosity and also the shape of the result more closely resemble the nonprobabilistic solution. Nonetheless, both methods deviate by about one order of magnitude from the nonprobabilistic solution at low μ values. For increasing μ, these relative differences on average decrease. Our results indicate that using the SPS method results in an improved estimate of the transmitted flux. Nonetheless, the overall underestimation of flux shows that the SPS method does not solve the problem of high optical depths.
Fig. 8.
Results for the infinite planeparallel slab problem with transverse optical depth τ_{max} = 75. Different advanced MCRT methods (solid lines) are compared to results from a nonprobabilistic numerical method (dotted line). 
The SePS method and the simple SePS method achieved a significantly better result. The SePS method in particular shows a low relative error for all directions μ, and is thus a suitable method for transverse optical depths of up to ∼75. The simple SePS method, however, slightly overestimates the resulting flux. Toward higher values of μ this effect is anticipated, first because the complete emission of the sphere is represented by a single photon package when using the simple SePS method, and second because this peeloff photon package is launched at a position close to the pole of the sphere that is facing the detector. Therefore, radiation that would otherwise originate from points of the sphere that are more deeply embedded in the slab are now experiencing a lower attenuation, resulting in a higher flux level. This effect can be expected to intensify for larger spheres and additionally lead to an underestimation of flux at lower μ values for similar reasons. For this reason the simple SePS method fails to produce proper flux estimates at higher transverse optical depths.
Whether a method eventually succeeds or fails in the context of high optical depth MCRT simulations depends on the computation time it requires and the accuracy it achieves. The performance gain due to the SePS method then becomes apparent when considering the computation time required. All the presented simulations were performed serially (i.e., with a single core^{1}) using PYTHON. The calculation for the SePS method took ∼2 days compared to the ∼52 days of equivalent serial computation time reported by Camps & Baes (2018), who were using the S&S method with 10^{13} photon packages. The latter simulation could be expected to require a far greater number of photon packages to reach the same error as our SePS based simulation, which in addition used 10^{4} times fewer photon packages. As a result, the SePS method saved > 95% of the simulation time and achieved a significantly better result. Even though a direct comparison between the two simulations is difficult, as expected, the SePS method clearly outperforms the S&S method at medium to high optical depths and proves to be the preferred method for this type of MCRT simulations.
To study the feasibility of performing MCRT simulations at higher optical depths, Fig. 9 shows the results of four different MCRT simulations for three transverse optical depths τ_{max} ∈ {75,150,300}. In every simulation 10^{9} photon packages were used, which is the same number as used previously for the results in Fig. 8. In addition to the SPS and SePS method, we also consider the forward SePS and the forward SPS method. These two methods differ from the SePS method and the SPS method, respectively, only in that the Stretch method is applied if and only if the direction of the photon package after scattering points toward the detector (i.e., in forward direction). As a consequence, methods using the forward label result in detected photon packages experiencing on average a reduced number of scattering events and having traveled a shorter distance before leaving the slab. For all three values of τ_{max} the total detected luminosities of methods that use the forward label exceed their respective bidirectional methods. Already at an optical depth of τ_{max} = 150 the SePS method struggles to produce the correct result. The forward SePS method, however, is closer to the nonprobabilistic solution.
Fig. 9.
Results for the infinite planeparallel slab problem for transverse optical depths τ_{max} ∈ {75,150,300}. Different advanced MCRT methods (solid lines) are compared to results from a nonprobabilistic numerical method (dotted lines). The forward label indicates a variant of the SePS method that differs only in its application of the Stretch method; it is applied if and only if the photon package scatters toward the detector. 
The SePS method particularly underestimates the flux toward smaller μ values, which mainly originated from photon packages that interacted close to the bottom boundary of the slab. A lack of flux may therefore indicate a lack of representative scattering events close to this boundary. A potential solution to this problem may involve the introduction of a lower limit of scattering events that a photon package has to undergo before it can leave the slab. This is closely related to the problem of simulating the path of photon packages to sufficiently high scattering orders. It may very well be that too many photon packages are leaving the slab too early while still carrying a significant weight, both during the infinite slab simulation or already during the precalculation of sphere spectra. The fact that photon packages are allowed to leave the model space may at this point become a limiting factor as it systematically leads to an underestimation of flux. However, this raises the question of the order of the minimum number of considered scattering events, which may be the topic of future studies. Additionally, the results shown in Fig. 2 also imply an upper limit for relevant scattering orders; this means that photon packages that exceed this number of scatter orders may safely be removed from the model space without significantly changing the final result for the estimated fluxes. Another potential source for the underestimation of flux may reside in spectra that have not been precalculated with sufficient accuracy. To ensure a sufficiently high accuracy of the precalculated spectra, and thus the possibility to control this additional source of error, it is highly recommended to use statistical tests for convergence, such as those described by Camps & Baes (2018). In this study, however, the number of simulated photon packages is fixed.
At τ_{max} = 300, none of the considered methods is capable of reproducing the nonprobabilistic solution. These three plots also suggest that for increasing transverse optical depth the resulting flux level obtained by the SPS method decreases more strongly than the flux level of the SePS method, meaning that it becomes increasingly relevant to use an advanced method like the SePS method at higher optical depths. Nonetheless, at τ_{max} = 300 the estimated flux based on the SePS method is orders of magnitude too low. Only at μ ≲ 1 does the forward SePS method predict fluxes that are compatible with nonprobabilistic solutions for all three considered τ_{max} ≤ 300.
Therefore, we conclude that despite the success of the SePS method at an optical depth of τ_{max} = 75, it is nonetheless not fully capable of solving the problem of MCRT simulations at high optical depths and it is questionable whether this problem can be handled by purely MC based methods. However, we find that the SPS method and the SePS method both benefit MCRT simulations in the regime of medium to high optical depths where fluxes may be severely underestimated as they significantly alleviate this problem. The main reason for their success lies in the coverage of a wide range of scattering orders, which is especially effective in the case of the SePS method, designed to reach high scattering orders already during early stages of the life cycle of every photon package.
The role of the albedo. As described in Sect. 2.3, the albedo plays a decisive role in MCRT simulations at high optical depths. Similar to the effect of an increased optical depth τ_{max}, which can be seen in Fig. 2, an increased value of A shifts the scattering order with the highest contribution toward higher values. Simultaneously, the range of highly contributing scattering orders also increases. However, the probability of simulating a path that contributes significantly decreases with its scattering order. Consequently, keeping the number of simulated photon packages fixed, it becomes less likely to generate a sufficient number of paths that are contributing significantly to the measured flux. In conclusion, higher values of the albedo A require a greater number of simulated photon packages for a proper flux determination. Therefore, increasing the albedo results in a lower value of τ_{max} that can be simulated with acceptable effort; in turn, a lower albedo results in a higher limit for τ_{max}. This can be seen in Figs. B.1 and B.2, which show examples of the results of MCRT simulations for transmission spectra through a planeparallel infinite slab for an increased albedo of A = 0.9 and a decreased albedo of A = 0.05, respectively. Figure B.1 illustrates that the limit for τ_{max} even drops below τ_{max} = 75 when increasing the albedo to A = 0.9, as the determined flux is strongly underestimated. This can be explained by the fact that for A = 0.5 the limit for τ_{max} lies between τ_{max} = 75 and 150. Increasing the albedo significantly reduces this limit. However, Fig. B.2 shows the results of MCRT simulations for a reduced value of the albedo of A = 0.05, which is a typical value for dust in the micrometer wavelength range. It shows the results for transmission spectra through slabs of four different optical depths τ_{max} ∈ {75,150,300,600}. For this low value of A the τ_{max} limit increases strongly. In particular, only at an optical depth of τ_{max} = 600 do the results of the MCRT simulation using the SePS method start to show slight differences from the nonprobabilistic solution. In conclusion, we find that a high albedo of A = 0.9 has its limit at τ_{max} < 75, while a low albedo of A = 0.05 allows simulations up to τ_{max} ≳ 600. We conclude that increasing the albedo shifts the onset of the scattering order problem toward lower values of τ_{max}; in other words, it reduces the τ_{max} limit, which is in agreement with the findings in Sect. 2.3.
4.2. Extension to the SePS method: Anisotropic scattering and polarization
In this study, we introduced the extended peeloff method and specifically applied it to the infinite planeparallel slab problem studied by Camps & Baes (2018), which is based on purely isotropic scattering. In the following we discuss anisotropic scattering and polarization and how it can be integrated in the extended peeloff method.
HenyeyGreenstein (HG) scattering (Henyey & Greenstein 1941) is a widely used oneparameter model for anisotropic scattering. By integrating it into the extended peeloff method, as described in Sect. 3.2, the resulting spectra for precalculated spheres become in general dependent on the initial direction of the photon package. Depending on the type of the detector, there are a few possibilities for successfully implementing this dependency. In the case of a curved detector, as is used in the infinite slab problem, patches were defined on the surface of precalculated spheres, which were oriented in the coordinate system of the detector. To implement HG scattering, the emission spectra of these patches have to be calculated separately for the full sampled range of initial directions of the photon package. In particular, the initial μ = cos(θ) value of the photon package needs to be sampled. While in the isotropic case the emission spectra of any point of the surface depend solely on the parameter μ_{s} = cos(θ_{s}), the anisotropic case leads to a dependence on a second angle, ϕ_{s}, which describes the direction in the tangent plane of that point. Assuming a sampling of the angular parameters θ, θ_{s}, and ϕ_{s} in steps of 1° as well as single precision (4 bytes per entry), the required memory per wavelength and patch is on the order of only ∼200 megabytes. In the case of a plane detector at a distant location, these patches can instead be defined in the coordinate system of the scattered photon package. This removes the necessity of sampling the initial direction of the photon package and, as a consequence, significantly reduces the memory consumption to ∼1 megabyte per wavelength and patch. The implementation for both types of detectors is therefore readily feasible.
The inclusion of polarization to the extended peeloff method is possible as well, but leads to higher memory demands. Polarization is an important quantity that can be simulated in MCRT simulations. A proper treatment of polarization is often achieved by the use of Miescattering, whose outcome of scattering events depends on and affects the polarization state of a photon. It opens up many possibilities for indepth studies of astrophysical systems, providing complementary insights to pure fluxes. The SPS method can utilize polarization in a straightforward manner and already constitutes a crucial part of many MCRT codes like Mol3D (Ober et al. 2015) and POLARIS (Reissl et al. 2016). Compared to the inclusion of HG scattering, precalculated sphere spectra additionally depend on the initial polarization state of the scattered photon package. Precalculated spectra are meant to be calculated for neutral photon packages, where the weight for the initial photon package is associated with I_{in}, which is the first component of the Stokes vector. The remaining Stokes components Q_{in}, U_{in}, and V_{in}, however, need to be sampled, and they increase the memory consumption of precalculated spectra. Since memory demand may quickly pose a problem, the choice of sampled Q_{in}, U_{in}, and V_{in} is critical. In general, these parameters are limited by the inequality
and sampling them may be realized for instance by using a grid of cartesian or spherical coordinates. This means that whenever the extended peeloff method is applied the threedimensional bin corresponding to the Stokes vector of the scattered photon package is determined and its precalculated spectrum subsequently used. Assuming a cartesian sampling of these parameters in steps of 10% thus results in another factor of ∼4000 for the memory consumption. Moreover, photon packages that cross the rim of the sphere are also polarized, resulting in another factor of 4 in the memory consumption to store all components of their Stokes vectors. Overall, this results in ∼16 gigabytes per wavelength and patch under the assumption of a flat detector at great distance. Considering how coarse the grid in the Stokes components is, it suggests that a proper treatment of polarization using the extended peeloff method is still questionable as it is currently limited by its enormous memory consumption. However, given the rapid development of computer hardware, it can be expected that a full treatment of polarization based on this approach will become feasible in the future.
5. Summary and conclusions
In this work, we studied the problems of high optical depths in MCRT simulations and various potential solutions. For this purpose, we adopted from Camps & Baes (2018) the numerical setup of a threedimensional infinite planeparallel homogeneous slab and studied the transmission of radiation through it (see Sect. 2.1). To explore and extend possible solutions to this problem, we initially turned to a further simplified onedimensional version of the slab experiment (Sect. 2.3). This analysis was based on a solution of a differential equation for photon distribution eigenstates in the slab and a decomposition of an arbitrary distribution in terms of these eigenstates. Doing so we derived an expression for the contribution of different scattering orders n to the transmission spectrum through a slab in Appendix A.5, the results of which are displayed in Fig. 2. We demonstrated that proper flux estimations for increasing optical depths involve not only the simulation of high scattering orders, but also wide ranges of different scattering orders. The scattering order problem is accompanied by a pathfinding problem; the weight of the path of a photon package and its probability may strongly deviate from one another, in particular in the regime of high optical depths. We found that the albedo is a crucial variable that determines the severity of the two problems. In particular, an increasing value of the albedo increases the relevant range of scattering orders, and thus decreases the pace at which highly contributing photon packages need to travel through the model space. Based on these findings, and despite its performance boost, we argue that the S&S method will inevitably fail at high optical depths as it is designed to significantly reduce the scattering order of detected photon packages.
In Sect. 3 we identify the peeloff method as a promising candidate for solving the scattering order problem. When applied, the basic peeloff method repeatedly generates peeloff photon packages at every scattering point on the path of a representative photon package. Since peeloff photon packages are sent on an interactionfree path toward the detector, the detected set of peeloff photon packages comprises a range of scattering orders n ∈ {0,1,…,n_{max}}, where n_{max} is the number of scattering events the original photon package has undergone before leaving the model space. In Fig. 8 the SPS method, which combines the basic peeloff method with the S&S method, shows a slight improvement compared to the S&S method in the τ_{max} = 75 slab experiment. Nonetheless, just like the S&S method, the SPS method also fails to quickly and properly calculate transmission spectra at high optical depths for a few reasons. First, the full range of relevant scattering orders has to be calculated, meaning that increasing values of τ_{max} require increasing computational effort. Second, as seen in Fig. 2, detected photon packages of low scattering order become increasingly futile with higher optical depth. This means that a potentially high amount of computation time has to be spent on the calculation of mostly negligible scattering orders before stronger contributing scattering orders can even be calculated.
For these reasons, we proposed a novel method that aims to solve the scattering order problem by taking into account environmental information of the path of a photon package: the extended peeloff method. Rather than generating a peeloff photon package at a point of interaction and sending it to the detector, the interacting photon package is placed in the center of a homogeneous sphere, and a number of peeloff photon packages are then sent from the surface of the sphere to the detector. Subsequently, the original photon package is immediately moved to and emitted from the surface of that sphere. As a consequence, the extended peeloff method emits photon packages that are composed of contributions from multiple scattering orders. The emission properties of the sphere need to be precalculated and stored before this method can be used (see Sect. 3.2). This approach tackles multiple problems at once. First, starting from its first application the extended peeloff method generates high scattering order photon packages, and thus the range of weakly contributing scattering orders is quickly surpassed and the range of highly contributing scattering orders reached. Second, the immediate displacement of the photon package skips the calculation of a potentially high number of scattering events during a MCRT simulation. Figure 8 compares the SePS method, which combines the S&S method with the extended peeloff method, with the SPS method. Additionally, it shows results of the simple SePS method, a variant in which the sphere is represented by a single peeloff photon package. The figure clearly shows that the SePS method results in the most accurate transmission spectrum. Moreover, this simulation took only ∼2 days, and thus saved > 95% of the computation time compared to the equivalent serial computation time reported by Camps & Baes (2018), while at the same time achieving a significantly lower error (see Sect. 4 for details).
Furthermore, we discuss the benefits of precalculating emission spectra compared to using spectra that are used for the MRW procedure. The two main reasons are first, that the directional emission properties of the surface of a sphere are important for a proper peeloff approach, and second, that MRW is not capable of simulating polarization and Miescattering. In Sect. 4.2, the inclusion of anisotropic scattering or even polarization and Miescattering in the SePS method was described. While the inclusion of anisotropic scattering, for instance HG scattering, is possible, a full treatment of polarization is currently limited by the memory consumption of precalculated spectra and feasible only for rather coarsely sampled grids for the polarization states when using stateoftheart computers.
Furthermore, the limitations of the SePS method in terms of the transverse optical depth τ_{max} have been studied. Figure 9 shows results of the SePS method for three different optical depths τ_{max} ∈ {75,150,300}. In this context we also tested an additional variant of the SePS method that deviates only in that it applies the Stretch method if and only if the photon package scatters toward the detector, hence the name forward SePS method. We find that the forward SePS method results in overall higher transmitted luminosities, as was expected, but that both variants have difficulties properly simulating the spectra at τ_{max} ≥ 150 when using 10^{9} photon packages. We note that it is possible to simulate higher values of optical depths at the cost of increased numbers of photon packages (i.e., with increased simulation time). However, in this regime of τ_{max} the number of photon packages that is needed for a proper flux determination already reaches values that are too high to be handled with acceptable computational effort. We conclude that of all the tested methods, the SePS method showed the best performance for all studied transverse optical depth. This makes it the preferred MCRT method that has been tested in the regime of high optical depths even though it is not fully capable of solving all problems that arise in this regime. The benefits it grants to MCRT simulations, however, significantly alleviate the issue of underestimated flux. We also find that its success is mainly based on the coverage of scattering orders that it already provides during the early stages of the life cycle of every photon package.
According to the findings in Sect. 2.3, which are supported by the results in Sect. 4.1, it seems that the probability of generating significantly contributing paths decreases with τ_{max} and with A, meaning that the portion of paths with negligibly low weight is found to increase. As a result, the chances of generating highly contributing photon package paths also decreases. This implies that any uncoordinated MC approach with a limited number of photon packages is out of reach of being able to sample the paths sufficiently well as the optical depth increases. Moreover, it suggests that the higher the optical depth, the more the choice of the path of a photon package has to be coordinated and the application of some form of biasing may be crucial for this undertaking. For the determination of flux in a setup of high optical depths, it thus seems inevitable to equip MCRT simulations with pathfinding techniques. We believe that the problem of high optical depths may pave the way for a new bottom up approach to MCRT simulations; we believe that this new approach may lie in the use of artificial intelligence to guide photon packages on their way through the model space. Whether a method like reinforcement learning has its place in this undertaking remains to be seen. The simultaneous solution to the pathfinding and scattering order problem is the prerequirement to a full solution to the optical depth problem of MCRT simulations, and with this paper we have presented a method to solve the latter problem.
Interestingly, this result can also be obtained by discretizing the problem and describing the evolution of an arbitrary mdimensional initial photon distribution by a m × m transition matrix T. Due to the simplicity of the problem, the matrix is a Toeplitz matrix that is positive definite, which is a consequence of its dominant diagonal entries. As a result, due to the PerronFrobenius theorem, T has a positive largest (dominant) eigenvalue with a corresponding onedimensional eigenspace, which is spanned by a positive eigenvector.
Acknowledgments
We thank all the members of the Astrophysics Department Kiel for helpful discussions and remarks. We acknowledge the support of the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets (WO 857/171)”.
References
 Adamson, D. 1975, The Role of Multiple Scattering in Onedimensional Radiative Transfer [Google Scholar]
 Baes, M., & Dejonghe, H. 2001, MNRAS, 326, 722 [Google Scholar]
 Baes, M., Gordon, K. D., Lunttila, T., et al. 2016, A&A, 590, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Camps, P., & Baes, M. 2018, ApJ, 861, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Camps, P., & Baes, M. 2020, Astronomy and Computing, 31, 100381 [CrossRef] [Google Scholar]
 di Bartolomeo, A., Barbaro, G., & Perinotto, M. 1995, MNRAS, 277, 1279 [Google Scholar]
 Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC3D: A Multipurpose Radiative Transfer Tool (Astrophysics Source Code Library) [Google Scholar]
 Fleck, J. A., Jr., & Canfield, E. H. 1984, J. Comput. Phys., 54, 508 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gordon, K. D., Baes, M., Bianchi, S., et al. 2017, A&A, 603, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Henyey, L. G., & Greenstein, J. L. 1941, ApJ, 93, 70 [Google Scholar]
 Ishimaru, A. 1978, Wave Propagation and Scattering in Random Media. Volume 1  Single Scattering and Transport Theory (New York: Academic Press), 1 [Google Scholar]
 Krieger, A., & Wolf, S. 2020, A&A, 635, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mie, G. 1908, Ann. Phys., 330, 377 [Google Scholar]
 Min, M., Dullemond, C. P., Dominik, C., de Koter, A., & Hovenier, J. W. 2009, A&A, 497, 155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ober, F., Wolf, S., Uribe, A. L., & Klahr, H. H. 2015, A&A, 579, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reissl, S., Wolf, S., & Brauer, R. 2016, A&A, 593, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Roberge, W. G. 1983, ApJ, 275, 292 [CrossRef] [Google Scholar]
 Robitaille, T. P. 2010, A&A, 520, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Solution to the onedimensional slab problem
In this section the problem of high optical depths of Monte Carlo radiative transfer (MCRT) simulations is studied. To this end, we consider the onedimensional slab problem, derive its solutions in terms of eigenstates of the formulated problem, and from the solution draw conclusions about the underlying difficulties a MCRT simulation encounters.
A.1. Onedimensional slab problem
We consider a onedimensional connected region of total optical depth τ_{tot} > 0 embedded in a vacuum. Let v(x) be the distribution of photons, where x ∈ ℝ is the (extinction) optical depth coordinate. Furthermore, we assume that v(x) may only be nonzero for x ∈ [−b,b], where b = τ_{tot}/2 defines the border of the region. In an iterative process, photons are released isotropically from their corresponding positions and experience their succeeding events of interaction at new positions. A portion of these photons, however, leave the region through one of its borders and do not interact any longer. The randomly traversed individual distance Δx that a photon travels follows a distribution p(Δx) given by
At the location of interaction a photon package either scatters and continues its path through the medium or is absorbed. The probability for scattering and absorption is given by A and 1 − A, respectively, where A is the albedo. The albedo is defined as the fraction between the scattering cross section C_{sca} and the extinction cross section C_{ext}. The expected resulting distribution of interacting photons can then be written as follows:
In the following we derive eigenstates of this operation and formula to expand v(x) in terms of these eigenstates.
A.2. Derivation of eigenstates
Eigenstates v(x) of the operation in Eq. (A.2) with their corresponding eigenvalues λ satisfy the equation
By writing the exponential function in Eq. (A.3) as
where H(x) is the Heaviside step function, and taking the second derivative of Eq. (A.3), we obtain the following secondorder differential equation:
Here v″ denotes the second derivative of v with respect to x. For 0 < λ < 1, the solutions for this common differential equation are given by the antisymmetric function sin(σx) and the symmetric function cos(σx), where
However, when inserting these solutions into Eq. (A.3) we find conditions for σ that have to be satisfied for these eigenstates. The condition for the symmetric and the antisymmetric eigenstates are
where the asterisk indicates that is belongs to the antisymmetric solution. We define the sets of positive solutions for σ to these conditions as Σ and Σ^{*}, respectively the symmetric and antisymmetric solutions, which are thus given by
and
The values of σ_{i} ∈ Σ and increase for increasing i ∈ ℕ_{> 0}, and the sets Σ and Σ^{*} are disjoint. Moreover, we define
Let ℱ be the set of eigenstates
We define the inner product ⟨⋅ ,⋅⟩ as
with for all f ∈ ℱ, where λ_{f} is the eigenvalue corresponding to f obtained by Eq. (A.11). Using these definitions we find that the set ℱ forms an orthonormal system of eigenstates:
Next, we expand v(x) in terms of the eigenstates ℱ as
where the coefficients are given by
A.3. Evolution of eigenstates
Let v_{0}(x) be the initial nonzero distribution of photons. After launch, the resulting distribution of interacting photons is partly absorbed, while the remaining part continues its path through the medium. The distribution of launching photons after the nth event of scattering satisfies
Let a_{k} and b_{k} with k ∈ ℕ_{> 0} be the expansion coefficients of v_{0}. Subsequently, the distribution v_{n} can be written as
where λ_{k} and are the eigenvalues corresponding to σ_{k} and , respectively, which are linked by Eq. (A.6). Given that all corresponding eigenvalues satisfy 0 < λ < 1, we find that the absolute value of the amplitude of every mode of the expansion decreases exponentially with increasing number n and vanishes only in the limit n → ∞. Additionally, all eigenvalues are pairwise different; in particular, they satisfy the relation for all k ∈ ℕ_{> 0}. Thus, any nonzero distribution v(x) has a dominant (i.e., highest) eigenvalue λ_{max}, that either corresponds to a dominant σ_{min} ∈ Σ or σ_{min} ∈ Σ^{*}. Therefore, as n → ∞, the distribution asymptotically converges:
Furthermore, in the case of a nonzero nonnegative initial distribution, the dominant eigenvalue is λ_{1}. Therefore, the shape of this photon distribution converges to the shape of the dominant eigenstate cos(σ_{1}x), independent of the specifics of the initital distribution, which is positive over the full domain of the slab^{2}.
A.4. Observer
In this subsection we derive the relevant equations that let us study the problem of high optical depth and help us to find tools to overcome this problem. To do this we place an observer at x → ∞ who detects photons that leave the optically thick region through its border at x = b. Given a distribution of launching photons v_{0}(x), the measurement of observed photons during a single step is given by
Comparing Eq. (A.20) with (A.2) we find
meaning that after a single step of photon transfer the measurement of the observer equals the photon density at the border of the region that is facing the observer, prior to its interaction with the medium. This equation holds for arbitrary distributions v(x). Thus, after the nth event of interaction and assuming A ≠ 0, the measurement I^{(n)} can be written as
Using the expansion coefficients a_{k} and b_{k} of v_{0}(x), we find the general expression of the measurement to be given by
For Eq. (A.24) we made use of Eqs. (A.35) and (A.38). The measurement of the observer increases with the number of considered interactions to a partial sum I^{( ≤ n)} given by
Since Aλ_{k} < 1 and , we can take the limit of the partial sums and arrive at the total measured quantity I^{tot}:
A.5. Application to a MCRT simulation
In this section we apply the results from Appendix A.4 to the case of a MCRT simulation, in particular, we derive a solution for the transmitted energy for the case of an onedimensional slab.
During a MCRT simulation isolated photon packages originate from different points in the model space and travel through it. The photon distribution thus corresponds to a delta distribution that launches at a position x = a:
Here E is the energy carried by the photon package. It then generally interacts multiple times before leaving its cell. As a consequence, on average, the energy of the photon package spreads over the entire region according to
Depending on the initial position a ∈ [−b,b], the energy that the observer detects can be calculated using Eq. (A.27). Assuming an albedo of A = 1 and using the relations from Appendix A.6, the total observed intensity is given by
which in the case of a = 0 equals exactly E/2; in the case a = b it reaches its maximum value, which asymptotically converges to E as b → ∞; and finally in the case of a = −b it reaches its minimum value of E/(2b + 2), which vanishes as b → ∞.
Transmission through a onedimensional slab. A slab, as described in Appendix A.1, that is illuminated from the left (i.e., from x → −∞) can be described by the initial distribution v_{0}(x) = 2δ(x + b). As a result, the energy that enters the slab is normalized to 1. Using Eqs. (A.16) and (A.26) we arrive at the contribution I^{(n)} to the total transmitted given by
We note that I^{(n)} can also be calculated by the ndimensional integral
with τ_{0} ≐ 0 and τ_{n + 1} ≐ τ_{tot}. The sum in Eq. (A.31) tends to converge slowly for low scattering orders n, which is a consequence of the slow decrease in the eigenvalues λ_{k} with k. The calculation of I^{(n)} based on its integral form, however, is simple in this case, while becoming impractical for increasing scattering order. Therefore, the integral form can be solved analytically up to a certain n and used a measure for convergence of the sum in Eq. (A.31).
The overall observed energy of photons is the sum of all contribution I^{(n)}. The average number of scattering events a photon package has undergone before its detection (N^{mean}) depends on the optical depth and is given by
It is obvious that N^{mean} increases with the albedo, and thus in the case of A = 1 reaches its highest value (). Using equations from Appendix A.6, it is given by
For 0 < b ≪ 1 we find that and for high optical depths (i.e., b ≫ 1, ). Thus, in order to simulate high optical depths properly, MCRT simulations are required to simulate photon packages of increasingly high scattering order.
A.6. Useful relations
When using the results from Appendix A.2 and expanding v(x) ∈ {1,x,x^{2},x^{3}} in terms of the eigenstates in ℱ, we can derive multiple relations, some of which are listed in the following. We note that specific variable names (i.e., σ,λ) marked with an asterisk correspond to antisymmetric eigenstates, while those without the asterisk correspond to symmetric eigenstates.
Appendix B: Method descriptions and transmission spectra
Fig. B.1.
Results for the transmitted intensity I_{trans} through an infinite planeparallel slab with an albedo of A = 0.9 and a transverse optical depth of τ_{max} = 75, where μ = cos(θ) is the cosine of the penetration angle θ. The SePS method (solid line) is compared against a result from a nonprobabilistic numerical method (dotted line). For details, see Sect. 4.1. 
Fig. B.2.
Results for the transmitted intensity I_{trans} through an infinite planeparallel slab with an albedo of A = 0.9 and for transverse optical depths τ_{max} ∈ {75,150,300,600}, where μ = cos(θ) is the cosine of the penetration angle θ. The SePS method (solid lines) is compared against results from a nonprobabilistic numerical method (dotted lines). For details, see Sect. 4.1. 
Overview of the MCRT methods used in this paper.
All Tables
All Figures
Fig. 1.
Sketch of the infinite planeparallel slab setup. Radiation enters the slab of transverse optical depth τ = τ_{max} through the top layer. The observer detects radiation that is transmitted through the slab and is located below the bottom layer of the slab. The initial direction of a photon package is described by the penetration angle θ_{0}. If the photon package scatters, its direction is chosen randomly according to its corresponding PDF and the new direction is described by the angle θ. 

In the text 
Fig. 2.
Detection probability for transmitted photon packages of scattering order n through an onedimensional slab. Plotted are the results for five different transverse optical depths τ_{max} ∈ {5,10,20,40,80} under the assumptions A = 1 and using isotropic scattering. Vertical dashed lines mark the mean scattering order of the corresponding optical depths. For details, see Appendix A.5. 

In the text 
Fig. 3.
Sketch of a radiating sphere at radial position τ_{r} embedded in a larger sphere with radius τ^{sph}. Radiation originating on the surface of the embedded sphere penetrates the surface of the larger sphere about the penetration angle θ_{s}, and thus contributes to its emission spectrum. For details, see Sect. 3.2.1. 

In the text 
Fig. 4.
Total transmitted peeloff contribution δL_{tot} of radiating embedded spheres or points, inside a larger sphere of size = 10^{1.5}. The orange and the blue curve show the contribution when using the extended and the basic peeloff method, respectively. The three gray vertical lines show the transition radii τ_{r}, where a differently sized sphere can be embedded. In particular, the left gray line gives the transition from using a sphere of size to a sphere of size . The middle line separates the τ^{sph} region from the region where the basic peeloff is used. Close to the rim of the sphere (i.e., to the right of the last gray line for radii τ_{r} ≲ ) none of the precalculated spheres fit and the basic peeloff method is applied. 

In the text 
Fig. 5.
Precalculated emission spectra for spheres of optical depth radius for k ∈ {1,2,3,4,5}. The internal source of luminosity is neutral and isotropically radiating. It is placed in the center of the sphere and propagates through the sphere until it leaves it through its rim. The intensity I_{SePS} leaves the surface of the sphere and depends on the cosine of the penetration angle μ_{s} = cos(θ_{s}). The intensity is highest in the direction radially outward and overall decreases toward a smaller μ_{s}. 

In the text 
Fig. 6.
Ratio of emission spectra I_{SePS} to I_{SPS}, which use the SePS method and the SPS method, respectively, during the precalculation. The SPS method reaches increasingly lower intensities for increasing sphere sizes (i.e., for higher values of k). Shown in brackets is the total loss of emitted luminosity when using the SPS method compared to the SePS method. For details, see Sect. 3.2.1. 

In the text 
Fig. 7.
Emitted luminosity per patch size dL^{p, det}/dμ_{p} as a function of the position μ_{p} = cos(θ_{p}) on the surface of the sphere. The area beneath the curve corresponds to the total emitted luminosity of the sphere. Values of μ_{p} increase when nearing a faceon view of the detector; resulting in a rising extended peeloff contribution. 

In the text 
Fig. 8.
Results for the infinite planeparallel slab problem with transverse optical depth τ_{max} = 75. Different advanced MCRT methods (solid lines) are compared to results from a nonprobabilistic numerical method (dotted line). 

In the text 
Fig. 9.
Results for the infinite planeparallel slab problem for transverse optical depths τ_{max} ∈ {75,150,300}. Different advanced MCRT methods (solid lines) are compared to results from a nonprobabilistic numerical method (dotted lines). The forward label indicates a variant of the SePS method that differs only in its application of the Stretch method; it is applied if and only if the photon package scatters toward the detector. 

In the text 
Fig. B.1.
Results for the transmitted intensity I_{trans} through an infinite planeparallel slab with an albedo of A = 0.9 and a transverse optical depth of τ_{max} = 75, where μ = cos(θ) is the cosine of the penetration angle θ. The SePS method (solid line) is compared against a result from a nonprobabilistic numerical method (dotted line). For details, see Sect. 4.1. 

In the text 
Fig. B.2.
Results for the transmitted intensity I_{trans} through an infinite planeparallel slab with an albedo of A = 0.9 and for transverse optical depths τ_{max} ∈ {75,150,300,600}, where μ = cos(θ) is the cosine of the penetration angle θ. The SePS method (solid lines) is compared against results from a nonprobabilistic numerical method (dotted lines). For details, see Sect. 4.1. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.