© 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, RADMC-3D (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 three-dimensional 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 plane-parallel 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 peel-off 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 peel-off 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 peel-off method, which may greatly boost the computation speed of a MCRT simulation, that is based on the basic peel-off 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 peel-off method as well as the extended peel-off 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 peel-off method. In particular, we look into its capabilities of performing proper high optical depth simulations with the use of Mie-scattering 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 plane-parallel homogeneous slab setup is described, and a MC based probabilistic numerical implementation and a non-probabilistic 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 plane-parallel slab

To study the performance and capabilities of MCRT simulations at high optical depths, we simulate an infinite plane-parallel 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

$$\begin{array}{c}\hfill p_{{\mu }_0}({\mu }_0)=2{\mu }_0,\mathrm{with}0\le {\mu }_0\le 1,\end{array}$$(1)

where μ₀ = cos(θ₀) is the cosine of the initial penetration angle θ₀; in particular, μ₀ = 0 and μ₀ = 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

$$\begin{array}{c}\hfill p_{\tau }(\tau )=exp(-\tau ),\mathrm{with}\tau \ge 0.\end{array}$$(2)

Fig. 1. Sketch of the infinite plane-parallel 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

$$\begin{array}{c}\hfill p_{\mu }(\mu )=\frac{1}{2},\mathrm{with}-1\le \mu \le 1.\end{array}$$(3)

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 composite-biasing (Baes et al. 2016), which is given by

$$\begin{array}{c}\hfill q_{\tau }(\tau )=\frac{1}{2}[exp(-\tau )+\alpha exp(-\alpha \tau )],\mathrm{with}\alpha =\frac{1}{1+T},\end{array}$$(4)

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

$$\begin{array}{c}\hfill w^{\mathrm{stretch}}=\frac{p_{\tau }(\tau )}{q_{\tau }(\tau )}.\end{array}$$(5)

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 equal-sized bins for the directions μ that cover the interval [0,1]. The jth detector bin of the detector measures the intensity I_j given by

$$\begin{array}{c}\hfill I_j=\frac{{\sum }_{n=1}^N{w}_j^{(n)}}{N}\frac{M}{2{\mu }_j},\end{array}$$(6)

where $w_j^{(n)}$ is the contribution (i.e., weight) of the nth photon package and μ_j the central value of the jth detector bin.

2.2. Non-probabilistic numerical implementation

The conceptually simple problem of determining the transmission of radiation through an infinite plane-parallel 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 non-probabilistic 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

$$\begin{array}{c}\hfill \mu \frac{\partial I}{\partial \tau }(\tau ,\mu )=I(\tau ,\mu )-\frac{A}{2}{\int }_{-1}^1\mathrm{d}{\mu }^{\prime }I(\tau ,{\mu }^{\prime })\mathrm{\Psi }(\mu ,{\mu }^{\prime }),\end{array}$$(7)

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 non-probabilistic 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 plane-parallel 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 one-dimensional slab in Appendix A as a new perspective may have the potential to lead to new insights.

Detected photon packages traveled through the one-dimensional slab, scattered in general several times, and eventually left the slab toward the detector. The probability p⁽ⁿ⁾ that a photon package undergoes n scattering events before its detection is given by

$$\begin{array}{c}\hfill p^{(n)}=\frac{I^{(n)}}{I^{\mathrm{tot}}},\end{array}$$(8)

where I⁽ⁿ⁾ 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⁽ⁿ⁾ for the case where A = 1 (i.e., for an absorption-free 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⁹, 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⁽ⁿ⁾ 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⁽ⁿ⁾ with Aⁿ 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 one-dimensional 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 composite-biasing 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: $N_{\text{int}}~{\tau }_{\text{sph}}^2$ (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)²), 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 $\mathrm{\Delta }T\sim f_x^{-1}$ meaning that it is reduced by the factor f_x. As a result the composite-biasing 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 peel-off 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 peel-off 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 peel-off photon package, at every point of interaction of the original photon package. The peel-off photon package is then sent immediately (i.e., interaction-free) 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 peel-off photon package $w_j^{\text{po}}$ which is detected by the jth detector bin is given by

$$\begin{array}{c}\hfill w_j^{\mathrm{po}}=\frac{w_{\mathrm{org}}}{2M}exp(-\frac{{\tau }_{\mathrm{tv}}}{{\mu }_j}),\end{array}$$(9)

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 peel-off method is referred to as the basic peel-off method. In Sect. 4 we present results for a comparison between the S&S method and the Split & Peel-off & Stretch (SPS) method.

3.2. Extended peel-off

According to the findings in Sect. 2.3 it is clear that the basic peel-off 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 peel-off method. Later we combine this method with the S&S method and compare this Split & extended Peel-off & Stretch (SePS) method against the SPS method. Instead of generating a single peel-off photon package at a point of interaction, this method creates several peel-off 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 peel-off 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 peel-off 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 peel-off photon package contains information of all orders n_po ∈ {n,n+1,n+2,…}, unlike the basic peel-off, 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 peel-off 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 peel-off method, on the contrary, contain contributions of arbitrarily high scattering orders already from the first application of the method. As a result, the extended peel-off 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 peel-off method. The MRW is similar to the extended peel-off 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 peel-off 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 Mie-scattering (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 peel-off 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 peel-off method.

3.2.1. Precalculation: Emission spectra

We precalculate spectra of spheres with optical depth radii ${\tau }^{\text{sph}}\in \left\{{\tau }_k^{\text{sph}}\mid k\in \left\{1,2,3,4,5\right\}\wedge {\tau }_k^{\text{sph}}={10}^{k/2}\right\}$. Transverse optical depths below ${\tau }^{\text{sph}}={\tau }_1^{\text{sph}}$ are well represented by basic peel-off. 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 peel-off 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 peel-off 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⁸) 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 $\mathrm{\Delta }{\tau }_{\text{min}}\ge {\tau }_1^{\text{sph}}$. During the precalculation this condition becomes ${\tau }^{\text{sph}}-{\tau }_r\ge {\tau }_1^{\text{sph}}$. 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 peel-off 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 peel-off 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 peel-off emission of the photon package in the direction ν contributes to the larger sphere’s μ_s-bin of direction

$$\begin{array}{c}\hfill {\mu }_{\mathrm{s}}(\nu )=\sqrt{1-\frac{{\tau }_r^2}{{{\tau }^{\mathrm{sph}}}^2}(1-{\nu }^2)}.\end{array}$$(10)

This implies that a scattering event at radius τ_r results in a non-zero contribution for all ${\mu }_{\mathrm{s}}\ge \sqrt{1-\frac{{\tau }_r^2}{{{\tau }^{\text{sph}}}^2}}$. 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

$$\begin{array}{c}\hfill J(\nu )=\frac{1}{2}{\displaystyle \int \mathrm{d}\nu }exp(-x(\nu )),\end{array}$$(11)

where $x(\nu )=-{\tau }_r\nu +\sqrt{{{\tau }^{\text{sph}}}^2-{\tau }_r^2(1-{\nu }^2)}$ 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

$$\begin{array}{c}\hfill J(\nu )=\{\begin{array}{cc}\frac{1}{2}exp(-{\tau }^{\mathrm{sph}})\nu \hfill & \mathrm{if}{\tau }_r=0,\hfill \\ \frac{1}{4{\tau }_r}[exp(-x(\nu ))+({{\tau }^{\mathrm{sph}}}^2-{\tau }_r^2)\mathbf{·}\cdots \hfill & \\ \cdots \mathbf{·}(\frac{exp(-x(\nu ))}{x(\nu )}+\mathrm{Ei}(-x(\nu )))]\hfill & \mathrm{if}0<{\tau }_r<{\tau }^{\mathrm{sph}},\hfill \\ \frac{1}{2}\nu \hfill & \mathrm{if}{\tau }_r={\tau }^{\mathrm{sph}}\wedge \nu >0,\hfill \\ \frac{1}{4{\tau }^{\mathrm{sph}}}[exp(2{\tau }^{\mathrm{sph}}\nu )-1]\hfill & \mathrm{if}{\tau }_r={\tau }^{\mathrm{sph}}\wedge \nu \le 0,\hfill \end{array}\end{array}$$(12)

where Ei(x) is the exponential integral defined as

$$\begin{array}{c}\hfill \mathrm{Ei}(x)=-{\displaystyle {\int }_{-x}^{\infty }}\mathrm{d}s\frac{exp(-s)}{s}.\end{array}$$(13)

The SePS method. If the distance between a scattering position and the surface of the sphere exceeds ${\tau }_1^{\text{sph}}$, 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 multi-dimensional 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−${\tau }_1^{\text{sph}}$,τ^sph], however, no sphere can be embedded, and instead the transmission spectrum is given by the basic peel-off 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 non-zero 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 ${\tau }_2^{\text{sph}}$, ${\tau }_3^{\text{sph}}$, ${\tau }_4^{\text{sph}}$, and ${\tau }_4^{\text{sph}}$ are 97, 161,199, and 221, respectively.

In Fig. 4, we show as an example for τ^sph = ${\tau }_3^{\text{sph}}$ the results for the total contribution δL_tot to the emission spectrum. For this calculation we used 10⁸ 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 peel-off 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 = ${\tau }_3^{\text{sph}}$. Compared to the basic peel-off, the extended peel-off, 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 peel-off. 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 ${\tau }_2^{\text{sph}}$ can be embedded. The second vertical gray line marks the highest radius for the embedment of a sphere of size ${\tau }_1^{\text{sph}}$. Past this radius none of the precalculated spheres fit inside and the basic peel-off method is used.

Fig. 4. Total transmitted peel-off contribution δLtot of radiating embedded spheres or points, inside a larger sphere of size ${\tau }_3^{\text{sph}}$ = 101.5. The orange and the blue curve show the contribution when using the extended and the basic peel-off 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 ${\tau }_2^{\text{sph}}$ to a sphere of size ${\tau }_1^{\text{sph}}$. The middle line separates the τsph region from the region where the basic peel-off is used. Close to the rim of the sphere (i.e., to the right of the last gray line for radii τr ≲ ${\tau }_3^{\text{sph}}$) none of the precalculated spheres fit and the basic peel-off 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, non-interacting 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 ${\tau }_k^{\text{sph}}$ 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 peel-off 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., ${\tau }_4^{\text{sph}}$ = 100) already reaches ∼68% and for k = 5 (i.e., ${\tau }_4^{\text{sph}}$ ≈ 310) even > 99%.

