bf heating NAN: macro-only problem

[jhmatthews]

Version discovered in: 77a_dev
Description: heat_photo is NAN in some cells in the macro atom CV models.
How to Reproduce: run ~/Dropbox/Python/bug_pf_files/sv_macronan.pf

Discussion:

In the recent CV models with macro atoms I'd discovered that bf heating was being calculated as not a number in macro atom runs. I've tracked the problem to the function q_recomb, which works out the collisional recombination coefficient- i.e. three body recombination. This is called by the routine macro_bf_heating, and is the bf process other than photoionization which can heat the plasma.

The piece of code for q_recomb is as follows

/* q_recomb. This returns the collisional recombination co-efficient
Calculated from inverse of q_ioniz.
*/

double
q_recomb (cont_ptr, electron_temperature)
     struct topbase_phot *cont_ptr;
     double electron_temperature;
{
  double coeff;
  double root_etemp;
  double q_ioniz ();

  root_etemp = sqrt (electron_temperature);
  coeff = 2.07e-16 / (root_etemp * root_etemp * root_etemp);
  coeff *= exp (cont_ptr->freq[0] * H_OVER_K / electron_temperature);
  coeff *= config[cont_ptr->nlev].g / config[cont_ptr->uplev].g;
  coeff *= q_ioniz (cont_ptr, electron_temperature);

  return (coeff);
}

The reason we get NANs is due to the exponential- here we have electron temperature of 880K and a frequency of 1.31e16, which means we try to calculate exp(715), which fails.

It occurs to me that I'm not sure I understand the above equation for q_recomb, which is:

and haven't as yet been able to find this in the literature or my version of Mihalas.

[jhmatthews]

So when we consider collisional excitation and de-excitation we get the relationship

q12 / q21 = g_2 / g_1 * exp (-hnu/kT)

so I think it's fairly obvious that you might expect the form

if your upper level is a continuum level, but I don't quite understand where the T_e dependence comes in- presumably because your continuum state is not really a quantised state but a Maxwellian velocity distribution.

I will discuss this with @lazygun37 and try and get it straight and check we understand what is being done by the code. Note that this also fits in with what Nick has been doing regarding collisional ionization.

In terms of a bugfix- firstly I should check why temperatures are getting that low, if it is physical and whether there should be checks to ensure it doesn't happen- note there are checks in calc_te but I feel like there should also be some here. Do we then just need a floor for this number? Presumably the temperature can never be this low if the collisional recombination heating is ever significant, as for this recombination heating to be significant you would need a highly ionized gas and thus a significant amount of heating would have come from photoionization meaning the temperature would have to be a lot higher.

[ssim]

The relation between q_recomb and q_ioniz just follows from considering TE: these are thermodynamic inverse processes and so you require that, in TE, n_l n_e q_ioniz = n_u n_e^2 q_recomb - so that q_recomb = q_ioniz (n_l / n_u / n_e)_TE. Then it's just the Saha-Boltzmann equation to evaluate n_l / n_u /n_e in TE:

n_l / n_u /n_e = (h^2 / 2 /pi /m_e / k / T)^3/2 g_l / 2 / g_u * exp( delta_E / k T)

The T is clearly the kinetic temperature (not radiative) since these are purely collisional processes.

[jhmatthews]

Ok, thanks Stuart.

I've been looking into why the wind gets so cool in the outer regions, and of course it is due to adiabatic cooling. I was also getting very confused and I think the thing I wrote earlier is completely wrong.

This opens up a a lot of questions/issues:

1) The numerical fix. The expression for q_recomb is (where u = h nu / kT, and the h nu is just the energy gap between the continuum state and the level in question.

 q_recomb = 2.07e-16 /  (T**3/2)   *  gl / gu   *  exp(u)  *  q_ioniz              (1)

where

q_ioniz = 1.55e13 /  (T**1/2)  * gaunt * sigma * exp(-u) / u                         (2)

so I think we can just eliminate the need for q_recomb to call q_ioniz and therefore eliminate the exponential. Combining 1 and 2 gives:

q_recomb = 3.2085e-3   *   (T**-2)   * gaunt * sigma  * gl/gu    / u

         = 3.2085e-3  * gaunt * sigma  * gl/gu  * (k/h nu)  / T

Which I think I've done right. This would be quicker and more numerically stable and is quite an elegant solution. It also makes me wonder if we can do similar things elsewhere.

2) Should adiabatic cooling be on / do we trust results in regimes where it is driving down the temperature?

