Continuum radiative transfer

Understanding the properties and evolutionary stage of protoplanetary disks requires to be able to derive physical quantities from multi-wavelength and multi-technique observations of these disks. Because of the complexity of the radiative transfer, we must rely on sophisticated numerical codes to do so. In this chapter, we present the radiative transfer problem in dusty media and discuss why this represents a difficult numerical problem to solve. We briefly describe the various methods that can be used to solve the radiative transfer equation and discuss their relative merits and drawbacks.


Introduction
Dust grains represent only a small fraction of the mass of protoplanetary disks, of the order of 1 %.But they play a central role in the disk thermal structure because as they completely dominate the continuum opacities.At short wavelengths, dust grains efficiently absorb, scatter, and polarise the starlight while at longer wavelengths dust re-emits the absorbed radiation.How much radiation is scattered and absorbed is a function of both the geometry of the circumstellar environment and the properties of the dust.In turn, the amount of absorbed radiation sets the temperature of the dust (and gas) and defines the amount of radiation that is re-emitted at longer, thermal wavelengths.
With the advent of high-angular resolution and high-contrast instruments, the basic structural properties (e.g., size, inclination, and surface brightness) of the circumstellar environments of the nearest and/or largest objects -disks and envelopes around young stars in nearby star-forming regions and around more distant evolved stars -are now under close scrutiny.With this unprecedented wealth of high-resolution data, from optical to radio, detailed studies of the dust content become possible and sophisticated radiative transfer (RT) codes are needed to fully exploit the data.
Because circumstellar disks remain optically thick over a broad wavelength range, a given observation only probes a very limited part of the disk.Figure 1 illustrates the various regions of a protoplanetary disk probed using different observational techniques.To obtain a global picture of disks, it is critical to combine all these observations and to model them simultaneously.

Observables and radiative quantities
For an astronomer looking at a source through a telescope, there are 4 observables for radiation: energy flux, direction, frequency or wavelength and the polarization.

Radiative flux
The radiative flux is defined as the energy going through an elemental surface per units of time per unit of surface area.In CGS units this has the dimension of erg s −1 cm −2 .
The flux F is called the bolometric flux, but most of the time we are interested in measuring the flux in a specific wavelength or frequency range, and we can defined a monochromatic flux F λ or F ν .The bolometric and monochromatic fluxes are related by Note that because dλ/ν = −λ/dν, F λ F ν but λF λ = νF ν .Several units are used in astronomy for the flux.In particular, it is often defined in Jansky, where 1 Jy = 10 −23 erg s −1 cm −2 Hz −1 .Note also that the flux is a vectorial quantity.The flux F we have defined is the component of the flux vector − → F that is perpendicular to dA: where − → n is a unitary vector normal to dA.

Specific intensity
The flux is a measure of the energy carried by all the photons, i.e. from all directions, passing through a surface.In radiative transfer theory, we generally work with a connected quantity, the specific intensity, which describes the amount of energy carried in a given direction.To define the specific intensity accurately in a given direction, we construct an infinitesimal area dA normal to the chosen direction and we consider all the rays whose direction is within a solid angle dΩ of the chosen direction.The energy crossing dA, in a time interval dt and a frequency range dν is given by where I ν is the specific intensity.I ν is a 6 dimension quantity, it depends on the frequency, the 3 spatial dimensions and 2 angular directions.It has units of erg s −1 cm −2 Hz −1 ster −1 .The specific intensity is connected to the flux by the following relation: Note that if I ν is isotropic, the corresponding net flux is zero.

Moments of specific intensity
Because the specific intensity contains an extremely rich information, it is often sufficient to work with quantities that are angularly averaged.In particular, it is convenient to define the zeroth, first and second tensor moment of the radiation field: The zeroth moment J ν is called the mean intensity and is the angular average of I ν .This quantity is of particular interest as it determines the heating (of dust and gas), the ionization and population levels.The first moment is proportional to the flux and describes the net flux of energy, while the second moment describes the radiation pressure.
To completely describe the radiation field, an infinite number of moments is required, but in most cases it is a reasonable approximation to only use the first 3 moments, together with an closure relation.
Summer School "Protoplanetary Disks: Theory and Modeling Meet Observations‰