Fig. 5. Precalculated emission spectra for spheres of optical depth radius ${\tau }_k^{\text{sph}}$ 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 ISePS 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 ISePS to ISPS, which use the SePS method and the SPS method, respectively, during the precalculation. The SPS method reaches increasingly lower intensities for increasing sphere sizes ${\tau }_k^{\text{sph}}$ (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 non-constant normal vector. In this regard, the detector of this specific problem is comparably complex.

First we describe a method for propagating extended peel-off 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 ${\tau }_k^{\text{sph}}$ 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 peel-off 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 non-zero density, different points on its surface may result in different detected peel-off 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₁ and p₂ in a patch and their respective normal vectors n₁ and n₂ fulfill n₁ · n₂ ≤ C_n for some chosen constant C_n >  0. The second criterion can be satisfied by imposing the condition |p₁ − p₂| ≤ 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 ${\stackrel{\sim }{\mathit{p}}}_1$ and ${\stackrel{\sim }{\mathit{p}}}_2$ with equal azimuthal position within a patch fulfill

$$\begin{array}{c}\hfill |{\stackrel{\sim }{\mathit{p}}}_1-{\stackrel{\sim }{\mathit{p}}}_2|\le \mathrm{\Delta }{\tau }_{\max }^p,\end{array}$$(14)

where we choose $\mathrm{\Delta }{\tau }_{\text{max}}^p=1$. The condition Eq. (14) results in circular patches at the poles and ring-like 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 ${\tau }_1^{\text{sph}}$, ${\tau }_2^{\text{sph}}$, ${\tau }_3^{\text{sph}}$, ${\tau }_4^{\text{sph}}$, and ${\tau }_4^{\text{sph}}$ the number of patches is N_p = 12, 20, 36, 64, and 112, respectively.

It is possible to precalculate the luminosity $L_j^p$ 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 peel-off method rather than intensities since the weight of a photon package is associated with its carried luminosity. In a MCRT simulation this peel-off luminosity is additionally attenuated by the material on the path between the patch and the detector. The luminosity $L_j^{p,\text{det}}$ that the jth detector bin measures due to the peel-off emission of the patch p, is given by

$$\begin{array}{c}\hfill L_j^{p,\det }=L_j^{p}exp(-\frac{{\tau }_p}{{\mu }_j}),\end{array}$$(15)

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 ${\tau }_3^{\text{sph}}$ and its total patch luminosity $L^{p,\text{det}}={\sum }_j{L}_j^p$ 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 $L^{\text{det}}={\sum }_{j,p}{L}_j^p$.

Fig. 7. Emitted luminosity per patch size dLp, 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 face-on view of the detector; resulting in a rising extended peel-off contribution.

In order to precalculate the projected luminosity of any patch onto the detector we use a MC method in which we generate 10⁷ 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 peel-off 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 peel-off method, we test a simplified version of the extended peel-off method, later referred to as the simple SePS method. In this variant, the whole sphere is represented by a single peel-off photon package, rather than N_p peel-off 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

$$\begin{array}{c}\hfill L_j=\frac{L_{\mathrm{sph}}^{\mathrm{tot}}}{2M},\end{array}$$(16)

where $L_{\mathrm{sph}}^{\mathrm{tot}}$ is the total luminosity emitted by the sphere. The peel-off 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 plane-parallel 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 peel-off 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 non-probabilistic numerical calculation (see Sect. 2.2) and is used for comparison with the results of different MC based methods. A total of 10⁹ photon packages are used in all MCRT simulations. The plot clearly shows that the three simulations based on the peel-off 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 peel-off based simulation sends individual peel-off 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 non-probabilistic solution. Nonetheless, both methods deviate by about one order of magnitude from the non-probabilistic 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 plane-parallel slab problem with transverse optical depth τmax = 75. Different advanced MCRT methods (solid lines) are compared to results from a non-probabilistic 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 peel-off 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¹) 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¹³ 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⁴ 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⁹ 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 bi-directional 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 non-probabilistic solution.

