Material prepared for the 52nd Saas-Fee Advanced Course on "The circum-galactic medium across cosmic time: an observational and modelling challenge"¶

by Alejandro Benitez-Llambay & Michele Fumagalli¶

Contact for questions or comments: alejandro.benitezllambay@unimib.it

Calculating and Understanding Cooling Functions with CLOUDY¶

Excercise 1 - Cooling by collisional ionization in a constant temperature hydrogen cloud¶

This exercise aims to become familiar with CLOUDY by studying the cooling rate arising from the collisional ionization of a hydrogen gas cloud as a function of density and temperature.

Tasks¶

1) Write a CLOUDY input parameter file to calculate the cooling rate of a pure hydrogen gas cloud in thermal equilibrium as a function of the cloud temperature and density.

Sample the temperature in the range, $4 \le \rm log_{10} (T / K) \le 8$, using a grid spacing $\rm \Delta log(T)=0.05$. Sample the density in the range $-4 \le \rm log_{10} (n_{H} / cm^{-3}) \le 0$, using a grid spacing $\rm \Delta log (n_{H})=1$.

Although not strictly necessary, running CLOUDY on the desired grid is convenient using the grid command. Check the documentation for details (Hazy1.pdf, Chapter 18).

2) Plot the cooling rate for each density as a function of temperature. What features do you see in the cooling curve? Discuss the physics behind the features observed in the plot

3) Does the cooling rate depend on density? Consider now the so-called cooling function $\Lambda(T)=\mathcal{C}(n_{\rm H}, T)/n_{H}^2$. Does the cooling function depend on density? Does the shape of the cooling curve depend on density?

Example Solution¶

CLOUDY parameter file for this excersice

hden -2 vary //Density of hydrogen in atoms per cubic centimetre
grid, range from -4 to 0 with 1.0 dex steps
constant temperature 4 vary // Temperature of the cloud
grid, range from 4 to 8 with 0.05 dex steps
abundances PRIM //We set the abundances to primordial. Equivalent to this would be to change PRIM by primordial.abn"
abundances all -20 // Set the abundance of all but hydrogen to 10^-20
set dr 0 //Thickness of the shell is 0
stop zone 1 // Stop after calculating one zone.
save cooling "cooling_only.cool" no hash //Save the desired data
save grid "cooling_only.grd" no hash
save overview "cooling_only.ovr" no hash
In [1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

nHmin = -4.0
nHmax = 0.0
nHsteps = 1.0

Tmin = 4.0
Tmax = 8.0
Tsteps = 0.05

NnHbins = 0
nH = nHmin
nHbins = []
while(nH < nHmax):
    nH = nHmin + nHsteps * NnHbins
    nHbins.append(nH)
    NnHbins += 1

NTbins = 0
T = Tmin
Tbins = []
while(T < Tmax):
    T = Tmin + Tsteps * NTbins
    Tbins.append(T)
    NTbins += 1

cooling_rate = np.zeros([NnHbins, NTbins])
heating_rate = np.zeros([NnHbins, NTbins])

with open('tables/cooling_only_H/cooling_only.cool','r') as input_file:
    lines = input_file.readlines()
    for i in range(NnHbins):
        for j in range(NTbins):
            cooling_rate[i,j] = float(lines[i*NTbins+j+1].split()[3])
          
        
fig = plt.figure(1, figsize=(10,5), tight_layout=True)
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)


for i in range(len(nHbins)):
    c = (nHbins[i]-nHmin)/(nHmax-nHmin)
    color = cm.jet(c)                
    ax1.plot(Tbins, np.log10(cooling_rate[i,:]), label=r'${%.1f}$'% nHbins[i], color=color)
    ax2.plot(Tbins, np.log10(cooling_rate[i,:])-2.0 * nHbins[i], label=r'${%.1f}$'% nHbins[i], color=color)

    
