Add two-photon continuum (#260)

* create two_photon continuum function

* fix spline

* add logic to observed vs theoretical wavelength

* remove units from A_sum

* add 2p method for IonCollection

* allow wavelength, temperature, density to be arrays

* bug fixes for array handling

* error handling for ions that aren't H- or He-like

* add needs_dataset decorator

* fix variable description in IDL tests

* fix code to correct units of EM and add missing 1/wvl

* change name of A_sum and add note in fiasco.io about its origin

* add include_protons option to two_photon calculation, default to True

* fix He-like level index

* change normalization

* revert to 4539f17

* fix normalization, now with current main tree

* add description of the function and references

* refactor EM, and move ioneq and abund to ioncollection

* fix LateX in description of 2p function

* add tests for ion and collections

* fix ion collection test

* add he_2.wgfa to test file list

* add he_2.scups to test file list

* typo

* test fix for test_ion.py -- roundoff error?

* text fix for test_collection.py -- roundoff error?

* add tests for helium-like ion and an ion missing input dataset

* fix bug in test_ion.test_two_photon

* add elvlc for c4 and c5 to test list

* add wgfa and scups for c4 and c5 to test list

* set include_protons to False by default

* update value for C V test

* add c_5.reclvl to test database

* require dbase_Version <= 8.0.7 for two_photon test

* minor refactoring of 2p continuum

* change comparison with == to np.allclose

* allclose -> isclose, oops

* add comment to docstring about no photons beneath rest wavelength

---------

Co-authored-by: Jeffrey Reep <reep@atrcl174.IfA.Hawaii.edu>
Co-authored-by: Jeffrey Reep <reep@kahiku.local>
Co-authored-by: Jeffrey Reep <reep@atrcw9.IfA.Hawaii.Edu>
Co-authored-by: Jeffrey Reep <reep@atrcw10.ifa.hawaii.edu>
Co-authored-by: Jeffrey Reep <reep@atrcw11.IfA.Hawaii.Edu>
Co-authored-by: Jeffrey Reep <reep@atrcw17.IfA.Hawaii.Edu>
Co-authored-by: Jeffrey Reep <reep@atrcw26.IfA.Hawaii.Edu>
Co-authored-by: Jeffrey Reep <reep@atrcw27.IfA.Hawaii.Edu>
Co-authored-by: Will Barnes <will.t.barnes@gmail.com>
Co-authored-by: Jeffrey Reep <reep@hawaii.edu>

docs/references.bib

@@ -157,6 +157,18 @@ @article{dere_ionization_2007
   doi = {10.1051/0004-6361:20066728}
 }
 
+@ARTICLE{drake_twophoton_1986,
+  title = "{Spontaneous two-photon decay rates in hydrogenlike and heliumlike ions}",
+  author = {{Drake}, G.~W.~F.},
+  journal = {Phys. Rev. A},
+  year = {1986},
+  month = oct,
+  volume = {34},
+  number = {4},
+  pages = {2871-2880},
+  doi = {10.1103/PhysRevA.34.2871}
+}
+
 @article{feldman_potential_1992,
   title = {The Potential for Plasma Diagnostics from Stellar Extreme-Ultraviolet Observations},
   author = {Feldman, U. and Mandelbaum, P. and Seely, J. F. and Doschek, G. A. and Gursky, H.},
@@ -181,6 +193,16 @@ @article{fontes_fully_1999
   doi = {10.1103/PhysRevA.59.1329}
 }
 
+@ARTICLE{gronenschild_twophoton_1978,
+  author = {{Gronenschild}, E.~H.~B.~M. and {Mewe}, R.},
+  title = "{Calculated X-radiation from optically thin plasmas. III. Abundance effects on continuum emission.}",
+  journal = {Astronomy and Astrophysics Supplement Series},
+  year = {1978},
+  month = may,
+  volume = {32},
+  pages = {283-305}
+}
+
 @article{itoh_relativistic_2000,
   title = {Relativistic {{Thermal Bremsstrahlung Gaunt Factor}} for the {{Intracluster Plasma}}. {{II}}. {{Analytic Fitting Formulae}}},
   author = {Itoh, Naoki and Sakamoto, Tsuyoshi and Kusano, Shugo and Nozawa, Satoshi and Kohyama, Yasuharu},

fiasco/collections.py

@@ -167,6 +167,48 @@ def free_bound(self, wavelength: u.angstrom, **kwargs):
                 free_bound += fb * abundance * ioneq[:, np.newaxis]
         return free_bound
 
+    @u.quantity_input
+    def two_photon(self, wavelength: u.angstrom, electron_density: u.cm**-3, **kwargs):
+        r"""
+        Compute the two-photon continuum emission.
+
+        .. note:: Both abundance and ionization equilibrium are included here.
+
+        The combined two-photon continuum is given by
+
+        .. math::
+
+            P_{2p}(\lambda,T,n_{e}) = \sum_{X,k}\mathrm{Ab}(X)f(X_{k})C_{2p, X_k}(\lambda,T,n_{e})
+
+        where :math:`\mathrm{Ab}(X)` is the abundance of element :math:`X`,
+        :math:`f(X_{k})` is the ionization equilibrium of the emitting ion :math:`X_{k}`,
+        and :math:`C_{fb, X_k}(\lambda,T)` is the two-photon emission of the
+        ion :math:`X_k` as computed by `fiasco.Ion.two_photon`.
+        The sum is taken over all ions in the collection.
+
+        Parameters
+        ----------
+        wavelength : `~astropy.units.Quantity`
+        electron_density: `~astropy.units.Quantity`
+
+        See Also
+        --------
+        fiasco.Ion.two_photon
+        """
+        wavelength = np.atleast_1d(wavelength)
+        electron_density = np.atleast_1d(electron_density)
+        two_photon = u.Quantity(np.zeros(wavelength.shape + self.temperature.shape + electron_density.shape),
+                                'erg cm^3 s^-1 Angstrom^-1')
+        for ion in self:
+            try:
+                tp = ion.two_photon(wavelength, electron_density, **kwargs)
+            except MissingDatasetException as e:
+                self.log.warning(f'{ion.ion_name} not included in two-photon emission. {e}')
+                continue
+            else:
+                two_photon += tp * ion.abundance * ion.ioneq[:, np.newaxis]
+        return two_photon
+
     @u.quantity_input
     def spectrum(self, density: u.cm**(-3), emission_measure: u.cm**(-5), wavelength_range=None,
                  bin_width=None, kernel=None, **kwargs):

fiasco/io/sources/continuum_sources.py

@@ -252,15 +252,22 @@ def to_hdf5(self, hf, df, **kwargs):
 
 
 class HSeqParser(GenericParser):
-    """
+    r"""
     Parameters for calculating two-photon continuum for hydrogen-like ions
+
+    Notes
+    -----
+    * The parameter :math: `\psi_{\text{norm}}` (called :math: `A_{\text{sum}}` in CHIANTI) is a normalization factor of
+    the integral of the spectral distribution function :math:`\psi(y)` from 0 to 1, such
+    that :math: `\frac{1}{\psi_{\text{norm}}} \int_{0}^{1} \psi(y) dy = 2`.  This normalization is only used For
+    hydrogenic ions.
     """
     filetype = 'hseq_2photon'
     dtypes = [int, float, int, float, float, float]
-    units = [None, u.dimensionless_unscaled, None, 1/u.s, 1/u.s, u.dimensionless_unscaled]
-    headings = ['Z', 'y', 'Z_0', 'A', 'A_sum', 'psi']
+    units = [None, u.dimensionless_unscaled, None, 1/u.s, u.dimensionless_unscaled, u.dimensionless_unscaled]
+    headings = ['Z', 'y', 'Z_0', 'A', 'psi_norm', 'psi']
     descriptions = ['atomic number', 'fraction of energy carried by one of the two photons',
-                    'nominal atomic number', 'radiative decay rate', 'summed radiative decay rate',
+                    'nominal atomic number', 'radiative decay rate', 'normalization of the integral of psi from 0 to 1',
                     'spectral distribution function']
 
     def __init__(self, filename, **kwargs):

fiasco/ions.py

@@ -6,7 +6,7 @@
 import numpy as np
 
 from functools import cached_property
-from scipy.interpolate import interp1d, PchipInterpolator, splev, splrep
+from scipy.interpolate import CubicSpline, interp1d, PchipInterpolator, splev, splrep
 from scipy.ndimage import map_coordinates
 
 from fiasco import proton_electron_ratio
@@ -1418,7 +1418,7 @@ def free_bound(self,
             C_{fb}(\lambda, T) = \frac{2}{hc^3(k_B m_e)^{3/2}\sqrt{2\pi}}\frac{E^5}{T^{3/2}}\sum_i\frac{\omega_i}{\omega_0}\sigma_i^{\mathrm{bf}}\exp{\left(-\frac{E-I_i}{k_BT}\right)}
 
         where :math:`E` is the energy of the outgoing photon,
-        :math:`\omega_i,\omega_0` are the statastical weights of the
+        :math:`\omega_i,\omega_0` are the statistical weights of the
         :math:`i`-th level of the recombined ion and the ground level of the recombining ion, respectively,
         :math:`\sigma_i^{\mathrm{bf}}` is the free-bound cross-section,
         and :math:`I_i` is the energy required to ionize the recombined ion from level :math:`i`.
@@ -1527,3 +1527,90 @@ def _karzas_cross_section(self, photon_energy: u.erg, ionization_energy: u.erg,
         cross_section = prefactor * ionization_energy**2 * photon_energy**(-3) * gaunt_factor / n
         cross_section[np.where(photon_energy < ionization_energy)] = 0.*cross_section.unit
         return cross_section
+
+    @needs_dataset('elvlc')
+    @u.quantity_input
+    def two_photon(self,
+                   wavelength: u.angstrom,
+                   electron_density: u.cm**(-3),
+                   include_protons=False) -> u.Unit('erg cm3 s-1 Angstrom-1'):
+        r"""
+        Two-photon continuum emission of a hydrogenic or helium-like ion.
+
+        .. note:: Does not include ionization equilibrium or abundance.
+                  Unlike the equivalent IDL routine, the output here is not
+                  expressed per steradian and as such the factor of
+                  :math:`1/4\pi` is not included.
+
+        In hydrogen-like ions, the transition :math:`2S_{1/2} \rightarrow 1S_{1/2} + h\nu` cannot occur
+        as an electric dipole transition, but only as a much slower magnetic dipole transition.
+        The dominant transition then becomes :math:`2S_{1/2} \rightarrow 1S_{1/2} + h\nu_{1} + h\nu_{2}`.
+
+        In helium-like ions, the transition from :math:`1s2s ^{1}S_{0} \rightarrow 1s^{2}\ ^{1}S_{0} + h\nu`
+        is forbidden under quantum selection rules since :math:`\Delta J = 0`.
+        Similarly, the dominant transition becomes
+        :math:`1s2s ^{1}S_{0} \rightarrow 1s^{2}\ ^{1}S_{0} + h\nu_{1} + h\nu_{2}`.
+
+        In both cases, the energy of the two photons emitted equals the energy difference of the two levels.
+        As a consequence, no photons can be emitted beneath the rest wavelength for a given transition.
+        See the introduction of :cite:t:`drake_1986` for a concise description of the process.
+
+        The emission is given by
+
+        .. math::
+
+            C_{2p}(\lambda, T, n_{e}) = hc \frac{n_{j}(X^{+m}) A_{ji} \lambda_{0} \psi(\frac{\lambda_{0}}{\lambda})}{\psi_{\text{norm}}\lambda^{3}}
+
+        where :math:`\lambda_{0}` is rest wavelength of the (disallowed) transition,
+        :math:`A_{ji}` is the Einstein spontaneous emission coefficient,
+        :math:`\psi` is so-called spectral distribution function, given approximately by
+        :math:`\psi(y) \approx 2.623 \sqrt{\cos{\Big(\pi(y-\frac{1}{2})\Big)}}` :cite:p:`gronenschild_twophoton_1978`,
+        :math:`\psi_{\text{norm}}` is a normalization factor for hydrogen-like ions such
+        that :math:`\frac{1}{\psi_{\text{norm}}} \int_{0}^{1} \psi(y) dy = 2` (and 1 for helium-like ions),
+        and :math:`n_{j}(X^{+m})` is the density of ions m of element X in excited state j, given by
+        :math:`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}`.
+
+        Parameters
+        ----------
+        wavelength : `~astropy.units.Quantity`
+        electron_density : `~astropy.units.Quantity`
+        include_protons : `bool`, optional
+            If True, use proton excitation and de-excitation rates in the level population calculation.
+        """
+        wavelength = np.atleast_1d(wavelength)
+        electron_density = np.atleast_1d(electron_density)
+
+        final_shape = wavelength.shape + self.temperature.shape + electron_density.shape
+        if self.hydrogenic:
+            A_ji = self._hseq['A']
+            psi_norm = self._hseq['psi_norm']
+            cubic_spline = CubicSpline(self._hseq['y'], self._hseq['psi'])
+            config = '2s'  # Get the index of the 2S1/2 state for H-like
+            J = 0.5
+        elif self.helium_like:
+            A_ji = self._heseq['A']
+            psi_norm = 1.0 * u.dimensionless_unscaled
+            cubic_spline = CubicSpline(self._heseq['y'], self._heseq['psi'])
+            config = '1s.2s'  # Get the index of the 1s2s 1S0 state for He-like:
+            J = 0
+        else:
+            return u.Quantity(np.zeros(final_shape),  'erg cm^3 s^-1 Angstrom^-1')
+        level_index = np.where((self._elvlc['config'] == config) & (np.isclose(self._elvlc['J'], J)) )[0][0]
+
+        E_obs = self._elvlc['E_obs'][level_index]
+        E_th = self._elvlc['E_th'][level_index]
+        E_2p = E_obs if E_obs > 0.0 else E_th
+        rest_wavelength = 1 / E_2p
+
+        psi_interp = cubic_spline((rest_wavelength / wavelength).decompose())
+        psi_interp = np.where(wavelength < rest_wavelength, 0.0, psi_interp)
+
+        energy_dist = (A_ji * rest_wavelength * psi_interp) / (psi_norm * wavelength**3)
+
+        level_population = self.level_populations(electron_density, include_protons=include_protons)
+        level_population = level_population[..., level_index]
+
+        level_density = level_population / electron_density
+        matrix = np.outer(energy_dist, level_density).reshape(*final_shape)
+
+        return const.h * const.c * matrix

fiasco/tests/idl/helpers.py

@@ -46,7 +46,7 @@ def run_idl_script(idl_env, script, input_args, save_vars, file_name, version, f
     file_name: `str`
         Name of the IDL results file
     version: `packaging.version.Version`, `str`
-        Version of IDL used to generate these test results
+        Version of CHIANTI used to generate these test results
     format_func: `dict`, optional
         Functions to use to format output from the IDL function.
         This is most useful for adding the necessary units to any of the outputs.

fiasco/tests/test_collections.py

@@ -106,6 +106,19 @@ def test_free_bound(another_collection, wavelength):
     index_t = 24  # This is approximately where the ioneq for Fe V peaks
     assert u.allclose(fb[index_t, index_w], 3.057781475607237e-36 * u.Unit('erg cm3 s-1 Angstrom-1'))
 
+@pytest.mark.parametrize('wavelength', [wavelength, wavelength[450]])
+def test_two_photon(collection, wavelength, hdf5_dbase_root):
+    # add Li III to the test to include an ion that throws a MissingDatasetException
+    collection = collection + fiasco.Ion('Li III', collection.temperature, hdf5_dbase_root=hdf5_dbase_root)
+    tp = collection.two_photon(wavelength, electron_density = 1e10 * u.cm**(-3))
+    if wavelength.shape:
+        assert tp.shape == wavelength.shape + temperature.shape + (1, )
+    else:
+        assert tp.shape == (1, ) + temperature.shape + (1, )
+    index_w = 450 if wavelength.shape else 0
+    index_t = 30
+    # This value has not been checked for correctness
+    assert u.allclose(tp[index_w, index_t, 0], 3.48580645e-27 * u.Unit('erg cm3 s-1 Angstrom-1'))
 
 def test_radiative_loss(collection):
     rl = collection.radiative_loss(1e9*u.cm**(-3))

fiasco/tests/test_ion.py

@@ -32,6 +32,16 @@ def fe10(hdf5_dbase_root):
     return fiasco.Ion('Fe 10', temperature, hdf5_dbase_root=hdf5_dbase_root)
 
 
+@pytest.fixture
+def c4(hdf5_dbase_root):
+    return fiasco.Ion('C IV', temperature, hdf5_dbase_root=hdf5_dbase_root)
+
+
+@pytest.fixture
+def c5(hdf5_dbase_root):
+    return fiasco.Ion('C V', temperature, hdf5_dbase_root=hdf5_dbase_root)
+
+
 @pytest.fixture
 def c6(hdf5_dbase_root):
     return fiasco.Ion('C VI', temperature, hdf5_dbase_root=hdf5_dbase_root)
@@ -323,6 +333,23 @@ def test_free_bound(ion):
     # This value has not been tested for correctness
     assert u.allclose(emission[0, 0], 9.7902609e-26 * u.cm**3 * u.erg / u.Angstrom / u.s)
 
+# The two-photon test currently fails for dbase_version >= 9 because it is missing c_5.reclvl
+@pytest.mark.requires_dbase_version('<=','8.0.7')
+def test_two_photon(c4, c5, c6):
+    # test both H-like and He-like ions, and one that doesn't have two-photon emission
+    c4_emission = c4.two_photon(200 * u.Angstrom, electron_density = 1e10* u.cm**(-3))
+    c5_emission = c5.two_photon(200 * u.Angstrom, electron_density = 1e10* u.cm**(-3))
+    c6_emission = c6.two_photon(200 * u.Angstrom, electron_density = 1e10* u.cm**(-3))
+    c6_emission_protons = c6.two_photon(200 * u.Angstrom, electron_density = 1e10* u.cm**(-3), include_protons=True)
+    assert c4_emission.shape == (1, ) + c4.temperature.shape + (1, )
+    assert c5_emission.shape == (1, ) + c5.temperature.shape + (1, )
+    assert c6_emission.shape == (1, ) + c6.temperature.shape + (1, )
+    assert c6_emission_protons.shape == (1, ) + c6.temperature.shape + (1, )
+    assert u.allclose(c4_emission[0, 30, 0], 0.0 * u.cm**3 * u.erg / u.Angstrom / u.s)
+    # These values have not been tested for correctness
+    assert u.allclose(c5_emission[0, 30, 0], 4.04634243e-25 * u.cm**3 * u.erg / u.Angstrom / u.s)
+    assert u.allclose(c6_emission[0, 30, 0], 8.25316887e-26 * u.cm**3 * u.erg / u.Angstrom / u.s)
+    assert u.allclose(c6_emission_protons[0, 30, 0], 6.79615958e-29 * u.cm**3 * u.erg / u.Angstrom / u.s)
 
 def test_free_bound_no_recombining(h1):
     # This is test the case where there is no data available for the recombining

fiasco/util/data/test_file_list.json

@@ -169,12 +169,19 @@
     "c_4.diparams",
     "c_4.drparams",
     "c_4.easplups",
+    "c_4.elvlc",
     "c_4.fblvl",
     "c_4.rrparams",
+    "c_4.scups",
+    "c_4.wgfa",
     "c_5.diparams",
     "c_5.drparams",
+    "c_5.elvlc",
     "c_5.fblvl",
+    "c_5.reclvl",
     "c_5.rrparams",
+    "c_5.scups",
+    "c_5.wgfa",
     "c_6.diparams",
     "c_6.drparams",
     "c_6.elvlc",
@@ -811,6 +818,8 @@
     "he_2.elvlc",
     "he_2.fblvl",
     "he_2.rrparams",
+    "he_2.scups",
+    "he_2.wgfa",
     "he_3.rrparams",
     "heseq_2photon.dat",
     "hseq_2photon.dat",