Reaction Rate Calculations

The CALCULATE_REACTION_RATES subroutine evaluates the rate coefficients for all reactions in the chemical network at the local gas temperature, dust temperature, radiation field, visual extinction, and column densities. The resulting rates are used by the chemical solver to advance the abundances.

Reaction rates depend on local physical conditions: gas temperature, dust temperature, radiation field, extinction, column densities, electron density, and cosmic-ray ionization rate.
Multiple entries for the same reaction are allowed in the rate file (DUPLICATE flag) and are activated based on their valid temperature ranges.
Photodissociation of H₂ and CO, and photoionization of C and S, are treated explicitly with shielding and ray tracing, rather than via simple analytic fits.

Inputs and Outputs

Key inputs

TEMPERATURE: gas temperature (K)
DUST_TEMPERATURE: dust temperature (K)
RAD_SURFACE(J): incident radiation field along ray J
AV(J): visual extinction along ray J
COLUMN_NH2, COLUMN_NHD, COLUMN_NCO, COLUMN_NC, COLUMN_NS: column densities along rays
ALPHA, BETA, GAMMA: Arrhenius-type reaction parameters
RTMIN, RTMAX: temperature validity range
ZETALOCAL: local cosmic-ray ionization rate
nelectron, density: electron and gas densities

Outputs

RATE(I): rate coefficient for reaction I
Stored indices for key reactions: NRGR, NRH2, NRHD, NRCO, NRCI, NRSI

Thermal Gas-Phase Reactions

Most two-body gas-phase reactions are computed using a modified Arrhenius form:

\[k(T) = \alpha \left(\frac{T}{300}\right)^\beta \exp\!\left(-\frac{\gamma}{T}\right)\]

Important implementation details:

Reactions with large negative activation energies (\(\gamma < -200\)) are suppressed below their minimum valid temperature.
For duplicated reactions, the code selects the appropriate entry based on the current temperature.
Rates are capped at unity (except for grain reactions) to maintain numerical stability.

H₂ Formation on Grains

The formation of molecular hydrogen on dust grains is treated separately and does not use the standard Arrhenius form.

Depending on compile-time options, the rate follows:

Cazaux & Tielens (2002); see detailed description here,
a simplified temperature-dependent prescription given by \(3\times10^{-18}\sqrt{T_{\rm gas}}e^{-\frac{T_{\rm gas}}{10^3 {\rm K}}}\), or
the rate of Röllig et al. (2007) given by \(3\times10^{-18}\sqrt{T_{\rm gas}}\).

The selected rate depends explicitly on both gas and dust temperatures.

The corresponding reaction index is stored as:

NRGR — H₂ grain formation

Suprathermal Ion–Neutral Reactions

Optionally, suprathermal chemistry is included following Visser et al. (2009).

For ion–neutral reactions at low visual extinction, the effective temperature is increased by a contribution proportional to the Alfvén speed:

\[T_{\rm eff} = T + \Delta T_{\rm sup}\]

This enhancement applies only when the local extinction is below a critical value Av_crit.

Photoreactions

Photoreaction rates are computed by explicit ray integration:

\[k = \sum_J \alpha \, G_J \, e^{-\gamma A_{V,J}}\]

where the sum runs over all rays J.

Special Cases

The following reactions are treated with dedicated shielding functions:

H2 photodissociation Computed using H2PDRATE and self-shielding by H₂.
HD photodissociation Treated analogously to H₂.
CO photodissociation Includes shielding by CO and H₂ via COPDRATE.
C and S photoionization Computed with CIPDRATE and SIPDRATE, including temperature dependence and column-density shielding.

The indices of these reactions are stored for later use (NRH2, NRHD, NRCO, NRCI, NRSI).

Cosmic-Ray Ionization

Primary cosmic-ray ionization rates are proportional to the local ionization rate:

\[k = \alpha \, \zeta_{\rm local}\]

Duplicate reactions are again selected by temperature range.

Cosmic-Ray–Induced Photoreactions

Secondary UV photons generated by cosmic rays drive additional photoreactions. Their rates follow:

\[k = \alpha \, \zeta_{\rm local} \left(\frac{T}{300}\right)^\beta \frac{\gamma}{1 - \omega}\]

where \(\omega\) is the dust albedo.

Freeze-Out onto Dust Grains

Neutral species and singly charged ions can freeze onto dust grains.

The rate depends on:

the thermal velocity of the species,
Coulomb focusing (for ions),
a fixed sticking probability (0.3).

The general scaling is:

\[k \propto \sqrt{T} \, C_{\rm ion} \, S\]

Desorption Processes

Cosmic-Ray Desorption

Desorption due to transient grain heating by cosmic rays follows Roberts et al. (2007), with a fixed cosmic-ray flux and a temperature- dependent yield.

Photodesorption

Photodesorption rates depend on:

dust temperature (via the yield),
attenuated FUV flux along each ray.

Thermal Desorption

Thermal evaporation from grains follows Hasegawa, Herbst & Leung (1992) and depends exponentially on the dust temperature:

\[k \propto \exp\!\left(-\frac{E_b}{T_{\rm dust}}\right)\]

Grain-Surface Reactions

Grain-mantle reactions are treated with constant rates provided directly in the reaction file:

\[k = \alpha\]

Grain-Assisted Recombination

Optionally, grain-assisted recombination of H⁺, He⁺, and C⁺ is included following Gong et al. (2017).

These rates depend on:

electron density,
gas density,
local FUV radiation field,
dust charging parameter \(\Psi\).

Numerical Safeguards

To ensure numerical stability:

Negative rates trigger a fatal error.
Rates below \(10^{-99}\) are set to zero.
Gas-phase rates are capped at unity.
Grain-surface and desorption reactions are allowed to exceed unity.

Summary

The reaction rate module in 3D-PDR:

Supports a wide range of gas-phase, grain-surface, photo-, and cosmic-ray–driven reactions;
Includes detailed ray-based attenuation and self-shielding;
Allows multiple temperature-dependent rate entries per reaction;
Incorporates optional suprathermal and grain-assisted processes;
Ensures physically consistent and numerically stable rate coefficients across extreme PDR conditions.