ax1.set_xlabel(r'$\rm log_{10} \left ( T / K \right ) $')
ax1.set_ylabel(r'$\rm log_{10} \left ( \mathcal{C} \right ) \ [erg/cm^3/s] $')
ax1.legend(ncol=2, title=r'$\rm log(n_{H}/cm^{-3})$', loc='lower right')
ax2.legend(ncol=2, title=r'$\rm log(n_{H}/cm^{-3})$', loc='lower right') 

ax2.set_xlabel(r'$\rm log_{10} \left ( T / K \right ) $')
ax2.set_ylabel(r'$\rm log_{10} \Lambda = log_{10} \left ( \mathcal{C} / nH^2 \right ) $')

fig.suptitle('Radiative cooling of a Hydrogen+Helium mixure')
plt.show()

Excercise 2 - Cooling by collisional ionization in a constant temperaturcmhydrogen + helium cloud¶

Tasks¶

1) Modify the CLOUDY input file used to solve the previous exercise to calculate the cooling rate arising from the collisional ionization of primordial hydrogen + helium(*) mixture. Can you identify the new features in the cooling curve arising from Helium? Does the cooling function, $\Lambda$, depend on density?

2) On top of the total cooling function, $\Lambda$, plot the contribution of Hydrogen and Helium. What determines the peak in the cooling rate at around $3\times 10^4 \ K$ and $10^5 \ K$?

(*) Note that in CLOUDY, the abundances are always specified by number relative to hydrogen, not by mass or relative to silicon. Therefore, the abundance of helium should be set using the command element abundance helium -1 (see hazy1.pdf, Chapter 7)

NOTE: Although we propose to use CLOUDY to calculate the cooling function in the two previous exercises, the cooling function can be calculated numerically using fitting functions and assuming ionization equilibrium. The same applies to gas subject to an external ionizing radiation field.

Example solution¶

CLOUDY parameter file for this excersice
hden -2 vary //Density of hydrogen in atoms per cubic centimetre
grid, range from -4 to 0 with 1.0 dex steps
constant temperature 4 vary // Temperature of the cloud
grid, range from 4 to 8 with 0.05 dex steps
abundances PRIM #We set the abundances to primordial. Equivalent to this would be to change PRIM by primordial.abn"
abundances all -20 // Set the abundance of all but hydrogen to 10^-20
element abundance helium -1
set dr 0 //Thickness of the shell is 0
stop zone 1 // Stop after calculating one zone.
save cooling each "cooling.cool" no hash //Save the desired data
save grid "cooling.grd" no hash
save overview "cooling.ovr" no hash
In [3]:
import numpy as np
import matplotlib.pyplot as plt

nHmin = -4.0
nHmax = 0.0
nHsteps = 1.0

Tmin = 4.0
Tmax = 8.0
Tsteps = 0.05

NnHbins = 0
nH = nHmin
nHbins = []
while(nH < nHmax):
    nH = nHmin + nHsteps * NnHbins
    nHbins.append(nH)
    NnHbins += 1

NTbins = 0
T = Tmin
Tbins = []
while(T < Tmax):
    T = Tmin + Tsteps * NTbins
    Tbins.append(T)
    NTbins += 1

cooling_rate_total = np.zeros([NnHbins, NTbins])
cooling_rate_H = np.zeros([NnHbins, NTbins])
cooling_rate_He = np.zeros([NnHbins, NTbins])

with open('tables/cooling_only_H_HE/cooling.cool','r') as input_file:
    lines = input_file.readlines()
    for i in range(NnHbins):
        for j in range(NTbins):
            cooling_rate_total[i,j] = float(lines[i*NTbins+j+1].split()[2])
            cooling_rate_H[i,j] = float(lines[i*NTbins+j+1].split()[3])
            cooling_rate_He[i,j] = float(lines[i*NTbins+j+1].split()[4])
        