Specific intensity is constant along a ray in free space
One of the fundamental properties of the specific intensity is that, if we are following a ray in free space (i.e.vacuum), the intensity in the direction of the ray remains constant along that ray: where ds is a differential element of length along the ray.Note that this is not in conflict with the fact that the flux F ν from an astronomical object is inversely proportional to the distance squared.The intensity is a quantity defined per unit of solid angle, while the flux is not.Because the angular scale of an object scales as ΔΩ ∝ 1/r 2 , so the does the flux.

Absorption and emission
When a ray passes through matter, energy can be added or removed and the specific intensity will not remain constant along the ray.
The energy emitted by a volume dV is defined by the spontaneous emission coefficient or emissivity j ν : dE = j ν dV dΩ dt dν (9) where j ν has dimension of erg s −1 cm −3 Hz −1 ster −1 .When moving a distance ds, a beam of cross section dA travels through a volume dV = dsdA and the intensity added to the beam is then: The loss of energy in a beam traveling a distance ds is described by where α ν is called the extinction coefficient and has units of cm −1 .Eq. 11 means that the intensity along a ray decreases as the exponential of the absorption coefficient integrated along the line of sight: This equation can easily understood if you consider a simple model of randomly distributed particles with a particle density n, each of the particles having a cross section σ ν (in cm 2 ).If you now consider a beam of radiation of surface dA along a distance ds, the number of absorbing particles seen by the beam is ndAds.The energy absorbed by the particles is: dE = −dI ν dA dΩ dt dν = I ν (n σ ν dA ds) dΩ dt dν (13) and then dI ν = −n σ ν I ν ds ( 14) The opacity can also be defined per unit of mass: κ ν = α ν /ρ (in cm 2 g −1 ), where ρ is the local density in g cm −3 .Note that depending if ρ is the gas or dust density, the opacity will be defined per gram of gas or per gram of dust.
The optical depth between 2 points s 1 and s 2 can be defined as: The extinction coefficient is connected to the photon mean free path, i.e. the mean distance a photon can travel through a medium before being absorbed or scattered.The mean optical depth is given by With τ ν = α ν l ν = 1, we find that the mean free path is l ν = 1/α ν .The optical depth τ ν of a medium can thus be understood as the number of mean free paths the photon travels through the medium before been absorbed.A medium is called optically thick at a given frequency if τ ν > 1 and optically thin if τ ν < 1.

The radiative transfer equation and its formal solution
If we include the effects of emission and absorption in the differential equation for the specific intensity along a ray, we obtain : This equation can be formally integrated between 2 points s 1 and s 2 giving:

Kirchhoff's law
Suppose you have a cavity filled with an emitting material of temperature T , in thermal equilibrium with the cavity.Then, the intensity should be everywhere I ν = B ν (T ).If we write the radiative transfer equation, we obtain and then This relation is called Kirchhoff's law.It means that the emissivity of a medium in thermal equilibrium is defined only by its temperature and the absorption coefficient α ν .It applies everywhere the medium is in local thermodynamic equilibrium (LTE).

Source function
As a generalization of Kirchhoff's law, we can define the source function S ν as the ratio of the emissivity by the absorption coefficient: allowing us to rewrite the radiative transfer equation: In the case of LTE, the source function is equal to the Planck function S ν = B ν (T ).Eq. 24 is very interesting in order to understand the behaviour of the radiation inside a medium, independently of the shape of S ν .The intensity always tends to approach exponentially the source function on spatial scales of the order of the mean free path.If S ν is constant along the ray, this appears clearly when integrating formally the radiative transfer equation: EPJ Web of Conferences 00006-p.6

Remark on time dependence
So far we have assumed that the light propagation is much faster than any timescale at which the object we are interested in can change.This is usually the case when studying disks, but this is not always true, as physical phenomena can occur on timescales of hours or minutes close to the star, i.e. shorter than the time needed for the light to reach the outer disk.
In those cases, the steady state radiative transfer equation must be replaced by the time dependent equation: For an example of numerical implementation and discussion of the effects of the time dependence of the radiative transfer in disks, we refer the reader to Harries (2011).

Radiative transfer in dust
Until now, all the equations and discussion were valid for radiative transfer in any medium, including line transfer in molecular and atomic gas and continuum transfer in dusty medium.From now, we will focus only on the continuum radiative transfer in this chapter, and refer the reader to Kamp (2015) for a discussion of line radiative transfer.

