© The Authors 2023

1 Introduction

Radiative transfer (RT) simulations provide a crucial tool in astrophysics across various areas of application. These simulations allow us to predict and analyze the wavelength-dependent appearance of astrophysical objects, linking it to the underlying physical properties of these systems (e.g., protoplanetary disks, exoplanet atmospheres, filaments). To enable the simulation of complex environments, modern RT simulations often apply the Monte Carlo (MC) approach, which relies on the probabilistic simulation of so-called photon packages (PPs) and their randomly generated paths through a model space. Within the last two decades, various MCRT codes have been developed, including MC3D (Wolf et al. 1999; Wolf et al 2003), MCFOST (Pinte et al. 2006, Pinte et al. 2022), MCMax (Min et al. 2009), RADMC-3D (Dullemond et al. 2012), Mol3D (Ober et al. 2015), POLARIS (Reissl et al. 2016), and SKIRT 9 (Camps & Baes 2020). These codes have been optimized with regard to their computation time by using such methods as a locally divergence free continuous absorption (Lucy 1999), immediate reemission scheme according to a temperature corrected emission spectrum (Bjorkman & Wood 2001), partial diffusion approximation (Min et al. 2009), modified random walk (Min et al. 2009; Robitaille 2010), biasing techniques (Baes et al. 2016), precalculated sphere spectra (Krieger & Wolf 2020), and an extended peel-off method (Krieger & Wolf 2021).

Nonetheless, the analyses adopting these methods may remain limited due to a lack of required computational resources, in particular, the computation time. This is especially the case for the simulations of the RT in systems of high optical depth, for which we have to consider high numbers of simulated interactions a PP undergoes before leaving the model space and, consequently, the computation time rises. For such systems, it has been found that an overly low number of simulated PPs or an insufficient number of simulated interactions may even lead to the so-called scattering order problem (Krieger & Wolf 2021), resulting in severely underestimated flux values (e.g., Baes et al. 2016; Gordon et al. 2017; Camps & Baes 2018; Krieger & Wolf 2020,Krieger & Wolf 2021). We note that the underestimation is a consequence of the needed restriction of simulated PPs and a consequence of very high MC noise combined with non-Gaussian statistics. In that regard, the scattering order problem generally states that the calculation of a reliable flux estimate requires a proper representation of simulated scattering orders, with the latter being the number of interactions a PP has undergone prior to its detection. However, the corresponding scattering order distribution has been shown to widen and shift toward larger scattering orders as the transverse optical depth of a system increases, quickly leading to infeasible demands for MCRT simulations and computational limits that can be difficult to overcome.

Additionally, the problem of unreliable flux estimates has also been reported for systems with optical depths between 10 and 30 (Gordon et al. 2017). As a result, it can even affect simulations of embedded radiation sources, such as young accreting planets (Krieger & Wolf 2022). Moreover, three-dimensional (3D) simulations of fainter sources can be extremely challenging, if the number of simulated sources is high and the feasible number of simulated PPs per source is small and becomes a limiting factor. In particular, the simulation of flux and polarization maps of protoplanetary disks resulting from self-scattering (i.e., dust grains thermally emitting photons that scatter off other dust grains) is a computationally extremely demanding task. However, its effect can be crucial in terms of the inference of disk properties based on observations (e.g., Brunngräber & Wolf 2020, Brunngräber & Wolf 2021).

Therefore, it is generally highly desirable to develop the tools and understanding that would allow us to expand the scope and complexity of simulated objects, while maintaining the expected reliability of the results. Motivated by this goal, we explore the difference of two fundamental MCRT frameworks. Even though both frameworks are generally well known, we find that the most popular MCRT codes currently stick to either one of them with no consideration of the other. We show that these frameworks, however, may differ significantly with regard to their required computation speed and quality of flux estimates. To that end, we use a setup composed of an infinite plane-parallel slab, illuminated from one side, as a testbed, and we compare the results of the MCRT simulations of the transmitted intensity performed in both frameworks. In Sect. 2, we briefly summarize the frameworks and introduce the setup. Subsequently, we perform MCRT simulations to compare the performance of the frameworks in Sect. 3. Lastly, we discuss implications regarding the choice of framework depending on the albedo and optical depth of the system in Sect. 4, before presenting our final conclusions.

2 Methods

