© ESO 2020

1. Introduction

Nowadays, the most widely used method to perform radiative transfer simulations in astrophysics is the Monte Carlo approach. The area of application covers a wide range of astrophysical objects, such as protoplanetary or circumbinary disks, Bok globules, filaments, or planetary atmospheres, which are subjects of studies of various codes, for example, MC3D (Wolf et al. 1999), Mol3D (Ober et al. 2015), and POLARIS (Reissl et al. 2016). During these simulations, the photon packages that are generated travel through the simulated model space and interact with matter, for example dust, by scattering, absorption, and reemission. The particular path of the photon package as well as its interactions are determined in a probabilistic way. Throughout this process, thermal energy is deposited in the dust, which increases its temperature. At every point of interaction, the photon package has a chance to be absorbed or scattered and thus continue its path through the model space. Therefore, every photon package has to be tracked from its first emission to the moment it leaves the simulated space. The number of interactions it undergoes increases strongly with the optical depth of the medium, resulting in a high computational demand (e.g., Baes et al. 2016; Gordon et al. 2017; Camps & Baes 2018). This problem is usually handled by imposing a maximum number of interaction whose exceedance results in the removal of the photon package. This generally leads to an underestimation of the ambient temperature and consequently to an altered simulated observational appearance. Additionally, the estimation of its caused error is difficult at the very least.

Being able to simulate optically thick regions is of great relevance, in particular, if the radiation stems from these regions. For example in the case of protoplanetary disks, the heating of dust grains by deeply embedded protoplanets offers a great challenge to the Monte Carlo approach. These young accreting protoplanets are embedded in an optically thick circumplanetary disk or envelope, which effectively shields their radiation from distant regions as well as observers. In order to infer properties of these dense regions, a proper temperature calculation is crucial. Multiple approximative methods have been proposed and implemented in order to solve the problem of high optical depths, for example, the partial diffusion approximation (Min et al. 2009), the modified random walk (Min et al. 2009; Robitaille 2010), or biasing techniques (Baes et al. 2016), each of which has its merits and drawbacks.

The diffusion approximation corrects the temperature of an optically thick region after the Monte Carlo simulation. It assumes that the probability for a photon package to leave this region without an interaction is zero and it obtains a temperature by solving the time-independent radiative diffusion equation using the temperature at the boundary of this region (Wehrse et al. 2000; Min et al. 2009). By using this method, a low number of photon packages can already lead to a smooth temperature distribution, whereas otherwise, in a regular Monte Carlo simulation a low photon count results in a high level of noise. This method gives approximated temperature values and is typically used if the smoothness of the temperature distribution is required for further analysis, for example, when iteratively solving for a vertical hydrostatic density distribution of a circumstellar disk (Min et al. 2009).

A different method to handle the problem of high optical depths is the modified random walk. In a region with constant density and optical properties, a sphere with a certain radius with the position of the photon package marking its center can be defined. If the optical depth of the sphere is high, the photon package scatters multiple times before it leaves the sphere for the first time. The modified random walk uses the results from the work of Fleck & Canfield (1984) for a random walk procedure representing a large number of scattering, absorption, and reemission events. Thus, it is possible to immediately move the photon package to the rim of the sphere and skip the time consuming calculation of the scattering events that the photon package would undergo on its way toward the rim. This method can greatly enhance the performance of Monte Carlo simulations (e.g., Min et al. 2009; Robitaille 2010). It was originally developed for isotropic scattering (Fleck & Canfield 1984) and can be extended to nonisotropic scattering assuming that the phase function for scattering only depends on the angle between the incident and outgoing particle direction (Ishimaru 1978). Unfortunately, there is no way to apply this method to Mie scattering (Mie 1908), which is anisotropic and whose outcome of a scattering event depends on the polarization state of the incident photon. Additionally, the exact properties of a photon package, which is emitted at the rim of the sphere, are not defined; this includes the package’s stokes vector and direction of emission, in particular.

Lastly, the biasing technique can solve the problem of high optical depth by favoring high optical depths and in turn reducing the weight of the photon package by a corresponding factor (Baes et al. 2016). This method is very powerful since it gives the possibility of altering the pathway of a photon package such that it caters to the individual goals of a simulation. This freedom comes with the difficulty of quantifying the actual error of the result and thus should be applied with caution. If the photon package scatters multiple times, the weighting procedure has to be applied after every scattering event, which can result in particularly high weighting factors. Furthermore, this method clearly demotes certain paths that the photon package would have traveled along by using unbiased Monte Carlo radiative transfer (Baes et al. 2016).

In the following, we describe a new method to solve the problem of high optical depth, which is (i) Monte Carlo based, (ii) can be used during a regular Monte Carlo simulation, (iii) does not interfere with multiple other prominent optimization techniques (e.g., continuous absorption of photon packages, immediate reemission, peel-off), (iv) does not introduce an additional bias, and (v) is (theoretically) applicable to arbitrarily high optical depths. Section 2 describes the basic principles of a Monte Carlo radiative transfer code. It follows a derivation of the necessary variables needed for the implementation of our method as well as a proof of their scalability, which is crucial since it ensures the viability of the method. We proceed to describe the method and discuss its limitations as well as upgrades before we measure its performance in a testbed. In Sect. 3 we apply our method to a simple model of a viscously heated circumstellar disk. We compare the results of the basic Monte Carlo method with the results of our method and study the impact of a limit on the number of interactions on the calculated disk temperature. We then point out the necessity of using a method that is capable of dealing with high optical depths. We discuss the short as well as the long range impact of an inaccurate treatment of high optical depths on the resulting disk temperature and thus on simulated observations of state-of-the-art instruments, such as ALMA (Kurz et al. 2002), PIONIER (Le Bouquin et al. 2011), and MATISSE (Lopez et al. 2014). In Sect. 4, we summarize the results of this paper.

2. Monte Carlo random walk

For the numerical implementation and performance test of the new method, we used the three-dimensional Monte Carlo radiative transfer Code Mol3D (Ober et al. 2015), which is a successor of MC3D (Wolf et al. 1999; Wolf 2003). It achieves a high computational efficiency by combining a locally divergence free continuous absorption of photon packages (Lucy 1999) as well as immediate reemission according to a temperature corrected emission spectrum (Bjorkman & Wood 2001). Details about the code can be found in Ober et al. (2015).

In this work, we describe a method that has been implemented in Mol3D, which uses precalculated results of Monte Carlo simulations of photon packages traveling through a homogeneous sphere of constant density, temperature, and optical properties in order to solve the problem of high optical depths in a dusty medium. We discuss the applicability of the method to arbitrarily high optical depths in arbitrary media as well as the feasibility in deriving an error estimate.

The basic concept of this method is to precalculate the random walk of a photon package traveling from the center of a sphere to its rim while keeping track of essential variables that are needed during a radiative transfer simulation. Since this process is repeated multiple times during a Monte Carlo radiative transfer simulation, we can choose a randomly selected outcome from our precalculation. The difficulty with this approach is to find a meaningful and efficient way to store the data of possible outcomes for this process that are additionally quickly accessible. The most intuitive way is to store the energy $E_{\mathrm{abs}}^{\mathrm{tot}}$ deposited in the sphere as well as the complete list of properties of the photon package at the moment it penetrates the rim of the sphere. This would include the wavelength λ, the Stokes vector S = (I, Q, U, V)^T, as well as the angle of penetration θ_out. Certainly, the outcome of this simulation depends on the temperature as well as the radius and density of the sphere. These variables a priori depend on each other, thus, they cannot be stored for each variable separately. When binning the outcome and storing their emerging probability distributions, the memory demand easily becomes enormous depending on the accuracy we intend to achieve, thus, making a different approach desirable.

2.1. General aspects of a Monte Carlo radiative transfer simulations

During a Monte Carlo radiative transfer simulation photon packages travel through a grid composed of cells for which each has constant physical properties. After the emission of a photon package, the direction and the optical depth it travels is chosen in a probabilistic manner determining its next location of interaction. Between two events of interaction, energy is stored in every dust grain inside a traversed cell according to the absorption probability of the photon package by a single grain. This probability depends on the path length l the photon package travels inside the cell. The deposited energy is given by (Lucy 1999)

$$\begin{array}{c}\hfill \mathrm{\Delta }E=E_{\gamma }\frac{C_{\mathrm{abs}}l}{V_{\mathrm{cell}}},\end{array}$$(1)

where C_abs is the wavelength-dependent absorption cross-section, V_cell is the volume of the cell, and E_γ is the energy the photon package carries. Throughout this process, the energy of the photon package is kept constant. When traveling through an optically thick cell, the photon package may interact multiple times before crossing one of the boundaries of the cell. Immediately after an absorption event, the photon package is reemitted isotropically with a randomly assigned direction and wavelength (λ) according to its ambient dust temperature (T_i) dependent spectrum. Under the assumption of a radiative equilibrium, the probability density function (PDF) P(λ | T_i) for the wavelength is given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}P(\lambda |T_i)}{\mathrm{d}\lambda }=\frac{C_{\mathrm{abs}}}{K_1}{\left(\frac{\mathrm{d}{B}_{\lambda }}{\mathrm{d}T}\right)}_{T=T_i},\end{array}$$(2)

in which B_λ is the Planck function and $K_1={\int }_0^{\infty }{C}_{\text{abs}}(\text{d}{B}_{\lambda }/\text{d}T)\text{d}\lambda$ is a normalization factor (Bjorkman & Wood 2001). For better performance during a simulation, this PDF and its corresponding cumulative distribution function (CDF) are precalculated in Mol3D for a set of wavelengths and temperatures. Despite the discrete sampling of reemission spectra, the resulting temperature of dust in a cell is a continuous function. We note that throughout the course of this paper, we refer to this general description of Monte Carlo radiative transfer as the basic method. Furthermore, since we do not make use of any of the aforementioned biasing techniques, we also refer to this method as unbiased.

2.2. Precalculated random variables

