Add two-photon continuum

[jwreep]

Fixes #43

Not functional yet -- just introducing a placeholder for now.

A few questions have popped up:

1. We need the rest wavelength of the 2S1/2 (hydrogenic) and 1s2s S0 (helium-like) states. Should we use the theoretical or observed wavelength? (i.e. self._elvlc['E_th'] or self._elvlc['E_obs'])

2. Where and why does A_sum enter the calculation? It's not in the Equations in Chianti papers so far as I can see.

3. Should the output units be u.erg * u.cm**3 / u.s / u.angstrom as with f-f and f-b?

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 92.59%. Comparing base (93a85d8) to head (05b47fe).

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #260      +/-   ##
==========================================
+ Coverage   92.32%   92.59%   +0.27%     
==========================================
  Files          38       38              
  Lines        2788     2863      +75     
==========================================
+ Hits         2574     2651      +77     
+ Misses        214      212       -2     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[wtbarnes]

1. We need the rest wavelength of the 2S1/2 (hydrogenic) and 1s2s S0 (helium-like) states. Should we use the theoretical or observed wavelength? (i.e. self._elvlc['E_th'] or self._elvlc['E_obs'])

In the case of the FB continuum, I used the theoretical energies to fill in where the observed energies were missing. I would assume we should do something similar here?

fiasco/fiasco/ions.py

Lines 1385 to 1390 in 65130bb

	# Fill in observed energies with theoretical energies
	E_obs = self._fblvl['E_obs']*const.h*const.c
	E_th = self._fblvl['E_th']*const.h*const.c
	level_fb = self._fblvl['level']
	use_theoretical = np.logical_and(E_obs==0*u.erg, level_fb!=1)
	E_fb = np.where(use_theoretical, E_th, E_obs)

2. Where and why does A_sum enter the calculation? It's not in the Equations in Chianti papers so far as I can see.

I'm not sure about the "why" but see my comment on #43 for the "where"

3. Should the output units be u.erg * u.cm**3 / u.s / u.angstrom as with f-f and f-b?

That seems sensible as you often want to add the continuum contributions together

[jwreep]

Three outstanding issues (I think):
-It works for single values of temperatures/densities/wavelengths, but needs to be extended to arrays.
-It needs to be validated against Chianti (IDL and/or ChiantiPy)
-Need to add descriptions of the new functions

>>> import fiasco
>>> import astropy.units as u
>>> iron = fiasco.Element('Fe', temperature=1e7*u.K)
>>> iron.two_photon(10*u.angstrom, 1e10*u.cm**-3)
WARNING: No proton data available for Fe 25. Not including proton excitation and de-excitation in level populations calculation. [fiasco.ions]
WARNING: No proton data available for Fe 26. Not including proton excitation and de-excitation in level populations calculation. [fiasco.ions]
<Quantity [[3.28559482e-25]] cm3 erg / (Angstrom s)>

[jwreep]

Array handling added. Still need to validate against Chianti.

@wtbarnes, this will return the two-photon as a function of (wavelength, temperature, electron_density). This is an extra dimension over the f-f and f-b calculations, which hopefully isn't an issue?

I've also changed the logic from Chianti slightly to ensure that this function returns 0 for wavelengths < rest_wavelength. In Chianti, it only calculates the intensity for wavelengths > rest_wavelength. This maintains the size of the input array. The rest wavelength in the example below is 1.78 A, for example.

