Source geometries

The calculations of the escape probability and the emerging flux depend on the adopted geometry. Below we give an overview of the geometries available in pythonradex. For all geometries, we assume a homogeneous medium (i.e. constant number density). For convenience, the escape probabilities are expressed as functions of optical depth \(\tau\) rather than the absorption coefficient \(\alpha\). The specific intensities are expressed using the source function \(S_\nu\). In absence of dust and overlapping lines, it is simply given by the Planck function evaluated at the excitation temperature of the transition: \(S_\nu=B_\nu(T_\mathrm{ex})\).

Note that RADEX always uses the formula for a slab geometry to calculate the specific intensity, regardless of the adopted geometry. This is incorrect for spherical geometries. See Differences between pythonradex and RADEX for more details.

Static sphere

A homogeneous, static sphere. The escape probability is given in [Osterbrock74]:

(1)\[\beta(\tau) = \frac{3}{2\tau}\left(1-\frac{2}{\tau^2}+\left(\frac{2}{\tau}+\frac{2}{\tau^2}\right) e^{-\tau}\right)\]

where \(\tau\) is the optical depth of the diameter of the sphere. A mean specific intensity (in [W/m²/Hz/sr]) can also be calculated from [Osterbrock74] and expressed as:

(2)\[I_\nu = \frac{2S_\nu}{\tau^2}\left(\frac{\tau^2}{2}-1+(\tau+1)e^{-\tau}\right)\]

where \(S_\nu\) is the source function. Due to the spherical geometry, the specific intensity is actually not constant across the source, unless it is completely optically thick (in general, the specific intensity towards the center is higher compared to the regions close to the sphere edge). The specific intensity given above corresponds to a mean, in the sense that if multiplied by the solid angle of the sphere (given by \(\Omega=R^2\pi/d^2\) with \(R\) the radius of the sphere and \(d\) the distance of the source), it gives the correct flux density (in [W/m²/Hz]). This mean intensity (or the corresponding brightness temperature) is what pythonradex returns to the user.

Derivation of escape probability and flux density by Osterbrock (1974)

Since the book by [Osterbrock74] is not easily accessible online, we present here the derivation of the escape probability and flux density for a static sphere. Consider the following sketch showing a sphere with optical depth \(\tau\) along its diameter:

Static sphere with optical depth \(\tau\) along the diameter. We consider a ray that makes an angle \(\theta\) with respect to the outward normal. The optical depth along that ray is \(\tau\cos\theta\).

Consider a ray making an angle \(\theta\) with respect to the outward normal. The optical depth along that ray is \(\tau\cos\theta\). Using the fact that the sphere is uniform (i.e. the emission coefficient \(j_\nu\) is constant), the specific intensity in the direction of the ray, at the surface of the sphere, is given by (see [Rybicki04], equation 1.30)

\[I_\nu(\theta) = \frac{j_\nu}{\alpha}(1-e^{-\tau\cos\theta})\]

where \(\alpha\) is the absorption coefficient. Next, we compute the outward flux density at the surface of the sphere (energy per unit area, time and frequency) by integrating over all outward directions. Here we use the solid angle element \(\mathrm{d}\Omega=\sin\theta\mathrm{d}\theta\mathrm{d}\phi\):

\[F_\nu^s = \int_0^{2\pi}\int_0^{\pi/2}I_\nu(\theta)\cos\theta\sin\theta\mathrm{d}\theta\mathrm{d}\phi = \frac{2\pi j_\nu}{\alpha\tau^2}\left( \frac{\tau^2}{2}-1+(\tau+1)e^{-\tau} \right)\]

In this expression, the cosine is necessary to take into account the projection of the surface of the sphere with respect to the outward direction. To calculate the flux density observed at the telescope, we multiply \(F_\nu^s\) with the surface of the sphere (\(4\pi R^2\)) and divide by \(4\pi d^2\) where \(d\) is the distance to the source:

\[F_\nu = \frac{2\pi j_\nu}{\alpha\tau^2}\left( \frac{\tau^2}{2}-1+(\tau+1)e^{-\tau} \right)\frac{4\pi R^2}{4\pi d^2}\]

Finally, \(F_\nu\) is divided by the solid angle of the sphere \(\Omega=\frac{R^2\pi}{d^2}\) to derive the observed specific intensity. Identifying the source function \(S_\nu=\frac{j_\nu}{\alpha}\), we recover Eq. 2. To calculate the escape probability, we consider the flux density at the surface of the sphere in the optically thin limit where all photons escape. It is simply given by the total emission within the spherical volume, divided by the surface of the sphere:

\[F_\nu^{s,\mathrm{thin}} = \frac{\frac{4}{3}R^3\pi4\pi j_\nu}{4\pi R^2} = \frac{4\pi R j_\nu}{3}\]

where \(4\pi j_\nu\) is the emission coefficient integrated over all directions (solid angles). Thus, the escape probability is given by