Temperature of a dust grain
Until now, we have assumed we knew the temperature of the dust grains.In practice, the temperature of the dust grains is set by the radiation field and the dust temperature and radiative transfer must be solved self-consistently.
To determine the temperature of a dust grain in equilibrium with the radiation field, we must compute the heating and cooling rates.The heating rate is determined by adding all the energy absorbed per second from the photons at all frequencies: Assuming the dust grain has reached a temperature T , the cooling rate is defined by the energy emitted by this grain per second: In radiative equilibrium, we have the thermal balance Q + = Q − , which leads to this equation fixing the temperature of the dust grain: It is interesting to notice that for a given radiation field J ν , the temperature of a dust grain is set by the shape of the opacity law.In particular, in a radiation field dominated by the stellar radiation, the dust temperature is mainly set by the ratio of opacities at the wavelengths of the peak of the photosphere and at the peak of the dust emission.For instance, a small sub-micronic grain has a steep opacity law between the optical and mid-infrared regimes (see Min 2015).They are not good emitters in the mid-infrared.To re-emit the absorbed energy, they will then need to reach high temperatures.Larger, millimetre-sized, grains on the contrary will have a flat opacity law.They will absorb less stellar radiation and re-emit the absorbed energy more easily in the infrared.As a consequence, they will have a much lower temperature than the small grains.The temperature decreases with distance as the stellar flux is reduced by 1/r 2 .The temperature roughly decreases with grain size, as the grain becomes a better emitter.For grains much larger than the wavelength of emission, the temperature becomes constant as the opacity becomes Grey over the wavelength range where there is significant absorption and emission.

Non-equilibrium dust grains
For very small grains in a diluted radiation field, like PAHs in the outer regions of a disk for instance, it is impossible to define a temperature.Small grains have small heat capacities and the absorption of a single UV photon can increase the grain temperature by a substantial amount.These small grains never reach an equilibrium temperature but undergo temperature fluctuations: they heat up very quickly when they absorb a photon and cool down progressively between 2 absorptions.As a consequence they emit on a much broader wavelength range than the one they would emit on if they were at equilibrium.For such transiently heated dust grains, the emissivity (Eq.21) can be written: where P(T ) is the probability density of the dust grain to be at a temperature T. Note that in practice, the integral stops at a maximum temperature above which the dust grain is sublimated.This probability density depends on the dust grain opacity and heat capacity as well as the hardness and dilution of the radiation field.Several methods exist to compute this probability density (e.g.Desert et al. 1986;Guhathakurta & Draine 1989;Leger & Puget 1984;Siebenmorgen & Kruegel 1992).

Scattering
Including light scattering strongly increases the complexity of the radiative transfer.Scattering is simultaneously a sink and a source term.Energy is scattered away from the beam, while energy is also scattered into the beam from all the other directions.
The radiative transfer equation becomes an integro-differential equation which relates the radiation fields at all positions and in all different directions: where α ext = α abs ν + α sca ν , α abs ν and α sca ν are the extinction, absorption and scattering cross section respectively.The ratio of the scattering and extinction cross sections is called the albedo : ) is the scattering phase function, which describes the anisotropy of scattering.It gives the probability that a photon incoming from the direction − → n will be scattered towards the direction − → n .This means that we need to solve the radiative transfer equation along all rays at the same time, as we need to know I ν (s, − → n ) for all the directions − → n in order to solve I ν (s, − → n ).

Radiative transfer of polarised light
The specific intensity only describes unpolarised light.Scattering on dust grains produces polarisation and aligned dust grains result in polarisation by emission and dichroïc absorption.Several formalisms exist to describe polarised light, one of the most commonly used uses the Stokes vector S = (I, Q, U, V), where I is the total specific intensity, Q and U describe the linearly polarised intensity and V the circularly polarised intensity.The radiative transfer equation can be easily extended in this formalism to account for polarisation: The modification from the non-polarised radiative transfer equation are two-fold: • the absorption and extinction coefficient are becoming vectorial to account for the polarised emission and dicroïc extinction from preferentially aligned grains • the phase function is replaced by the Mueller matrix (see a full description in Bohren & Huffman 1983;van de Hulst 1981 andMin 2015).