Fig. 9. Results for the infinite plane-parallel slab problem for transverse optical depths τmax ∈ {75,150,300}. Different advanced MCRT methods (solid lines) are compared to results from a non-probabilistic 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 non-probabilistic 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 non-probabilistic 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 plane-parallel 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 non-probabilistic 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 peel-off method and specifically applied it to the infinite plane-parallel 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 peel-off method.

Henyey-Greenstein (HG) scattering (Henyey & Greenstein 1941) is a widely used one-parameter model for anisotropic scattering. By integrating it into the extended peel-off 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 peel-off 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 Mie-scattering, whose outcome of scattering events depends on and affects the polarization state of a photon. It opens up many possibilities for in-depth 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

$$\begin{array}{c}\hfill I_{\mathrm{in}}^2\le Q_{\mathrm{in}}^2+U_{\mathrm{in}}^2+V_{\mathrm{in}}^2,\end{array}$$(17)

and sampling them may be realized for instance by using a grid of cartesian or spherical coordinates. This means that whenever the extended peel-off method is applied the three-dimensional 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 peel-off 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 three-dimensional infinite plane-parallel 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 one-dimensional 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 peel-off method as a promising candidate for solving the scattering order problem. When applied, the basic peel-off method repeatedly generates peel-off photon packages at every scattering point on the path of a representative photon package. Since peel-off photon packages are sent on an interaction-free path toward the detector, the detected set of peel-off 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 peel-off 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 peel-off method. Rather than generating a peel-off 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 peel-off 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 peel-off 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 peel-off 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 peel-off 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 peel-off 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 peel-off approach, and second, that MRW is not capable of simulating polarization and Mie-scattering. In Sect. 4.2, the inclusion of anisotropic scattering or even polarization and Mie-scattering 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 state-of-the-art 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⁹ 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.

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/17-1)”.