For our method, we precalculated the outcome of a simulation of photon packages traveling through a sphere of specific (constant) properties and used it to generate multidimensional binned CDFs, which account for the different possible outcomes of the simulation. During the actual Monte Carlo radiative transfer simulation, we then made use of the precalculated CDFs in order to quickly move photon packages through optically thick regions. The quantities we tracked during the precalculation are the total absorbed energy of the sphere ($\mathrm{\Delta }{E}_{\mathrm{abs}}^{\mathrm{tot}}$), the location of the last event of absorption before the photon package leaves the sphere for the first time (r_la), as well as its wavelength when penetrating the rim of the sphere (λ). For a medium of constant and homogeneous optical properties, the distribution of these variables solely depends on the radius (r_sph), temperature (T), and number density (n_dust) of the sphere. The photon package is assumed to be unpolarized after reemission at the last location of absorption. Therefore, the distribution of polarization states is trivial and does not have to be tracked. However, over the course of our precalculations, we made use of Mie scattering including the proper treatment of the polarization states of photon packages. In the following, we introduce three random variables that are crucial for our method, the effective optical depth ${\widehat{\tau }}_{\mathrm{ext}}$, the minimal missing absorption optical depth ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$, and the relative total absorbed energy X.

2.2.1. Effective optical depth ${\widehat{\tau }}_{\mathsf{ext}}^{\mathsf{sph}}$

We performed the calculations for spheres of specific optical depths (${\widehat{\tau }}_{\mathrm{ext}}$) and the set of temperatures that have already been used previously to precalculate the PDF from Eq. (2). For the simulated spheres, we chose radii that correspond to optical depths of

$$\begin{array}{c}\hfill {\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}={\widehat{C}}_{\mathrm{ext}}{r}_{\mathrm{sph}}{n}_{\mathrm{dust}}\in {\mathcal{T}}_{\mathrm{sph}}\equiv \{{10}^i\mid i\in \{1,1.5,2,2.5,3\}\},\end{array}$$(3)

where n_dust is the dust number density and Ĉ_ext is a temperature-dependent effective extinction cross-section. We define the latter by

$$\begin{array}{c}\hfill {\widehat{C}}_{\mathrm{ext}}(T_i)=\frac{1}{K_2}{\displaystyle {\int }_0^{\infty }}{C}_{\mathrm{ext}}{\left(\frac{\mathrm{d}{B}_{\lambda }}{\mathrm{d}T}\right)}_{T=T_i}\mathrm{d}\lambda ,\end{array}$$(4)

where C_ext is the wavelength-dependent extinction cross-section and $K_2={\int }_0^{\infty }(\text{d}{B}_{\lambda }/\text{d}T)\text{d}\lambda$ is a normalization factor.

For illustrative purposes and to stress the meaning of the optical depth ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$, Fig. 1 displays properties of reemitted photon packages and their experienced optical properties. The top plot displays the difference spectrum dB_λ/dT that is normalized to a maximum value of 1 for each temperature value T_dust. It shows that at higher temperatures, the peak of the spectrum moves toward shorter wavelengths. In particular, the maximum of this spectrum follows a law that is similar to Wien’s displacement law given by λ_max T ≃ 2410 μm K, see Appendix A.3 for a short derivation. Compared to Wien’s displacement law, the maximum of the difference spectrum is red-shifted, which is a necessary condition for the temperature correction technique for which it is used. Additionally, in the bottom plot, we show the temperature dependence of the effective extinction cross-section Ĉ_ext defined in Eq. (4). Starting from T_dust = 2.73 K, the cross-section Ĉ_ext increases strongly with the dust temperature, except for a local maximum with a subsequent dip, which is a consequence of a silicate feature¹. Since Ĉ_ext is proportional to the optical depth ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$, see Eq. (3), both quantities share the same temperature dependence. To illustrate the corresponding value of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ in this plot (right axis), we chose the arbitrary values n_dust = 7.8 × 10⁷ m⁻³ and r_sph = 1 au for the dust number density and the sphere radius, respectively. Altogether, both plots show that as the dust temperature rises in a radiative transfer simulation, the probability of reemitting photon packages with shorter wavelengths increases and, as a consequence, the effective optical depth, that is, the opaqueness of the medium, which is experienced by these photon packages, increases as well. Thus, in order to achieve a certain level of opaqueness at lower temperatures, the number density needs to be much higher compared to the case of high temperatures. Overall, the reason for using this temperature-dependent extinction optical depth ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ rather than an optical depth for a specific wavelength is that it is a much better measure for the spheres optical thickness with regard to the radiation that it produces at its particular temperature, that is, similar values of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ at different temperatures correspond to a similar computational effort.

Fig. 1. Visualization of the temperature dependence of reemitted photon packages when using the (difference) spectrum dBλ/dT. Top: difference spectrum dBλ/dT that is normalized for each temperature value Tdust to a maximum value of 1 is displayed. Bottom: effective extinction cross-section Ĉext (see Eq. (4); left axis) and the optical depth ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ (see Eq. (3); right axis) as a function of the dust temperature Tdust. For details, see Sect. 2.2.

We used 501 logarithmically sampled temperatures T ranging from 2.7 K to 3000 K and 132 wavelengths λ between 50 nm to 2 mm and kept the radius of the sphere fixed at the arbitrary value r_sph = 1 au. In Sect. 2.3 we show the scalability of our simulations with regard to the sphere radius and dust density.

2.2.2. Minimal missing absorption optical depth ${\tau }_{\mathsf{abs}}^{\mathsf{mis}}$

The radius of the last event of absorption r_la strongly depends on the optical depth of the sphere. Higher values of ${\widehat{\tau }}_{\mathrm{ext}}$ push the last location of absorption closer toward the rim of the sphere. We chose the minimal missing absorption optical depth ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$ that the photon package has to overcome in order to leave the sphere as a measure for r_la. It is given by

$$\begin{array}{c}\hfill {\tau }_{\mathrm{abs}}^{\mathrm{mis}}=(r_{\mathrm{sph}}-r_{\mathrm{la}})C_{\mathrm{abs}}{n}_{\mathrm{dust}}.\end{array}$$(5)

We note that C_abs depends on the wavelength λ the photon package has when escaping the sphere, which also corresponds to the wavelength it acquires at its last event of reemission. The only means of interaction between these two events is scattering, which leaves the wavelength of the photon package unaffected. We used 100 logarithmically sampled bins with a maximum value of ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}=30$. We note that the use of an analogously defined missing extinction optical depth ${\tau }_{\mathrm{ext}}^{\mathrm{mis}}$ instead of the missing absorption optical depth ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$ is conceivable; however, ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$ shows to be a reliable measure and we do not expect any advantages in using ${\tau }_{\mathrm{ext}}^{\mathrm{mis}}$.

In Fig. 2 the probability distributions (upper row) for three different temperatures 2.73 K, 45 K, and 1508 K and different values of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ are presented. The underlying probability distributions for ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$ show, in accordance with our expectation, striking similarities among all temperatures and spheres of sufficiently high values of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$. In the case of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}=10$, the resulting distributions differ the most. This effect stems from the wavelength-dependent optical properties of the dust. The last location of absorption can only be as deep as the sphere’s center, resulting in a strict upper limit for ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$, which falls below ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}=10$ in certain wavelength regimes and thus causes the differences seen in Fig. 2. After performing ∼10¹¹ simulations for different sphere sizes and temperatures, we find that not a single photon package originated from an optical depth of ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}\ge 30$. Since the probability distribution is a bin size dependent function, we calculated its corresponding probability density displayed in the second row of Fig. 2. Photon packages leaving the sphere have the highest probability of being absorbed and reemitted from an optical depth of about ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}=0.1$. The average optical depth, nonetheless, is close to ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}=1$. Due to the method of immediate reemission (Bjorkman & Wood 2001) and the fact that the photon package starts its path at the center of a sphere, high values of ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$ are favored. A photon package may only stem from a low value of ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$ if all previous reemission events from further inside the sphere failed to transport it through the rim of the sphere. However, the fact that a photon package may scatter counteracts this trend as the geometrical reach of a photon package from one event of emission to a subsequent event of absorption decreases.

Fig. 2. Result for the distribution ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$ using 106 simulations per temperature and a value of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$. Three different temperatures 2.73 K, 45 K, and 1508 K are displayed in the first, second, and third columns, respectively; the probability is shown in the first row and the probability density is illustrated in the second row. Different sphere sizes ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ exhibit notable similarities with regard to the shape and average value of their corresponding ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$ distributions. This shows that radiation, which is emitted by dense spheres of different sizes but with identical temperatures, thus, originates from a similar distance from the rim of each sphere.

2.2.3. Relative total absorbed energy X

The total absorbed energy of the sphere $\mathrm{\Delta }{E}_{\mathrm{abs}}^{\mathrm{tot}}$ due to a single photon package depends on the pathway segments the photon package travels along inside the sphere and is given by

$$\begin{array}{c}\hfill \mathrm{\Delta }{E}_{\mathrm{abs}}^{\mathrm{tot}}=E_{\gamma }{n}_{\mathrm{dust}}{\displaystyle \sum _i}{C}_{\mathrm{abs}}({\lambda }_i)l_i,\end{array}$$(6)

where λ_i denotes the wavelength of the photon package after its ith event of interaction and l_i is its corresponding pathway segment.

Figure 3 shows the average relative energy a photon package deposits in a sphere for different temperatures and sphere sizes. All depicted sphere sizes feature a similar behavior, especially toward higher values of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$. Furthermore, due to the temperature-dependent definition of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$, each curve is located at a distinct level. Denser spheres, that is spheres with a higher value of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$, generally result in higher absorbed energies and show a nearly constant shift in the log-log plot. In Fig. 4 we can see that the gray area, marking the region for possible values of $\mathrm{\Delta }{\overline{E}}_{\mathrm{abs}}^{\mathrm{tot}}/E_{\gamma }$, exhibits a clear dependence of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$, in particular,

