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

Excercise 1: The onset of self-shielding

This exercise aims to study the density above which a parcel of gas with primordial abundance becomes self-shielded from an external photoionizing radiation field.

Tasks:

1) Modify the CLOUDY input parameter used previously to calculate the cooling rate of a Hydrogen + Helium gas to run CLOUDY in the density range $\rm -6 \le log_{10} \left ( n_{H} / cm^{-3} \right ) \le 1$, and the temperature range, $\rm 3 \le log_{10} \left ( T / K \right ) \le 6$. Include the presence of the external HM12 UVB.

2) Plot the fraction of neutral hydrogen and helium as a function of density and temperature. Discuss the results. Above which density does the warm ISM become self-shielded from the external UVB?

3) Assuming ionization equilibrium, i.e., that the number of ionization balances the number of recombination events, it is possible to find an approximation to the neutral hydrogen fraction, $f_{\rm HI}$, present in the gas. This approximation is expressed as a function of the hydrogen recombination rate, $\alpha_{\rm HII}(T)$, the hydrogen collisional ionization rate, $\Gamma_{\rm H, coll}$, and the hydrogen photoionization rate, $\Gamma_{\rm H, UVB}$.

The photoionization rate is:

$\Gamma_{\rm H, UVB}$ = $\displaystyle\int_{\nu_{\rm th}}^{\infty} \displaystyle\frac{J_{\nu}}{h \nu} \sigma_{\rm H}(\nu) d\nu $,

where $J(\nu)$ is the shape of the HM12 spectrum, $\nu_{\rm th}$ is the threshold frequency for ionization of hydrogen, and $\sigma_{H}(\nu)$ is the photoionization cross-section for hydrogen (which depends on density).

$\alpha_{\rm HII}(T)$ and $\Gamma_{\rm H, coll}$ are functions of the temperature, and can be taken from, e.g., Cen (1992) (https://ui.adsabs.harvard.edu/abs/1992ApJS...78..341C/abstract)

Your task is to find out the fraction of neutral hydrogen, $f_{\rm HI}$, as a function of temperature, in the optically-thin limit, assuming the contribution of the ionized helium to the electron number density is negligible. In this case, $f_{\rm HI}$ depends on density and temperature.

4) Show that in the optically-thick limit ($\Gamma_{\rm H, UVB} \to 0$), the approximation becomes:

$\rm f_{\rm HI} = \displaystyle\frac{\alpha_{HII}}{\alpha_{HII}+\Gamma_{H, coll}}$,

which is a function of temperature only. Plot this solution on top of the neutral hydrogen fraction vs Temperature that results from CLOUDY. Does this solution agree with the CLOUDY results? Where do the discrepancies originate from?

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 1 with 1.0 dex steps
constant temperature 4 vary // Temperature of the cloud
grid, range from 3 to 6 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

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

nHmin = -6.0
nHmax = 2.0
nHsteps = 1.0

Tmin = 3.0
Tmax = 6.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 tables
with open('tables/self-shielding/cooling.ovr','r') as input_file:
    lines = input_file.readlines()
    header = lines[0].split()

    ovr_dic = {}
    for i in header:
        ovr_dic[i] = np.zeros([NnHbins, NTbins])
        
    for i in range(NnHbins):
        for j in range(NTbins):
            for k in range(len(header)-1):
                ovr_dic[header[k]][i,j] = np.asarray(lines[i*NTbins+j+1].split()[k])

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(NnHbins):
    c = (nHbins[i]-nHmin)/(nHmax-nHmin)
    color = cm.jet(c)                
    
    n_HII = ovr_dic['HII'][i,:]
    n_HI  = ovr_dic['HI'][i,:]
    
    n_HeIII  = ovr_dic['HeIII'][i,:]
    n_HeII = ovr_dic['HeII'][i,:]
    n_HeI  = ovr_dic['HeI'][i,:]

    ax1.plot(Tbins, np.log10(n_HI/(n_HI+n_HII)), color=color, label=r'${%.1f}$'% nHbins[i])
    ax2.plot(Tbins, np.log10(n_HeI/(n_HeI+n_HeII+n_HeIII)), color=color, label=r'${%.1f}$'% nHbins[i])
    
n_HII = ovr_dic['HII'][:,20]
n_HI  = ovr_dic['HI'][:,20]
n_HeIII  = ovr_dic['HeIII'][:,20]
n_HeII = ovr_dic['HeII'][:,20]
n_HeI  = ovr_dic['HeI'][:,20]
    
ax3.plot(nHbins, np.log10(n_HI/(n_HI+n_HII)), 'o-', label='HI fraction')    
ax3.plot(nHbins, np.log10(n_HeI/(n_HeI+n_HeII+n_HeIII)), 'o-', label='HeI fraction')


ax1.set_ylim(-7,0.1)
ax1.set_xlim(3.0,6.0)

ax2.set_ylim(ax1.get_ylim())
ax2.set_xlim(ax1.get_xlim())

ax1.set_title('Hydrogen')
ax2.set_title('Helium')

ax1.set_xlabel(r'$\rm log_{10} \left ( T / K \right )$', size=15)
ax1.set_ylabel(r'$\rm log_{10} \left ( n_{HI} / n_{H} \right )$', size=15)

ax2.set_xlabel(r'$\rm log_{10} \left ( T / K \right )$', size=15)
ax2.set_ylabel(r'$\rm log_{10} \left ( n_{HeI} / n_{He} \right )$', size=15)

ax3.set_xlabel(r'$\rm log_{10} \left ( n_{H} / cm^{-3} \right )$', size=15)
ax3.set_ylabel(r'$\rm log_{10} \left ( n_{0} / n_{total} \right )$', size=15)

ax2.legend(title=r'$\rm log_{10} \left ( n_{H} / cm^{-3} \right )$')
ax3.legend()                
                

<matplotlib.legend.Legend at 0x7fe826db3700>

To answer question 4), we consider the following code:

def alpha_HII(T):
    T3 = T*1e-3
    T6 = T*1e-6
    return 8.4e-11 * T**(-0.5) * T3**(-1./5.) / (1.0+T6**(0.7))

def Gamma_H_coll(T):
    T5 = T*1e-5
    return 5.85e-11 * T**(0.5) / (1.0+T5**0.5) * np.exp(-157809.1/T)

def fHI_self_shielding(T):
    return alpha_HII(T) / (alpha_HII(T)+Gamma_H_coll(T))

fig = plt.figure(1, figsize=(5,5), tight_layout=True)
ax1 = fig.add_subplot(111)

for i in range(NnHbins):
    c = (nHbins[i]-nHmin)/(nHmax-nHmin)
    color = cm.jet(c)                
    
    n_HII = ovr_dic['HII'][i,:]
    n_HI  = ovr_dic['HI'][i,:]
    
    n_HeIII  = ovr_dic['HeIII'][i,:]
    n_HeII = ovr_dic['HeII'][i,:]
    n_HeI  = ovr_dic['HeI'][i,:]

    ax1.plot(Tbins, np.log10(n_HI/(n_HI+n_HII)), color=color)
    
T = 10**np.array(Tbins)
ax1.plot(Tbins, np.log10(fHI_self_shielding(T)), ls='--', color='k', label='approximation')

ax1.set_xlabel(r'$\rm log_{10} \left ( T / K \right )$', size=15)
ax1.set_ylabel(r'$\rm log_{10} \left ( f_{HI} \right )$', size=15)

ax1.legend();