Why is radiative transfer difficult ?
In an optically thin medium, for instance a debris disk, solving the radiative transfer and dust temperature is easy.The radiation field is dominated at any point by the stellar radiation field, which also sets the dust temperature.When the medium is optically thick however, like in a protoplanetary disk, optical depths effects strongly increase the complexity of the problem: In some part of the volume, the stellar radiation will be strongly extinct and the specific intensity can be dominant, depending on the position and wavelength, either by the scattering term or the thermal emission term, which both depend on the radiation field at any other points in the volume.
We can summarize this extra complexity by rewriting the previously described equations: and Summer School "Protoplanetary Disks: Theory and Modeling Meet Observations‰ The intensity at any point I ν (s) depends on the source function S ν (s), which in turn depends on the specific intensity in any other direction, and coming from any other point in the model, both for emission and scattering source function.Equations 34 and 35 must be solved simultaneously.

Methods to solve the continuum radiative transfer
In this section, we briefly list the various methods that have been used to solve the continuum radiative transfer.We focus more on the Monte Carlo method as it seems to be become used more and more widely thanks to its flexibility and the development of new algorithms strongly improving its efficiency.
There are 2 main methods to solve the continuum radiative transfer, which differ in the way they deal with the angular dependence of the radiation field.Discrete ordinate (or ray-tracing) methods use a set of fixed directions on which they solve the radiative transfer equation, and Monte Carlo randomly samples all the directions by following photon packages.
The last decade has seen considerable progress in the development of RT techniques, going from 1+1D models (Chiang & Goldreich 1997), with only vertical transport of radiation and often grey opacities, to more complex models, 2 or 3D with more realistic opacities and treatment of scattered light.Monte Carlo (MC) codes, thanks to their high flexibility, seem to become the rule in the modeling of continuum emission, in particular in the field of protoplanetary disks.Only a couple of raybased have been developed to study disks: RADICAL (Dullemond et al. 2002), ProDiMo (Woitke et al. 2009) HO-CHUNK (Whitney & Hartmann 1992), while there is now a relatively large number of MC-based radiative transfer codes in this field: MC3D (Wolf 2003), TORUS (Harries et al. 2004), MOCASSIN (Ercolano et al. 2005), MCFOST (Pinte et al. 2009(Pinte et al. , 2006)), MCMax (Min et al. 2009).

Lambda iteration
The Lambda iteration is a conceptually simple method to solve the radiative transfer.Basically, it iterates between a global calculation (the radiative transfer equation) and a local calculation (the calculation of the source function) until convergence is reached.
The formal solution of the radiative transfer equation 34 can be written via the linear Λ operator: As we discussed it before, the formal solution of the RT equation and hence the Λ operator is not a local operation but involves integral over the entire volume.The simple idea between the Λ iteration is to iterate between equations 35 and 38: first, a guess on the source function is made, the formal radiative transfer equation is solved on a set of rays, allowing to compute J ν from S ν1 , finally a new source function at any point is computed from the local specific intensity, and the whole procedure is iterated.
Each iteration transports the radiation further by about a mean photon free path, meaning that the Λ iteration requires N iter τ 2 to reach convergence.
To overcome this difficulty, several accelerations can be used, in particular the Accelerated Λ Iteration (ALI, Cannon 1973;Collison & Fix 1991;Hubeny 2003) which splits the Λ operator into an approximate Λ * operator, which is easy or easier to invert, and the difference between the true and approximate Λ operators: A typical choice for Λ * is the diagonal or tridiagonal of Λ.This allows to separate the self-coupling of the cells with themselves (diagonal operator) and with their close neighbors (tridiagonal operator), which is now solved by a matrix inversion, from the global calculation.This results in a dramatic speed-up of the convergence of the iteration procedure.Furthermore, ALI methods can be coupled with additional acceleration such as the Ngacceleration (Dullemond & Dominik 2004;Ng 1974).