In this section, we describe two frameworks, which we call the extinction framework (EF) and the scattering framework (SF), which are both widely used in modern MCRT codes when estimating flux values. The underlying procedure of both frameworks is based on the simulation of PPs on randomly generated paths, which follow certain probability density functions. On its path, a PP may interact multiple times with the medium before eventually leaving the model space. The interactions usually include absorption and scattering, both of which are a cause for extinction. As a result, the weight carried by a PP, which corresponds to an intensity or energy, may be reduced prior to its detection by an observer. To properly simulate this process, the optical properties of the interacting medium have to be taken into account. These include the cross-sections for scattering, C_sca, absorption, C_abs, and extinction, C_ext = C_sca + C_abs, as well as the (single scattering) albedo, A = C_sca/C_ext. These quantities additionally determine the optical depth τ along a straight path of length Δl through the medium with τ = C_intρΔl, where C_int is either of the aforementioned cross-sections and ρ the (average) number density along the path. During a MCRT simulation, the path lengths between two consecutive events of interaction are determined on the basis of randomly chosen optical depths, which follow an exponential distribution.

2.1 MCRT frameworks

The considered frameworks differ with regard to the procedure of path determination. In the EF, the path length between two consecutive events of interaction is determined on the basis of extinction optical depths, τ_ext, and the weight of an interacting PP is reduced at each point of interaction by a factor that equals the albedo of the medium. At its core, this procedure first randomly selects a location for the subsequent interaction event and then selects its interaction type. For improved performance, the PP is split at the location of interaction into an absorbed and a scattered part, of which only the latter is traced during the flux determination, carrying a weight that is reduced by a weighting factor that corresponds to the albedo, A. In the SF, on the contrary, the process of absorption is assumed to occur continuously along the PP path, which can be interpreted as leaving out the probabilistic determination of the absorption process. In other words, this framework allows us to “passify” the absorption process. As a consequence, only the scattering process is actively simulated in a probabilistic manner, such that the scattering optical depth is used as a basis to determine the path length between two consecutive scattering events and the weight is reduced passively. We note that MCRT simulations, which contrary to this study trace the absorbed radiation energy, have been shown to significantly benefit from the usage of a continuous absorption procedure (Lucy 1999). However, throughout this study, we only consider the effect of its choice on derived flux estimates.

In particular, the implementation of these frameworks encompasses three aspects. First, there is the determination of a pseudo random number r ∈ [0, 1) according to a uniform distribution. Second, there is the calculation of the (corresponding) optical depth, τ = − ln (1 − r), and, third, the framework-dependent determination of the thus resulting interaction-free length, Δl, and weighting factor, w. In the EF (SF), the latter are given by Δl = τ/ρ/C_ext (Δl = τ/ρ/C_sca) and w = A ($w=e^{-\text{Δ}l\rho C_{\text{abs}}}$). We note that these equations only apply if the PP experiences another interaction at its subsequent location of interaction. If, on the contrary, it leaves the interacting medium after traversing a distance of only Δl’ < Δl, the weighting factor is not applied within the EF and a changed weighting factor of $w={\text{e}}^{-\text{Δ}{l}^{\prime }\rho C_{\text{abs}}}$ is instead used in the SF. This procedure ensures that the PP leaves the medium with a proper weight.

In the limit of infinite simulated PPs, both frameworks yield the same result in terms of estimated flux values. However, their differences may show in the form of: 1) different convergence speeds in terms of the number of required PPs or simulation time; 2) the computation time per PP; and 3) the level of variance of flux estimates (MC noise) as a function of the number of simulated PPs.

We note, that (in principle) there are many more conceivable frameworks that could be used, however, the EF and SF are the most prominent. Moreover, these frameworks can be transformed into each other by introducing the correct constant stretching factor. As a result of this biasing procedure, the weighting factor would be affected, as the path determination and weighting scheme are interdependent. Eventually, whether the photons are removed along the photon path (SF) or only at the location of the scattering event (EF) depends solely on the chosen stretching factor. A brief derivation of the transformation that links both frameworks is presented in Appendix A.

2.2 Setup

To explore the differences of these frameworks, we adopted the setup of Camps & Baes (2018), which is composed of an infinite plane-parallel 3D homogeneous slab with a total transverse (extinction) optical depth of τ_max, which is embedded in a vacuum. This setup is often used as a testbed for MCRT simulations, in which the slab is illuminated by an isotropic (monochromatic) radiation source from one side and the intensity transmitted through the slab is measured, as a function of µ = cos θ, by a simulated observer on the other side of the slab (for details regarding the setup, see Camps & Baes 2018; Krieger & Wolf 2021). Here, the quantity θ describes the penetration angle of the transmitted PPs with θ = 0 corresponding to the direction perpendicular to the surface of the slab. In total, we use M_bin = 41 detector bins, which are linearly sampled regarding the direction µ ∈ [0, 1]. The transmitted intensity of the jth detector bin is then given by: $$I_{\text{trans}}({\mu }_j)={\displaystyle \sum _{i=1}^N\frac{w_{i,j}}{N}}\frac{M}{2{\mu }_j},$$

