© The Authors 2024

1 Introduction

Monte Carlo radiative transfer (MCRT) simulations are a widely used tool that enables the simulation and study of the physical conditions prevalent in various astrophysical systems. Among other applications, they allow us to calculate temperature distributions of these systems on the basis of numerical models, as well as to determine their appearance in terms of wavelength-dependent intensity maps. Due to its universality, this Monte Carlo (MC) based approach even allows for simulations of complex systems, such as filaments, protoplanetary or circumbinary disks, or planetary atmospheres. The basic principle behind MCRT codes such as RADMC-3D (Dullemond et al. 2012), POLARIS (Reissl et al. 2016), or SKIRT 9 (Camps & Baes 2020) is the simulation of photon packages on randomly determined probabilistic paths. These photon packages are emitted from radiation sources and travel through the simulated model space, in which they interact potentially multiple times with the traversed medium before eventually leaving it.

However, this approach faces a fundamental problem when simulating regions of high optical depths, as the required number of interactions and with it the computational demand rises. This ultimately puts a practical limit on the maximum optical depth that can be simulated. This well-known problem can result in noisy temperature distributions and flux values that are severely underestimated (e.g., Baes et al. 2016; Gordon et al. 2017; Camps & Baes 2018; Krieger & Wolf 2020). Although it has already been shown that this issue may already arise for systems with optical depths between 10 and 30, where it has led to significant discrepancies in the outcome of different MCRT codes (Gordon et al. 2017), the severity of it seems to increase strongly with the optical depth of the system, making the regime of high optical depths well-suited for the search for potential solutions.

Attempts have been made to alleviate this issue by the usage of advanced MCRT methods such as the partial diffusion approximation (Min et al. 2009), the modified random walk (Min et al. 2009; Robitaille 2010), biasing techniques (Baes et al. 2016), precalculated sphere spectra (Krieger & Wolf 2020), and the extended peel-off method (Krieger & Wolf 2021). Each of these methods has their respective use cases, nonetheless, the core problem remains unsolved. For this reason, Camps & Baes (2018) carried out a study, focussing on the phenomenon of underestimated flux values in the regime of high optical depths. They analyzed the limit of optical depth values that can be simulated on the basis of the numerical setup of an infinite plane-parallel slab that is illuminated from one side, for which they determined the transmitted intensity using different MCRT methods. According to this study, the most reliable estimate for the transmitted intensity was achieved, when utilizing a combination of different MCRT methods. Their best method, in particular, used a split method, which reduces the weight of a photon package at each point of interaction, but in return enables the enforcement of the interaction to be a scattering event rather than an absorption event. Additionally, a stretching method was applied, which biases the distance between two consecutive points of interaction toward larger optical depths, i.e., it stretches the path of the photon package.

Krieger & Wolf (2021) then performed an in-depth study on the core problems of MCRT simulations in the regime of high optical depths and identified the scattering order problem, which is accompanied by a pathfinding problem. The scattering order problem can be summarized as follows: The total simulated flux is composed of contributions from simulated photon packages of different scattering orders; the latter refers to the number of scattering events a photon package has undergone prior to its simulated detection. A reliable flux estimate, thus, relies on a proper representation of the scattering order distribution (SOD), namely, the distribution of contributions to the total flux with regard to the different scattering orders. At higher optical depths, however, this is becoming increasingly difficult, as the range of relevant scattering orders may broaden and shift toward extremely large values (Krieger & Wolf 2021), resulting in the need for the simulation of a large quantity of photon packages of potentially extremely high scattering orders. This problem was shown to be alleviated with the additional utilization of the peel-off method (Yusef-Zadeh et al. 1984) or the extended peel-off method (Krieger & Wolf 2021); nonetheless, (transverse) optical depth values of τ_max ≳ 150 still resulted in a notable underestimation of transmitted flux values when restricting the simulation time to a feasible extent. This breakdown at high optical depths could, yet again, be attributed to an inadequate representation of the SOD, stressing the need for more efficient MCRT methods. Thus, apart from a potential improvement of the method that is used to determine photon package paths, the question of the smallest number of required simulated scattering events was raised. In MCRT simulations, the total flux strictly increases with the number of considered scattering events per photon package, assuming the simultaneous usage of a peel-off based method; flux values may therefore be severely underestimated, if the simulated photon packages do not undergo a certain minimum number of scattering events.

The concept of such a minimum scattering order is the subject of the present study. In Sect. 2, we derive an approximate SOD for the testbed of a one-dimensional (1D) slab and we explore two cases that shape the kind of difficulty that arises in the regime of high optical depth: the albedo-dominated case and the optical depth-dominated case. Based on that, we present a simple formula in Sect. 3 that estimates the required minimum scattering order and we propose a simple method that allows for its utilization. In Sects. 3.1 and 3.2 the proposed method is applied to assess its benefits based on MCRT simulations of two different systems, a three-dimensional (3D) infinite plane-parallel slab and a protoplanet embedded in an edge-on circumstellar disk, respectively. In Sect. 4, the limits of the method are discussed, along with the implications of our findings regarding the effectiveness of stretching methods and a potential way for their assessment. Finally, our results are summarized in Sect. 5.

2 Derivation

In order to explore the possibility for improving the capabilities of modern MCRT simulations at high optical depths, we studied the setup of a homogeneous 1D slab, which is illuminated from one side by a monochromatic radiation source and we analyzed the SOD of photons, that are represented by photon packages in a MCRT simulation, transmitted through the slab. This study builds upon the results of our previous work (Krieger & Wolf 2021), which we summarize in Sect. 2.1. In Sect. 2.2, we derived an approximate solution for the underlying SOD of the studied setup, which we use to infer various of its properties. In Sect. 2.3, we discuss these findings with regard to the implications regarding the complexity of MCRT simulations, before proposing a simple to implement method that aims at solving the problems at high optical depths in Sect. 3.

2.1 Setup

