What is this nonequilibrium cooling nonsense?

Jonathan and I have met a few times now to discuss nonequilibrium cooling since he is implementing the molecular hydrogen routines and variable gamma. While trying to understand how NEQCooling modifies DMCooling, we realized that it was not being modified correctly. First, let me define some of the main parts in the different types of cooling. 

Hex_eq = equilibrium hydrogen excitation
Hex_neq = nonequilibrium hydrogen excitation
DM_em = metal excitation via electron/metal collisions (from DM)
DM_Hm = metal excitation via hydrogen/metal collisions (from DM)
Z_em = metal excitation via electron/metal cololisions (from Pat)

Now, here is how the different types of cooling are related to each other: 

DM = Hex_eq + DM_em + DM_Hm
NEQ = DM - Hex_eq + Hex_neq
Z = NEQ - DM_em + Z_em

NEQ

For the NEQ case, I think there was some confusion in the past and this was not handled correctly. Hydrogen excitation was being double counted in a weird way. There was a routine for an equilibrium ionization fraction table, which should have been used to get Hex_eq, but it was not used and I could not find a table. To fix this, I implemented a function to calculate the equilibrium ionization fraction: 

! solves for equilibrium ionization fraction of hydrogen (derived from Saha equation)
FUNCTION EqIonFrac(Temp, nH)
   REAL(KIND=qPREC) :: EqIonFrac, Temp, nH, b

   b = (2d0*Pi*emass*Boltzmann*Temp/planck**2d0)**(1.5d0)*EXP(-13.6d0*ev/(Boltzmann*Temp))/nH
   IF (b>10d10) THEN
      EqIonFrac = 1d0
   ELSE
      EqIonFrac = (-b + SQRT(b**2d0 + 4d0*b))/2d0
   END IF

END FUNCTION

The equation for EqIonFrac is quadratic, and b is used with the quadratic formula to get EqIonFrac. The IF statement is in there because I learned a valuable lesson about limits and large numbers. You can see as T becomes large, the exponential in b goes to 1 and b scales as T1.5. So b also becomes large as T becomes large. Due to rounding, the computer takes the square root part to be approximately b as b becomes large, and thus EqIonFrac to be zero. However, if you take the limit of EqIonFrac as b goes to infinity you find that it is 1 (which makes more sense physically). I found that 1010 is a good limit to use for b in that IF statement.


Z

The implementation of Z is slightly more complicated. Pat's cooling tables have a relatively limited range. Outside of this range, Z should revert back to NEQ. This transition needs to be smooth, so DM_em is not completely subtracted from NEQ and Z_em is not completely added to NEQ. These two components are weighted depending on the temperature.

Another issue is that in the code, we do not have a good way of separating DM_em and DM_Hm from the DM curve. So for Z cooling, Z_em is actually replacing DM_em + DM_Hm. This is okay as long as DM_Hm is small in the range of Pat's table. According to DM, DM_Hm is only important when the ionization fraction is low (approximately < 0.001). In Pat's table, the ionization fraction never goes below 0.01, so DM_Hm is not important and can safely be ignored.

Below is an example of how the fixes affected the temperature profile behind a particular shock. DM is black, NEQ is blue, Z is red…the old version is dashed, and the new version is solid. vs = 30 km/s, namb = 100 cm-3, B = 0 (hydro run)

Attachments (1)

Download all attachments as: .zip

Comments

No comments.