© A. Bensberg and S. Wolf 2022

1 Introduction

Variability is known as a characteristic feature of young stellar objects (YSOs) since their first observations (Joy 1945). The study of these time-dependent phenomena can reveal different underlying physical mechanisms such as spots on the stellar surface caused by magnetic fields, variable mass accretion onto the star, and obscuration by circumstellar material along the line of sight to the observer (Herbst et al. 1994). Since all of these mechanisms have an impact on the measured flux of the YSO, observed light curves of variable YSOs show different shapes ranging from periodical to linear or curved (Park et al. 2021). The time scales of the underlying event can be on the order of hours to years (e.g., Grankin et al. 2007; Cody & Hillenbrand 2014; Günther et al. 2014). In the case of the extreme accretion events that can be found in FUor and EXor objects, the outburst can even last decades (e.g., Hartmann & Kenyon 1996; Herbig 2008). The variability of the central object does also affect the illumination and thus the heating of the dust of the circumstellar disk. Therefore, variability can also be studied in the scattered light and traced in the thermal reemission radiation of the surrounding dust (e.g., Carpenter et al. 2001; Johnstone et al. 2013).

Radiative transfer simulations play an important role in the interpretation of astronomical observations. Because of its flexibility, the Monte Carlo method is often applied to simulate the radiative transfer in different astrophysical objects (e.g., Wolf et al. 1999; Harries 2000; Pinte et al. 2006; Dullemond et al. 2012; Ober et al. 2015). However, an underlying assumption in most of these solutions is that of radiative equilibrium. Moreover, while there is an implementation of Monte Carlo radiative transfer that uses gas properties to approach time-dependent heating and cooling (Harries 2011), there is currently no solution available focusing on the thermal properties of the dust.

In this paper, we present an implementation of an algorithm for full time-dependent 3D Monte Carlo radiative transfer (MCRT), which includes simulations of the scattered light, dust thermal reemission radiation, and temperature distribution with a focus on the latter. For this purpose, we extend the publicly available 3D Monte Carlo radiative transfer code POLARIS¹ (Reissl et al. 2016), which is used for the simulation of various astrophysical objects ranging from galaxies (Pellegrini et al. 2020; Reissl et al. 2020) and Bok globules (Brauer et al. 2016; Zielinski et al. 2021) to circumstellar disks (Brauer et al. 2019; Brunngräber & Wolf 2020) and exoplanetary atmospheres (Lietzow et al. 2021).

This paper is organized as follows: The computational method is described in Sect. 2. In Sect. 2.1, we focus on the theoretical background of the relevant thermal processes. The basic algorithm for the time-dependent simulation of temperature distributions is presented in Sect. 2.2. The procedures for the calculation of high-resolution images of the thermal reemission radiation and scattered light, that is to say algorithms for time-dependent ray tracing as well as time-dependent scattering, are briefly described in Sects. 2.3 and 2.4. An overview of different tests of the implemented routines is given in Sect. 3. The impact of the temporal step width and the number of photon packages per time step is discussed in Sects. 3.1 and 3.2. The influence of the optical depth on the temperature distribution simulations for 1D and 2D models are examined in Sects. 3.3 and 4.1, respectively. In Sect. 4.2, we discuss the impact of a luminosity outburst of a central star embedded in a circumstellar dust disk on the resulting dust temperature distribution, as well as on the scattered light and thermal reemission radiation. Finally, a summary as well as an outlook on further applications and potential improvements of the method are provided in Sect. 5.

2 Computational method

In the following, we outline the theoretical background of time-dependent continuum radiative transfer. After discussing the thermal processes and properties of the dust in Sect. 2.1, we describe the implementation of a time-dependent method of radiative transfer simulations for temperature, dust reemission and scattering in Sect. 2.2. For this, we extend the publicly available 3D Monte Carlo radiative transfer code POLARIS (Reissl et al. 2016).

2.1 Description of heating and cooling

We first focus on the thermal processes and corresponding properties of the dust. The ith photon package with wavelength λ emitted by a radiation source with luminosity L that is emitting N_ph photon packages in the time-interval Δt carries the energy $$e_{i,\lambda }=\frac{L\text{Δ}t}{N_{\text{ph}}}.$$(1)

Along its path through the model space, the dust absorbs energy of the photon package in relation to the wavelength-dependent mass absorption coefficient κ_λ. The total absorption rate $\dot{A}$ in a region with mean radiation intensity J_λ is then given by: $$\dot{A}=4\pi {\displaystyle {\int }_0^{\infty }{\kappa }_{\lambda }{J}_{\lambda }\text{d}\lambda \text{.}}$$(2)

Following Lucy (1999), the mean intensity J_λ of a discrete number of photon packages with energy e_i,λ within a wavelength range dλ traveling a path length δl_i in a volume V in the time-interval Δt can be estimated by: $$J_{\lambda }\text{d}\lambda =\frac{1}{4\pi }\frac{1}{V\text{Δ}t}{\displaystyle \sum _i{e}_{i,\lambda }\delta l_i}.$$(3)

An estimator for the absorption rate in Eq. (2) is then derived by substitution of Eq. (3) and discretization: $$\dot{A}=\frac{1}{V\text{Δ}t}{\displaystyle \sum _i{\kappa }_{\lambda }{e}_{i,\lambda }\delta l_i}.$$(4)

Due to the absorption process the dust will heat up to a temperature T, thus the enthalpy u(T) of the dust grain will increase. According to Guhathakurta & Draine (1989) the enthalpy u(T) of an N atom dust grain with volume V_d and heat capacity per volume C_v can be calculated via: $$u\left(T\right)=\left(1-2\text{/}N\right)V_{\text{d}}{\displaystyle {\int }_0^T{C}_{\text{v}}\left(T\right)\text{d}T}.$$(5)

Under the assumption of local thermodynamic equilibrium (LTE), the dust will emit radiation with an emission rate $\dot{E}$ corresponding to its emission spectrum B_λ at temperature T scaled with the absorption coefficient κ_λ: $$\dot{E}\left(T\right)=4\pi {\displaystyle {\int }_0^{\infty }{\kappa }_{\lambda }{B}_{\lambda }\left(T\right)\text{d}\lambda \text{.}}$$(6)

Considering the difference between absorption and emission rate within one time step, that is combining Eqs. (4) and (6), we now can estimate the enthalpy u_n+1 of the dust after n + 1 time steps with step width Δt: $$u_{n+1}=u_n+\left(\dot{A}-\dot{E}\right)\text{Δ}t.$$(7)

However, this holds only true if we assume the absorption and emission rate is constant over the step width Δt. We can constrain the temporal step width by considering the time it takes to emit the total enthalpy of a dust grain at temperature T, that is the cooling time t_c(T). For the case of emission without absorption, we find $$t_{\text{c}}\left(T\right)=\frac{u\left(T\right)}{\dot{E}\left(T\right)}$$(8)