We adopted the setup of Krieger & Wolf (2021), which is composed of a homogeneous slab with a total (extinction) optical depth of τ_max that is embedded in a vacuum. For any given distribution v₀(x) of photons (or photon packages) inside the slab, where x ∈ ℝ describes the position in terms of optical depth units with x = 0 the center of the slab, the radiative transfer process can be described by an iterative process. During its n-th step, photons are launched at locations x, according to their current distribution v_n(x), and traverse the slab until either leaving it through one of its borders or interacting at a different location within the slab. This results in the distribution v_n+1 (x) of the next step. The randomly selected traversed distances between two consecutive points of interaction follow an exponential distribution, resulting in a spreading of the photon distribution during each step, while the directions are chosen randomly according to an isotropic distribution. The event of interaction itself is simulated by rescaling the whole photon distribution by a factor of A, which corresponds to the single scattering albedo of the slab. The process of spreading can be described in terms of eigenstates (Krieger & Wolf 2021), which themselves are photon distributions that remain unchanged during this process, apart from some eigenstate-specific constant factor λ, namely, the eigenvalue. To distinguish between symmetric and anti-symmetric eigenstates, we marked variables corresponding to the latter with an asterisk. For more details, we refer to the appendix in Krieger & Wolf (2021).

The total intensity of photons that escape the slab on one side can then be described in terms of a weighted sum of all eigenstates of the slab with eigenvalues λ_k and ${\lambda }_k^{*}$ with k ∈ ℕ and the albedo A. The contribution I⁽ⁿ⁾ of the nth scattering order to the total intensity is given by (Krieger & Wolf 2021) $$I^{(n)}=A^n{\displaystyle \sum _{k=1}^{\infty }{(-1)}^{k+1}}\left({\lambda }_k^{n+1}\sqrt{1-{\lambda }_k}{a}_k+{\lambda }_k^{*n+1}\sqrt{1-{\lambda }_k^{*}}{b}_k\right).$$(1)

The coefficients a_k and b_k of this expansion can be derived from an initial photon distribution, hereby given as: $$v_0(x)=2\delta (x+b).$$(2)

This particular distribution was chosen since it represents a scenario in which the slab is illuminated from one side by a monochromatic radiation source. To this end, the source is placed on one side of the slab at the optical depth coordinate x = −b with b = τ_max/2. Furthermore, it is normalized due to a factor of 2 in Eq. (2), which ensures that the totality of all simulated photons that traverse the slab after the first event of emission carry a weight of one, as the underlying phase function follows an isotropic distribution. The coefficients a_k and b_k in Eq. (1) are then given by (Krieger & Wolf 2021) $$a_k=\langle \cos \left({\sigma }_{k}x\right),v_0(x)\rangle =\frac{2\mathrm{ cos }\left({\sigma }_{k}b\right)}{{\lambda }_k+b}={(-1)}^{k+1}\frac{2\sqrt{1-{\lambda }_k}}{{\lambda }_k+b}$$

and $$b_k=\langle \sin \left({\sigma }_k^{*}x\right),v_0(x)\rangle =\frac{-2\mathrm{ sin }\left({\sigma }_k^{*}b\right)}{{\lambda }_k^{*}+b}={(-1)}^k\frac{2\sqrt{1-{\lambda }_k^{*}}}{{\lambda }_k^{*}+b}.$$

In combination with Eq. (1) this results in the following expression for I⁽ⁿ⁾: $$I^{(n)}=2A^n{\displaystyle \sum _{k=1}^{\infty }\frac{{\lambda }_k^{n+1}\left(1-{\lambda }_k\right)}{{\lambda }_k+b}}-\frac{{\lambda }_k^{*n+1}\left(1-{\lambda }_k^{*}\right)}{{\lambda }_k^{*}+b}.$$(3)

This expression can then be used to calculate the probability $\Delta p^{(n)}=I^{(n)}/{\displaystyle {\sum }_{n=0}^{\infty }{I}^{(n)}}$ for a transmitted photon package to have scattered n times before leaving the slab, and with it the SOD, namely, the distribution of scattering orders n of transmitted photon packages. Figure 1 exemplarily shows the SODs ∆p⁽ⁿ⁾/∆log₁₀ (n) for three different values of the optical depth τ_max (upper plot) and for three different value of the albedo A (lower plot), which were calculated as described in Krieger & Wolf (2021). These plots clearly show a shift of the SOD toward higher values of the scattering order n for increasing values of optical depth and albedo, which will be addressed in greater detail later on in this study. In the following sections, we derive an approximate solution for the SOD and discuss its implications for the complexity of MCRT simulations in the regime of high optical depths and the related scattering order problem. The results of these derivations will form the basis for a novel proposed method that aims at solving or at least alleviating this problem and should thus improve MCRT simulations in the regime of high optical depths.

Fig. 1 SOD ∆p(n)/∆ log10 (n) as function of the scattering order n. Top: impact of τmax using A = 1. Bottom: impact of A using τmax = 40.

2.2 The scattering order distribution

Every eigenvalue λ is related to its corresponding (positive) σ value, in accordance with the equation λ = 1/(1 + σ²), which can be obtained by solving the equations (Krieger & Wolf 2021) $$b{\sigma }_k=\pi (k-1)+b\text{Δ}{\sigma }_k\text{\hspace{1em} and}$$(4) $$b{\sigma }_k^{*}=\pi (k-0.5)+b\text{Δ}{\sigma }_k^{*}$$(5)

for k ∈ ℕ, where $$b\text{Δ}{\sigma }_k={\tan }^{-1}\left({\sigma }_k^{-1}\right)\text{ and }b\text{Δ}{\sigma }_k^{*}={\tan }^{-1}\left({\sigma }_k^{*-1}\right)$$(6)

It can be shown, that the inequalities 0 < b∆σ_k+1 < b∆σ_k < π/2 and $0<b\text{Δ}{\sigma }_{k+1}^{*}<b\text{Δ}{\sigma }_k^{*}<\pi /2$ hold, meaning that both b∆σ_k and $b\text{Δ}{\sigma }_k^{*}$ decrease with increasing value of k. We note that previous equations imply that quantities with asterisks are related to the same quantities without asterisks according to ${\lambda }_k^{*}={\lambda }_{k+0.5}$ and ${\sigma }_k^{*}={\sigma }_{k+0.5}$. However, in order to use these relations, we have to extend σ_k and ${\sigma }_k^{*}$ as well as λ_k and ${\lambda }_k^{*}$ to real numbers of k. Using Eqs. (4)–(6) and performing derivatives with respect to k, the extension is defined by the derivative $$\frac{\text{d}\sigma }{\text{d}k}=\frac{\pi }{{\left(1+{\sigma }^2\right)}^{-1}+b}$$

for $\sigma \in \left\{{\sigma }_k,{\sigma }_k^{*}\right\}$ as well as $$\frac{\text{d}\lambda }{\text{d}k}=\frac{-2\pi \sqrt{1-\lambda }{\lambda }^{3/2}}{\lambda +b}$$(7)