$$\begin{array}{c}\hfill \mathrm{\Delta }{\overline{E}}_{\mathrm{abs}}^{\mathrm{tot}}\sim {\left({\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\right)}^2⟹\mathrm{\Delta }{\overline{E}}_{\mathrm{abs}}\sim n_{\mathrm{dust}}{r}_{\mathrm{sph}}^2.\end{array}$$(7)

Fig. 3. Average relative deposited energy using 104 simulations per ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ and temperature value.

Fig. 4. Extinction optical depth ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ dependence of the average total relative absorbed energy per photon package. The red line exemplarily shows the behavior for a T ≈ 1500 K sphere. The gray area marks the region for $\mathrm{\Delta }{\overline{E}}_{\mathrm{abs}}^{\mathrm{tot}}/E_{\gamma }$ found for the complete range of temperatures between 2.73 and 3000 K.

The latter relation in Eq. (7) describes the average energy a single dust grain in the sphere acquires through a single escaping photon package. In order to mitigate the problem of binning potentially arbitrary high values of $\mathrm{\Delta }{E}_{\mathrm{abs}}^{\mathrm{tot}}$, we instead binned the quantity X, which is defined as follows:

$$\begin{array}{c}\hfill X\equiv \frac{\mathrm{\Delta }{E}_{\mathrm{abs}}^{\mathrm{tot}}}{E_{\gamma }{\left({\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\right)}^2}=\frac{n_{\mathrm{dust}}}{{\left({\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\right)}^2}{\displaystyle \sum _i}{C}_{\mathrm{abs}}({\lambda }_i)l_i.\end{array}$$(8)

This particular definition of the quantity X serves a few purposes: (i) It is a measure for the total deposited energy in a sphere, (ii) it is independent of the photon energy E_γ (for a comparison, see Eq. (6)), and (iii) its factor of $1/{\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}{}^2$ results in a rather weak dependence of the sphere size (for a comparison, see Eq. (7)).

In Fig. 5 we show the result of performing 10⁶ simulations per temperature and sphere size for the normalized distribution of the quantity X. We used 400 logarithmically sampled bins for X ranging from 10⁻⁷ to 10 with an additional bin for X = 0. We find that in all cases, the probability density of X peaks strongly at around X ≈ 10⁻²…10⁻¹ and only in the case of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}=10$, the gray curve, does the distribution exhibit significant contributions from particularly low values of X. In addition to that, Fig. 5 shows a log-normal fit that we performed. We find that the distribution of X is well described by a log-normal distribution for all cases of high optical depth. Additionally, toward higher values of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$, the distributions show indications of convergence for each individual temperature. A remarkable common property among all distributions is the lack of a X >  10 contribution, that is, for all preformed simulations we find

$$\begin{array}{c}\hfill X<X^{\max }=10.\end{array}$$(9)

Fig. 5. Results for the distributions of the variable X using 106 simulations per temperature and ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ value. Outcomes for six different temperatures are displayed, which correspond to the resulting distributions (histograms) as well as their log-normal distribution fits (dashed lines) using corresponding colors. For details, see Sect. 2.2.

Overall, per calculated sphere size ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$, the PDF and the CDF each have a total of 501 × 400 × 100 × 132 ($\doteq M_T\times M_X\times M_{{\tau }_{\mathrm{abs}}^{\mathrm{mis}}}\times M_{\lambda }$ where M_x is the number of bins for a random variable x) entries that have to be stored and quickly accessed during an actual Monte Carlo radiative transfer simulation. Using single precision (4 byte per entry), this corresponds to ∼10 gigabytes per sphere size ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$. In realistic simulations, the dust density of cells differs and thus also the spatial extent of spheres with equal values of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ varies on a cell basis. To avoid costly precalculations for each individual cell, which would lead to arbitrarily high memory requirements, we made use of the scalability of our simulations, which is the topic of the following subsection. Additionally, we performed a series of further tests to make sure that the number of bins, their logarithmic scaling, and full range result in acceptable coverage of the specifics of all distributions of X and ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$. We find that this specific sampling (i) gives reliable results and good coverage for all distributions while (ii) also keeping the memory requirements reasonably low. In particular, firstly, using a linear sampling instead increases the error of the distributions averages, secondly, a higher number of bins does not justify the resulting memory consumption, and thirdly, a lower number of bins leads to a a higher error of the average value of distributions of both random variables X and ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$.

2.3. Scalability

In a Monte Carlo radiative transfer simulation, the path a photon package travels in a homogeneous environment, with nonzero density, is determined by and thus equivalent to a sequence of random numbers. The optical depth τ a photon package travels is a random variable that solely depends on one of these random numbers. During our precalculations, we kept the optical properties of the environment fixed. Therefore, the path length l that a photon package travels only depends on the random variable τ and the ambient dust number density n_dust according to l ∼ τ/n_dust. Since it does not appear elsewhere in the path determination, the number density of dust only effects the path by scaling its length. Thus, given a realization of the aforementioned sequence of random numbers, the derived quantity l ⋅ n_dust takes the same value independent of n_dust. Consequently, the precalculated distributions of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$, ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$, and X (see Eqs. (3), (5), and (8)) are also independent of n_dust. This is one crucial property that makes our method viable in the first place.

In Fig. 6 we show results of a scalability test for the distribution of X. For the calculation of the distribution, we used sphere sizes of r_sph = 1 au and adjusted the density such that ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ reached the specific values from Eq. (3). We then performed the same simulations, but we altered r_sph values by a factor of 100 and 0.01 while keeping ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ fixed and readjusting the number density n_dust. Figure 6 exemplarily shows the probability differences between simulations with a different sphere size and density for three different temperatures 2.73 K, 45 K, and 1508 K. The first row compares the simulations using r_sph = 100 au with the simulations using r_sph = 1 au. The second and third rows compare r_sph = 1 au with r_sph = 0.01 au and r_sph = 100 au with r_sph = 0.01 au, respectively. We find that the differences are on the order of ∼1%. Additionally, for every subplot and value of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$, we calculated the relative differences of the average value of X. The maximum value of the relative difference is displayed in each subplot and is also on the order of ∼1%. Considering the stochastic nature of this process and the fact that the number of counts per bin follows a Poisson distribution, this value equals the expected S/R of $1/\sqrt{{10}^4}=1\%$. The differently scaled models are thus in agreement with our expectation. As an additional test, we increased the number of simulated photon packages by two orders of magnitude and we obtain a higher level of agreement for the distributions of one order of magnitude, which can be seen in Fig. A.1. This shows that the error of the precalculated distributions behaves as expected and that scalability is given.

Fig. 6. For three different temperatures 2.73 K, 45 K, and 1508 K in the first, second, and third columns, respectively, show the probability differences between distributions of X using simulations with different sphere sizes rsph but equal values of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$. For every temperature and sphere size, 104 simulations were performed. Text boxes contain information about the compared sphere sizes. For example, “100–1” denotes the comparison of the distribution using a sphere size of rsph = 100 au with the distribution using a sphere size of rsph = 1 au. In addition to that, the relative difference of the averages of two compared distributions is given by δ⟨X⟩ in the upper left corner of each plot. Low values of δ⟨X⟩ indicate a higher level of agreement.

2.4. Method

During a Monte Carlo radiative transfer simulation, a photon package travels through a grid composed of cells with homogeneous physical properties. Succeeding an event of emission and given that the photon package is in a sufficiently large distance from the cells boundaries, we can make use of the results of our precalculations. For the highest performance of this method, we have to find the largest possible precalculated sphere that (i) fits into the cell and (ii) leaves the cell’s physical conditions unchanged when the photon is transfered from the sphere’s center to the last point of absorption. In order to meet the first (spatial) requirement, we calculated the minimal distance l_min to the next cell border. Then the corresponding optical depth ${\widehat{\tau }}_{\mathrm{S}}^{\max }$ was calculated as follows:

$$\begin{array}{c}\hfill {\widehat{\tau }}_{\mathrm{S}}^{\max }={\widehat{C}}_{\mathrm{ext}}{l}_{\min }{n}_{\mathrm{dust}}.\end{array}$$(10)

To meet the requirements, the chosen size of the sphere has to satisfy the condition ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\le {\widehat{\tau }}_{\mathrm{S}}^{\max }$. The second condition requires the constancy of all physical properties that effect the path of a photon package during its way to the sphere’s rim. The physical conditions inside the cell (optical properties, reemission spectrum, dust number density) are homogeneous but not necessarily constant since the cell temperature changes during the run of a simulation and as do its spectrum of reemitted photon packages and effective optical properties, see for example Eq. (4). However, the effective optical properties were precalculated for each of the 501 temperature values and they only change when reaching the next temperature bin. Therefore, there is a maximum energy $\mathrm{\Delta }{E}_{\mathrm{T}}^{\uparrow }$ that can be deposited in a single dust grain before the next temperature bin is reached. Using Eqs. (1) and (8), the second (temperature) condition thus requires

$$\begin{array}{c}\hfill {\left({\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\right)}^2\stackrel{!}{<}{V}_{\mathrm{cell}}{n}_{\mathrm{dust}}\frac{\mathrm{\Delta }{E}_{\mathrm{T}}^{\uparrow }}{E_{\gamma }}\frac{1}{X}·\end{array}$$(11)

In making use of the result from Eq. (9), that is the value of X has an upper limit, this requirement implies a maximum optical depth ${\widehat{\tau }}_{\mathrm{T}}^{\max }$ that must not be exceeded by a chosen sphere in order to leave the cell’s optical properties unaffected. The maximum optical depth ${\widehat{\tau }}_{\mathrm{T}}^{\max }$ is thus given by

$$\begin{array}{c}\hfill {\widehat{\tau }}_{\mathrm{T}}^{\max }={\left(\frac{V_{\mathrm{cell}}{n}_{\mathrm{dust}}}{X^{\max }}\frac{\mathrm{\Delta }{E}_{\mathrm{T}}^{\uparrow }}{E_{\gamma }}\right)}^{1/2}·\end{array}$$(12)

Therefore, for the highest efficiency possible, we chose a precalculated sphere size ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ given as follows:

$$\begin{array}{c}\hfill {\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}=\underset{\tau \in {\mathcal{T}}_{\mathrm{sph}}}{max}\{\tau |\tau \le min({\widehat{\tau }}_{\mathrm{S}}^{\max },{\widehat{\tau }}_{\mathrm{T}}^{\max })\}.\end{array}$$(13)

If no such sphere has been calculated, it either means that the photon package is too close to a cell boundary or that the next temperature level is close. This implies that the method provides the greatest performance boost in the case of optically thick cells, that is, when the efficiency of the basic Monte Carlo method is low.

Given a value for ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ and the cells temperature, we used the precalculated multidimensional CDF to randomly choose values for X, ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$, and λ. The value of X was used to immediately increase the deposited energy of every dust grain in the corresponding cell. Using ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ and ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$, the distance r_la can be determined. Since this method was used after an event of isotropic reemission occurs, the direction in which the photon package is displaced was chosen isotropically. During the precalculation, the quantity r_la marks the last point of absorption before the photon package left the sphere. Therefore, we had to launch the photon package at that position and transport it through the sphere’s rim without another absorption on its way. In a real simulation, this is done by basic Monte Carlo transfer. This process has to be restarted until the photon package leaves the sphere without being absorbed. Only after a photon package is successful can the transfer of it continue. If the optical density of the cell is high, this method can be used multiple times to displace the photon package before eventually leaving the cell. The most time-consuming part of the method is its hit-or-miss nature. Even though it is already launched close to the rim, it can take many attempts before the photon package manages to leave it. Nonetheless, every time the method can be used prevents having to calculate a potentially high number of interactions and thereby computation time. If it were possible to remove either of the method’s limitations, including grid induced or temperature induced limits, it would greatly benefit the effectiveness of our method. The grid induced restrictions are difficult to overcome without the loss of spatial resolution of the simulation. The temperature induced limitation, though, may potentially be removed by adjusting the reemission spectrum (see Eq. (2)) of previously absorbed photon packages during our precalculations. However, this is not part of this paper and may very well be part of a follow-up study in the future.

2.5. Optimization – performance boosts

One step of this method involves the absorption-free transfer of a photon package from its last location of absorption through the rim of the sphere. If the photon package is absorbed before leaving the sphere, it has to be reset to the previous location of reemission. Two major issues arise that significantly lower the chance of a successful escape of the photon package: (i) the optical depth and (ii) the initial direction of the photon package. The probability distribution of the distance a photon package travels before its absorption is given by an exponential function. It becomes increasingly difficult to transport a photon package over a large distance, which can lead to a high number of failed attempts. The direction of emission of the photon package is chosen randomly according to an isotropic distribution. The chance for a successful escape, however, may strongly depend on the direction of emission. Thus, a large portion of photons are initially sent in a direction with a very low escape rate. In general, it is crucial to implement optimization techniques that alleviate this problem. In fact, the viability of our method is largely grounded in their utilization. However, there are simple solutions that help to speed up this costly procedure significantly. In the following, we present three different techniques that serve this purpose.

2.5.1. Optical depth problem – splitting scattering and absorption length determination

A photon package with wavelength λ is emitted and follows a path along which it may encounter several points of interaction. The type of interaction (absorption or scattering) is determined randomly by using the albedo A = C_sca/C_ext. Thus, if A ≠ 1, it may happen at any point of interaction during the photons path that it is absorbed and does not escape the sphere successfully. Therefore, a large portion of photons may not even have a chance to escape the sphere as their traveled distance before absorption does not reach the minimum distance between the sphere’s rim and the location of emission. In order to avoid sending out these photon packages, we split the determination of scattering and absorption events. Instead of using the extinction cross-section C_ext as a measure for the interaction rate and subsequently deciding the type of interaction, we used the absorption and scattering cross-sections C_abs and C_sca, respectively, and determined the traveled distances that lead to events of absorption and scattering independently. A short proof for the viability of this approach can be found in Appendix A.1.

During the simulation, the photon package was emitted inside the sphere. Its location of last absorption was determined according to its absorption optical depth ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$. Subsequently, another absorption optical depth τ_abs was randomly chosen according to the absorption cross-section and dust number density. Escaping photon packages with ${\tau }_{\mathrm{abs}}<{\tau }_{\mathrm{abs}}^{\mathrm{mis}}$ can immediately be discarded. Thus, rather then limiting the probability density of τ_abs to values $\ge {\tau }_{\mathrm{abs}}^{\mathrm{mis}}$, we instead used the sum of ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$ and τ_abs as a measure for the propagated path length l_abs. In doing so, no photon package was created at the point of last absorption, which does not have a chance of leaving the sphere successfully. Until traversing the distance l_abs, the only means of interaction is scattering. If the photon package reaches l_abs before leaving the sphere, it has to be reset. Additionally, whenever a photon package arrives at a scattering point and before it scatters, it is determined whether the photon package can potentially reach the sphere’s rim after the forthcoming scattering event. If its remaining path length before absorption does not suffice in overcoming the remaining minimal distance to the rim of the sphere, the photon package is immediately reset to its point of emission. The importance of this step becomes apparent when considering a high optical depth of, for example, ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}=3$. In this case, the probability for ${\tau }_{\mathrm{abs}}<{\tau }_{\mathrm{abs}}^{\mathrm{mis}}$, which are certainly unsuccessful attempts, is ≳95%. This step of optimization is crucial for our method.

2.5.2. Initial direction problem – precalculated direction distribution

When the photon package inside the sphere is reemitted after an absorption event, the direction of emission follows an isotropic distribution. If the position of the last location of emission r_la coincides with the center of the sphere described by r_c, the escape probability p_esc for the photon package, which generally depends on the direction emission, is the same in any direction. If this is not the case, that is, if r_la ≠ r_c, the setup is axially symmetric along the axis defined by r_la − r_c. Thus, the direction-dependence of p_esc reduces to a dependence of the launching angle θ alone. The launching angle θ is the angle about which the direction of emission differs from the direction of the shortest distance to the rim of the sphere and is defined by the following equation:

$$\begin{array}{c}\hfill cos(\theta )=\frac{({\mathit{r}}_{\mathrm{la}}-{\mathit{r}}_{\mathrm{c}})·{\mathit{v}}_{\gamma }}{\Vert {\mathit{r}}_{\mathrm{la}}-{\mathit{r}}_{\mathrm{c}}\Vert \Vert {\mathit{v}}_{\gamma }\Vert },\end{array}$$(14)

where v_γ is the direction of emission of the photon package. When launching a photon package close to the rim of a sphere with a relatively high value of ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$, the probability may peak strongly at θ = 0 since shorter distances have a higher chance of being traversed without absorption. When randomly choosing θ according to an isotropic distribution, a large portion of launched photon packages therefore have low escape probabilities. This problem can be solved by precalculating the probability distribution for a successful escape for every possible sphere size, wavelength, and value ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$. We used 181 bins for the launching angle θ, ranging from 0 (forward emission) to π (backward emission), and set up simulations to determine the probability distribution. Compared to the size of the already precalculated CDF, which is described in Sect. 2.2, this CDF increases the memory consumption by ∼0.1%. In order to calculate the CDF properly, it is important to consider that the escape probability for every θ-bin has to be weighted with the emission probability of a photon in this bin. In the case of isotropic emission, the PDF for cos(θ) is constant.

During the radiative transfer simulation, the above distribution is used whenever the photon package is emitted at its last point of absorption r_la. It is important to keep the randomly chosen value of θ fixed during all escape attempts. In order to avoid performance drops, we restarted this process with newly assigned random values when a high number of attempts was exceeded. Since this number was barely reached, we were not able to find a measurable impact of this additional parameter on the resulting temperature as well as on the emission properties of the sphere. Nontheless, the distributions that are obtained from precalculated spheres are completely generated without it.

2.5.3. Infinite beam splitting

Whenever our method is used and a photon package is transported from the center of the sphere r_c to its last location of absorption r_la, a value for X is randomy chosen. By doing so, a single randomly chosen value for the amount of deposited energy is determined and increases the deposited energy of the corresponding cell. Instead of choosing a random value for X, it is possible to choose the mean ⟨X⟩ instead. This corresponds to a scenario in which the photon package is split at the center of the sphere into an infinite number of photon packages, which all follow their individual paths through the sphere until they eventually meet and combine at the location of last absorption. Since we assume that Monte Carlo radiative transfer simulations converge to the proper solution when increasing the number of photon packages, this approach decreases the level of noise of a Monte Carlo simulation. While the computation speed is not necessarily boosted significantly, the memory requirements for the results of our precalculations decrease strongly, in our case, by a factor of 400.

2.6. Performance test

In this section, we showcase the performance of this method and compare it to the basic Monte Carlo method, see Sect. 2.1, before we eventually apply it to the exemplary case of a circumstellar disk in Sect. 3. The testbed that we use consists of a dense isothermal sphere with a radius of r_sph = 1 au and a constant density corresponding to an effective extinction optical depth of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}={10}^{3.5}$ at a temperature of T = 1500 K. The bottom plot of Fig. 1 shows the resulting temperature dependence of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$. We chose this particular setup for the following reasons. Firstly, the value of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}={10}^{3.5}$ exceeds the optical depth of our precalculated spheres (see Eq. (3)). Secondly, in a realistic simulation, the precalculated sphere sizes ideally cover the full range of potentially encountered optical densities. Thus, this value of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ only exceeds the largest precalculated sphere size by half an order of magnitude.