fig = plt.figure(1, figsize=(15,5), tight_layout=True)
ax1 = fig.add_subplot(131)
ax2 = fig.add_subplot(132)
ax3 = fig.add_subplot(133)


for i in range(len(nHbins)):
    c = (nHbins[i]-nHmin)/(nHmax-nHmin)
    color = cm.jet(c)            
    ax1.plot(Tbins, np.log10(cooling_rate_total[i,:]), label=r'${%.1f}$'% nHbins[i], color=color)
    ax2.plot(Tbins, np.log10(cooling_rate_total[i,:])-2.0 * nHbins[i], label=r'${%.1f}$'% nHbins[i], color=color)

#Any curve will work, as the cooling function is independent of density

ax2.set_ylim(-25,-21.5)

ax3.plot(Tbins, np.log10(cooling_rate_H[0,:])-2.0 * nHbins[0], label='Hydrogen', ls='-.')
ax3.plot(Tbins, np.log10(cooling_rate_He[0,:])-2.0 * nHbins[0], label='Helium', ls='--')
ax3.plot(Tbins, np.log10(cooling_rate_total[0,:])-2.0 * nHbins[0], label='Hydrogen', ls='-', color='k', lw=2)
ax3.set_ylim(ax2.get_ylim())

ax1.set_xlabel(r'$\rm log_{10} \left ( T / K \right ) $')
ax1.set_ylabel(r'$\rm log_{10} \left ( \mathcal{C} \right ) \ [erg/cm^3/s] $')
ax1.legend(ncol=3, title=r'$\rm log(n_{H}/cm^{-3})$', loc='lower right')
ax2.legend(ncol=3, title=r'$\rm log(n_{H}/cm^{-3})$', loc='lower right') 

ax2.set_xlabel(r'$\rm log_{10} \left ( T / K \right ) $')
ax2.set_ylabel(r'$\rm log_{10} \Lambda = log_{10} \left ( \mathcal{C} / nH^2 \right ) $')

ax3.set_xlabel(r'$\rm log_{10} \left ( T / K \right ) $')
ax3.set_ylabel(r'$\rm log_{10} \Lambda = log_{10} \left ( \mathcal{C} / nH^2 \right ) $')

fig.suptitle('Radiative cooling of a Hydrogen+Helium mixure')
plt.show()

Excercise 3 - Cooling and heating by collisional ionization in a constant temperature hydrogen + helium cloud, subject to photoionization by an external radiation field¶

Tasks¶

1) Modify the CLOUDY input file used to solve the previous exercise to calculate the cooling and heating rates arising from the collisional ionization of a primordial hydrogen + helium mixture subject to the external present-day Haardt \& Madau (2012)(*) ultraviolet background radiation.

2) Plot, as usual, the cooling rate as a function of temperature for each density. Can you identify the new features in the cooling curve arising from the presence of the external UVB? Does the cooling function, $\Lambda$, depend now on density? If so, why?

3) Similarly to the cooling rate, plot the heating rate as a function of temperature for each density.

4) Study the regimes where cooling and heating dominates. Why does cooling dominate over heating at high densities? Why heating dominates over cooling at low densities?

5) What would happen if a gas cloud is exposed to the external UVB? How will the final temperature of the cloud depend on its initial temperature?

(*) You can specify the spectrum of the incoming radiation to be that of the HM2012 by using the command: table hm12 redshift 0.

Example solution¶

CLOUDY parameter file for this excersice

table hm12 redshift 0
hden -2 vary //Density of hydrogen in atoms per cubic centimetre
grid, range from -6 to 0 with 1.0 dex steps
constant temperature 4 vary // Temperature of the cloud
grid, range from 4 to 8 with 0.05 dex steps
abundances PRIM //We set the abundances to primordial. Equivalent to this would be to change PRIM by primordial.abn"
abundances all -20 // Set the abundance of all but hydrogen to 10^-20
element abundance helium -1
set dr 0 //Thickness of the shell is 0
stop zone 1 // Stop after calculating one zone.
save cooling "cooling.cool" no hash //Save the desired data
save grid "cooling.grd" no hash
save overview "cooling.ovr" no hash
In [4]:
import numpy as np
import matplotlib.pyplot as plt

nHmin = -6.0
nHmax = 0.0
nHsteps = 1.0

Tmin = 4.0
Tmax = 8.0
Tsteps = 0.05

NnHbins = 0
nH = nHmin
nHbins = []
while(nH < nHmax):
    nH = nHmin + nHsteps * NnHbins
    nHbins.append(nH)
    NnHbins += 1

NTbins = 0
T = Tmin
Tbins = []
while(T < Tmax):
    T = Tmin + Tsteps * NTbins
    Tbins.append(T)
    NTbins += 1

cooling_rate = np.zeros([NnHbins, NTbins])
heating_rate = np.zeros([NnHbins, NTbins])

with open('tables/cooling_heating_H_HE/cooling.cool','r') as input_file:
    lines = input_file.readlines()
    for i in range(NnHbins):
        for j in range(NTbins):
            cooling_rate[i,j] = float(lines[i*NTbins+j+1].split()[3])
            heating_rate[i,j] = float(lines[i*NTbins+j+1].split()[2])
          
        
fig = plt.figure(1, figsize=(15,5), tight_layout=True)
ax1 = fig.add_subplot(131)
ax2 = fig.add_subplot(132)
ax3 = fig.add_subplot(133)


for i in range(len(nHbins)):
    c = (nHbins[i]-nHmin)/(nHmax-nHmin)
    color = cm.jet(c)
    ax1.plot(Tbins, np.log10(cooling_rate[i,:]),color=color)
    ax2.plot(Tbins, np.log10(cooling_rate[i,:])-2.0 * nHbins[i], color=color)
    ax1.plot(Tbins, np.log10(heating_rate[i,:]), color=color, ls='--')
    ax2.plot(Tbins, np.log10(heating_rate[i,:])-2.0 * nHbins[i], color=color, ls='--')
    ax3.plot(Tbins, np.log10(abs(cooling_rate[i,:]-heating_rate[i,:]))-2.0 * nHbins[i], label=r'${%.1f}$'% nHbins[i], color=color)

ax2.set_ylim(-25,-21.5)
    
ax1.set_xlabel(r'$\rm log_{10} \left ( T / K \right ) $')
ax1.set_ylabel(r'$\rm log_{10} \left ( \mathcal{C} \right ) \ [erg/cm^3/s] $')
ax3.legend(ncol=1, title=r'$\rm log(n_{H}/cm^{-3})$', loc='lower right') 

ax2.set_xlabel(r'$\rm log_{10} \left ( T / K \right ) $')
ax2.set_ylabel(r'$\rm log_{10} \Lambda = log_{10} \left ( \mathcal{C} / nH^2 \right ) $')

ax3.set_xlabel(r'$\rm log_{10} \left ( T / K \right ) $')
ax3.set_ylabel(r'$\rm log_{10} \Lambda_{\rm eff} = log_{10} \left ( |\mathcal{C}-\mathcal{H}| / nH^2 \right ) $')


fig.suptitle('Radiative cooling and Heating of a Hydrogen+Helium mixure + external HM12 UVB')
plt.show()

Excercise 4 - Cooling and heating by collisional ionization in a constant temperature cloud with solar metallicity, subject to photoionization by an external radiation field¶

Tasks¶

1) Modify the CLOUDY input file used to solve the previous exercise to calculate the cooling and heating rates arising from the collisional ionization of a mixture at solar metallicity, assuming the solar abundance pattern. Also, assume that the cloud is subject to the external present-day Haardt \& Madau (2012)(*) ultraviolet background radiation.

2) What new features do you identify now in the cooling curve?

