Escape Probability and Line Cooling

This routine computes molecular line cooling and the mean radiation field using the Large Velocity Gradient (LVG) escape probability approximation along multiple rays (see description here). It is designed for use within the 3D-PDR framework and operates on a given grid cell, using level populations obtained from the statistical equilibrium solver.

The routine performs the following tasks:

Computes the optical depth \(\tau_{ij}\) for each radiative transition \(i \rightarrow j\) along all rays.
Computes the ray-by-ray escape probability \(\beta_{ij}^{(r)}\).
Averages the escape probability over rays to obtain \(\beta_{ij}\).
Evaluates the source function \(S_{ij}\).
Computes line cooling rates.
Updates the mean radiation field \(\langle J_{ij} \rangle\).
Constructs the radiative transition matrix used by the level population solver.

The LVG approximation assumes that velocity gradients Doppler-shift photons out of resonance, allowing radiative transfer to be treated locally via escape probabilities.

Physical Assumptions

Statistical equilibrium holds for all molecular levels.
Line transfer is treated locally using escape probabilities.
Velocity broadening includes both thermal and microturbulent components.
The radiation field includes contributions from:
- Cosmic Microwave Background (CMB)
- Thermal dust emission
Strong maser action is regulated to avoid numerical divergence.

Key Quantities

Velocity Dispersion

The effective inverse velocity dispersion entering the LVG optical depth is

\[\phi^{-1} = \frac{1}{\sqrt{ \frac{2 k_B T}{m} + v_{\mathrm{turb}}^2 }}\]

where:

\(T\) is the gas temperature
\(m\) is the molecular mass
\(v_{\mathrm{turb}}\) is the microturbulent velocity

This factor is stored in the code as frac2.

Einstein Prefactor

For each transition \(i \rightarrow j\), the LVG optical depth prefactor is

\[f_{ij} = \frac{A_{ij} c^3}{8 \pi \nu_{ij}^3}\]

implemented as frac1.

Optical Depth Calculation

The optical depth for a transition \(i \rightarrow j\) along ray \(r\) is computed as a sum over ray segments:

\[\tau_{ij}^{(r)} = \sum_k f_{ij} \, \phi^{-1} \left[ \frac{n_j g_i - n_i g_j}{g_j} \right]_k \Delta s_k\]

where:

\(n_i\), \(n_j\) are the level populations
\(g_i\), \(g_j\) are the statistical weights
\(\Delta s_k\) is the path length of segment \(k\)
The populations are evaluated at each ray integration point

In the code, this is accumulated as tau_ij(j).

One-dimensional models

When the ONEDIMENSIONAL flag is enabled:

Only one ray (ray index 6) contributes to radiative transfer.
All other rays are assigned extremely large optical depths:

\[\tau_{ij}^{(r \neq 6)} \rightarrow 10^{50}\]

This effectively enforces 1D slab geometry while preserving the ray-based infrastructure.

Escape Probability

The escape probability for each ray is computed using the standard LVG form:

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

Special Cases

To ensure numerical stability and physical realism, the following limits are applied:

Strong masing (\(\tau < -5\)):

\[\beta = \frac{1 - e^{5}}{-5}\]
Weak masing (\(-5 \le \tau < 0\)):

Standard LVG expression is used.
Optically thin limit (\(|\tau| < 10^{-8}\)):

\[\beta = 1\]

These limits prevent floating-point overflow and unphysical cooling rates.

Mean Escape Probability

The total escape probability for a transition is the average over all rays:

\[\begin{split}\beta_{ij} = \begin{cases} \sum_r \beta_{ij}^{(r)} & \text{(ONEDIMENSIONAL)} \\ \frac{1}{N_{\mathrm{rays}}} \sum_r \beta_{ij}^{(r)} & \text{(3D)} \end{cases}\end{split}\]

Source Function

The line source function is computed as

\[\begin{split}S_{ij} = \begin{cases} \dfrac{2 h \nu_{ij}^3}{c^2} \left[ \dfrac{1}{\dfrac{n_j g_i}{n_i g_j} - 1} \right] & n_i \neq 0 \\ 0 & n_i = 0 \end{cases}\end{split}\]

This follows the formulation used in UCL_PDR and is valid for arbitrary population inversions.

Background Radiation Field

The background radiation field consists of:

CMB Contribution

\[B_\nu(T_{\mathrm{cmb}}) = \frac{2 h \nu^3}{c^2} \left[ \frac{1}{e^{h\nu / k_B T_{\mathrm{cmb}}} - 1} \right]\]

Dust Emission

Thermal dust emission is included as

\[B_\nu(T_{\mathrm{dust}}) \, \epsilon_\nu\]

with emissivity

\[\epsilon_\nu = \rho_{\mathrm{grain}} \, n_{\mathrm{grain}} \left( 0.01 \frac{1.3 \nu}{3 \times 10^{11}} \right)\]

The total background field is

\[B_{ij} = B_\nu(T_{\mathrm{cmb}}) + B_\nu(T_{\mathrm{dust}})\,\epsilon_\nu\]

Mean Radiation Field

The angle-averaged radiation field entering the rate equations is

\[\langle J_{ij} \rangle = (1 - \beta_{ij}) S_{ij} + \beta_{ij} B_{ij}\]

This interpolates between optically thick trapping and optically thin escape.

Line Cooling Rate

The cooling rate contributed by transition \(i \rightarrow j\) is

\[\Lambda_{ij} = A_{ij} h \nu_{ij} \, n_i \, \beta_{ij} \frac{S_{ij} - B_{ij}}{S_{ij}}\]

The total cooling rate is the sum over all downward transitions:

\[\Lambda = \sum_{i>j} \Lambda_{ij}\]

Radiative Transition Matrix

The routine updates the radiative transition matrix:

\[R_{ij} = A_{ij} + B_{ij} \langle J_{ij} \rangle + C_{ij}\]

This matrix is used in the statistical equilibrium solver to update level populations.

Memory Management

Temporary arrays allocated within the routine:

tau_ij – optical depth per ray
field – mean radiation field matrix

These arrays are deallocated before returning.

Summary

This routine provides a self-consistent LVG-based treatment of line transfer, including dust and CMB backgrounds, ray-based anisotropy, and numerical safeguards for maser action. It is suitable for both 3D and 1D geometries and forms a core component of the 3D-PDR radiative cooling framework.