© M. Baes et al. 2022

1 Introduction

Radiative transfer is the general process that describes the interaction and transfer of energy when electromagnetic radiation passes through a medium. In astronomy, where we need to extract virtually all of the information on the Universe from radiation, radiative transfer is relevant in every field, from exoplanets to cosmology. General reviews on astrophysical radiative transfer include the seminal works by Chandrasekhar (1960) and Rybicki & Lightman (1979).

The Monte Carlo method is probably the most widely used approach to deal with radiative transfer problems (e.g., Whitney 2011; Steinacker et al. 2013; Noebauer & Sim 2019). The essence of the Monte Carlo radiative transfer (MCRT) technique is to represent the radiation field as the flow of a huge number of individual photon packets. The life cycle of each photon packet is followed individually through the medium, whereby each step in this cycle is determined in a probabilistic way by generating random numbers from the appropriate probability density function (PDF). The radiation field is constructed from the statistical analysis of the trajectory of all the different photon packets. Important advantages of the MCRT method are its flexibility and its inherent 3D nature, which makes the method particularly attractive for complex geometries.

The physical processes of emission, absorption, and scattering are easily incorporated in the MCRT framework. These processes are the most important ones for continuum radiative transfer, for example dust radiative transfer (Steinacker et al. 2013). An additional process that needs to be considered when performing non-LTE line radiative transfer calculations is stimulated emission. By this process, a photon of a specific wavelength interacts with an atom or molecule in an excited state, causing it to drop to a lower energy level and liberating the internal energy in the form of a new photon. In many circumstances, stimulated emission can easily be incorporated in the MCRT framework ‘under the hood’, as absorption and stimulated emission are reverse processes proceeding at different rates, and one can combine them into a single process (e.g., Hogerheijde & van der Tak 2000; van der Tak et al. 2007; Brinch & Hogerheijde 2010). In most media, the absorption rate exceeds the stimulated emission rate because there are more electrons in the lower energy states than in the higher energy states. The combination of both processes then results in a net absorption cross section.

There are, however, situations where this assumption is not valid. When a population inversion is present, the stimulated emission rate exceeds the absorption rate, resulting in a net stimulated emission cross section, or equivalently, a formally negative absorption cross section. The result is that the intensity of the radiative field increases exponentially when it passes through the medium, a physical phenomenon known as optical amplification. In astrophysics, optical amplification is mainly important for masers. Traditional masers are produced by rotational or vibrational transitions in molecules such as OH, H₂0 or SiO, when the density is considerably higher than in the diffuse ISM and a luminous source acts as a pumping mechanism to generate a population inversion (Reid & Moran 1981; Elitzur 1992; Lo 2005). Another class consists of hydrogen recombination line masers, where the recombination lines are amplified by stimulated emission in the Rayleigh-Jeans limit, even at low optical depths (Krolik & McKee 1978; Martin-Pintado et al. 1989; Báez-Rubio et al. 2018).

The presence of net stimulated emission, or absorption with a negative cross section, does not fit the standard MCRT method. Indeed, one of the core ingredients in MCRT is the conversion of a randomly generated optical depth to a physical path length (Niccolini et al. 2003; Baes et al. 2003, 2019; Whitney 2011; Steinacker et al. 2013). If the cross section along a portion of the path is negative, the optical depth is no longer a monotonic function of the distance along the path, rendering the conversion impossible - or at least ambiguous - and this fundamental MCRT ingredient breaks down. One could attempt to incorporate net stimulated emission in MCRT using an iterative approach. First perform an MCRT phase that ignores stimulated emission. When this is finished, launch additional photon packets corresponding to stimulated emission in a second phase, and iterate this phase until convergence is reached. This approach has substantial disadvantages. The second iterative phase adds a significant layer of complexity and convergence might be slow because of the expected exponential growth of the radiation field intensity. Moreover, this approach requires storing the radiation field including directional information at every location, which is extremely expensive for 3D simulations.

In this paper, we present a technique, which we call explicit absorption, that allows us to efficiently handle net stimulated emission, or equivalently, absorption with negative cross sections, within a single MCRT phase. The text is organised as follows. We first describe the method in Sect. 2 and then verify its correctness using a straightforward 1D problem that can be solved analytically in Sect. 3. Subsequently, we discuss the implementation of the explicit absorption technique in our 3D MCRT code SKIRT (Camps & Baes 2020) in Sect. 4, verify its operation using an idealised 3D configuration and published benchmark problems in Sect. 5, and study its efficiency compared to the regular MCRT method in Sect. 6. Finally, we summarise and conclude in Sect. 7.

2 Monte carlo radiative transfer with stimulated emission

2.1 The radiative transfer equation

The standard time-independent monochromatic radiative transfer equation can be written as $$\frac{\text{d}I}{\text{d}s}\left(x,n\right)=j\left(x,n\right)-n\left(x\right)C_{\text{abs}}\left(x\right)I\left(x,n\right)-n\left(x\right)C_{\text{sca}}\left(x\right)I\left(x,n\right)+n\left(x\right)C_{\text{sca}}\left(x\right){\displaystyle \int I\left(x,n^{\prime }\right)\text{Φ}\left(n,n^{\prime }\right)\text{d}{\text{Ω}}^{\prime }}.$$(1)

In this expression, x and n represent the position and the propagation direction of the radiation, respectively. I represents the specific intensity of the radiation, j the emissivity, n the number density of the medium, C_abs and C_sca the absorption and scattering cross sections, and Φ the scattering phase function. In most cases, the second and third terms on the right-hand side of this equation are combined in a single extinction term, with the extinction cross section $$C_{\text{ext}}=C_{\text{abs}}+C_{\text{sca}}$$(2)

simply the sum of the absorption and scattering cross sections. We prefer to write them explicitly to stress the two processes at work. In the case of stimulated emission, the absorption cross section is the sum of the positive contribution of actual absorption and the negative contribution of stimulated emission. When stimulated emission dominates over absorption, we find that Cabs < 0.

For the sake of simplicity, we assume in this paper that the absorption and scattering cross sections are independent of the intensity of the radiation field. In real astrophysical situations, the level populations of the atoms and molecules in the medium, and thus the absorption and scattering coefficients, do depend on the radiation field strength. In particular, a population inversion will be reduced when the radiation field intensity is sufficiently high, ultimately leading to a saturation of the radiation field (e.g., Reid & Moran 1981; Elitzur 1990). The simplifying assumption we make comes down to assuming that we remain in the unsaturated regime.

2.2 Basic monte carlo radiative transfer

The simplest MCRT calculation considers radiation at a single wavelength. Throughout the calculation, the state of each individual photon packet is characterised by its weight W (equivalent to the number of photons within the packet), its position x, and its propagation direction n. Once the photon packet is launched into the medium, its life cycle is determined by a loop. This loop consists of three steps, during which the properties of the packet are updated.

The first step in the photon packet life cycle (PLC) consists of randomly determining whether the photon packet will interact with the medium, and if so, where this interaction will take place. The PDF that describes the free path length s before an interaction is most conveniently described in optical depth space, where it has an exponential distribution. In the light of the remainder of this section, it is useful to revisit the reasoning behind this. For a photon packet moving through an attenuating medium with density n(s) and extinction cross section C_ext(s), the radiative transfer equation reads $$\frac{\text{d}I}{\text{d}s}\left(s\right)=-n\left(s\right)C_{\text{ext}}\left(s\right)I\left(s\right),$$(3)

where we note that no additional emission processes are relevant for this particular photon packet. The solution of this simple ordinary differential equation is readily found, $$I\left(s\right)=I\left(0\right){\text{e}}^{-{\tau }_{\text{ext}}\left(s\right)},$$(4)

with the monotonically increasing optical depth scale τ_ext(s) defined as $${\tau }_{\text{ext}}\left(s\right)={\displaystyle {\int }_0^{s}n\left(s^{\prime }\right)C_{\text{ext}}\left(s^{\prime }\right)\text{d}{s}^{\prime }}.$$(5)