as upper limit for the step width. Since both emission rate and enthalpy depend on the size of the grain, the cooling time is size-dependent. An exemplary plot of the cooling time of astronomical silicate (see, Sect. 3) with different grain sizes is shown in Fig. 1. It can be seen that the cooling times of grains with sizes up to 1.5 µm are of the same order of magnitude. Larger grains are cooling slower, because they are less efficient emitters due to a smaller ratio of their emitting surface to their volume. In both cases, the cooling time is rapidly decreasing with increasing temperature. The choice of the temporal step width Δt is thus crucial for the goodness of the estimators described in Eqs. (4) and (7). A detailed discussion of the corresponding numerical consequences can be found in Sect. 4.

Fig. 1 Cooling time described in Eq. (8) for different temperatures of grains with different radii consisting of astronomical silicate.

2.2 Numerical implementation

Based on the concepts described above, we implemented a time-dependent continuum radiative transfer algorithm in our 3D MCRT code POLARIS. Therefore, the temporal step width is added as parameter to the stationary radiative transfer scheme for the simulation of the temperature distribution (hereafter: stationary method). A detailed description of the stationary scheme can be found for example in Wolf et al. (1999) and Reissl et al. (2016). POLARIS allows the use of several grid geometries to build the model space including cylindrical, spherical, octree, and voronoi. The number of grid cells as well as their spacing can be set individually for every dimension. Multiple stellar sources can be set in the grid at any position. However, for the sake of simplicity we are considering the case of a single central star (spherical, isotropic irradiation) for the following description. A schematic flowchart of the whole procedure for a single time step is shown in Fig. 2. The algorithm for the simulation of the time-dependent temperature distribution can be divided into a series of four steps that are repeated every time step until the total simulation time, that is the time at which the simulation is stopped, is reached.

Fig. 2 Flowchart of the numerical procedure for the simulation of time-dependent temperature distributions described in Sect. 2.2 (the four main steps outlined there are color coded: 1 – orange, 2 – green, 3 – blue and 4 – red).

1 Emission from dust or star

At the beginning of each time step, a fixed number of photon packages N_ph is emitted from the star or the dusty environment. Therefore, we have to determine the luminosity of the dust of every cell. This is done by reading out the current temperature of the cell from the grid and using Eq. (6) to calculate the respective emission rate of a single dust grain. We can now calculate the luminosity of the ith cell by multiplying the emission rate ${\dot{E}}_i$ of the cell with its number density of the dust grains n_i and volume $V_i:L_i={\dot{E}}_i{V}_i{n}_i$. The total dust luminosity L_tot is then obtained by summing up the luminosity of all cells. The ratio of total dust luminosity to the combined luminosity of star and dust gives the probability of a photon package to be emitted from the corresponding sources. In the case of thermal emission from within the dust distribution, the emitting cell is determined by the cumulative probability density function p_i for the ith cell given by: $$p_i=\left({\displaystyle \sum _{k=0}^i{L}_k}\right)\text{/}{L}_{\text{tot}}.$$(9)

Once the source is determined, a random direction of emission is drawn and a wavelength is assigned to the photon package. The respective cumulative probability density function $p_{{\lambda }_i}$ of a wavelength distribution with bin size Δλ and a source luminosity $L_{{\lambda }_i}$ at wavelength λ_i is given by: $$p_{{\lambda }_i}=\left({\displaystyle {\int }_{{\lambda }_i}^{{\lambda }_i+\text{Δ}\lambda }{L}_{{\lambda }_i}}\right)\text{/}{L}_i.$$(10)

Next the energy of the photon package is assigned using Eq. (1) and the first point of interaction at optical depth τ is determined by using a random number ζ ∈ [0, 1) (see also Reissl et al. 2016): $$\tau =-\log \left(1-\zeta \right).$$(11)

Lastly, the photon package is stored by adding it to a stack (hereafter: photon stack) that is carried over from one time step to another. The emission process of a photon package can be divided into the following steps:

1. Determination of the source of emission by using the probability density function determined by the luminosity ratio of the emitting sources.

2. Calculation of the emitting cell in the case of thermal reemission by dust using Eq. (9).

3. Computation of the random direction of emission.

4. Determination of the wavelength using the cumulative. probability density function of the emission spectrum of each source.

5. Calculation of the energy of the photon packages according to Eq. (1).

6. Computation of the optical depth to first interaction via Eq. (11).

7. Storing of the photon package (photon stack).

This procedure is repeated until all N_ph photon packages are emitted.

It should be noted that the emission of photon packages with equally distributed energy and wavelength determined by the cumulative probability density function of the emission spectrum discussed above will overestimate the energy at the most probable wavelength if the number of photons is not high enough to be statistically significant. While this will affect the local radiation field, it will not significantly affect the resulting dust temperature distribution, since the local deviations average out on a global scale. However, to make sure that the local radiation field is reproduced correctly, steps 4 and 5 can be replaced by another emission method: Instead of drawing the wavelength of photon packages using a probabilistic approach, a fixed number of photon packages N_λ is emitted for every wavelength bin and the energy is assigned using the respective value of the emission spectrum. The energy e_*,j of a photon package with wavelength λ for a stellar source with radius R_* and effective temperature T is then given by: $$e_{*,j}=4{\pi }^2{R}_{*}^2\frac{B_{\lambda }\left(T\right)\text{Δ}t}{N_{\lambda }}.$$(12)

In the case of thermal reemission radiation by the dust, the emitting cells are determined by the luminosity ratio described by the probability density function in Eq. (9). This is done to save computation time by neglecting the thermal reemission of dust with low temperatures that will not significantly contribute to the radiation field. The wavelength distribution of the emitted photon packages depends on the temperature of the dust of the corresponding cell. Since photon packages will only be emitted from cells with high dust luminosity, it must be ensured that the total luminosity L_tot of the dust of all cells is conserved. To do this for the emission of photon packages of every wavelength of the chosen cells, Eq. (1) is modified by scaling with the emission spectrum ${\tilde{B}}_{\lambda }\left(T\right)$ of the respective cell. Here, the emission spectrum is normalized such that the integral over the wavelength is equal to 1. The energy e_{dust, j} of a photon package emitted by the dust is then calculated via: $$e_{\text{dust,}j}=\frac{L_{\text{tot}}\text{Δ}t}{N_{\lambda }}·{\tilde{B}}_{\lambda }\left(T\right).$$(13)

However, this method is computationally expensive because N_λ photon packages need to be emitted per wavelength, emitting source and time step. Thus, it is only used if the detailed, local radiation field is of interest.

2 Continuous absorption

In the next step, all photon packages of the photon stack are propagated through the grid until either the optical depth of interaction (see, Eq. (11)) or the end of the time step is reached. This is determined by integrating the path length and the optical depth of each photon package. While traveling through a cell, a photon package will deposit energy that adds to the estimated absorption rate of the cell, which is determined using Eq. (4). The photon package is not loosing energy on its way through the grid, that is, the absorbed energy is not subtracted. However, there are other approaches to handle absorption that deliver equivalent results, but require higher computational times. A more detailed discussion of these different approaches can be found in Appendix A.