3) Plot the contribution of H, He, C, Ne, Fe, and Mg to the cooling function. What are the primary coolants of gas at solar metallicity, assuming the same solar abundance pattern?

Example solution¶

CLOUDY parameter file for this excersice

table hm12 redshift 0
hden -2 vary //Density of hydrogen in atoms per cubic centimetre
grid, range from -6 to 0 with 1.0 dex steps
constant temperature 4 vary // Temperature of the cloud
grid, range from 4 to 8 with 0.05 dex steps
abundances "solarC10.abn" 
set dr 0 //Thickness of the shell is 0
stop zone 1 // Stop after calculating one zone.
save cooling each "cooling.cool" no hash //Save the desired data
save grid "cooling.grd" no hash
save overview "cooling.ovr" no hash
In [5]:
import numpy as np
import matplotlib.pyplot as plt

nHmin = -6.0
nHmax = 0.0
nHsteps = 1.0

Tmin = 4.0
Tmax = 8.0
Tsteps = 0.05

NnHbins = 0
nH = nHmin
nHbins = []
while(nH < nHmax):
    nH = nHmin + nHsteps * NnHbins
    nHbins.append(nH)
    NnHbins += 1

NTbins = 0
T = Tmin
Tbins = []
while(T < Tmax):
    T = Tmin + Tsteps * NTbins
    Tbins.append(T)
    NTbins += 1

#Load the cooling rates for all elements
with open('tables/cooling_heating_Zsun/cooling.cool','r') as input_file:
    lines = input_file.readlines()
    header = lines[0].split()

    #dima is reported but not included explicitly in the table, so we do not need to read it? 
    cooling_rates = {}
    for i in header:
        cooling_rates[i] = np.zeros([NnHbins, NTbins])
        
    for i in range(NnHbins):
        for j in range(NTbins):
            for k in range(len(header)-1):
                cooling_rates[header[k]][i,j] = np.asarray(lines[i*NTbins+j+1].split()[k])
          
        
fig = plt.figure(1, figsize=(10,5), tight_layout=True)
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

for i in range(len(nHbins)):
    c = (nHbins[i]-nHmin)/(nHmax-nHmin)
    color = cm.jet(c)
    ax1.plot(Tbins, np.log10(cooling_rates['Ctot(erg/cm3/s)'][i,:])-2.0 * nHbins[i], label=r'${%.1f}$'% nHbins[i], color=color)

ax2.plot(Tbins, np.log10(cooling_rates['Ctot(erg/cm3/s)'][-1,:])-2.0 * nHbins[-1], color='k', lw=2)

#Now we plot the contribution to the cooling rate of H, He, C, Ne, Fe, and Mg
elements_list = ['H', 'He', 'O', 'N', 'Ne', 'Si', 'Fe', 'Mg']
for i in elements_list:
    ax2.plot(Tbins, np.log10(cooling_rates[i][-1,:])-2.0 * nHbins[-1],lw=1.0, label='%s'%i)
    
    
ax1.set_ylim(-24.0,-21)
ax2.set_ylim(ax1.get_ylim())

ax1.set_xlabel(r'$\rm log_{10} \left ( T / K \right ) $')
ax1.set_ylabel(r'$\rm log_{10} \Lambda $')

ax2.set_xlabel(r'$\rm log_{10} \left ( T / K \right ) $')
ax2.set_ylabel(r'$\rm log_{10} \Lambda $')

ax2.legend(ncol=2)

fig.suptitle('Radiative cooling and Heating of a cloud with solar matallicity + external HM12 UVB')
plt.show()

From now on, you can complicate these calculations further and construct cooling functions as a function of density, temperature, metallicity, and redshift. By doing so, you will build cooling function grids that can later be ported to simulation codes (e.g., Gadget, Arepo) to simulate the formation of galaxies and the IGM.

Congratulations, aided by CLOUDY, you have made an essential ingredient of simulations of galaxy formation and a tool to understand the interstellar medium.