The exponential decline of the intensity in reality means that individual photons are continuously disappearing from the photon packet due to interactions with the medium as the photon packet propagates. Ideally, we should model all of these interactions individually. The essence of the Monte Carlo method consists of replacing this continuous behaviour by stochastically sampling a single interaction location for the entire photon packet. The challenge is to find the appropriate PDF p(s) from which we should sample a physical path length before the next interaction. This PDF can be found by noting that the probability that a photon does not interact along the path between 0 and s is equal to I(s)/I(0). The cumulative distribution P(s) corresponding to the desired PDF p(s) is thus $$P\left(s\right)=1-\frac{I\left(s\right)}{I\left(0\right)}=1-{\text{e}}^{-{\tau }_{\text{ext}}\left(s\right)},$$(6)

or in optical depth space, $$P\left({\tau }_{\text{ext}}\right)=1-{\text{e}}^{-{\tau }_{\text{ext}}}.$$(7)

The PDF is found by differentiation, and turns out to be a simple exponential distribution, $$p\left({\tau }_{\text{ext}}\right)={\text{e}}^{-{\tau }_{\text{ext}}}.$$(8)

Generating a random τ_ext from an exponential distribution is straightforward. If this random optical depth exceeds the total optical depth along the path, $${\tau }_{\text{ext,\hspace{0.17em}path}}={\displaystyle {\int }_0^{\infty }n\left(s\right)C_{\text{ext}}\left(s\right)\text{d}s},$$(9)

the photon packet survives the journey along the path and exits the medium. If, on the other hand, the random optical depth is smaller than the total optical depth, we can calculate the physical path length s to the next interaction location by inverting the relation in Eq. (5). For complex 3D models, this inversion is often the most time-consuming part of the entire PLC. It is usually accomplished through the use of discretisation on a spatial grid, combined with efficient grid traversal algorithms (e.g., Baes et al. 2003; Camps et al. 2013; Saftly et al. 2013, 2014; Hubber et al. 2016). With the physical path length s determined, the photon packet can be propagated to the next interaction location.

The second step in the PLC, at least if the photon packet has not left the system before the interaction, is the determination of the nature of the interaction at the interaction site. The photon packet can either be absorbed or scattered; the appropriate PDF is hence not a continuous but a discrete one with only two possible values. If the absorption and scattering cross sections at the location x are C_abs(x) and C_sca(x), respectively, the discrete PDF that describes the probabilities for both types is $$p\left(\text{abs}\right)=\frac{C_{\text{abs}}\left(x\right)}{C_{\text{ext}}\left(x\right)}=1-\omega \left(x\right),$$(10a) $$p\left(\text{sca}\right)=\frac{C_{\text{sca}}\left(x\right)}{C_{\text{ext}}\left(x\right)}=\omega \left(x\right),$$(10b)

where ω represents the scattering albedo. Drawing a random interaction type from this discrete PDF is trivial. If the event is an absorption event, the photon packet’s life cycle terminates. In the case of a scattering event, its journey continues, and we move to the third step in the cycle.

The third and final step in the PLC is the actual simulation of a scattering event at the location of the interaction. This comes down to generating a new propagation direction for the photon packet, based on the scattering phase function. The details of this step are not very relevant for the discussion in the remainder of this paper.

Each of these three steps are repeated until the photon packet either escapes from the system or is absorbed by the medium. We can then continue with the next photon packet and repeat the entire procedure. This entire process is repeated for all the photon packets until the last one has left the medium. The ultimate goal of a MCRT simulation is usually the determination of the radiation field, either at locations within the model space or corresponding to an external viewpoint. This can be achieved by statistically binning the photon packets according to location, propagation direction and weight. The details of these steps are also less relevant for the remainder of the discussion.

2.3 Challenges in the case of net stimulated emission

The basic MCRT algorithm discussed in the previous section assumed that both the absorption and the scattering cross sections are positive. As discussed in the Introduction, this is not necessarily the case if we consider stimulated emission as a negative contribution to the absorption. In the case that the stimulated emission rate dominates over the absorption rate, we have C_abs < 0, and we encounter two problems for the MCRT method, depending on the ratio between C_abs and C_sca.

If stimulated emission dominates over absorption, but the scattering cross section is sufficiently large such that the total extinction cross section is still positive, we are in what we call the weak net stimulated emission regime. In this regime, the first step in the PLC can proceed as described above. We do have an issue with the second step, the determination of the nature of the interaction, however. Indeed, since C_abs < 0, we formally have from Eqs. (10) that p(abs) < 0 and p(sca) > 1. This is obviously not a valid discrete PDF from which we can stochastically sample.

If the stimulated emission rate is so large or the scattering cross section so low that C_abs is larger in magnitude than C_sca, we enter the strong net stimulated emission regime. In this case the total extinction cross section becomes negative. Concerning the second step in the PLC we have from Eqs. (10) that p(abs) > 1 and p(sca) < 0, which leads to a similar problem as before. Even more troublesome is that we now encounter a fundamental problem in the first step of the cycle. If the extinction cross section is negative in some regions of the computational domain, we find that the optical depth along each path will no longer be a monotonically increasing function of physical path length, which means that we can no longer interpret Eq. (6) as a valid cumulative distribution function. It can even occur that the optical depth is completely negative, which completely breaks the procedure.

2.4 The absorption-scattering split

The way to deal with the challenges posed by net stimulated emission, or negative absorption cross sections, is to introduce deterministic elements in the Monte Carlo PLC. The introduction of deterministic elements is, next to biasing or variance reduction, usually done to accelerate the algorithm (for an overview, see Steinacker et al. 2013, Sect. 5.2). As we will show here, they can also be used to extend the application range of MCRT. A simple example of the introduction of a deterministic element in the PLC is the so-called absorption-scattering split (e.g., Mattila 1970; Witt 1977). The second step in the standard PLC consists of randomly deciding whether the photon packet will absorb or scatter at the interaction site. In reality, there will be both absorption and scattering: a fraction of the photons in the packet will be absorbed, and the remaining fraction will scatter. Rather than stochastically determining the nature of the interaction, we can split the photon packet into two parts: one that is absorbed and one that scatters. We know exactly which fraction should be absorbed and scattered: this is given by 1 - ω and ω, respectively. So at the location of the interaction, we replace the original photon packet with weight W that is about to interact by two new photon packets: one with weight (1 - ω) W that will be absorbed, and one with weight ω W that will scatter. In practice, since the absorbed energy disappears, we are only concerned with the photon packet that scatters and continues its journey.

This very simple change to the PLC replaces a stochastic choice of the interaction nature by a deterministic procedure. Since stochasticity always introduces Poisson noise, this straightforward adjustment helps to reduce noise. In other words, by introducing the absorption-scattering split, we need fewer photon packets in our simulation to reach a given level of signal-to-noise compared to a simulation that does not incorporate it.

The same method can, however, also be used to allow for stimulated emission in the weak net stimulated emission regime. In this regime we have 0 < −C_abs < C_sca, so formally ω > 1. If a photon packet with weight W is to interact, we create two new photon packets with weights (1 - ω) W and ω W, respectively. The former packet has a negative weight, but this does not jeopardise the PLC. As it will be absorbed, this photon package does not travel through the medium and and does not contribute to the radiation field; the negative weight is thus not an issue. In case the absorption rate of a cell in the medium, and subsequently the mean intensity of the radiation field, the level populations, and the emission rate, are calculated by summing the weights of all absorbed photon packets in that cell (e.g., Gordon et al. 2001; Bjorkman & Wood 2001; Baes et al. 2005), the resulting absorption rate will be negative. A negative absorption rate simply reflects the fact that we have net stimulated emission: since the absorption cross section is negative as well, the standard formulae to determine the mean radiation field intensity can still be used¹. The result is that we can continue the PLC with the second photon packet, which has a positive weight that is larger than the original one. In this way we can simulate net stimulated emission, at least in the weak net stimulated emission regime.

2.5 Explicit absorption

The absorption-scattering split can be used to deal with stimulated emission in the weak net stimulated emission regime, but not in the strong net stimulated emission regime. In this latter regime we have the more fundamental problem of negative extinction cross sections. We can, however, introduce another deterministic element in the Monte Carlo PLC that can overcome this problem. We call this approach explicit absorption.