The initial temperature of the dense sphere is 2.73 K. Throughout the simulation, photon packages were emitted isotropically at the center of the sphere. The initial wavelengths were randomly chosen according to the latest temperature during the simulation run. Launched photon packages move inside the dense sphere until they penetrate its rim. Since the temperature of the sphere is low at the start of the simulation, the spectrum of emitted photons initially peaks at a relatively long wavelength and thus these photon packages encounter a rather low dust extinction cross-section (for a comparison, see Fig. 1). While the temperature increases with the number of photon packages that were created, the dense sphere becomes increasingly opaque for the newly created photon packages on average. Eventually, it reaches its maximum effective optical depth at a temperature of 1500 K. During the whole process, we tracked the properties of escaping photon packages in terms of X, ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$, and λ as we did with the precalculated spheres. Additionally, we tracked the number of photon packages needed to heat up the dense sphere from its initial temperature to any temperature up to 1500 K. The photon package energy E_γ together with the tracked number of photon packages corresponds to a central heating luminosity, which is required to sustain the corresponding sphere temperature. This simulation was performed using the basic Monte Carlo method with E_γ = 10⁻⁵ as well as by using our method for highly optically thick regions with E_γ ∈ {10⁻⁵,10⁻⁷,10⁻⁹}, where E_γ is given in units of L_⊙s. For the latter case, we made use of the performance boosts described in Sects. 2.5.1 and 2.5.2, involving absorption-scattering splitting and an improved initial direction, respectively. Results of our comparison are displayed in Fig. 7.