(base) reep@atrcl174 fiasco % python
Python 3.11.7 (main, Dec 15 2023, 12:09:56) [Clang 14.0.6 ] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import fiasco
>>> import astropy.units as u
>>> T = [1e6,1e7,1e8]*u.K
>>> wvl = [1,3,5,10,25,100]*u.angstrom
>>> ne = [1e9,1e10]*u.cm**-3
>>> iron = fiasco.Element('Fe',temperature=T)
>>> iron.two_photon(wvl,ne)
WARNING: No proton data available for Fe 25. Not including proton excitation and de-excitation in level populations calculation. [fiasco.ions]
WARNING: No proton data available for Fe 26. Not including proton excitation and de-excitation in level populations calculation. [fiasco.ions]
<Quantity [[[0.00000000e+00, 0.00000000e+00],
            [0.00000000e+00, 0.00000000e+00],
            [0.00000000e+00, 0.00000000e+00]],

           [[0.00000000e+00, 0.00000000e+00],
            [1.01445614e-19, 1.01445611e-17],
            [2.32387744e-15, 2.32387740e-13]],

           [[0.00000000e+00, 0.00000000e+00],
            [3.64770442e-20, 3.64770431e-18],
            [8.34977178e-16, 8.34977164e-14]],

           [[0.00000000e+00, 0.00000000e+00],
            [7.31047937e-21, 7.31047915e-19],
            [1.67241168e-16, 1.67241165e-14]],

           [[0.00000000e+00, 0.00000000e+00],
            [6.45079443e-22, 6.45079424e-20],
            [1.47745478e-17, 1.47745476e-15]],

           [[0.00000000e+00, 0.00000000e+00],
            [9.67610084e-24, 9.67610056e-22],
            [2.23657453e-19, 2.23657449e-17]]] cm3 erg / (Angstrom s)>
>>> iron.two_photon(wvl,ne).shape
WARNING: No proton data available for Fe 25. Not including proton excitation and de-excitation in level populations calculation. [fiasco.ions]
WARNING: No proton data available for Fe 26. Not including proton excitation and de-excitation in level populations calculation. [fiasco.ions]
(6, 3, 2)

[wtbarnes]

@jwreep Thanks! Sorry for the unresponsiveness the past few days. This is looking great. Regarding the difference in shape due to density, I don't think the shape mismatch is an issue. In general, to add along temperature and wavelength with the ff and fb, you have to assume something about the density for the two photon anyway which here you can do be either a) assuming a single density or b) assuming a constant pressure by using the couple_density_to_temperature keyword (which actually will work for any temperature-density dependency). We just need to make sure that keyword argument can propagate through from calling the two-photon continuum to the level populations call.

Regarding the IDL comparisons, I've developed a way to systematically test against the IDL results and include it as part of the test suite. You can see an example of how this works in fiasco/tests/idl if you want to take a crack at it or I can as well. As part of #259, I'm improving the docs on writing these kinds of tests so we could wait on merging that (nearly done) and then write the IDL comparison tests.

[jwreep]

I'm having trouble getting hissw installed on Windows, and trouble getting an IDL license server to work on Mac . . . so it might be a bit of time before I can do the comparisons until I sort out other issues.

[jwreep]

The function needs error handling for ions missing data, which will cause issues when trying to calculate the total two-photon emission from all ions.

>>> lith = fiasco.Element('Li', temperature=1e7*u.K)
>>> lith[2]._elvlc
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/reep/Documents/Forks/fiasco/fiasco/base.py", line 157, in property_template
    return DataIndexer.create_indexer(self.hdf5_dbase_root, data_path)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/reep/Documents/Forks/fiasco/fiasco/io/datalayer.py", line 32, in create_indexer
    return DataIndexerHDF5.create_indexer(*args)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/reep/Documents/Forks/fiasco/fiasco/io/datalayer.py", line 71, in create_indexer
    raise KeyError(f'{top_level_path} not found in {hdf5_path}')
KeyError: 'li/li_3/elvlc not found in /Users/reep/.fiasco/chianti_dbase.h5'

[wtbarnes]

Thankfully, this is exactly what the needs_dataset decorator is for. If the listed datasets are not available for that ion, then a MissingDatasetException is raised before doing anything in the function. Have a look at some of the other methods on Ion to see how to use this decorator.

[jwreep]

Clever design! That fixes the issue.

[jwreep]

The attached plot shows a comparison for Fe XXVI at 10^7 K. (Wavelength vs intensity)

