Add continuum contributions to IonCollection.radiative_loss (#282)

* update docstring for radiative_loss * update radiative loss method * add Li III to test to throw missing data exception * fix bug in test * add function to calculate total f-f gaunt factor * add f-f radiative loss function * start adding f-b rad loss * add note to f-b Gaunt factor - f1 is analytic * add g_fb_total -- needs to be checked * add fb rad loss + fix fb Gaunt * add new functions to full loss * small bug fixes * add units tests, remove collection test for the moment * pin numpy < 2.0 * fix recombined/recombining in fb radiative losses * add some tests * small edit to gaunt_fb * fix missing ions in f-b * fix up tests * add one more test * remove numpy version pin * refactor rad loss into three separate functions * handle possible overflows in gf exp term * add idl comparison tests for ff and fb rad losses * fix gf test * Update fiasco/ions.py Co-authored-by: Will Barnes <will.t.barnes@gmail.com> * Update fiasco/ions.py Co-authored-by: Will Barnes <will.t.barnes@gmail.com> * Update fiasco/ions.py Co-authored-by: Will Barnes <will.t.barnes@gmail.com> * Update fiasco/ions.py Co-authored-by: Will Barnes <will.t.barnes@gmail.com> * Update fiasco/ions.py Co-authored-by: Will Barnes <will.t.barnes@gmail.com> * Update fiasco/ions.py Co-authored-by: Will Barnes <will.t.barnes@gmail.com> * fix up small issues * fix zeta0 tests * add topic guide for f-b Gaunt factor * add collection tests + doc strings * fix typos in tests * remove li 3 from f-f and f-b tests -- doesn't do anything * fix refs; fix topic guide formatting * change method names to be consistent with ion methods * move two-photon description to topic guide; fix up refs * fix the tests --------- Co-authored-by: Jeffrey Reep <reep@atrcw30.ifa.hawaii.edu> Co-authored-by: Jeffrey Reep <reep@atrcw40.ifa.hawaii.edu> Co-authored-by: Jeffrey Reep <reep@kahiku.local> Co-authored-by: Will Barnes <will.t.barnes@gmail.com> Co-authored-by: Jeffrey Reep <reep@atrcw4.IfA.Hawaii.Edu>
wtbarnes · Aug 13, 2024 · 8ee514b · 8ee514b
1 parent 7cba8b8
commit 8ee514b
Show file tree

Hide file tree

Showing 13 changed files with 564 additions and 69 deletions.
diff --git a/docs/references.bib b/docs/references.bib
@@ -157,16 +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},
+@article{drake_spontaneous_1986,
+  title = {Spontaneous Two-Photon Decay Rates in Hydrogenlike and Heliumlike Ions},
+  author = {Drake, G. W. F.},
   year = {1986},
   month = oct,
+  journal = {Physical Review A},
   volume = {34},
   number = {4},
-  pages = {2871-2880},
-  doi = {10.1103/PhysRevA.34.2871}
+  pages = {2871--2880},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRevA.34.2871},
+  urldate = {2024-08-12}
 }
 
 @article{feldman_potential_1992,
@@ -193,16 +195,19 @@ @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},
+@article{gronenschild_calculated_1978,
+  title = {Calculated {{X-radiation}} from Optically Thin Plasmas. {{III}} - {{Abundance}} Effects on Continuum Emission},
+  author = {Gronenschild, E. H. B. M. and Mewe, R.},
   year = {1978},
   month = may,
+  journal = {Astronomy and Astrophysics Supplement Series},
   volume = {32},
-  pages = {283-305}
+  pages = {283--305},
+  issn = {0365-0138},
+  urldate = {2017-04-17}
 }
 
+
 @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},
@@ -311,6 +316,19 @@ @article{landi_chiantiatomic_2013
   doi = {10.1088/0004-637X/763/2/86}
 }
 