Moment method
In section 2.3, we have introduced the first 3 moments of the specific intensity.If we integrate the radiative equation over all the angles and then multiply by the directional vector and integrate again, we obtain the 2 first moment equations: These equations represent a non-closed system, with more unknowns than equations and we would need an infinite number of equations to end up with a system equivalent to the full angle-dependent radiative transfer equation.In practice, this system is truncated by introducing a closure equation, for instance the Eddington approximation: which is a good approximation for an isotropic radiation field.The equation system can then be reduced to one equation, which is easy to solve numerically with standard numerical methods: This is the diffusion approximation method.The Variable Eddington Tensor (VET) method (Mihalas & Mihalas 1984) generalizes such an approach, by computing the coefficients f i, j,ν in rather than fixing them.The calculation of the coefficients f i, j,ν is done by an procedure iterating the calculations of J ν and f i, j,ν and converges much more rapidly than an ALI method.Such a method is used by the code RADICAL (Dullemond & Dominik 2004).It remains extremely difficult to implement in 3 dimensions and for anisotropic scattering, which has limited its use in the field of protoplanetary disks so far.

Monte Carlo method
An alternative method to solve the radiative transfer problem is called the Monte Carlo method.Monte Carlo methods are numerical algorithms which can simulate random processes (Markov chains, transport model, finance, . . . ) but also completely deterministic problems, by building a probabilistic model.It is for instance well suited for multi-dimensional integrals.The generally accepted birthdate for this method is 1949 when Metropolis and Ulam published an article "The Monte Carlo" method.
In the case of radiative transfer, the basic idea is to propagate many photon packets by randomly sampling from probability distribution functions for directions, frequencies, path lengths, interactions.Essentially, a Monte Carlo radiative transfer code generates packets and follows them individually, one interaction event after another until they exit the system and reach a detector.All the parameters of the photon random walk are obtained by drawing random numbers following probability distribution function that reflects the local properties (optical depth, albedo, scattering phase function, temperature, . . .).This procedure is repeated for a large number of photon packets until the desired convergence is reached.
The Monte Carlo scheme estimates physical quantities by statistical means, which potentially leads to noisy results when the number of packets sampling some regions and/or directions in the model becomes low.Several variance reduction techniques have been developed to improve the sampling of the Monte Carlo method: forced first scattering (Cashwell & Everett 1959), peel-off techniques (Yusef-Zadeh et al. 1984), estimation of the mean specific intensity (Lucy 1999), immediate re-emission (Bjorkman & Wood 2001), polychromatic packets (Ercolano et al. 2005), importance weighting schemes (Juvela 2005;Pinte et al. 2006), partial diffusion approximation and modified random walk (Min et al. 2009), estimation of the scattering source function and combination with ray-tracing (Min et al. 2009;Pinte et al. 2009).These various techniques allowed the Monte Carlo method to progressively become competitive against grid-based methods, and it is now more and more commonly used to solve the continuum RT problem.
The Monte Carlo has also the strong advantage that it mimics the propagation of photons (Fig. 4), making them intrinsically 3D, and thus making it easier to implement additional physics, whereas the complexity of a discrete ordinate code increases dramatically with dimension.

How to choose random variables from any probability distribution function ?
To compute the various events of the packet's random walk, one needs to be able to sample an eventually complex probability distribution function.One way to do so is to use the fundamental equation of the Monte Carlo method.If a set of random variables A i follows a uniform distribution, the random variables X i = F −1 (A 1 ) follow the distribution f , with One of the main quantities to compute in a Monte Carlo is the distance that a packet will travel up to its next interaction with a dust grain.The fraction of packets that will interact between τ and τ + dτ is given by P(τ) dτ = e −τ − e −(τ+dτ) dτ.The probability distribution for the optical depth τ is then given by: The optical depth to the next interaction by is then: The free path l until the next interaction is then obtained by integrating: This equation is solved by propagating the packet in the medium and adding the opacity along the path until the optical depth τ is reached.This is usually the most computing intensive part of a Monte Carlo code and this integration must be performed with care to insure accurate numerical results.For most probability distributions, an analytic integration, like Eq. 47, is usually not possible and a numerical integration of Eq. 45 is performed instead as illustrated in Fig. 5.

Drawbacks of the basic Monte Carlo method
In principle, the Monte Carlo method is very simple and easy to implement, but it can become extremely slow in a variety of conditions: Summer School "Protoplanetary Disks: Theory and Modeling Meet Observations‰ The difference can be approximated as a differential as δT will be small if a large enough number of packets is used.This procedure is illustrated in Figure 6.The general mathematical basis for the Bjorkman & Wood (2001) calculation are discussed in details in Baes et al. (2005).