The Python version as coded currently is missing a factor of 1/wavelength. I don't know how the units could then work out . . .

It's also not clear to me why it's a factor of 10^28 off from IDL. I don't see the assumption about EM in the IDL code.

[jwreep]

I suspect the equation in Young et al and the CHIANTI user's guide is not correct. The IDL coding (and that in ChiantiPy) does not match that equation.

I had assumed the code used a column EM, but it should be a volume EM. A_sum still appears to be unitless, and I somehow missed a factor of 1/n_e.

[wtbarnes]

The Python version as coded currently is missing a factor of 1/wavelength. I don't know how the units could then work out

Yeah, I'm confused as well. Where does that 1/Å factor come from? Unless the units of the spectral distribution function are 1/Å? I've looked at the Goldman and Drake papers and I cannot find anything that obviously indicates that. In the IDL code, I'm also confused as to how to reconcile that with what is in Landi et al. (1999), but maybe there's not point as what is implemented now in the IDL routines is clearly from Young et al. (2003).

It's also not clear to me why it's a factor of 10^28 off from IDL. I don't see the assumption about EM in the IDL code.

I noticed looking through the IDL code that, like the ff and fb methods, there is a 1e40 scaling that is being applied. Do you account for that? Could it be that there is an errant 1e12 missing that combined with 1e40 gives you your magical 1e28?

A_sum still appears to be unitless,

This is still worrying to me. Looking in the Parpia and Johnson paper, I cannot find any value that looks like it corresponds to A_sum. Since that posting on the CHIANTI mailing list isn't getting much traction. Maybe we should post an issue to ChiantiPy directly since it implements these in the exact same way.

I had assumed the code used a column EM, but it should be a volume EM.

Why? In the documentation for the IDL routines, it seems to be in units of column emission measure and I cannot see any obvious place where there is a conversion to volume EM.

In general, when it comes to implementing this in fiasco, I'd prefer to not even bother with providing the EM as an input. As long as the outputs are not integrated over temperature, I'd prefer to just allow the user to deal with the EM themselves as it reduces the number of inputs and simplifies the logic of the code. Currently in fiasco, the intensity method is the only place where EM is an input and that's because there is an integration over temperature.

[jwreep]

I'm trying to follow their implementation as closely as possible, and they have an extra factor of 1/wavelength that I missed (along with 1/n_e). There is an emission measure built into CHIANTI (with a value of 1, so scale up as you wish), and from dimensional analysis I assumed it had to be a column EM. Fixing my mistakes (1/(n_e * wvl) ~ 1/cm^2), it should actually be a volume EM. The latest commit should have the desired units (erg cm^3 s^-1 A^-1).

Yes, the factor of 1e40 is present in their two_photon routine, so that's part of the huge missing factor. Part of it was the missing 1/n_e, but I don't think that fixes it all.

Where does that 1/Å factor come from? Unless the units of the spectral distribution function are 1/Å?

As to the exact equation, I'm still a bit stumped, but the Gronenschild & Mewe paper I linked in #43 does have the extra factor of 1/wavelength (Equation 12 - rest wavelength/wavelength^3). They've clearly manipulated this equation or a similar one somehow, but I can't find any papers that match the form they use, and can't figure out exactly what they've done. If you look closely at the IDL code, it does actually have rest wavelength/wavelength^3 in the output:

 y=wvl0/wvl[w_ind]