where N is the number of simulated PPs and w_{i, j} is the total weight that the ith PP contributes to the detector bin, corresponding to the direction µ_j. Unless mentioned otherwise, we assume an albedo of A = 0.5 and simulate scattering events using an isotropic phase function. For this setup, the result for the transmitted intensity can also be calculated on the basis of a non-probabilistic method (see Camps & Baes 2018; Krieger & Wolf 2021), which we use to assess the quality of the obtained MCRT-based results.

Fig. 1 Transmission curves for τmax = 20 using Ntrans = 104 transmitted PPs per MC-based simulation.

3 Results

To compare the performance of the SF with that of the EF (see Sect. 2.1), we performed MCRT simulations to obtain the transmitted intensity for the slab setup (see Sect. 2.2). We assumed an (extinction) optical depth of τ_max = 20 and ran each simulation until a total of N_trans = 10⁴ PPs were successfully transmitted through the slab. Consequently, the required number of simulated PPs N for both frameworks may differ. For each of these simulations, we additionally calculated the corresponding transmission curve that would emerge if the peeloff method (Yusef-Zadeh et al. 1984) was applied, which is a commonly used MCRT method, that has shown to significantly boost the performance of MCRT simulations (e.g., Krieger & Wolf 2021). The obtained transmission intensity is shown in Fig. 1, together with the result of the non-probabilistic reference solution. We find, that MCRT simulations that use the peel-off method (labeled “peel-off”) generally lead to better results compared to simulations that use the basic approach (labeled “original”), which are prone to higher levels of noise and show a stronger underestimation of flux values, typical of simulations in the regime of high optical depths (e.g., Camps & Baes 2018; Krieger & Wolf 2021). Furthermore, we find that the SF clearly outperforms the EF for this task, as the estimated transmission is much closer to the reference solution in terms of shape and magnitude. It is important to note that these differences occur despite the fact that both simulations were performed for a fixed number of transmitted PPs. In fact, the EF required N = 157 738 simulated PPs of which only N_trans = 10⁴ were transmitted, which is ~83% higher compared to the SF. However, simulation times scale with the number of simulated interactions, N_int, rather than the number of simulated PPs, which in the case of the EF amount to N_int = 6 476 737, exceeding the case of the SF by ~258%. Consequently, using the SF rather than the EF results in a reduction in the computation time of ~72%, while providing better results.

The superiority of the SF over the EF regarding this task can be explained as follows. Since the PP path determination takes place in optical depth coordinates and the albedo satisfies A < 1, the experienced optical depth of the system is smaller within the SF. In this case, the step length of simulated PP path segments is larger in the SF. As a result, PPs are transported with a faster pace through the model space in terms of traversed units of the extinction optical depth per simulated scattering event. This, in turn, leads to a wider spatial spread of simulated scattering locations, which seems to be crucial for the improvement of flux estimates. Additionally, the complexity of the problem is reduced, indicated by the fact that fewer simulated scattering events lead to better flux estimates. In other words, switching to the SF has practically lowered the number of PPs that need to be simulated for a sufficiently well-sampled representation of the distribution of scattering orders (SOs; Krieger & Wolf 2021). Here, the SO of a PP describes its total number of scattering events, n, that have occurred. However, we note in passing that in the limiting case of A → 1, these differences among both frameworks disappear.

3.1 Impact on scattering order distribution

For the detector bin of the direction µ ≈ 0, Fig. 2 shows the SO dependent number of counts C⁽ⁿ⁾ with the corresponding scattering order distribution (SOD) for detected PPs. The SOD is shown in terms of $\text{Δ}{I}_{\text{trans}}^{\left(n\right)}/\text{Δ\hspace{0.17em}}{\log }_{10}\left(n\right)$, where $I_{\text{trans}}^{\left(n\right)}$ is the transmission intensity for SO n. These results stem from simulations, in which the peel-off method was applied. As a consequence, the distributions of the number of counts are the same for all detector bins, since during each simulated interaction one peel-off PP is sent to every detector bin. The SODs, however, vary between different detector bins. The plot clearly shows that within the SF less interactions were simulated, which is reflected by the smaller area below the SF-based count distribution. This is especially the case for lower SOs. Despite that, the corresponding SOD exceeds that of the EF, namely, paths that were generated in the SF are overall of greater relevance for the task of determining transmitted intensity values. Moreover, in this framework, PPs often leave the model space with a lower SO, which is beneficial, since the SOs far from the peak of the SOD (n ≫ n_peak ≈ 6) barely contribute to the transmission intensity. However, these SOs make up the largest portion of the invested computation time. Comparing the SODs, we find that in the EF the transmission intensity is underestimated across all relevant scattering orders n > 0, suggesting that the simulation of highly contributing PP paths is less likely in this framework; furthermore, this even fails for the first simulated scattering event. The difference in performance between both frameworks at τ_max = 20 is striking and can be expected to, firstly, already exist at lower optical depths, and secondly, to increase at higher optical depths.