This is possibly the most important question.

3) I have not yet implemented adiabatic cooling in macro atoms.
This means that we are actually violating conservation of energy in these regions of the wind, because there is no kpkt destruction process for adiabatic cooling yet, although I have started working on this and it won't be difficult to sort out. The question that may be slightly tricky is sorting it out properly so that it can also serve as a net input of kpkts, something I'm less keen on, if we include an adiabatic heating term.

4) Inconsitent temperature checks and floors

Here, the temperature that is used to normalize the MC macro atom estimators is that which was calculated earlier (882K), which is less than TMIN of 2000K (which should not be allowed?), and also is reported different in py_wind as 1300. We need to check how this is done and why.

[kslong]

James,

I am in favor of maintaining adiabatic cooling as an option. I don't really see what this would not be operative in real CVs. I would like us to be able to turn it off, since I think it's an interesting question how big an effect it is.

Note that if you want to worry about something. You might worry about the question of whether the flow times in CV winds are so short that the ionization is frozen in. I know some of the X-ray folks have concluded that it might be.

Knox

Knox S. Long
Space Telescope Science Institute
1-410-338-4862

[jhmatthews]

I've implemented my suggested change to q_recomb in my fork and I'm testing. I will commit that to dev when I'm happy it works and gives the answers it should.

I suggest that my first job then is to sort out adiabatic cooling as a kpkt destructions process, which is pretty easy, I just need to work out how to break out of the photon loop properly. We can then have a proper think about what to do in these regions.

I guess my point was not so much 'we should turn it off' -- more that in these regimes where it is very important and driving the temperature down do we trust our treatment (I'm pretty sure @ssim has raised this point). I was perhaps advocating something more like having it on, but cutting off the model once we hit this region, but I guess that has no real benefit. Either it's not an important region in which case cutting it off has no benefit, or it is in which case cutting it off is wrong.

I feel a little uncomfortable about modeling regions with supposed temperatures of 800K. I will have a think about this, and the flow times point and try to articulate this a little better.

[kslong]

Hi james,

Our treatment is I believe correct. The problem in the case of the hydro models was that we had converging flows which would drive shocks. For cooling, which is what occurs in an expanding wind, I do not believe there is a problem.

To zeroth order, the cooling time scale is just the thermal energy content / (total rate for energy losses from the gas). Note that since we balance heating, one could equivalently put heating into the denominator.

If the cooling (or heating) time is short compared to flow time then there could be no problem.

All this is clearly beyond what we would be concerned about in a paper, but it's useful to know about the possibility.

Cheers,
Knox

Knox S. Long
Space Telescope Science Institute
1-410-338-4862

[jhmatthews]

I've tested and my numerical fix seems to work fine, and results in the models converging to 96% with no NANs. There is some difference in wind temperature far out and I'm looking into why.

Calculations seems to suggest that cooling timescales and flow timescales are rather similar in fast moving, cool, expanding regions in the outer regions of the wind where adiabatic cooling is dominant

t_cooling ~ NkT / cool_adiab ~  n_h kT / cool_adiab ~ 500s

t_flow ~ z_cell / v_z ~ 2e11 / 2.8e8 ~ 700s

These are the worst cells as they have the lowest energy content, are fast moving, and have the most significant adiabatic cooling. This does suggest that we need to think about this though- perhaps I should move this discussion elsewhere now it is diverging somewhat from the issue (no pun intended).

[jhmatthews]

This specific bug has been fixed in pull request #65

The question of flow timescales is very much open.

Adiabatic cooling implementation is discussed in pull request #66