Fig. 7. Viability and performance test – comparison between the basic Monte Carlo method and our method. Setup: a dense sphere with a radial effective extinction optical depth ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}(T=1500\mathrm{K})={10}^{3.5}$ is heated by a central source from 2.73 K to 1500 K. Upper left: heating curves of the sphere show high level of agreement. Upper right: performance comparison. Bottom: emitted photon statistic shows high level of agreement.

We find that the heating curves in the upper left plot of Fig. 7 greatly overlap. The required heating power to reach the final temperature among all four different simulations is in good agreement and only differs by <0.1%. Generally, relative differences between the heating curves at arbitrary temperatures are within the expected level of inherent noise due to the use of the Monte Carlo method. The three lower plots show the distributions of the tracked quantities X, ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$, and λ from left to right, respectively. Distributions with higher values for E_γ require a lower number of emitted photon packages to reach any temperature level, which in turn results in a higher level of noise. Nonetheless, it becomes apparent that for all tracked quantities, the distributions overlap, thus also showing agreement between our method and the basic Monte Carlo method. We note that the shape of these normalized distributions is strongly affected by the particular bin sizes that have been used. For the sampling of wavelengths, for example, we used a catalog which samples certain domains with higher density than others, resulting in this particular shape. These tests show agreement between the heating process as well as agreement between the spectrum of the spheres and the properties of emitted photon packages. We note that in both cases, we were able to simulate the full radiative transfer, including the changing polarization state of photon packages due to scattering. The only striking difference between both methods can be seen in the upper right plot of Fig. 7. This plot shows the number of photons that we were able to transport from the center of the dense sphere outward and through its rim per second of real computation time. Certainly, the specific numbers in this plot vary for different computational architectures; the overall trends, nonetheless, are representative. The orange curve corresponds to the basic Monte Carlo method and achieves an overall rate of emitted photon packages of ∼1 photon package per second. The three other curves correspond to simulations with our newly introduced method. We find that for these three cases, the overall rate of emitted photon packages reaches ∼10³ photon packages per second, thus a boost of three orders of magnitude. Except for the overall boost in performance, we find that this method has the positive effect of stabilizing the number of transferred photon packages during the whole run. As the temperature gradually increases, the performance of the basic method decreases. The reason for this effect is the overall increase in effective optical thickness of the sphere with temperature and its accompanying increased number of interactions. Our method, however, simply transfers the photon package close to a precalculated sphere’s rim, independent of the temperature. In other words, we transport the photon package at a consistent pace in terms of spatial units while the basic Monte Carlo method transports in terms of units of the extinction optical depth.

One interesting feature of this plot is the decreased efficiency at low temperatures for high values of E_γ. The reason for this drop is the fact that our method is limited in its choice of usable precalculated sphere sizes, in this case, by the criterion of constancy of the dense spheres optical properties. The higher the energy of the photon package, the more energy is deposited in the cell, thus, resulting in a quicker increase in temperature and change of optical properties of the sphere (for a comparison, see Eq. (12)). Since the temperature levels are sampled logarithmically, this leads to a drop of the efficiency at low temperatures. Nonetheless, the method outperforms the basic Monte Carlo approach significantly, especially toward higher temperatures when dust becomes increasingly opaque. Additionally, most of the computation time is usually spent at higher temperatures, that is when our method is at its peak performance. Another feature of these curves can be seen at about 800 K. Here we find a peak in performance for all three curves. This is due to the fact that when starting from this temperature, the effective extinction optical depth exceeds ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\approx {10}^3$ and thus the previously unutilized largest precalculated sphere, ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}={10}^3$, is used.

3. Application

In this section we apply our method to the situation of the radiative transfer in a viscously heated circumstellar disk and compare the results to a simulation based on the basic Monte Carlo method. Subsequently, we present a study of the impact of N_int on Monte Carlo radiative transfer simulations of systems with high optical depth.

3.1. Model description

We use the following three-dimensional parametrized circumstellar disk density distribution (Shakura & Sunyaev 1973)

$$\begin{array}{c}\hfill \rho (r,z)=\frac{\mathrm{\Sigma }(r)}{\sqrt{2\pi }h(r)}·exp[-\frac{1}{2}{\left(\frac{z}{h(r)}\right)}^2],\end{array}$$(15)

where r is the radial distance from the z-axis at the center and z is the spatial coordinate perpendicular to the midplane. The surface density distribution Σ(r) and density scale height h(r) are given by (Lynden-Bell & Pringle 1974; Hartmann et al. 1998)

$$\begin{array}{c}\hfill \mathrm{\Sigma }(r)={\mathrm{\Sigma }}_0{\left(\frac{r}{r_0}\right)}^{\beta -\alpha }·exp[-{\left(\frac{r}{r_0}\right)}^{2+\beta -\alpha }]\mathrm{and}h(r)=h_0{\left(\frac{r}{r_0}\right)}^{\beta },\end{array}$$(16)

where Σ₀ (h₀) is the surface density (scale height) at the reference radius r₀. The dimensionless parameters α and β are model parameters that shape the radial density profile and flaring of the disk, respectively. The total mass of the disk M_disk is obtained by the integration of Eq. (15) and is composed of gas and dust with a mass ratio of M_gas : M_dust = 100 : 1. In dense regions of a circumstellar disk, viscosity may become an important source of heating. Assuming a geometrically thin, optically thick accretion disk, the dissipated energy per unit area Ė_ν generated through viscosity is given by (Hartmann 1998)

$$\begin{array}{c}\hfill \frac{{\dot{E}}_{\nu }}{2}=\frac{3G{M}_{\ast }{\dot{M}}_{\ast }}{8\pi r^3}[1-{\left(\frac{R_{\ast }}{r}\right)}^{1/2}]\stackrel{R_{\ast }\ll r}{\approx }\frac{3G{M}_{\ast }{\dot{M}}_{\ast }}{8\pi r^3},\end{array}$$(17)

where R_*, M_*, and Ṁ are the stellar radius, mass, and its accretion rate, respectively, and G is the gravitational constant. In order to show the applicability of our method, we performed radiative transfer simulations of solely viscously heated disks, meaning that the radiation of the central star does not contribute.

We used a spherical grid with 300 logarithmically sampled cells in the r-direction and 91 linearly sampled cells in the θ-direction as well as one cell in the ϕ-direction. The disk has an inner radius r_in, which is defined by the sublimation of dust at its sublimation temperature T_subl. The outer radius r_out was chosen such that the disk extents to optically thin regions that are far from the star where the density has dropped significantly.

In our simulation, dissipated energy originates purely from the midplane cell layer and is determined by integrating Eq. (17) over the intersection of every cell with the midplane. For most model parameters, we chose commonly observed values that are typical of circumstellar disks around T Tauri stars (Galli et al. 2015; Andrews et al. 2009), for the specific values see Table 1. To generate a model with efficient viscous heating, we chose a rather high value for the dust mass M_dust. This model results in an overall dissipated luminosity of ∼1.58 × 10⁻³ L_⊙. In our model, the sources of radiation are dust grains, whose spectrum is temperature-dependent. We assume spherical dust grains composed of 62.5% silicate and 37.5% graphite with a ratio of parallel : ortho = 1 : 2 (Draine & Malhotra 1993) and use their corresponding mixed optical properties (Draine & Lee 1984; Laor & Draine 1993; Weingartner & Draine 2001). Dust grains have a bulk density of ρ_bulk = 2.5 g cm⁻³ and radii a, ranging from 5 nm to 250 nm following a grain size distribution given by dn ∼ a^−3.5da (Mathis et al. 1977). Additionally, we make use of mean optical properties of dust grains (Wolf 2003). The process for scattering is described by the Mie theory.

Finally, two important parameters that limit the accuracy and have a significant effect on the computation time are the number of photon packages N_γ and the maximum number of interactions N_int. The first parameter determines the energy of a single photon package E_γ = L_*/N_γ Δt for the stellar radiation. Increasing the number of photon packages N_γ leads to a lower level of noise in the temperature distribution of the circumstellar disk, but this comes with an increase in computation time. The parameter N_int limits the number of interactions a single photon package can encounter before being removed from the model space. In the best case, this number is set to infinity, thus, having no undesired impact on the simulation. We note that the code was designed to ensure energy conservation and that, in the case of a viscously heated disk, every dissipating cell emits an integer number of photon packages. Consequently, the energy per photon package depends on the emission rate of the cell from which it originates and is chosen such that it is closest to the previously mentioned value, that is, E_γ ≲ L_*/N_γ Δt. The way we treat radiation that is produced through viscosity is by first determining the required dust temperature of each cell in order to provide its specific luminosity. We use that temperature as an initial cell temperature. During the run of a simulation, emitted photon packages heat up the disk. During the simulation, photon packages are created according to the gradually increasing latest temperature of their corresponding cells.

3.2. Results