The first step in the standard PLC consists of picking a random path length for the next interaction location based on the interpretation of the intensity in Eq. (4) as a cumulative distribution in Eq. (6). Since the extinction cross section is just the sum of the absorption and scattering cross sections, we can also rewrite Eq. (4) as $$I\left(s\right)=\left[I\left(0\right){\text{e}}^{-{\tau }_{\text{abs}}\left(s\right)}\right]{\text{e}}^{-{\tau }_{\text{sca}}\left(s\right)},$$(11)

with the absorption and scattering optical depth scales defined as: $${\tau }_{\text{abs}}\left(s\right)={\displaystyle {\int }_0^{s}n\left(s^{\prime }\right)C_{\text{abs}}\left(s^{\prime }\right)\text{d}{s}^{\prime }},$$(12) $${\tau }_{\text{sca}}\left(s\right)={\displaystyle {\int }_0^{s}n\left(s^{\prime }\right)C_{\text{sca}}\left(s^{\prime }\right)\text{d}{s}^{\prime }}.$$(13)

Based on this equation, we can change the PLC. Rather than simulating the entire extinction stochastically, we split the absorption and scattering contributions. We treat the absorption contribution explicitly and only the scattering contribution stochastically. We implement the explicit absorption by continuously decreasing the intensity, or the weight, of the photon packet as it travels along the path according to: $$W\left(s\right)=W\left(0\right){\text{e}}^{-{\tau }_{\text{abs}}\left(s\right)}.$$(14)

On the other hand, we randomly determine the next scattering location by picking a random scattering optical depth τ_sca from an exponential PDF, $$p\left({\tau }_{\text{sca}}\right)={\text{e}}^{-{\tau }_{\text{sca}}},$$(15)

and we compare it to the total scattering optical depth along the path, $${\tau }_{\text{sca,\hspace{0.17em}path}}={\displaystyle {\int }_0^{\infty }n\left(s\right)C_{\text{sca}}\left(s\right)\text{d}s}.$$(16)

If τsca > τ_{sca, path}, the photon packet leaves the system. If, on the other hand, τ_sca < τ_sca,path, the photon packet will scatter along this path. The scattering location is determined by inverting the relation in Eq. (13) for s. In both cases, we should not forget to alter the photon packet’s weight according to expression in Eq. (14).

We note that this alteration of the Monte Carlo method reduces the PLC to two steps. Indeed, the second step in the original method, the random determination of the nature of the interaction, disappears. The remaining steps are the determination of and the propagation to the next scattering location, and the simulation of the new direction after the scattering. Since we eliminate a stochastic element from the PLC, this alteration should in principle reduce the noise. On the other hand, it also makes the first step more computationally demanding: we need to determine both the absorption and the scattering optical depth along the path, and take into account the continuous change of the weight of the photon packet.

The most interesting aspect of this explicit treatment of absorption is that it can also be applied in the case of negative absorption cross sections. The only effect of absorption is the continuous change of the weight of a photon packet as it moves through the medium according to Eq. (14). Whether τ_abs(s) is positive or negative, and hence whether the weight decreases or increases as a function of distance along the path, does not really matter. For the determination of the scattering location, on the other hand, it is crucial that the optical depth is a non-decreasing function of distance along the path. Since C_sca is always non-negative, this condition is satisfied.

3 A two-stream radiative transfer problem

3.1 Problem Definition

We test our new approach using a simple 1D radiative transfer problem for which the solution can be calculated analytically. We consider radiation that propagates through a finite column of matter of length x_max = 1. We indicate the distance within the column as x, and we denote the intensity of the radiation field moving in positive and negative directions as I⁺(x) and I⁻(x), respectively. A single source of radiation with unit intensity that emits in the positive direction is placed on the left side of the column. The density of matter and the optical properties within the column are uniform, and we arbitrarily set n = 1. The scattering phase function is such that half of the radiation scatters in the forward direction and the remaining half scatters in the backward direction. The problem is depicted graphically in Fig. 1.

With these assumptions, the radiative transfer Eq. (1) can be transformed into a set of two coupled differential equations for I⁺(x)aηd I⁻(x), $$\frac{\text{d}{I}^{\pm }}{\text{d}x}\left(x\right)=\mp \left(C_{\text{abs}}+\frac{1}{2}{C}_{\text{sca}}\right)I^{\pm }\left(x\right)\pm \frac{1}{2}{C}_{\text{sca}}{I}^{\mp }\left(x\right).$$(17)

The boundary conditions are: $$I^{+}\left(0\right)=1,$$(18) $$I^{-}\left(x_{\max }\right)=0.$$(19)

Fig. 1 Graphical illustration of the two-stream radiative transfer problem considered in Sect. 3.

3.2 Analytical solution

For any value of the cross sections C_abs and C_sca, this set of coupled equations can be solved analytically. We note that, as discussed in the previous sections, the scattering cross section is always positive, whereas the absorption cross section can also be negative, as it combines the contribution of actual absorption and stimulated emission. We assume that C_sca > 0 and we consider five different cases, depending on the value of C_absnnet absorption (C_abs > 0), pure scattering (C_abs = 0), weak net stimulated emission (−C_sca < C_abs < 0), critical net stimulated emission (−C_sca = C_abs < 0), and strong net stimulated emission (C_abs < −C_sca < 0). These five different regimes are indicated on Fig. 2, which shows the (C_abs, C_sca) parameter space for this simple radiative transfer problem.

3.2.1 Net absorption

In the regime of net absorption, thus with C_abs > 0, the set of coupled differential Eq. (17) has two real eigenvalues. Given the boundary conditions in Eqs. (18)–(19), we can write the full solution as: $$I^{+}\left(x\right)=\frac{\xi \cosh \xi \left(1-x\right)+\zeta \sinh \xi \left(1-x\right)}{\xi \cosh \xi +\zeta \sinh \xi },$$(20) $$I^{-}\left(x\right)=\frac{\frac{1}{2}{C}_{\text{sca}}\sinh \xi \left(1-x\right)}{\xi \cosh \xi +\zeta \sinh \xi },$$(21)

where we have defined the constants $$\xi =\sqrt{C_{\text{abs}}\left(C_{\text{abs}}+C_{\text{sca}}\right)},$$(22) $$\zeta =C_{\text{abs}}+\frac{1}{2}{C}_{\text{sca}}.$$(23)

Both these constants are strictly positive, and so are the solutions in Eqs. (20)–(21) for all values of C_abs and C_sca in this regime.

Fig. 2 (Cabs, Csca) parameter space for the 1D radiative transfer model discussed in Sect. 3.3. The green shaded area in this plot corresponds to the region with physically viable solutions. The three major regimes (net absorption, weak net stimulated emission, and strong net stimulated emission) are indicated. The orange and red lines indicate the boundaries between these major regimes, and correspond to pure scattering and critical net stimulated emission, respectively. The coloured dots indicate the models shown explicitly in Fig. 3. The white area corresponds to non-physical solutions. The solid green line that indicates the boundary between physical and non-physical solutions is given by the Eqs. (30) and (34).

3.2.2 Pure scattering

When C_abs = 0, the differential Eq. (17) simplify to: $$\frac{\text{d}{I}^{\pm }}{\text{d}x}\left(x\right)=-\frac{1}{2}{C}_{\text{sca}}\left[I^{+}\left(x\right)-I^{-}\left(x\right)\right].$$(24)

Subject to the boundary conditions in Eqs. (18)–(19), the solution becomes: $$I^{+}\left(x\right)=\frac{2+C_{\text{sca}}\left(1-x\right)}{2+C_{\text{sca}}},$$(25) $$I^{-}\left(x\right)=\frac{C_{\text{sca}}\left(1-x\right)}{2+C_{\text{sca}}}.$$(26)

The same solution can also be found by taking the limit C_abs → 0 in the expressions in Eqs. (20)–(23). It is immediately clear that this solution is positive, and thus physically valid, for all values Of C_sca.