Computing the mean intensity
In the basic Monte Carlo method, the thermal balance in a given cell is performed by equating the energy to re-emit to sum of the energies of all the absorbed packets in the cell.
A much more effective method, introduced by Lucy (1999), is to directly use Eq. 29 by computing the mean intensity from the Monte Carlo run: where γ indicates a packet, γ is the packet luminosity and Δl γ is the length the packet propagated through in the cell i.Each packet crossing a cell is then contributing to the local mean intensity, as illustrated in Fig. 7.In particular, this allows to very efficiently compute the temperature in optically thin cells where absorption is very unlikely to happen.

Modified random-walk and partial diffusion approximation
In very optically thick regions, solving the radiative transfer via a Monte Carlo method becomes extremely challenging.Photon packets can interact many times (several millions for the optical depths encountered in protoplanetary disks) before exiting the system, increasing the computational times by several orders of magnitude.Min et al. (2009) introduced two complementary methods to speed up the calculations in high optical depth regions: the modified random walk (MRW), simplified by Robitaille (2010), and the partial diffusion approximation (PDA).The idea behind the MRW is to replace several thousands or even millions of interactions by a single propagation of the packet using the solution of diffusion approximation.Let's consider a packet "trapped" in a very optical thick cell.If the optical depth to any edge of the cell is much larger than unity, we can set up a sphere of constant density and temperature where the packet will propogate Summer School "Protoplanetary Disks: Theory and Modeling Meet Observations‰ Figure 7. Information created by a packet during its propagation.Left: In the basic Monte Carlo method, information is created only in cells where the packet interacts with the medium.Right: using the concept of specific and/or mean intensity allows creating information in all the cells crossed by the packet, resulting in a much lower number of packets required to reach a given convergence, in particular in the optically thin regions of the model.following a random walk.Using the diffusion approximation, we can then move the packet to the edge of the sphere in a single step.The actual distance propagated by the packet4 can be sampled from a probability distribution, which depends on the radius of the sphere and on the diffusion coefficient.The photon is then re-emitted from a random position on the sphere, at a wavelength sampled from the source function, i.e. from the Planck function (as we are in a very optically thick region and LTE can be assumed), or from its temperature derivative if the Bjorkman & Wood (2001) method is used.The procedure is repeated until the packet gets close enough to one edge of the cell, and a classical Monte Carlo propagation step can be used.
The PDA can additionally be used after the Monte Carlo run to compute the temperature structure in cells where a small number of packets have travelled.The temperature structure can then be solved using: ∇.(D∇ ) = 0 (52) where D = 1/3ρκ R is the local diffusion coefficient with κ R the local Rosseland opacity, and = 4σT 4 /c is the energy density, or more simply: This is a system of linear equation which can be solved very quickly, using the Monte Carlo temperature as a boundary condition.Note that continuum observations usually do not require this step, as it is used in regions which are not seen by the observer.The main goal of the PDA is to to obtain a reliable temperature structure at any point in the disk, for instance to compute the vertical hydrostatic equilibrium or the chemistry in the disk midplane.