for $\lambda \in \left\{{\lambda }_k,{\lambda }_k^{*}\right\}$. For a fixed value of λ, the derivative in Eq. (7) is on the order of 1/b in the regime of b ≫ 1 > λ, while the second derivative $\frac{{\text{d}}^2\lambda }{\text{d}{k}^2}$ is on the order of 1/b². In this regime, the difference ∆λ_k = λ_k+1 − λ_k between two consecutive λ values can thus be approximated by $$\text{Δ}{\lambda }_k\approx \frac{\text{d}\lambda }{\text{d}k}\text{Δ}k,$$

with ∆k = 1. Using this relation, the series in Eq. (3) becomes $$I^{(n)}=-{\displaystyle \sum _{k=1}^{\infty }{U}^{(n)}}\left({\lambda }_k\right)\text{Δ}{\lambda }_k-U^{(n)}\left({\lambda }_k^{*}\right)\text{Δ}{\lambda }_k^{*}\text{, }$$(8)

with $$U^{(n)}(\lambda )\equiv \frac{A^n}{\pi }{\lambda }^{n-0.5}\sqrt{1-\lambda },$$

such that I⁽ⁿ⁾ converges in the limit of b → ∞ to the difference of two Riemann integrals, one performed over λ_k and another over ${\lambda }_k^{*}$. In particular, the upper boundary of the integration intervals of both integrals would differ by $\text{Δ}\lambda ={\lambda }_1-{\lambda }_1^{*}$ , while the lower boundary for both these integrals would be zero and thus equivalent. By approximating the right-hand side of Eq. (8) with this difference, we arrive at $$I^{(n)}\approx U^{(n)}\left({\lambda }_1\right)\left({\lambda }_1-{\lambda }_1^{*}\right).$$(9)

Since the scattering order, n, appears as an exponent in the function U⁽ⁿ⁾, it can be shown that the approximation in Eq. (9) is underestimating the actual intensity at very high scattering orders. Moreover, in the case of very low scattering orders, we find that the approximation overestimates the flux, suggesting that the approximate SOD is generally shifted toward lower scattering orders compared to its exact solution. This suggests that any derived expression for a minimum scattering order can be expected to be rather underestimated. This is discussed in the following sections.

Finally, the determination of I⁽ⁿ⁾ requires the calculation of λ₁ and ${\lambda }_1^{*}$, which can be deduced from their corresponding σ values, σ₁ and ${\sigma }_1^{*}$, respectively. In the regime of high optical depths, the conditions 0 < σ₁ ≪ 1 and $0<{\sigma }_1^{*}\ll 1$ are satisfied. Thus, we perform a series expansion of Eqs. (4) and (5) about σ = 0 to deduce both values as functions of the optical depth b. This results in equations of the form: $$b\sigma =f\pi -\sigma +\frac{1}{3}{\sigma }^3+O\left({\sigma }^5\right),$$

where f = 0.5 for σ = σ₁ and f = 1 for σ = σ₁. In the considered optical depth regime, this cubic equation has three real roots, σ_A, σ_B, and σ_C, which satisfy ${\sigma }_A<0<{\sigma }_B<\sqrt{b+1}<{\sigma }_C$, with σ_B being the searched-for solution. It can be shown that its solution is given by $${\sigma }_B=\frac{f\pi }{b}-\frac{f\pi }{b^2}+\frac{f\pi }{b^3}+\frac{f\pi \left({(f\pi )}^2/3-1\right)}{b^4}+O\left(\frac{1}{b^5}\right).$$

Making use of the relation between σ and λ values, we find $$\begin{array}{l}{\lambda }_1=1-\frac{{\pi }^2}{4b^2}+\frac{{\pi }^2}{2b^3}+O\left(\frac{1}{b^4}\right)\text{\hspace{1em}and}\\ {\lambda }_1^{*}=1-\frac{{\pi }^2}{b^2}+\frac{2{\pi }^2}{b^3}+O\left(\frac{1}{b^4}\right).\end{array}$$

Finally, using these results in combination with Eq. (9) and neglecting higher order contributions, we arrive at an approximate solution for the SOD for photons, which traversed the slab in the regime of high optical depths, that is given by $$I^{(n)}=\frac{3{\pi }^2}{8b^3}{\left(A{\lambda }_1\right)}^n.$$(10)

2.3 Complexity of radiative transfer simulations

The approximate solution for I⁽ⁿ⁾ shown in Eq. (10) for the regime of high optical depths can then be used to calculate the position of the maximum of the SOD, n_max, which satisfies $${\frac{\text{d}}{\text{d}n}\left(n{I}^{(n)}\right)|}_{n=n_{\max }}=0$$

and is thus given by $$n_{\max }={\frac{-I^{(n)}}{\text{d}{I}^{(n)}/\text{d}n}|}_{n=n_{\max }}=\frac{-1}{\ln A+\ln {\lambda }_1}.$$(11)

Equations (10) and (11) suggest that while in the regime of high optical depths the transmitted flux decreases strongly with the optical depth, the complexity of the radiative transfer process can still be limited by the albedo rather than the optical depth.

Here, we refer to the complexity in terms of the scattering order ranges that need to be simulated to properly represent the radiative transfer process, which is particularly important for deriving reliable flux estimates. This implies that the most suitable MCRT method for solving a specific problem should always consider the complexity that is introduced due to the albedo and the optical depth.

In particular, if the complexity is dominated by the optical depth, namely, if | ln A| ≪ | ln λ₁|, the position of the SOD becomes practically independent of the albedo and is determined by the optical depth according to $$n_{\max }\approx -1/\ln {\lambda }_1\to 4b^2/{\pi }^2+O(b)\text{\hspace{1em} (optical depth-dominated). }$$(12)

Comparing this result with the SODs in the upper plot of Fig. 1, we find that the actual position of the maximum exhibits an approximately quadratic dependency on the optical depth, as predicted by Eq. (12). However, we also find, that the derived approximation underestimates the scattering order of the maxima, as was expected as a result of a shift of the approximate solution toward lower scattering orders (see Sect. 2.2). Furthermore, we note that the special case of A = 1 is particularly interesting, in the sense that it always results in an optical depth-dominated scenario, independent of the actual value of τ_max, and thus the complexity of the problem increases drastically with the optical depth.

On the contrary, if the complexity of the radiative transfer process is albedo-dominated, namely, if | ln λ₁| ≪ | ln A|, the position of the maximum is mostly determined by the albedo and given by $$n_{\max }\approx -1/\ln A\text{\hspace{1em} (albedo-dominated). }$$(13)

