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 shape of the UVB

This exercise aims to study how the external radiation field's shape impacts a gas cloud's ionization state.

Tasks:

1) Modify some of the CLOUDY input parameters used previously for CLOUDY to consider a gas at a constant temperature, $T = 10^3 \rm K$, and density, $\rm n_{\rm H} = 10^{-2} \ cm^{-3}$, and solar metallicity.

Consider gas illuminated by an external radiation field. Instead of using the given HM12 UVB spectrum, consider a power law continuum and vary the power law index in the range $-3 \le n \le 0$. Also, vary the continuum "intensity" in the range $-6 \le \rm log_{10} \left ( 4 \pi J / cm^{-2} s^{-1} \right ) \le -1$ (see Sec. 2.4 in Hazy1.pdf).

Use this setup to study (plot) how the ionization state of metals (e.g., CII, CIII, CIV, and CV) change as a function of the intensity of the incoming radiation field and the shape of its spectrum.

• Discuss why the ionization state of the metals depends on the intensity of the external radiation field.
• Discuss why the metal's ionisation state depends on the spectrum's shape at fixed intensity.

Repeat this exercise for other metals and ionic species.

Example solution

CLOUDY parameter file for this excersice

table power law 0
intensity -3 vary 
grid -6 -1 1
hden -2 //Density of hydrogen in atoms per cubic centimetre
constant temperature 4 
abundances "solarC10.abn" 
set dr 0 //Thickness of the shell is 0
stop zone 1 // Stop after calculating one zone.
save overview "ubv_shape_0.ovr" no hash
save incident continuum "ubv_continnum_0.cont" no hash

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

Imin = -6
Imax = -1
Isteps = 1.0

NIbins = 0
I = Imin
Ibins = []
while(I < Imax):
    I = Imin + Isteps * NIbins
    Ibins.append(I)
    NIbins += 1


#Load the tables
slopes = [0,1,2,3]
input_files = ['tables/ubv_shape/ubv_shape_%d.ovr'%i for i in slopes]

fig = plt.figure(1, figsize=(7,7), tight_layout=True)
ax1 = fig.add_subplot(221)
ax2 = fig.add_subplot(222)
ax3 = fig.add_subplot(223)
ax4 = fig.add_subplot(224)

for input_file, slope in zip(input_files, slopes):
    
    with open(input_file,'r') as fo:
        lines = fo.readlines()
        header = lines[0].split()

        ovr_dic = {}
        for i in header:
            ovr_dic[i] = np.zeros([NIbins])

        for i in range(NIbins):
            for k in range(len(header)-1):
                ovr_dic[header[k]][i] = np.asarray(lines[i+1].split()[k])

    total_carbon = (ovr_dic['C1']+ovr_dic['C2']+ovr_dic['C3']+ovr_dic['C4'])
    ax1.plot(Ibins, np.log10(ovr_dic['C1'] / total_carbon), 'o-', label='-%d'%slope)
    ax2.plot(Ibins, np.log10(ovr_dic['C2'] / total_carbon), 'o-', label='-%d'%slope)
    ax3.plot(Ibins, np.log10(ovr_dic['C3'] / total_carbon), 'o-', label='-%d'%slope)
    ax4.plot(Ibins, np.log10(ovr_dic['C4'] / total_carbon), 'o-', label='-%d'%slope)
    
    
ax1.set_ylim(-6,0.5)
ax2.set_ylim(ax1.get_ylim())
ax3.set_ylim(ax1.get_ylim())
ax4.set_ylim(ax1.get_ylim())

axes = [ax1, ax2, ax3, ax4]
[ax.set_xlabel(r'$\rm log_{10} \left ( Intensity \right )$', size=15) for ax in axes]
[ax.set_ylabel(r'$\rm log_{10} \left ( f_{C%d} \right )$'%i, size=15) for ax,i in zip(axes, [1,2,3,4])]

ax2.legend(title='continuum index')

<matplotlib.legend.Legend at 0x7fb38aa48760>

This is easier to understand if we plot the continuum

#Read Sec. 6.12.9 in Hazy1.pdf
input_files = ['tables/ubv_shape/ubv_continnum_%d.cont'%i for i in range(4)]

N = 6
data = np.loadtxt(input_files[0])
Ndata = len(data)

fig = plt.figure(1, figsize=(10,10), tight_layout=True)
axes = [fig.add_subplot(3,2,i+1) for i in range(6)]

slopes = [0,1,2,3]

for i,ax in zip(range(6), axes):
    for input_file, slope in zip(input_files, slopes):
        n1 = int(Ndata/6) * i
        n2 = int(Ndata/6) * (i+1)
        data = np.loadtxt(input_file)
        
        ax.plot(np.log10(data[n1:n2,0]*13.6), np.log10(data[n1:n2,1]/data[n1:n2, 0]), label="%-d"%slope)
    ax.axvline(np.log10(11.2), ls='--')
    ax.set_ylim(-10,6)
    ax.legend()
    ax.grid()
    ax.set_xlabel(r'$\rm log_{10} \left ( h \nu / eV \right )$')
    ax.set_ylabel(r'$\rm log_{10} \left ( 4 \pi J_{\nu} / h \right ) [cm^{-2} s^{-1}$')
    
plt.show()