Photochemistry

JAFF provides first-class support for photodissociation and photoionisation reactions. This page explains the physical model, how JAFF reads cross-section data, and how rate coefficients are derived from a discretised radiation field.

Overview

A photochemical reaction has the form

The reaction rate per unit volume is

where \(\sigma(E_\nu)\) is the reaction cross section (\(\text{cm}^2\)), \(c\) is the speed of light (\(\text{cm s}^{-1}\)), and \(n_\gamma\) is the photon number density (\(\text{cm}^{-3}\)). The subscript \(\gamma\) signifies the local radiation field and a subscript of \(\nu\) signifies the frequency of the photon. JAFF uses high-resolution cross-section tables from the Leiden Observatory Database and the NORAD Database for photo-dissociaton and photo-ionization reactions respectively to store and calculate band averaged cross-sections.

When the radiation energy density per band is calculated, JAFF replaces the continuous integral with a sum over discrete radiation frequency bands:

where \(\rho_i\) is the radiation energy or photon density in band \(i\) and \(k_i\) is a pre-computed rate coefficient for that band.

Radiation Field Discretisation

The radiation field is divided into contiguous energy bands. You specify the band boundaries in jaffgen.toml (see Configuration File):

[network.radiation]
bands             = [13.6, "inf"]   # band edges in eV; "inf" for open upper bound
profile_index   = 0               # photon-number spectrum index α (or one per band, e.g. [0, 1])
mode    = "nph"           # photon number density ("nph") or energy density ("u")
rsl               = 2.99792458e10   # speed of light (cm/s)

Photon spectrum

The photon number density spectrum is assumed to be a power law in photon energy:

where \(\alpha\) is profile_index. Setting \(\alpha = 0\) gives \(n(E) \propto E^{-2}\), i.e. equal energy per logarithmic bin which is the default assumption.

profile_index may also be a list with one \(\alpha_i\) per band (e.g. profile_index = [0, 1] for two bands). Each band \(i\) then uses \(n(E) \propto E^{\alpha_i - 2}\) and is normalised on its own, so the spectrum is a histogram and is discontinuous at band edges by design.

Valid band / index combinations

RadiationProps checks the configuration when it is constructed:

Only the first band can reach \(E = 0\) and only the last can reach \(E = \infty\), so the divergence checks use those bands' indices; interior bands always converge.

Band-averaged cross section

For band \(i\) spanning \([E_\text{lo}, E_\text{hi}]\):

Rate coefficient

In photon-density mode (mode = "nph"):

In energy-density mode (mode = "u"):

where \(\langle E\rangle_i\) is the photon-number–weighted average energy in the band.

Generated Rate Expressions

After band integration, each photochemical reaction contributes one term per band to the ODE right-hand side:

The generated code uses the radiation density array den[i] as the runtime variable for \(\rho_i\). The ODE for the radiation field itself can be obtained from Network.sradodes().

Marking a Reaction as Photochemical

In the PRIZMO format, append PHOTO followed by the energy threshold (eV) instead of a rate expression:

H -> H+ + E    []    PHOTO, 13.60

JAFF will look up the matching cross section in its bundled database and integrate it over the configured radiation bands.

Cross-Section Data

JAFF bundles cross sections from three sources. At network-load time the serialized reaction key (Reactant1.Reactant2__Product1.Product2, the serialized form of the reaction; photo reactions carry the _PHOTON agent, e.g. H._PHOTON__H+.e-) is looked up in jaff.db, and the cross-section arrays are attached to the reaction's xsecs_dict.

Storage layout

Tabulated cross sections are stored as two collapsed HDF5 files, one group per serialized reaction, all datasets co-sorted by ascending photon energy. These files are not bundled in the package: on first network load JAFF downloads them (via pooch) from a remote mirror into the local data/xsecs directory, then reuses the cached copies on subsequent runs:

src/jaff/data/xsecs/leiden.hdf5      # one group per reaction (absorption + decay)
src/jaff/data/xsecs/norad.hdf5       # one group per reaction (ionization decay)
src/jaff/data/xsecs/verner_1996.csv  # analytic-fit parameters

