notebook.community

Edit and run

Calculations related to the IGM [v1]

These will be presented in a forthcoming paper:  Simha & Prochaska (2019)




In [1]:



    


#%matplotlib notebook



    











In [2]:



    


# imports
from importlib import reload
import numpy as np
import os

from pkg_resources import resource_filename

from matplotlib import pyplot as plt

from astropy import units
from astropy.table import Table
from astropy.cosmology import Planck15

from frb import igm

Stellar mass (baryons locked up)





In [3]:



    


#stellar_mass_file = resource_filename('frb', 'data/IGM/stellarmass.dat')
stellar_mass_file = resource_filename('frb', os.path.join('data','IGM','stellarmass.dat'))



    











In [4]:



    


rho_mstar_tbl = Table.read(stellar_mass_file, format='ascii')



    











In [5]:



    


rho_mstar_tbl[0:5]



    


















    
Out[5]:






Table length=5



zt_Gyrrho_Mstar


float64float64float64


0.013.48576600000.0


0.112.18560400000.0


0.211.05542400000.0


0.310.06522900000.0


0.49.194501900000.0

Method





In [6]:



    


zval = np.linspace(0., 4., 100)
rho_Mstar = igm.avg_rhoMstar(zval, remnants=False)



    











In [7]:



    


plt.clf()
ax = plt.gca()
ax.plot(zval, np.log10(rho_Mstar.value))
# Label
ax.set_xlabel('z')
ax.set_ylabel(r'$\log \, \rho_* \; [M_\odot/ \rm Mpc^3]$ ')
plt.show()

Following Fukugita 2004 (Table 1)





In [8]:



    


M_sphere = 0.0015
M_disk = 0.00055
M_WD = 0.00036
M_NS = 0.00005
M_BH = 0.00007
M_BD = 0.00014



    











In [9]:



    


f_remnants = (M_WD+M_NS+M_BH+M_BD) / (M_sphere+M_disk)
f_remnants



    


















    
Out[9]:








0.30243902439024395

















In [10]:



    


rho_Mstar_full = igm.avg_rhoMstar(zval, remnants=True)



    











In [11]:



    


plt.clf()
ax = plt.gca()
ax.plot(zval, np.log10(rho_Mstar.value), label='No Remnants')
ax.plot(zval, np.log10(rho_Mstar_full.value), label='With Remnants')
# Label
ax.set_xlabel('z')
ax.set_ylabel(r'$\log \, \rho_* \; [M_\odot/ \rm Mpc^3]$ ')
# Legend
legend = plt.legend(loc='upper right', scatterpoints=1, borderpad=0.2,
                       handletextpad=0.1, fontsize='large')
plt.show()

ISM

$z=0$ -- Fukugita 2004





In [12]:



    


M_HI = 0.00062
M_H2 = 0.00016
M_ISM = M_HI + M_H2



    











In [13]:



    


M_ISM/(M_sphere+M_disk)



    


















    
Out[13]:








0.38048780487804873

$z>0$ -- Could use DLAs and [silly] K-S relation

For now, assume M_ISM = M* at z>1

Do it





In [14]:



    


rhoISM = igm.avg_rhoISM(zval)



    











In [15]:



    


plt.clf()
ax = plt.gca()
ax.plot(zval, np.log10(rhoISM.value), label='ISM')
ax.plot(zval, np.log10(rho_Mstar_full.value), label='M*')
# Label
ax.set_xlabel('z')
ax.set_ylabel(r'$\log \, \rho \; [M_\odot/ \rm Mpc^3]$ ')
# Legend
legend = plt.legend(loc='upper right', scatterpoints=1, borderpad=0.2,
                       handletextpad=0.1, fontsize='large')
plt.show()

$Y$ -- Helium

https://arxiv.org/abs/1807.09774




In [16]:



    


#He_file = resource_filename('frb', 'data/IGM/qheIII.txt')
He_file = resource_filename('frb', os.path.join('data','IGM','qheIII.txt'))



    











In [17]:



    


qHeIII = Table.read(He_file, format='ascii')



    











In [18]:



    


qHeIII



    


















    
Out[18]:






Table length=833



zQ_HeIII_18Q_HeIII_18_lQ_HeIII_18_uQ_HeIII_21Q_HeIII_21_lQ_HeIII_21_u


float64float64float64float64float64float64float64


12.01e-101e-101e-101e-101e-101e-10


11.9884.374927e-092.218632e-101.508848e-071.860788e-092.667322e-101.104288e-08


11.9768.743113e-093.519406e-103.034705e-073.67851e-094.431452e-102.215664e-08


11.9641.320629e-084.904414e-104.578895e-075.554472e-096.294819e-103.344437e-08