3.2.3 Weak net stimulated emission

In the weak net stimulated emission regime, we can still use the expressions in Eqs. (20)–(23) as the formal solution of the coupled set of radiative transfer Eq. (17). However, since C_abs < 0 and C_abs + C_sca > 0, ξ is an imaginary number. A more elegant way to write the solution is: $$I^{+}\left(x\right)=\frac{{\xi }^{\prime }\cos {\xi }^{\prime }\left(1-x\right)+\zeta \sin {\xi }^{\prime }\left(1-x\right)}{{\xi }^{\prime }\cos {\xi }^{\prime }+\zeta \sin {\xi }^{\prime }},$$(27) $$I^{-}\left(x\right)=\frac{\frac{1}{2}{C}_{\text{sca}}\sin {\xi }^{\prime }\left(1-x\right)}{{\xi }^{\prime }\cos {\xi }^{\prime }+\zeta \sin {\xi }^{\prime }},$$(28)

where we have set: $${\xi }^{\prime }=\sqrt{-C_{\text{abs}}\left(C_{\text{abs}}+C_{\text{sca}}\right)}>0.$$(29)

Interestingly, instead of the standard exponential functions usually encountered in radiative transfer problems, we now have trigonometric functions, as typically encountered in solutions of coupled differential equations with imaginary eigenvalues.

Furthermore, it turns out that this formal solution does not always correspond to a physically sound solution: only when both I⁺(x) and I−(x) are positive at all positions x along the column, the formal solution can be physical. This limits the region in the (C_abs, C_sca) parameter space that corresponds to physical solutions. For C_sca ⩽ 2, the solutions are always physical in the weak net stimulated emission regime. For every value of C_sca > 2, however, there is a limit value of C_abs below which the solutions in Eqs. (27)–(28) become negative. This limit value is reached when the denominator of these expressions becomes zero, that is when: $$\frac{\tan {\xi }^{\prime }}{{\xi }^{\prime }}=\frac{1}{\zeta }.$$(30)

For each C_sca > 2, this equation can be solved numerically to determine the limit value of C_abs. The boundary between the physical and non-physical models in the weak net stimulated emission regime is indicated as the upper part of the thick green line in Fig. 2.

3.2.4 Critical net stimulated emission

The critical net stimulated emission case corresponds to the situation where absorption, scattering and stimulated emission are perfectly in balance such that their combined cross section is zero, that is, C_abs + C_sca = 0. For this special case, we can rewrite the differential Eq. (17) as: $$\frac{\text{d}{I}^{\pm }}{\text{d}x}\left(x\right)=\pm \frac{1}{2}{C}_{\text{sca}}\left[I^{+}\left(x\right)+I^{-}\left(x\right)\right].$$(31)

Taking into account the boundary conditions in Eqs. (18)–(19), the solution is: $$I^{+}\left(x\right)=\frac{2-C_{\text{sca}}\left(1-x\right)}{2-C_{\text{sca}}},$$(32) $$I^{-}\left(x\right)=\frac{C_{\text{sca}}\left(1-x\right)}{2-C_{\text{sca}}}.$$(33)

Also in this case, the solutions are not physically valid for all values of C_sca. In fact, it is easy to see that they are only physical when C_sca < 2, or equivalently, when C_abs > −2.

3.2.5 Strong net stimulated emission

Finally, we consider the strong net stimulated emission regime where the net stimulated emission cross section is so large that it dominates over the combined absorption and scattering cross sections. In this case, the solution of the radiative transfer problem is formally exactly the same as in the net absorption regime, that is, the expressions in Eqs. (20)–(21).

Since both C_abs and C_abs + C_sca are negative in this regime, we have that ξ > 0. The main difference with the case of net absorption is that ζ is negative in this case, which implies that the denominator will become negative for certain combinations of C_abs and C_sca. For every value of C_sca < 2, there is a limit value of C_abs below which the solutions in Eqs. (20)–(21) become negative. This limit value is reached when: $$\frac{\tanh \xi }{\xi }=-\frac{1}{\zeta }.$$(34)

For C_sca > 2, the radiative transfer problem has no physical solutions in the strong net stimulated emission regime. The boundary between physical and non-physical solutions in the strong net stimulated emission regime smoothly connects to the similar boundary in the weak net stimulated emission regime, and is indicated as the thick green line in Fig. 2.

3.2.6 No scattering

To streamline the discussion of the various regimes in the previous subsections, we assumed that C_sca > 0. We now consider the case where there is no scattering, i.e., C_sca = 0. The radiative transfer Eq. (17) simplify drastically and are no longer coupled. Taking into account the boundary conditions in Eqs. (18)–(19), the trivial solution is: $$I^{+}\left(x\right)=\exp \left(-C_{\text{abs}}x\right),$$(35) $$I^{-}\left(x\right)=0.$$(36)

This solution is physically valid for all finite values of C_abs. The pure scattering, weak net stimulated emission, and critical net stimulated emission regimes collapse into the single point where C_sca = C_abs = 0. This no extinction regime is represented in Fig. 2 by the intersection of the orange and red lines.

3.3 Qualitative discussion of the solution

In Fig. 3, we illustrate the systematic change of the solution of the two-stream radiative transfer problem as we vary the model parameters. We have fixed C_sca = 1, and the different lines in this figure represent solutions of the RTE when we gradually decrease the absorption cross section C_abs from 2 to −2.5 in steps of ∆C_abs = −0.25. In Fig. 2, we show the position of each of these models in the (C_abs, C_sca) parameter space. This set of parameters covers each of the five regimes discussed in the previous subsection.

For positive values of C_abs, we are obviously in the net absorption regime. The largest value of C_abs corresponds to the strongest attenuation of the intensity in the positive direction. The intensity profile I⁺(x) is almost a perfect exponentially decaying function, as expected when absorption strongly dominates over scattering. The intensity in the negative direction, which is entirely due to scattered radiation, is also a monoton-ically decreasing function of x. When C_abs decreases, the total extinction cross section decreases as well. The result is that I⁺(x) gradually becomes flatter and less exponential-like, and that I⁺(x) gradually increases with decreasing C_abs.

When C_abs becomes zero, the model contains only scattering. Compared to the models with net absorption, the intensity profile in the positive direction is even flatter. As shown by expression in Eq. (25), the I⁺(x) profile is now a monotonically decreasing linear function of x rather than an exponential function. The I⁻(x) profile is also a linearly decreasing function of x that is globally increased compared to the solutions of the net absorption regime.

As soon as C_abs decreases even more and hence becomes negative, we enter the net stimulated emission regime. We first enter the weak net stimulated emission regime, in which the scattering cross section dominates over the negative absorption cross section. As C_abs gradually becomes more negative, the intensity profile in the forward direction becomes flatter and flatter, as a result of the competition between stimulated emission and scattering. At a certain point, I⁺(x) becomes a non-monotonic function of x: it first slightly increases as a function of x and subsequently slightly decreases. More specifically, this starts at C_abs = −0.2824. A nice illustration of such a model is the case C_abs = −0.3799: for this specific value, we find that I⁺(1) = 1, that is, the forward intensity at the start and the end point of the column has exactly the same value. It does reach a maximum value of 1.0302 at x = 0.5. If C_abs becomes even more negative, there is no longer attenuation and the intensity at the end point of the column exceeds the input intensity. As soon as C_abs < −0.5,I⁺(x) becomes a monotonically increasing function of x, with a strong initial increase and a gradual flattening at large values of x. Over the entire range of C_abs values, the negative intensity profile remains a decreasing function of x, and the magnitude steadily increases as C_abs becomes more negative.

When C_abs = −1, we end up in the situation where the absorption, scattering and stimulated emission cross sections are perfectly in equilibrium, that is, we are in the critical net stimulated emission regime. As can be seen from Eq. (35), I⁺(x) is now a monotonically increasing linear function of x. The intensity in the negative direction is a monotonically decreasing linear function of x.