3 Scattering or absorption

When the optical depth of interaction determined by Eq. (11) is reached, the probability for scattering and absorption is determined by the albedo of the dust grains. The type of interaction is then obtained by using a random number drawn from a uniform distribution. Since the process of choosing the type of interaction and the treatment of scattering are not different from the stationary case, we refer to Reissl et al. (2016) for a detailed description. In the case of an absorption event, the photon package is marked for deletion without adding additional energy to the absorbing cell. If the photon package is scattered, a new optical depth for interaction is calculated using Eq. (11) and the photon package is propagated further through the grid, that is steps 2 and 3 are repeated until the total path length of the photon package reaches the maximum path length corresponding to the temporal step width. Photon packages leaving the grid within the time step will also be marked for deletion.

4 Update of enthalpy and temperature

At the end of a time step, the photon packages marked for deletion are removed from the stack and the temperature of each cell is updated. Therefore, the new enthalpy u(T) of the dust of each cell is calculated by evaluating Eq. (7). The corresponding temperature is determined by interpolating precalculated values of the enthalpy using Eq. (5). After the temperatures are updated, the algorithm proceeds with the next time step.

2.3 Time-dependent ray tracing

The algorithm presented in Sect. 2.2 allows producing images of the scattered stellar radiation and the thermal reemission radiation of the dust. However, the simulation of high resolution images using this algorithm is computationally very expensive since every photon package needs to be stored on and loaded from the memory in every time step. Therefore, we implemented two fast, time-dependent methods to calculate scattering images and dust emission maps. For the latter, a time-dependent ray tracing algorithm is used. The time-dependent algorithm for the simulation of images of the scattered light is presented in Sect. 2.4.

The ray tracer implemented in POLARIS follows the path of parallel rays through a model space, which is divided into multiple grid cells. A detector with a two-dimensional array of pixels is set outside of the model space. The single rays are sent out from points on the opposite side of the computational grid, opposite to the detector pixels at which they are to be observed. The intensity of the ray at the start is either zero or set by the background radiation. When crossing a grid cell, the intensity of a ray is reduced corresponding to the optical depth τ of the dust in the cell by a factor e^−τ and increased by the thermal reemission radiation of the dust. This is done by solving the corresponding radiative transfer equation using the Runge–Kutta–Fehlberg method of order 4(5). When all cells along the path are processed, the total intensity of the ray is saved on the corresponding detector pixel. To ensure the applicability of the method for high optical depth, a recursive refinement for the integration step as well as the detector pixel size is applied (see, Ober et al. 2015).

The time-dependent ray tracer is following this procedure, but saves the intensity of the rays for different time steps. This is done by stopping the calculation of the intensity along the ray after a path length δl corresponding to the speed of light c and the step width of a time step Δt = δl/c. The current total intensity of each ray is then scaled with the remaining optical depth towards the observer and projected onto the detector. For the next time step, the detector is replaced by a new detector on which the total intensity of the rays after the new time step is saved, and the intensity represented by the rays is again obtained by using the Runge–Kutta–Fehlberg method. This procedure is repeated for as many time steps as needed to account for the contributions of all volume elements along the path of each ray.

It should be noted that the temperature and thus the thermal reemission of each cell is changing with time. Therefore, one ray tracing simulation is only showing the contribution of the thermal reemission radiation of the temperature distribution of one time step. Consequently, every image calculated with the ray tracer only contains the contribution of the dust with the temperature distribution for that particular time step. Thus, some of the resulting images are showing the thermal reemission radiation at the same time steps, e.g., the first image simulated using the second temperature distribution shows the dust emission of the same time step as the second image of the simulation using the first temperature distribution. To get the final images, the images corresponding to the same time steps are added up. A sketch of the complete time-dependent ray tracing procedure for two different temperature distributions is shown in Fig. 3. A simple test of the corresponding light travel time effects for a spherical dust distribution can be found in Appendix B.

2.4 Time-dependent scattering

As long as the dust density distribution of a model is not changed, for example, due to dust sublimation, the distribution of the scattered stellar light does not depend on the temperature of the dust. Time-dependent scattering simulations can thus be performed independently of the temperature simulations. For this purpose, we modified the stationary Monte Carlo method for scattering simulations with peel-off first presented in Yusef-Zadeh et al. (1984).

Here, the basic handling of stellar emission and scattering are similar to the procedure described for the simulation of time-dependent temperature distributions in Sect. 2.2: For the wavelength of interest, a fixed number of photons packages are emitted from the source (e.g., star). The energy of the photon packages is assigned according to the spectrum of that source (e.g., Eq. (12) for a stellar source). Next, the random directions of the photon packages and the optical depth of the first interaction are determined. The photon packages are then propagated through the model space until the optical depth of the first interaction is reached. In the case of a scattering event, a copy of the respective photon package (hereafter: peel-off photon) is sent to the detector. The energy of the peel-off photon is diminished in relation to the optical depth as well as the scattering probability to the observer. The scattered photon package is then propagated to the next optical depth of interaction. If a photon package leaves the model space, its energy is stored on a detector.

For the time-dependent scattering simulations, the path length of each photon package within the model space is integrated. In addition, the different distances to the detector due to different viewing angles have to be considered. This is done by adding the remaining distance to the total path length of a photon package after it leaves the model space. Peel-off photons inherit the path length of the interacting photon package at the point of interaction and can thus be treated equally. The single detector used for the stationary simulations is replaced by multiple detectors which cover different intervals in time. The energy of a photon package leaving the model space will thus be saved on the detector that covers the time interval corresponding to the total path length of the photon package. Thus, using the speed of light c, we determine the index n of a detector corresponding to the time step t_n with step width Δt for a photon package with total path length l via: t/c ∈ [t_n: t_n + Δt].

Changes of the properties of the emitting source are treated likewise by defining multiple sources which are emitting at different time steps. Since the time step currently assigned to a photon package is determined by the path length, the total path length of a photon package emitted by the ith source at time ti has an offset of l₀ = t_i · c. The procedure described above is then repeated until all photon packages of all sources left the model space.

Fig. 3 Sketch of the time-dependent ray tracing procedure for two different temperature distributions. The straight arrows denote the path of integration of each ray. The ray tracing simulation is performed based on the temperature distributions at time steps t0 and t1 separately. Detector images, i.e., images at the position of the observer obtained for the same time step are added up.

3 Tests

In the method described above, we introduced two parameters to handle time-dependent radiative transfer within the framework of an existing Monte Carlo approach, namely the temporal step width Δt and the number of photon packages per time step N_ph. In the following, we describe the influence of these parameters as well as the influence of the optical depth on the temperature distribution. This is done by simulating the heating process of a simple one-dimensional dust distribution around a centrally embedded illuminating and heating star for different values of each of the above parameters.