\[\beta(\tau) = \frac{F_\nu^s}{F_\nu^{s,\mathrm{thin}}} = \frac{3}{2R\alpha\tau^2}\left( \frac{\tau^2}{2}-1+(\tau+1)e^{-\tau} \right) = \frac{3}{\tau^3}\left( \frac{\tau^2}{2}-1+(\tau+1)e^{-\tau} \right)\]

Here we used the definition of the optical depth (\(\tau=2R\alpha\)). One easily verifies that the above expression is the same as Eq. 1.

Static slab

A homogeneous, static slab. The escape probability is given by (e.g. [Elitzur92]):

\[\beta(\tau) = \frac{\int_0^1 (1-e^{-\tau/\mu})\mu\mathrm{d}\mu}{\tau}\]

The specific intensity (in [W/m²/Hz/sr]) is given by:

(3)\[I_\nu = S_\nu(1-e^{-\tau})\]

To get the flux density in [W/m²/Hz], \(I_\nu\) needs to be multiplied by the solid angle of the emitting region.

Derivation of escape probability for a homogeneous static slab

The escape probability for a static slab can be derived following a method similar to the static sphere. Consider a slab in the \(x\)-\(y\) plane with a thickness in the \(z\)-direction of \(L\). It is assumed that the slab’s extension in the \(x\) and \(y\) directions is much larger compared to its thickness. Consider a ray making an angle \(\theta\) with the surface normal of the slab.

Static slab with depth \(L\) and optical depth \(\tau\). We consider a ray that makes an angle \(\theta\) with respect to the surface normal. The optical depth along that ray is \(\tau/\cos\theta\).

The specific intensity is given by

\[I_\nu(\theta) = \frac{j_\nu}{\alpha}(1-e^{-\tau/\cos\theta})\]

where \(\tau\) is the optical depth of the slab along the surface normal. As in the case of the static sphere, we derive the flux density at the surface of the slab:

\[F_\nu^s = 2\pi\int_0^{\pi/2}I_\nu(\theta)\cos\theta\sin\theta\mathrm{d}\theta = 2\pi\frac{j_\nu}{\alpha}\int_0^{\pi/2}(1-e^{-\tau/\cos\theta})\cos\theta\sin\theta\mathrm{d}\theta\]

Setting \(\mu=\cos\theta\), this transforms to

\[F_\nu^s = 2\pi\frac{j_\nu}{\alpha} \int_0^1(1-e^{-\tau/\mu})\mu\mathrm{d}\mu\]

Now consider the flux density in absence of absorption. Like for the sphere, it is given by the total emission within the volume of the slab, divided by its surface. Assume that one side of the slab (front or back side) has a surface area \(A\) (thus the volume is \(A L\)). Then

\[F_\nu^{s,\mathrm{thin}} = \frac{A L 4\pi j_\nu}{2A} = 2\pi j_\nu L\]

where \(2A\) is the total surface area of the slab (front and back). Like for the static sphere, the escape probability is equal to the ratio \(F_\nu^s/F_\nu^{s,\mathrm{thin}}\):

\[\beta(\tau) = \frac{F_\nu^s}{F_\nu^{s,\mathrm{thin}}} = \frac{\int_0^1(1-e^{-\tau/\mu})\mu\mathrm{d}\mu}{\alpha L} = \frac{\int_0^1(1-e^{-\tau/\mu})\mu\mathrm{d}\mu}{\tau}\]

where the definition of optical depth was used for the last step.

LVG sphere

The Large Velocity Gradient (LVG) approximation is applicable if the characteristic flow velocity along the line of sight is much larger than the local random (e.g. thermal) velocities (e.g. [Scoville74], [Elitzur92]). In other words, all photons escape due to Doppler shifting, unless absorbed locally, that is, close to the emission location. The local region where an emitted photon can be absorbed is sometimes referred to as the radiatively connected region.

We use the model by [Goldreich74]: A homogeneous sphere with a constant, radial velocity gradient \(\mathrm{d}v/\mathrm{d}r=V/R\) where \(V\) is the velocity at the sphere surface and \(R\) is the radius of the sphere. The escape probability is given by (see e.g. Eq. 5 in [Goldreich74])

(4)\[\beta(\tau) = \frac{1-e^{-\tau}}{\tau}\]

\(\tau\) is the optical depth of the radiatively connected region. When looking through the centre of the sphere, the optical depth equals \(\tau\) in the interval \([-V,V]\), and is zero outside. The observed specific intensity can be derived by using an approach similar to [deJong75]. Consider the following sketch:

Sketch of an LVG sphere. We consider a radiatively connected region. It is a thin slab perpendicular to the line of sight, at a distance \(p\) from the center of the sphere. The velocity projected onto the line of sight is constant in that slab.

The radiatively connected region along the line of sight is a thin slab perpendicular to the line of sight. Indeed, consider such a slab at a distance \(p\) from the center of the sphere. The velocity projected onto the line of sight is indeed constant in the slab:

\[v_p = v(x)\cos\mu = v(x)\frac{p}{x} = \frac{V}{R}x \frac{p}{x} = \frac{V}{R}p\]