Finally, when C_abs gets even more negative, we enter the strong net stimulated emission regime. The smaller C_abs, the stronger the forward intensity profile I⁺(x) will grow as a function of x. It is interesting to see that the slope of the I⁺(x) curve is strongest at small x and gradually flattens at larger x. The intensity profile in the negative direction remains a decreasing function of x, but its magnitude increases steadily as C_abs grows more negative. At some point, the exponential growth rate of I⁺(x) at small x is so strong that the intensity diverges before it reaches the other side of the column. This is the point where we reach the green curve in Fig. 2, corresponding to the relation in Eq. (34). This is where the solution in Eqs. (20)–(21) becomes unphysical. For the current choice of parameters, the critical boundary between physical and unphysical intensity profiles is C_abs = −2.7340.

Fig. 3 Intensity of the radiation field in the positive (leftpanel) and negative (rightpanel) direction for the 1D radiative transfer model discussed in Sect. 3.3. All models have the same scattering cross section Csca = 1 and the different lines correspond to models with different absorption cross sections Cabs, ranging from 2 to −2.5 in steps of −0.25. The thick orange and red curves correspond to pure scattering and critical net stimulated emission, respectively.

3.4 Monte Carlo simulations

For the simple ID radiative transfer problem defined in the previous section, we constructed a simple Monte Carlo code. We implemented three different variants: one version with the standard implementation as described in Sect. 2.2, one version with the absorption-scattering split (Sect. 2.4), and one version with explicit absorption (Sect. 2.5). In each case, the radiation field is recorded at 50 locations spread uniformly along the column.

In Fig. 4, we show the results of our simulations. In each panel, the solid lines are the analytical results also shown in Fig. 3 and the dots are the results of the MCRT calculations. The top panels show the results for the standard implementation. This procedure only works for C_abs ⩾ 0 and it clearly reproduces the analytical results perfectly.

The panels on the middle row show the results of MCRT simulations including the absorption-scattering split. For the models with Cabs ⩾ 0, the analytical results are still reproduced accurately. As argued in Sect. 2.4, the introduction of the absorption-scattering split also allows us to extend the applicability of the MCRT method to the regime of weak net stimulated emission. We see that, indeed, the MCRT results also reproduce the analytical results for all values C_abs > −C_sca = −1.

Finally, the bottom panels show the results for the MCRT simulations with explicit absorption. In this case, the MCRT simulations can be applied for all parameter settings, including the regimes of critical and strong net stimulated emission. The MCRT results reproduce the analytical results accurately over the entire range of C_abs values. Only for the most negative values of C_abs, we see deviations from the analytical result. The deviation is particularly notable for C_abs = −2.5, the most negative value considered. This value is close to the limit value −2.7340 that C_abs can have before the problem diverges and becomes unphysical.

4 Implementation in SKIRT

Having demonstrated the power of the explicit absorption technique using a simple 1D problem, we now describe its implementation in a state-of-the-art 3D code.

Fig. 4 Monte Carlo simulation results of the 1D radiative transfer problem discussed in Sect. 3.3. The three rows correspond to three different variants of the MCRT method: the standard implementation (top row), the version with the absorption-scattering split (middle row), and the version with explicit absorption (bottom row). Solid lines are the analytical results, the dots are the results of the MCRT calculations.

4.1 The SKIRT code

SKIRT² (Baes et al. 2011; Camps & Baes 2015, 2020) is a feature-rich 3D MCRT code that is commonly being applied to construct synthetic observations of astrophysical models. Recent examples include Di Mascia et al. (2021a,b); Granato et al. (2021); Vandenbroucke et al. (2021); Kapoor et al. (2021); Lin et al. (2021); Whitney et al. (2021); Mosenkov et al. (2021); Olsen et al. (2021); Victoria-Ceballos et al. (2022); Popping et al. (2022); Shen et al. (2022); Vijayan et al. (2022); Flores-Freitas et al. (2022). SKIRT offers an extensive library of built-in geometries (Baes & Camps 2015) as well as the capability of importing models from hydro-dynamical simulation snapshots (Saftly et al. 2015; Camps et al. 2016). Sources can be assigned spectra from built-in or user-provided stellar or stellar population template libraries (Camps & Baes 2015). Media types and physical processes include dust extinction and emission using turn-key or highly configurable dust models (Camps et al. 2015), polarisation by dust grains (Peest et al. 2017; Vandenbroucke et al. 2020), Lyman-α resonant line transfer (Camps et al. 2021), electron Thomson and Compton scattering, photo-ionisation and fluorescence by neutral atoms (Vander Meulen et al., in prep.), and kinematics (Baes et al. 2003; Camps & Baes 2020). Supported wavelengths span the X-ray, UV to sub-millimetre, and radio regimes.

On the numerical side, SKIRT offers several spatial discretisation options, including 1D spherical and 2D cylindrical grids for models with the corresponding symmetries, and advanced hierarchical and unstructured grids for high-dynamic-range 3D models (Saftly et al. 2013, 2014; Camps et al. 2013). The configuration for a given model can include any mixture of source and media types. The Monte Carlo PLC supports primary emission by the sources, absorption and scattering by the media, and secondary emission by eligible media, while tracking the radiation field and recording any requested observables. Self-consistent iteration can be enabled to properly handle processes such as dust destruction near high-intensity sources (during primary emission) or dust self-absorption in regions with high optical depth (during secondary emission). The PLC runs in parallel on multiple execution threads and/or in multiple processes, allowing the code to efficiently scale from laptops to multi-node high-performance clusters (Verstocken et al. 2017).

SKIRT uses several acceleration and variance-reduction techniques. These include peel-off towards the instruments at emission and scattering events for recording observables (Yusef-Zadeh et al. 1984), the absorption-scattering split described in Sect. 2.4 (Mattila 1970; Witt 1977), continuous absorption along the photon packet path to improve sampling of the radiation field (Lucy 1999; Niccolini et al. 2003), forced scattering to improve sampling of optically thin regions (Cashwell & Everett 1959), path length biasing to more easily penetrate optically thick regions (Baes et al. 2016), and wavelength biasing to improve signal-to-noise in wavelength ranges of interest. Many of these options can be configured by the user. For example, forced scattering can be turned off to speed up resonant line transfer simulations, where each photon packet undergoes a potentially large number of scattering interactions with very short free path lengths (often within the same grid cell). In addition to all this, the code switches to optimised versions of some procedures depending on the configuration. For example, if a photon packet’s wavelength cannot change along its path because the model has no kinematics and optical media properties are constant across space, the absorption and scattering cross sections can be calculated once at the start of the path instead of repeating the calculation for each crossed cell.

4.2 Adding explicit absorption

In view of the complexity described in the previous subsection, adding the explicit absorption technique to the SKIRT PLC seems daunting. The most important question is whether explicit absorption can co-exist with the other acceleration and variance reduction techniques, most notably peel-off, continuous absorption, forced scattering, and path length biasing. A secondary concern is whether explicit absorption can be added as a user-configurable run-time option without impacting performance of the existing PLCs, while limiting the amount of code duplication. It turns out that, with the proper implementation choices, co-existence with the other techniques is not a problem and supporting the various types of PLCs in a single code is quite feasible.

In the updated SKIRT version³, we allow enabling or disabling forced scattering and/or explicit absorption independently of each other. This yields four PLC variants, each of which has specific advantages or drawbacks.

Forced scattering tends to reduce noise for simulations with low optical depths. However, for each scattering event, it requires the calculation of the geometry and optical depth of the full path up to the model boundary, regardless of the location of the scattering event along the path. In models with very intensive scattering, such as for Lyman-α line transfer, the PLC without forced scattering is often the better choice, because it avoids calculating the path geometry and optical depth beyond the scattering location. However, our implementation does not support storing the radiation field. The main reason for this choice is that some additional optimisations can be applied. For example, there is no need to store the calculated path segment information. As a consequence, forced scattering cannot be turned off for simulations that require the radiation field to calculate secondary emission by media.

