Estimating line ratios

[fsciortino]

This is a short demo aimed at answering questions by @Didou09, related to the issue/post in #45.

Here we demonstrate how to get line ratios, particularly focusing on the k and w ratios of the K-alpha Ar spectrum. We will show here how to reproduce Fig. 8 in J E Rice et al 2015 J. Phys. B: At. Mol. Opt. Phys. 48 144013 (https://iopscience.iop.org/article/10.1088/0953-4075/48/14/144013) using theoretical Photon Emissivity Coefficients (PECs). We will make use of PECs that are not distributed by ADAS (taken originally from atomDB), but that can be requested by emailing me. The same Aurora functionality can be used for any spectra for which ADAS data is available. Important: it is up to the user to assess the atomic data quality. Ask an expert for help if you're in doubt.

We first identify the k and w lines from the PEC files. We will only consider the dominant processes that populate upper levels of these lines: dielectronic recombination from the Li-like state of Ar for the k line, and excitation from the He-like state of Ar for the w line. Looking through the ADF15 (PEC) files mentioned above, we can identify the 2 lines of interest as the following

lam_Ar_k = 3.990000 # A, 1s 2p^2 {}^2D_{3/2} - 1s^2 2p {}^2P_{1/2} -- ISEL=1107 (DRSAT)
lam_Ar_w = 3.949075 # A, 1s 2p {}^1P_1 - 1s^2 {}^1S_0 -- ISEL=30 (EXCIT)

We next load all the PECs in the chosen files

# the k-like is only populated via dielectronic recombination from Li-like Ar, the w line mostly by excitation from He-like Ar
# PEC files also contain other minor contributions, which could be summed for completeness
filepath_he = 'pec#ar16.dat'
filepath_li = 'pec#ar15.dat'
log10pec_dict_he = aurora.read_adf15(filepath_he)
log10pec_dict_li = aurora.read_adf15(filepath_li)

After looking either at the ADF15 files or the Aurora log10pec_dict_* keys, you can plots the rates for excitation and recombination components using the plot_lines argument. The following, for example, will show that the Li-like 3.99A line is much stronger than the 3.9877A one:

log10pec_dict_li = aurora.read_adf15(filepath_li, plot_lines=[3.99,3.9877])

To look at line ratios, one can either identify specific lines of interest (typically, looking at the spectral notation at the bottom of ADF15 files), or add all lines within a central wavelength range near the line of interest. If most lines except the one of interest are actually very weak, then only the strong one will matter.

Here, we plot only the ratio of the 2 lines identified above. Let's load the interpolanting functions specifically for the 2 lines of interest:

log10pec_k_interp = aurora.read_adf15(filepath_li)[lam_Ar_k]['drsat']
log10pec_w_interp = aurora.read_adf15(filepath_he)[lam_Ar_w]['excit']

We then create Te grid of interest, to reproduce the Te range of Rice et al. J. Phys. B: At. Mol. Opt. Phys. 48 (2015) 144013:

Te_grid = np.linspace(1e3, 3e3, 1000) 

and we evaluate the intensity of the 2 lines at a single value of ne:

ne_cm3 = 1e14
k_intens = log10pec_k_interp(np.log10(ne_cm3), np.log10(Te_grid))[0,:]
w_intens = log10pec_w_interp(np.log10(ne_cm3), np.log10(Te_grid))[0,:]

Finally, we can plot the results and compare to Fig. 8 in the aforementioned paper:

fig,ax = plt.subplots()
ax.plot(Te_grid/1e3, 10**k_intens/10**w_intens)
ax.set_xlabel(r'$T_e$ [keV]')
ax.set_ylabel(r'k/w emissivity ratio')
plt.tight_layout()

This will give the following plot:

[Didou09]

There is an issue with this plot:
it is inaccurate by ~3 orders of magnitude, that's because we are not using the k line intensity but rather the intensity of another line with identical wavelength also found in pec#ar15.dat, as explained in a dedicated issue #48.

[fsciortino]

@Didou09 this has now been addressed in commit 9c87d01, which is already in the master branch. The (updated) documentation should explain well the change, but let me clarify it here as well. You can choose 2 options to deal with the issue that you pointed out:

1. Just run as you were doing, and note that the second line appearing in the ADF15 file with the same wavelength will have a 1e-6 A added to it. This is not ideal/clear, but at least allows dictionary elements not to be overwritten.

2. Add the keyword argument index_lines=True to each call of read_adf15. This will then load the log10pec_dict dictionary with keys that are equal to the ADAS "ISEL" indices that you can find in files. Note that there will be different indices for different populating processing, i.e. there may be an index for the excitation component, another index for the recombination component, etc.. The option index_lines=False, identifying lines by their wavelengths, is probably clearer/easier for most users, but can run into issues in some cases, as you pointed out.

I have not changed any dependencies on aurora.read_adf15. Let me know if you find anything that doesn't work as expected! Note that the change is in the master branch, but not yet on PyPi or conda.

[fsciortino]

Quick update for future readers, following PR #51 , which has now changed the way ADF15 (PEC) files are processed, now using pandas DataFrames: the following example reproduces the behavior discussed above with the new syntax.

import numpy as np
import matplotlib.pyplot as plt
plt.ion()
from scipy.interpolate import interp1d
import aurora

##### Plot k/w ratio over temperature at one value of density
filepath_he = 'pec#ar16.dat'
filepath_li = 'pec#ar15.dat'
isel_Ar_k = 1107 # 1s 2p^2 {}^2D_{3/2} - 1s^2 2p {}^2P_{1/2} -- ISEL=1107 (DRSAT)
isel_Ar_w = 30  # 1s 2p {}^1P_1 - 1s^2 {}^1S_0 -- ISEL=30 (EXCIT)

# the k-like is only populated via dielectronic recombination from Li-like Ar, the w line mostly by excitation from He-like Ar
# PEC files also contain other minor contributions, which could be summed for completeness
trs_he = aurora.read_adf15(filepath_he)
trs_li = aurora.read_adf15(filepath_li)

# select w and k lines of the Ar K-alpha spectrum
tr_Ar_k = trs_li.iloc[isel_Ar_k-1]  # "-1" because of python indexing (ISEL indices start at 1)
tr_Ar_w = trs_he.iloc[isel_Ar_w-1]

assert tr_Ar_k['lambda [A]'] == 3.990000 # A
assert tr_Ar_w['lambda [A]'] == 3.949075 # A

# select a single transition and plot it as a function of ne and Te:
tr = trs_he.iloc[isel_Ar_w]
aurora.plot_pec(tr)

# to look at line ratios, one can either identify specific lines of interest (typically, looking at
# the spectral notation at the bottom of ADF15 files), or add all lines
# within a central wavelength range near the line of interest. If most lines except the one of interest
# are actually very weak, then only the strong one will matter

# chosen density (units of :math:`cm^{-3}`)
ne_cm3 = 1e14

# create Te grid of interest (units of eV)
Te_grid = np.linspace(1e3, 3e3, 1000) # range of Rice et al. J. Phys. B: At. Mol. Opt. Phys. 48 (2015) 144013

k_intens = tr_Ar_k['log10 PEC fun'].ev(np.log10(ne_cm3), np.log10(Te_grid))
w_intens = tr_Ar_w['log10 PEC fun'].ev(np.log10(ne_cm3), np.log10(Te_grid))

# plot line ratio
fig,ax = plt.subplots()
ax.plot(Te_grid/1e3, 10**(k_intens-w_intens))
ax.set_xlabel(r'$T_e$ [keV]')
ax.set_ylabel(r'k/w emissivity ratio')
plt.tight_layout()