Similarly to the optical depth-dominated case, we find that this approximate position of the maximum also underestimates the scattering order of the actual position of the maximum for the SODs shown in the lower plot of Fig. 1. This suggests that the derived expressions for n_max in Eqs. (12) and (13) can be used as lower boundaries for the actual position of the maximum of an SOD in both the optical depth-dominated and the albedo-dominated cases, respectively.

Another property that strongly impacts the complexity is the mean number of scattering orders n_mean, which is given by $$n_{\text{mean }}=\frac{{\displaystyle {\sum }_{n=0}^{\infty }n}{I}^{(n)}}{I^{\text{tot }}}=\frac{A{\lambda }_1}{1-A{\lambda }_1}\text{, }$$

where I^tot is the total transmitted intensity¹ with: $$I^{\text{tot}}={\displaystyle \sum _{n=0}^{\infty }{I}^{(n)}}.$$

We find, that n_mean is closely related to the position of the peak of the SOD n_max, which satisfies Aλ₁ ≤ n_mean/n_max ≤ 1. In the optical depth-dominated case, the position of the peak of the SOD therefore approximately traces the mean number of scattering orders.

Furthermore, the width of the SOD also strongly impacts the complexity of the radiative transfer problem. Using the approximation for I⁽ⁿ⁾ in Eq. (10), we find that the apparent width of the SOD in the double logarithmic plot of the SOD (see Fig. 1) indeed converges to a constant value in the regime of high optical depths, as it has been conjectured in Krieger & Wolf (2021). Considering the results of this study regarding the width of the SOD and position of its maximum as a function of the albedo and optical depth, we therefore find more evidence in support of the importance of the scattering order problem (Krieger & Wolf 2021), suggesting that it indeed poses one of the fundamental problems of MCRT simulations in the regime of high optical depths.

3 Method: Minimum scattering order

In the regime of high optical depths, it is crucial to simulate photon packages of high scattering orders to achieve a proper representation of the actual SOD and, thus, a reliable flux estimate. Photon packages that leave the simulated model space with a very small scattering order cannot contribute toward a sufficient coverage of the SOD and may thus lead to a misrepresentation of the distribution and the underestimation of flux (Krieger & Wolf 2021).

One possibility to counteract this is to enforce a minimum scattering order, m, which we recommend coupling with the usage of a peel-off method, meaning, photon packages may only leave the model space after their number of scattering events reached or even exceeded this threshold. Using Eq. (10), we arrive at the following relation between the coverage p ∈ [0, 1) of the total transmitted intensity and the approximate scattering order m, that needs to be simulated in order to reach it: $$p=\frac{{\displaystyle {\sum }_{n=0}^m{I}^{\left(n\right)}}}{I^{\text{tot}}}=1-{\left(A{\lambda }_1\right)}^{m+1}$$

Hence, in the regime of high optical depths the minimum scattering order, as a function of the albedo A > 0, optical depth, and coverage p, is given by $$m=\frac{\ln (1-p)}{\ln \left(A{\lambda }_1\right)}-1\stackrel{\text{ Eq. (11) }}{=}-\ln (1-p)n_{\max }-1.$$(14)