The dust mixture used for all models corresponds to the typical composition of the ISM with 62.5% astronomical silicate and 37.5% graphite, for which we apply the $\frac{1}{3}-\frac{2}{3}$ relation for parallel and perpendicular orientations (see, Draine & Malhotra 1993). The grain size distribution follows Mathis et al. (1977) and can thus be described by a power law n(s) ∝ s^−3.5 for grain sizes between 5 nm and 250 nm. The optical properties of the dust are derived using the wavelength-dependent refractive indices by Draine & Lee (1984) and Laor & Draine (1993). The calorimetric data is obtained by applying the models for silicate and graphite presented in Draine & Li (2001). To save computational time, we assume mean values for all dust grain parameters. Following Wolf (2003) they are weighted with respect to the abundances of the chemical components and the grain number density as well as the grain size distribution. To justify this approach, we compared temperature simulations using multiple and single grain sizes to the respective simulations using the stationary method. We found that the temperature distributions calculated with weighted mean values are deviating less than 1% from temperature distributions simulated with unaveraged values.

3.1 Temporal step width Δt

In Sect. 2.1, we have shown that the temporal step width is a crucial parameter for the goodness of the estimator of the enthalpy and thus for the simulation of the temperature distribution. To illustrate the effect of different values of this parameter, we simulated the temperature distribution of a simple one-dimensional dust distribution with constant density (inner radius R_in = 10 AU, outer radius R_out = 100 AU, optical depth τ_V = 0.1) and an embedded T Tauri-like star (hereafter: 1D model) and compared it to the resulting temperature distribution of the stationary method of POLARIS. Since the stellar parameters had been fixed during the simulation, the dust in the grid heated up until the equilibrium temperature was reached. The simulation was started with a dust temperature of 2.7 K for all cells and stopped when the temperature of the dust in each cell is no longer changing within a ΔT given by the bin size of the temperature grid (logarithmically distributed between 2.7 K and 2500 K). To ensure that the time interval is long enough to reach this equilibrium, the simulations covered at least twice the light traveling time through the grid. Thus, a total simulation time of 10⁵ s was chosen for the given model. To rule out the influence of the number of photon packages per time step on the simulations, we fixed N_ph at 10³ for all simulations. Thus, the total number of photon packages emitted over the simulation time for the smallest temporal step width was 10⁸. Therefore, this number of photon packages was used for the simulation of the stationary reference temperature distribution. An overview of all model parameters can be found in Table 1. The resulting temperature distributions of the simulations with step widths Δt = 1 s, 2 s, 4 s and 6 s are shown in Fig. 4.

While the temperature distributions obtained for the temporal step widths Δt = 1 s and 2 s converged to the stationary temperature distribution, the temperature distributions calculated with Δt = 4 s and 6 s are oscillating in the inner region. This can be explained by calculating the cooling time t_c for the cells with the highest dust temperature using Eq. (8). If the temporal step width is close to the cooling time, these cells will cool almost to the ground temperature of 2.7 K in a single time step. Thus, the estimated emission rate of the dust will be low compared to the absorption rate. Since absorption and emission rates are considered constant within one time step, this leads to an overestimation of the dust enthalpy given by Eq. (7) and thus an overestimated dust temperature. These cells with higher dust temperatures will cool even faster. The dust temperatures of the affected cells are then oscillating between very high and low values. Since the dust is emitting photon packages in every time step, and overestimated dust temperatures lead to an overestimated energy of the emitted photon packages, the oscillation propagates through the grid. The larger the temporal step width compared to the cooling time, the stronger the oscillation. For this reason, the step width should always be chosen according to the cooling time of the estimated maximum dust temperature of the simulation.

Fig. 4 Resulting temperature distribution (top) and relative difference (bottom) of the 1D model discussed in Sect. 3.1 with different temporal step widths At for the time-dependent (dashed) and stationary (solid) method. The number of photon packages per step is set to 103, the total number of photon packages used for the stationary simulation is 108, respectively. Lines for radii larger than 20 AU are overlapping and are therefore not shown.

3.2 Number of photons per time step N_ph

The number of photon packages is a crucial parameter for every numerical approach to the radiative transfer problem that makes use of the Monte Carlo method. In the case of time-dependent simulations, this has to be taken into account for each individual time step. The effect of an increase of the number of photon packages on the computational time of each simulation is stronger than in the stationary case because all photon packages need to be stored and loaded for every time step. Furthermore, the total number of photon packages processed in every time step will increase until the number of photon packages leaving the grid equals the number of photon packages that are emitted. The computation time of a single time step is thus increasing linearly. It is therefore crucial to evaluate the effect of the number of photon packages per time step on the resulting temperature distribution.

For this purpose, we performed a series of simulations for the simple 1D model presented in Sect. 3.1 with a fixed temporal step width of Δt = 2 s and different number of photon packages per time step N_ph = 10, 10², 10³, 10⁴. The temperature distribution obtained with the stationary approach using 10⁸ photon packages in Sect. 3.1 was used as reference. The resulting temperature distributions and relative differences are shown in Fig. 5. It can easily be seen, that the temperature distribution gets smoother for higher numbers of photon packages. Furthermore, a value of N_ph = 10 leads to an overestimation of the dust temperature in the inner cells.

Since the probability of an absorption event in the case of τ_V = 0.1 is sufficiently small, it can be assumed that most of the photon packages will not be absorbed, thus their energy will be carried through all cells until the photon package leaves the model space. However, some photon packages will undergo an absorption event, that is their further transfer through the grid will be stopped. Following the approach to handle absorption discussed in Sect. 2.2 a photon package is only contributing energy to the enthalpy of the dust by continuous absorption along its path through the grid. Thus, after an absorption event, the radiation field in the following cells is weakened. In the case of a very low number of photon packages per time step (N_ph < N_Cells) nearly no photon package will be absorbed and the whole energy emitted by the source in every time step contributes to the dust enthalpy of every cell. Thus, the dust temperature of each cell is overestimated. However, in the case of a single absorption event, the impact on the local radiation field will be larger (e.g., in the case of N_ph = 10 roughly 10% of the energy distributed in one time step is taken out), leading to large variations of the temperature distribution. In the case of higher optical depth, that is if absorption is more probable, the dust temperature will be underestimated. The same applies for the wavelength selection, which is also done by using a probability density function (see, Eq. (10)). The flux at wavelength where the emission probability is highest will be overestimated, while the flux at wavelength where the emission probability is lowest will be underestimated.

Fig. 5 Resulting temperature distribution (top) and relative difference (bottom) of the 1D model discussed in Sect. 3.2 with different number of photon packages per time step Nph, for the time-dependent (dashed) and stationary (solid) method. The temporal step width is set to 2 s, the total number of photon packages used for the stationary simulation is 108, respectively.

3.3 Optical depth