Fig. 2 Distributions of counts C(n) of PPs (right axis; solid lines) and corresponding transmitted intensities $I_{\text{trans}}^{\left(n\right)}$ (left axis; dashed lines) for scattering order, n, measured by the detector bin for direction, µ ≈ 0 (perpendicular to the surface of the slab).

3.2 Onset of flux underestimation

In the following, we investigate the conditions that lead to the onset of flux underestimation, meaning that we analyze its dependence on the optical depth, albedo, and number of simulated PPs. In contrast to our previous simulations, here we keep the number of PP send out per simulation fixed rather than the number of transmitted PPs. Furthermore, from now on we make use of the peel-off method (unless mentioned otherwise), as it has shown to reliably and significantly boost the quality and performance (e.g., Krieger & Wolf 2021), see for instance Fig. 1.

For this purpose, we analyze the resulting bias and MC noise of estimates of the transmitted intensity depending on different key parameters: the extinction optical depth 2 ≤ τ_ext ≤ 20, the albedo A ∈ {0.1, 0.5, 0.9}, and the number of simulated PPs 10 ≤ N ≤ 10⁶. Additionally, every simulation was repeated 100 times each, allowing for the determination of medians and interquartile ranges (IQRs) to estimate the bias and MC noise of the obtained intensity values, respectively. Figure 3 shows the result of the 100 performed simulations assuming τ_max = 10, A = 0.5, and N = 10⁴. Dashed lines represent the median transmission intensities for both frameworks, the hatched areas mark the regions in the plot within which all determined transmission curves of the corresponding color lie, and the shaded regions account for the central 50% of simulated values for each direction, µ. The displayed median curves are suitable for determining the bias of expectable intensity estimates, which are associated with the usage of a limited number of simulated PPs, when performing a MCRT simulation. Similarly, the IQR, which here represents the central 50% of simulated intensity values, can be used as a measure for the expectable spread of obtained intensity estimates caused by the inherent MC noise of the simulations.

Comparing the results between the SF and EF, we find that the SF outperforms the EF in the task of determining both with regard to the bias, as its median curve $I_{\text{trans}}^{\text{median}}$ is a better estimator for the reference solution (purple solid line) and the MC noise (vertical width of shaded area) is significantly smaller for all directions, µ.

Fig. 3 Transmission curves for τmax = 10. Results of MCRT simulations comprise median transmission curves, each based on 100 individual simulations using N = 104 PPs. The hatched areas mark regions in the plot within which all 100 simulated transmission curves of the corresponding color lie, while the inner shaded regions display the central 50% of simulated values.

3.2.1 Dependence on optical depth

Furthermore, Fig. 4 shows the bias (upper plot) and MC noise (lower plot) as a function of the extinction optical depth of the slab exemplarily for the direction µ = 0.5, which corresponds to radiation leaving the slab with a penetration angle of θ = 60°. In particular, the bias shows in the plot in the form of underestimated intensity values that lead to $I_{\text{trans}}^{\text{median}}/I_{\text{trans}}^{\text{ref}}<1$, where $I_{\text{trans}}^{\text{ref}}$ represents the (linearly interpolated) reference solution. This plot suggests that for τ_ext ≲ 8 both frameworks lead to a rather small underestimation, however, for higher optical depths, the bias associated with the EF lead to a significant underestimation. On the contrary, we find that the bias of the SF, originating from the limitation of PPs, is much smaller. Hence, it results in significantly better intensity estimates, especially for τ_ext > 8, and shifts the onset of the problem toward larger optical depths of τ_ext ≳ 14. The shape of the displayed shaded areas indicates that the spread of expectable intensity values highly depends on the optical depth. In fact, in the case of the EF, the range of the central 50% of simulated intensity values does not even include the reference values for τ_ext ≳ 12. To better illustrate the MC noise as a function of the extinction optical depth of the slab, the shaded region in the lower plot in Fig. 4 shows the central 50% data points for the ratio $I_{\text{trans}}/I_{\text{trans}}^{\text{median}}$ as a function of τ_ext. Here, the vertical width of the shaded region is a measure of the spread of the distribution of intensity estimates. We find, that for increasing optical depth, the noise generally increases in both frameworks. However, this plot clearly suggests that the MC noise associated with the SF is much smaller for all simulated values of τ_ext. As a result, it provides a much more reliable estimate for the transmitted intensity than the EF.

Fig. 4 Bias (upper plot) and MC noise (lower plot) as a function of τmax. These results are based on 100 MCRT simulations each, and comprise median transmission curves (dashed lines) using N = 104 PPs and µ = 0.5. The shaded regions mark the region within which the central 50% of simulated values lie. For details, see Sect. 3.2.1.