distr=y * spl_interp(y0,psi0[*,iz-1],y2,y)/wvl[w_ind]
distr1=rescale*factor*1d8*avalue[iz-1]*this_abund* $
                       (distr/**wvl[w_ind]**) * $
                       (ioneq1[i]*pop[pop_idx]/edens[t_ind[i]]) * dem_arr[t_ind[i]]

[jwreep]

Ah, I think I should clarify: the function is per unit EM (erg cm^+3 s^-1 A^-1), so the output of this function can still be multiplied by a column EM to give intensity units, just like f-f and f-b.

[wtbarnes]

Ah, I think I should clarify: the function is per unit EM (erg cm^+3 s^-1 A^-1), so the output of this function can still be multiplied by a column EM to give intensity units, just like f-f and f-b.

Right, I think that's the best way to think about it. Effectively then, the EM doesn't have to be part of the function at all.

[jwreep]

For completeness, the equation, as far as I can tell from the coding in IDL, is:

and A_sum is just 1 for helium-like.

I'll post an issue to ChiantiPy for clarification.

[jwreep]

Fe XXVI (edit: it should be the hydrogenic Fe XXVI, the plots are mislabeled) looks alright, but there are issues with other ions. Not clear at the moment what the issue is since I haven't broken it down by ion yet. The factor of 1e30 is likely 1e40 / n_e, since I'm assuming a density of 10^10 cm^-3 here.

Other temperatures have similar issues.

@wtbarnes, does the IDL version have proton excitation/de-excitation data for any ions? When I calculate this for everything, I get these warnings for all of the ions:

WARNING: No proton data available for O 7. Not including proton excitation and de-excitation in level populations calculation. [fiasco.ions]
WARNING: No proton data available for O 8. Not including proton excitation and de-excitation in level populations calculation. [fiasco.ions]

[jwreep]

Carbon is the (primary?) culprit:

In my previous plot, the bump is He II (304!), but I don't think helium is a problem. I think its contribution to the total is just getting swamped by carbon.

I have no idea why carbon is an issue. Have the ionization equilibria and level populations been thoroughly checked against the IDL? If I had to guess, I would presume the ionization balance between C V and VI are off.

[jwreep]

Might need to recheck this floating point math. It might be a deeper issue. I've added tests to fill in the CodeCov gaps (He-like ion and one with a missing data set). The comparison for C V is failing, and gives a different value from what I get when calculating by hand. The different OSes/Python versions are giving slightly different answers as well.

This is frustrating . . .

[wtbarnes]

Sorry for being so unresponsive for several weeks!

The comparison for C V is failing, and gives a different value from what I get when calculating by hand. The different OSes/Python versions are giving slightly different answers as well.

Regarding the differences in OS/Python versions, having a cursory look through the CI, I see that the test that is failing for every OS/Python version is generating an output of 4.020540e-25 consistent out to those digits with some differences out past the sixth decimal place while the expected answer in the test is approximately 4.04e-25. The relatively small differences between the different CI runs is not as concerning to me as the difference between the expected value and what the tests on the CI are returning for all OS/Python versions.

These small differences between versions could be down to different versions of numpy and could be a result of solving the level populations the way I am. Especially for small level populations, there can be relatively large fluctuations due to numerical stability.

I also noticed that some of the tests are issuing warnings about missing proton data. Should the comparisons be run with or without the protons? My understanding of the above discussion was that we wanted to exclude the proton rates by default in the two-photon calculation.

but I'd like to request that you perhaps could write the IDL comparison test since I'm still a bit shaky about using hissw and whether it's working properly with CHIANTI for me.

Yes, I'd be happy to do that! Once we get the above issues ironed out, I can either push directly to this PR, or we can open up a new PR to do it.

Thanks for your persistence. This has been much more of a slog than I anticipated though to put it in perspective, it took me years (!) to sort out what the issues were with the free-free and free-bound continua!

[jwreep]

Thanks, Will! I assumed you were busy with eclipse meetings.

Yes, this is definitely a lot more of a slog than I wanted it to be. Nearing the end point, though, I hope!

while the expected answer in the test is approximately 4.04e-25.

Yep, this is what I'm getting when doing it by hand:

>>> temperature = np.logspace(5, 8, 100)*u.K
>>> c5 = fiasco.Ion('C V', temperature,hdf5_dbase_root=fiasco.defaults['hdf5_dbase_root'])
>>> c5_emission = c5.two_photon(200 * u.Angstrom, electron_density = 1e10* u.cm**(-3))
>>> c5_emission[0,30,0]
<Quantity 4.04634243e-25 cm3 erg / (Angstrom s)>

Is there anything in the smaller test database versus full database that might cause the difference?

Edit: yes! I rebuilt with the test database and get a different result.

>>> files = get_test_file_list()
>>> build_hdf5_dbase(fiasco.defaults['ascii_dbase_root'], fiasco.defaults['hdf5_dbase_root'],files=files)
>>> temperature = np.logspace(5, 8, 100)*u.K
>>> c5 = fiasco.Ion('C V', temperature,hdf5_dbase_root=fiasco.defaults['hdf5_dbase_root'])
>>> c5_emission = c5.two_photon(200 * u.Angstrom, electron_density = 1e10* u.cm**(-3))
>>> c5_emission[0,30,0]
<Quantity 4.02054043e-25 cm3 erg / (Angstrom s)>

I don't know whether that's a bug or expected, though!

There's no difference for C VI:

>>> c6 = fiasco.Ion('C VI', temperature,hdf5_dbase_root=fiasco.defaults['hdf5_dbase_root'])
>>> c6_emission = c6.two_photon(200 * u.Angstrom, electron_density = 1e10* u.cm**(-3))
>>> c6_emission[0,30,0]
<Quantity 8.25316887e-26 cm3 erg / (Angstrom s)>

which is the same as with the full database. The test in test_collections (H I + He II) also gives the same result.

[wtbarnes]

Well this is mysterious. And now the test for the full database fails (as expected given your earlier findings)! https://github.com/wtbarnes/fiasco/actions/runs/8655635096/job/23735165230?pr=260#step:10:178

I'm confused as to why this value should be dependent on the presence of other files in the database. My first thought was the proton rates, but C V has no psplups files so it can't be that.

[wtbarnes]

Ok I think I figured it out. The reason was due to the absence of the c_5.reclvl file in the test database. In the level populations calculation, there is optional correction of level resolved ionization and recombination that is applied if the reclvl or cilvl files are present. If it is missing, it is skipped as it is assumed there isn't enough information there to calculate that correction.

I've just pushed a commit that adds this file to the test database and I reverted your test back to the old value. The test and full database tests should now pass 🤞.

[jwreep]

Thanks! Would've taken me a week to figure that one out.

[wtbarnes]

It is extremely subtle. Also, its probably my fault for not adding a warning if the code tries to calculate the correction, but can't. The reason its excluded is because most ions don't have these files so you'd be getting warnings all the time unless you set the keyword argument to false.

[wtbarnes]

Ha well now the v9 test is failing because C V has no reclvl file in that version since the level-resolved ionization and recombination is treated differently.

I think the right way to resolve that is to just skip the two-photon test for v9, at least until we fully support v9 anyway and then figure out a better solution then. To do this, you can use the decorator as is done here:

fiasco/fiasco/tests/test_ion.py

Line 245 in 93a85d8

@pytest.mark.requires_dbase_version('<=','8.0.7')

[wtbarnes]

Now that the tests are all passing, I'll do a review of the code later today and we can hopefully get this merged in the next day or so. I'm now inclined to just leave the IDL test to a separate PR.

Thanks again for your persistence!

[wtbarnes]

@jwreep I've just pushed a commit that does some minor refactoring of the logic in the 2p method on Ion and IonCollection. Have a look at it and let me know what you think.

If the tests pass here, I'm happy to merge. We'll need to open up an issue to track the IDL comparisons for this method first.

[jwreep]

All of these changes look reasonable and looks like cleaner code than mine. I think we should also open an issue to check / verify whether dielectronic ions impact this calculation, as discussed above. If so, they need to be properly incorporated.

[wtbarnes]

Everything looks good to me and the tests are all passing and the docs look good! The only remaining bit is the clarification on the intensity comment.

[wtbarnes]

Thanks for this massive effort @jwreep!