pynonthermal

pynonthermal is a Python solver for the Spencer-Fano equation, which describes the energy distribution of non-thermal (fast) electrons slowing down in a plasma. When high-energy leptons — such as the Compton, photoelectric, and pair-production electrons and positrons produced by radioactive decay in supernova ejecta — are injected into a partially ionised gas, they lose energy through three competing channels: Coulomb heating of the free thermal electrons, collisional ionisation, and collisional excitation of bound states.

Given a set of ions (with number densities) and an energy deposition rate, pynonthermal computes:

• the degradation spectrum y(E) of the non-thermal electron population,
• the fraction of deposited energy going to heating, ionisation, and excitation (per channel and per ion),
• non-thermal ionisation rate coefficients for each ion and excitation rate coefficients for individual bound-bound transitions, ready to be used in non-LTE plasma modelling.

These quantities are important, for example, in modelling the late-time spectra and light curves of Type Ia and core-collapse supernovae, where non-thermal ionisation can dominate over photoionisation. The solver follows the method of Kozma & Fransson (1992) (see Method background for details and further references) and ships with the atomic data needed to run out of the box: ionisation cross sections for a wide range of ions, and level/transition data for bound-bound excitation.

Contents

• Installation
• Quick start
• Usage guide
• Where the ion populations come from
• Complete example: pure-oxygen plasma
• Units and conventions
• Method background
• Cross-section datasets
• Advanced usage: custom cross sections
• Citing pynonthermal
• License

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

Run the test suite with:

uv run -- python3 -m pytest

Quick start

Make a solver, add the elements of the gas, solve, then read the results:

import pynonthermal

with pynonthermal.SpencerFanoSolver() as sf:
    # O II (ion_stage=2, i.e. charge +1) at a number density of 1e8 cm^-3
    sf.add_element(8, ion_densities={2: 1.0e8})

    sf.solve(deposition_ev_per_s_per_cm3=1.0e8)  # the rate of energy deposition per volume

    print("heating fraction:", sf.get_frac_heating())
    print("ionisation fraction:", sf.get_frac_ionisation_tot())
    print("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 solver is a builder: each call adds ions, ionisation channels, or excitations to the matrix, and solve() solves it. The with block is optional and only scopes the solver.

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

Usage guide

1. Make the solver

sf = pynonthermal.SpencerFanoSolver(emin_ev=0.1, emax_ev=16000.0, npts=4096)

• emin_ev, emax_ev: the bounds of the uniform energy grid in eV (defaults 0.1 and 16000.0). An electron that degrades below emin_ev is taken to have thermalised, so its energy counts as heating. Every ionisation potential must lie above emin_ev, and a ValueError says which lower emin_ev to use. The default emax_ev covers the K-shell ionisation of iron at about 7 keV; lower it to 3000 eV to match Kozma & Fransson (1992). emin_ev also sets the highest free electron density that the solver accepts, because the Coulomb logarithm of the loss function must stay positive at the bottom of the grid: about 7e16 cm^-3 at emin_ev=0.1, and about 7e19 cm^-3 at emin_ev=1.
• npts: the number of energy grid points (default 4096). More points cost memory and time; check get_frac_sum().
• verbose: print the setup, each added ionisation channel, and a per-ion, per-shell breakdown.
• use_ar1985: use the original Arnaud & Rothenflug (1985) ionisation cross sections (see Cross-section datasets).
• heating_only_approximation: leave the excitation and ionisation loss terms out of the matrix and solve with the heating loss alone. The rates still follow from that approximate solution, so the fractions do not sum to one.

2. Set the temperature and the atomic data

sf.set_temperature(6000)  # K, for the LTE populations of the excitation levels
sf.set_atomic_data(use_collstrengths=True, maxnlevelslower=5, maxnlevelsupper=250)

The solver has one temperature. Set it before any excitation with LTE level populations; it is not needed otherwise.

set_atomic_data() chooses the level data for the excitations and how to build their cross sections. adata_polars takes your own level/transition table in the format of artistools.atomic.get_levels(); the other options default to the ARTIS values (collision strengths where available, and transitions from the lowest 5 levels up to the lowest 250). Without this call the solver uses the internal database with those defaults.

3. Add the elements

sf.add_element(8, 1.0e10, ion_fractions={1: 0.99, 2: 0.01}, excitation=True)
sf.add_element(26, 1.0e9, recomb_ratecoeffs={2: 1.0e-11, 3: 1.5e-11, 4: 3.0e-11, 5: 6.0e-11})

add_element() takes:

• Z: the atomic number, and n_elem: the number density of the element in cm^-3, summed over its ion stages. Every rule needs n_elem except ion_densities, which gives the densities themselves.
• the population rule: exactly one of ion_densities, ion_fractions, or recomb_ratecoeffs — see the next section.
• excitation: also add the bound-bound excitations of every ion stage that has level data, with LTE level populations at the temperature of the solver. Every stage gets the built-in ionisation cross sections either way.
• builtin_channels: set it to False to leave the built-in ionisation cross sections out and give every channel yourself (see custom cross sections).

To add one ion at a time instead, call add_ionisation(Z, ion_stage, n_ion) for the built-in shells of that ion, and add_ion_excitation(Z, ion_stage, n_ion) for its bound-bound excitations with LTE level populations. Both take n_ion=None for an ion whose population the solver already holds.

4. Solve

sf.solve(deposition_ev_per_s_per_cm3=1.0e8)

• deposition_ev_per_s_per_cm3: the rate of energy deposition per volume in eV s^-1 cm^-3 (positive and finite). With fixed populations the energy fractions do not depend on it and the rate coefficients scale linearly with it; with recomb_ratecoeffs the populations depend on it too.
• balance_tol: the relative tolerance of the population ratios of an element with recomb_ratecoeffs (default 1e-4).

The free electron density comes from the ion charges of ionpopdict, so it counts only the electrons of the ions that the solver holds. Give it yourself with sf.override_n_e(n_e) before solve(), for example when species that are not in the solver also give electrons:

sf.override_n_e(2.5e6)  # cm^-3; None takes it from the ion charges again

It works with every population rule. With recomb_ratecoeffs it replaces charge neutrality: solve() then finds the populations at your density, and they do not have to be neutral with it. The value holds until another call changes it, and a call after solve() discards the solution, so solve() must run again. (The override_n_e argument of solve() is deprecated.) The iron notebook shows what a given density does to a balance.

5. Read the results

Energy fractions, as shares of the deposited energy:

sf.get_frac_heating()  # to heating of the thermal electrons
sf.get_frac_ionisation_tot()  # to ionisation, over all ions
sf.get_frac_excitation_tot()  # to excitation, over all ions
sf.get_frac_sum()  # the sum; ~1.0 when the grid resolves every ionisation and excitation
sf.get_frac_ionisation_ion(Z, ion_stage)  # one ion's share
sf.get_frac_excitation_ion(Z, ion_stage)

Populations, densities, and rates:

sf.get_ion_fractions(Z)  # {ion_stage: fraction of the element}
sf.ionpopdict  # {(Z, ion_stage): number density [cm^-3]}
sf.get_n_e()  # free (thermal) electron density [cm^-3]
sf.get_n_e_nt()  # non-thermal electron density [cm^-3]
sf.get_n_ion_tot()  # total nuclei [cm^-3]
sf.balance_iterations  # iterations that a recomb_ratecoeffs balance took

sf.get_ionisation_ratecoeff(Z, ion_stage)  # [s^-1]
sf.get_excitation_ratecoeff(Z, ion_stage, transitionkey)  # [s^-1]
sf.excitationlists[(Z, ion_stage)]  # the transitions of that ion, keyed by transition key
sf.get_eff_ionpot(Z, ion_stage)  # effective ionisation potential [eV]

Multiply get_ionisation_ratecoeff() by the ion number density for ionisations per second per cm^3, and get_excitation_ratecoeff() by the lower level's population density for excitations per second per cm^3. For the excitations of add_element(excitation=True) and of add_ion_excitation() the key is (lower_level_index, upper_level_index), for example (0, 8).

Each value of excitationlists[(Z, ion_stage)] is a pynonthermal.ExcitationTransition, which holds the lower level population levelnumberdensity in cm^-3, the cross sections xs_vec in cm^2 on sf.engrid, and the transition energy epsilon_trans_ev in eV.

The solution itself is sf.yvec over sf.engrid, both read-only arrays. With verbose=True the solver prints the per-ion and per-shell breakdown as it analyses the solution.

6. Plot the solution

sf.plot_yspectrum()  # degradation spectrum y(E)
sf.plot_channels(xscalelog=True)  # energy going to each channel vs electron energy
sf.plot_spec_channels("channels.pdf")  # both panels in one figure, saved to file

Each method shows the figure interactively, or saves it when outputfilename is given; plot_yspectrum() and plot_channels() also accept a Matplotlib axis to draw into an existing figure.

Where the ion populations come from

Every add_element() call gives exactly one of three rules. A future non-LTE rule will be a fourth keyword.

ion_densities

sf.add_element(26, ion_densities={2: 3.0e5, 3: 7.0e5})

The number density in cm^-3 of each ion stage, keyed by ion stage. n_elem is their sum, so do not give it as well. This is usually what a plasma code already holds.

ion_fractions

sf.add_element(26, 1.0e6, ion_fractions={2: 0.3, 3: 0.7})

The same populations as a share of n_elem. The fractions must lie between 0 and 1 and sum to one.

recomb_ratecoeffs

sf.add_element(8, 1.0e10, recomb_ratecoeffs={2: 3.0e-13, 3: 3.0e-12, 4: 1.0e-11})

The recombination rate coefficients in cm^3 s^-1, keyed by the ion stage that recombines. For each pair of adjacent stages i and i+1 the balance is n_i Gamma_i = n_{i+1} n_e alpha_{i+1}, where Gamma_i is the non-thermal ionisation rate coefficient of stage i from the Spencer-Fano solution and alpha_{i+1} is the coefficient you give. The chain runs from one below the lowest key to the highest key, so the example is O I to O IV.

A coefficient outside 1e-16 to 1e-8 cm^3 s^-1 raises a warning (RECOMB_RATECOEFF_MIN_WARN and RECOMB_RATECOEFF_MAX_WARN). The published radiative and dielectronic fits stay inside that range, so a value outside it is nearly always a unit error: a coefficient in m^3 s^-1 is 1e-6 of the same coefficient in cm^3 s^-1. The value is still used.

The solution depends on the ion densities, so solve() iterates: it solves the equation, updates the densities from the balance and the free electron density from charge neutrality, and repeats until the population ratios agree to balance_tol. Typical cases converge in about 5 to 10 iterations; a RuntimeError reports a balance that did not converge within 100. sf.balance_iterations says how many it took.

Points to note:

• The balance includes only non-thermal ionisation and the recombination that you give. It does not include thermal collisional ionisation, photoionisation, or charge exchange. The ion fractions therefore depend on the deposition rate density, unlike the fixed-population case.
• The top stage of the chain is a sink: its ionisation is an energy loss in the matrix, but the ions it makes have no stage to go to. A warning is raised if the ionisation rate out of the top stage exceeds 1 % of the total ionisation rate of the element, because about that fraction of the element then belongs in a higher stage. Extend the chain with a rate coefficient for the next stage.
• The free electron density comes from charge neutrality, unless override_n_e() gives it.

The functions behind the balance are in pynonthermal.ionbalance: get_ion_fractions(), solve_charge_neutral_n_e_ratios(), and the general root find solve_charge_neutral_n_e(), which takes any charge density function that does not increase with the free electron density.

The Saha equation as a comparison

The Saha equation is not a population rule, because it describes a gas whose ionisation is thermal, and a gas with non-thermal ionisation is not in local thermodynamic equilibrium. It is still worth comparing against, so pynonthermal.ionbalance gives it as a function of its own:

fractions = pynonthermal.ionbalance.get_saha_ion_fractions(26, [1, 2, 3, 4, 5], 6000.0, n_elem=1.0e6)

This needs no Spencer-Fano solution. For each pair of adjacent stages, n_{i+1} n_e / n_i = 2 (U_{i+1} / U_i) (2 pi m_e k_B T / h^2)^(3/2) exp(-chi_i / (k_B T)), with the ionisation potentials chi_i from the NIST table. n_elem gives the charge-neutral free electron density of that element alone; n_e instead fixes it, for the comparison at the density of a solution. The partition functions U_i come from the LTE level populations of the level data; the built-in data covers He, O, and Fe. For other elements give them as partfuncs={ion_stage: U, ...}, or supply a level table in adata_polars. The bare nucleus (ion_stage = Z + 1) has a partition function of 1. get_saha_factor() gives one pair's ratio coefficient by itself.

The iron ionisation balance notebook is a worked example: the ion fractions of iron in the core of a Type Ia supernova at 250 days, with the deposition rate from the 56Co decay, a comparison with the Saha equation, and the evolution from 150 to 400 days.

Complete example: pure-oxygen plasma

This reproduces Figure 2 of Kozma & Fransson (1992): a pure-oxygen plasma with electron fraction x_e = 0.01, including both ionisation and excitation channels. With verbose=True the solver prints its setup and a per-ion, per-shell breakdown as it runs.

import pynonthermal

n_e = 1e8  # free electron density [cm^-3]
x_e = 1e-2  # ionisation fraction n_OII / (n_OI + n_OII)
n_oxygen = n_e / x_e

# the grid is theirs rather than the default: emin_ev=1 is their low-energy cutoff E_0,
# and emax_ev=3000 is the top of their energy range.
with pynonthermal.SpencerFanoSolver(emin_ev=1, emax_ev=3000, npts=4096, verbose=True) as sf:
    sf.set_temperature(6000)  # K, for the LTE level populations
    sf.add_element(8, n_oxygen, ion_fractions={1: 1 - x_e, 2: x_e}, excitation=True)

    # with fixed ion densities, any positive deposition rate works here: the energy fractions
    # are independent of it (with recomb_ratecoeffs they would not be).
    sf.solve(deposition_ev_per_s_per_cm3=2950.49 * n_oxygen)

    sf.plot_channels(xscalelog=True)

The resulting plot shows the energy distribution of contributions to ionisation, excitation, and heating; the area under each curve gives the fraction of deposited energy in that channel:

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).
• deposition_ev_per_s_per_cm3 is in eV s^-1 cm^-3.
• get_ionisation_ratecoeff() and get_excitation_ratecoeff() both return rates in s^-1.
• The recomb_ratecoeffs of add_element() are in cm^3 s^-1, keyed by the ion stage that recombines.

Method background

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, ApJ, 390, 602), based on the electron slowing-down equation of Spencer and Fano (1954, Phys. Rev., 93, 1172). A similar approach is used in CMFGEN.

The integral form of the Kozma and Fransson degradation equation (their equation 7) is discretised on a uniform energy grid as an upper-triangular matrix equation and solved by back-substitution from the highest energy downward. The SpencerFanoSolver class docstring maps each term of the equation to the method that implements it, and the code comments cite the specific Kozma and Fransson equations at each site. The secondary-electron energy distribution follows Opal, Peterson and Beaty (1971) as applied by Kozma and Fransson, and the energy loss rate to thermal electrons uses their Coulomb-logarithm prescription (after Schunk and Hays 1971).

If internal level/transition data are used (for example, via add_ion_excitation()), they are imported from the CMFGEN atomic data compilation (see the source data files for references), with excitation cross sections computed from the tabulated collision strengths (Li, Dessart and Hillier 2012, equation 11) or, for permitted transitions without one, from the oscillator strength via the van Regemorter (1962) approximation with the g-bar factor of Mewe (1972), as described in Shingles et al. (2020, section 2.5).

Cross-section datasets

Ionization cross sections from H (Z=1) to Ni (Z=28) use the shell-resolved analytical fits compiled by Arnaud and Rothenflug (1985, A&AS, 60, 425), with updates to Fe from Arnaud and Raymond (1992, ApJ, 398, 394). For heavier elements (Z>28) and any other ions missing from the fit data, the approximation of Axelrod (1980, PhD thesis, Eq. 3.38) is used — the high-energy limit of the Lotz (1967, Z. Phys., 206, 205) formula with relativistic corrections — with subshell binding energies from Lotz (1970, J. Opt. Soc. Am., 60, 206).

Passing use_ar1985=True to the solver selects the original Arnaud and Rothenflug (1985) compilation without the Fe updates, which can be useful for comparison with older published results.

Advanced usage: custom cross sections

Give a cross section as a function of an array of energies in eV that returns cross sections in cm^2 (the type pynonthermal.CrossSectionFunc). The solver calls it on its own grid, and between the grid points where it needs to, so the plasma does not depend on the energy grid of the solution. An array at every energy of sf.engrid is accepted too, but then the plasma is tied to that grid and the solver can only interpolate between the points.

A custom cross section follows the same path through the solver as a built-in one, so the matrix, the energy fractions, and the rate coefficients stay consistent.

import numpy as np
import pynonthermal


def my_ionisation_xs(en_ev):
    return np.interp(en_ev, my_en_ev, my_xs_cm2, left=0.0, right=0.0)


with pynonthermal.SpencerFanoSolver() as sf:
    sf.add_element(8, ion_densities={2: 1.0e8})
    # this keeps the built-in shells and adds one channel; builtin_channels=False replaces them
    sf.add_ionisation_channel(8, 2, None, ionpot_ev=35.0, xs_vec=my_ionisation_xs, channelkey="mine")
    sf.add_excitation(8, 2, None, my_excitation_xs, epsilon_trans_ev=20.0, transitionkey=(0, 3), levelpopfrac=0.9)

add_ionisation_channel() takes:

• Z and ion_stage: the ion that the channel ionises, and n_ion: its number density in cm^-3, or None for an ion whose population the solver already holds.
• ionpot_ev: the ionisation potential in eV. It must lie between emin_ev and emax_ev, and the cross section must be zero at and below it. Any value in that range is allowed, so a channel need not be a subshell of the built-in table; a total ionisation cross section for the ion works too.
• xs_vec: the cross section, non-negative and finite.
• channelkey: any key that is unique within the ion, for the verbose output.

add_excitation() takes:

• Z, ion_stage, and levelnumberdensity: the population of the lower level in cm^-3. Give None instead, with levelpopfrac between 0 and 1, for an ion whose population the solver already holds. The level population then follows the ion population, whether it is fixed or comes from a balance.
• epsilon_trans_ev: the transition energy in eV. It must be positive and no greater than emax_ev, since no electron the solver represents could otherwise drive the transition. Transitions below emin_ev are allowed here, but add_element(excitation=True) drops them: Kozma and Fransson (1992) take every electron below emin_ev to have thermalised, so that energy is accounted for as heating instead.
• xs_vec, and transitionkey: the key to pass to get_excitation_ratecoeff().

sf.calculate_N_e() integrates over a domain just above the ionisation potential that is narrower than one grid cell. A cross section given as a function is called there; an array can only be interpolated, so resolve that region with npts if it matters for your ion. The term it feeds is the energy that thermalises below emin_ev, which is a small part of the heating fraction.

The solver keeps the Lorentzian secondary-electron distribution of Kozma and Fransson (1992, equation 4), whose width comes from pynonthermal.collion.get_J(). The matrix fill integrates that distribution analytically, so its shape is not adjustable.

Citing pynonthermal

If you use pynonthermal, please cite it via the Zenodo record. Please also consider citing the papers describing the method: Kozma and Fransson (1992) and Shingles et al. (2020).

License

Distributed under the MIT license. See LICENSE for details.