The PLC without explicit absorption uses the extinction (sum of scattering and absorption) along a photon packet’s path to locate the next interaction point. This requires the cumulative extinction optical depth to be a non-decreasing function of path length. It is thus not possible to handle negative extinction cross sections. The PLC with explicit absorption instead uses the scattering optical depth to locate the next interaction point. While the scattering cross section still must be non-negative, this allows the extinction cross section to be negative. As a drawback, this technique requires calculating both the scattering and absorption optical depths for the photon packet path.

4.3 The photon packet

At the heart of the Monte Carlo PLC in SKIRT is the C++ class used for representing a photon packet. A photon packet tracks all relevant properties during its life cycle, including wavelength, weight (i.e., number of photons), polarisation state, position, and propagation direction. It can also store details on the path that traverses the simulation’s spatial grid in the packet’s current propagation direction and starting at its most recent interaction or emission location. Specifically, for each path segment n crossing a given spatial grid cell, this information includes the cumulative distance s_n and the cumulative optical depth τ_n up to the point where the path exits the cell. Once calculated and stored, these path segment properties can be used to determine the next interaction location and to add the packet’s contribution to the radiation field in each crossed cell.

To support explicit absorption, we adjusted the path segment data structure so that it can now hold two distinct cumulative optical depth values, ${\tau }_n^{\text{ext}|\text{sca}}$ and ${\tau }_n^{\text{abs}}.$ When explicit absorption is enabled, the PLC code calculates and stores the scattering and absorption optical depths ${\tau }_n^{\text{sca}}$ and ${\tau }_n^{\text{abs}}$ in these fields. When explicit absorption is not used, the code instead calculates and stores the extinction optical depth ${\tau }_n^{\text{ext}}$ in the first field, and the second field is set to zero. As a result of this scheme, the procedure for locating the next interaction point can always use the ${\tau }_n^{\text{ext}|\text{sca}}$ field, which will store the proper value depending on whether explicit absorption is enabled or not, and the procedure for adjusting the radiation field can always obtain the extinction optical depth through a simple addition, i.e. ${\tau }_n^{\text{ext}|\text{sca}}+{\tau }_n^{\text{abs}}.$ With this ‘trick’, we were able to substantially limit the changes to the code in areas that are not specifically concerned with the explicit absorption technique.

4.4 The photon packet life cycle

The top-level function driving the SKIRT PLC already had two separate code branches, one with and one without forced scattering. To allow selecting explicit absorption at run-time as well, we added an extra code path at just two places in each branch. When calculating the optical depths for the segments along the current path, we invoke a different function to ensure that the appropriate type of optical depth values are stored in the photon packet (see Sect. 4.3). And when applying the bias factor to the packet’s weight for a scattering interaction, we either use the scattering albedo (without explicit absorption) or the absorption extinction up to the interaction location according to Eq. (14) in Sect. 2.5 (with explicit absorption). The required if-then-else statements are performed only once per photon packet path segment and thus do not affect performance.

For most simulation configurations, a substantial fraction of the execution time is spent in the function calculating the path segments and the corresponding optical depths for the photon packet path. It is therefore justified to optimise this portion of the code for speed. There now is a separate version of this function for each of the four variants of the PLC. Furthermore, each of these functions have separate code paths depending on the overall configuration of the simulation. For example, there is a version optimised for configurations that allow calculating the absorption and scattering cross sections just once at the start of each path. In total, these various implementations occupy under 600 lines of code in a single class, including comments and documentation, resulting in an acceptable trade-off between complexity and performance.

No other changes to the code were needed, as we will verify in the next section. When explicit absorption is enabled, the photon packets follow different random paths because the interaction points are determined based on the scattering optical depth rather than the extinction optical depth. However, the photon weight is adjusted differently to compensate, so that the other aspects of the simulation continue to work properly. Specifically, the radiation field is properly tracked and the peel-off photon packets properly sample the fluxes observed by the configured instruments. Path length stretching, when enabled, also continues working as expected in the adjusted context.

5 Verification of the SKIRT implementation

5.1 Regression tests

Over the last decade of SKIRT development, we continuously maintained and expanded our set of functional test cases used for regression testing. We now have about 700 test cases that cover most features and code paths and can be run in 10 min or so on a present-day desktop computer. Each test simulation runs on a single execution thread and uses a fixed pseudo-random number sequence. This ensures that output files are binary identical between runs, enabling automated verification. Using this test suite, we easily confirmed that the updated implementation described in Sect. 4 does not break the already present PLCs.

5.2 The 1D two-stream problem

As an initial test for the explicit absorption technique, we set up a 3D SKIRT model that is equivalent to the 1D two-stream problem discussed in Sect. 3. The model includes a narrow tube of material along the x-axis and a point source injecting a laser beam into the tube at one end. The material has configurable absorption and scattering cross sections at the laser’s wavelength. When a photon packet scatters, the interaction reverses the packet’s direction or leaves it unchanged with fifty per cent probability for each action. The simulation is further configured to record the radiation field along the tube in regularly spaced grid cells. This automatically enables forced scattering (see Sect. 4.2). Because SKIRT does not record directional information for the radiation field, the resulting quantity is equivalent to the sum of the intensities, I⁺(x) + I⁻(x), in the 1D problem.

Running this model for a set of relevant cross section values sampling the shaded area in Fig. 2, we confirmed that SKIRT properly reproduces the analytic solutions listed in Sect. 3. We note that the model used here can be constructed with standard SKIRT components except for the special material properties, which we defined through a trivial custom class. Because of its 1D scattering behaviour, this material definition is not very useful for general 3D modelling. A more useful variation described in the next section has been included as a standard SKIRT component.

5.3 A spiral galaxy model

To test our implementation with an actual 3D geometry for which the results are still easily interpreted, we use the smooth spiral galaxy model presented by Peest et al. (2017) in their Sect. 6 with properties listed in their Table 2. In short, the model includes a stellar bulge and disk with an older star population, a second stellar disk with a younger star population, and a dust disk. The density distributions of the stellar and dust disks are perturbed by a smooth spiral arm structure. The left panels of Fig. 5 illustrate the observed V-band surface brightness as calculated by SKIRT for near face-on and near edge-on viewing angles.

We replace the dust in the model by an idealised medium for which we can directly configure the optical properties at the single, arbitrary wavelength used in our tests. We configure a slight preference for forward scattering (<cos 9> = 0.25), set the scattering cross section C_sca to a fixed value, and let the absorption cross section C_abs vary as a multiple of C_sca, i.e. C_abs = kC_sca. The total mass of the medium is normalised such that the extinction optical depth for a configuration with k = 1 is similar to that for the model with dust. As we vary k, the extinction optical depth varies as well.

The panels in the middle column of Fig. 5 show the average surface brightness calculated by SKIRT along a central horizontal slice for each of the two viewing angles. Results are shown for optical properties in three regimes, respectively with dominating absorption (green, k = 4), no absorption (orange, k = 0), and net stimulated emission (blue, k = −1.5). As expected, the surface brightness increases as k decreases, and the features of the curves trace the spiral arms in the model. These effects are more prominent for the near edge-on view because this line of sight has longer paths through the medium. For the first two regimes (green and orange), the calculations were performed without and with explicit absorption, shown as solid and dashed curves, respectively. It is apparent from the figure that the results are essentially identical. The regime with net stimulated emission (blue) requires explicit absorption to be enabled, so only this result can be shown. Even if no direct comparison is possible, the progression seen in the curves for the three regimes does offer comfort that this latter result is correct as well.

The right panels ofFig. 5 show the mean intensity of the radiation field inside each of these model variations, averaged along the line of sight and along a horizontal slice of the view. These results confirm that the implementation with explicit absorption also properly records the radiation field.

All of the results shown in Fig. 5 were calculated with forced scattering enabled because of the need for recording the radiation field during the simulation. We separately confirmed, however, that the PLC without forced scattering also reproduces the surface brightness results of the middle panel.

