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[bot]]

Codecov Report

All modified and coverable lines are covered by tests :white_check_mark:

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

Additional details and impacted files

```diff @@ 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 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: 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? https://github.com/wtbarnes/fiasco/blob/65130bb53b42a9d5606470d0ed85eccde87057c9/fiasco/ions.py#L1385-L1390

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.

[wtbarnes]

Hmm do C V and VI have a trparams file? About a week and a half ago I pushed a bug fix that fixed a bug where those files weren't being read which did change the ioneq quite a bit. After that fix, I checked the ioneq for all ions and found that it agreed with what was in Chianti.

However, if you're just reading from the ioneq files directly and not calculating it separately, then this is not the issue.

[wtbarnes]

Re your question about the proton rates, those are only used when there is a psplups file for that ion. If there is one, the proton excitation and de-excitation rates will be included though these don't have a huge effect usually.

[jwreep]

Neither ion has a .trparams file. I'm just using self.ioneq for these values, which I presume is just reading the files? It looks comparable to CHIANTI IDL's plot_ioneq output. The level populations also look similar to what IDL is giving.

Now I'm really confused!

[jwreep]

The revision history of the IDL code has the following comment:

;       v.21, 19-Jan-2018, Peter Young
;            This is a major update. The common blocks have been
;            removed, and I now call ch_setup_ion to load the atomic
;            data for the ions. In the process of doing this I found
;            three bugs in the original code (i) proton rates were not
;            loaded previously, yet they were important for c_6; (ii)
;            the upper levels for the He-like ions are different for
;            most of the dielectronic ions, so these are now defined
;            separately (heseq_lvl2d); (iii) for He-like dielectronic
;            ions, the routine was incorrectly using atomic data from
;            the previous ion in masterlist. I've also added a
;            check on the upper levels of the He-like ions to make
;            sure they really do correspond to the two-photon
;            transition.

Fiasco appears to be including the proton rates in the calculation for C VI, but I don't see where this is in the IDL calculation, likely because it's nested in other functions? Something like this might be the source of the problems, perhaps.

For another possibility, I don't think the function I've written has any reference to dielectronic ions. I don't see where heseq_lvl2d files are stored, either (not in the continuum folder or individual element folders).

[wtbarnes]

Re proton rates, you can turn off the proton rate calculation by passing include_protons=False to level_populations.

As far as heseq_lvl2d, isn't there some sort or hardcoded array in the IDL software itself that is named something like this? Maybe I'm misremembering.

[wtbarnes]

Also ignore my previous comment about the trparams stuff. That's definitely not the culprit. If you're just using the ioneq attribute then it's just reading ionization fractions from the file and interpolating them to the temperature array of that ion. I've carefully checked that these outputs match those in IDL.

[jwreep]

Alright, so, the proton rates appear to be the issue with C VI:

It's curious that it makes so much difference. I'm using CHIANTI v8.0.7, which predates Peter's update that added in proton rates, so I guess we need to be careful about comparisons here to make sure we have the same versions. In any case, the proton rates need to be added to this implementation (as a kwarg?) since it makes orders of magnitude difference.

C V doesn't have any proton data at all, so that can't be its issue.

The heseq_lvl2d is indeed hard-coded into later versions of two_photon.pro, but also appears to have been removed prior to the current version of CHIANTI. I'm very confused.

; PRY, 18-Jan-2018
; Here the level index of upper level of the two-photon transition of
; the helium sequence ions is hard-coded. Note that -1 indicates the
; ion does not have a model.
; Note that the index can be different for the dielectronic ions.
; PRY, 12-Nov-2020
; These are no longer used.
;
heseq_lvl2 = [-1, 3,-1,-1,-1,5,6,6,-1,6, 6,6,5,5, 3,5, 3,5, 3,5,-1,-1,-1,-1,-1,4,-1,4,-1, 4]
heseq_lvl2d =[-1,-1,-1,-1,-1,3,3,3,-1,3,-1,3,3,3,-1,3,-1,3,-1,3,-1,-1,-1,-1,-1,3,-1,3,-1,-1]

[jwreep]

On a side note, I don't see the p/e ratio in the IDL code. I wonder if that's why there's a $\frac{1}{n_{e}}$ factor in there, and perhaps why I'm getting 1e30 rather than 1e40 off from the IDL. Might not have the equation fully correct quite yet, but getting there slowly . . .

[wtbarnes]

Re: the C VI proton rates, I wonder if its related to this bug: https://groups.google.com/g/chianti/c/9dvCwFUReZ0. I'm confused though because as I understand it, that was a problem with the database, so it should manifest the same in both the IDL and Python versions. However, it seems too coincidental...

This is likely a dumb question, but in your testing, you're using the v8.0.7 of the database for both the IDL and the Python cases and you're using v8.0.7 (or some other later minor version) of the IDL software?

[jwreep]

Not a dumb question -- since it's hard to keep track!

For IDL, I'm using both the 8.0.7 data and routines. (I have the 10.1 version also installed, but I specifically downgraded SSW for the moment.)

For Python, I'm using the bundle with fiasco, which looks like it's 8.0.6.

Since the bug fix in that e-mail thread is from a year ago, presumably the data (proton rates) were only fixed in the latest versions of Chianti. So, both Python and IDL should have the same issue in my plot.

However, since fiasco uses proton rates by default in level populations, the C VI in my first plot was off. I guess the proton rate data was available but not used in the IDL two_photon.pro prior to Peter's bug fix in 2018. (Version 8 was released in 2015)

[wtbarnes]

However, since fiasco uses proton rates by default in level populations, the C VI in my first plot was off. I guess the proton rate data was available but not used in the IDL two_photon.pro prior to Peter's bug fix in 2018. (Version 8 was released in 2015)

Ah, I see. I guess then it makes sense to just turn the proton rates off for all ions in the two-photon calculation by default since it likely makes little difference for most ions and clearly gives wrong answers in the case of C VI.

[jwreep]

Alright, getting somewhere! Found a major mistake -- I was using the wrong index for the He-like levels, which I assumed was just the second transition as it is for the H-like case. I've added a line to grab the correct index by accessing the actual level configuration and angular momentum quantum number rather than making any assumptions.

The results are much more encouraging:

There are still small differences. I'm not entirely sure where to add and/or deal with dielectronic ions (as in Peter's note to the IDL code), so that might likely be the major remaining difference.

[jwreep]

I'm guessing #25 might also affect the end result of the two-photon calculation as well. It's possible that might need to be implemented prior to getting perfect agreement between Fiasco and Chianti.

[wtbarnes]

I'm guessing #25 might also affect the end result of the two-photon calculation as well. It's possible that might need to be implemented prior to getting perfect agreement between Fiasco and Chianti.

Hmm, I don't think that's likely the cause. For these effects to be included, you have to specify a particular radial distance and radiation field. I think by default in IDL, those are both zero and thus those effects are not included by default so I don't think that would play a role here.

[wtbarnes]

Those plots are encouraging! For the second plot, are the labels switched? I thought that not including the protons made the comparison better? Or am I reading this plot wrong?

[jwreep]

Oh, thanks. That's good to hear.

Yes, the labels are switched accidentally.

[wtbarnes]

Ok, good 😅. If I'm reading that right then, disregarding the factor of 1e30, it looks like w/o proton ratio is close-ish to 1 over the whole wavelength rate (though looking at the top plot I see there is some constant offset) with some small deviation near 300 Å. I wonder what is causing that offset...

I'm not entirely sure where to add and/or deal with dielectronic ions (as in Peter's note to the IDL code), so that might likely be the major remaining difference.

That's possible. I've not dealt with the satellite lines at all in fiasco. In v9 and above, the manner in which these are dealt with changes dramatically (see #234 and I'm working on that in #254) by coupling two ions together in the level populations calculation (the so-called "two-ion model"), but in this version I cannot exactly recall how this is handled.

You can get the dielectronic ion (if it exists) by just appending "dielectronic" to the front of the usual dataset name, e.g. to get the satellite line energy information for Fe IX,

import fiasco
import astropy.units as u
fe_9 = fiasco.Ion('Fe IX', 1*u.MK)
fe_9._dielectronic_elvlc

These correspond to the files in the directories that are prepended by a "_d" in the database.

[jwreep]

Thanks, I'll test to see if I can get the dielectronic emission included using that.

[jwreep]

Might be a dumb question, but how do I get the level populations for the dielectronic ion? I don't see a dielectronic_level_populations() function or a switch to tell the level_populations() function to use the dielectronic data. Am I missing something?

[wtbarnes]

Not a dumb question. There's currently not a way to do this in fiasco. I thought that it was a matter of including the satellite lines within the main level populations calculation but maybe I'm mistaken?

If we're going to go down the road of chasing this dielectronic issue down, I'd suggest doing this in a separate PR. When it comes to this PR, we could either a) merge this without that change once we're sure that's the cause of the lingering difference in Chianti and IDL or b) hold off on merging this until we've resolved the dielectronic issue.

One way to diagnose whether this is the culprit would be to compute the 2-photon emission for an ion that doesn't have an accompanying "_d" directory. If the difference persists, then something else must be the issue.

[jwreep]

He II is still showing small differences (that's the bump at 304). It doesn't have any dielectronic data or proton rates, so it must be something else. Back to the drawing board for now!

[wtbarnes]

Well that's progress! I wonder how the level populations compare for this ion. I'd guess that's where the problem lies.

[jwreep]

There are errors on the level of 0.5% in the level populations. I think you might have better insight about the reason.

For He II, at $10^{10}\ \text{cm}^{-3}$, the population of $2S_{1/2}$:

IDL, $10^5$ K: 1.763e-04 Fiasco, $10^5$ K: 1.75277910e-04

IDL, $10^6$ K: 3.486e-03 Fiasco, $10^6$ K: 3.50703844e-03

[wtbarnes]

I'm definitely not worried about 0.5% in the level populations. The IDL code and fiasco use very different approaches to solving the set of equations so I'm perfectly happy with that level of agreement.

OK, so the problem in the two-photon is elsewhere. What about the abundances? I've had issues in the past with getting the IDL code to use the abundance set I want. In fiasco by default, it will use the sun_coronal_1992_feldman_ext set. Are you specifying this when you run the IDL code?

[jwreep]

That's a really good question. It's not clear to me if I am. The patch notes from the latest IDL version have this note:

;       v.22, 12-Aug-2019, Peter Young
;            The abund_file and ioneq_file inputs were not working, so
;            this has been fixed.

since v8.0.7 is older than this patch, it's possible I'm not using the abundance file I think I am. The IDL default when I open SSW is showing as sun_photospheric_1998_grevesse. When I use that in fiasco, the discrepancy is reduced, but still present.

So, yes, possibly the wrong abundance. I'll have to look at this more carefully.

[jwreep]

Not the abundance. If it were, the IDL case appears to have a helium abundance 10.93 < IDL's < 10.986, and none of the abundance sets have values in that range, in either the old or current version of the database.

[wtbarnes]

If it were, the IDL case appears to have a helium abundance 10.93 < IDL's < 10.986,

Likely a dumb question, but is it possible we're missing an order of magnitude somewhere else and the variation is actually 0.93 to 0.986?

[jwreep]

Sorry, should've made another update, but I don't fully understand what's going on yet.

The remaining discrepancy is the p/e ratio, but I don't understand why it's absent from IDL or ChiantiPy.

The density of excited ions can be written: $n{j}(X^{+m}) = \frac{n{j}(X^{+m})}{n(X^{+m})} \frac{n(X^{+m})}{n(X)} \frac{n(X)}{n{H}} \frac{n{H}}{n{e}} n{e}$ which could be rewritten in terms of either the column or volume EM $n{j}(X^{+m}) = \frac{n{j}(X^{+m})}{n(X^{+m})} \frac{n(X^{+m})}{n(X)} \frac{n(X)}{n{H}} \frac{EM{V}}{V n{e}}$ or $n{j}(X^{+m}) = \frac{n{j}(X^{+m})}{n(X^{+m})} \frac{n(X^{+m})}{n(X)} \frac{n(X)}{n{H}} \frac{EM_{c}}{s n_{e}}$

This appears to be the coding in IDL and ChiantiPy (sort of -- don't actually see a volume or column depth), but hides the temperature dependence of p/e into the assumed volume or column depth, which I don't like. I also don't understand how the units could then work out properly (seems to be missing a density-squared).

ChiantiPy:

rate[it, goodWvl] = f*pop[it, l2]*distr*ab*thisIoneq[it]*self.Em[it]/eDensity[it]

IDL:

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]]

tl;dr my coding is wrong because I was trying to make the units work out, but they've also cleverly hidden p/e ratio because they wanted to be able to use any form of EM as input.

[wtbarnes]

Huh ok. The approach you're describing in the IDL/ChiantiPy case feels a bit wonky to me and I'm definitely more inclined to use the approach we've taken here, using the proton/electron ratio explicitly. I'm also still confused about how the units actually work out in their case but maybe I'm just thinking about it wrong.

When you say "whichever form of the EM they want" do you mean a choice between column and volume EM or do you mean whether the EM represents $n_en_e$ or $n_en_H$?

[jwreep]

The former. In IDL, you can use the column EM, volume EM, or DEM as input, so I guess the function is written to accommodate that. I still can't quite wrap my head around the units, which is why I haven't committed a fix yet.

[jwreep]

The crux of my confusion is that the IDL code doesn't seem to have a volume/column depth built into it, but it seems like it should be there?

$F{2p} \propto (hc) * (A{ji} / \lambda^{2}) * (n_{j}(X^{+m}))$

where the units are $(\text{erg}\ Å) \times (\frac{1}{\text{s}\ Å^{2}}) \times (\text{cm}^{-3}) = \frac{\text{erg}}{\text{s}\ Å\ \text{cm}^{3}}$

Writing $n{j}(X^{+m}) \propto \frac{EM{V}}{V n{e}}$, then $F{2p} \propto (hc) (A_{ji} / \lambda^{2}) (\frac{1}{n{e}}) * \frac{EM{V}}{V} $

which, in terms of units is $F{2p} \propto \frac{\text{erg}\ \text{cm}^{3}}{\text{s}\ Å} * \frac{EM{V}}{V} $

or similarly, $F{2p} \propto \frac{\text{erg}\ \text{cm}^{3}}{\text{s}\ Å} * \frac{EM{c}}{h}$ for $h$ the column depth

It's a bit odd at first glance, but I think that's what going on. In short, normalizing to $\frac{EM{V}}{V} = \frac{EM{c}}{h} = 1\ \text{cm}^{-6}$.

I suppose we could drop out the p/e ratio as they've done, which will then give the same result as CHIANTI (besides small differences like level populations), but then that also drops the temperature dependence of p/e.

I think we could keep p/e(T), but we would still have to normalize by an EM, or, divide the whole function by $\frac{EM{V}}{V} = \frac{EM{c}}{h} = 1\ \text{cm}^{-6}$. (Maybe? How else could you get to units of $\frac{\text{erg}\ \text{cm}^{3}}{\text{s}\ Å}$?) We would then just have to be aware in comparisons that the function is off by this factor.

The weird part is that this implies the normalization in CHIANTI depends on temperature . . .

I think that's the gist of it.

[wtbarnes]

I think the confusion here may be due to the fact that the equation 11 in Young et al. already has an assumed factor of $EM_C/h$ (or $EM_V/V$) built into it. That is, the two photon intensity in the case of a LOS column integration is,

$$ I_{2p} = \int\mathrm{d}h \frac{d\epsilon}{d\lambda} $$

which has units of erg cm-2 sr-1 Å-1. To be consistent with what is already done in fiasco for both the contribution function and the free-free and free-bound continua, I suggest what we do is express the two-photon continuum as,

$$ P_{2p} = \frac{d\epsilon}{d\lambda}\frac{4\pi}{ne^2} = \frac{hc}{\lambda}A{ij}\frac{n_j(X^{+m})}{n(X^{+m})}\frac{n(X^{+m})}{n(X)}\frac{n(X)}{n(H)}\frac{n(H)}{n_e}\phi(\lambda_0/\lambda)\frac{1}{n_e} $$

which has units of erg cm3 s-1 Å-1. I've also dropped the $1/4\pi$ factor to be consistent as well such that,

$$ I_{2p} = \int\mathrm{d}h ne^2 P{2p} \approx EMC P{2p} $$

Writing this all out, I'm now realizing that we could just compute the contribution function (i.e. literally use Ion.contribution_function), select the relevant transition, and multiply by $\phi\lambda_0/\lambda$. Actually, we wouldn't want to do this because the A values are different here and also to be consistent with what is done in free_free and free_bound, I guess we should also drop the abundance and ionization fraction components as well.