3.2.2 Dependence on albedo

We performed the same analysis for a lower (A = 0.1, Fig. C.1) and a higher (A = 0.9, Fig. C.2) albedo and find that the differences between both frameworks increase for lower values of the albedo. This results in an earlier onset of the problem of intensity underestimation when using the EF. In particular, we find that as low values as τ_ext > 6 may lead to errors at the order of ~10% (see Fig. C.1 in the appendix). Considering the amplification of the bias induced underestimation at lower values of µ (e.g., see Fig. 3), even higher deviations can be expected. On the contrary, the SF performs reasonably well for all considered optical depth values up to τ_ext = 20 even with regard to the associated MC noise, which is significantly lower than for the same slab simulations using the EF. Conversely, for A = 0.9 the performance of both frameworks improves both with regard to the bias and the MC noise (see Fig. C.1). For this case, we report only a small benefit from the SF compared to the EF, as can be expected, since both frameworks become more similar as A → 1.

Fig. 5 Bias (upper plot) and MC noise (lower plot) as a function of the number of simulated PPs, N. These results are based on 100 MCRT simulations each and comprise median transmission curves (dashed lines) using N = 104 PPs and µ = 0.5. The shaded regions mark the region within which the central 50% of simulated values lie. For details, see Sect. 3.2.3.

3.2.3 Dependence on number of photon packages

Moreover, the number of send out PPs N was varied for a simulated slab of thickness τ_ext = 10 and its effect on the bias and MC noise of thereby obtained intensity estimates was analyzed. The results of these simulations are displayed in Fig. 5, assuming µ = 0.5. As expected, the bias in the upper plot decreases for higher values of N, which is indicated by the ratio $I_{\text{trans}}^{\text{median}}/I_{\text{trans}}^{\text{ref}}$ approaching a value of 1 and, simultaneously, the noise shown in the lower plot decreases. Additionally, the onset of the problem of flux underestimation therefore shifts toward smaller optical depths as N decreases. A comparison of the results of both frameworks clearly suggests that the number of simulated PPs that is required to achieve a low bias is much smaller for the SF than for the EF. In particular, the number of simulated PP that is needed to achieve a median intensity of ~80% of the reference solution is about 10–100 times higher for the EF than for the SF in this example. Considering the fact that the SF furthermore leads to a smaller number of simulated interactions per PP, this corresponds to a significant potential for saving computation time, already at optical depths as small as τ_ext = 10, which can be expected to further increase in the regime of higher optical depths.

4 Discussion and conclusions

In the previous section, we compare the results of MCRT simulations that use different frameworks, the EF or the SF, for determining the transmitted intensity through an infinite plane-parallel slab that is illuminated from one side. Our findings clearly suggest that the SF is generally better suited for this task than the EF, as the shape and magnitude of those obtained transmission curves are a better match to the non-probabilistic solution, while at the same time requiring a significantly smaller computation time. The benefits not only entail a smaller bias associated with the usage of the SF, but also greater reliability of transmitted flux estimates due to a smaller MC noise. The observed superiority of the SF may stem from a reduced problem complexity, which seems to be determined in this case by the scattering optical depth rather than the extinction optical depth, with the former being smaller by a factor of A.

Generally, the albedo plays a crucial role in determining the complexity of a MCRT problem. For a source that is deeply embedded inside a region of high optical depth, for instance, the number of interaction events scales approximately quadratically with the optical depth (e.g. Krieger & Wolf 2020,Krieger & Wolf 2021). In the SF, the computation time can therefore be expected to also scale with ~A² when simulating deeply embedded sources of radiation (e.g., protostars, protoplanets, visous heating in a protoplanetary disk, etc.). Considering the results for the SODs shown in Fig. 2, the improved quality manifests itself for all n > 0, suggesting, that the SF is promoting highly contributing paths at a much higher rate than the EF. Moreover, it is known that for higher values of A; firstly, the SOD shifts toward larger scattering orders and, secondly, broader SO intervals have to be simulated to ensure a reliable transmitted intensity estimation Krieger & Wolf (2021). The path determination in the SF, interestingly, follows exactly these trends, while the EF path determination utilizes a procedure that is completely independent of the albedo. Therefore, it can be expected that the SF outperforms the EF for all values of A, which is supported by our results (presented in Sect. 3).

Our analysis of the onset of flux underestimation (Sect. 3.2) has shown that this issue can manifest itself already at optical depths as small as τ_ext ≈ 6 or lower, depending on the number of feasibly simulated PPs per source. Such an early onset of this problem may very well explain the high level of MC noise present in MCRT simulations of embedded planets (Krieger & Wolf 2022) or of the process of self-scattering (e.g., Brunngräber & Wolf 2020), which were both performed within the EF. Meaning, a shift to the SF may significantly benefit such simulations in the form of a reduced simulation times and less MC noise.