When assuming a coverage of p = 99.9%, this results in the following expressions for the minimum scattering order for both the albedo-dominated and the optical depth-dominated cases: $$m\approx \{\begin{array}{cc}\frac{28b^2}{{\pi }^2}-1& \text{ (optical depth-dominated; }A=1)\\ \frac{-7}{\ln A}-1& \text{ (albedo-dominated; }A<1\text{ and }b\to \infty ).\end{array}$$(15)

In order to apply our proposed method of enforced minimum scattering order during a MCRT simulation, simulated photon packages are kept inside the model space for a minimum of m scattering events, where m is chosen at least as high as stated in Eq. (15). During a standard MCRT simulation, an arbitrarily large optical depth τ is chosen after every single event of interaction randomly according to its corresponding probability distribution to decide the next location of interaction. However, when using our proposed method, instead the interval is limited to values of τ < τ_los, where τ_los is the total optical depth in the direction of the original photon package. This enforces an interaction within the model space and to compensate for that, the weight of the resulting photon package needs to be adjusted by a factor of $w_{\text{enforce }}=1-e^{-{\tau }_{\text{los }}}$. This factor corresponds to the probability with which the original photon package interacts at an optical depth of τ < τ_los.

We note that this method furthermore requires the simulation of a secondary, non-interacting photon package to compensate for the loss of weight carried by the original photon package. Its direction corresponds to that of the original photon package, while its weight is equal to that of the original package, reduced by a factor of $w_{\text{sec }}=e^{-{\tau }_{\text{los }}}$. The simulation of such a secondary photon package, however, is no longer required if a peel-off method is used simultaneously, since in this case its weight must not add to the flux estimate, rendering its simulation unnecessary. Additionally, if the peel-off method is used, every simulated event of interaction contributes to the SOD, leading to an improvement of the overall coverage of scattering orders.

3.1 Application: Slab experiment

In this section, we apply our proposed method in a simulation of an infinite plane-parallel 3D homogeneous slab. This setup is often used as a testbed, particularly for the study of MCRT simulations in the regime of high optical depths (for details regarding the setup, see Camps & Baes 2018; Krieger & Wolf 2021). The slab is illuminated from one side by an isotropic monochromatic radiation source, where the slab has a total transverse optical depth of τ_max = 30. In this setup, we simulate the transmitted intensity as a function of µ = cos θ, which is the cosine of the penetration angle θ. For this test, the single scattering albedo is set to A = 0.7 and the scattering of photon packages is simulated according to an isotropic phase function. To estimate the overall performance of the MCRT simulations, we additionally calculate the expected transmitted intensity based on a non-probabilistic solution. To this end, we adopted the procedure used by Camps & Baes (2018) and Krieger & Wolf (2021) that is based on work done by Roberge (1983) and di Bartolomeo et al. (1995).

Figure 2 shows results for the determined transmitted intensity I_trans(µ) for the non-probabilistic solution (purple solid line) as well as for four different types of MC based radiative transfer simulations. The latter differ with regard to their utilization of our proposed minimum scatter order method (MSO) as well as regarding the usage of the composite-biasing technique (Camps & Baes 2018) for stretching photon package paths. For each of these types of simulations, we performed 10⁴ simulations with 10³ photon packages per simulation for our further analysis. These data were subsequently divided into ten batches of 10³ simulations each, across which the average transmitted intensity was calculated. This effectively yields intensity curves corresponding to N = 10⁶ simulated photon packages each. Then, the four displayed dashed curves represent the respective median curves of all ten thereby obtained intensity curves, while the hatched colored areas mark regions in the plot, where all individual simulated curves lie. We note, that it is crucial to consider median curves to study the underestimation of transmitted intensity values, since they are a measure of the most likely outcome of a simulation for a fixed number of simulated photon packages. Comparing the non-probabilistic solution in Fig. 2 with the MC based intensity curves, we find that the simulated optical depth may result in severely underestimated values for the intensity, particularly if the stretching method is not used. In this case, even the usage of the method of enforced minimum scattering orders does not yield results that are close to the reference solution without simulating a higher number of photon packages. On the contrary, our results suggest, that using the stretching method leads to significantly better estimates and that the additional utilization of our proposed method (orange dashed curve) results in a further improvement of their quality. Altogether, enforcing a minimum scattering order seems to benefit the simulations twofold: First, the flux estimates are higher across all penetration angles. Second, the noise of the intensity estimates, which is inherent to the MC based approach, is smaller, which results in the overall narrower hatched orange region as compared to the blue region in the plot.

Figure 3 shows the ratio of the median intensity curves $I_{\text{trans}}^{+\text{MSO}}/I_{\text{trans}}^0$ of simulations that utilize our proposed method $\left(I_{\text{trans }}^{+\text{MSO}}\right)$ and those that do not make use of it $\left(I_{\text{trans }}^0\right)$. In this plot, values above one (gray dotted line) indicate a gain of intensity associated with the usage of our method. The upper plot in Fig. 3 shows the case when the stretching method is not used, and the lower plot that when it is used. Generally, we find that the intensity ratio is typically greater than one, indicating a gain of intensity associated with the usage of our proposed method. This is especially the case, when the stretching method is used, occasionally even leading to intensity ratios above 5. When stretching is not used, the gain still amounts to values between ~20–90%. In further tests, where we reduced the number of simulated photon packages per simulation to 10⁵ or even to 10⁴, we found, that the severity of the underestimation of flux estimates clearly increases with lower photon package numbers. Likewise, an increase of the number of simulated photon packages has led to further improvement. Overall, the enforcement of the minimum scattering order therefore appears to be most effective, when the number of simulated photon packages generally does not suffice to ensure a proper representation of relevant photon package paths. In turn, this implies that its utilization can therefore greatly reduce the required number of simulated photon packages.

The previous analysis has demonstrated that the application of our proposed method can lead to improved flux estimates, when keeping the number of simulated photon packages constant. However, the additional simulation of interaction events is also accompanied by an increase of computation time, which roughly scales with the number of simulated events. To assess the effect of its utilization on the computational demand, Fig. 4 shows the dependence of calculated intensity estimates as a function of the number of simulated events of interaction N_interact. In particular, the upper plot shows median curves of the ratio $I_{\text{trans}}/I_{\text{trans}}^{\text{ref}}$ for the direction µ ≈ 0.99, which roughly corresponds to the direction perpendicular to the slab surface, while assuming the usage of the stretching method. Here, $I_{\text{trans}}^{\text{ref}}$ is the respective linearly interpolated non-probabilistic reference solution. Additionally, the color-coded shaded areas in the plot mark the regions where the central 50% (also called the interquartile range) of corresponding data points lie. We note that to generate these plots, the results of the previously mentioned 10⁴ performed simulations have been gathered in batches, assuming various sizes between 1 and 10⁴, and subsequently averaged. Additionally, for each batch, the corresponding numbers of simulated interactions of all simulations have been added up. The averaged intensity curves obtained in this way, with the corresponding numbers of interactions, were binned using 20 logarithmically sampled bins with regard to N_interact, allowing for the determination of the interquartile ranges, as shown in this plot.

The upper plot exhibits a clear increase in the median ratio $I_{\text{trans}}/I_{\text{trans}}^{\text{ref}}$ for both types of simulations. Simulations that utilize our method, though, achieve higher flux estimates for the same number of simulated events and are therefore closer to the reference solution. We note that data points in the area of the plot where N_interact ≳ 10⁷ originate from the averaging of large batches of simulations, of which there are only few. As a result, this procedure can be expected to lead to comparably poor statistics at very high values of N_interact. This also explains the decrease in slope at N_interact ≈ 10⁷ of both curves, which is an artifact that is expected for this analysis procedure. Apart from the clear improvement of the median curve, also the interquartile range of intensity values is much closer to the reference solution (see orange shaded area). Its width in the logarithmic plot, however, appears to be similar for both simulation types. The lower plot shows similar results, but now assuming a direction of µ = 0.5, which corresponds to a penetration angle of 60°. Both plots show similar trends, but the intensity is significantly more underestimated for the direction µ = 0.5 for all considered values of N_interact. Nonetheless, the underestimation becomes less severe as N_interact increases, as expected, which is similar to the case of µ ≈ 0.99. This again suggests that for increasing penetration angles, the number of interactions events that needs to be simulated to achieve a certain level of quality increases. While the actual solution will be reached with or without the utilization of our method, given a sufficient number of photon packages are simulated, its usage appears to significantly reduce the required number of simulated interactions.

Moreover, N_interact scales linearly with the computational time and is independent of the specific computer hardware, which makes it a suitable quantity for measuring the computational demand of MCRT simulations. Based on that, the plot also suggests that the usage of our method decreases the computational demand needed to achieve a certain level of quality of flux estimates. In particular, for a constant ratio of we find that fewer interactions (by about two to three times) have to be simulated for achieving the same quality, which corresponds to a likewise reduction of computation time, when enforcing our proposed minimum number of scattering orders.

Fig. 2 Results for the infinite plane-parallel slab experiment for an optical depth of τmax = 30 and an albedo of A = 0.7. Curves represent transmission intensities obtained by either a non-probabilistic method (purple) or MC based methods (dashed curves). The latter are additionally divided into four types, depending on their usage of advanced MCRT methods, such as the stretching method (stretching) or the enforcement of a minimum scattering order (MSO). The underlaid colored areas mark the regions of all simulated intensity curves, and dashed curves of the same color represent the corresponding medians.

Fig. 3 Ratios of transmission intensity curves $I_{\text{trans}}^{+\text{MSO}}/I_{\text{trans}}^0$ as function of the direction of penetration µ, where $I_{\text{trans}}^{+\text{MSO}}\left(I_{\text{trans}}^0\right)$ denotes the curve corresponding to simulations applying (not applying) the enforcement of a minimum scattering order. Values above one (gray dotted line) indicate a gain of intensity associated with the method. The upper (lower) plot shows the results for simulations that do not (do) additionally make use of a stretching method.

Fig. 4 Ratio $I_{\text{trans}}/I_{\text{trans}}^{\text{ref}}$ a function of the number of simulated interaction events Ninteract, where $I_{\text{trans}}^{\text{ref}}$ denotes the value of the reference solution. The lines correspond to median curves, and the (color-coded) shaded areas encompass the central 50% of corresponding data points. The upper (lower) plot shows the results for a penetration direction of µ ≈ 0.99 (µ = 0.5).

3.2 Application: Embedded protoplanet in an edge-on disk

In order to assess the benefits of our proposed method in an astrophysical setting, we implemented the method in the 3D Monte Carlo radiative transfer Code Mol3D (Ober et al. 2015), a successor of MC3D (Wolf et al. 1999; Wolf 2003). Mol3D uses a locally divergence free continuous absorption scheme of photon packages (Lucy 1999), applies a method of immediate reemission according to a temperature corrected emission spectrum (Bjorkman & Wood 2001), and is particularly optimized for systems of extremely high optical depth by utilizing a large database of precalculated photon package paths (Krieger & Wolf 2020). Furthermore, Mol3D uses the composite-biasing method (Camps & Baes 2018), with which it has been equipped particularly for the purpose of this analysis. Details regarding Mol3D can be found in Ober et al. (2015).

The simulated system is composed of an edge-on circumstellar disk at a distance of d = 140 pc that harbors an embedded protoplanet. The disk density distribution is described by $$\rho (r,z)=\frac{\text{Σ}(r)}{\sqrt{2\pi }h(r)}\cdot \exp \left[-\frac{1}{2}{\left(\frac{z}{h(r)}\right)}^2\right],$$(16)

where r is the radial distance to the center of the disk perpendicular to the z-axis and z describes the distance to the mid-plane of the disk. Σ(r) is the surface density and h(r) the density scale height, which are given by (Lynden-Bell & Pringle 1974; Hartmann et al. 1998) $$\text{Σ}(r)={\text{Σ}}_0{\left(\frac{r}{r_0}\right)}^{\beta -\alpha }\cdot \exp \left[-{\left(\frac{r}{r_0}\right)}^{2+\beta -\alpha }\right]$$(17)

and $$h(r)=h_0{\left(\frac{r}{r_0}\right)}^{\beta },$$(18)

where variables with subscript 0 denote reference values evaluated at the reference radius r₀, α is the compactness parameter, and β the flaring parameter of the disk. The total disk mass M_disk = 10⁻² M_⊙ includes the mass of dust and gas. Generally, the chosen values of all model parameters of the circumstellar disk correspond to commonly observed values that are typical for disks around T Tauri stars (Andrews et al. 2009). The dust is composed of spherical dust grains which themselves consist of 62.5% silicate and 37.5% graphite, the latter with a ratio of parallel:ortho = 1:2 (Draine & Malhotra 1993), for which we use mixed optical properties (Draine & Lee 1984; Laor & Draine 1993; Weingartner & Draine 2001). The grains have a bulk density of ρ_bulk = 2.5 g cm⁻³ and radii a between 5 and 250 nm, that follow the distribution dn ~ a^−3.5da (Mathis et al. 1977). Mean optical properties of these dust grains were calculated (Wolf 2003) and used throughout our simulations. In these simulations, we assume Mie theory for the description of simulated scattering events.

The embedded protoplanet orbits the center of the edge-on circumstellar disk at a distance of r_p = 20 au and is placed between the observer and the disk center inside the mid-plane. It has a surface temperature of 2000 K and a luminosity of L_p = 10⁻³ L_⊙, which is consistent with evolutionary tracks of planets Burrows et al. (1997). A list of all relevant model parameters can be found in Table 1.

Here, we consider two types of simulations for which we determine the observed total (specific) flux density that originates from the embedded planet. One type uses our proposed method of enforced minimum scattering order, while the other is not. Since the planet is embedded deeply inside the disk, we can expect the optical depth to be high in all directions, while the albedo value is rather small for the considered wavelength range. Therefore, we use the expression for the minimum scattering order m in Eq. (15) for the albedo-dominated case for our simulations. We note, that due to the general lack of a non-probabilistic reference solution for such a system, it is not possible to access the performance of MCRT simulations in the same manner as has been done for the case of the infinite plane-parallel slab (see Fig. 2). Therefore, we only focus on the flux gain associated with the usage of our proposed method.

In total, we performed 400 simulations per type, each for 45 wavelengths λ in the mid-infrared range between ~3 and 50 µm using 10⁵ photon packages per simulation and wavelength. Subsequently, we determined the ratio of the two resulting median flux density curves F^+MSO /F⁰, where F^+MSO (F⁰) denotes the curve obtained from simulations using (not using) the method of enforced minimum scattering order. The results of these simulations are shown in Fig. 5. The plot shows, that for the majority of considered wavelengths, the usage of our proposed method leads to an increase of the determined flux estimate by a factor that partly even exceeds the value of 100. Similar to the findings in Sect. 3.1, we additionally find in this setup that the flux gain associated with the usage of our proposed method seems to increase as the expected level of flux decreases. Further tests were conducted in which we varied the surface temperature and luminosity of the planet as well as its orbital position, both radially and azimuthally, which only provided further evidence in support of this conclusion.

Fig. 5 Results from simulations of an embedded protoplanet in an edge-on circumstellar disk, showing the ratio of two median flux density curves (F+MSO /F0) corresponding to radiation originating from the planet. Here, F+MSO (F0) denotes the flux density curve obtained from simulations using (not using) the method of enforced minimum scattering order. Values above one (gray dotted line) indicate a gain in flux density associated with the usage of the method.

4 Discussion and implications

In this section, we give the reasons for using our proposed method as a simple way to implement the means of providing a minimum quality during the step of flux determination in a MCRT simulation. We additionally address potential issues concerning the applicability of the method and the interpretation of the results in Eqs. (14) and (15). Altogether, we argue in favor of using a minimum scattering order m that is at least as high as the expression derived in Eq. (15) for MCRT simulations of all systems which, for instance, contain regions of high optical depths, or systems that contain a radiation source, which is embedded inside an interacting medium.

While the actual SOD of a simulated system depends on its particular properties, our method of the enforced minimum scattering order is aimed at avoiding the consideration of the full complexity and gives an estimate for it in terms of a minimum scattering order m, which can readily be used throughout a simulation. To clarify the role of m in such a simulation, we next discuss three changes to the setup of the 1D slab, which would lead to a shift and reshaping of the SOD; while these changes are specific, the drawn conclusions are generally valid. First, we discuss the effects of placing the radiation source inside the slab, which shifts the SOD toward lower scattering orders. Second, we consider the difference between the approximate solution for the SOD (Eq. (10)) and the actual solution, which is shifted toward higher scattering orders. Third, we consider the transition from the 1D to the 3D case, which leads to a shift and reshaping of the SOD. We discuss these changes to the SOD while assuming the case of A < 1, as this applies to the majority of materials in astrophysics, which, in accordance with Eq. (15), correspond to a minimum scattering order in the albedo limited case of m = −7/ ln A − 1. The reasoning, however, can be applied to the optical depth-dominated case as well.

By placing the radiation source inside the slab rather than outside, on the side that is not facing the detector, the optical depth between the detector and the source of radiation decreases. One crucial consequence of the reduced optical depth is that photon packages of lower scattering orders gain in relevance regarding their contribution to the total detected intensity, tendentially leading to a shift of the SOD toward lower scattering orders. This could potentially result in a situation whereby m is overestimated, such that the coverage of scattering orders would exceed p > 99.9%. When making use of our proposed method, this by itself is not an unfavorable result, as it would simply increase the quality of the determined flux estimates. If, however, the required simulation time for this is expensive, the achieved gain in accuracy may not be worth the accompanied increase in simulation time. Therefore, a shift of the SOD toward lower scattering orders does not result in a lowering of the quality of estimated flux values, when using our proposed method of enforced minimum scattering order.

If, on the other hand, the SOD is shifted toward higher scattering orders, as can be expected (see Sect. 2.2) for the actual distribution, compared to its approximation in Eq. (10), the proposed value of m may be below what would be required to achieve the desired coverage of p = 99.9%. As a result, the derived flux values may still be underestimated, however, due to the enforced contribution of scattering orders up to the order of m, the level of underestimation can be expected to be lower. Therefore, the usage of the method would at the very least provide an improvement of quality of determined flux values (as has been shown in Sects. 3.1 and 3.2) and offset the scattering orders of photon packages closer to the relevant scattering order interval.

The same reasoning holds when transitioning from the 1D to the 3D scenario. This can be understood, for instance, by assuming a slightly changed setup and comparing the scenario of a central radiation source that is located inside a 1D slab of half transverse width b with that of a 3D sphere of radius b, both in units of optical depth. Since in the 1D case, the photon package path is confined, such that it always moves toward or away from the center, its movement in the radial direction leads to an average shift of its position of ∆τ_r = 1 between two consecutive events of interaction. In the 3D case and in the regime of high optical depths, the average displacement of the radial position of the photon package is comparably smaller. In particular, if the photon package is already at a radial optical depth position of τ_r ≫ 1, the average displacement reduces to ∆τ_r = 0.5. The reduction of average radial displacement is a result of the fact, that the direction of the photon package in the 1D case always is parallel to the radial direction, which it is not in general in the 3D case. As a result, the SOD in the 3D case can be expected to be, firstly, narrower compared to that of the 1D case, and secondly, shifted toward higher scattering orders, since the radial displacement of photon packages is delayed due to the smaller radial pacing at high optical depths. Therefore, we can conclude, that despite the SOD exhibiting a different shape and being shifted, using the minimum scattering order m from Eq. (15) for 3D MCRT simulations can still be expected to elevate the quality of determined flux values (compare with Sects. 3.1 and 3.2) and serve as a lower boundary for the measure of the complexity of the underlying SOD.

Our findings clearly suggest that using the quantity m from Eqs. (14) and (15) as a measure of the minimum scattering order for the proposed method of enforced minimum scattering order generally leads to more reliable and higher flux values in both the 1D and the 3D cases. This is particularly relevant in the regime of high optical depths, which can be explained as a result of the alleviation of the scattering order problem. This has implications not only for simulations in which photon packages themselves are used to construct observables, but also when relying on secondary results derived from their simulation. One such example would be the utilization of a raytracing scheme for deriving flux maps, which is based on a scattering source function that is derived during a MCRT simulation. For such a case, our findings clearly suggest that this source function would also be lacking in terms of the contribution of photon packages, which would ultimately manifest itself in underestimated flux values; the latter could likewise be improved by the usage of our proposed method. Moreover, we note that while it is possible to adjust the method for the estimation of m by assuming a different coverage p of the SOD, we generally recommend using high values close to unity. Ultimately, it is still necessary to limit its value.

Furthermore, our results have implications for the performance of other optimization methods used to improve MCRT in the regime of high optical depths. This is the case, for instance, in the context of biasing methods that affect the determination of photon package paths, such as stretching methods (e.g., Camps & Baes 2018), which increase the average distance between two consecutive points of interaction of a simulated photon package. As we discussed previously (Krieger & Wolf 2021), the usage of such a stretching method, despite its great potential benefits, has its limitations, which can be the result of a too strong decrease of sampled scattering orders of simulated photon packages. This can lead to a scenario in which the SOD is not sufficiently well represented as the relevant scattering order intervals are practically not reached due to the usage of the stretching method. This may very well explain the improvement of flux estimates that was associated with the usage of our method in the embedded protoplanet case described in Sect. 3.2. In further tests, where we performed the same simulations but did not apply the stretching method, we found no difference between simulations that enforced a minimum number of scattering orders and those that did not. This can be explained by the fact that simulated photon packages that are emitted by the planet cannot escape without a high number of interactions.

The findings of this study therefore suggest that the usage of a stretching method, which results in a sampling of scattering orders well below the minimum scattering order, m, can be considered as an indication for too high stretching factors. Moreover, the best suited stretching methods should take into account the albedo of the medium, as it strongly shapes the complexity of the simulation, see Sect. 2.2. However, our results do suggest that in the albedo-dominated case, namely, low albedo and high optical depth, the usage of a method that increases the overall pace with which a simulated photon package is transported, generally offers a great potential for MCRT simulations, if used correctly, as can be seen in Fig. 2. This implies that the best-suited method to overcome the problem of high optical depths and the flux underestimation problem seems to utilize a stretching technique combined with a method to ensure a sufficient number of simulated interactions events.

As a precautionary note, it is worth acknowledging that when dealing with more complex systems beyond the scope of our current consideration, it may be generally difficult to decide whether the problem is albedo-dominated or optical depth-dominated. This could for instance be the case, if the albedo is direction-dependent or the optical properties of the simulated system vary across the medium. However, as long as the optical depth is generally high, we believe that it is worth considering the usage of a stretching technique combined with the enforcement of at least as many scattering orders m as derived for the albedo-dominated case, see Eq. (15). Furthermore, if this is leading to higher flux estimates, it would be crucial to assess whether a further increase of m would yield even higher estimates.

Moreover, we note that in the case of a low albedo values and depending on the particular MCRT implementation, the typical number of simulated scattering events may also exceed the minimum scattering order for a large portion of simulated scattering orders, while estimated flux values are still underestimated. This may hint at another problem, which lies in a suboptimal method for path determination, which in this case leads to an over-sampling of paths of high scattering orders and an insufficient representation of paths of low scattering orders. Thus, it may be beneficial to make use of a different MCRT implementation and apply, for instance, better suited biasing techniques. Additionally, the question regarding a concept of a maximum scattering order may also arise, such that photon packages that have undergone a certain number of scattering events may safely be removed from the model space. The knowledge of such a maximum scattering order may then be used to increase the efficiency of MCRT simulations by spending less time on the simulation of insignificantly and more on the simulation of significantly contributing photon packages. Overall, we find that it would be highly desirable to design an adaptive stretching method in the future, which takes into account the findings of this study in order to adjust the stretching as needed to ensure the most reliable flux estimates.

5 Summary

In this study, we analyzed the scattering order distribution (SOD) of transmitted photon packages, which traversed a 1D slab that is illuminated from one side. On the basis of this testbed, we derived an approximate solution for the SOD and studied its implications for the complexity of Monte Carlo radiative transfer (MCRT) simulations, as described in Sects. 2.2 and 2.3. Here, we find further evidence in favor of the conjecture that the scattering order problem constitutes one of the issues that lies at the core of MCRT simulations, particularly in the regime of high optical depths. Moreover, we identified two different cases that can dominate the complexity of a MCRT simulation at high optical depths with notably distinct properties, i.e., the albedo-dominated case (| ln λ₁| ≪ | ln A|) and the optical depth-dominated case (| ln A| ≪ | ln λ₁|) On this basis, we derived the following three properties of the approximate SOD, all of which are in agreement with previously made predictions (see Krieger & Wolf 2021). First, the peak of the SOD at the scattering order n_max can be used as a measure of the average scattering order of the distribution. Second, the value of n_max scales approximately quadratically with the optical depth τ_max in the optical depth-dominated case, while being limited by n_max ≲ −1/ ln A independent of the value of the optical depth in the albedo-dominated case. Third, the apparent width of the SOD in the double logarithmic plot converges in the limit of τ_max → ∞, which stresses the severity of the scattering order problem, especially in the optical depth-dominated case.

Furthermore, we used the SOD approximation to derive an estimate for the minimum scattering order m, which is required to achieve a certain coverage of the SOD and the total transmitted intensity (see Sect. 3). We thus propose a method of enforced minimum scattering order that ensures a minimum scattering order of simulated photon packages of at least the derived value of m. Our proposed method aims at increasing the reliability of flux estimates, by alleviating the scattering order problem, which thus serves as a simple to implement means of providing a minimum quality both in 1D and 3D simulations.

We then tested the method by applying it in two different simulated systems: 1) an infinite plane-parallel slab (see Sect. 3.1) and 2) a protoplanet embedded in an edge-on circumstellar disk (see Sect. 3.2). In both systems, the enforcement of the minimum scattering order led to an alleviation of the underestimation of flux and its effectiveness seemed to increase for lower observed values of flux. Overall, we find that this concept bares great potential for an improved flux determination, particularly in the regime of high optical depths and in combination with a photon path length stretching technique, as described in Sect. 4. Finally, we have discussed the possibility of more advanced optimization techniques for the flux determination step that build on our findings regarding the SOD and minimum scattering order. With the help of such techniques, it might then be possible to overcome the core problems in the regime of high optical depths and as a result maximize the overall potential of MCRT simulations.

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