References

Appendix A: Solution to the one-dimensional 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 one-dimensional 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. One-dimensional slab problem

We consider a one-dimensional 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 non-zero 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

$$\begin{array}{c}\hfill p(\mathrm{\Delta }x)=\frac{1}{2}exp(-|\mathrm{\Delta }x|).\end{array}$$(A.1)

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 $\stackrel{\sim }{\mathit{v}}(x)$ can then be written as follows:

$$\begin{array}{c}\hfill \stackrel{\sim }{v}(x)=\frac{1}{2}{\displaystyle {\int }_{-b}^b\mathrm{d}yv(y)exp(-|x-y|).}\end{array}$$(A.2)

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

$$\begin{array}{c}\hfill \lambda v(x)=\frac{1}{2}{\displaystyle {\int }_{-b}^b\mathrm{d}yv(y)exp(-|x-y|).}\end{array}$$(A.3)

By writing the exponential function in Eq. (A.3) as

$$\begin{array}{c}\hfill exp(-|x-y|)=\mathrm{H}(x-y)exp(y-x)+\mathrm{H}(y-x)exp(x-y),\end{array}$$(A.4)

where H(x) is the Heaviside step function, and taking the second derivative of Eq. (A.3), we obtain the following second-order differential equation:

$$\begin{array}{c}\hfill \lambda v^{″}+(1-\lambda )v=0.\end{array}$$(A.5)

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 anti-symmetric function sin(σx) and the symmetric function cos(σx), where

$$\begin{array}{c}\hfill {\sigma }^2=\frac{1}{\lambda }-1.\end{array}$$(A.6)

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 anti-symmetric eigenstates are

$$\begin{array}{ccc}& cos\left(\sigma b\right)=\sigma sin\left(\sigma b\right)\hfill & \hfill (\mathrm{symmetric}\mathrm{states}),\end{array}$$(A.7)

