Theory

Basic radiative transfer

We briefly discuss basic theory of radiative transfer that is relevant for pyhonradex. A more detailed discussion can for example be found in [Rybicki04] or [Rybicki91].

The radiation field in every point of space can be described by the specific intensity \(I_{\nu}\), defined as the energy radiated per unit of time, surface, frequency and solid angle, i.e., \(I_{\nu}\) has units of W/m²/sr/Hz. The differential equation describing the change of the specific intensity along a spatial coordinate \(s\) is given by

\[\frac{\mathrm{d}I_\nu}{\mathrm{d}s} = -\alpha I_\nu + j_\nu\]

Here, \(\alpha\) is the absorption coefficient in m^-1. It describes how much is removed from the beam per unit length. On the other hand, the emission coefficient \(j_\nu\) is the energy emitted per unit time, solid angle, volume and frequency. Defining the optical depth as \(\mathrm{d}\tau=\alpha\mathrm{d}s\), we can rewrite the equation as

\[\frac{\mathrm{d}I_\nu}{\mathrm{d}\tau} = -I_\nu + S_\nu\]

with the source function \(S_\nu=\frac{j_\nu}{\alpha}\). Our basic goal is to solve this equation. For example, for a uniform medium (the emission and absorption coefficients do not depend on position) the solution reads \(I_\nu=I_\nu(0)e^{-\tau}+S_\nu(1-e^{-\tau})\). In summary, we need to know the emission coefficient and the absorption coefficient to solve the radiative transfer.

