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.

Connecting 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.

1D Jet Simulations

The set up for these simulations is very simple. The grid is initialized with an ambient object and a wind object. The goal is to study the effect of magnetic fields on the emission behind shocks in a pulsed jet.

Resolution Study

Resolution Study

Simulations involving cooling tend to require a lot of resolution. For this reason, I have been running a simplified version of the 1D Jet simulation to see how resolution affects the dynamics. Here are the initial parameters:

Ambient:
density = 1000 particles/cc
velocity = 0 km/s
temperature = 1000 K
By = 0 uG

Wind:
density = 10000 particles/cc
velocity = 50 km/s
temperature = 1000 K
By = 0 uG

Global Data:
Domain length = 30 AU
Final time = 1.7 computational units (approx 7 cooling times)

Note that these runs feature a stationary ambient, and they do not include MHD or jet perturbations. To get an idea of the form of the shock waves that are produced, here is a lineout of velocity:

This is from a high resolution run (100 base cells + 7 levels of AMR = 12800 effective). You can see that even with this resolution, the front shock is somewhat smooth instead of a sharp jump. Now in order to make some predictions based on the initial parameters, we need some equations…

Important Equations

A strong shock implies that the mach number M >> 1. We will also be assuming an ideal gas, so \gamma = \frac{5}{3}. The following equations can be derived from the usual shock jump conditions:

T_{ps} = \frac{3m_{H}v_{s}^{2}}{16k} \hspace{1 in} (1)

v_{s} = \frac{4}{3}v_{ps} \hspace{1in} (2)

n_{ps} = 4n_{a} \hspace{1in} (3)

Where s stands for shock, ps for post-shock, and a for ambient. T is temperature, v is velocity, n is number density, mH is the mass of a hydrogen atom, and k is Boltzmann's constant.

In order to figure out these post-shock values, we need to somehow relate them to the jet. This is where these equations come in:

v_{ps} = v_{j}(1 + \eta^{-\frac{1}{2}})^{-1} \hspace{1in} (4)

\eta = \frac{n_{j}}{n_{a}} \hspace{1in} (5)

Where j stands for jet, and \eta is basically the density contrast of the jet material to the ambient material. Now we can combine these equations and constants to come up with an equation for the post-shock temperature in terms of jet parameters:

T_{ps} \approx (40.41 \mathrm{\frac{K s^{2}}{km^{2}}}) v_{j}^{2} (1+\eta^{-\frac{1}{2}})^{-2} \hspace{1in} (6)

Where vj is more conveniently in units of km/s. Now we need equations for the characteristic cooling time and cooling length:

t_{cool} = \frac{3kT_{ps}}{n_{ps}\Lambda} \hspace{1in} (7)

l_{cool} = v_{ps} t_{cool} \hspace{1in} (8)

Where \Lambda is the cooling rate in units of cm³erg/s. Using equations (3) and (4), these cooling equations can be rewritten in terms of our known initial parameters, and the post-shock temperature which we will want to have calculated already:

t_{cool} = \frac{3kT_{ps}}{4n_{a}\Lambda} \hspace{1in} (9)

l_{cool} = v_{j}(1 + \eta^{-\frac{1}{2}})^{-1} \frac{3kT_{ps}}{4n_{a}\Lambda} \hspace{1in} (10)

Some Predictions

Based on all these equations, we can make some predictions about the post-shock temperature and what resolution we might want. Using the aforementioned initial parameters, you get the following predictions:

Tps (103 K)	tcool (107 s)	lcool (AU)
58.3123	1.26155	3.20345

A rate of 4.7863x10^-22 was used, which is fairly consistent with Dalgarno McCray. So if we have a domain that is 30 AU and we want 10 cells per cooling length, then we would need a resolution of about 94 cells. That doesn't seem too bad, so I will round up to 100 and start from there.

Results

