Saha Ionization Model
Athelas models ionization equilibrium with the non-degenerate Saha-Boltzmann
equations based on the formalism of Zaghloul et al. (2000).
This model is used to compute ionization fractions and the mean charge state of
each atomic species. The ionization state enters the equation of state
through contributions to the internal energy and pressure and is therefore
solved repeatedly as part of the temperature inversion.
It is a large computational expense.
The Saha-Boltzmann equations have the form
where \(n_s\) is the number density of atoms in the \(s\)-th ionization state, \(g_s\) is the associated statistical weight, \(m_e\) is the electron rest mass, \(h\) is Planck’s constant, \(k_B\) is Boltzmann’s constant, \(\chi_s\) is the ionization energy for the \(s \to (s+1)\) process, and \(\mathbb{Z}\) is the atomic number of the species.
The Saha-Boltzmann equations can be combined with the condition of charge neutrality
and conservation of nuclei
where \(n_k\) is the number density of nuclei of species \(k\), and formed into a single set of transcendental equations that can be solved for the mean charge \(\bar{\mathbb{Z}}\). Let us denote the fraction of atoms in the \(s\)-th ionization state as \(y_s = n_s/n_k\) and express the mean charge as \(\bar{\mathbb{Z}} = n_e/n_k\) we can, following Zaghloul , reexpress the above equations in the following forms
and observe the recurrence relation
Substituting this into (5) we get an expression for the neutral fraction
Finally, combining the relations for the ionization fractions \(y\) with the condition for conservation (4) of nuclei we arrive at a transcendental equation for \(\bar{\mathbb{Z}}\)
Equation (9) may be solved iteratively for \(\bar{\mathbb{Z}}\) using a any
root finding algorithm. Athelas implements a couple of method for finding
\(\bar{\mathbb{Z}}\).
Once \(\bar{\mathbb{Z}}\) is known then the constribution from this species to the electron number density can be tallied as \(n_e = \bar{\mathbb{Z}}n_k\).
Numerical Methods
Athelas implements two methods for solving equation (9).
Both solvers use a Newton-Raphson iteration so derivatives are necessary.
The first, which we refer to as the linear model solves Eq. (9)
directly. The second, referred to as log in Athelas, reformaulates
Eq. (9) by taking natural logs and expressing \(\prod f\) as
exponentials.
The linear Saha solver is faster but can be numerically fragile
for large \(\mathbb{Z}\) and intended for light elements.
For the simulations that Athelas is designed to carry out
this is usually fine – it is stable through CNO. The logarithmic
solver is robust for all species and thermodynamic conditions,
at increased cost.
Linear
The original Saha solver in Athelas, which we refer to as linear,
solves Eq (9) directly using a Newton-Raphson iteration.
Before writing down the derivative \(F'(\mathbb{\bar{Z}})\) we note that
this formulation can be subject to numerical instabilities. An attempt to
circumvent this, inspired by the implementation in SNEC, is to introduce thresholds on
prospective electron number densities such that they are ignored. This serves
both as a guard to instability and improves performance.
We rewrite Eqs (9) and (8) as
\(s_{\rm min}\) and \(s_{\rm max}\) are set based on thresholds
ZBARTOL and ZBARTOLINV as compared to \(f_s/(\bar{\mathbb{Z}} n_k)\).
If \(f_s/(\bar{\mathbb{Z}} n_k) <\) ZBARTOL then we do not solve for the
\((s+1)\)-th ionization state and assume that state and all following are empty.
Similarly, if \(f_s/(\bar{\mathbb{Z}} n_k) >\) ZBARTOLINV we assume that the
\(s\)-th ionization state and all preceeding are empty.
With this in hand we may write down the derivative \(F'(\bar{\mathbb{Z}})\)
where we have defined
Log
The products of exponentials in Eq. (9) can overflow for sufficiently high-\(\mathbb{Z}\) atoms. To improve numerical stability, Athelas can express the Saha ladder in logarithmic space. Begin by compactifying Eq (9) as
where \(\mathcal{N}\) and \(\mathcal{D}\) can be identified from Eq (9). Taking the natural logarithm of \(F(\bar{\mathbb{Z}})\) and rearranging, we notice that a root of \(F\) corresponds to a root of
In what follows we will exploit the “LogSumExp” trick borrowed from machine learning. The core idea is to notice that
All that remains is to express \(\mathcal{N}\) and \(\mathcal{D}\) in a useful form. Note that the above can be written into an inline algorithm that doesn’t require a-priori knowledge of \(a_{\rm max}\).
We define
The transcendental equation for \(\bar{\mathbb{Z}}\) is rewritten in terms of \(x = \ln \bar{\mathbb{Z}}\). Defining the sums
With manipulation these can be identified as the numerator and denominator of Eq (9). The mean charge equation then becomes
This form avoids explicit products and powers of \(\bar{\mathbb{Z}}\), and all summations are evaluated using stable log-sum-exp techniques. The equation \(G(x) = 0\) is solved using a Newton–Raphson iteration in logarithmic space. This requires the derivative which, with respect to \(x\), can easily be written down:
where
Then we may recover \(\bar{\mathbb{Z}} = e^{x}\). The iteration typically converges in one or two steps.
Note
The minimum/maximum ionization state tricks used above in the linear
solver are not yet implemented in the log solver.
Tip
We recommend the linear solver for most use cases. It is slightly
faster and stable for most species of concern in the context of
supernova ejecta.
Atomic Data
We take ionization potentials and statistical weights from the NIST database.
Ionization potentials are stored in athelas/data/atomic_data_ionization.dat
and statistical weights in athelas/data/atomic_data_degeneracy_factors.dat.
These are provided for convenience but the user can supply their own.
Attention
We assume, as is common, that the ground states are the dominant contributors in the Saha-Boltzmann equation. Therefore, there is no temperature dependence in the partition functions and we simply use the ground state statistical weights.
Configuration
The ionization physics has a few controls in the [ionization] block.
The paths to the ionization potentials and degeneracy factors must
be specified. The provided data are stored in athelas/data and can be set
with fn_ionization and fn_degeneracy. These are required when
ionization is enabled. Then, the number of species to solve Saha ionization for
may be set with the ncomps option. Finally, the solver method may be set
as solver = linear or solver = log.
As is inferred above, we do not require that we do Saha ionization solves for all species present on the domain. Instead, we can do ionization for any number of species.
Note
We currently assume full ionization for species which we do not do Saha solves. This can be extended with little work if desired.
Note
The number of Saha species must be less than or equal
to the total number of species set in the composition block.
If neutrons are present in the compositional profile then this
inequality is strict.
Input Deck
Ionization is controlled in the input deck as follows:
[physics]
ionization = true
[ionization]
fn_ionization = "../data/atomic_data_ionization.dat" [Required]
fn_degeneracy = "../data/atomic_data_degeneracy_factors.dat" [Required]
ncomps = 3 # must be <= [composition.ncomps] !! [Required]
solver = "linear" # or "log" [Optional]