To examine the behavior of the optical depth on the time-dependent temperature simulation, we performed simulations of models with different optical depths (τ_V = 0.01, 0.1, 1, 10 as seen from the central star). We fixed the temporal step width and the number of photons per time step to values that lead to well converged temperature distributions in Sects. 3.1 and 3.2, that is, Δt = 2 s and N_ph = 10³. For every optical depth, a reference temperature distribution was simulated using the stationary method for temperature distributions. The resulting temperature distributions of the time-dependent simulations and the relative differences to the respective reference temperature distributions after a total simulation time of 10⁵ s are shown in Fig. 6. It can be seen that the temperature distributions of the models up to τ_V = 1 are in very good agreement with the temperature distributions of the stationary simulations. In the case of the τ_V = 10 model, the dust temperature is always lower than the reference temperature distribution. This is caused by a slight discrepancy between the optical and calorimetric data used in the time-dependent temperature simulations, as well as some numerical constraints that are further discussed in Sect. 5. However, the temperature distribution of the time-dependent simulation is deviating less than 8% from the respective temperature distribution of the stationary approach and is therefore considered consistent.

Fig. 6 Resulting temperature distributions (top) and relative differences (bottom) of the 1D model for different optical depths τV discussed in Sect. 3.3 for the time-dependent (dashed) and stationary (solid) method. In the case of time-dependency, the temporal step width is fixed at 2 s and the number of photon packages per time step is 103. The total number of photon packages used for the stationary simulations is 108.

4 Examples

Since the goal of our method is the simulation of young stellar objects with circumstellar disks, we now consider the radiative transfer of a model of a very low-mass circumstellar disk. The model description as well as a comparison between time-dependent and stationary simulation of the temperature distribution is given in Sect. 4.1. Time-dependent simulations of the scattered light as well as the thermal reemission radiation for a simple luminosity outburst of the stellar source are presented in Sect. 4.2.

4.1 Circumstellar disk model

The density distribution of the disk is following the considerations of Shakura & Sunyaev (1973) and has been successfully used for modeling circumstellar disks in various studies (e.g., Pinte et al. 2008; Ratzka et al. 2009; Brunngräber et al. 2016). It is given by: $$\rho \left(r,z\right)={\rho }_0{\left(\frac{r}{R_0}\right)}^{-\alpha }\exp \left[-\frac{1}{2}{\left(\frac{z}{h\left(r\right)}\right)}^2\right],$$(14)

with the cylindrical coordinates (r, z), the density scaling parameter ρ₀ determined by the total disk mass, the radial density profile parameter α and h(r) the scale height, which can be written as $$h\left(r\right)=h_{\text{ref}}{\left(\frac{r}{R_0}\right)}^{\beta }.$$(15)

Here, R₀ is the radius of the reference scale height h_ref. The parameter β determines the disk flaring. We applied this density distribution for a model with an inner radius of R_in = 20 AU, outer radius R_out = 100 AU and moderate flaring (β = 1.125 and h_ref = 10 AU).

To make sure that the two-dimensional density structure has no effect on the time-dependent temperature distribution, we performed simulations of the disk model with midplane optical depths of τ_V = 0.01, 0.1, 1, and 10. We then compared these to the temperature distributions simulated with the stationary method. The temporal step width (Δt = 8 s), number of photons per time step (N_ph = 10^s) and the total simulation time (10⁵ s) were chosen according to the criteria discussed in Sect. 3. Corresponding to the total number of photon packages used for the time-dependent temperature simulations, the number of photon packages applied for the stationary simulation of each model was 10⁹. An overview of all parameters of the disk model can be found in Table 2.

A selection of vertical cuts through the temperature distribution of the time-dependent simulation for an optical depth τ_V = 1 for different time steps, that is different stages of the heating process, is depicted in Fig. 7. The resulting temperature distributions and relative differences of the midplane temperature for all models are shown in Fig. 8. The trend of the temperature distribution with different optical depth of the 1D model discussed in Sect. 3.3 can also be found for the disk model. The time-dependent temperature distributions with optical depths of τ_V = 0.01, 0.1, 1 are in good agreement with the temperature distributions of the stationary simulations. Since the optical depth is highest in the midplane of the disk model, we can assume that the temperatures of the entire disk are equally well fitting. In the case of τ_V = 10, the temperature distribution in the midplane is deviating by about 5% from that in the midplane of the stationary simulation. To check whether this also affects the dust temperature of cells above and below the midplane, a plot of the difference of the dust temperature distribution between stationary and time-dependent simulation for all grid cells is given in Fig. 9. It can be seen that the deviation of the dust temperature of the cells above and below the midplane is always below 2% and thus comparable to the level of the statistical noise of the simulation.

To make sure that the corresponding spectral energy distribution (SED) is not significantly affected by the underestimated dust temperature in the midplane, we performed simulations of the scattered light and thermal reemission radiation for the time-dependent temperature distribution as well as for the reference temperature distribution. The resulting SEDs as well as the relative difference between the SEDs calculated using the time-dependent and the stationary temperature distribution are shown in Fig. 10.

While the overall shape of the SED is well reproduced, the difference plot reveals a peak difference of about 15% for wavelengths of around 20 µm. This is close to the shortest wavelength at which the thermal reemission radiation of the dust for the given temperature distribution becomes significant. The difference plots in Figs. 8 and 9 show that the dust temperature of the cells in the inner few AU of the midplane is underestimated. These cells correspond to the cells with the highest dust temperature determined with the stationary method. The spectrum of the thermal reemission radiation of the dust with the time-dependent temperature distribution is thus shifted to longer wavelengths compared to the spectrum calculated with the stationary method. However, this only affects the radiation of wavelengths corresponding to the thermal reemission of the few cells with the highest dust temperature. Therefore, the flux of most of the wavelengths larger than 10 µm is showing not more than 10% difference to the reference SED.

Fig. 7 Vertical cuts through the temperature distribution of the time-dependent disk model after 20 000, 40 000 and 60 000 s described in Sect. 4.1 with τ = 1 and temporal step width 8 s and 105 photons per time step.

Fig. 8 Resulting midplane temperature distributions (top) and relative differences (bottom) of the disk model with different optical depth τV discussed in Sect. 4.1 for the time-dependent (dashed) and stationary (solid) method. In the case of time-dependency, the temporal step width is fixed at 8 s and the number of photon packages per time step is 105. The total number of photon packages used for the stationary simulations is 108.

Fig. 9 Difference plot of a vertical cut through the dust temperature distribution of the disk model with τ = 10 described in Sect. 4.1. The reference temperature distribution is calculated with the stationary method using 108 photon packages.

Fig. 10 Normalized SED (top) and relative difference (bottom) of the time-dependent temperature distribution of the disk model with τ = 10 discussed in Sect. 4.1 (dotted). The reference SED (solid) is calculated using the temperature distribution of the stationary method.

4.2 2D disk outburst

Lastly, we performed simulations of a model for an outburst of a star with a circumstellar disk. For this, we used the disk model described in Sect. 4.1 (see also Table 2). For the sake of simplicity and to save computational time, we chose a low optical depth of τ_V = 0.01. In this case the dust temperature of the cells are converged after about 50 000 s, that is roughly the time light needs to travel once through the model space. The luminosity of the central star was then increased by a factor of 4 for a duration of 10 000 s. The total simulation time was set to 120 000 s to make sure, that all grid cells return to their previous dust temperature after the outburst.

