Add continuum contributions to `IonCollection.radiative_loss`

[jwreep]

Fixes #232

I've added f-f and f-b , and 2-p continuum emission to the radiative loss calculation.

It appears to be working, but I'd like to do a bit more testing.

One potential issue: the calculation currently assumes a fixed bin width. It might be nice to generalize and do the integral properly to allow for non-fixed bins, particularly since it can be integrated over a very large range.

[codecov[bot]]

Codecov Report

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

Project coverage is 92.83%. Comparing base (f7e2560) to head (956a1c6).

Additional details and impacted files

```diff @@ Coverage Diff @@ ## main #282 +/- ## ========================================== + Coverage 92.38% 92.83% +0.44% ========================================== Files 38 37 -1 Lines 2904 3013 +109 ========================================== + Hits 2683 2797 +114 + Misses 221 216 -5 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[jwreep]

The output from all ions looks reasonable, except for log T > 7. There should be an increase with temperature from free-free emission (I thought?).

Edit: The calculation of radiative losses in CHIANTI uses a wavelength-averaged Gaunt factor. Perhaps this is the cause of the discrepancy.

[wtbarnes]

Hmm do you have a comparison of the IDL and Python cases? I definitely agree that you should see the uptick in the losses at high temperatures. I had forgotten that IDL computes the radiative loss contribution from the continuum in a slightly different way. I would think these should be consistent though...

Also, if you just plot the wavelength-averaged ff contribution, does that show an increase at high temperatures?

[jwreep]

I haven't had much time to come back to this since I posted. The CHIANTI User's guide shows the uptick, but that plot is from v6 of the database.

I don't know if it's the Gaunt factor causing the discrepancy, but I've tried changing the wavelength range and binning without any change. If there's a programming mistake here, it's a subtle one!

[wtbarnes]

I wonder if this is a consequence of not including the lower wavelengths that make the dominant contributions at high temperatures? For example, computing just the free-free contribution and averaging over the wavelength dimension,

import astropy.units as u
import fiasco
import numpy as np
import matplotlib.pyplot as plt
from astropy.visualization import quantity_support

temperature = np.logspace(5,8,100) * u.K
all_ions = fiasco.IonCollection(*[fiasco.Ion(i, temperature) for i in fiasco.list_ions()])
wavelength = np.logspace(-2,2,100) * u.AA
ff_short_wave = all_ions.free_free(wavelength).mean(axis=1)
wavelength = np.logspace(1,2,100) * u.AA
ff_long_wave = all_ions.free_free(wavelength).mean(axis=1)

with quantity_support():
    plt.plot(temperature, ff_short_wave, label=r'$\lambda_{min}=0.01$ Å')
    plt.plot(temperature, ff_long_wave, label=r'$\lambda_{min}=10$ Å')
plt.xscale('log')
plt.yscale('log')
plt.legend()

gives

[jwreep]

I tried to redo the function this morning, and while the f-f does have an uptick, I think the b-b is dominating the total losses.

Here's how I changed the continuum calculation:

        wavelength = np.logspace(-2,4,200) * u.angstrom

        ff = self.free_free(wavelength, **kwargs)
        fb = self.free_bound(wavelength, **kwargs)
        tp = self.two_photon(wavelength, density, **kwargs)

        for i in range(len(density)):
            for j in range(len(self.temperature)):
                rad_loss[j,i] += np.trapz(ff[j,:], wavelength)
                rad_loss[j,i] += np.trapz(fb[j,:], wavelength)
                rad_loss[j,i] += np.trapz(tp[:,j,i], wavelength)

and I broke it apart by components as well:

It looks like f-f is still small compared to the total, which might be because b-b is too large or one or all of the continua are too small. Alternatively, maybe this is right?

[wtbarnes]

That plot almost looks like the FF and FB are missing some scaling factor. I'm surprised that they are so far below the bound-bound losses as well as being comparable to the two-photon losses at most temperatures.

[jwreep]

I think I figured it out. The discrepancy might be the integration. Changing the code to ignore the bin width

        ff = self.free_free(wavelength, **kwargs).sum(axis=1)*u.angstrom
        fb = self.free_bound(wavelength, **kwargs).sum(axis=1)*u.angstrom
        tp = self.two_photon(wavelength, density, **kwargs).sum(axis=0)*u.angstrom

I get:

This isn't correct (fb > ff), but there's something about the integration that's not working.

In IDL ff_rad_loss.pro, the calculation is

rad_loss[k]=rad_loss[k]+this_abund*this_ioneq[k]*factor*sqrt(t[k])*z^2*gaunt[k]

which is just a straight sum over all the ions (no integration in wavelength space).

So, I guess the question is why does integration not seem to work? Shouldn't the total radiation for e.g. free-free be given by something like $R{ff} = \int{\lambda} I_{ff}(\lambda) d\lambda$ ?

[jwreep]

The f-f emissivity according to technical report 11 is

CHIANTI's technical report #9 addresses the radiative losses. They neglect two-photon's contributions to the total losses. The equation for free-free (2) doesn't appear to be a direct integration of the emissivity.

Note the exponential factor has disappeared. This is perhaps why a direct numerical integration is not producing the correct result. I have no idea where the exponential went, though!

Edit: The Gaunt factor depends on wavelength as well. Sutherland 1998 gives full derivations of both quantities. See section 2.4.

[jwreep]

Tentatively fixed. Still needs a few unit tests to be added, but the losses compare well to IDL.

This uses the same implementation as the IDL formulation, which is based on Mewe et al 1986 for f-b losses and Sutherland 1998 for f-f. Two-photon is neglected.

[jwreep]

@wtbarnes, this should be good to go!

I've followed the implementation in the IDL version, which doesn't integrate the f-f and f-b functions directly, but implements separate functions for radiative losses. Two-photon emission is neglected in this.

I pinned the numpy version < 2 for now.

One small note: the implementation of the f-b Gaunt factor in IDL uses a polynomial approximation to an infinite sum (Eq. 14 of Mewe et al 1986). Scipy has the appropriate function for the analytic solution, so I've used that instead. The plot below compares the approximation to the analytic solution. It should be a negligible difference but could cause a small discrepancy with IDL's output. We could make the approximation an option, if wanted.

[wtbarnes]

Thanks Jeff! This is great and sorry for being so unresponsive on the review front. I will try to get to this either today or early next week.

Re: your last comment, I think using the scipy implementation is better. I did something similar for the Gauss-Laguerre quadrature in the case of the ionization rates.

I'm happy to add an IDL comparison test here so that we can check the output as we do any refactoring. Would you mind pasting here snippet of IDL code you were running when you were comparing against the CHIANTI IDL results?

[jwreep]

No worries -- enjoy Tenerife! I'll add the IDL code when I'm in the office tomorrow.

[jwreep]

Using sun_coronal_2021_chianti.abund:

IDL> ff_rad_loss,temperature,ff 
IDL> fb_rad_loss,temperature,fb      
IDL> rad_loss,temperature,total   
IDL> plot,temperature,total,chars=2,xtit='Temperature',ytit='Loss Rate',/ylog,/xlog,yr=[1e-25,1e-21]
IDL> oplot,temperature,fb,line=1,thick=2     
IDL> oplot,temperature,ff,line=2    
IDL> oplot,temperature,(total-ff-fb),line=3