Need to Review?

Radiation has a very prominent role in astrophysics, both as the only observable and also as a dominant mechanism for energy transfer within and out of astrophysical systems. Consequently, its transport through a medium is one of the most fundamental processes that needs to be considered. Analyzing the radiation received from an object provides us with useful information not only about the radiative source itself but also the medium in between and surrounding the object and the observer. The radiative transfer problem is formally defined by the non-steady state radiative transfer equation (See lecture on radiative transport (Th, Sep. 5) for detailed description).

${\frac {dI_{\nu }(x,n,t)}{ds}}+{\frac {dI_{\nu }(x,n,t)}{cdt}}=j_{\nu }(n,t)-(\kappa _{\nu ,abs}(n,t)+\kappa _{\nu ,scat}(n,t))\rho (x,t)I_{\nu }(n,t)+\kappa _{\nu ,scat}(n,t)\rho (x,t)\int _{4\pi }\Phi (n,n^{/}prime,x,\nu )I_{\nu }(x,n^{/}prime)d\Omega ^{/}prime\,\!$ where s is the distance along the path defined by a position x and propagation direction n and $\kappa (s,\lambda )$ is the mass extinction coefficient ($\alpha (s,\lambda )=\kappa (s,\lambda )\rho (s)$ ). $\Phi (n,n^{/}prime,x,\nu )$ is the scattering phase function which describes the probability that a photon originally propagating in direction n’ and scattered at position x, will have n as its new propagation direction after the scattering event.

The challenging aspect of radiative transfer is that it is a six-dimensional problem (3 spatial + 2 directions ($\theta$ ,$\phi$ ) + frequency) with additional time dependence. Additionally, significant complexity is introduced as a consequence of the ability of radiation to affect the state of the medium which is the source of the radiation itself. Thus, to know the state of the matter, we need to know the radiation field. To know the radiation field, we need to know the photon emission, absorption and scattering rates. However, to estimate these rates, we in turn again need to know the state of the matter. Furthermore, the RT problem is nonlocal in space (photons propagate within the entire domain of interest), direction (scattering and absorption/re-emission can change the direction of photons), and wavelength (absorption at a particular wavelength gets re-emitted at another wavelength). As a consequence of these challenges, radiative transfer in astrophysics is frequently calculated using radiative transfer codes which can also account for 3D geometry and non-linear affects due to dust properties (See Lecture on dust grains).

A numerical algorithm for integrating the formal transfer equation

Note : To understand what follows, a basic knowledge of numerical methods is useful. If you are completely new to numerical methods, "Numerical Recipes" by Press, Teukolsky, Vetterling and Flannery is a good reference source.

The primary aspects of numerical radiative transfer can be illustrated using a 1D plane parallel atmosphere. In this model, the gas density, gas temperature etc in the atmosphere is assumed to depend only on height (z). Consequently, we only require one directional variable ($\mu =cos(\theta )$ ), with the steady state Radiative Transfer Equation (RTE) reducing to the following form (overall 3D problem - spatial + direction + wavelength) -

$\mu {\frac {I_{\nu }(z,\mu )}{dz}}=j_{\nu }(z)-\alpha _{\nu }(z)I_{\nu }(z,\mu )$ The first step in numerical solution techniques is to discretize the solution vector or the physical properties in the RTE. The quantities requiring discretization are the spatial coordinates, the directional coordinates, the wavelengths, and/or the dust properties. Let us divide up ’z’ into cells with indices i $=$ 1,2,3,... ,Nz where Nz is the chosen number of grid cells. The cell walls, which separate the cells, also have indices, which we will give half-numbers: i $=$ 1/2, 3/2, 5/2,... ,Nz + 1/2 (See Figure).

The next step in the technique is to choose a numerical integration technique. We will discuss the first and second order versions of the method of Olson and Kunasz (1987, J. Quant. Spectros. Radiat. Transfer 38, 325) - OK87 method - in this example, since it is generally stable and reliable. However, depending on the complexity of the problem, other more suitable techniques may be required to get correct results.

Next, we need to make an assumption for the functional form of j(z) and $\alpha$ (z) between the cell boundaries, depending on the expected physical processes that we expect to dominate the radiation. For instance, one might use the emission and absorption co-efficient for small dust grains high up in the atmosphere. With these chosen functionals, we can then solve the formal transfer equation exactly. We illustrate the first and second order methods here with the higher order methods following the same general approach.

First Order

The first order version of OK87 method assumes that the emissivity j and extinction coefficient $\alpha$ are constant within each cell, but can be arbitrarily different from one cell to the next. Thus the Optical Depth in the cell ’i’ would be

$\delta \tau _{i}=(s_{i+1/2}-s_{i-1/2})\alpha _{i}\,\!$ or

$\delta \tau _{i}=(z_{i+1/2}-z_{i-1/2})\alpha _{i}/\mu \,\!$ Thus the source function S$_{i}={\frac {j_{i}}{\alpha _{i}}}$ is constant and the exact solution for RTE is

$I_{i+1/2}=\exp ^{-\delta \tau _{i}}I_{i-1/2}+S_{i}(1-\exp ^{-\delta \tau _{i}})\,\!$ Thus, after choosing a direction $\theta$ (i.e. $\mu$ ) to calculate I$_{\nu }$ (z,$\mu$ ) measured by an observer, we can integrate systematically from one grid wall to the next. We start at i $=$ 1 and step up for $\mu >0$ whereas we initialize at the top and work our way down for $\mu <0$ .

Second Order

In the second order method, we assume that S and $\alpha$ vary linearly between the cell walls (instead of constant within the cell). For cell i, this implies (for $\mu >0$ )

{\begin{aligned}S(z)={\frac {z_{i+1/2}-z}{z_{i+1/2}-z_{i-1/2}}}S_{i-1/2}+{\frac {z-z_{i-1/2}}{z_{i+1/2}-z_{i-1/2}}}S_{i+1/2}\\\alpha (z)={\frac {z_{i+1/2}-z}{z_{i+1/2}-z_{i-1/2}}}\alpha _{i-1/2}+{\frac {z-z_{i-1/2}}{z_{i+1/2}-z_{i-1/2}}}\alpha _{i+1/2}\\I_{i+1/2}=e^{-\delta \tau _{i}}I_{i-1/2}+Q_{i}\\Q_{i}={\frac {1-(1+\delta \tau _{i})e^{-\delta \tau _{i}}}{\delta \tau _{i}}}S_{i-1/2}+{\frac {\delta \tau _{i}-1+\exp ^{-\delta \tau _{i}}}{\delta \tau _{i}}}S_{i+1/2}\end{aligned}}\,\! Boundary/Initial Conditions

When $\mu >0$ , we have to integrate from z $=$ 0 upward and choose an appropriate boundary condition for I$_{1/2}$ in the first cell. For instance, in the plane atmosphere case we have been discussing here, the Plank black-body function B$_{\nu }(T)$ is a reasonable choice for mid-infrared wavelengths whereas reflected light of the solar spectrum is a more appropriate choice for the optical wavelength regime. Similarly, for $\mu <0$ , we have to impose a boundary condition for I$_{Nz+1/2}$ . Generally, I$_{Nz+1/2}=0$ is a reasonable choice for astrophysical scenarios. However, due to irradiation by the Sun, I$_{Nz+1/2}\neq 0$ and it depends on both angles $\theta$ and $\phi$ . Since this breaks the rotational symmetry in x-y plane, we have to include another variable in our problem.

Note that for higher order schemes, we use more values of S(z). For example, in the third order integration scheme for obtaining I_i+1/2 for $\mu >0$ , we do not only use S_i-1/2 and S_i+1/2, but also S_i+3/2. The subgrid model for S(z) is thus a quadratic fit through three values.

Appropriate Grid Size

Some general points to consider for the right spatial resolution which both provides accurate results and is numerically efficient -

• Try to spatially resolve the photosphere of the object with sufficient number of grid points, because it is here where the observed spectrum is formed.
• Use a stable numerical integration scheme that also works properly when large steps in $\tau$ are taken.
• For high optical depths use second order/higher integration methods if possible.
• Regions that are at high optical depth at all wavelengths (i.e at peak of Plank function at that wavelength e.g. Synchrotron self absorption) can be mapped with optically thick grid spacing.
• Posteriori checks - the intensity function I$_{\nu }$ along a light ray should not make large jumps from one grid cell) to the next
• Always choosing $\delta \tau <1$ over each grid cell is not required or advised.

Higher Dimensional Problems

For radiative transfer problems with 2-D/3-D spatial dimensions, the basic principles of the numerical integration method can be directly generalized. One important difference between the 1D vs higher dimensional problems is that each light ray passing through the grid would will intersect with the cell walls multiple times and get divided into a number of ray segments. Each of these segment belongs to a cell that is being traversed and the segment length can be very irregular. Thus, proper consideration for the changing $\delta s_{i}$ (the ray segment length in cell i) is required while calculating the intensity integration along a light ray.

The difficult aspect of radiative transfer calculations is that we generally do not know the values of j$_{\nu }$ and $\alpha _{\nu }$ in advance since the radiation field can affect the medium in such a way as to modify j$_{\nu }$ and $\alpha _{\nu }$ . Additionally, we cannot separate the problem into individual rays (as we did in the simple example) because a change in extinction/emission coefficients due to a particular ray at a particular location (x) will affect the formal transfer equation for all rays passing through x. Besides this "ray coupling" challenge, we also have to consider the "radiative cell coupling" - emission from a volume element can travel to distant regions - which makes the problem non-local.

There are primarily three main classes of numerical techniques which can be used to solve such problems - (1) Monte Carlo methods, which simulate the multiple scattering/emission/absorption processes directly, (2)Discrete ordinate methods, which solve the problem by dividing all coordinates, including the angles and the frequency, into discrete grid points or grid cells, and (3) Moment methods, including the diffusion method, which treat the angular and/or frequency domain in terms of its moments.

Monte Carlo Method

The primary concept behind the Monte Carlo method is to follow the path of a photon from one scattering event to the next, and to use random numbers to decide in which direction the photon will proceed after each scattering event. We must also use a random number to find out where along the present path the next scattering event will take place. By repeating this process for many millions of such photons, we can get a good statistical sample of paths that the photons will follow. Since the optical depths $\tau$ that the photon will travel until the next scattering event is given by the probability distribution –

$p(\tau )=\exp ^{-\tau }\,\!$ , a uniformly distributed random number (w) is drawn between 0 and 1, and

$\tau =-ln(w)\,\!$ . Once we have fixed the optical depth where the next scattering event would take place, the light ray travels cell-by-cell until it arrives at the point where the combined optical depth between that point and the location of previous scattering event is equal to the chosen $\tau$ . Since the light ray/photon can also get absorbed instead of scattering, we have to calculate the probability that the interaction in a particular cell is a scattering event; it is equal to the dust albedo a=$\kappa _{sca}/\kappa _{ext}$ . Using a random number w, the nature of the interaction is easily determined: If w $<$ a, we have a scattering event, otherwise an absorption event. Once the photon is absorbed or it exits from the considered medium, the simulation moves to the next photon.

It is important to note that when a "photon" is modeled in the Monte Carlo method, it typically represents many photons at once. Thus, we are really considering the trajectories of photon packets, each with a large number of photons. The approximation we make is to assume that all the photons in a single photon packet follow the same path. Thus, instead of considering absolute absorption, a fraction of the photon packet luminosity is absorbed and it is stored in the interaction cell. This absorbed luminosity is used at a later stage to calculate the emission spectrum from the medium in that volume cell thus acting as the source function for another MC cycle.Once, we have a statistical representation of the overall radiation field, we can then use it to predict what an observer would detect based on his orientation and location. Typically, for an observer at infinity one would assign a solid angle around the direction of the observer, and collect all photons that end up in that direction. Since an observer has to collect a large number of photon packets to reduce the noise, Monte Carlo method is very numerically intensive.

A related method, which is often used to circumvent this issue is to continuously "peel off" radiation from a photon packet as it passes through the medium, and compute what the observer would see from that photon packet, including the extinction from each of these points to the observer. Another approach is called the scattering source function approach. In this method, we follow each photon packet as it moves through the cloud. As it does so, we add an appropriate contribution of this photon packet to the scattering emissivity function $j_{\nu }$ of each cell it passes. The length of the ray segment is taken into account, so that if a photon packet passes right through the middle of a cell, its contribution to that cell will be larger than if it just passes through the edge. Once we have launched all photon packets, the function $j_{\nu }(x)$ is known throughout the medium. We can thus get the final image by integrating the formal transfer equation along a pre-defined set of rays all ending at the observer (like in previous section). Additional variations of MC radiative transfer techniques have been developed to optimize numerical efficiency - See Steinacker et al. 2013 article (link at the top) for a more detailed discussion.

The Following figures shows the typical procedure of a radiative transfer calculation as well an example from a paper modelling a Proto-planetary Disk (Mathews et al. 2013)  