Use FAC data files for CR modeling

[cjperks7]

0) What's New

• Function under colradpy.convolve_EEDF.py to convolute given cross-section data (units [cm2]) with an electron energy distribution function (EEDF, units [1/eV]) to generate rate coefficients (units [cm3/s]). More details in the following sections:

  • 1. Cross-sections and their extrapolations
  • 2. EEDF convolution

• Function under colradpy.read_FAC.py and colradpy.read_FAC_utils.py to load data files for energy levels, transition probabilities, cross-sections, and rate coefficients produced by the Flexible Atomic Code (FAC) and prepare them for collisional-radiative modeling by the ColRadPy main class. More details in the following sections:

  • 3. Minimum working example of using FAC data
  • 4. Various sequence modeling validation (H-, He-, Ne-like + satellites)

• Note that FAC provides data for the following processes, all of which ColRadPy loads: - Energy levels - Spontaneous emission decay rates (up to all multipole orders) - 2-photon decay rates for H-like and 1S0 state of He-like - Collisional excitation cross-sections - Radiative recombination (and photoionization, unused here) cross-section - Dielectronic recombination strengths (units [cm2*eV]) as well as autoionization rates for doubly excited states - Collisional ionization cross-sections - Also know that the function pfac.fac.MawellRate() can convolve these cross-sections into rate coefficients and print to an ASCII file, but I was having trouble with non-physical results so preferred my own convolution within ColRadPy. Would see the rate go back up at high-Te so seems like a problem with numerical tricks employed

• FAC data files are assumed to have a common prefix and regular suffix tied to the process, i.e. Ar01a.en has common prefix Ar01 and suffix a.en denoting it's an ASCII-format energy level file

  • a.en is the ASCII-format energy level file
  • a.tr is the ASCII-format spontaneous emission rate file
  • a.ce, a.rr, a.ai, and a.ci are the expected suffixes for ASCII-format excitation, radiative recombination, autoionization/dielectronic recombination, and collisional ionization cross-section files
  • 2-photon emission rates are by an interpolation function included in FAC

• Note that I've also given the option to use Maxwell-average data files from pfac.fac.MaxwellRate

  • These then have the expected suffixes ce.mr, rr.mr, and ci.mr

• Note that I have also made some changes to colradpy.write_pecs_adf15() that I haven't committed to this branch

  • The changes include some quality-of-life conveniences and further documentation I like to have in ADF15 files
  • Do you want any of this?

1) Cross-sections and their extrapolation

• Within the cross-section data files produced by FAC, included are formulas for the asymptotic behavior of the cross-sections for collisional excitation, radiative recombination/photoionization, and collisional ionization - Note that the default simulation grid in FAC is from the transition energy to about 30x larger - From playing around with extending and refining that grid, I think it can be assumed you get non-physical results outside the default energy grid. Or am I wrong? - The below examples in this section will present cross-section from this refinement exercise, but note that cross-section data used in other sections are run with FAC in it's default grid

• Excitation for H-like Ar

  • The asymptotic behavior for collisional excitation cross-sections should follow the Bethe formula as illustrated in Eqns 1 and 2 in [R. Mewe, Astron&Astrop, vol 20 (1972)] for optically allowed and forbidden transitions respectively.
  • 
  • It can be seen that the FAC cross-section data converges to the Bethe formula at ~30x the transition energy before diverging. I assume that divergence is non-physical since it's flatlining.
  • In ColRadPy, it is assumed you ran FAC with appropriate energy grids, so as a function of the energy grid that the electron energy distribution function is calculated on, any energy points above the FAC simulation grid are filled with the Bethe formula. All lines labeled "ColRadPy" in this section were created using the default FAC grid
  • Doing this, I find that the trend looks good, but there's ~10% bump in magnitude as can be clearly seen in the gray and orange lines above. I think this is decent enough.

• Radiative recombination for H-like Ar

  • 
  • Here it can be seen that as E_electron->0, the cross-section strongly increases. As such, it's actually more important to include the low-energy behavior of the cross-section using the asymptotic formula provided in FAC
  • At the high-energy side, it can be seen that the trend of the formula doesn't make a major distinction between S and P states at fixed J even though they ultimately have different behavior. It gets S states really well
  • As such, it is assumed that for electron energies higher than the FAC simulation grid, the cross-section is zero. This can cause an under-estimation of the rate coefficients at high-Te

• Ionization for H-like Ar

  • 
  • As can be seen, the FAC data basically just follows the formula
  • The cross-section above the binding energy is set to zero. At energies above the FAC simulation grid, data is filled using the formula
  • In between FAC data points a linear interpolation is used which slightly underestimates the cross-section as can be seen

2) EEDF Convolution

• This function is a pretty simple integration given cross-section data

• Plotted is a validation exercise of this function

  • 
  • Top left is the inputted cross-section data as a function of the incident electron energy, top right is the resultant rate coefficients, bottom left is the velocity weight used, and bottom right is a span of electron energy distribution functions used

• Notice that I have found it important to include relativistic effects (denoted "rel") when doing this integration

  • This is because at Te~10keV you start to have a significant population of electrons with energy ~mc^2. As such, when using a Maxwellian distribution, the code automatically turns on relativistic effects for Te>10keV.
  • The strongest effect that I observed is in the velocity weight. Obviously you can't go faster than the speed of light
  • To a lesser extent, also your maximal entropy energy distribution is no longer Maxwell-Boltzmann, but Maxwell-Juttner
  • Without relativistic effects and using the same cross-section data as in [Mewe 1972], I get the same rate coefficient as the explicitly derived formula given in the paper.

• Right now, only Maxwellian distributions are supported. I have intentions for doing some non-Maxwellian CR modeling on WEST, so in another Issue I'll address adding maybe a data file. TBD.