Fig. 5 SKIRT results for a smooth spiral galaxy model with varying optical medium properties, viewed at inclinations of 20° (top row) and 80° (bottom row). The left column shows the surface brightness assuming a dust medium at V-band wavelengths, using an arbitrary colour scale. The middle column shows the average surface brightness along a central horizontal slice, assuming the model is observed at a distance of 10 Mpc. Results are shown for three types of media replacing the dust (green, orange and blue colours) and calculated without (solid curves) or with (dashed curves) explicit absorption. See Sect. 5.3 for details. The right column shows the mean intensity of the radiation field inside each of these model variations, averaged along the line of sight (weighted by medium mass) and averaged along a horizontal slice of the view.

5.4 Published benchmarks

In previous work, we employed several benchmark problems that have been published in the literature to test SKIRT with dust media (Camps & Baes 2020) and for Lyman-α line transfer (Camps et al. 2021). For the current work, to verify our updated implementation, we ran the benchmark problems listed below with explicit absorption enabled.

• Ivezić et al. (1997): a spherically symmetric (1D) circumstel-lar dust shell with various dust density profiles and optical depths.

• Pascucci et al. (2004): an axisymmetric (2D) circumstellar dust disk with various edge-on optical depths.

• Pinte et al. (2009): a disk geometry similar to that in the previous item, now including the effects of polarisation for anisotropic scattering by spherical dust grains.

• Gordon et al. (2017): a uniform dust slab (3D) externally illuminated by a star, testing dust extinction, heating, and emission for various optical depths of the slab.

• Dijkstra et al. (2006): analytical approximation for the spectrum emerging from a static, uniform hydrogen sphere with a central point source emitting at the Lyman-α line centre.

• Tasitsiomi (2006): numerical solutions for static, expanding and contracting hydrogen spheres with a central Lyman-α point source or with uniform emission throughout the sphere.

The dust benchmarks were run with forced scattering and the Lyman-α benchmarks without forced scattering. All results are within the uncertainty levels noted in previous work (Camps & Baes 2020; Camps et al. 2021). Even if all materials in these benchmark models have non-negative cross sections, we can still conclude that the explicit absorption technique is very robust, at least in this regime.

6 Efficiency

6.1 The figure of merit (FOM)

We noted in Sect. 2.5 that the explicit absorption technique should reduce Monte Carlo noise because it eliminates a stochastic element (the absorption/scattering split) but, on the other hand, requires the calculation of both the absorption and the scattering optical depth along the photon packet path. We now investigate how these differences influence the overall performance of the algorithm. To accomplish this, we need a mechanism to quantify and measure such performance.

Camps & Baes (2018) discuss quantitative statistical tests for the convergence of MCRT results, borrowed from the manual for the Monte Carlo N-Particle transport code MCNP (X-5 Monte Carlo Team 2003). We refer to these works and references therein for details. The method entails tracking the moments ${\displaystyle {\sum }_i{w}_i^k,k=1,2,3,4}$ of the individual photon packet contributions w_i to each result bin during the simulation. From these moments and the total number N of photon packets launched, two dimen-sionless quantities are calculated for each bin: the relative error, R, and the variance of the variance, VOV. If both values are below 0.1, the result can be considered to be converged and R can be interpreted as a true relative error on the result recorded for the bin. Moreover, one can then define the figure of merit, $$\text{FOM}=\frac{1}{R^{2}T},$$(37)

where T is the computer time spent on the simulation in seconds. Because the computer time is proportional to the number of photon packets N and the square of the relative error should scale as 1/N, the FOM is approximately constant once the solution has converged. The FOM thus provides a measure for the efficiency of a particular method for solving the problem at hand. Specifically, a higher FOM indicates that the simulation can reach a given level of signal-to-noise in less computer time.

SKIRT allows tracking the information underlying these statistics for the observed fluxes (Camps & Baes 2020). We construct a basic slab model similar to the geometry used by Gordon et al. (2017) and Camps & Baes (2020) for benchmarks and tests. In our version, a monochromatic source illuminates the top of a uniform slab of material with configurable optical properties and we observe the radiation emerging from the slab on the other side. The instrument views a central slice of the slab and has 150 bins along its length (each bin covers the width of the slice). After each simulation, we calculate mean 〈R〉 and 〈VOV〉 statistics averaged over these bins. For a given configuration of material properties and PLC variant, we run simulations with a successively larger number of photon packets until 〈R〉 < 0.1 and 〈VOV〉 < 0.1. Then we run two extra simulations with three times more photon packets and use the results from these three converged simulations to obtain an average FOM value and an estimate of its uncertainty caused by remaining Monte Carlo noise.

The material in our slab has uniform density and optical properties. For the purpose of calculating the penetrating flux we could therefore limit the spatial discretisation to a single grid cell enclosing the complete slab. However, we instead configure a regular grid with 75 cells in each of the horizontal directions and 100 cells in the vertical direction (the direction towards the source). This ensures that the calculation of the path geometries and optical depths for a fair number of cell segments are taken into account in the efficiency measure. The grid also allows the radiation field to be recorded with a fair spatial resolution.

This brings us to an important caveat. Our efficiency measure does not take into account the noise levels in the recorded radiation field, because SKIRT is not able to track the required statistics. It is possible that a given Monte Carlo method produces high-quality observed fluxes but records a more noisy radiation field. This will not be captured by our FOM metric.

Because of the dependency on the simulation time T, FOM values cannot be compared between computer systems. FOM ratios, however, should be stable. On the other hand, different model setups (i.e. other than our slab) and different implementations of the PLC (i.e. in codes other than SKIRT) might change relative FOM values in unpredictable ways. In spite of this, we believe the overall trends noted in the following subsection to be fairly robust.

6.2 Comparing algorithm variants

In Fig. 6, we show the FOM as a function of transverse extinction optical depth for the slab model described in the previous subsection. The figure shows three values of the scattering albedo (panels from left to right) and several variants of the PLC (colours and line styles). Similar to our earlier approach for the spiral galaxy model, we configure a fixed value for C_sca and a multiple of this value for the absorption cross section, C_abs = kC_sca. The total mass of the medium is then normalised such that the extinction optical depth (taking into account both scattering and absorption) across the slab in the vertical direction equals the required value, ranging from τ_ext = 0.1 to 10. We arbitrarily configure a fairly strong preference for forward scattering, with 〈cos θ〉 = 0.5.

We first consider the previously implemented PLCs, i.e. those without explicit absorption (solid curves). As expected, forced scattering (solid green) substantially improves performance over the basic method (solid orange) in all albedo regimes for low optical depths (τ ≲ 0.5). For higher optical depths, the situation is much less clear. As long as absorption dominates over scattering (middle and right panels), the forced scattering technique can be considered the method of choice, except perhaps for a narrow but possibly important intermediate optical depth range (0.5 ≲ τ ≲ 3) where it is slower by up to a factor of two. However, when scattering dominates (left panel), forced scattering slows down the simulation by more than an order of magnitude for τ ≳ 5. This can be attributed at least in part to the large number of peel-off photon packets that need to be sent in this situation. It is thus wise to disable forced scattering for media of this type. The third solid curve (blue) shows the performance of the forced scattering PLC when storing the radiation field. As can be expected, this additional task slows down the simulation by a constant factor.

More interesting in the context of this paper is the performance of the explicit absorption technique (dashed curves). Comparing the two PLCs without forced scattering (dashed and solid orange curves), we conclude that explicit absorption tends to slow down the simulation for lower optical depths (τ ≲ 2) and accelerate it for higher optical depths. The performance differences rise with decreasing scattering albedo, reaching an acceleration by three orders of magnitude at τ ≈ 10 when absorption dominates (right panel). This remarkable improvement seems to be a direct result of taking absorption out of the probabilistic part of the PLC and handling it separately in a deterministic, explicit manner. When forced scattering is enabled (dashed and solid green curves, or equivalently, dashed and solid blue curves), the performance differences caused by explicit absorption are much smaller. Even for a very low scattering albedo (right panel), the acceleration is limited to less than a factor of two.

When reviewing the performance of all PLC variants, one notable feature stands out. In the regime of balanced scattering and absorption (middle panel), enabling explicit absorption for the PLC without forced scattering (orange curves) boosts performance for τ ≳ 3 beyond that of the PLC with forced scattering (green curves). A similar but somewhat weaker boost can be seen for a medium with low scattering albedo (right panel).