The temperature simulation was performed using a temporal step width of 4 s and 10⁵ photon packages per time step. The temporal step width was chosen such that the cooling time of the expected highest dust temperature is not reached within one step. Three vertical cuts through the temperature distribution of the disk model for different time steps during the outburst are shown in Fig. 11. It can easily be seen that a heat wave corresponding to the outburst in the source luminosity is traveling through the grid in radial direction. Furthermore, the thermal response times of the dust are comparably small, that is the dust of a cell is heating or cooling to the equilibrium temperature within one time step, leading to a discrete form of the heat wave.

Furthermore, we investigated the scattered stellar light and thermal reemission radiation of the dust as observables of this system. In particular, we simulated the thermal reemission radiation of the N band (10 µm) and the scattered stellar light of the V band (550 nm) using the respective time-dependent algorithms for ray tracing and scattering described in Sect. 2.3 and Sect. 2.4. In the case that the observed disk is inclined, it is necessary to take the different light traveling times to the observer into account. To illustrate this effect, the intensity maps were derived for an inclination of 25° as well as for the face-on case (0°). Since the scattering and reemission simulations are not sensitive to the temporal step width, the step width was set to 1000 s and a total number of photon packages per time step of 10⁸ were used for scattering and reemission, respectively. In order to avoid an impact of the initial heating process, all simulations were started once a steady state was reached. To rule out an offset of the light curves of the scattered and reemitted radiation caused by additional light travel times due to different positioning of the detector in the scattering and ray tracing treatments, the start of the outburst was set to 60 000 s in both cases.

The resulting light curve for the scattered light is shown in Fig. 12. A sample of the corresponding scattering images for the inclined disk can be found in Fig. 13. There are three effects that can be seen in the light curves: First, the overall flux is higher in the case of the inclined disk. This is caused by the different probabilities of the scattering direction. The scattering of the stellar light towards the observer is more probable for the inclined than for the face-on disk. Secondly, the scattered light of the outburst is first visible in the light curve of the inclined disk. This shows the importance of the light traveling time toward the observer. Since the lower part of the disk in the images is closer to the observer, the outburst will be seen there first. This is also visible in the first image of Fig. 13. Furthermore, the outburst is stretched out compared to the face-on case, because parts of the disk with the same radial distance to the source will have different distances to the observer. This can be seen in the second of the scattering images. The third effect is a long afterglow caused by scattered photon packages with long distances to the observer. This wave of dim, scattered photon packages can be found in the upper half of the third scattering image.

The light curve of the simulation of the thermal reemission radiation is shown in Fig. 14. The corresponding images of the reemission radiation are given in Fig. 15. The main effect visible in the shape of the light curves is the light traveling time of the geometrical depth structure of the disk. While the heating and cooling time scales are short compared to the temporal step width, the heating of each layer of the disk in the reemission light curve will be visible with an offset relative to the light traveling time within the disk. Combined with the different light travel times due to the inclination of the disk, the light curve of the inclined disk shows a broader peak than the face-on light curve. It should be noted that the duration of the outburst relatively to the light travel times within the grid is a crucial factor for the shape of the emission peak, both in the case of the thermal reemission radiation and for the scattered radiation. There is also an afterglow visible in both light curves. However, the duration of the afterglow in the thermal reemission case is small compared to the one seen in the light curves of the scattered light. The reason for this is that the relative cold outer parts of the disk are not significantly contributing to the measured reemission radiation flux before or during the outburst. This also causes the emission peaks of the inclined and face-on disk to be closer together than in the case of scattering.

Fig. 11 Vertical cut through the temperature distribution of the outburst disk model after 65 000, 85 000 and 105 000 s described in Sect. 4.2. The temporal step width of the simulation is set to 4 s; 105 photon package are emitted per time step.

Fig. 12 Light curve of the scattered stellar light (550 nm) of the outburst model discussed in Sect. 4.2. The duration of the outburst is 10 000 s; the total simulation time is 120 000 s. During the outburst, the stellar luminosity is increased by a factor of 4. The flux is normalized with respect to the maximum of the integrated flux of all images of the different time steps.

4.3 2D disk outburst – High optical depth

To test the treatment of backwarming, that is the heating of regions of the disk by the thermal reemission radiation of different regions, we repeated the time-dependent temperature simulations of the disk model used in Sect. 4.2 with an optical depth of τ_V = 100 in the midplane. In this case, the outer parts of the midplane of the disk are shielded from the direct emission of the star and can thus only be heated by the thermal reemission of the surrounding dust. A vertical cut through the temperature distribution of the high optical depth disk model after 95 000 s during the outburst is shown in Fig. 16. It can be seen that there is a warm region in the midplane around 70 AU that is enclosed by colder regions. This means it was not heated by the stellar emission of the outburst, but by the thermal reemission of the dust above and below the midplane. It should also be noted, that the width and position of the warm spot do not correspond directly to the outburst in the optically thinner layers of the disk. This is caused by the delay between the heating by the stellar emission and the thermal reemission of the dust, as well as the density structure of the disk.

5 Discussion and conclusions

We presented an implementation of an algorithm for full time-dependent 3D Monte Carlo radiative transfer for the MCRT code POLARIS. It includes time-dependent methods for the simulation of temperature distributions and thermal reemission as well as scattering images. The method for the simulation of the temperature distribution focuses for the first time on the thermal properties of dust instead of gas. It can be used for the simulation of the radiative transfer of various astrophysical objects.

We discussed the influence of the newly introduced parameters temporal step width Δt and number of photon packages per time step N_ph. It was found that the temporal step width has to be small compared to the shortest cooling time of the systems, whilst the number of photon packages per time step has to be large enough for correct scattering and absorption treatment. Since all photon packages in the grid have to be loaded from and stored on a stack in each time step, the access speed and size of the CPU cache as well as the RAM are crucial parameters for the total computation time. Even though the treatment of the single photon packages in POLARIS is parallelized, the computation time is thus not necessarily decreasing linear with the number of threats. Temporal step width and number of photon packages per time step should therefore be chosen accordingly.

The time-dependent simulation of temperature distributions was successfully tested by simulating 1D and 2D models with different optical depth (τ_V = 0.01, 0.1, 1, 10, 100). A deviation from the reference temperature distributions calculated with the stationary approach of POLARIS was only seen for the case of τ_V = 10. The main reason for this is a slight inconsistency between the applied optical and calorimetric properties of the dust. While the optical properties are derived from measurements of crystalline olivine and fitted for astronomical purposes (see, Laor & Draine 1993), the thermal properties are derived from measurements of the heat capacity of basalt glass, obsidian glass and SiO₂ glass (see, Draine & Li 2001). A discrepancy between optical and thermal properties leads to an underestimation of the temperature, if the heat capacity is overestimated compared to the optical properties. Better fitting dust data could therefore improve our method.

