PyAtomDB: Extending the AtomDB Atomic Database to Model New Plasma Processes and Uncertainties

Abstract

The AtomDB project provides models of X-ray and extreme ultraviolet emitting astrophysical spectra for optically thin, hot plasma. We present the new software package, PyAtomDB, which now underpins the entire project, providing access to the underlying database, collisional radiative model calculations, and spectrum generation for a range of models. PyAtomDB is easily extensible, allowing users to build new tools and models for use in analysis packages such as XSPEC. We present two of these, the kappa and ACX models for non-Maxwellian and Charge-Exchange plasmas respectively. In addition, PyAtomDB allows for full open access to the apec code, which underlies all of the AtomDB spectra and has enabled the development of a module for estimating the sensitivity of emission lines and diagnostic line ratios to uncertainties in the underlying atomic data. We present these publicly available tools and results for several X-ray diagnostics of Fe L-shell ions and He-like ions as examples.

1. History of the AtomDB Project

The AtomDB project was originally designed to model X-ray emission from collisionally ionized, optically thin plasma in thermal equilibrium. It built on a history of development, starting with the Raymond–Smith code [1] in 1977 which created spectra for the most abundant 12 elements in a hot plasma using tabulated Gaunt factors and effective collision strengths. In 1985, this was updated by Brickhouse, Raymond, and Smith [2] to improve both the code and atomic data, particularly for the iron L- and K-shell ions. Effective collision strengths were updated and, significantly, the atomic data was separated from the code into independent data files.

In 2001, the first version of AtomDB was released [3]. This included the first release of the publicly accessible atomic database (APED: Astrophysical Plasma Emission Database), complete with defined filetypes as standardized FITS files [4]. In addition, the thermal plasma model, the Astrophysical Plasma Emission Code (APEC), was written and used to convert this database into line and continuum emissivities, which were and are still used extensively in spectral analysis tools such as XSPEC [5] and Sherpa [6]. This was accompanied by another large scale update to the atomic data in APED, in particular for modeling He-like and H-like ions.

In 2012, AtomDB version 2 was released [7]. This again included a large update to atomic data in APED: Ionization and recombination rates were updated, along with new excitation data for He-like and H-like ions. In addition, APED was expanded to include all the elements up to Nickel (previously it had covered only the 14 most abundant elements).

Four years later, AtomDB version 3 was released. This involved a major expansion of the code, the database and tools used to interpret the spectra to allow for the modeling of non-equilibrium ionization plasma, as well as updates to existing atomic data. In particular, there was a need to handle much larger data files (inner shell processes opened up many new levels), calculations for one ion at a time, and new input and output formats suitable for more flexible spectral analysis.