$$\begin{array}{ccc}\hfill & sin\left({\sigma }^{\ast }b\right)=-{\sigma }^{\ast }cos\left({\sigma }^{\ast }b\right)\hfill & \hfill (\mathrm{anti}\text{-}\mathrm{symmetric}\mathrm{states}),\end{array}$$(A.8)

where the asterisk indicates that is belongs to the anti-symmetric solution. We define the sets of positive solutions for σ to these conditions as Σ and Σ^*, respectively the symmetric and anti-symmetric solutions, which are thus given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathrm{\Sigma }& =\{{\sigma }_i|i\in {\mathbb{N}}_{>0}\wedge i-1<\frac{{\sigma }_{i}b}{\pi }<i\wedge \cdots \hfill \\ \hfill & \cdots \wedge cos({\sigma }_{i}b)={\sigma }_{i}sin({\sigma }_{i}b)\}\hfill \end{array}\end{array}$$(A.9)

and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathrm{\Sigma }}^{\ast }& =\{{\sigma }_i^{\ast }|i\in {\mathbb{N}}_{>0}\wedge i-\frac{1}{2}<\frac{{\sigma }_i^{\ast }b}{\pi }<i+\frac{1}{2}\wedge \cdots \hfill \\ \hfill & \cdots \wedge sin({\sigma }_i^{\ast }b)=-{\sigma }_i^{\ast }cos({\sigma }_i^{\ast }b)\}.\hfill \end{array}\end{array}$$(A.10)

The values of σ_i ∈ Σ and ${\sigma }_i^{*}\in {\mathrm{\Sigma }}^{*}$ increase for increasing i ∈ ℕ_{> 0}, and the sets Σ and Σ^* are disjoint. Moreover, we define

$$\begin{array}{c}\hfill {\lambda }_i=\frac{1}{1+{{\sigma }_i}^2}\wedge {\lambda }_i^{\ast }=\frac{1}{1+{{\sigma }_i^{\ast }}^2}\forall i\in {\mathbb{N}}_{>0}.\end{array}$$(A.11)

Let ℱ be the set of eigenstates

$$\begin{array}{c}\hfill \mathcal{F}=\{cos({\sigma }_{i}x)|i\in {\mathbb{N}}_{>0}\}\cup \{sin\left({\sigma }_i^{\ast }x\right)|i\in {\mathbb{N}}_{>0}\}.\end{array}$$(A.12)

We define the inner product ⟨⋅ ,⋅⟩ as

$$\begin{array}{c}\hfill \langle f,g\rangle ={\displaystyle {\int }_{-b}^b\mathrm{d}x\overline{f(x)}g(x),}\end{array}$$(A.13)

with $\overline{f(x)}\doteq f(x)/(b+{\lambda }_f)$ 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:

$$\begin{array}{c}\hfill \forall f,g\in \mathcal{F}:\langle f,g\rangle =\{\begin{array}{cc}1\hfill & \mathrm{if}f=g\hfill \\ 0\hfill & \mathrm{else}.\hfill \end{array}\end{array}$$(A.14)

Next, we expand v(x) in terms of the eigenstates ℱ as

$$\begin{array}{c}\hfill v(x)={\displaystyle \sum _{k=1}^{\infty }{a}_{k}cos({\sigma }_{k}x)+b_{k}sin({\sigma }_k^{\ast }x),}\end{array}$$(A.15)

where the coefficients are given by

$$\begin{array}{c}\hfill a_k=\langle cos({\sigma }_{k}x),v(x)\rangle \wedge b_k=\langle sin\left({\sigma }_k^{\ast }x\right),v(x)\rangle .\end{array}$$(A.16)

A.3. Evolution of eigenstates

Let v₀(x) be the initial non-zero distribution of photons. After launch, the resulting distribution of interacting photons ${\stackrel{\sim }{\mathit{v}}}_0(x)$ is partly absorbed, while the remaining part ${\mathit{v}}_1(x)\doteq A{\stackrel{\sim }{\mathit{v}}}_0(x)$ continues its path through the medium. The distribution of launching photons after the nth event of scattering satisfies

$$\begin{array}{c}\hfill v_n(x)=A{\stackrel{\sim }{v}}_{n-1}.\end{array}$$(A.17)

Let a_k and b_k with k ∈ ℕ_{> 0} be the expansion coefficients of v₀. Subsequently, the distribution v_n can be written as

$$\begin{array}{c}\hfill v_n(x)=A^n{\displaystyle \sum _{k=1}^{\infty }{{\lambda }_k}^n{a}_{k}cos({\sigma }_{k}x)+{{\lambda }_k^{\ast }}^n{b}_{k}sin({\sigma }_k^{\ast }x),}\end{array}$$(A.18)

where λ_k and ${\lambda }_k^{*}$ are the eigenvalues corresponding to σ_k and ${\sigma }_k^{*}$, 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 ${\lambda }_k>{\lambda }_k{}^{*}>{\lambda }_{k+1}>{\lambda }_{k+1}{}^{*}$ for all k ∈ ℕ_{> 0}. Thus, any non-zero 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:

$$\begin{array}{c}\hfill v_n(x){\underset{n\to \infty }{\overset{A}{\to }}}^n{\lambda }_{\max }^n·\{\begin{array}{cc}cos\left({\sigma }_{\min }x\right)\hfill & \mathrm{if}{\sigma }_{\min }\in \mathrm{\Sigma }\hfill \\ sin\left({\sigma }_{\min }x\right)\hfill & \mathrm{if}{\sigma }_{\min }\in {\mathrm{\Sigma }}^{\ast }.\hfill \end{array}\end{array}$$(A.19)

Furthermore, in the case of a non-zero non-negative initial distribution, the dominant eigenvalue is λ₁. Therefore, the shape of this photon distribution converges to the shape of the dominant eigenstate cos(σ₁x), independent of the specifics of the initital distribution, which is positive over the full domain of the slab².

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₀(x), the measurement of observed photons during a single step is given by

$$\begin{array}{c}\hfill I(v_0)=\frac{1}{2}{\displaystyle {\int }_{-b}^b\mathrm{d}y{v}_0(y)exp(y-b).}\end{array}$$(A.20)

Comparing Eq. (A.20) with (A.2) we find

$$\begin{array}{c}\hfill I(v_0)={\stackrel{\sim }{v}}_1(b),\end{array}$$(A.21)

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⁽ⁿ⁾ can be written as

$$\begin{array}{c}\hfill I^{(n)}=\frac{v_{n+1}(b)}{A}.\end{array}$$(A.22)

Using the expansion coefficients a_k and b_k of v₀(x), we find the general expression of the measurement to be given by

$$\begin{array}{cc}& I^{(n)}=A^n{\displaystyle \sum _{k=1}^{\infty }{{\lambda }_k}^{n+1}{a}_{k}cos({\sigma }_{k}b)+{{\lambda }_k^{\ast }}^{n+1}{b}_{k}sin({\sigma }_k^{\ast }b)}\hfill \end{array}$$(A.23)

$$\begin{array}{cc}\hfill & =A^n{\displaystyle \sum _{k=1}^{\infty }{(-1)}^{k+1}({{\lambda }_k}^{n+1}\sqrt{1-{\lambda }_k}{a}_k+{{\lambda }_k^{\ast }}^{n+1}\sqrt{1-{\lambda }_k^{\ast }}{b}_k).}\hfill \end{array}$$(A.24)

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

$$\begin{array}{cc}& I^{(\le n)}={\displaystyle \sum _{m=0}^n{I}^{(m)}}\hfill \end{array}$$(A.25)

$$\begin{array}{cc}\hfill & ={\displaystyle \sum _{k=1}^{\infty }{(-1)}^{k+1}({\lambda }_k\frac{1-{\left(A{\lambda }_k\right)}^{n+1}}{1-A{\lambda }_k}\sqrt{1-{\lambda }_k}{a}_k+\cdots }\hfill \\ \hfill & \cdots +{\lambda }_k^{\ast }\frac{1-{\left(A{\lambda }_k^{\ast }\right)}^{n+1}}{1-A{\lambda }_k^{\ast }}\sqrt{1-{\lambda }_k^{\ast }}{b}_k).\hfill \end{array}$$(A.26)

Since |Aλ_k| <  1 and $\text{|}A{\lambda }_k^{*}<1|$, we can take the limit of the partial sums and arrive at the total measured quantity I^tot:

$$\begin{array}{c}\hfill I^{\mathrm{tot}}={\displaystyle \sum _{k=1}^{\infty }{(-1)}^{k+1}(\frac{{\lambda }_k\sqrt{1-{\lambda }_k}}{1-A{\lambda }_k}{a}_k+\frac{{\lambda }_k^{\ast }\sqrt{1-{\lambda }_k^{\ast }}}{1-A{\lambda }_k^{\ast }}{b}_k).}\end{array}$$(A.27)

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 one-dimensional 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:

$$\begin{array}{c}\hfill v_0(x)=E\delta (x-a).\end{array}$$(A.28)

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

$$\begin{array}{c}\hfill v^{\mathrm{tot}}={\displaystyle \sum _{n=0}^{\infty }{v}_n(x).}\end{array}$$(A.29)

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

$$\begin{array}{c}\hfill I^{\mathrm{tot}}=\frac{E}{2}(1+\frac{a}{b+1}),\end{array}$$(A.30)

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 one-dimensional 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₀(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⁽ⁿ⁾ to the total transmitted given by

$$\begin{array}{c}\hfill I^{(n)}\doteq A^n{I}_{A=1}^{(n)}=2A^n{\displaystyle \sum _{k=1}^{\infty }\frac{{\lambda }_k^{n+1}(1-{\lambda }_k)}{b+{\lambda }_k}-\frac{{{\lambda }_k^{\ast }}^{n+1}(1-{\lambda }_k^{\ast })}{b+{\lambda }_k^{\ast }}.}\end{array}$$(A.31)

We note that I⁽ⁿ⁾ can also be calculated by the n-dimensional integral

$$\begin{array}{c}\hfill I^{(n)}={\left(\frac{A}{2}\right)}^n{{\displaystyle \int }}_{{[0,{\tau }_{\mathrm{tot}}]}^n}\mathrm{d}{\tau }_1\cdots \mathrm{d}{\tau }_{n}exp{\displaystyle \{-\sum _{k=0}^n|{\tau }_{k+1}-{\tau }_k|\},}\end{array}$$(A.32)

with τ₀ ≐ 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⁽ⁿ⁾ 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⁽ⁿ⁾. 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

$$\begin{array}{c}\hfill N^{\mathrm{mean}}=\frac{1}{I^{\mathrm{tot}}}{\displaystyle \sum _{n=0}^{\infty }n{I}^{(n)}.}\end{array}$$(A.33)

It is obvious that N^mean increases with the albedo, and thus in the case of A = 1 reaches its highest value ($N_{A=1}^{\mathrm{mean}}$). Using equations from Appendix A.6, it is given by

$$\begin{array}{c}\hfill N_{A=1}^{\mathrm{mean}}=\frac{2b^3+6b^2+3b}{3b+3}.\end{array}$$(A.34)

For 0 <  b ≪ 1 we find that $N_{A=1}^{\mathrm{mean}}\sim b$ and for high optical depths (i.e., b ≫ 1, $N_{A=1}^{\mathrm{mean}}\sim 2b^2/3$). 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²,x³} 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 anti-symmetric eigenstates, while those without the asterisk correspond to symmetric eigenstates.

$$\begin{array}{cc}& cos\left({\sigma }_{n}b\right)={(-1)}^{n+1}\sqrt{1-{\lambda }_n}\hfill \end{array}$$(A.35)

$$\begin{array}{cc}\hfill & sin\left({\sigma }_{n}b\right)={(-1)}^{n+1}\sqrt{{\lambda }_n}\hfill \end{array}$$(A.36)

$$\begin{array}{cc}\hfill & cos\left({{\sigma }_n}^{\ast }b\right)={(-1)}^n\sqrt{{{\lambda }_n}^{\ast }}\hfill \end{array}$$(A.37)

$$\begin{array}{cc}\hfill & sin\left({{\sigma }_n}^{\ast }b\right)={(-1)}^{n+1}\sqrt{1-{{\lambda }_n}^{\ast }}\hfill \end{array}$$(A.38)

$$\begin{array}{cc}\hfill & {\displaystyle \sum _{n=1}^{\infty }\frac{sin\left({\sigma }_{n}b\right)}{{\sigma }_n}\frac{cos\left({\sigma }_{n}a\right)}{b+{\lambda }_n}=\frac{1}{2}\forall a\in [-b,b]}\hfill \end{array}$$(A.39)

$$\begin{array}{cc}\hfill & {\displaystyle \sum _{n=1}^{\infty }\frac{{\lambda }_n}{b+{\lambda }_n}=\frac{1}{2}}\hfill \end{array}$$(A.40)

$$\begin{array}{cc}\hfill & {\displaystyle \sum _{n=1}^{\infty }\frac{{(-1)}^{n+1}}{b+{\lambda }_n}\frac{{\lambda }_n}{\sqrt{1-{\lambda }_n}}=\frac{1}{2}}\hfill \end{array}$$(A.41)

$$\begin{array}{cc}\hfill & {\displaystyle \sum _{n=1}^{\infty }\frac{1}{b+{\lambda }_n}\frac{{\lambda }_n^2}{1-{\lambda }_n}=\frac{b}{2}}\hfill \end{array}$$(A.42)

$$\begin{array}{cc}\hfill & {\displaystyle \sum _{n=1}^{\infty }\frac{cos\left({\sigma }_n^{\ast }b\right)}{{\sigma }_n^{\ast }}\frac{sin\left({\sigma }_n^{\ast }a\right)}{b+{\lambda }_n^{\ast }}=-\frac{1}{2}\frac{a}{b+1}\forall a\in [-b,b]}\hfill \end{array}$$(A.43)

$$\begin{array}{cc}\hfill & {\displaystyle \sum _{n=1}^{\infty }\frac{{\lambda }_n^{\ast }}{b+{\lambda }_n^{\ast }}=\frac{1}{2}\frac{b}{b+1}}\hfill \end{array}$$(A.44)

$$\begin{array}{cc}\hfill & {\displaystyle \sum _{n=1}^{\infty }\frac{sin\left({\sigma }_{n}b\right)}{{\sigma }_n^3}\frac{cos\left({\sigma }_{n}a\right)}{b+{\lambda }_n}=\frac{1}{4}(b^2+2b-a^2)\forall a\in [-b,b]}\hfill \end{array}$$(A.45)

$$\begin{array}{cc}\hfill & {\displaystyle \sum _{n=1}^{\infty }\frac{{(-1)}^{n+1}}{b+{\lambda }_n}\frac{{\lambda }_n^2}{{(1-{\lambda }_n)}^{\frac{3}{2}}}=\frac{b}{4}(b+2)}\hfill \end{array}$$(A.46)

$$\begin{array}{cc}\hfill & {\displaystyle \sum _{n=1}^{\infty }\frac{cos\left({{\sigma }_n}^{\ast }b\right)}{{{\sigma }_n^{\ast }}^3}\frac{sin\left({\sigma }_n^{\ast }a\right)}{b+{\lambda }_n^{\ast }}}\hfill \\ \hfill & =\frac{a^3(1+b)-a{b}^2(b+3)}{12{(1+b)}^2}\forall a\in [-b,b]\hfill \end{array}$$(A.47)

$$\begin{array}{cc}\hfill & {\displaystyle \sum _{n=1}^{\infty }\frac{1}{b+{\lambda }_n^{\ast }}\frac{{{\lambda }_n^{\ast }}^2}{1-{\lambda }_n^{\ast }}=\frac{1}{6}\frac{b^3}{{(1+b)}^2}}\hfill \end{array}$$(A.48)

Appendix B: Method descriptions and transmission spectra

Fig. B.1. Results for the transmitted intensity Itrans through an infinite plane-parallel 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 non-probabilistic numerical method (dotted line). For details, see Sect. 4.1.

Fig. B.2. Results for the transmitted intensity Itrans through an infinite plane-parallel 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 non-probabilistic numerical method (dotted lines). For details, see Sect. 4.1.

All Tables

All Figures

Fig. 1. Sketch of the infinite plane-parallel 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 one-dimensional 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 peel-off contribution δLtot of radiating embedded spheres or points, inside a larger sphere of size ${\tau }_3^{\text{sph}}$ = 101.5. The orange and the blue curve show the contribution when using the extended and the basic peel-off 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 ${\tau }_2^{\text{sph}}$ to a sphere of size ${\tau }_1^{\text{sph}}$. The middle line separates the τsph region from the region where the basic peel-off is used. Close to the rim of the sphere (i.e., to the right of the last gray line for radii τr ≲ ${\tau }_3^{\text{sph}}$) none of the precalculated spheres fit and the basic peel-off method is applied.

In the text

Fig. 5. Precalculated emission spectra for spheres of optical depth radius ${\tau }_k^{\text{sph}}$ 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 ISePS 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 ISePS to ISPS, which use the SePS method and the SPS method, respectively, during the precalculation. The SPS method reaches increasingly lower intensities for increasing sphere sizes ${\tau }_k^{\text{sph}}$ (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 dLp, 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 face-on view of the detector; resulting in a rising extended peel-off contribution.
In the text

	Fig. 8. Results for the infinite plane-parallel slab problem with transverse optical depth τmax = 75. Different advanced MCRT methods (solid lines) are compared to results from a non-probabilistic numerical method (dotted line).
In the text

Fig. 9. Results for the infinite plane-parallel slab problem for transverse optical depths τmax ∈ {75,150,300}. Different advanced MCRT methods (solid lines) are compared to results from a non-probabilistic 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 Itrans through an infinite plane-parallel 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 non-probabilistic numerical method (dotted line). For details, see Sect. 4.1.
In the text

	Fig. B.2. Results for the transmitted intensity Itrans through an infinite plane-parallel 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 non-probabilistic numerical method (dotted lines). For details, see Sect. 4.1.
In the text