To evaluate the performance of MCRT simulations in the regime of high optical depths, we additionally performed simulations for a slab of transverse optical depth τ_max = 75 for A ∈ {0.1, 0.5, 0.9} and determined the transmitted intensity (see Fig. B.1). Here, we analyzed the interplay between the framework, the albedo, and an additional usage of composite-biasing (Baes et al. 2016), a prominent method that affects the path determination. These simulations provided further support for the thus far presented conclusions (for details, see Appendix B). Overall, using the SF resulted in benefits for all tested albedos and optical depths, both independent of the additional usage of the composite-biasing method.

Despite the benefits the SF clearly offers, there are cases, in which the EF may outperform the SF. For instance, the back-scattered part of radiation may be misrepresented in the optically thick case when A ≪ 1. In this case, the corresponding flux values may be underestimated as a consequence of the less densely sampled distribution of scattering locations close to the illuminated side of the slab. Whether this is of concern, however, depends on the particular type of simulation and the specifics of the simulated system. In general, when deciding on the framework, it seems crucial to consider the distribution of simulated interaction points. If, for instance, a dense distribution close to the origin of radiation (in units of optical depth) is desired, the EF may be beneficial. Similarly, if the optical depth is very small, the EF may have benefits due to a denser sampling; however, this could be compensated within the SF by the usage of a method that enforces a minimum number of simulated scattering events per PP. A dedicated study that aims at addressing these topics would certainly be of great interest.

However, apart from the path determination procedure, both frameworks also exhibit striking differences regarding their weight determination schemes, which has implications with respect to their benefits and shortcomings. In general, the contribution of a simulated PP path to a flux estimate equals the product of its probability (density) and weight. Considering that the means by which a deviation from some most contributing path is penalized by these frameworks differs. In particular, the SF penalizes a deviation by a weight associated with the spatial deviation from the most contributing path. The EF, on the other hand, penalizes by attributing a lower probability for traversing large extinction optical depths. In the EF, there is always a residual possibility for a simulated PP to traverse any optical depth without losing weight, while in the SF the weight always decreases as the PP traverses an absorbing medium. Therefore, a PP transmitted within the SF always carries information about its complete previous path, as its total length passively affected its (energy) state due to the application of a corresponding weighting factor. As a result, the final weight may vary significantly stronger within the EF, resulting in an overall higher MC noise of the flux estimates (see Sect. 3) even at low optical depths, and notably in the regime of high optical depths. As a result, choosing the SF over the EF likely results in a variance reduction, which has been supported by our results. In other words, “passifying” the absorption process allows us to reduce the number of simulated scattering events while simultaneously decreasing the bias and inherent MC noise. This raises the question of whether it is possible to passify other probabilistic procedures in MCRT simulations or improve them by utilizing machine learning techniques such as reinforcement learning to generate probability density functions for random variables that depend on the simulated environment.

Based on the results of this study, we can conclude that the SF clearly outperforms the EF in the task of determining transmission intensities using MCRT simulations. We showed that this affects the quality of the results significantly in the form of a reduced bias, smaller MC noise, and an alleviation of the problem of underestimated flux values. Additionally, the computational demand to achieve a low bias within the SF can be multiple orders of magnitude lower, leading to an overall considerably faster convergence speed. Its implementation can therefore greatly benefit all modern MCRT codes that are aimed at determining transmitted flux values.

Acknowledgements

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-2)”.

Appendix A Framework transformation

In the following, we demonstrate how the EF and the SF are linked via the application of a biasing technique. By doing so, we furthermore show, that the discrete nature of removing PPs at the location of scattering in the EF thereby transforms into the continuous absorption scheme that is applied in the SF. In the EF, the weighting factor is w = A and the traversed extinction optical depth is selected by the determination of a pseudo random number r ∈ [0, 1) according to a uniform distribution, followed by its calculation: $${\tau }_{\text{ext}}=-\ln (1-r).$$(A.1)

Therefore, the probability density function for the random variable, τ_ext, is given by dr/dτ_ext = exp(–τ_ext). Next, we assume that a PP is emitted and follows a straight path until it interacts at a distance l ∈ [l*, l* + δl) with the medium, where δl << l* is small. In the EF, this interval corresponds to an interval of random numbers r ∈ [r_i, r_i + δr) and to an optical depth interval τ_ext ∈ [τ_e, τ_e + δτ_e), where τ_e = l*ρ C_ext and δτ_e = δlρ C_ext. The expected portion of its initial energy f_sca, that is scattered within this region, is then given by $$f_{\text{sca}}={\displaystyle {\int }_{r_i}^{r_i+\delta r}\text{d}}rA,$$(A.2) $$={\displaystyle {\int }_{{\tau }_e}^{{\tau }_e+\delta {\tau }_e}\text{d}}{\tau }_{\text{ext}}\left(\frac{\text{d}r}{\text{d}{\tau }_{\text{ext}}}\right)A,$$(A.3) $$={\displaystyle {\int }_{{\tau }_e}^{{\tau }_e+\delta {\tau }_e}\text{d}}{\tau }_{\text{ext}}{e}^{-{\tau }_{\text{ext}}}A,$$(A.4) $$=A\delta {\tau }_e{e}^{-{\tau }_e}.$$(A.5)