+@article{mewe_calculated_1986,
+  title = {Calculated {{X-radiation}} from Optically Thin Plasmas. {{VI}} - {{Improved}} Calculations for Continuum Emission and Approximation Formulae for Nonrelativistic Average {{Gaunt}} Factors},
+  author = {Mewe, R. and Lemen, J. R. and {van den Oord}, G. H. J.},
+  year = {1986},
+  month = sep,
+  journal = {Astronomy and Astrophysics Supplement Series},
+  volume = {65},
+  pages = {511--536},
+  issn = {0365-0138},
+  urldate = {2017-04-17}
+}
+
+
 @book{phillips_ultraviolet_2008,
   title = {Ultraviolet and {{X-ray Spectroscopy}} of the {{Solar Atmosphere}}},
   author = {Phillips, Kenneth J. H. and Feldman, Uri and Landi, Enrico},
@@ -447,3 +465,23 @@ @techreport{young_chianti_2019
   urldate = {2024-01-09},
   institution = {None}
 }
+
+@techreport{young_chianti_2019-1,
+  title = {{{CHIANTI Technical Report No}}. 9: {{The}} Radiative Loss Function},
+  author = {Young, Peter R.},
+  year = {2019},
+  month = may,
+  number = {9},
+  urldate = {2024-08-13},
+  institution = {None}
+}
+
+@techreport{young_chianti_2021,
+  title = {{{CHIANTI Technical Report No}}. 12: {{Computing}} the Free-Bound Continuum},
+  author = {Young, Peter R.},
+  year = {2021},
+  month = feb,
+  number = {12},
+  langid = {english},
+  institution = {None}
+}
diff --git a/docs/topic_guides/freebound_gaunt_factor.rst b/docs/topic_guides/freebound_gaunt_factor.rst
@@ -0,0 +1,34 @@
+.. _fiasco-topic-guide-freebound-gaunt-factor:
+
+The implementation of the total free-bound Gaunt factor in fiasco
+=================================================================
+
+The calculation of the wavelength-averaged total free-bound Gaunt factor, :math:`G_{fb}`,
+as suggested by :cite:t:`mewe_calculated_1986` and used in CHIANTI to compute the total free-bound
+radiative losses (:meth:`fiasco.Ion.free_bound_radiative_loss`) makes an approximation to an infinite
+sum that can be improved upon.
+Specifically, Equation 14 of :cite:t:`mewe_calculated_1986` has a simple analytic solution.
+They make the approximation,
+
+.. math::
+
+    f_{1}(Z, n, n_{0} ) = \sum_{1}^{\infty} n^{-3} - \sum_{1}^{n_{0}} n^{-3} = \zeta(3) - \sum_{1}^{n_{0}} n^{-3} \approx 0.21 n_{0}^{-1.5}
+
+where :math:`\zeta(x)` is the Riemann zeta function.
+However, the second sum is analytic,
+
+.. math::
+
+    \sum_{1}^{n_{0}} n^{-3} = \zeta(3) + \frac{1}{2}\psi^{(2)}(n_{0}+1)
+
+where :math:`\psi^{n}(x)` is the :math:`n`-th derivative of the digamma function (or polygamma function).
+As such, we can write the full solution of Equation 14 as,
+
+.. math::
+
+    f_{1}(Z, n, n_{0}) = \zeta(3) - \sum_{1}^{n_{0}} n^{-3} = - \frac{1}{2}\psi^{(2)}(n_{0}+1)
+
+The final expression is therefore simplified and more accurate than the approximation used by :cite:t:`mewe_calculated_1986`.
+
+Note also that, unlike in :cite:t:`mewe_calculated_1986`, the calculation of the total free-bound Gaunt factor,
+used by :meth:`~fiasco.Ion.free_bound_radiative_loss`, does not include the abundances and ionization equilibria.
diff --git a/docs/topic_guides/index.rst b/docs/topic_guides/index.rst
@@ -7,3 +7,5 @@ Topic Guides
   :maxdepth: 1
 
   database_format
+  freebound_gaunt_factor
+  two_photon_continuum
diff --git a/docs/topic_guides/two_photon_continuum.rst b/docs/topic_guides/two_photon_continuum.rst
@@ -0,0 +1,48 @@
+.. _fiasco-topic-guide-two-photon-continuum:
+
+Calculating the two-photon continuum emission
+=============================================
+
+This topic guide provides a detailed description of the two-photon continuum emission calculation
+as implemented in :meth:`fiasco.Ion.two_photon`.
+
+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_spontaneous_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)}}
+
+according to :cite:p:`gronenschild_calculated_1978` and :math:`\psi_{\text{norm}}` is a normalization factor such that
+
+.. math::
+
+    \frac{1}{\psi_{\text{norm}}} \int_{0}^{1} \psi(y) dy = 2
+
+for hydrogen-like ions and :math:`1` for helium-like ions.
+Finally, :math:`n_{j}(X^{+m})` is the density of ions with charge state :math:`m` of element :math:`X`
+in excited state :math:`j` and is 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}.
diff --git a/fiasco/collections.py b/fiasco/collections.py
@@ -224,7 +224,7 @@ def spectrum(self, density: u.cm**(-3), emission_measure: u.cm**(-5), wavelength
         wavelength_range : `~astropy.units.Quantity`, optional
             Tuple of bounds on which transitions to include. Default includes all
         bin_width : `~astropy.units.Quantity`, optional
-            Wavelength resolution to bin intensity values. Default to 1/10 of range
+            Wavelength resolution to bin intensity values. Default to 1/100 of range
         kernel : `~astropy.convolution.kernels.Model1DKernel`, optional
             Convolution kernel for computing spectrum. Default is gaussian kernel with thermal width
 
@@ -303,18 +303,105 @@ def spectrum(self, density: u.cm**(-3), emission_measure: u.cm**(-5), wavelength
         return wavelength, spectrum
 
     @u.quantity_input
-    def radiative_loss(self, density: u.cm**(-3), **kwargs):
+    def radiative_loss(self, density: u.cm**(-3), **kwargs) -> u.Unit('erg cm3 s-1'):
+        r"""
+        Calculate the total wavelength-integrated radiative loss rate including the
+        bound-bound, free-bound, and free-free emission contributions
+
+        .. note::  The calculation does not include two-photon continuum emission, which is also
+                   neglected in the CHIANTI IDL routines.
+
+        Parameters
+        ----------
+        density : `~astropy.units.Quantity`
+            Electron number density
+
+        Returns
+        -------
+        rad_loss_total : `~astropy.units.Quantity`
+            The total bolometric radiative loss rate
         """
-        Calculate radiative loss curve which includes multiple ions
+        rad_loss_bound_bound = self.bound_bound_radiative_loss(density, **kwargs)
+        rad_loss_free_free = self.free_free_radiative_loss()
+        rad_loss_free_bound = self.free_bound_radiative_loss()
+
+        rad_loss_total = (rad_loss_bound_bound
+                          + rad_loss_free_free[:, np.newaxis]
+                          + rad_loss_free_bound[:, np.newaxis])
+
+        return rad_loss_total
+
+    @u.quantity_input
+    def bound_bound_radiative_loss(self, density, **kwargs) -> u.Unit('erg cm3 s-1'):
+        r"""
+        Calculate the radiative loss rate from bound-bound emission (line emission)
+        integrated over wavelength.
+
+        Parameters
+        ----------
+        density : `~astropy.units.Quantity`
+            Electron number density
+
+        Returns
+        -------
+        rad_loss : `~astropy.units.Quantity`
+            The bolometric bound-bound radiative loss rate per unit emission measure
         """
         density = np.atleast_1d(density)
         rad_loss = u.Quantity(np.zeros(self.temperature.shape + density.shape), 'erg cm^3 s^-1')
         for ion in self:
             try:
                 g = ion.contribution_function(density, **kwargs)
             except MissingDatasetException as e:
-                self.log.warning(f'{ion.ion_name} failed to be added to the contribution function: {e}')
+                self.log.warning(f'{ion.ion_name} not included in the bound-bound emission. {e}')
                 continue
             rad_loss += g.sum(axis=2)
-
         return rad_loss
+
+    @u.quantity_input
+    def free_free_radiative_loss(self) -> u.Unit('erg cm3 s-1'):
+        r"""
+        Calculate the radiative loss rate from free-free emission (bremsstrahlung)
+        integrated over wavelength.
+
+        Returns
+        -------
+        rad_loss : `~astropy.units.Quantity`
+            The bolometric free-free radiative loss rate per unit emission measure
+        """
+        free_free = u.Quantity(np.zeros(self.temperature.shape), 'erg cm^3 s^-1')
+        for ion in self:
+            try:
+                ff = ion.free_free_radiative_loss()
+                abundance = ion.abundance
+                ioneq = ion.ioneq
+            except MissingDatasetException as e:
+                self.log.warning(f'{ion.ion_name} not included in free-free radiative loss. {e}')
+                continue
+            else:
+                free_free += ff * abundance * ioneq
+        return free_free
+
+    @u.quantity_input
+    def free_bound_radiative_loss(self) ->  u.Unit('erg cm3 s-1'):
+        r"""
+        Calculate the radiative loss rate from free-bound emission (collisional recombination)
+        integrated over wavelength.
+
+        Returns
+        -------
+        rad_loss : `~astropy.units.Quantity`
+            The bolometric free-bound radiative loss rate per unit emission measure
+        """
+        free_bound = u.Quantity(np.zeros(self.temperature.shape), 'erg cm^3 s^-1')
+        for ion in self:
+            try:
+                fb = ion.free_bound_radiative_loss()
+                abundance = ion.abundance
+                ioneq = ion.ioneq
+            except MissingDatasetException as e:
+                self.log.warning(f'{ion.ion_name} not included in free-bound radiative loss. {e}')
+                continue
+            else:
+                free_bound += fb * abundance * ioneq
+        return free_bound