We performed simulations of the viscous circumstellar disk introduced in Sect. 3.1 using our method, while imposing no restriction on the maximum number of interactions, that is, N_int → ∞ and N_γ = 10⁸ photon packages. The left plot of Fig. 8 shows the resulting vertical cut of the temperature distribution through the midplane. In this case, within the first 0.1 au, multiple cells reach temperatures above the sublimation limit. As a consequence, the dust sublimates and the inner cavity size increases. As can be seen, the midplane is heated stronger than regions in the disk at higher scale heights. Looking at even higher scale heights, however, the circumstellar disk becomes optically thin and is heated by the edge of the disk’s innermost region, which, overall, results in the specific temperature distribution seen in Fig. 8. This setup is symmetric with regard to the midplane with differences that arise, inherently, due to the stochastic nature of the Monte Carlo method.

Fig. 8. Comparison of the resulting temperature distributions of a viscously heated circumstellar disk. Left: temperature distribution for Nγ = 108 and Nint → ∞ using our method. Right: relative difference between the basic Monte Carlo method and our method using Nint = 109.

Since the basic Monte Carlo method relies on N_int <  ∞ for many problems of high optical density and in order to ensure comparability, we limit our method by using the same value of N_int = 10⁹ as well. To do that, we have to use an estimator for the number of interactions undergone by a photon package whenever we make use of our method, that is, when transferring a photon package to the rim of a precalculated sphere. This number is not tracked during the precalculation since it increases the required memory consumption. Furthermore, its only purpose is to make sure that the number of interactions of a photon package does not exceed the limiting value N_int. In Appendix A.2, we derive an estimator for the number of interactions that we use for the comparison between both methods. The results of these simulations are shown in the right plot of Fig. 8. The plotted quantity |δ_T| is the relative difference between the temperature distribution obtained by using our method T_MC+ and by using the basic Monte Carlo method T_MC. It is defined by

$$\begin{array}{c}\hfill {\delta }_T=\frac{T_{\mathrm{MC}+}-T_{\mathrm{MC}}}{T_{\mathrm{MC}+}+T_{\mathrm{MC}}}·\end{array}$$(18)

We find that the difference between both methods is on the order of ∼2% for the largest part of the plot, which is represented by dark red colors. The highest differences of ∼5%, corresponding to bright yellow regions in the plot, are only found in single cells and can be explained by a low photon count in these cells. We thus find that the results of both methods are in agreement with each other. For this specific setup, we achieved a boost in computation speed by about one order of magnitude.

3.3. Impact of imposing a maximum interaction number limit

In this section we study the impact of the parameter N_int on the resulting temperature distribution of a circumstellar disk, as described in Sect. 3.1. Due to its reliability and the achieved boost in computation speed, the following simulations were performed using our newly developed method. In Fig. 9 we see the results of a comparison between a reference simulation with N_int → ∞ and N_γ = 10⁸ with simulations of a finite number of N_int. In the left (right) plot, we limited the number of interactions to N_int = 10⁷ (N_int = 10⁹) and plotted the relative difference to the reference model. Positive values of δ_T correspond to bright regions in the plots where the simulation with N_int → ∞ reached higher temperatures. We find that the differences are larger for the lower limit N_int = 10⁷, which is in accordance with our expectation. If the number of interactions of a photon package reaches its maximum (N_int), the photon package was removed. This reduces the number of photons emitted from that cell, thus, leading to an underestimation of its radiation field and consequently its dust temperature. The effect of this reduced temperature is most prominent in the midplane, resulting in relative differences of up to ∼1/3. Regions in the plot with negative values, that is, higher resulting temperatures for models with a lowered value of N_int, are rare. The highest differences greatly exceed the level of noise seen in the right plot of Fig. 8. The right plot of Fig. 9, which has a much higher limit of N_int = 10⁹, shows the greatest differences in cells that are close to the star, which are slightly off the midplane cell layer. These differences stem from the fact that a high number of photon packages are created in the innermost part of the disk and, in the case of a limited number of interactions, a part of them are removed from the model space without entering the regions off the midplane. On the contrary, in the reference model, these regions are heated up partly with photon packages that have undergone a number of interactions that exceed N_int. Nonetheless, it can be seen that the differences in the largest part of the circumstellar disk are compareable to the noise level (see Fig. 8) and only exceed this level in the innermost region of the disk.

Fig. 9. Studying the impact of Nint. Comparison of temperature distributions of viscously heated circumstellar disks. The reference model uses Nint → ∞. All displayed results stem from simulations using our method. Left: comparison with a model using Nint = 107. Right: comparison with a model using Nint = 109.

To discuss the consequences and highlight the significance of our results, Fig. 10 shows the radial midplane temperature profile of these disks. This layer of cells produces radiation due to viscous heating. We show results for our reference model with N_int → ∞ and N_γ = 10⁸ as well as for three models, each of which has one modified parameter. We find that our reference model as well as the model with a reduced number of photon packages N_γ = 10⁷ and the model with a reduced value of N_int = 10⁹ all produce comparable results for the midplane temperature profile. Only in the case of a further restricted value of N_int = 10⁷ do we find strong deviations from all other simulations. We identify two major differences: (i) The innermost temperature profile and the sublimation radius are both underestimated and (ii) the temperature of the outer part of the midplane is strongly affected by the radiation originating from the inner part of the disk. The first difference occurs since the density in this part of the disk is high and photon packages may reach their interaction limit even before leaving the cell, which results in an underestimation of the temperature in this region. Consequently, a finite value of N_int results in an underestimation of the inner cavity size. This effect may change simulated near-infrared long baseline interferometry observations, for instance, by MATISSE and Pioneer, and therefore impacts the interpretation of protoplanetary disk observations. The second temperature difference seen at higher radii also arises due to the removal of photon packages originating from the inner region of the disk. In the case of stellar radiation, outer regions are effectively shielded from direct irradiation. On the contrary, a viscously heated disk generates radiation throughout the extent of the disk. Thus, regions with strong viscous dissipation are shielded less efficiently form viscous dissipation than they are from stellar radiation. Due to a combination of relatively high luminosity and reduced shielding, photon packages from warm regions of the midplane reach further out laying regions and lead to a significant increase in the temperature. A limit on the number of interactions therefore underestimates the resulting midplane temperature, even at higher radii. This effect may very well be crucial, for example, in the case of simulated millimeter and submillimeter observations (e.g., by ALMA) and impact our understanding and interpretation of real observations. It is generally difficult to estimate the impact of photon packages that are removed too early and our study shows that the impact of these photon packages may very well extend to regions of a circumstellar disk, which are far away from the optically thick region from which they were emitted.

Fig. 10. Impact of Nint and Nγ on the midplane temperature of a viscously heated circumstellar disk.

An investigation of the impact of stellar irradiation on the resulting temperature distribution shows a rather weak dependence of N_int, in particular, when choosing N_int = 10⁷ results in a midplane temperature profile that is already comparable to the case of N_int → ∞, independent of whether the basic method or our new method was used. Furthermore, the level of noise in the optically thick regions is dominated by the inherent noise of a Monte Carlo simulation. In order to reduce the level of noise, an increased number of emitted photon packages N_γ is required. However, the midplane temperature in our specific setup is significantly dominated by heating through viscous dissipation. Therefore, a proper treatment of optically thick regions is especially relevant when dealing with radiation that originates from deep inside these regions, such as in the case of viscous dissipation or young accreting protoplanets.

4. Summary

In this work, we present a newly developed method to solve the problem of high optical depths in Monte Carlo radiative transfer simulations. More precisely, this method reduces the computation time with which a photon package is transported through a dense dusty medium in order to calculate the temperature distribution. We assumed homogeneously distributed dust grains and used the radiative transfer code Mol3D as a testbed for this method. We made use of precalculated results from simulations of photon packages that travel from the center of a sphere to its rim. During the precalculation, we performed Monte Carlo radiative transfer simulations to determine properties of the emitted radiation originating from a homogeneous sphere in dependence of its temperature and size. We track the relative total deposited energy in the sphere (X), the last location of absorption (using ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$), as well as the wavelength of emitted photon packages (λ). The results are binned and stored in a multidimensional CDF. The CDF is then used to quickly transport photon packages over long distances in a probabilistic manner. We note that the precalculation was performed using basic radiative transfer, which includes wavelength-dependent optical properties, a temperature-dependent dust emission spectrum, as well as Mie-scattering with a proper treatment of the polariazation state of the photon package. Therefore, the results of these simulations are unbiased. Next, we provide a summary of our findings as well as key features of our method.

1. Precalculations were performed for a set of discrete values of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$, see Eq. (3). Apart from the wavelength λ, we chose to track the variables X and ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$, which depend on the product of the sphere radius r_sph and dust density n_dust and not on any one of these variables separately. In Sect. 2.3, we show that the results for distributions of any value of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$, as a consequence, can be scaled and adjusted to the problem at hand, making this method viable in the first place.

2. We find that the tracked distributions show hints for convergence for increasing ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$. Interestingly, this leads to a simple relation for the mean deposited energy per dust grain and simulated photon package: $\Delta E_{\text{abs}}~n_{\text{dust}}{r}_{\text{sph}}^2$. Furthermore, we find that among all of the performed simulations, the probability for X >  10 is very low. In fact, not even a single transported photon package among all of the performed simulations reached that value, see Sect. 2.2.3. This serves as a great constraint for the maximally deposited energy per photon package, which is dependent on the temperature, density, and size of the spheres.

3. In Sect. 2.4, we derived and discuss the limitations of our method with regard to the change of the medium’s optical properties as well as spatial limits due to a model grid, see Eq. (13). Subsequently, in Sect. 2.5 we also present upgrades to the method that significantly boost the effectiveness – the splitting of scattering and absorption length (see Sect. 2.5.1), the utilization of our prior knowledge of the initial direction (see Sect. 2.5.2), as well as the infinite beam splitting (see Sect. 2.5.3).