Here we look at what post-shock temperature the simulation gives as a function of the effective resolution. The effective resolution = (# cells) x 2^{(# AMR levels)}. Each run used a base grid of 100 cells, so the only thing that changed from run to run was the number of AMR levels. For some reason, I could not get it to run with more than 7 levels, so the 25600 run is actually a base grid of 200 cells with 7 levels. For Tps I will use the highest temperature that I get from the simulation. The highest temperature always occurs early on in the simulation and then the post-shock temperature slowly decreases over time, but stays relatively constant at some lower temperature.
Z

Effective Resolution	cells/lcool	Tps (103 K)	Relative Error (%)
100	10.68	34.37	41.06
200	21.36	37.11	36.36
400	42.71	40.21	31.04
800	85.43	43.68	25.09
1600	170.85	46.66	19.98
3200	341.70	49.25	15.54
6400	683.40	51.18	12.23
12800	1366.81	53.20	8.77
25600	2733.61	53.66	7.98

And in graphical form, the table looks like this where I have plotted the expected Tps as the green horizontal line:

So the big question is, why does it appear that we need way more resolution than is predicted from the cooling length? Perhaps it has to do with how the refinement is being triggered. So I will try a resolution of 400, but with no AMR and see if it does any better than the 400 effective resolution which used 2 levels of AMR. With this run I got Tps = 40.59 x 10³ K. So there really was not a significant difference from AMR to fixed grid.

My next idea was to revert back to just DM cooling instead of the new Z cooling. Perhaps, the Z cooling routine is not implemented quite right, and there is some "double counting" for cooling strength or something like this. So here is a data table where I used DM instead of Z. The improvement column is the decrease in relative error from the comparable Z cooling run.
DM

Effective Resolution	cells/lcool	Tps (103 K)	Relative Error (%)
100	10.68	25.71	55.91
200	21.36	31.63	45.76
400	42.71	41.07	29.57
800	85.43	47.65	18.28
1600	170.85	51.37	11.91
3200	341.70	53.30	8.60
6400	683.40	53.98	7.43
12800	1366.81	54.14	7.16
25600	2733.61	58.16	0.26

The DM cooling gets closer to what I expect but the Tps is still a little low. In both cases, Tps seems to be converging to some other value different from my expected value of 58.3123. However, when I do a run with no cooling, I get post-shock temperatures reaching 58.3123, and some even higher…closer to 60.97 with just 100 cells and no AMR. So I decided to take a closer look at the adiabatic case to see if I could figure out the discrepancy. Again, here is the data for no cooling:
None

Effective Resolution	cells/lcool	Tps (103 K)	Relative Error (%)
100	10.68	60.97	4.56
200	21.36	61.09	4.76
400	42.71	61.00	4.61
800	85.43	61.06	4.71
1600	170.85	61.76	5.91
3200	341.70	63.45	8.81
6400	683.40	64.17	10.04
12800	1366.81	65.20	11.81
25600	2733.61	64.16	10.02

You can see that these simulations actually give higher than expected post-shock temperatures. If I used these values as the "correct" values then the post-shock temperatures from the cooling simulations are even further off. Here is a plot that contains all the data:

Shock Structure

I've made some lineouts of density, velocity, pressure, and temperature to get a better understanding of how the simulation is interpreting the shock structure. Every image is taken at the middle of the simulation for the highest resolution run. The solid black line is density, dashed blue is velociity, dotted green is pressure, and dash/dot red is temperature. Everything is in computational units, and it's plotted on a log scale. The structure is more difficult to see in the cooling cases, so I included some magnified images for those two cases.

No Cooling

DM Cooling

Z Cooling

1D Radiative Shock Simulations

These simulations follow Ch. 4 of Delamarter '01 and his treatment of the 1D steady radiative shock problem. The purpose of this problem is to check that the cooling source terms are being handled correctly. Here are some useful links:

​AstroBEAR 1.0 test page

​Delamarter '01 pdf

Equations

The shock jump equations for a stationary shock are used to solve for the initial post-shock values. The post-shock velocity v₂ can be written as:

v_2 = v_1 \frac{(\gamma - 1)M^2 + 2}{(\gamma + 1)M^2}

where v₁ is the ambient velocity, and M is the ambient mach number. Remember that the mach number M = v₁/c where c is the ambient sound speed, and c = \sqrt{\frac{\gamma P_{1}}{\rho_{1}}} where P₁ is ambient pressure, and \rho_{1} is ambient density. The post-shock density and pressure (\rho_{2} and P₂ respectively) can be found by using mass flux and momentum flux conservation across the shock:

\rho_1 v_1 = \rho_2 v_2

\rho_1 v_1 ^2 + P_1 = \rho_2 v_2 ^2 + P_2

These post-shock values become the boundary conditions for the fluid equations in the cooling region:

\rho v = \rho_2 v_2

\rho v ^2 + P = \rho_2 v_2 ^2 + P_2

\frac{d}{dx}(\frac{1}{2}\rho v ^3 + \frac{\gamma}{\gamma - 1} P v) = -\Lambda

where \rho, v, and P represent the density, velocity, and pressure in the cooling region as functions of x, and \Lambda is the cooling rate.

Initial Parameters

n₁ = 60 particles/cc

v₁ = 10⁷ cm/s

T₁ = 10⁴ K

For analytic cooling,

\Lambda = n^{2} \alpha T^{\beta}

where n and T and the number density and temperature in the post-shock region respectively.

For these simulations,

\beta = 2

\alpha = 1.23786 * 10^-34 erg*cm³/s/K²

cell length = 2.5 x 10¹⁵ cm

problem domain = 400 cells

final time ~ 4000 years

Results

Results

A run with no cooling was done to make sure the post-shock values were correct. With no cooling (adiabatic), the hydrodynamic quantities should jump discontinuously at the shock, and then remain constant while the shock itself remains stationary. Only the first 100 cells of the simulation are shown here.

Adiabatic
Pressure Movie	Temperature Movie
Density Movie	Velocity Movie

Then, a run with the aforementioned analytic cooling parameters was done. The hydrodynamic quantities will change discontinuously at the shock as before, but then change continuously according to the cooling rate. The ambient temperature is set to be the floor temperature, so that no cooling occurs at or below this temperature. Once the post-shock temperature reaches this floor temperature, the hydrodynamic quantities will remain constant. This is called the quiescent region. Only the first 100 cells of the simulation are shown here.

Analytic Cooling
Pressure Movie	Temperature Movie
Density Movie	Velocity Movie

For Dalgarno-McCray (DM) cooling, some of the initial parameters had to be changed. The initial velocity used in the previous simulations would give rise to a post-shock temperature in the unstable region of the DM curve. This led to the radiative instability, and thus unsteady shocks. To avoid this, the initial velocity was lowered to 80 km/s. Also, the cooling rate from the DM curve is much higher than that of the analytic form. In order to see the cooling region, the grid had to be shrunk, resulting in a cell length of 10¹³ cm.

DM Cooling
Temperature Movie

For Non-equilibrium (NEQ) cooling, the length scale had to be adjusted. The shock velocity is the same as the DM case, but the cell length is now 7.75e11 cm.

NEQ Cooling
Temperature Movie

Z Cooling has the same initial parameters and scaling as the NEQ case.

Z Cooling
Temperature Movie

Here is a fun plot that shows the temperature profiles for DM, NEQ, and Z cooling. This gives you a good idea of the temperature ranges that each type of cooling works in. Also, you can see the differences in the cooling length which is why I would normally change the length scale depending on the type of cooling.

Including MHD

Making these simulations work with MHD was no simple task. Having a magnetic field adds another jump condition, and it makes the equations within the cooling region more complicated. Here are what the equations look like now: