Create a class to represent ionization state distributions

[namurphy]

There are a variety of situations for which plasma physicists need to keep track of ionization state distributions (e.g., what fraction of iron is neutral, singly charged, doubly charged, etc.). We should create a base class that is able to store and access this information. This can probably be related to the Particle class.

This class should be able to be extended to have functions that can, e.g., predict the resulting spectrum. There should probably also be a standard way to include information such as the temperature (or a more general particle distribution function) and number density, as well as relative abundances of elements.

We might want to create a way to more easily iterate across charge states and elements.

@lemmatum and @wtbarnes - This seems like a great thing to collaborate on among SpectroscoPy, fiasco, PlasmaPy, and perhaps other existing spectroscopy packages. I've been meaning to look more deeply into how fiasco and specutils store these data.

Cross-reference: https://github.com/NEI-modeling/NEI/issues/18

[wtbarnes]

When using the CHIANTI database, one can obtain the ionization fractions 1 of 2 ways,

• Access pre-computed ionization fractions stored in CHIANTI tabulated at different temperatures (there are multiple datasets available for the same elements/ions)
• Compute the ionization fractions yourself using the ionization and recombination rates

Note that in both cases, it is assumed that the rates (and by extension the fractions) are not density-dependent and that the ionization and recombination rates are in equilibrium (though of course the rates can be used to compute the time-dependent ionization states).

In fiasco, you can compute the the ionization fractions as a function of temperature using either method, e.g. for iron,

import numpy as np
import astropy.units as u
import fiasco
import matplotlib.pyplot as plt
iron = fiasco.Element('iron', np.logspace(4,9,100)*u.K)
# Compute ionization fractions using ionization and recombination rates
ioneq = iron.equilibrium_ionization()
for i,ion in enumerate(iron):
    # Method 1--precomputed data from database
    plt.plot(ion.temperature, ion.ioneq, f'C{i%10}', ls='--')
    # Method 2--computed from rates
    plt.plot(ion.temperature, ioneq[:,i], f'C{i%10}', ls='-', label=ion.ion_name)
plt.xscale('log')
plt.show()

[lemmatum]

@wtbarnes When computing rates is this done simply via Saha-Boltzmann or do you have a more sophisticated model? Have these been checked against other codes/databases?

Another way to get ionization rates would be to query NIST. I believe they have a temperature and density sensitive method for getting these rates. In any case, if density sensitive rates aren't yet implemented in fiasco, perhaps we can make that an issue in either fiasco or SpectroscoPy (or both and then compare) since density sensitivity is quite important... especially for the plasmas I'm interested in :smile_cat: .

[wtbarnes]

@lemmatum CHIANTI is primarily used for diffuse, high-temperature optically-thin astrophysical plasmas, in particular the solar corona, where the assumption of local thermodynamic equilibrium is not likely to hold. Because of this, the ionization fractions must be determined by a balance between the ionization and recombination rates,

Section 4 of this paper provides a nice overview of how the ionization fractions are generally treated in the low-density, high-temperature solar corona.

In CHIANTI (and by extension in fiasco), the following processes are included when computing the ionization equilibrium: direct ionization, excitation-autoionization, radiative recombination, and dielectronic recombination (all functions of temperature) such that the balance equation becomes,

Thus, finding the populations is just a matter of inverting this Z+1-by-Z+1 (as a function of T!) matrix. The tricky part of course is actually calculating these rates properly. fiasco (and originally the CHIANTI IDL software and ChiantiPy) employs several specialized methods from the literature to do this. Sections 2 and 3 of this paper summarize the different methods used for calculating each of these terms.

Presumably some density-dependence could be added to these and I know there is some interest in doing so (particularly in the dielectronic recombination term).

Sorry this is a bit long winded. The TL;DR is that yes, CHIANTI (and thus fiasco) employs more specialized methods for calculating the ionization fractions, though these methods may not generalize well to other plasmas as they are methods derived for somewhat specific parameters (and in some cases specific ions).

I believe that some informal comparisons have been made between CHIANTI and other databases (e.g. ADS), though I don't know about any thorough cross-checks being done. @namurphy and I discussed this a bit back in the summer. A common API between multiple databases would certainly make this much less painful.

For what it is worth, I've compared the output of fiasco with the IDL procedures and ChiantiPy and they are consistent.

[lemmatum]

@wtbarnes thanks for the details! They'll be useful in figuring out the scope of what needs to be done in SpectroscoPy without duplicating effort. I've opened an issue in SpectroscoPy to deal with this feature.

[namurphy]

Closed by #404.