Similarly, in the SF, the interval corresponds to an optical depth interval τ_sca ∈ [τ_s, τ_s + δτ_s), where τ_s = l*ρ C_sca = A τ_e and δτ_s = δlρ C_sca = Aδτ_e. In this framework, the scattered portion of energy thus equals: $$f_{\text{sca}}={\displaystyle {\int }_{{\tau }_s}^{{\tau }_s+\delta {\tau }_s}\text{d}}{\tau }_{\text{sca}}{e}^{-{\tau }_{\text{sca}}}{e}^{-{\tau }_{\text{abs}}},$$(A.6)

where $e^{-{\tau }_{\text{abs}}}$ is the weighting factor accounting for the continuous absorption of energy, which eventually yields the same result as the calculation in the EF.

Next, we assume that the determined PP path is stretched in Eq. (A.1) by a constant factor of α, namely, τ_ext = −α In (1 – r), such that: $$\left(\frac{\text{d}r}{\text{d}{\tau }_{\text{ext}}}\right)=\frac{1}{\alpha }{e}^{-{\tau }_{\text{ext}}/\alpha }.$$(A.7)

We can then use this result and re-write Eq. (A.4) as follows: $${\displaystyle {\int }_{{\tau }_e}^{{\tau }_e+\delta {\tau }_e}\text{d}}{\tau }_{\text{ext}}{e}^{-{\tau }_{\text{ext}}}A,$$(A.8) $$={\displaystyle {\int }_{{\tau }_e}^{{\tau }_e+\delta {\tau }_e}\mathrm{d}}{\tau }_{\text{ext}}\frac{1}{\alpha }{e}^{-{\tau }_{\text{ext}}/\alpha }\left(\alpha e^{-{\tau }_{\text{ext}}+{\tau }_{\text{ext}}/\alpha }A\right),$$(A.9)

where the term in front of the bracket corresponds to the probability density function for stretched PP paths and the term inside the bracket is the resulting corresponding weighting factor. When switching from the EF to the SF, the stretching factor is α = 1/A, which, when plugging it into the previous equation, yields $${\displaystyle {\int }_{{\tau }_e}^{{\tau }_e+\delta {\tau }_e}\mathrm{d}}{\tau }_{\text{ext}}{e}^{-{\tau }_{\text{ext}}}$$(A.10) $$={\displaystyle {\int }_{{\tau }_e}^{{\tau }_e+\delta {\tau }_e}\mathrm{d}}{\tau }_{\text{ext}}A{e}^{-{\tau }_{\text{ext}}A}\left(e^{-{\tau }_{\text{ext}}+{\tau }_{\text{ext}}A}\right)$$(A.11) $$={\displaystyle {\int }_{{\tau }_s}^{{\tau }_s+\delta {\tau }_s}\mathrm{d}}{\tau }_{\text{sca}}{e}^{-{\tau }_{\text{sca}}}{e}^{-{\tau }_{\text{abs}}}.$$(A.12)

Here, we used a coordinate transformation and substituted τ_ext with τ_sca/A to arrive at Eq. (A.12). As can be seen, this is the same expression as used in the SF (see Eq. (A.6)). Meaning, by selecting τ_ext values and stretching them to match τ_sca, it is necessary to apply the weighting factor $w=e^{-{\tau }_{\text{abs}}}$ which can be interpreted as the continuous absorption of photons along the path. As a result, the EF, in which no stretching is applied (i.e., α = 1), is the only case in which the weighting factor is independent of the selected path length and there is no continuous change of energy along the path.

Appendix B Regime of high optical depths

In order to evaluate the performance of MCRT simulations in the regime of high optical depths, we additionally performed simulations for a slab of transverse optical depth of τ_max = 75 and albedo A = 0.5, and determined the transmitted intensity, using N = 10⁶ send out PPs. However, since the obtained transmission curves exhibit a high level of noise, which is a consequence of the high optical depth, each MC-based simulation was repeated 100 times. Apart from the calculation of a non-probabilistic solution for the problem, we performed three types of MC-based simulations, which all included the peel-off method. One simulation type is based on the SF and does not apply any further MCRT methods (besides the peel-off method) to boost its performance. The other two types additionally use composite-biasing (Baes et al. 2016), which is a MCRT biasing technique used for stretching PP paths and both types differ regarding their utilized frameworks. Stretching paths is particularly interesting, as it has the potential to greatly improve transmitted flux estimates in the regime of high optical depths, (Camps & Baes 2018; Krieger & Wolf 2021) and its performance has not previously been tested in the SF.