This large scale restructuring of the code exposed inflexibility in the APEC code, which had been designed to run from start to finish and produce one complete output file for the APEC model. To allow for future flexibility, the entire project was shifted to python with the release of the open source PyAtomDB (https://github.com/AtomDB/pyatomdb), which now underpins all aspects of the AtomDB project.

During and subsequent to this, we discovered that new data formats allowed for a significant improvement in AtomDB’s capabilities, enabling a wide range of additional models to be produced. In this paper, we will outline the PyAtomDB project, the new data formats created for non-equilibrium modeling, the new models which this has enabled for charge exchange (CX), non-Maxwellian electrons (Kappa), and the ability to perform uncertainty studies.

1.1. The Astrophysical Plasma Emission Database (APED)

The AtomDB database consists of a selection of data for every ion from H to Ni, all stored as binary FITS (http://fits.gsfc.nasa.gov/fits_standard.html) table files. The data are stored in several different files for each ion, separated by process. Throughout the database, all ions are indexed by their ion charge plus one, so for example He_1 is neutral Helium, and files related to it are stored in the directory $ATOMDB/APED/he/he_1/. There are up to eight different file types for each ion, as listed in

Table 1. File types for per-ion atomic data in the Astrophysical Plasma Emission Database (APED).

Type	Short Description	Details
IR	Ionization & recombination	Ionization from and recombination into the ion.
LV	Energy levels	Energy levels and photoionization cross section data.
LA	Lambda & A-values	Wavelengths and Einstein-A probabilities of transitions.
EC	Electron collisions	Electron-impact excitation rate coefficients or effective
		collision strengths.
PC	Proton collisions	Proton-impact excitation rate coefficients or effective
		collision strengths.
DR	Dielectronic recombination	DR Satellite line wavelengths and intensities.
PI	Photoionization	Photoionization cross sections for ions with data from XSTAR.
AI	Autoionization	Final state selective autoionization transition probabilities.

. In addition, APED contains several multi-element files, which hold data covering many ions together, as described in

Table 2. File types for global atomic data in APED. A full description of the data files and formats is included in Appendix A.

Type	Short Description	Details
abund	Elemental abundance	Abundance relative to solar for elements $Z\le 30$
ionbal	Ionization balance	Equilibrium ionization balance and non-equilibrium
		eigenvectors for each element
filemap	List of filenames	The data files used to create the latest release of the
		Astrophysical Plasma Emission Code (APEC) model

.

Data in the APED is stored as close to the original data format as possible. As a result, there is no single standard for values that are frequently represented as functions, such as collision strengths. Instead there are 10 or more different formats ranging from temperature/value pairs to fit coefficients, each of which can then be interpreted by the APEC code to produce the actual rate coefficient or other value required. This can be confusing for an end user, so part of the motivation for PyAtomDB was to provide a unified interface allowing users to request a piece of atomic data and get the actual rate or coefficient back, instead of the raw fit parameters, etc., which require further interpretation. This accessibility improvement makes using data in other models significantly easier.

1.2. The Astrophysical Plasma Emission Code (APEC)

APEC reads APED and processes the data using a collisional-radiative model to produce line and continuum emissivities on a set temperature and density grid. The plasma is assumed to be optically thin. For each ion at each grid node the excitation rates are calculated and combined with the spontaneous emission coefficients to form a collisional radiative matrix connecting the populations of all levels. Solving this gives level populations through simple multiplication and line emissivities per ion. Ionization and recombination into excited levels is also included in this matrix.

These line emissivities are then multiplied by the elemental abundance (by default [8]), and the ion population in equilibrium to produce the emissivity ${\epsilon }_{ij}$ for each line:

$${\epsilon }_{ij}=A_{ij}{N}_i=A_{ij}\frac{N_i}{N_{z1}}\frac{N_{z1}}{N_Z}\frac{N_Z}{N_H}{N}_H$$

(1)

where $A_{ij}$ is the spontaneous emission coefficient from level i to j, $N_i$ is the number density of ions in excited state i, $N_{z1}$ is the number density of the ion, $N_Z$ is the number density of the element, and $N_H$ is the number density of hydrogen. Therefore these terms from left to right are the spontaneous emission coefficient, the level population, the ion fraction, the elemental abundance relative to hydrogen, and the hydrogen density.

APEC then produces two output files for the plasma, storing two kinds of emissivity data. The line file contains the ${\epsilon }_{ij}$, wavelength, ion and transition identification for each line above a set ${\epsilon }_{ij}$ (by default, ${10}^{-20}$ ph cm³ s⁻¹). The coco files contain the continuum emissivities from bremsstrahlung, radiative recombination and two photon decay, in ph cm³ keV⁻¹ s⁻¹, compressed onto an emissivity vs energy grid by an algorithm [9] which ensures that when the table is linearly interpolated the result is within 1% of the original emissivity value at every point. This is stored on an element by element basis, assuming equilibrium ionization fractions. The coco file also contains the sum of the emission from weak lines not contained in the line file, known as the pseudocontinuum, which is compressed in the same way as the continuum. Within these two files, both of which are standard FITS files, the first Header/Data Unit (HDU) is the list of temperatures, then each subsequent HDU contains the line or continuum emission for one of these temperatures. This data can then be rapidly read by fitting codes such as XSPEC and used to analyze data, interpolating between neighboring temperatures in $log\left(T\right)-log\left(\epsilon \right)$ space.

With AtomDB 3, the introduction of the non-equilibrium ionization (NEI) model means that the ionization fraction, $N_{z1}/NZ$, could not be assumed to be fixed. There are therefore new nei versions of the same data files that do not have this factor included in their emissivity, and for each line and continuum component specify the driving ion which created the feature, as opposed to just the element or ion of the line. For example, the forbidden line of He-like O can originate from the excitation of O⁶⁺, ionization of O⁵⁺, or recombination of O⁷⁺. Creating a model spectrum then involves calculating the ionization fraction for the plasma in question and then multiplying the tabulated emissivities by this to the result to obtain the final emissivity for the model.

1.3. PyAtomDB

The introduction of AtomDB 3 created a vastly larger database to handle all of the underlying atomic processes involved in non-equilibrium processes. In particular, inner-shell excitation and ionization, along with excitation-autoionization, were included in the database for the first time and it was necessary to include ions that are not significant emitters in the X-ray region in equilibrium plasmas, such as M-shell ions.

Initially, APEC was modified to handle these data sets, although it became clear that this was not an efficient method for handling this. The original code started by loading all of the atomic data in the database, then iterating over each temperature, density, element, and ion to produce the spectra. Due to the code’s structure, changing the memory storage alone (required due to the entire database no longer fitting in a reasonable machine’s memory) required a complete re-write of the code. While doing this, it was also converted to Python due to the ease of development, support across different systems, widespread existing use in astronomy, and existing relevant packages such as Astropy [10,11] which PyAtomDB now uses. This conversion enables several long term goals of the project: Users can now run the code themselves, direct access to the APED is possible for those wishing to integrate atomic data lookups into their code, and spectral models based on APED can be developed by users and integrated easily into analysis tools. The resulting six different modules of the PyAtomDB package are listed in

Table 3. Modules within PyAtomDB and their purpose. A full description is provided in the code documentation (https://atomdb.readthedocs.io).

Module	Purpose
atomdb	Accessing the database, returning atomic data
apec	Thermal plasma model
spectrum	Creating spectra from apec outputs
atomic	Basic atomic data (element names, etc.)
const	Physical and code related constants
util	Utilities and helper routines

.

A full manual of the capabilities of PyAtomDB is beyond the scope of this paper. We will instead highlight some of the capabilities of the spectrum module, creating spectra of interest to end users, and extensions to it for modeling different, non-standard plasma types and investigating the effects of uncertainties on atomic data.

2. Spectral Modeling with PyAtomDB

The spectrum module in PyAtomDB contains routines for creating collisional ionization equilibrium (CIE) spectra of thermal plasmas, complete with applying instrument responses. These tools are also provided with wrappers, allowing them to be used in PyXspec (https://heasarc.gsfc.nasa.gov/xanadu/xspec/python/html/).

2.1. Non-Equilibrium Ionization Modeling

In addition to simple CIE plasmas, there are extensive models for NEI plasmas within PyAtomDB. These models represent a plasma where the electron temperature has rapidly changed, however the charge state distribution of the ions has not yet had time to reach a new equilibrium at this temperature. As the timescale for ionization and recombination is significantly longer than that for individual level populations to change, there is significant emission from ions at temperatures where they would not normally be observed in CIE. Such plasmas can be found in any place with shock heating—for example young Supernova Remnants (SNR) where the reverse shock rapidly heats the electrons leaving the plasma in an ionizing state until the plasma reaches equilibrium again e.g., [12]. Similarly, a recombining plasma, where the electrons have undergone a rapid cooling, can be found in many Mixed Morphology SNR, where rapid expansion leads to rapid cooling of the electrons and a subsequent recombination dominated plasma [13].

Figure 1 shows an equilibrium, ionizing from $k{T}_e=0.07$ keV and recombining from $k{T}_e=7.0$ keV spectrum for a $k{T}_e=0.7$ keV plasma at an ionization timescale of ${10}^{10}$ cm⁻³s. The ionizing case is dominated by He-like lines of Si and Ne, while the recombining case shows significant emission from He- and H-like Fe, as well as the H-like Si and Mg ions.

This model is identical to the existing rnei model in XSPEC, but has additional flexibility such as allowing the user to identify strong lines in a selected wavelength region of a non-equilibrium plasma. There is a similar model for a plane-parallel shock (pshock) [12].

2.2. Non-Maxwellian Modeling

Non-Maxwell–Boltzmann electron energy distributions can occur in a range of astrophysical plasmas such as the Earth’s magnetosphere [14], supernova shock waves (e.g., [15]), or solar flares (e.g., [16]). In these cases, the electron energy distribution is not well represented by a Maxwellian distribution due to an excess of high energy electrons. These can be modeled by a $\kappa$ distribution:

$$f_{\kappa }\left(E;\kappa ,T\right)=A_{\kappa }\frac{2}{\sqrt{\pi }}{\left(\frac{1}{k_{B}T}\right)}^{\frac{3}{2}}\sqrt{E}{\left[1+\frac{E}{\left(\kappa -\frac{3}{2}\right)k_{B}T}\right]}^{-\left(\kappa +1\right)}$$

(2)

where:

$$A_{\kappa }=\frac{\Gamma (\kappa +1)}{\Gamma \left(\kappa -\frac{1}{2}\right){\left(\kappa -\frac{3}{2}\right)}^{\frac{3}{2}}}.$$

(3)

When $\kappa$ approaches its lower limit of 3/2, the distribution is highly non-Maxwellian, when $\kappa =\infty$ it is once again Maxwellian.

Most electron energy dependent atomic data is, however, collated assuming the electron energies follow a Maxwellian distribution, as it saves substantial disk space, speeds computation of emissivities, and is relevant to the majority of applications. In theory the parameters such as electron collision strengths could be recalculated from first principles for non-Maxwellian plasmas, however the underlying cross section data either no longer exists for the vast majority of ions or, if they are, are too unwieldy for practical use. As a result, representing the $\kappa$ distribution as a sum of Maxwellian plasmas presents an attractive compromise. We have implemented [17] the decomposition method of [18] to create spectra for non-Maxwellian plasma, wherein several different Maxwellian temperatures are added to closely approximate a true $\kappa$ energy distribution.

Initially this was to be performed as a complete rewrite of the underlying model, however with the flexibility of PyAtomDB and the new data formats, we are able to achieve this using a simple 2-step process: The effective ionization and recombination rates are calculated by appropriately summing the Maxwellian rates as shown in Figure 2, the charge state distribution is calculated, and then the non-equilibrium emissivities can be combined with these to provide a non-Maxwellian spectrum.

2.3. Charge Exchange

Charge exchange occurs when a donor ion or atom with one or more electron still attached interacts with a recombining ion, which is missing at least one electron, resulting in the transfer of the electron from the donor to the receiver:

$$R^{y+}+D^{z+}\to R^{(y-1)+}+D^{(z+1)+}$$

(4)

where R is the recombining ion, and D is the donor. In most situations where charge exchange is postulated (e.g., the solar wind), the donor is neutral hydrogen as the only element with a significant enough abundance for the interaction to be observable. Typically, the exchanged electron is captured into an excited state of the ion, and these are often significantly higher n and l shells than those populated by electron impact excitation. As the ion relaxes to its ground state, the line emission is therefore significantly different from that observed in a typical thermal plasma. In some cases this can allow for the direct identification of CX, but in many more cases the foreground CX emissions adds an additional nuisance component to the spectrum of interest.

We released version 1 of the AtomDB Charge Exchange (ACX) model in 2014. This used a combination of analytical nl shell capture distribution formulae with a cascade to ground calculated through the APED database to create spectra useful for identifying CX emission [19], and successfully used it to model the spectrum of the hot and cool solar wind components.

However, the model had several shortcomings. On the physics side, the CX process was modeled without any velocity dependence, with the user selecting from one of four analytical cross section models to see which worked best. On the implementation side, the model relied on a computationally expensive process combining the AtomDB database to produce the spectra and line lists of interest. This meant that as the APED database was updated in the future, the changes in line wavelength and A-values were not reflected in the ACX model.

To expand the model, for ACX version 2 theoretical charge exchange cross sections (as opposed to the analytical $nl$ distributions), including velocity dependent effects, were obtained from the Kronos database [20,21,22]. PyAtomDB was then used to create spectra for capture into each $nl$ shell, and subsequent cascade to the ground state. By treating each of these as independent spectra, as if they were separate ions in the NEI model, spectra can be rapidly created in two steps: First, the cross section for capture into each $nl$ shell is calculated, then the spectrum for each $nl$ shell is multiplied by this cross section and then summed:

$${\epsilon }_{CX}=\sum ^{nl}{\sigma }_{nl}\left(E\right)v{\epsilon }_{nl}$$

(5)

where ${\epsilon }_{CX}$ is the total CX emissivity, $nl$ is the shell captured into, E and v are the energy and velocity of the center of mass of the donor-recombining ion, and ${\epsilon }_{nl}$ is the emissivity for capture into that one specific $nl$.

2.4. Cross Section Data Selection

The Kronos database contains only nlS selective cross sections for collisions with a hydrogen donor and a H-like recombining ion, and also for the bare recombining ion of C, N, O, and Ne. For all other bare ions, and for He as the donor, data is only resolved to the n shell.

In the case of nlS resolved data, the total capture into all the states which match the nlS is split evenly. For n resolved data the l shell distribution is handled by one of the analytical distribution functions described for ACX version 1 (the user may choose which one). For all other ions, where there is no data in the Kronos database, the model falls back to the ACX 1 models, with the total cross section fixed at $3\times {10}^{-15}$ cm².

Within the Kronos database, there are often different data sets present for the same ion. As the authors of the database themselves advised, we choose the data which exists from the best available technique if no specific set is recommended. The techniques are ranked from quantum molecular orbital close coupling (QMOCC) [23], atomic orbital close coupling (AOCC) [24], classical trajectory monte carlo (CTMC) [25,26], to multi-channel Landau Zener (MCLZ) [27].

2.5. Spectral Generation

For each nl or nlS to which capture is possible given the cross section in the Kronos database, we calculate the spectrum by creating a radiative matrix assuming no further collisional excitation as the electron cascades to the ground state. The A-values are taken from the AtomDB database for all ions where they exist.

For some heavier ions, the n shell for capture ends up being very high, e.g., 16 for the donor Fe²⁵⁺. This is beyond the AtomDB holdings, which stop at $n=10$ for He- and H-like ions and at $n\le$ 7 for all other ions. For these cases we have topped up the existing AtomDB A-values and wavelengths using calculations from Autostructure [28]. Where possible, the wavelengths for these calculations have been shifted to match the NIST ASD [29] values. The new energy levels were matched to the existing AtomDB levels using energy-order, symmetry, and degeneracy matching. For the Li-, He-, and H- like ions these matches were relatively trivial, though in the middle of the Fe L shell ions this matching becomes more problematic. Fortunately, the inaccuracies introduced by poor wavelengths at high energies in these complex ions can be largely discounted due to the fact that the sheer number of lines prevents any one line from producing a strong spectral signature. By contrast, a simple ion such as O⁷⁺ has a straightforward series of lines with well established wavelengths, and therefore discrepancies would be more pronounced.

Constructing the model in PyAtomDB has two distinct advantages: It allows the model to be rapidly developed (as it reuses most of the tools for NEI plasma) and it enables users to deploy it using the same commands.

3. Estimating the Impact of Uncertainties with PyAtomDB

Varying Atomic Data

Estimating the magnitude of the uncertainties on rates for fundamental atomic processes is a difficult task with no obvious correct solution. The source data comes from a wide range of theoretical and experimental sources, often lacking sufficient information on the associated uncertainties to make a meaningful estimate. Nevertheless, there is significant effort being applied to estimating atomic data uncertainties: Most atomic data publications in the past decade have some effort to make reasonable uncertainty assessments, and there are several groups e.g., [30,31] that are working on trying to estimate uncertainties on the existing atomic data.

Given the magnitude of this task, being able to direct the efforts by identifying the atomic data which line emission and ratio diagnostics are most sensitive to is essential for more accurate data analysis. We have developed the variableapec (https://github.com/AtomDB/variableapec) module, which makes use of routines in PyAtomDB to identify which atomic processes most strongly affect the resulting spectra. The user can vary any fundamental atomic data in the AtomDB database by a set amount and re-run emissivity calculations to identify which lines are sensitive to the changed parameter. This allows users to investigate the sensitivity of line and line ratio emissivities over a range of temperatures and densities. The module also checks for other lines affected from perturbing a single line.

Varying an atomic rate of a single transition may have cascading effects for other transitions, depending on the strength of the line varied or the levels of the transition. The user can input a range of wavelengths and cutoff threshold for the fractional emissivity change and plot lines sensitive to the transitions changed. This is an important diagnostic of lines that are sensitive to cascade effects as well as a diagnostic of correlated errors. Throughout this work we provide errors as the mean value of the positive and negative uncertainties for the fractional emissivity change. Line ratio errors are half the range of the maximum and minimum ratios obtained when applying both positive and negative uncertainties to the individual atomic data.

Figure 3 shows Fe xvii lines sensitive in the range of 10–20 Å to a ±30% change in A value in key emission lines at low density and the temperature where Fe xvii peaks in abundance. Varying the A value for the 3E transition alone resulted in a greater than 0.01% change in emissivity for not only the 3E transition but also four other lines in the range of 10–20 Å. Note that at a low density, the emissivities of the lines are almost only sensitive to the rate of population into the level (i.e., the effective collision strengths) as opposed to the rate out (i.e., the A values). 3E also has the greatest fractional change in emissivity of the six key transitions with A values perturbed by 30%. The value of 30% chosen for the uncertainty on the A-value is deliberately large for demonstration purposes, although it is in line with the uncertainties on the 3C and 3D line ratios reported in Bernitt et al. [32], and NIST Atomic Spectral Database uncertainties for these lines are class C ($\pm 18\%$) and D ($\pm 50\%$).

The new variability module also allows the user to vary the atomic rates of two separate transitions and recalculate diagnostics for that line ratio, such as the Fe xvii 3C/3D ratio shown in Figure 4. Line ratios can be sensitive to either temperature or density, making them useful diagnostics of estimating a plasma’s temperature or density. It is more informative to add uncertainties to two lines individually and recalculate line ratios rather than adding a blanket uncertainty to the final line ratio.

variableapec can also apply uncertainties for more complicated line ratios, such as the G and R ratios in the $n=2\to 1$ transition lines of He-like ions [33]. These are the ratios of the forbidden plus intercombination lines to the resonance lines (G) and the forbidden to intercombination (R) lines, which are temperature and density dependent respectively. It adds the same uncertainty to the forbidden, intercombination, and resonance lines and then calculates the error propagation on the individual rate changes to produce an error on the final line ratio. Figure 5 and Figure 6 show the effect of perturbing these key emission lines for the O vii R ratio and Fe xxv G ratio respectively. We applied a constant error over a temperature or density range, however the errors on the underlying data can be significantly temperature dependent, e.g., effective collision strengths. Future steps include recalculating line emissivities and line ratios with representative temperature-dependent errors.

The tool allows the user to recalculate line diagnostics at multiple temperatures or densities. It is informative to estimate the fractional error on the resulting line ratio after varying the A value or direct excitation rate of a density- or temperature-sensitive line. The O vii R ratio at three different plasma temperatures in Figure 5 shows that the magnitude of fractional error varies for different plasmas.

When a pair of lines is blended, it is more useful to add the uncertainty onto the three lines individually and recalculate the blended line ratio. Figure 7 shows the blended line ratio of (209.62 + 209.92) Å/213.77 Å for Fe xiii, a density-sensitive line ratio often used in solar physics [34].

4. Summary

The new PyAtomDB package had been released as open source code. It is now possible for the astronomy community to obtain all of the data, thermal plasma codes, and spectral models based on the project, and to develop extensions based on it.

We have developed modules for charge exchange and non-Maxwellian electrons, which are compatible with PyAtomDB and analysis tools such as PyXspec. The code can easily be extended to include new models the community may contribute.

In addition, we have used the newly available codes to develop a tool for estimating the sensitivity of spectral diagnostics to uncertainties in the underlying atomic data, which is immediately applicable to a range of plasmas for upcoming high resolution missions such as XRISM.