which is independent of \(x\). Here we used the assumption that the radial velocity gradient is constant and given by \(V/R\). Note that \(v_p\) is proportional to \(p\) (and \(p<0\) if the slab is located in the half sphere closer to the observer), so each slab will have a different projected velocity. We again see that any photon escapes, unless absorbed locally (within the radiatively connected region).

To calculate the flux density at velocity \(v_p\), we can simply integrate over the solid angle of the slab as seen by the observer. At each point of the slab, the specific intensity is given by \(S_\nu(1-e^{-\tau})\), where it is understood that the frequency \(\nu\) corresponds to the velocity \(v_p\), and \(\tau\) is the optical depth of the slab. The element of solid angle is given by \(\mathrm{d}\Omega_p=\frac{r\mathrm{d}\phi\mathrm{d}r}{d^2}\) where \(r\) and \(\phi\) are polar coordinates on the slab surface and \(d\) is the distance to the observer. Thus, the flux density is given by

\[F_\nu = \int_0^{2\pi}\int_0^{\sqrt{R^2-p^2}}S_\nu(1-e^{-\tau})\frac{r}{d^2}\mathrm{d}r\mathrm{d}\phi = \frac{2\pi S_\nu}{d^2}(1-e^{-\tau})\int_0^{\sqrt{R^2-p^2}}r\mathrm{d}r\]

Using \(p=\frac{v_p R}{V}\) and the solid angle of the sphere \(\Omega=\frac{R^2\pi}{d^2}\), this evaluates to

\[F_\nu = S_\nu(1-e^{-\tau})\Omega\left(1-\frac{v^2}{V^2}\right)\]

where we dropped the index \(p\) from the velocity \(v\). The latter can be converted to frequency as \(v=c(1-\nu/\nu_0)\) with \(c\) the speed of light and \(\nu_0\) the rest frequency. Note that the above equation applies for \(|v|\leq V\). For \(|v|>V\), the flux vanishes. Finally, the specific intensity is simply expressed as

\[I_\nu = F_\nu/\Omega = S_\nu(1-e^{-\tau})\left(1-\frac{v^2}{V^2}\right)\]

As in the case of the static sphere, also for the LVG sphere RADEX incorrectly uses the formula for a slab geometry to calculate the specific intensity (see Differences between pythonradex and RADEX for more details), although note that at \(v=0\) (line center), the two formulas agree. So in terms of emission at the line center, the two codes agree.

LVG slab

Here we consider a homogeneous slab with a constant velocity gradient \(\mathrm{d}v/\mathrm{d}z\), where \(z\) is the coordinate perpendicular to the slab surface, along the line of sight. The escape probability is given by (see Section II of [Scoville74] or Section III of [deJong75])

\[\beta(\tau) = \frac{1-e^{-3\tau}}{3\tau}\]

Here again, \(\tau\) is the optical depth of the radiatively connected region. It is constant over the velocity interval with a width given by \(\frac{\mathrm{d}v}{\mathrm{d}z}Z\) with \(Z\) the total depth of the slab, and zero outside. The specific intensity is given by the same formula as for the static slab.

Geometries emulating RADEX

Mainly for test and legacy purposes, pythonradex offers two additional geometries that emulate the behaviour of RADEX.

Static sphere RADEX

This geometry uses the same formula (Eq. 1) for the escape probability as the static sphere, but uses the formula (Eq. 3) for the static slab to calculate the specific intensity.

LVG sphere RADEX

The paper describing RADEX ([vanderTak07]) says that RADEX uses the formula given in Eq. 4 to calculate the escape probability of an LVG sphere. However, by inspecting the RADEX source code, it is found that at least the current version uses a different formula that seems to be based on [deJong80]. It is given by

\[\beta(\tau) = \frac{1}{\tau\sqrt{\ln(\tau/(2\sqrt{\pi}))}} \qquad \text{if } \tau\geq 7\]

and

\[\beta(\tau) = \frac{4-4e^{-2.34\tau/2}}{4.68\tau} \qquad \text{if } \tau< 7\]

Thus, the geometry “LVG sphere RADEX” uses this formula. For the flux density, it uses the same formula (Eq. 3) as for the static slab (despite the spherical geometry).

The line width parameter width_v

It is important to understand the different interpretations of the input parameter width_v used by pythonradex (see API of the Source class). For static geometries, this refers to the local emission width (typically the thermal width). pythonradex allows two kinds of local emission profiles (parameter line_profile_type): “Gaussian” (in which case width_v refers to the FWHM) and “rectangular”. On the other hand, for the LVG geometries (“LVG sphere” and “LVG slab”), width_v refers to the global velocity width of the source. For the “LVG sphere”, width_v is equal to \(2V\) (with \(V\) the velocity at the sphere surface). For the “LVG slab”, width_v equals \(\mathrm{d}v/\mathrm{d}zZ\) (with \(Z\) the depth of the slab and \(\mathrm{d}v/\mathrm{d}z\) the constant velocity gradient). For these geometries, the parameter line_profile_type needs to be set to “rectangular”. This ensures that the optical depth is calculated correctly.