Finally, we tested a possible application of the code by simulating an outburst of a star with a circumstellar disk. Here, the temperature distribution, scattering images and light curve as well as thermal reemission images and light curve were calculated. The scattering and emission images and light curves were simulated for a face-on disk as well as for a disk with an inclination of 25°. The time-dependent ray tracing as well as the time-dependent scattering simulation are fast and can be performed independently of the simulation of the temperature distribution.

The effect of the heating of the disk by reemission radiation of the dust was shown by repeating the temperature simulation of the outburst model with an optical depth of τ_V = 100 in the midplane.

It should be noted, that the presented method for the simulation of temperature distributions is most efficient at time scales from hours to weeks. To simulate variable events in circumstellar systems on longer time scales, it is needed to include the physics of the thermal coupling of the gas and dust phase. While this is out of the scope of this paper, the presented methods provide a valuable basis on which additional physical effects can be implemented. Since the methods for the time-dependent simulation of scattering and thermal reemission images are in principle independent of the temperature calculation, they can already be used on larger time scales. In the case of the thermal reemission simulations, the corresponding time-dependent temperature distributions could also be provided from hydrodynamics simulations.

Even though the presented method for the simulation of temperature distributions is in principle not dependent on the structure of the used dust distribution model, the computational afford increases with its complexity and optical depth. In order to represent the underlying optical depth structure for models with a high optical depth, an increased number of computational grid cells is needed. As discussed in Sect. 3.2 a high number of grid cells requires a higher number of photon packages per time step. Steep gradients in the dust density structure will also lead to high temperature gradients that require smaller temporal step widths. For example, in the case of a circumstellar disk with optical depth τ_V > 10^s the temperature simulation would require at minimum 10⁴ photon packages per time step with a temporal step with of around 0.01 s. Depending on the radial extent of the model, the total number of photon packages processed per time step are then about 10⁸ to 10⁹. This corresponds roughly to the computational afford of a single stationary radiative transfer simulation and should therefore be considered when setting up time-dependent simulations.

Future applications are aiming at the simulation of images and light curves of observed variable illumination and heating sources embedded in dust distributions, such as young stellar objects embedded in their parental molecular cloud or surrounded by circumstellar disks, to study their geometrical structure and dust properties. Potential improvements of the numerical implementation cover a faster memory management as well as the possibility to start the temperature simulation from a steady state instead of including the simulation of the initial heating process of the disk.

Fig. 13 Images of the scattered light of the outburst model after 65 000, 85 000 and 105 000 s for 550 nm described in Sect. 4.2 (temporal step width: 1000 s; 108 photon packages per time step). The flux is normalized with respect to the maximum flux of all images of the different time steps.

Fig. 14 Light curve of the thermal reemission radiation (10 µm) of the outburst model discussed in Sect. 4.2. The duration of the outburst is 10 000 s; the total simulation time is 120 000 s. During the outburst, the stellar luminosity is increased by a factor of 4. The flux is normalized with respect to the maximum of the integrated flux of all images of the different time steps.

Fig. 15 Images of the thermal reemission radiation of the outburst model after 65 000, 85 000 and 105 000 s for 10 µm described in Sect. 4.2 (temporal step width: 1000 s; 108 photon packages per time step). The flux is normalized with respect to the maximum flux of all images of the different time steps.

Fig. 16 Vertical cut through the temperature distribution of the outburst disk model with high optical depth after 95 000 s described in Sect. 4.3. The temporal step width of the simulation is set to 4s: 105 photon package are emitted per time step. The contour lines mark three levels of temperatures. The red arrow outlines a region in which the backwarming effect can be seen.

Appendix A Absorption estimator

Following the discussion of the probability density functions of absorption optical depth τ_abs and interaction optical depth τ_ext in Lucy (1999; further discussed in Krieger & Wolf 2020), we are subsequently focusing on the different ways to treat the absorption of a photon package without explicitly taking scattering into account. We consider the energy of a photon package deposited while passing an optical depth interval of [τ₁,τ₂). The absorption probability for this interval can be written as ${\displaystyle {\int }_{{\tau }_1}^{{\tau }_2}\text{d}\tau e^{-\tau }}$. For simplicity, we are normalizing the energy of a photon package to 1. It will be shown, that the same energy is deposited in all three approaches to handle absorption:

Case 1: No continuous absorption, only single absorption events In the case of only single absorption events, all energy is deposited in the grid and the photon package is deleted after the event. The energy ε deposited in the grid is then given by: $$\epsilon ={\displaystyle {\int }_{{\tau }_1}^{{\tau }_2}\text{d}\tau e^{-\tau }·1=e^{-{\tau }_1}-e^{-{\tau }_2}.}$$(A.1)

Case 2: Continuous absorption, photon package is loosing energy We are now considering a photon package that is continuously loosing energy which is deposited in the grid, that is ($e^{-{\tau }_1}-e^{-\tau }$). In the case of an absorption event, the photon package leaves all remaining energy in the grid and is deleted. The deposited energy can then be calculated with: $$\begin{array}{lll}\epsilon \hfill & =\hfill & {\displaystyle {\int }_{{\tau }_1}^{{\tau }_2}\text{d}\tau e^{-\tau }\left(\left(e^{-{\tau }_1}-e^{-\tau }\right)+\left(1-e^{-\tau }\right)\right)}\hfill \\ \hfill & \hfill & +{\displaystyle {\int }_{{\tau }_2}^{\infty }\text{d}\tau e^{-\tau }\left(e^{-{\tau }_1}-e^{-{\tau }_2}\right)=e^{-{\tau }_1}-e^{-{\tau }_2}.}\hfill \end{array}$$(A.2)

Case 3: Continuous absorption, photon energy stays constant The last case is similar to case 2, except that the energy of the photon package stays constant and no energy is deposited in the grid when the photon is deleted. Thus, a photon package traveling a distance δτ leaves the energy δτ in the grid: $$\epsilon ={\displaystyle {\int }_{{\tau }_1}^{{\tau }_2}\text{d}\tau e^{-\tau }\left(\tau -{\tau }_1\right)+{\displaystyle {\int }_{{\tau }_2}^{\infty }\text{d}\tau e^{-\tau }\left({\tau }_2-{\tau }_1\right)=e^{-{\tau }_1}-e^{-{\tau }_2}.}}$$(A.3)

It can be seen that the deposited energy is the same in all three cases. However, in a numerical simulation, we do not always reach this limit. In the simulations presented in this paper, we chose case 3, in which the energy is distributed evenly in the grid even for small numbers of photon packages.

Appendix B Exemplary study: Dust shell outburst model

In Sect. 2.3 we presented a time-dependent ray tracing method for the simulation of thermal reemission images. A simple test for the correct treatment of light travel times is an outburst of a central source embedded in a spherical dust distribution. We thus set up a one-dimensional spherical model containing five dust shells with constant dust density (inner radius R_in = 0.1 AU, outer radius R_out = 10 AU, optical depth τ_V = 55) and an initial dust temperature of 50 K. An overview of all model parameters can be found in Table B.1.