Because the relative performance of the PLC variants so strongly depends on the regime imposed by the medium properties, there does not seem to be a single optimal algorithm. The specific model geometry and many other factors (see the caveats discussed at the end of Sect. 6.1) will also further influence the comparison. Nevertheless, based on Fig. 6 and our analysis, we propose as an overall rule of thumb to enable forced scattering except for media that are strongly dominated by scattering and with optical depths τ ≳ 0.5, and to always enable explicit absorption.

With this rule the simulation may perform sub-optimally (by a factor of less than two) in the balanced scattering/absorption regime (middle panel). To get optimal performance, one would have to configure different options depending on optical depth (the solid and dashed orange curves cross over at τ ≈ 2). And evidently, explicit absorption would always be required if the medium exhibits net stimulated emission. While a generic rule of thumb should avoid such complexity, it might in some cases be worthwhile to measure the FOM and select the appropriate algorithm, for example, when processing a large set of similar models.

Fig. 6 Figure of merit (FOM, defined in Sect. 6.1) as a function of the transverse extinction optical depth of our slab model, τext = 0. 1, 0. 2, 0. 5, 1, 2, 3, 5, 8, 10. The simulations were run on a basic multi-core desktop computer. A higher FOM value indicates a more efficient algorithm. Results are given for several variants of the photon packet life cycle (PLC): colours indicate without (orange) or with (green) forced scattering, or with forced scattering and also recording the radiation field (blue). The three panels correspond to different values of the scattering albedo. The left panel (ω = 0.9) represents the strong scattering regime and is applicable to dust attenuation at UV wavelengths. The central panel (ω = 0.5) corresponds to a situation where absorption and scattering are equally important and corresponds, for example, to dust attenuation at optical wavelengths. The right panel (ω = 0.1) corresponds to weak scattering, as applicable to dust attenuation at infrared wavelengths or radiative transfer in molecular and atomic lines. Line styles indicate without (solid) or with (dashed) explicit absorption. Error bars indicate the estimated uncertainty on the FOM values caused by the remaining Monte Carlo noise in the converged simulations.

7 Summary and conclusions

Most 3D radiative transfer codes employ the probabilistic Monte Carlo method, which follows the flow of photon packets through the obscuring medium rather than attempting to directly solve the radiative transfer equation. In conjunction with appropriate variance reduction and acceleration techniques, the approach is very powerful and allows one to emulate a wide range of physics. However, the classic Monte Carlo radiative transfer (MCRT) photon packet life cycle (PLC) requires a non-negative extinction cross section of the medium at every location in the simulated domain. This prohibits the straightforward handling of media with net stimulated emission, where the energy emitted at a given wavelength dominates the energy absorbed at that wavelength, so that the effective cross section becomes negative.

In this paper, we presented the explicit absorption technique that allows an adjusted MCRT PLC to efficiently handle net stimulated emission. The new technique uses the scattering optical depth along a photon packet’s path to randomly select the next interaction location, as opposed to using the extinction (absorption plus scattering) optical depth, and offers a separate, deterministic treatment of the absorption. While the scattering cross section still must be non-negative, a requirement that is easily met, the absorption and extinction cross sections are allowed to be negative. As a result, the method can handle net stimulated emission without any problems.

After describing the explicit absorption technique in detail, we verified its correctness and applicability in several ways. We first formulated a simple 1D radiative transfer problem and provided analytic solutions in various regimes, including net absorption and net stimulated emission. We verified that a special-purpose MCRT code incorporating explicit absorption is indeed capable of recovering these solutions in all regimes (see Fig. 4). We then reported on the implementation of explicit absorption in SKIRT, a fully featured state-of-the-art 3D MCRT code. This turned out to be a rather straightforward endeavour. Specifically, explicit absorption was easily combined with the variance reduction and acceleration techniques already incorporated in SKIRT, including peel-off and forced scattering. We verified the operation of the SKIRT implementation in the net absorption regime through regression tests and published benchmarks, and in the net stimulated regime through the same ID problem mentioned earlier and using a smooth 3D spiral galaxy model (see Fig. 5).

In the final section, we studied the efficiency of the explicit absorption technique using a basic setup for observing the radiation penetrating a 3D slab of material with varying optical properties (see Fig. 6). We found that explicit absorption slightly accelerates the simulation in most regimes as long as forced scattering is appropriately enabled or disabled. Based on these results, we recommended to always include explicit absorption, even if it is not strictly necessary for simulations in the net absorption regime (see the rule of thumb at the end of Sect. 6.2). Future experience with the technique will tell whether these trends hold for more realistic and complex models as well.

Acknowledgements

K.M. is a PhD Fellow of the Flemish Fund for Scientific Research (FWO-Vlaanderen) and acknowledges the financial support provided through grant number 1169822N.

References

All Figures

	Fig. 1 Graphical illustration of the two-stream radiative transfer problem considered in Sect. 3.
In the text

Fig. 2 (Cabs, Csca) parameter space for the 1D radiative transfer model discussed in Sect. 3.3. The green shaded area in this plot corresponds to the region with physically viable solutions. The three major regimes (net absorption, weak net stimulated emission, and strong net stimulated emission) are indicated. The orange and red lines indicate the boundaries between these major regimes, and correspond to pure scattering and critical net stimulated emission, respectively. The coloured dots indicate the models shown explicitly in Fig. 3. The white area corresponds to non-physical solutions. The solid green line that indicates the boundary between physical and non-physical solutions is given by the Eqs. (30) and (34).

In the text

Fig. 3 Intensity of the radiation field in the positive (leftpanel) and negative (rightpanel) direction for the 1D radiative transfer model discussed in Sect. 3.3. All models have the same scattering cross section Csca = 1 and the different lines correspond to models with different absorption cross sections Cabs, ranging from 2 to −2.5 in steps of −0.25. The thick orange and red curves correspond to pure scattering and critical net stimulated emission, respectively.

In the text

Fig. 4 Monte Carlo simulation results of the 1D radiative transfer problem discussed in Sect. 3.3. The three rows correspond to three different variants of the MCRT method: the standard implementation (top row), the version with the absorption-scattering split (middle row), and the version with explicit absorption (bottom row). Solid lines are the analytical results, the dots are the results of the MCRT calculations.

In the text

Fig. 5 SKIRT results for a smooth spiral galaxy model with varying optical medium properties, viewed at inclinations of 20° (top row) and 80° (bottom row). The left column shows the surface brightness assuming a dust medium at V-band wavelengths, using an arbitrary colour scale. The middle column shows the average surface brightness along a central horizontal slice, assuming the model is observed at a distance of 10 Mpc. Results are shown for three types of media replacing the dust (green, orange and blue colours) and calculated without (solid curves) or with (dashed curves) explicit absorption. See Sect. 5.3 for details. The right column shows the mean intensity of the radiation field inside each of these model variations, averaged along the line of sight (weighted by medium mass) and averaged along a horizontal slice of the view.

In the text

Fig. 6 Figure of merit (FOM, defined in Sect. 6.1) as a function of the transverse extinction optical depth of our slab model, τext = 0. 1, 0. 2, 0. 5, 1, 2, 3, 5, 8, 10. The simulations were run on a basic multi-core desktop computer. A higher FOM value indicates a more efficient algorithm. Results are given for several variants of the photon packet life cycle (PLC): colours indicate without (orange) or with (green) forced scattering, or with forced scattering and also recording the radiation field (blue). The three panels correspond to different values of the scattering albedo. The left panel (ω = 0.9) represents the strong scattering regime and is applicable to dust attenuation at UV wavelengths. The central panel (ω = 0.5) corresponds to a situation where absorption and scattering are equally important and corresponds, for example, to dust attenuation at optical wavelengths. The right panel (ω = 0.1) corresponds to weak scattering, as applicable to dust attenuation at infrared wavelengths or radiative transfer in molecular and atomic lines. Line styles indicate without (solid) or with (dashed) explicit absorption. Error bars indicate the estimated uncertainty on the FOM values caused by the remaining Monte Carlo noise in the converged simulations.

In the text