Scattered specific intensity and ray-tracing
The basic Monte Carlo builds images and SEDs by putting detectors in a given direction.The probability of a packet to reach this detector is small.This is compensated by running large (a few 10 to 100 millions) packets to ensure that enough of them will reach the observer.But for a specific configuration and wavelength, for instance an optically edge-on disk in the mid-IR, the probability of a packet reaching the observer is so small that it would require prohibitively long CPU-time to compute an emerging flux.
To avoid this, ray-tracing has for instance been used in combination with Monte Carlo methods to produce SEDs and emission maps in the infrared and millimetre regimes, where scattering can be considered isotropic in some cases (Dullemond & Dominik 2004;Wolf 2003).When the scattering is isotropic, only the mean specific intensity and temperature structure are required to calculate the source function.These two quantities can be easily estimated with a Monte Carlo method (Lucy 1999) as discussed before and do not require a large amount of memory to be stored.After an initial Monte Carlo run computing the total mean intensity and temperature structure in the disk, SEDs and/or maps can then be produced by integrating the source function of rays originating from the observer.Which such a method, the Monte Carlo method is only used to estimate the specific intensity, but not the images and/ or SEDs.The resulting noise in the observables is much lower as it only reflects the noise in the mean specific intensity.
This method of combining Monte Carlo and ray-tracing can be extended to any wavelength if the angular dependence of the scattering component of the source function is preserved in the calculations (Min et al. 2009;Pinte et al. 2009).The Monte Carlo method produces all the information needed to perform such calculations, as it can give an estimate of the specific intensity, with its complete angular dependence, and not only of the mean specific intensity.
As a packet crosses a cell, we add its contribution to the specific intensity in a given direction for the given cell.Note that we need the full angular dependence of the specific intensity for each cell, i.e.I ν (x, y, z, θ, φ).This requires a large amount of memory and is only possible in practice for models with an azimuthal symmetry.For 3D models, an alternative is to set in advance the directions from where you wish to observe the disk, and store the scattered emissivity instead of the specific intensity.This makes the problem much more tractable but requires to rerun a Monte Carlo calculation for each new direction.
In both cases, at the end of the Monte Carlo calculations, i.e. once the specific intensity or scattered emissivity is known, as well as the dust temperature in any cell of the disk, it is possible to formally integrate the radiative transfer equation along the rays originating from the observer.

Conclusion
The continuum radiative transfer represent a complex numerical problem, in particular in protoplanetary disks where extremely high optical depths can be reached.The Monte Carlo method has become the method of choice to solve this problem in the last few years.With the development of recent algorithms, modern radiative transfer codes have become highly flexible, extremely robust numerically, and reasonably fast, allowing the exploration of a significant fraction of the parameter space when comparing model with observations.

Figure 1 .
Figure 1.Left: origin of emission in the disk as a function of wavelength (Figure from R. Lachaume).Right: Schematic view of the complementarity between observations at various wavelengths.Various observing techniques allow probing distinct parts of the disk.In anti clock-wise direction from top right: scattered light images probe the disk surface, going deeper as the wavelength increases, PAHs' emission probes a very superficial layer excited by stellar UV photons, near-IR interferometry probes the very central parts of the disk, mid-IR spectroscopy allows studying the properties of the grains up to a few 10s of AUs, and observations at millimeter wavelengths allow probing the bulk of the disk mass.

Figure 2 .
Figure 2. Shape of the radiation field J ν as a function of the position in the disk.The central plot shows the dust temperature structure of a typical disk surrounding a T Tauri star.Left column: the radiation field is dominated by the stellar radiation in the central parts, with a thermal component at longer wavelengths.At a few AUs in the midplane, the stellar radiation is completely attenuated and only the thermal emission remains (bottom right).At large distances, the radiation field is dominated by the diluted stellar light + thermal emission of the disk (top right), or the scattered light + thermal emission (middle right).Figure kindly provided by P. Woitke.

Figure 3 .
Figure3.Temperature of a single dust grain as a function of grain size and distance.The temperature decreases with distance as the stellar flux is reduced by 1/r 2 .The temperature roughly decreases with grain size, as the grain becomes a better emitter.For grains much larger than the wavelength of emission, the temperature becomes constant as the opacity becomes Grey over the wavelength range where there is significant absorption and emission.

Figure 4 .
Figure 4. Schematic view of the path of photon packets in a MC code.Packets perform a random walked with absorption and scattering events and eventually escape from the circumstellar medium.Full line represents the path of a stellar packet while the dashed line displays the path of a thermal packet emitted by the disk.

Figure 5 .
Figure5.Illustration of the sampling of a probability distribution in a MC code.Left panel: phase function to be sampled.Note that this is the phase of a single grain computed using the Mie theory.This is not a realistic representation of a phase function of dust grains in a disk and is just chosen for illustration here.Right panel: integrated probability distribution function.A random number A is chosen, the corresponding scattering angle of the value θ for which θ 0 S 11 (θ ) sin θ dθ / π 0 S 11 (θ ) sin θ dθ = A is obtained.

Figure 6 .
Figure6.Selection of the wavelength in the immediate re-emission algorithm.Emissivity is plotted before (lower curve) and after the absorption of an individual packet (upper curve).The spectral distribution of the previously distributed packets is represented by the white area.To correct the spectral distribution, the new packet must be emitted following the difference of the 2 curves, i.e. following the shaded area.