• Note that this function also convolves dielectronic recombination strengths assuming the energy-dependence is just a delta function

3) Minimum working example of using FAC data

• The following is what's generally needed to run ColRadPy using FAC data

# Modules
import sys, os
from colradpy import colradpy

# Ion
sp = 'Ar'
nele = 1
Znuc = 18

# Common FAC atomic data files name
file = os.path.join(
    '/home/cjperks/2008_SPARCxray',
    'FAC/ATOM',
    sp,
    sp+'%02d'%(nele)
    )

# Simulation grids
temp_arr = np.logspace(np.log10(8.62), np.log10(8.5e4), 41) # [eV]
dens_arr = np.logspace(np.log10(1e10), np.log10(1e20), 21) # [cm^-3]
meta_arr = np.array([0])

# Load atomic data
crm = colradpy(
    file,
    meta_arr,
    temp_arr,
    dens_arr,
    use_recombination=True,
    use_recombination_three_body = False,
    use_ionization=True,
    suppliment_with_ecip=False,
    # New material
    atomic_data_type = 'FAC',       # Flag data source
    ele = sp,                       # Ion species
    nele = nele,                    # Number of electrons
    Zele = Znuc,                    # Nuclear charge
    EEDF = 'Maxwellian',            # Electron energy distriubtion function
    atomic_physics = 'incl_all',    # Atomic data files to search for
    )

# Runs collisional-radiative model
crm.solve_cr()

• The additions are back-compatible if you just don't include anything under # New material
• The file must be the path to the common prefix of all the FAC data files, i.e. /path/to/data/Ar01
• atomic_data_type default is adf04 for back-compatibility
  • Controls what read data file functions are used
• ele, nele, and Zele were found to be needed when constructing crm.data['atomic'] to sort out ionization stages stored in the energy level data file as well as some expected documentation from the ADF04 format
• EEDF is a flag for what can of electron energy distribution function is used
  • if EEDF = 'Maxwellian_mr' you're signalling that you've used the pfac.fac.MaxwellRate function so the code searches for files with the appropriate suffix, i.e. ce.mr
  • if EEDF = Maxwellian you're signalling that you want ColRadPy to Maxwell-average cross-section data so the code searches for files with the appropriate suffix, i.e. a.ce
• atomic_physics is a flag or list of the atomic physics data files to use
  • if not incl_all, needs to be a list that a subset of ['en', 'tr', 'ce', 'rr', 'ai', 'ci']
• suppliment_with_ecip does't work with FAC data right now since I set all zpla and zpla1 = -1 for lack of something better to do

4) Various sequence modeling validation (H-, He-, Ne-like + satellites)

• H-like Ar - I have been given high-quality PEC data files for Argon made by Adam Foster (Harvard) for Francesco Sciortino's (MIT) thesis work. From what Adam told me, he created the excitation data using FAC. Recombination and ionization were made by some guy named Nigel and he doesn't know what he did. So, to benchmark my work I should be able to reproduce his data. In the following legends, his data is marked as "AF". - Note that during this exercise, I used Gaussian, Voigt, and 2-photon emission profiles using the Aurora code as detailed in the following PR: https://github.com/fsciortino/Aurora/pull/88 - - - Per [Mewe 1972], we expect the 1s-2p PEC's to have 6% correction to the excitation rate from ground as we do find (all that work to basically get coronal equilibrium...) - Good agreement is found. At high-Te for the excitation rates, the discrepancy is due to including relativistic effects in the integration. At high-Te for the recombination rates, the discrepancy is because I didn't implement a high-energy extrapolation to the cross-section but the rate is relatively low enough this seems like an acceptable issue. I think the discrepancy in magnitude can be chalked up to numerics since which is greater swaps around with Te - Lines have been corrected to experimental wavelengths - Note no ionization PECs. To be addressed in Issue #17 - Note He-like satellites were not included in the above plot - Note that I don't think I get the doublet ratio exactly right as per "J.E. Rice et al, J. Phys. B: AMO, 48 (2015)". I think this is more so a problem with FAC underestimating for the 2S_1/2 -> 2P_1/2 collisional transitions and such. The results above are close enough so a different battle for another day.

• He-like Ar

  • 
  • 
  • In the second plot, you can see again the same discrepancies observed in the H-like Ar data. Mostly my magnitudes are slightly smaller which maybe be mostly numerics, but I get the overall temperature dependence the same
  • In the left panel of the first plot, you can see a synthetic spectrum for the He-like Ar K-alpha sequence comparing data I generated against Adam Foster's (AF). Note that I didn't correct my wavelengths to experimental values yet but I did for Adam's. Good agreement for the excitation and recombination components are found.
  • Voigt profiles were used for the w,x,y,z lines
  • On the right panel of the first, you can see a synthetic spectrum for the He-like satellites to the H-like Ar Lyman-alpha sequence. - The shape of my curve is in good agreement with that published in Figure 2 of "J.E. Rice et al, J. Phys. B: AMO, 48 (2015)" which was also calculated using FAC but also with the wavelengths experimentally adjusted which was about 1mA consist with my results. - The references within explain calculations of line factors for the satellites that I haven't compared against but should some day. - Recall Adam's recombination PECs were not calculated using FAC but some other source which could explain the observed discrepancy. I have checked that I've modeled all the same singly and doubly excited states (more actually) as Adam's, so it might not be because a process is missing - Anyway, the satellites illustrate that excitation-autoionization and dielectronic recombination is included and working

• Note that per how ColRadPy is written, I found it more convenient to add together the rates from radiative and dielectronic recombination due the preparation step so when I write "recomb" it includes both processes