11.9521.776445e-086.373657e-106.141419e-077.488672e-098.257424e-104.490605e-08


11.942.24176e-087.927134e-107.722276e-079.48111e-091.031927e-095.654171e-08


11.9282.716575e-089.564846e-109.321466e-071.153179e-081.248034e-096.835133e-08


11.9163.200887e-081.128679e-091.093899e-061.36407e-081.474066e-098.033491e-08


11.9043.694699e-081.309297e-091.257485e-061.580786e-081.710021e-099.249246e-08


11.8924.19801e-081.498339e-091.422904e-061.803325e-081.955901e-091.04824e-07


.....................


2.1141.01.01.01.01.01.0


2.1021.01.01.01.01.01.0


2.091.01.01.01.01.01.0


2.0781.01.01.01.01.01.0


2.0661.01.01.01.01.01.0


2.0541.01.01.01.01.01.0


2.0421.01.01.01.01.01.0


2.031.01.01.01.01.01.0


2.0181.01.01.01.01.01.0


2.0061.01.01.01.01.01.0

Plot





In [19]:



    


plt.clf()
ax=plt.gca()
ax.plot(qHeIII['z'], qHeIII['Q_HeIII_18'])
#
ax.set_xlabel('z')
ax.set_ylabel('f(HeIII)')
#
plt.show()

DM -- Piece by piece (as coded)

$\rho_b = \Omega_b \rho_c (1+z)^3$

$\rho_{\rm diffuse} = \rho_b - (\rho_* + \rho_{\rm ISM})$

$\rho_*$ is the mass density in stars

$\rho_{\rm ISM}$ is the mass density in the neutral ISM

Number densities

$n_{\rm H} = \rho_{\rm diffuse}/(m_p \, \mu)$

$\mu \approx 1.3$ accounts for Helium

$n_{\rm He} = n_{\rm H}/12$

$n_e = n_{\rm H} [1-f_{\rm HI}] + n_{\rm He} Y$

$f_{\rm HI}$ is the fraction atomic Hydrogen [value betwee 0-1]

$Y$ gives the number of free electrons per He nucleus [value between 0-2]

Integrating

$DM = \int \frac{n_e \, dr}{1+z} = \frac{c}{H_0} \int \frac{n_e \, dz}{(1+z)^2 \sqrt{(1+z)^3 \Omega_m +

                                                                             \Omega_\Lambda}}$

DM -- Altogether (using the code)





In [20]:



    


reload(igm)
DM = igm.average_DM(1.)
DM



    


















    
Out[20]:







$1053.0827 \; \mathrm{\frac{pc}{cm^{3}}}$

Cumulative plot





In [21]:



    


DM_cumul, zeval = igm.average_DM(1., cumul=True)



    











In [22]:



    


# Inoue approximation
DM_approx = 1000. * zeval * units.pc / units.cm**3



    











In [23]:



    


plt.clf()
ax = plt.gca()
ax.plot(zeval, DM_cumul, label='JXP')
ax.plot(zeval, DM_approx, label='Approx')
# Label
ax.set_xlabel('z')
ax.set_ylabel(r'${\rm DM}_{\rm IGM} [\rm pc / cm^3]$ ')
# Legend
legend = plt.legend(loc='lower right', scatterpoints=1, borderpad=0.2,
                       handletextpad=0.1, fontsize='large')
plt.show()



    


















    
























In [24]:



    


plt.clf()
ax = plt.gca()
ax.plot(zeval, DM_approx/DM_cumul, label='Approx/JXP')
#ax.plot(zeval, DM_approx, label='Approx')
# Label
ax.set_xlabel('z')
ax.set_ylabel(r'Ratio of ${\rm DM}_{\rm IGM} [\rm pc / cm^3]$ ')
# Legend
legend = plt.legend(loc='upper right', scatterpoints=1, borderpad=0.2,
                       handletextpad=0.1, fontsize='large')
plt.show()



    


















    
























In [25]:



    


DM_cumul[0:10]



    


















    
Out[25]:







$[0.083956468,~0.16792254,~0.25189822,~0.3358835,~0.41987839,~0.50388286,~0.58789694,~0.6719206,~0.75595385,~0.83999669] \; \mathrm{\frac{pc}{cm^{3}}}$
















In [26]:



    


DM_approx[0:10]



    


















    
Out[26]:







$[0.10001,~0.20002,~0.30003,~0.40004,~0.50005001,~0.60006001,~0.70007001,~0.80008001,~0.90009001,~1.0001] \; \mathrm{\frac{pc}{cm^{3}}}$
















In [27]:



    


zeval[0]



    


















    
Out[27]:








0.00010001000100010001

Development

Content source: FRBs/FRB

Similar notebooks:

notebook.community | gallery | about