BackLinksMenu

Z Cooling

The Physics of Z Cooling

Z cooling is an energy loss due to forbidden line emission. The routines use Pat Hartigan's cooling table Phcooling.tab. The table does not include recombination or permitted line emission.

Implementation into AstroBEAR 2.0

PHcooling.tab

Understanding the format of this table is crucial to interpolating correctly. There are a couple of properties that this table must have:

1.) It must be ordered as log(ne), T, log(X), (9 species cooling fractions), cooling rate. Where ne is the electron density, T is temperature, X is ionization fraction, and the 9 species cooling fractions are the fractions of how much each species contributes to the cooling rate. Note that here log refers to the base 10 logarithm. The 9 species are ordered as follows: OI, OII, NI, NII, SII, FeII, SiII, MgII, CII.

2.) The step sizes of log(ne), T, and log(X) must remain constant. In other words, each column is evenly spaced.

Also, it is important to note that the cooling rate has units of erg * cm^-3 * nH^-2. So before the rate can be added to the change in energy, dqdt(iE), it must be multiplied by nH² where nH is the number density of hydrogen.

Zcooling.tab

PHcooling.tab contains a lot of data that is not needed. Reading in this data one line at a time is inefficient, especially when astrobear is being run on multiple processors. This is why a separate script called ph2zformat.s^* is used to reformat PHcooling.tab into Zcooling.tab before running astrobear. Remember that PHcooling.tab still requires the aforementioned properties in order to be reformatted correctly. Zcooling.tab has the following format:

nDensities
nTemps
nXs
lognemin
lognemax
tempmin
tempmax
logxmin
logxmax
coolingrate(1,1,1)
coolingrate(2,1,1)
...
coolingrate(nDensities,nTemps,nxs)

Where nDensities, nTemps, and nXs are the number of unique values for log(ne), T, and log(X) respectively. So if these parameters are equal to 40, 30, and 21, then there are a total of 25,200 cooling rates (40*30*21 = 25200). lognemin is the minimum value of the log(ne) column in PHcooling.tab, and the other minimum and maximum values are defined in the same way. Now the data does not need to be read in one line at a time for 25,200 or so values. The entire array of cooling rates can be read in with a single READ statement. No iteration or nested DO loops required.

^*The only thing that ph2zformat.s assumes is that it is being run from the TABLES directory, and that the file you want to reformat is named PHcooling.tab.

Preparing for the Interpolation

There are several steps that need to occur before the actual interpolation. When Zcooling.tab is read in, the cooling rates are assigned to a properly allocated 3D array called Zcoolingtab(i,j,k). The indices i, j, k of Zcoolingtab are integers ranging from 1:nDensities, 1:nTemps, and 1:nXs respectively. These indices can also be thought of as coordinates. Before the cooling rate can be calculated, values for ne, T, and X are calculated. The electron density ne and ionization fraction X are converted to log(ne) and log(X). Then these values for log(ne), T, and log(X) can be scaled to the integer ranges of the coordinates i, j, and k. This is where the minimum and maximum values in Zcooling.tab are used.

In order to do a 3D interpolation, we need to find the 8 closest points to our newly scaled point of interest. You can think of these 8 points as the vertices of a cube or rectangular prism, and our point of interest is somewhere inside this box. First, we find the 8 closest coordinates. Since the coordinates i, j, k are all integers, the 8 closest coordinates are easily found by using the FLOOR and CEILING functions. Using these functions, you get a "low" and "high" value for each coordinate. Since there are 3 coordinates, 2*2*2 = 8 unique sets of coordinates.

The next step is to match these 8 closest coordinates with their corresponding cooling rates. This is done easily with 8 simple assignments. For example, rate1 = Zcoolingtab(lognelow, templow, logxlow) would be the rate that corresponds to the lower bound coordinates. Now we have the 8 closest points.

The final step before interpolating is to do another rescaling of the coordinates. We take the 8 closest coordinates to be the coordinates of a unit cube aligned with the origin. The coordinates of our point of interest are then scaled accordingly. This makes the interpolation very simple. Remember that in all these rescaling steps, nothing happens to the cooling rate because this is the value that we are interested in.

Trilinear Interpolation

Now we have 8 points that are positioned on the vertices of a unit cube aligned with the origin. Each point has a unique cooling rate. Our point of interest is inside this cube, but we do not yet know what its cooling rate should be. Trilinear interpolation is the method used to find this rate. The formula used is from Paul Bourke's website: ​http://paulbourke.net/miscellaneous/interpolation/. In cooling.f90, this is done on 8 lines to avoid having one lengthy, confusing line…

    ZCoolingRate = rate1*(1d0-lognein)*(1d0-tempin)*(1d0-logxin)
    ZCoolingRate = ZCoolingRate + rate2*lognein*(1d0-tempin)*(1d0-logxin)
    ZCoolingRate = ZCoolingRate + rate3*(1d0-lognein)*tempin*(1d0-logxin)
    ZCoolingRate = ZCoolingRate + rate4*(1d0-lognein)*(1d0-tempin)*logxin
    ZCoolingRate = ZCoolingRate + rate5*lognein*(1d0-tempin)*logxin
    ZCoolingRate = ZCoolingRate + rate6*(1d0-lognein)*tempin*logxin
    ZCoolingRate = ZCoolingRate + rate7*lognein*tempin*(1d0-logxin)
    ZCoolingRate = ZCoolingRate + rate8*lognein*tempin*logxin

Where rate1,rate2,…,rate8 are assigned as aforementioned, and the "in" values represent the coordinates of our point of interest inside the unit cube.

Connection to NEQ Cooling

One component of the cooling rate calculated by NEQ cooling is called metal_loss. This is the energy loss due to metal excitation, and it is based on the ionization fraction and DM cooling curve. Z cooling passes the parameter iCooling to the subroutine Cool_Derivatives. Now, if iCooling==ZCool, metal_loss will be subtracted from the cooling rate, and it is going to be replaced by a different rate.

Back in the main cooling file cooling.f90, ZCoolingRate calculates the rate due to forbidden line emission. It is based on temperature, electron density, and ionization fraction. In the ZCoolingStrength subroutine, the metal_loss and ZCoolingRate are combined to give a weighted cooling rate. This is done by a parabolic weight function which I have called Zweight. The purpose of this weighted rate is to make the transition from NEQ cooling to Z cooling smooth.

The table for ZCoolingRate is currently from 2000 to 16,500 K. Regular NEQ cooling will only cool effectively to about 8000 K or so. So you need to use Z cooling if you need to look at cooling below about 1e4 K.

---

Z cooling is related to ticket #224, milestone Implementing Multi-Physics, and it was tested with the 1D radiative shock problem.