4. We carried out a performance test using, as a testbed, a dense sphere of optical depth ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\approx 3162$ and achieve a boost of about three orders of magnitude in computational speed, see upper right plot in Fig. 7. In even denser environments, for instance, around young planets or stars, this value is easily exceeded and our method is expected to perform even better than that.

5. Simulations for increasingly high values of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ are time-consuming when using the basic Monte Carlo method. In the bottom plots of Fig. 7, we show that the statistic of emitted photon packages as well as the statistic of the deposited energy of the basic method and our method are in remarkable agreement. Therefore, the calculation of larger spheres can be sped up significantly when using the set of already precalculated smaller spheres. Its error has been shown to be in agreement with the variance of a Poisson distribution.

6. We applied our method to a simple model of an optically thick viscously heated circumstellar disk and show the validity of the method. We study the impact that the imposed maximum number of interactions N_int has on the resulting temperature profile of the disk. We find that the choice of N_int is decisive throughout the full extent of the disk and not only in the optically thick regions where it is generated, see Sect. 3.3. The impact of a value of N_int that is too low is generally difficult to predict.

High optical depths pose a well known problem to Monte Carlo radiative transfer simulations. This study presents a novel method that solves this problem of Monte Carlo radiative transfer simulations and thus enables us to properly predict temperatures of objects with optically thick regions.

Acknowledgments

We thank all the members of the Astrophysics Department Kiel for helpful discussions and remarks. We acknowledge the support of the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets (WO 857/17-1)”.

References

Appendix A: Derivation of expected optical depths

A.1. Interactions and absorption-scattering

The following derivation shows that the determination of absorption lengths and scattering lengths can be done independent of each other. The use of extinction lengths as well as interactions with a randomly assigned type of interaction (absorption or scattering) is equivalent to that.

A.1.1. Splitting absorption and scattering

The probability density for an absorption optical depth τ_abs is given by f(τ_abs) = exp(−τ_abs). The probability density for reaching a distance interval l is then given by

$$\begin{array}{c}\hfill p(l)=\underset{\mathrm{\Delta }l\to 0}{lim}\frac{P(l\in [l,l+\mathrm{\Delta }l])}{\mathrm{\Delta }l}={\displaystyle {\int }_{{\stackrel{\sim }{C}}_{\mathrm{abs}}l}^{{\stackrel{\sim }{C}}_{\mathrm{abs}}(l+\mathrm{\Delta }l)}}f({\tau }_{\mathrm{abs}})\mathrm{d}{\tau }_{\mathrm{abs}}.\end{array}$$

By expanding f(τ_abs) at τ_abs = C̃_absl where C̃_abs is the specific cross-section for absorption, we find that

$$\begin{array}{c}\hfill p(l)={\stackrel{\sim }{C}}_{\mathrm{abs}}exp(-{\stackrel{\sim }{C}}_{\mathrm{abs}}l).\end{array}$$(A.1)

A.1.2. Using interactions

The probability density for an interaction optical depth τ_ext is given by f(τ_ext) = exp(−τ_ext). Subsequent to a reemission event, a photon package is moved along a straight line until it reaches a randomly selected extinction optical depth. At that point, another random number determines whether the subsequent interaction is an event of absorption or scattering. The probability for a scattering event is given by C̃_sca/C̃_ext where C̃_sca (${\stackrel{\sim }{C}}_{\mathrm{ext}}$) is the specific cross-section for scattering (extinction).

