atomic77 data gives anomalous recombination cooling results for hydrogen and others.

[Higginbottom]

This behaviour was first seen by NSH when running a loop of thin shell models with a flat power law spectrum to compare with cloudy.
The recombination cooling vs IP curve has a step change, at IP about -3 > -2.

Further investigation reveals that below the threshold, hydrogen gives zero recombination cooling, and above the threshold it dominates the recombination cooling.
Hydrogen is mostly ionised in both cases.
Further investigation revealed that switching from standard77 to standard73 fixes the problem.
Presumably this means it is a numerical issue with the integration scheme in integ_fb with the new much larger frequency range over which the cross sections are defined.

Is there any way of attaching a file to this? If so I'll attach a .pf file that reproduces this behaviour. If not - anyone who wants it just ask and I'll mail it..

[Higginbottom]

Further testing shows that if, by hand, I restrict the range of the integrals in xinteg_fb, then the answers for fb emissivity return to the values of the standard73.
Thus, this is a numerical issue of doing an integration from (for example) 3e15 to 3e18 Hz where the integrand is effectively zero for most of the range.
We need a scheme to prevent qromb going wrong…

[kslong]

Nick,

I'm glad your are pursuing this.

We should take a look at what the actual integrals are. If we are doing a piece-wise power law for the cross section we may be able to calculate the results directly.

The alternative is to find an upper frequency bound to carry out the integration over. My guess is that this could be estimated from a hydrogenic approximation.

Knox

[Higginbottom]

So, this is a problem I've encountered before - when computing the ground state recombination rate via the milne relation. It comes about because there is a factor of e^(-hnu/kt) in the integral.
So, indépendant of what the cross section is doing, this totally swamps the integration, and if the plasma is cold enough that there is h(nu-nu0) << kT at some critical point, the qromb integration fails because the integrand is effectively zero for a significant part of the integral.

I fixed it very simply, by the following code

  ntmin = ion[nion].ntop_ground;            get the relevant transition         
  fb_xtop = &phot_top[ntmin];             
  fthresh = fb_xtop->freq[0];                     get the threshold frequency (i.e. ionisation potential of the ion)
  fmax = fb_xtop->freq[fb_xtop->np - 1];     get the maximum frequency
  dnu = 100.0 * (fbt / H_OVER_K);           set the variable such that we know that exp (-h (fthresh+dnu)/kT) is real
  if (fthresh + dnu < fmax)                    if the fmax we have from the cross section is too high 
{
  fmax = fthresh + dnu;                    set fmax to a value that will not underflow the integral
}
  rate = qromb (fb_topbase_partial, fthresh, fmax, 1e-5);   do the integration
}

This solution could be used here - it is independent of what cross section we are integrating over - even if it is dropping as 1/nu^3, it is the exponential that causes the problem…

[Higginbottom]

Another little note:

I was slightly confused that the integration which computes the number of recombinations was working, but the emissivity was failing. I now realise that this was because the integrand of the emissivity is multiplied by the energy released by the recombination. So, for low temperatures, the only point returning a non zero value was the intrgrand evaluated at the threshold frequency. Of course, this is fine for the number, but since the energy released by a recombination of an electron with just the threshold energy is zero, so the only non zero value, gets multiplied by zero!
Now going to test the solution above…

[Higginbottom]

I think I've fixed this particular issue, see figure below

The free bound cooling curve (free crosses) no longer has a discontinuity.
Now, why is cloudy so much higher than python - probably due to ionisation state.

[jhmatthews]

Commit 86f1d17 closes #62