pynonthermal

pynonthermal is a Spencer-Fano equation solver for non-thermal electron energy deposition in plasmas. It computes how deposited energy is partitioned into heating, ionisation, and excitation, and provides non-thermal ionisation and excitation rate coefficients.

Contents

• Quick start
• Installation
• Usage
• Units and conventions
• Example output
• Method background
• Cross-section dataset choice
• Citing pynonthermal
• License

Quick start

import pynonthermal

sf = pynonthermal.SpencerFanoSolver(emin_ev=1.0, emax_ev=3000.0, npts=4096)

# Add ions that can be non-thermally ionised.
# Here: O II (ion_stage=2, i.e. charge +1) with number density in cm^-3.
sf.add_ionisation(Z=8, ion_stage=2, n_ion=1.0e8)

# Solve for a deposition rate density in eV s^-1 cm^-3.
sf.solve(depositionratedensity_ev=1.0e8)

print("heating fraction:", sf.get_frac_heating())
print("total ionisation fraction:", sf.get_frac_ionisation_tot())
print("total excitation fraction:", sf.get_frac_excitation_tot())
print("sum of fractions:", sf.get_frac_sum())
print("ionisation rate coeff [s^-1]:", sf.get_ionisation_ratecoeff(Z=8, ion_stage=2))

The quickstart notebook contains a fuller worked example, and can be launched on Binder:

Installation

Released package (recommended for most users):

pip install pynonthermal

Development install with uv:

git clone https://github.com/lukeshingles/pynonthermal.git
cd pynonthermal
uv sync --frozen
source ./.venv/bin/activate
uv pip install --editable .
prek install

Usage

Typical solver workflow:

1. Create SpencerFanoSolver with an energy grid (emin_ev, emax_ev, npts).
2. Add ionisation channels with add_ionisation(Z, ion_stage, n_ion).
3. Optionally add excitation channels with add_excitation(...) or add_ion_ltepopexcitation(...).
4. Call solve(depositionratedensity_ev=..., override_n_e=...).
5. Query outputs such as get_frac_heating(), get_frac_ionisation_tot(), get_frac_excitation_tot(), get_ionisation_ratecoeff(Z, ion_stage), and get_excitation_ratecoeff(Z, ion_stage, transitionkey).

Units and conventions

• Energies are in eV.
• Number densities are in cm^-3.
• Cross sections are in cm^2.
• ion_stage = charge + 1 (for example, Fe I has ion_stage=1, Fe II has ion_stage=2).
• depositionratedensity_ev in solve() is in eV s^-1 cm^-3.
• get_ionisation_ratecoeff() and get_excitation_ratecoeff() both return rates in s^-1, and both scale linearly with depositionratedensity_ev.

Example output

The following plot shows the energy distribution of contributions to ionisation, excitation, and heating for a pure oxygen plasma (electron fraction 0.01), reproducing Figure 2 of Kozma and Fransson (1992). The area under each curve gives the fraction of deposited energy in that channel.

This figure is generated from the same solver setup demonstrated in the quickstart workflow.

Method background

When high-energy leptons (electrons and positrons) are injected into a plasma, they lose energy through ionisation, excitation, and Coulomb interactions with thermal electrons. Tracking these rates is important, for example, in late-time Type Ia supernova modelling.

The numerical solver is similar to the Spencer-Fano implementation in the ARTIS radiative transfer code (Shingles et al. 2020), itself an independent implementation of Kozma and Fransson (1992), based on Spencer and Fano (1945). A similar approach is used in CMFGEN.

Impact ionisation cross sections use fits from Arnaud and Rothenflug (1985) and Arnaud and Raymond (1992) for Z=1 to 28 (H to Ni). Heavier elements use the approximation of Axelrod (1980), Eq. 3.38 with Lotz (1967).

If internal level/transition data are used (for example, via add_ion_ltepopexcitation()), they are imported from the CMFGEN atomic data compilation. See source data files for references.

Cross-section datasets

Ionization cross sections from H (Z=1) to Ni (Z=28) are sourced from the analytical fits to data from M. Arnaud & R. Rothenflug (1985, A&AS, 60, 425), with updates to Fe from M. Arnaud & J. Raymond (1992, ApJ, 398, 394). For heavier elements (Z>28), the approximation of Axelrod (1980, PhD thesis, Eq. 3.38) is used, with Lotz (1967, Z. Phys., 206, 205) for the required parameters.

Advanced usage: custom excitation cross sections

You can supply your own excitation cross section table:

sf.add_excitation(Z, ion_stage, levelnumberdensity, xs_vec, epsilon_trans_ev, transitionkey=(lower, upper))

• Z: atomic number.
• ion_stage: one more than ion charge.
• levelnumberdensity: population density of the lower level (cm^-3), non-negative.
• xs_vec: NumPy array of cross sections (cm^2), non-negative, defined at every energy in sf.engrid (eV).
• epsilon_trans_ev: transition energy (eV). Must be positive and no greater than emax_ev, since no electron the solver represents could otherwise drive the transition.
• transitionkey: any unique key within the ion, used to retrieve the excitation rate coefficient.

Transitions below emin_ev are allowed here, but add_ion_ltepopexcitation() drops them: Kozma and Fransson (1992) take every electron below emin_ev to have thermalised, so that energy is accounted for as heating instead.

Retrieve the non-thermal excitation rate coefficient [s^-1] with:

nt_exc = sf.get_excitation_ratecoeff(Z, ion_stage, transitionkey)

Citing pynonthermal

If you use pynonthermal, please cite it via the Zenodo record.

License

Distributed under the MIT license. See LICENSE for details.