For the sake of simplicity, we model a heatwave by raising the dust temperature of every shell to 75 K, beginning from the innermost shell according to the light travel time through the entire model. After each time step, the dust temperature is reset to 50 K and the temperature of the next cell is set to 75 K. The temporal step width of Δt = 1000 s is chosen such that one time step corresponds to the light traveling time through one of the dust shells in radial direction.

The position of the reemitting, heated dust as seen from the observer at a time t after which the outburst is first detected can be described analytically by a parabola known for example from supernova light echoes (see, Wright 1980): $$z=\frac{x^2}{2ct}-\frac{1}{2}ct,$$(B.1)

where z is the distance to the observer, x the distance parallel to the observering plane, and c the speed of light. The corresponding parabolas for the outburst model described above are shown in Fig. B.1.

Fig. B.1 Dust shell outburst model: Regions of dust reemission radiation seen by an observer for different time steps after an outburst for the model described in Sect. B calculated using Eq. B.1. The dashed line indicates the outer radius of the dust distribution. The position of the central source is indicated by an asterisk.

Fig. B.2 Images of the thermal reemission radiation of the dust shell outburst model after 1000, 2000 and 3000 s at a wavelength of 10 µm as described in Sect. B (temporal step width: 1000 s). The flux is normalized with respect to the maximum flux of the image after 1 000 s. The white dotted lines indicate the maximum radius of the parabola depicted in Fig. B.1.

These analytic predictions can now be compared to the images of the thermal reemission radiation obtained using the time-dependent ray tracing algorithm described in Sect. 2.3. The thermal reemission images for three time steps after the first detection of the outburst can be found in Fig. B.2. The contour lines represent the maximum radius of the parabolas shown in Fig. B.1. The resulting images obtained with the time-dependent ray tracing algorithm match the analytic predictions perfectly.

References

All Tables

All Figures

	Fig. 1 Cooling time described in Eq. (8) for different temperatures of grains with different radii consisting of astronomical silicate.
In the text

	Fig. 2 Flowchart of the numerical procedure for the simulation of time-dependent temperature distributions described in Sect. 2.2 (the four main steps outlined there are color coded: 1 – orange, 2 – green, 3 – blue and 4 – red).
In the text

	Fig. 3 Sketch of the time-dependent ray tracing procedure for two different temperature distributions. The straight arrows denote the path of integration of each ray. The ray tracing simulation is performed based on the temperature distributions at time steps t0 and t1 separately. Detector images, i.e., images at the position of the observer obtained for the same time step are added up.
In the text

Fig. 4 Resulting temperature distribution (top) and relative difference (bottom) of the 1D model discussed in Sect. 3.1 with different temporal step widths At for the time-dependent (dashed) and stationary (solid) method. The number of photon packages per step is set to 103, the total number of photon packages used for the stationary simulation is 108, respectively. Lines for radii larger than 20 AU are overlapping and are therefore not shown.

In the text

	Fig. 5 Resulting temperature distribution (top) and relative difference (bottom) of the 1D model discussed in Sect. 3.2 with different number of photon packages per time step Nph, for the time-dependent (dashed) and stationary (solid) method. The temporal step width is set to 2 s, the total number of photon packages used for the stationary simulation is 108, respectively.
In the text

Fig. 6 Resulting temperature distributions (top) and relative differences (bottom) of the 1D model for different optical depths τV discussed in Sect. 3.3 for the time-dependent (dashed) and stationary (solid) method. In the case of time-dependency, the temporal step width is fixed at 2 s and the number of photon packages per time step is 103. The total number of photon packages used for the stationary simulations is 108.

In the text

	Fig. 7 Vertical cuts through the temperature distribution of the time-dependent disk model after 20 000, 40 000 and 60 000 s described in Sect. 4.1 with τ = 1 and temporal step width 8 s and 105 photons per time step.
In the text

Fig. 8 Resulting midplane temperature distributions (top) and relative differences (bottom) of the disk model with different optical depth τV discussed in Sect. 4.1 for the time-dependent (dashed) and stationary (solid) method. In the case of time-dependency, the temporal step width is fixed at 8 s and the number of photon packages per time step is 105. The total number of photon packages used for the stationary simulations is 108.

In the text

	Fig. 9 Difference plot of a vertical cut through the dust temperature distribution of the disk model with τ = 10 described in Sect. 4.1. The reference temperature distribution is calculated with the stationary method using 108 photon packages.
In the text

	Fig. 10 Normalized SED (top) and relative difference (bottom) of the time-dependent temperature distribution of the disk model with τ = 10 discussed in Sect. 4.1 (dotted). The reference SED (solid) is calculated using the temperature distribution of the stationary method.
In the text

	Fig. 11 Vertical cut through the temperature distribution of the outburst disk model after 65 000, 85 000 and 105 000 s described in Sect. 4.2. The temporal step width of the simulation is set to 4 s; 105 photon package are emitted per time step.
In the text

	Fig. 12 Light curve of the scattered stellar light (550 nm) of the outburst model discussed in Sect. 4.2. The duration of the outburst is 10 000 s; the total simulation time is 120 000 s. During the outburst, the stellar luminosity is increased by a factor of 4. The flux is normalized with respect to the maximum of the integrated flux of all images of the different time steps.
In the text

	Fig. 13 Images of the scattered light of the outburst model after 65 000, 85 000 and 105 000 s for 550 nm described in Sect. 4.2 (temporal step width: 1000 s; 108 photon packages per time step). The flux is normalized with respect to the maximum flux of all images of the different time steps.
In the text

	Fig. 14 Light curve of the thermal reemission radiation (10 µm) of the outburst model discussed in Sect. 4.2. The duration of the outburst is 10 000 s; the total simulation time is 120 000 s. During the outburst, the stellar luminosity is increased by a factor of 4. The flux is normalized with respect to the maximum of the integrated flux of all images of the different time steps.
In the text

	Fig. 15 Images of the thermal reemission radiation of the outburst model after 65 000, 85 000 and 105 000 s for 10 µm described in Sect. 4.2 (temporal step width: 1000 s; 108 photon packages per time step). The flux is normalized with respect to the maximum flux of all images of the different time steps.
In the text

	Fig. 16 Vertical cut through the temperature distribution of the outburst disk model with high optical depth after 95 000 s described in Sect. 4.3. The temporal step width of the simulation is set to 4s: 105 photon package are emitted per time step. The contour lines mark three levels of temperatures. The red arrow outlines a region in which the backwarming effect can be seen.
In the text

	Fig. B.1 Dust shell outburst model: Regions of dust reemission radiation seen by an observer for different time steps after an outburst for the model described in Sect. B calculated using Eq. B.1. The dashed line indicates the outer radius of the dust distribution. The position of the central source is indicated by an asterisk.
In the text

	Fig. B.2 Images of the thermal reemission radiation of the dust shell outburst model after 1000, 2000 and 3000 s at a wavelength of 10 µm as described in Sect. B (temporal step width: 1000 s). The flux is normalized with respect to the maximum flux of the image after 1 000 s. The white dotted lines indicate the maximum radius of the parabola depicted in Fig. B.1.
In the text