Now considering a gas, the absorption coefficient for a transition between energy levels \(l\) and \(l'\) (where \(E_l>E_{l'}\), hereafter indicated by \(l\succ l'\)) of a molecule is given by [Rybicki04]

\[\alpha_{ll'} = \frac{h\nu}{4\pi}(n_{l'}B_{l'l}-n_lB_{ll'})\phi_\nu^{ll'}\]

Here \(h\) is the Planck constant, \(n_l\) the number density of molecules in energy level \(l\), \(B_{l'l}\) and \(B_{ll'}\) are the Einstein coefficients for absorption and stimulated emission respectively, and \(\phi_\nu\) is the normalised line profile (i.e. \(\int\phi_\nu\mathrm{d}\nu=1\) and \(\phi_\nu\) describes how the energy is distributed over frequency).

The emission coefficient is given by (where again \(l\succ l'\))

\[j_{ll'} = \frac{h\nu}{4\pi}n_lA_{ll'}\phi_\nu^{ll'}\]

with \(A_{ll'}\) the Einstein coefficient for spontaneous emission.

The total absorption and emission coefficients are simply given by the sum over all transitions, plus contributions from the dust continuum (\(\alpha_c\) and \(j_c\)):

\[ \begin{align}\begin{aligned}\alpha = \sum_{l\succ l'}\alpha_{ll'}(\nu) + \alpha_c(\nu)\\j_\nu = \sum_{l\succ l'}j_{ll'}(\nu) + j_c(\nu)\end{aligned}\end{align} \]

We see that we need to know the fractional level population \(n_l\) to calculate \(\alpha\) and \(j_\nu\), which are needed to solve the radiative transfer. How can we do that? There are two kinds processes that can excite or de-excite a molecular level: radiative processes (emission or absorption of photons) or collisions. If the density of colliders (e.g. H₂ or electrons) is high enough, the energy levels become thermalised and we are in local thermodynamic equilibrium (LTE). In this case, the level population is simply given by the Boltzmann distribution:

\[n_l = n\frac{e^{-E_l/(kT_\mathrm{kin})}}{Q}\]

where \(Q=\sum_{l'} e^{-E_{l'}/(kT_\mathrm{kin})}\) is the partition function, \(T_\mathrm{kin}\) is the kinetic temperature of the gas, and \(n\) is the number density of the molecule. Thus, in the LTE case, the radiative transfer can be solved easily. In the more general non-LTE case, the level population is not known a priori and needs to be calculated. In that case, the level population is often characterised by an excitation temperature \(T_\mathrm{ex}\) defined by

\[\frac{n_l}{n_{l'}} = \frac{g_l}{g_{l'}}e^{-\Delta E/(kT_\mathrm{ex})}\]

where \(g\) is the statistical weight and \(\Delta E\) the energy difference between levels \(l\) and \(l'\) (note that in LTE, \(T_\mathrm{kin}=T_\mathrm{ex}\)).

Now, how can we calculate the level populations if LTE does not apply? We assume statistical equilibrium (SE). In other words, we assume that the rate of processes that populates a level equals the rate of processes that depopulates a level:

\[\sum_{l'}n_{l'}(C_{l'l}+R_{l'l}) = \sum_{l'}n_l(C_{ll'}+R_{ll'})\]

Here \(C_{l'l}\) and \(R_{l'l}\) are the rates per volume of collisional and radiative transitions respectively, from level \(l'\) to level \(l\). Thus, the left-hand side is the total rate of transitions into level \(l\), while the right-hand side is the total rate of transitions out of level \(l\). By writing down the statistical equilibrium for each level and solving the system of equations, the level populations can be calculated and the radiative transfer solved.

However, there is a problem: while the collisional rates \(C_{l'l}\) are known, the radiative rates depend, as one might expect, on the radiation field. If \(l\succ l'\), we have

\[R_{ll'} = A_{ll'} + B_{ll'}\bar{J}\]

while if \(l\prec l'\)

\[R_{ll'} = B_{ll'}\bar{J}\]

Here, \(\bar{J}\) is the radiation field averaged over frequency and solid angle: \(\bar{J}=\frac{1}{4\pi}\int I_\nu\phi_\nu^{ll'}\mathrm{d}\Omega\mathrm{d}\nu\). Thus, in order to solve the SE equations, we need to know the radiation field, which is what we were after in the first place… In order to solve this chicken and egg problem, we need to adopt an iterative technique.

Accelerated Lambda Iteration (ALI)

For the following discussion, we introduce the Lamda Operator, which essentially is a way of writing down the formal solution of the radiative transfer:

\[I_\nu = \Lambda[S_\nu]\]

So, the Lambda operator computes the radiation field \(I_\nu\) for a given source function \(S_\nu=\frac{j_\nu}{\alpha}\), the latter being completely determined by the level population (and the known contribution from the dust continuum). The simplest iteration scheme (the so-called Lambda Iteration scheme) is very straightforward: one starts with an initial guess for the level population and solves the radiative transfer (as formalised by the above equation). The solution \(I_\nu\) is then inserted into the statistical equilibrium equations, resulting in an updated level population. This procedure is repeated until convergence is established.

However, the Lambda iteration scheme can suffer from extremely slow convergence in optically thick systems (see e.g. the lecture notes by Dullemond or [Rybicki91]). An alternative scheme, known as Accelerated Lambda Iterations (ALI), provides much better convergence. The idea is to introduce an approximate Lambda operator \(\Lambda^*\) and to write

\[I_\nu = \Lambda^*[S_\nu] + (\Lambda-\Lambda^*)[S_\nu^\dagger]\]

where the \(\dagger\) indicates quantities from the previous iteration. This is inserted into the equations of statistical equilibrium, which can then be solved for an updated level population (and thus updated \(S_\nu\)). See for example [Rybicki91] or [Hubeny03] for more details about ALI.

The ALI method by Rybicki & Hummer (1992)

pythonradex implements a variation of the ALI scheme presented by [Rybicki92]. The method is capable of treating overlapping lines and full continuum. In contrast, RADEX can only treat non-overlapping lines without continuum.

pythonradex implements the “Full preconditioning strategy” presented in section 2.3 of [Rybicki92]. Instead of a Lambda operator, the method presented in [Rybicki92] uses a Psi operator that acts on the emission coefficient:

\[I_\nu = \Psi[j_\nu]\]

The approximate iteration scheme is then based on \(I_\nu=\Psi^*[j_\nu] + (\Psi-\Psi^*)[j_\nu^\dagger]\), which is inserted into the equations of statistical equilibrium.

Escape probability

We still need to specify the formal solution of the radiative transfer we adopt via the operator \(\Psi\). Same as RADEX, we use an escape probability method. We consider the probability \(\beta\) of a newly created photon to escape the source. This probability depends on the geometry of the source and the absorption coefficient (or, equivalently, optical depth). If the source is completely optically thick (\(\beta\approx 0\)), we expect the radiation field to equal the source function \(S_\nu=\frac{j_\nu}{\alpha}\). Thus, we write

\[I_\nu = \Psi[j_\nu] = \beta(\alpha^\dagger) I_\mathrm{ext} + (1-\beta(\alpha^\dagger))\frac{j_\nu}{\alpha^\dagger}\]

Here \(I_\mathrm{ext}\) is an external radiation field that irradiates the source from the outside (for example the CMB). If the source is completely optically thick, external radiation cannot penetrate the source and the corresponding term vanishes.

For the approximate Psi operator, we choose

\[\Psi^*[j_\nu] = (1-\beta(\alpha^\dagger))\frac{j_\nu}{\alpha^\dagger}\]

Please see the section about source geometries for a list of all geometries available in pythonradex with the corresponding formulas for the escape probability.

Ng-acceleration

pyhonradex employs Ng-acceleration [Ng74] to further accelerate convergence. Ng-acceleration uses the last three iteration steps to compute the next step. See the lecture notes by Dullemond (Sect. 4.4.7) for more details.