""" Module for IGM calculations
"""
from __future__ import print_function, absolute_import, division, unicode_literals
import numpy as np
import os
from IPython import embed
import importlib_resources
from scipy.interpolate import interp1d
from scipy.interpolate import InterpolatedUnivariateSpline as IUS
from astropy import units
from astropy.table import Table
from astropy.utils import isiterable
from astropy import constants
from frb.halos import hmf as frb_hmf
from frb import mw
from frb import defs
[docs]
def fukugita04_dict():
"""
Data from Fukugita 2004, Table 1
Returns:
f04_dict (dict): data dict.
"""
f04_dict = {}
f04_dict['M_sphere'] = 0.0015
f04_dict['M_disk'] = 0.00055
f04_dict['M_HI'] = 0.00062
f04_dict['M_H2'] = 0.00016
f04_dict['M_WD'] = 0.00036
f04_dict['M_NS'] = 0.00005
f04_dict['M_BH'] = 0.00007
f04_dict['M_BD'] = 0.00014
# Return
return f04_dict
[docs]
def average_fHI(z, z_reion=7.):
"""
Average HI fraction
1 = neutral
0 = fully ionized
Args:
z (float or ndarray): redshift
z_reion (float, optional): redshift
of reionization.
Returns:
fHI (float or ndarray): float or ndarray
"""
z, flg_z = z_to_array(z)
fHI = np.zeros_like(z)
#
zion = z > z_reion
fHI[zion] = 1.
# Return
if flg_z:
return fHI
else:
return fHI[0]
[docs]
def average_He_nume(z, z_HIreion=7.):
"""
Average number of electrons contributed by He as a function of redshift
per He nucleus
Following Kulkarni, Worseck & Hennawi 2018
https://arxiv.org/abs/1807.09774
Args:
z (float or ndarray): Redshift
z_HIreion (float, optional): Helium reionization
redshift.
Returns:
neHe (float or ndarray): Number of free electrons
per Helium nucelus.
"""
z, flg_z = z_to_array(z)
# Load Kulkarni Table
He_file = importlib_resources.files('frb.data.IGM')/'qheIII.txt'
qHeIII = Table.read(He_file, format='ascii')
# Fully re-ionized
first_ionized = np.where(qHeIII['Q_HeIII_18'] >= 1.)[0][0]
z_HeIIreion = qHeIII['z'][first_ionized]
#
fHeI = np.zeros_like(z)
fHeII = np.zeros_like(z)
# HeI ionized at HI reionization
zion = z > z_HIreion
fHeI[zion] = 1.
# HeII ionized at HeII reionization
zion2 = (z > z_HeIIreion) & (z < z_HIreion)
fi_HeIII = interp1d(qHeIII['z'], qHeIII['Q_HeIII_18'])
fHeII[zion2] = 1. - fi_HeIII(z[zion2])
# Combine
neHe = (1.-fHeI) + (1.-fHeII) # No 2 on the second term as the first one gives you the first electron
# Return
if flg_z:
return neHe
else:
return neHe[0]
[docs]
def z_from_DM(DM, cosmo=defs.frb_cosmo, coord=None, corr_nuisance=True):
"""
Report back an estimated redshift from an input IGM DM
Any contributions from the Galaxy and/or host need to have been 'removed'
Args:
DM (Quantity): Dispersion measure.
cosmo (Cosmology, optional): Cosmology
of the universe. LambdaCDM with the Repo cosmology
used by default.
coord (SkyCoord, optional): If provided, use it to remove the ISM
corr_nuisance (bool, optional): If True, correct for the MW Halo
and the host with 100 DM units.
Returns:
z (float): Redshift
"""
if coord is not None:
DM_ISM = mw.ismDM(coord)
DM_use = DM - DM_ISM
else:
DM_use = DM
# Correct
if corr_nuisance:
DM_use -= 100 * units.pc/units.cm**3
# Calculate DMs
all_DM, zeval = average_DM(5., cosmo=cosmo, neval=20000, cumul=True)
# Interpolate
fint = interp1d(all_DM.value, zeval)
# Evaluate
z = fint(DM_use.to('pc/cm**3').value)
# Return
return z
[docs]
def f_diffuse(z, cosmo=defs.frb_cosmo,
return_rho:bool=False,
perturb_Mstar:float=None):
"""
Calculate the cosmic fraction of baryons
in diffuse gas phase based on our empirical
knowledge of baryon distributions and their
ionization state.
Note that the default values use the standard
values from Madau & Dickinson (2014) and Fukugita (2004).
The former use a Salpeter IMF for rho_* which is no longer
in fashion.
Args:
z (float or ndarray): Redshift
cosmo (Cosmology, optional): Cosmology of
the universe.
return_rho (bool, optional): If true,
the diffuse gas density
is returned too.
perturb_Mstar (float, optional):
If provided, scale rho_Mstar by this value.
Useful for exploring the uncertainty in f_diffuse
Returns:
f_diffuse (float, ndarray): Diffuse gas baryon fraction.
rho_diffuse (Quantity, optional): Physical diffuse gas density.
Returned if return_rho is set to true.
"""
# Get comoving baryon mass density
rho_b = cosmo.Ob0 * cosmo.critical_density0.to('Msun/Mpc**3')
# ##
# Dense components
# Stars
rho_Mstar = avg_rhoMstar(z, remnants=True)
if perturb_Mstar is not None:
rho_Mstar *= perturb_Mstar
# ISM
rho_ISM = avg_rhoISM(z, cosmo=cosmo, perturb_Mstar=perturb_Mstar)
# Diffuse gas fraction
f_diffuse = 1 - ((rho_Mstar+rho_ISM)/rho_b).value
if not return_rho:
return f_diffuse
else:
return f_diffuse, rho_b*f_diffuse*(1+z)**3
[docs]
def sigma_fd(z, rel_err_Mstar):
"""
Calculate the uncertainty in the diffuse fraction at a given redshift.
Args:
z (float or np.ndarray): The redshift value.
rel_err_Mstar (float): The relative error in the stellar mass.
Returns:
sigma_fd (float or np.ndarray): The uncertainty in the diffuse fraction.
"""
# Calculate the 3 values
f_d_low = f_diffuse(z, perturb_Mstar=1-rel_err_Mstar)
f_d_high = f_diffuse(z, perturb_Mstar=1+rel_err_Mstar)
# Calculate the sigma
s_fd = np.abs(f_d_high - f_d_low)/2.
# Return
return s_fd
[docs]
def ne_cosmic(z, cosmo = defs.frb_cosmo, mu = 4./3):
"""
Calculate the average cosmic electron
number density as a function of redshift.
Args:
z (float or ndarray): Redshift
cosmo (Cosmology, optional): Cosmology in
which the calculations are to be performed.
mu (float): Reduced mass
Returns:
ne_cosmic (Quantity): Average physical number
density of electrons in the unverse in cm^-3.
"""
# Get diffuse gas density
_, rho_diffuse = f_diffuse(z, cosmo=cosmo, return_rho=True)
# Number densities of H and He
n_H = (rho_diffuse/constants.m_p/mu).to('cm**-3')
n_He = n_H / 12. # 25% He mass fraction
# Compute electron number density
ne_cosmic = n_H * (1.-average_fHI(z)) + n_He*(average_He_nume(z))
return ne_cosmic
[docs]
def average_DM(z, cosmo = defs.frb_cosmo, cumul=False, neval=10000, mu=4/3):
"""
Macquart Relation
Calculate the average cosmic DM 'expected' based on our empirical
knowledge of baryon distributions and their ionization state.
This includes both the IGM and galactic halos, i.e. any and all diffuse gas
Args:
z (float): Redshift
mu (float): Reduced mass correction for He when calculating n_H
cumul (bool, optional): Return the DM as a function of z
Returns:
DM (Quantity or Quantity array): DM values evaluated at
the required redshifts. An array is returned only if
cumul is True.
zeval (ndarray): evaluation redshifts. Only returned if
cumul is True.
"""
# Init
zeval, dz = np.linspace(0., z, neval,retstep=True)
# Get n_e as a function of z
n_e = ne_cosmic(zeval, cosmo=cosmo)
# Cosmology -- 2nd term is the (1+z) factor for DM
denom = cosmo.H(zeval) * (1+zeval) * (1+zeval)
# Time to Sum
DM_cum = (constants.c * np.cumsum(n_e * dz / denom)).to('pc/cm**3')
# Return
if cumul:
return DM_cum, zeval
else:
return DM_cum[-1]
[docs]
def average_DMhalos(z, cosmo = defs.frb_cosmo, f_hot = 0.75, rmax=1.,
logMmin=10.3, logMmax=16., neval = 10000, cumul=False):
"""
Average DM_halos term from halos along the sightline to an FRB
Args:
z (float): Redshift of the FRB
cosmo (Cosmology): Cosmology in which the calculations
are to be performed.
f_hot (float, optional): Fraction of the halo baryons in diffuse phase.
rmax (float, optional): Size of a halo in units of r200
logMmin (float, optional): Lowest mass halos to consider
Cannot be much below 10.3 or the Halo code barfs
The code deals with h^-1 factors, i.e. do not impose it yourself
logMmax (float, optional): Highest halo mass. Default to 10^16 Msun
neval (int, optional): Number of redshift values between
0 and z the function is evaluated at.
cumul (bool, optional): Return a cumulative evaluation?
Returns:
DM_halos (Quantity or Quantity array): One value if cumul=False
else evaluated at a series of z
zeval (ndarray): Evaluation redshifts if cumul=True
"""
zeval, dz = np.linspace(0, z, neval, retstep = True)
# Electron number density in the universe
ne_tot = ne_cosmic(zeval, cosmo = cosmo)
# Diffuse gas mass fraction
f_diff = f_diffuse(zeval, cosmo = cosmo)
# Fraction of total mass in halos
zvals = np.linspace(0, z, 20)
fhalos = frb_hmf.frac_in_halos(zvals, Mlow = 10**logMmin, Mhigh = 10**logMmax, rmax = rmax)
fhalos_interp = IUS(zvals, fhalos)(zeval)
# Electron number density in halos only
ne_halos = ne_tot*fhalos_interp*f_hot/f_diff
# Cosmology -- 2nd term is the (1+z) factor for DM
denom = cosmo.H(zeval) * (1+zeval) * (1+zeval)
# DM halos
DM_halos = (constants.c * np.cumsum(ne_halos * dz / denom)).to('pc/cm**3')
# Return
if cumul:
return DM_halos, zeval
else:
return DM_halos[-1]
[docs]
def average_DMIGM(z, cosmo = defs.frb_cosmo,
f_hot = 0.75, rmax=1.,
logMmin=10.3, neval = 10000,
cumul=False, return_DMhalos=False):
"""
Estimate DM_IGM in a cumulative fashion
Args:
z (float): Redshift of the FRB
cosmo (Cosmology, optional): Cosmology in which
the calculations are to be performed. LambdaCDM
with the Repo cosmology assumed by default.
f_hot (float, optional): Fraction of the halo
baryons in diffuse phase.
rmax (float, optional):
Size of a halo in units of r200
logMmin (float, optional):
Lowest mass halos to consider. Cannot be much below
10.3 or the Halo code barfs. The code deals with
h^-1 factors, i.e. do not impose it yourself
neval (int, optional): Number of redshift values between
0 and z the function is evaluated at.
cumul (bool, optional):
Return a cumulative evaluation?
return_DMHalos (bool, optional): Also return avgDM_halos?
Returns:
float or list:
DM_IGM (Quantity or Quantity array): One value if cumul=False
else evaluated at a series of z
zeval (ndarray, optional): Evaluation redshifts if cumul=True
DM_halos (ndarray, optinal):
"""
# DM cosmic
DM_cosmic, zeval = average_DM(z, cosmo = cosmo, cumul=True, neval=neval)
# DM_halos
DM_halos, _ = average_DMhalos(z,cosmo = cosmo, logMmin = logMmin,
f_hot=f_hot, cumul = True, rmax = rmax, neval = neval)
# Subtract the two
DM_IGM = DM_cosmic - DM_halos
# Return
if cumul:
ret_val = [DM_IGM, zeval]
else:
ret_val = DM_IGM[-1]
if return_DMhalos:
if not isinstance(ret_val, list):
ret_val = [ret_val]
ret_val += [DM_halos]
# Finally
return ret_val
[docs]
def avg_rhoISM(z, cosmo=defs.frb_cosmo,
perturb_Mstar:float=None):
"""
Co-moving Mass density of the ISM
Interpolates from z=0 values to z=1 where
we assume M_ISM = M* and also for z>1
Args:
z (float or ndarray): Redshift
cosmo (Cosmology, optional): Cosmology in which
the calculations are to be performed. LambdaCDM
with defs.frb_cosmo parameters assumed by default.
perturb_Mstar (float, optional):
If provided, scale rho_Mstar by this value.
Useful for exploring the uncertainty in f_diffuse
Returns:
rhoISM (Quantity): Units of Msun/Mpc^3
"""
# Init
z, flg_z = z_to_array(z)
# Mstar
rhoMstar = avg_rhoMstar(z, remnants=False)
if perturb_Mstar is not None:
rho_Mstar *= perturb_Mstar
# z=0 (Fukugita+ 2004)
f04_dict = fukugita04_dict()
M_ISM = f04_dict['M_HI'] + f04_dict['M_H2']
f_ISM_0 = M_ISM/(f04_dict['M_sphere']+f04_dict['M_disk'])
# Assume M_ISM = M* at z=1
f_ISM_1 = 1.
# Ages
t0 = cosmo.age(0.).to('Gyr').value
t1 = cosmo.age(1.).to('Gyr').value
t1_2 = (t0+t1)/2.
tval = cosmo.age(z).to('Gyr').value
# Interpolate
f_ISM = interp1d([t0, t1_2, t1], [f_ISM_0, 0.58, f_ISM_1], kind='quadratic',
bounds_error=False, fill_value=1.)
# Calculate
rhoISM_unitless = f_ISM(tval) * rhoMstar.value
# Finish
rhoISM = rhoISM_unitless * units.Msun / units.Mpc**3
#
return rhoISM
[docs]
def avg_rhoMstar(z, remnants=True):
"""
Return mass density in stars as calculated by
Madau & Dickinson (2014).
Note: these authors assumed the Salpeter IMF
which is no longer in fashion.
Args:
z (float or ndarray): Redshift
remnants (bool, optional): Include remnants in the calculation?
Returns:
rho_Mstar (Quantity): Units of Msun/Mpc^3
"""
# Init
z, flg_z = z_to_array(z)
# Load
stellar_mass_file = importlib_resources.files('frb.data.IGM')/'stellarmass.dat'
rho_mstar_tbl = Table.read(stellar_mass_file, format='ascii')
# Output
rho_Mstar_unitless = np.zeros_like(z)
# Extrema
highz = z > rho_mstar_tbl['z'][-1]
rho_Mstar_unitless[highz] = rho_mstar_tbl['rho_Mstar'][-1]
# Interpolate
fint = interp1d(rho_mstar_tbl['z'], rho_mstar_tbl['rho_Mstar'], kind='cubic')
rho_Mstar_unitless[~highz] = fint(z[~highz])
# Finish
rho_Mstar = rho_Mstar_unitless * units.Msun / units.Mpc**3
# Remnants
if remnants:
# Fukugita 2004 (Table 1)
f04_dict = fukugita04_dict()
f_remnants = (f04_dict['M_WD'] + f04_dict['M_NS'] + f04_dict['M_BH'] + f04_dict['M_BD']) / (
f04_dict['M_sphere'] + f04_dict['M_disk'])
# Apply
rho_Mstar *= (1+f_remnants)
# Return
if flg_z:
return rho_Mstar
else:
return rho_Mstar[0]
[docs]
def avg_rhoSFR(z):
"""
Average SFR density
Based on Madau & Dickinson (2014)
Parameters
----------
z: float or ndarray
Redshift
Returns
-------
SFR: Quantity
Units of Msun/yr/Mpc^3
"""
rho_SFR_unitless = 0.015 * (1+z)**2.7 / (1 + ((1+z)/2.9)**5.6)
rho_SFR = rho_SFR_unitless * units.Msun / units.yr / units.Mpc**3
# Return
return rho_SFR
[docs]
def z_to_array(z):
"""
Convert input scalar or array to an array
Parameters
----------
z: float or ndarray
Redshift
Returns
-------
z: ndarray
flg_z: int
0 -- Input was a scalar
1 -- Input was an array
"""
# float or ndarray?
if not isiterable(z):
z = np.array([z])
flg_z = 0
else:
flg_z = 1
# Return
return z, flg_z