photon_energy is stored in eV and every cross-section dataset (photoabsorption / photodecay) in cm², where photodecay is the reaction's single ionization-or-dissociation channel.

These assets feed two SQLite tables in jaff.db, which is what JAFF actually queries at runtime (see Codebase Structure):

• photo_reaction_cross_sections — one row per reaction, with a photo_absorption flag, a decay_type ("ionization" / "dissociation") and a file.hdf5::<group> pointer into the Leiden or NORAD HDF5 file.
• verner_cross_sections — the Verner analytic σ(E) expression as a SymPy-parseable string (symbol E, photon energy in eV, σ in cm²).

What lands on the reaction

Which database supplies a photoionization reaction is governed by pi_database (norad by default, or verner / leiden, globally or per reaction). The lookup order is the chosen database, then NORAD, then Verner, then Leiden (duplicates skipped); any fallback is logged once as a summary after loading. Photodissociation always uses Leiden.

For tabulated sources (Leiden / NORAD), reaction.xsecs_dict is an XsecsProps dict carrying the photon_energy grid (eV) plus photo_absorption and the single photodecay channel (cm² arrays, or None); _equations["decay_type"] records whether that channel is ionization or dissociation. The radiation integrator integrates these arrays numerically over each band.

Verner fits are kept symbolic: the dict then carries the SymPy Piecewise in xsecs_dict["photodecay_expr"] (symbol E, eV, σ in cm²) and the arrays are None; the integrator integrates the expression analytically per band. The database actually used is recorded in xsecs_dict["database"] ("norad", "verner" or "leiden").

Shielding

Photo-rates computed from cross sections assume an unattenuated radiation field. In a real cloud the photons that drive a reaction are absorbed by intervening gas — including the dissociating species shielding itself ("self-shielding"). JAFF models this with a dimensionless shielding factor \(f_\text{sh} \in [0, 1]\) that multiplies the rate coefficient:

The factor is applied to every band of a reaction (it multiplies each band's \(k_i\) before the band sum), so a single \(f_\text{sh}\) attenuates the whole reaction rate. \(f_\text{sh}\) is a symbolic expression in runtime quantities — column densities (ncol_<species>) and, for H2, the velocity dispersion (vdisp) — that the host code supplies at integration time.

When a shielding function depends on several shielding species, the per-species factors are multiplied together:

Enabling shielding

Shielding is opt-in per reaction, declared in jaffgen.toml under [network.reactions."<serialized>".shielding] (see the configuration reference). The reaction must be a photo-reaction; the type key selects the shielding function (default "leiden"). The serialized key contains . separators and so must be quoted in the TOML header.

# Leiden tabulated line shielding for CO photodissociation
[network.reactions."CO._PHOTON__C.O".shielding]
type        = "leiden"
radiation   = "ISRF"           # radiation-field subgroup (default "ISRF")
shielded_by = ["self", "H2"]   # shielding species; "self" = the reactant (CO)

# H2 self-shielding via the Hartwig et al. (2015) fit
[network.reactions."H2._PHOTON__H.H".shielding]
type      = "hg2015"
min_ncol  = 1.0e-35            # optional floors (see below)
min_vdisp = 1.0e-20

Shielding types

Leiden tabulated shielding (type = "leiden")

The shielding factor is read from the collapsed Leiden tables (data/shielding/leiden.hdf5, one group per reaction, downloaded on first use alongside the cross sections). Parameters:

For each species in shielded_by, JAFF emits a per-reaction shielding_<reaction>.hdf5 table next to the generated code and one interpolation call interp_<index>_shielding_<species>(ncol_<species>); the total factor is their product. "self" resolves to the reaction's reactant, so it interpolates over that species' own column density (e.g. ncol_CO for CO._PHOTON__C.O).

H2 self-shielding (type = "db1996" / "hg2015")

Both apply only to the H2._PHOTON__H.H (H2 → H + H) reaction and evaluate the standard three-term analytic fit

with \(x = N_{\text{H2}}/5\times10^{14}\,\text{cm}^{-2}\) and \(b_5 = b/10^5\,\text{cm s}^{-1}\), where the Doppler parameter \(b = \sqrt{2}\,\sigma_v\). The exponent \(\alpha\) is the only difference between the two: \(\alpha = 2\) for db1996, \(\alpha = 1.1\) for hg2015. The runtime inputs are the H2 column density ncol_H2 and velocity dispersion vdisp.

Self-consistent photoelectric field (chi_pe)

Grain photoelectric heating is driven by far-UV photons in the photoelectric band — from the grain work function (6 eV) up to the hydrogen ionisation edge (13.6 eV). Rate and heating expressions parametrise this by chi_pe: the local radiation field in that band, scaled to a reference (Draine/ISRF) field.

The band edges are configurable via DustProps(pe_threshold_low=..., pe_threshold_high=...) (defaults 6.0 and 13.6 eV), or the matching pe_threshold_low / pe_threshold_high keys in the [network.dust] table.

Whereas the plain chi symbol is a runtime input (an assumed external field), chi_pe is computed from the network's own radiation bands, so photoelectric heating stays consistent with the radiation transport actually being solved.

Definition

chi_pe is the ratio of two energy densities in the [6, 13.6] eV band:

• numerator — the energy density carried by the network's radiation bands, summing over each band the fraction of its energy that overlaps the photoelectric band. In photon-density mode each band's photden is multiplied by the band-average photon energy eavg to get an energy density.
• denominator — the energy density of the reference background_field in the same band, obtained from its tabulated photon-flux spectrum.

Both are energy densities in erg cm⁻³, so chi_pe is dimensionless.

Enabling it

chi_pe needs radiation transport and the dust module both on. From Python:

from jaff.physics import RadiationProps, DustProps

net = Network(
    "networks/GOW/GOW.jet",
    radiation_props=RadiationProps(
        bands=[6.0, 11.2, 13.6, "inf"],   # band edges in eV
        background_field="draine",         # reference field for the scaling
    ),
    dust_props=DustProps(),                # enables the dust module
)

From a jaffgen.toml add a [network.dust] table and a background_field key under [network.radiation]. If either radiation or dust is missing when a rate references chi_pe, generation aborts with a ParserError.

Under a reduced speed of light (rsl / RadiationProps(c="c_hat")), chi_pe is scaled by c_reduced / c so the photoelectric field stays consistent with the reduced-speed radiation transport (commit cddf25d).

Reference the symbol in any rate or .jfunc expression exactly like chi; JAFF substitutes the computed expression when the network is built.

Python API

from jaff import Network
from jaff.physics import RadiationProps

# Enable photochemistry by declaring radiation bands
net = Network(
    "networks/h_photoionization/h_photo.jet",
    radiation_props=RadiationProps(
        bands=[13.6, "inf"],   # band edges in eV
        profile_index=0,
        mode="nph",
    ),
)

# One spectral index per band (here two bands: 13.6–20 eV and 20–100 eV)
net_pb = Network(
    "networks/h_photoionization/h_photo.jet",
    radiation_props=RadiationProps(bands=[13.6, 20.0, 100.0], profile_index=[0, 1]),
)
[grp.profile_idx for grp in net_pb.radiation.groups]   # [0, 1]

# Inspect photochemical reactions
photo = net.reactions.photo_reactions()
for rxn in photo:
    print(rxn.verbatim)

# Cross sections are attached at load time
rxn = photo[0]
rxn.xsecs_dict["photon_energy"]     # eV grid
rxn.xsecs_dict["photodecay"]        # cm² array (or None)
rxn.plot_xsecs()                    # visualise σ(E)

# Symbolic radiation ODEs
print(net.sradodes())

Cross-section lookup is also exposed directly via the jaff.physics.Photochemistry class. Constructing it downloads the cross-section and line-shielding data files on first use (cached thereafter), so instantiate once and reuse:

from jaff.physics import Photochemistry

photo = Photochemistry()
photo.get_xsec(rxn)         # XsecsProps from the tabulated databases
photo.get_verner_xsec(rxn)  # analytic Verner σ(E) (sympy) or None