In the scattering free case, the probability density for traveling the distance l after traversing the extinction optical depth τ_ext = C̃_extl is, in analogy to the derivation of Eq. (A.1), given by C̃_extexp(−C̃_extl. For the probability density of a subsequent event of absorption, we thus find

$$\begin{array}{c}\hfill p_0={\stackrel{\sim }{C}}_{\mathrm{ext}}exp(-{\stackrel{\sim }{C}}_{\mathrm{ext}}l)·(1-\frac{{\stackrel{\sim }{C}}_{\mathrm{sca}}}{{\stackrel{\sim }{C}}_{\mathrm{ext}}})={\stackrel{\sim }{C}}_{\mathrm{abs}}exp(-{\stackrel{\sim }{C}}_{\mathrm{ext}}l).\end{array}$$(A.2)

In order to calculate the probability for n >  0 events of scattering before the eventual event of absorption at distance l, we have to integrate over the probability density of all possible paths of this kind. Let l_i ≥ 0 denote the path length the photon package travels before its ith event of scattering. It follows that after its nth and final scattering event, it travels a distance of $l_a=l-\sum _{i=1}^n{l}_i$ before being absorbed. For the ith event of scattering, the probability density is, in analogy to the derivation of Eq. (A.2), given by C̃_scaexp(−C̃_extl_i. Due to the fact that the individual path lengths are constrained by $\sum _{i=1}^n{l}_i\le l$, the probability density p_n is given by

$$\begin{array}{cc}\hfill p_n& ={\displaystyle {\int }_0^l}d{l}_n{\stackrel{\sim }{C}}_{\mathrm{sca}}exp(-{\stackrel{\sim }{C}}_{\mathrm{ext}}{l}_n){\displaystyle {\int }_0^{l-l_n}}d{l}_{n-1}{\stackrel{\sim }{C}}_{\mathrm{sca}}exp(-{\stackrel{\sim }{C}}_{\mathrm{ext}}{l}_{n-1})\cdots \hfill \\ \hfill & \cdots {\displaystyle {\int }_0^{l-{\sum }_{i=2}^n{l}_i}}d{l}_1{\stackrel{\sim }{C}}_{\mathrm{sca}}exp(-{\stackrel{\sim }{C}}_{\mathrm{ext}}{l}_1)·{\stackrel{\sim }{C}}_{\mathrm{abs}}exp(-{\stackrel{\sim }{C}}_{\mathrm{ext}}{l}_a).\hfill \end{array}$$(A.3)

In integrating Eq. (A.3), one arrives at

$$\begin{array}{c}\hfill p_n=\frac{l^n{\stackrel{\sim }{C}}_{\mathrm{sca}}^n}{n!}{\stackrel{\sim }{C}}_{\mathrm{abs}}exp(-{\stackrel{\sim }{C}}_{\mathrm{ext}}l).\end{array}$$

This implies, that with the highest probability, the photon package scatters n_max = ⌊l · C̃_sca⌋ times. Eventually, we have to take the sum of p_n over all possible numbers of scattering events and find for the probability density:

$$\begin{array}{cc}\hfill p(l)& ={\displaystyle \sum _{n=0}^{\infty }}{p}_n\hfill \\ \hfill & ={\displaystyle \sum _{n=0}^{\infty }}\frac{l^n{\stackrel{\sim }{C}}_{\mathrm{sca}}^n}{n!}{\stackrel{\sim }{C}}_{\mathrm{abs}}exp(-{\stackrel{\sim }{C}}_{\mathrm{ext}}l)\hfill \\ \hfill & ={\stackrel{\sim }{C}}_{\mathrm{abs}}exp(-{\stackrel{\sim }{C}}_{\mathrm{ext}}l){\displaystyle \sum _{n=0}^{\infty }}\frac{l^n{\stackrel{\sim }{C}}_{\mathrm{sca}}^n}{n!}\hfill \\ \hfill & ={\stackrel{\sim }{C}}_{\mathrm{abs}}exp(-{\stackrel{\sim }{C}}_{\mathrm{ext}}l)exp\left({\stackrel{\sim }{C}}_{\mathrm{sca}}l\right)\hfill \\ \hfill & ={\stackrel{\sim }{C}}_{\mathrm{abs}}exp(-{\stackrel{\sim }{C}}_{\mathrm{abs}}l).\hfill \end{array}$$

Since this result and the result obtained in Eq. (A.1) are the same, we can split the calculation of absorption and scattering length. ▫

A.2. Estimator for N_int

For comparison purposes between the basic Monte Carlo method and our newly developed method, we kept the maximum number of interactions N_int equal in both simulations. In order to keep track of the number of interactions when using the new method, we estimated it every time the method was used. It is possible to derive an exact estimate when assuming $⟨X⟩\stackrel{{\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\to \infty }{\to }{X}^{\ast }>0$ for the case of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\to \infty$, where ⟨ ⋅ ⟩ denotes the average of all possible paths of the photon package:

$$\begin{array}{c}\hfill ⟨X⟩\stackrel{\mathrm{eq}.8}{=}\langle \frac{n_{\mathrm{dust}}}{{\left({\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\right)}^2}{\displaystyle \sum _i}{C}_{\mathrm{abs}}({\lambda }_i)l_i\rangle =\frac{n_{\mathrm{dust}}}{{\left({\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\right)}^2}\langle {\displaystyle \sum _a}{C}_{\mathrm{abs}}({\lambda }_a)l_a\rangle .\end{array}$$

For the latter equation, we changed the sum over interaction path lengths to a sum over absorption path lengths. This is possible since the wavelength of the photon package does not change between two consecutive absorption events. Since ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\to \infty$ in the limit, any two consecutive absorption events are independent of each other, we can exchange the average of the sum by the sum of averages. Let ⟨⋅⟩_τabs, λ denote the average of all possible paths between two consecutive absorption events. It follows that

$$\begin{array}{c}\hfill ⟨X⟩=\frac{n_{\mathrm{dust}}}{{\left({\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\right)}^2}{\displaystyle \sum _a}{\langle C_{\mathrm{abs}}l\rangle }_{{\tau }_{\mathrm{abs}},\lambda }=\frac{N_{\mathrm{abs}}}{{\left({\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\right)}^2}{\langle {\tau }_{\mathrm{abs}}\rangle }_{{\tau }_{\mathrm{abs}},\lambda }.\end{array}$$

For the latter equation, we used τ_abs = n_dustC_absl and replaced the sum by the total number of absorption events N_abs. The average is simply ⟨τ_abs⟩_τabs, λ = 1. Therefore, we arrive at

$$\begin{array}{c}\hfill ⟨X⟩=\frac{N_{\mathrm{abs}}}{{\left({\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\right)}^2}\Rightarrow N_{\mathrm{abs}}=⟨X⟩{\left({\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\right)}^2.\end{array}$$(A.4)

Fig. A.1. For three different temperatures 2.73 K, 45 K, and 1508 K in the first, second, and third columns, respectively, show the probability differences between distributions of X using simulations with different sphere sizes rsph but equal values of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$. For every temperature and sphere size, 106 simulations were performed. Text boxes contain information about the compared sphere sizes. For example, “100–1” denotes the comparison of the distribution using a sphere size of rsph = 100 au with the distribution using a sphere size of rsph = 1 au. In addition to that, the relative difference of the averages of two compared distributions is given by δ⟨X⟩ in the upper left corner of each plot. Low values of δ⟨X⟩ indicate a higher level of agreement.

In this limiting case, the total number of interactions N_int = N_abs + N_sca converges to $N_{\mathrm{int}}\to {\tau }_{\mathrm{ext}}^{\mathrm{tot}}$ since one interaction occurs per unit of τ_ext on average. Analogously to the derivation of Eq. (A.4) and using C_ext = C_abs + C_sca, we thus arrive at

$$\begin{array}{c}\hfill N_{\mathrm{int}}=\langle {\displaystyle \sum _i}{n}_{\mathrm{dust}}{C}_{\mathrm{ext}}({\lambda }_i)l_i\rangle =N_{\mathrm{abs}}(1+{\langle \frac{C_{\mathrm{sca}}(\lambda )}{C_{\mathrm{abs}}(\lambda )}{\tau }_{\mathrm{abs}}\rangle }_{{\tau }_{\mathrm{abs}},\lambda }).\end{array}$$

Since λ and τ_abs are independent random variables, we can split the average into a product of two separate averages. Thus,

$$\begin{array}{c}\hfill {\langle \frac{C_{\mathrm{sca}}(\lambda )}{C_{\mathrm{abs}}(\lambda )}{\tau }_{\mathrm{abs}}\rangle }_{{\tau }_{\mathrm{abs}},\lambda }={\langle \frac{C_{\mathrm{sca}}(\lambda )}{C_{\mathrm{abs}}(\lambda )}\rangle }_{\lambda },\end{array}$$

where ⟨⋅⟩_λ denotes the average over all possible wavelengths using the corresponding temperature-dependent PDF. With Eq. (A.4), we find that

$$\begin{array}{c}\hfill N_{\mathrm{int}}=⟨X⟩{\left({\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\right)}^2{\langle \frac{C_{\mathrm{ext}}(\lambda )}{C_{\mathrm{abs}}(\lambda )}\rangle }_{\lambda }·\end{array}$$

Therefore, we estimate the number of undergone interactions during the transport of a photon package from the center of a sphere to its last location of absorption with

$$\begin{array}{c}\hfill \mathrm{\Delta }{N}_{\mathrm{int}}=\lceil X·{\left({\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}\right)}^2{\langle \frac{C_{\mathrm{ext}}(\lambda )}{C_{\mathrm{abs}}(\lambda )}\rangle }_{\lambda }\rceil ·\end{array}$$

A.3. Displacement law of the difference spectrum

In this section we study the difference spectrum dB_λ/dT introduced by Bjorkman & Wood (2001). We find that (i) dB_λ/dT is positive for all λ, T >  0, (ii) it is unimodal, that is, it has exactly one (global) maximum, and (iii) the maximum follows a law similar to Wien’s displacement law, that is, λ_maxT=const. The difference spectrum is given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}{B}_{\lambda }}{\mathrm{d}T}=\frac{2h^2{c}^3}{k_{\mathrm{B}}{T}^2{\lambda }^6}{e}^{\frac{\mathit{hc}}{k_{\mathrm{B}}T\lambda }}{\left(e^{\frac{\mathit{hc}}{k_{\mathrm{B}}T\lambda }}-1\right)}^{-2},\end{array}$$

which is strictly positive for all λ, T >  0. For a fixed temperature T, its derivative with respect to the wavelength λ satisfies

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}\lambda }\frac{\mathrm{d}{B}_{\lambda }}{\mathrm{d}T}=s(T,\lambda )(-6x-1+2\frac{e^{1/x}}{e^{1/x}-1}),\end{array}$$

where $x=\frac{k_{\mathrm{B}}T\lambda }{\mathit{hc}}$ and s(T, λ) > 0 for all λ, T >  0. Thus, an extremal point requires

$$\begin{array}{c}\hfill 6x+1=2\frac{e^{1/x}}{e^{1/x}-1}·\end{array}$$

This equation has exactly one positive solution x_max >  0 since (i) the left-hand side (l.h.s.) is smaller than the right-hand side (r.h.s.) as x → 0⁺ and (ii) the slope of the l.h.s (m_l.h.s.) is given by m_l.h.s. = 6, while the slope of the r.h.s (m_r.h.s.) satisfies 0 <  m_r.h.s. <  2 for all x >  0. Furthermore, we find that

$$\begin{array}{c}\hfill \frac{\mathrm{d}{B}_{\lambda }}{\mathrm{d}T}\stackrel{\lambda \to 0^{+}}{⟶}0\wedge \frac{\mathrm{d}{B}_{\lambda }}{\mathrm{d}T}\stackrel{\lambda \to \infty }{⟶}0,\end{array}$$

thus requiring the extremal point at x_max to be a maximum and, consequently, dB_λ/dT to be unimodal. A numerical calculation shows that the maximum is located at

$$\begin{array}{c}\hfill x_{\max }\simeq 0.1675.\end{array}$$

Therefore, for a given temperature T, the unique (global) maximum of the difference spectrum dB_λ/dT, which is located at the wavelength λ_max, satisfies:

$$\begin{array}{c}\hfill {\lambda }_{\max }T\simeq 2410\mu \mathrm{m}\mathrm{K}.\end{array}$$

All Tables

All Figures

Fig. 1. Visualization of the temperature dependence of reemitted photon packages when using the (difference) spectrum dBλ/dT. Top: difference spectrum dBλ/dT that is normalized for each temperature value Tdust to a maximum value of 1 is displayed. Bottom: effective extinction cross-section Ĉext (see Eq. (4); left axis) and the optical depth ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ (see Eq. (3); right axis) as a function of the dust temperature Tdust. For details, see Sect. 2.2.

In the text

Fig. 2. Result for the distribution ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$ using 106 simulations per temperature and a value of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$. Three different temperatures 2.73 K, 45 K, and 1508 K are displayed in the first, second, and third columns, respectively; the probability is shown in the first row and the probability density is illustrated in the second row. Different sphere sizes ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ exhibit notable similarities with regard to the shape and average value of their corresponding ${\tau }_{\mathrm{abs}}^{\mathrm{mis}}$ distributions. This shows that radiation, which is emitted by dense spheres of different sizes but with identical temperatures, thus, originates from a similar distance from the rim of each sphere.

In the text

	Fig. 3. Average relative deposited energy using 104 simulations per ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ and temperature value.
In the text

Fig. 4. Extinction optical depth ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ dependence of the average total relative absorbed energy per photon package. The red line exemplarily shows the behavior for a T ≈ 1500 K sphere. The gray area marks the region for $\mathrm{\Delta }{\overline{E}}_{\mathrm{abs}}^{\mathrm{tot}}/E_{\gamma }$ found for the complete range of temperatures between 2.73 and 3000 K.

In the text

	Fig. 5. Results for the distributions of the variable X using 106 simulations per temperature and ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$ value. Outcomes for six different temperatures are displayed, which correspond to the resulting distributions (histograms) as well as their log-normal distribution fits (dashed lines) using corresponding colors. For details, see Sect. 2.2.
In the text

Fig. 6. For three different temperatures 2.73 K, 45 K, and 1508 K in the first, second, and third columns, respectively, show the probability differences between distributions of X using simulations with different sphere sizes rsph but equal values of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$. For every temperature and sphere size, 104 simulations were performed. Text boxes contain information about the compared sphere sizes. For example, “100–1” denotes the comparison of the distribution using a sphere size of rsph = 100 au with the distribution using a sphere size of rsph = 1 au. In addition to that, the relative difference of the averages of two compared distributions is given by δ⟨X⟩ in the upper left corner of each plot. Low values of δ⟨X⟩ indicate a higher level of agreement.

In the text

Fig. 7. Viability and performance test – comparison between the basic Monte Carlo method and our method. Setup: a dense sphere with a radial effective extinction optical depth ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}(T=1500\mathrm{K})={10}^{3.5}$ is heated by a central source from 2.73 K to 1500 K. Upper left: heating curves of the sphere show high level of agreement. Upper right: performance comparison. Bottom: emitted photon statistic shows high level of agreement.

In the text

	Fig. 8. Comparison of the resulting temperature distributions of a viscously heated circumstellar disk. Left: temperature distribution for Nγ = 108 and Nint → ∞ using our method. Right: relative difference between the basic Monte Carlo method and our method using Nint = 109.
In the text

	Fig. 9. Studying the impact of Nint. Comparison of temperature distributions of viscously heated circumstellar disks. The reference model uses Nint → ∞. All displayed results stem from simulations using our method. Left: comparison with a model using Nint = 107. Right: comparison with a model using Nint = 109.
In the text

	Fig. 10. Impact of Nint and Nγ on the midplane temperature of a viscously heated circumstellar disk.
In the text

Fig. A.1. For three different temperatures 2.73 K, 45 K, and 1508 K in the first, second, and third columns, respectively, show the probability differences between distributions of X using simulations with different sphere sizes rsph but equal values of ${\widehat{\tau }}_{\mathrm{ext}}^{\mathrm{sph}}$. For every temperature and sphere size, 106 simulations were performed. Text boxes contain information about the compared sphere sizes. For example, “100–1” denotes the comparison of the distribution using a sphere size of rsph = 100 au with the distribution using a sphere size of rsph = 1 au. In addition to that, the relative difference of the averages of two compared distributions is given by δ⟨X⟩ in the upper left corner of each plot. Low values of δ⟨X⟩ indicate a higher level of agreement.

In the text