References

All Tables

All Figures

	Fig. 1 SOD ∆p(n)/∆ log10 (n) as function of the scattering order n. Top: impact of τmax using A = 1. Bottom: impact of A using τmax = 40.
In the text

Fig. 2 Results for the infinite plane-parallel slab experiment for an optical depth of τmax = 30 and an albedo of A = 0.7. Curves represent transmission intensities obtained by either a non-probabilistic method (purple) or MC based methods (dashed curves). The latter are additionally divided into four types, depending on their usage of advanced MCRT methods, such as the stretching method (stretching) or the enforcement of a minimum scattering order (MSO). The underlaid colored areas mark the regions of all simulated intensity curves, and dashed curves of the same color represent the corresponding medians.

In the text

Fig. 3 Ratios of transmission intensity curves $I_{\text{trans}}^{+\text{MSO}}/I_{\text{trans}}^0$ as function of the direction of penetration µ, where $I_{\text{trans}}^{+\text{MSO}}\left(I_{\text{trans}}^0\right)$ denotes the curve corresponding to simulations applying (not applying) the enforcement of a minimum scattering order. Values above one (gray dotted line) indicate a gain of intensity associated with the method. The upper (lower) plot shows the results for simulations that do not (do) additionally make use of a stretching method.

In the text

Fig. 4 Ratio $I_{\text{trans}}/I_{\text{trans}}^{\text{ref}}$ a function of the number of simulated interaction events Ninteract, where $I_{\text{trans}}^{\text{ref}}$ denotes the value of the reference solution. The lines correspond to median curves, and the (color-coded) shaded areas encompass the central 50% of corresponding data points. The upper (lower) plot shows the results for a penetration direction of µ ≈ 0.99 (µ = 0.5).

In the text

Fig. 5 Results from simulations of an embedded protoplanet in an edge-on circumstellar disk, showing the ratio of two median flux density curves (F+MSO /F0) corresponding to radiation originating from the planet. Here, F+MSO (F0) denotes the flux density curve obtained from simulations using (not using) the method of enforced minimum scattering order. Values above one (gray dotted line) indicate a gain in flux density associated with the usage of the method.

In the text