The results of these in total 300 MC-based simulations and the non-probabilistic solution are shown in Fig. B.1. We find that simulations that do not use the stretching method fail to properly estimate the transmission curve for all penetration directions, μ. Therefore, using the SF does not by itself suffice to solve the problem of high optical depths, and further techniques are needed to efficiently sample highly contributing paths. When combining the stretching method with either of the two frameworks, the results significantly improve and both types of simulations seem to perform similarly well. This suggests that the usage of the composite-biasing method has practically alleviated the downsides of the EF enough, such that the performance of both frameworks appears to be equivalent. However, it must be stressed, that this is found to be the case for an albedo of A = 0.5.

Fig. B.1 Transmission for τmax = 75. MC-based results comprise median transmission curves, each based on 100 peel-off-based individual simulations with N = 106 send out PPs per simulation. The hatched areas mark regions in the plot within which all 100 simulated transmission curves of the corresponding color lie. For details, see Sect. B.

In further tests, we found that the performance of these three types of simulations strongly depends on the albedo. In the case of a small albedo (A = 0.1), it seems that switching to the SF suffices to properly estimate the transmission, such that no stretching method was needed. Using the EF without stretching, however, resulted in significantly underestimated intensities. This may be explained by the fact that the scattering optical depth of the system equals only 7.5, which is not particularly high. Therefore, the complexity of the MCRT simulations within the SF is much lower compared to the case of EF-based simulations. For an albedo of A = 0.5, the SF without stretching has shown to improve the quality beyond similar simulations performed in the EF, which also further increases when additionally using the composite-biasing technique. Finally, if the albedo is relatively large (A = 0.9), the SF still performs better than the EF when stretching is not used. However, the reduction of computation time was not particularly high in this case. As expected, though, we find that using the composite-biasing method degrades the quality of intensity estimates, which is likely caused by an insufficient coverage of the underlying SOD, in accordance with the findings of Krieger & Wolf (2021).

In summary, we conclude that using the SF rather than the EF resulted in overall lower computation times and consistently better transmitted intensity estimates. This is independent of the fact of whether the stretching method was applied or not.

Appendix C Supplementary material

Fig. C.1 Similar to Fig. 4, but assuming A = 0.1. For details, see Sect. 3.2.2.

Fig. C.2 Similar to Fig. 4, but assuming A = 0.9. For details, see Sect. 3.2.2.

References

All Figures

	Fig. 1 Transmission curves for τmax = 20 using Ntrans = 104 transmitted PPs per MC-based simulation.
In the text

	Fig. 2 Distributions of counts C(n) of PPs (right axis; solid lines) and corresponding transmitted intensities $I_{\text{trans}}^{\left(n\right)}$ (left axis; dashed lines) for scattering order, n, measured by the detector bin for direction, µ ≈ 0 (perpendicular to the surface of the slab).
In the text

	Fig. 3 Transmission curves for τmax = 10. Results of MCRT simulations comprise median transmission curves, each based on 100 individual simulations using N = 104 PPs. The hatched areas mark regions in the plot within which all 100 simulated transmission curves of the corresponding color lie, while the inner shaded regions display the central 50% of simulated values.
In the text

	Fig. 4 Bias (upper plot) and MC noise (lower plot) as a function of τmax. These results are based on 100 MCRT simulations each, and comprise median transmission curves (dashed lines) using N = 104 PPs and µ = 0.5. The shaded regions mark the region within which the central 50% of simulated values lie. For details, see Sect. 3.2.1.
In the text

	Fig. 5 Bias (upper plot) and MC noise (lower plot) as a function of the number of simulated PPs, N. These results are based on 100 MCRT simulations each and comprise median transmission curves (dashed lines) using N = 104 PPs and µ = 0.5. The shaded regions mark the region within which the central 50% of simulated values lie. For details, see Sect. 3.2.3.
In the text

	Fig. B.1 Transmission for τmax = 75. MC-based results comprise median transmission curves, each based on 100 peel-off-based individual simulations with N = 106 send out PPs per simulation. The hatched areas mark regions in the plot within which all 100 simulated transmission curves of the corresponding color lie. For details, see Sect. B.
In the text

	Fig. C.1 Similar to Fig. 4, but assuming A = 0.1. For details, see Sect. 3.2.2.
In the text

	Fig. C.2 Similar to Fig. 4, but assuming A = 0.9. For details, see Sect. 3.2.2.
In the text