Zeus21 with Population II and III Stars
This quick tutorial covers all the new features in the updated public version of Zeus21 (see our paper Cruz et al. 2024) for more information.
In this new version, we have included the influence of Population III stars as well as the suppressive feedback effects from Lyman-Werner radiation and the relative velocity between CDM and baryons.
[1]:
import math
import numpy as np
from matplotlib import pyplot as plt
#import Zeus
import zeus21
#set up the CLASS cosmology
from classy import Class
#and the user parameters
UserParams = zeus21.User_Parameters(precisionboost=1.2)
Step 1: Set up Cosmology
We begin by running CLASS, where the associated parameters can be altered below. Then we save the cosmo parameters, the correlation functions, and the halo mass function at all desired z and M.
[2]:
# cosmo inputs from table 2 last column of 1807.06209, as 21cmFAST
ombh2 = 0.02242
omch2 = 0.11933
tau_re = 0.0544
hLittle = 0.6766
ns = 0.9665
As = np.exp(3.047)*10**(-10.)
[3]:
#set up user parameters
UserParams = zeus21.User_Parameters(FLAG_FORCE_LINEAR_CF=False,zmin_T21=10.)
#define all cosmology (including derived) parameters, and save them to the CosmoParams structure
CosmoParams = zeus21.Cosmo_Parameters(UserParams=UserParams, omegab = ombh2, omegac = omch2, h_fid = hLittle, As = As, ns = ns, tau_fid = tau_re, USE_RELATIVE_VELOCITIES = True, Flag_emulate_21cmfast=False)
print('CLASS has run, we store the cosmology.')
# Compute the HMF structure that stores associated quantities
HMFintclass = zeus21.HMF_interpolator(User_Parameters=UserParams,Cosmo_Parameters=CosmoParams)
print('HMF interpolator built. This ends the cosmology part -- moving to astrophysics.')
CLASS has run, we store the cosmology.
HMF interpolator built. This ends the cosmology part -- moving to astrophysics.
Step 2: Astrophysics
After that we set up the astro parameters, calculate the SFRD and with it all the global signal and related quantities. In principle one can re-run the astrophysics part only if you’re certain of the cosmology parameters.
With the current implementation, the computation of global quantities below takes less than a second on a laptop.
For more details on parameters, please refer to those used in Table 1 of our paper linked above.
NOTE: Pop III stars and LW feedback are turned on using the USE_POP_III and USE_LW_FEEDBACK flags below.
[ ]:
### POP II quantities first
################################
### Model Parameters
accretion_model = "exp" # Accretion model. "exp" for exponential, "EPS" for EPS. "RP16" for the dynamically averaged fitting function in Rodríguez-Puebla+16. Default is "exp"
################################
### SFR(Mh) Parameteres
alphastar = 0.69 # alphastar powerlaw index for low masses, default 0.5
betastar = -1.68 # betastar powerlaw index for high masses, default -0.5
epsstar = 10**-1.11 # epsilonstar = fstar at Mc
Mc = 10**11.93 # Pivot mass at which the power law cuts for model 0, default Mc = 3e11
dlog10epsstardz = -0.08 # dlog10epsilonstar/dz, default 0
################################
### Escape fraction parameters
fesc10 = 0.1 # fesc(M) parameter. Power law normalized (fesc10) at M=1e10 Msun with index alphaesc
alphaesc = 0.0
L40_xray = 10**0.5 # L40_xray: soft-band (E<2 keV) lum/SFR in Xrays in units of 10^40 erg/s/(Msun/yr)
E0_xray = 500. # E0_xray: minimum energy in eV
alpha_xray = -1.0 # Xray SED power-law index
Emax_xray_norm=2000 # max energy in eV to normalize SED. Keep at 2000 eV normally
################################
### LyA parameters
N_alpha_perbaryon_II = 9690 # number of Pop II photons between LyA and Ly Cont. per baryon (from BL05)
N_alpha_perbaryon_III = 17900 # number of Pop III photons between LyA and Ly Cont. per baryon value of 17900 is from Klessen & Glover 2023 (2303.12500), table A2
################################
### MTURN Parameters:
Mturn_fixed = None # Mturn_fixed: None if use Matom(z) at each z, Some value if fixed Mturn
FLAG_MTURN_SHARP= False # Mturn_sharp: False if regular exponential cutoff, True if sharp cutoff, active only if Mturn_fixed is on
################################
# Pop III Quantities
alphastar_III = 0
betastar_III = 0
epsstar_III = 10**-2.5709708032788794/3 #AV sugmaUV = 1.3
Mc_III = 1e7
dlog10epsstardz_III = 0.0
fesc7_III = 10**(-1.35)
alphaesc_III = -0.3
L40_xray_III = 10**0.5
alpha_xray_III = -1.0
USE_POPIII = True
USE_LW_FEEDBACK = True
A_LW = 2.0
beta_LW = 0.6
A_vcb = 1.0
beta_vcb = 1.8
AstroParams = zeus21.Astro_Parameters(CosmoParams=CosmoParams,
quadratic_SFRD_lognormal=False,
FLAG_USE_PSD=False,
accretion_model = accretion_model,
alphastar = alphastar,
betastar = betastar,
epsstar = epsstar,
Mc = Mc,
dlog10epsstardz = dlog10epsstardz,
fesc10 = fesc10,
alphaesc = alphaesc,
L40_xray = L40_xray,
E0_xray = E0_xray,
alpha_xray = alpha_xray,
Emax_xray_norm = Emax_xray_norm,
N_alpha_perbaryon_II = N_alpha_perbaryon_II,
N_alpha_perbaryon_III = N_alpha_perbaryon_III,
Mturn_fixed = Mturn_fixed,
FLAG_MTURN_SHARP = FLAG_MTURN_SHARP,
USE_POPIII = USE_POPIII,
USE_LW_FEEDBACK = USE_LW_FEEDBACK,
alphastar_III = alphastar_III,
betastar_III = betastar_III,
epsstar_III = epsstar_III,
Mc_III = Mc_III,
dlog10epsstardz_III = dlog10epsstardz_III,
fesc7_III = fesc7_III,
alphaesc_III = alphaesc_III,
L40_xray_III = L40_xray_III,
alpha_xray_III = alpha_xray_III,
A_LW = A_LW,
beta_LW = beta_LW,
A_vcb = A_vcb,
beta_vcb = beta_vcb,
)
CoeffStructure = zeus21.get_T21_coefficients(UserParams=UserParams, CosmoParams=CosmoParams ,AstroParams=AstroParams, HMFinterp=HMFintclass)
SFRD_class = zeus21.sfrd.SFRD_class(UserParams, CosmoParams, AstroParams, HMFintclass)
zlist= CoeffStructure.zintegral
print('SFRD and coefficients stored. Move ahead.')
[ ]:
### POP II quantities first
################################
### Model Parameters
accretion_model = "exp" # ACCRETION MODEL: 0 for exponential, 1 for EPS, default EXP
################################
### SFR(Mh) Parameteres
alphastar = 0.5 # alphastar powerlaw index for low masses, default 0.5
betastar = -0.5 # betastar powerlaw index for high masses, default -0.5
epsstar = 10**-1. # epsilonstar = fstar at Mc
Mc = 3e11 # Pivot mass at which the power law cuts for model 0, default Mc = 3e11
dlog10epsstardz = 0.0 # dlog10epsilonstar/dz, default 0
################################
### Escape fraction parameters
fesc10 = 0.1 # fesc(M) parameter. Power law normalized (fesc10) at M=1e10 Msun with index alphaesc
alphaesc = 0.0
L40_xray = 10**0.5 # L40_xray: soft-band (E<2 keV) lum/SFR in Xrays in units of 10^40 erg/s/(Msun/yr)
E0_xray = 500. # E0_xray: minimum energy in eV
alpha_xray = -1.0 # Xray SED power-law index
Emax_xray_norm=2000 # max energy in eV to normalize SED. Keep at 2000 eV normally
################################
### LyA parameters
N_alpha_perbaryon_II = 9690 # number of Pop II photons between LyA and Ly Cont. per baryon (from BL05)
N_alpha_perbaryon_III = 17900 # number of Pop III photons between LyA and Ly Cont. per baryon value of 17900 is from Klessen & Glover 2023 (2303.12500), table A2
################################
### MTURN Parameters:
Mturn_fixed = None # Mturn_fixed: None if use Matom(z) at each z, Some value if fixed Mturn
FLAG_MTURN_SHARP= False # Mturn_sharp: False if regular exponential cutoff, True if sharp cutoff, active only if Mturn_fixed is on
################################
ZMIN = 10.0 # down to which z we compute the evolution
################################
# Pop III Quantities
alphastar_III = 0
betastar_III = 0
epsstar_III = 10**(-3.0)
Mc_III = 1e7
dlog10epsstardz_III = 0.0
fesc7_III = 10**(-1.35)
alphaesc_III = -0.3
L40_xray_III = 10**0.5
alpha_xray_III = -1.0
USE_POPIII = True
USE_LW_FEEDBACK = True
A_LW = 2.0
beta_LW = 0.6
A_vcb = 1.0
beta_vcb = 1.8
#set up your astro parameters too, here the peak of f*(Mh) as an example
AstroParams = zeus21.Astro_Parameters(CosmoParams=CosmoParams,
quadratic_SFRD_lognormal=False,
FLAG_USE_PSD=False,
accretion_model = accretion_model,
alphastar = alphastar,
betastar = betastar,
epsstar = epsstar,
Mc = Mc,
dlog10epsstardz = dlog10epsstardz,
fesc10 = fesc10,
alphaesc = alphaesc,
L40_xray = L40_xray,
E0_xray = E0_xray,
alpha_xray = alpha_xray,
Emax_xray_norm = Emax_xray_norm,
N_alpha_perbaryon_II = N_alpha_perbaryon_II,
N_alpha_perbaryon_III = N_alpha_perbaryon_III,
Mturn_fixed = Mturn_fixed,
FLAG_MTURN_SHARP = FLAG_MTURN_SHARP,
USE_POPIII = USE_POPIII,
USE_LW_FEEDBACK = USE_LW_FEEDBACK,
alphastar_III = alphastar_III,
betastar_III = betastar_III,
epsstar_III = epsstar_III,
Mc_III = Mc_III,
dlog10epsstardz_III = dlog10epsstardz_III,
fesc7_III = fesc7_III,
alphaesc_III = alphaesc_III,
L40_xray_III = L40_xray_III,
alpha_xray_III = alpha_xray_III,
A_LW = A_LW,
beta_LW = beta_LW,
A_vcb = A_vcb,
beta_vcb = beta_vcb)
CoeffStructure = zeus21.get_T21_coefficients(UserParams=UserParams, CosmoParams=CosmoParams ,AstroParams=AstroParams, HMFinterp=HMFintclass)
SFRD_class = zeus21.sfrd.SFRD_class(UserParams, CosmoParams, AstroParams, HMFintclass)
zlist= CoeffStructure.zintegral
print('SFRD and coefficients stored. Move ahead.')
Step 3: Plotting globally-averaged quantities
[ ]:
plt.figure(figsize = (12, 6.75))
plt.plot(zlist, CoeffStructure.T21avg, color = "#E64D4E", linewidth = 3)
plt.xlim([10, 35])
plt.ylim([-120, 30])
plt.xlabel(r'$z$', fontsize = 30)
plt.ylabel(r'$T_{21}$ [mK]', fontsize = 30)
plt.xticks(fontsize=30)
plt.yticks(fontsize=30)
plt.tick_params(which='major', length=12, width=2, direction='in', top = True, bottom = True, left = True, right = True)
plt.tick_params(which='minor', length=5, width=2, direction='in', top = True, bottom = True, left = True, right = True)
plt.tick_params(axis="y", labelsize=30, pad = 10)
plt.tick_params(axis="x", labelsize=30, pad = 10)
plt.tight_layout()
plt.show()
[ ]:
plt.figure(figsize = (12, 6.75))
plt.plot(zlist,CoeffStructure.T_CMB , color = "#EE9063", linewidth = 3, linestyle = 'dotted', label = r'$T_\mathrm{CMB}$')
plt.plot(zlist,CoeffStructure.Tk_ad , color = "#0B92B1", linewidth = 3, label = r'$T_\mathrm{ad}$')
plt.plot(zlist,CoeffStructure.Tk_xray , color = "#E64D4E", linewidth = 3, label = r'$T_x$')
plt.plot(zlist,CoeffStructure.Tk_avg, color = "#665191", linewidth = 3, label = r'$T_k$')
plt.plot(zlist,CoeffStructure._invTs_avg**-1 , color = "#77AC54", linewidth = 3, label = r'$T_s$')
plt.xlim([10, 35])
plt.ylim([0, 100])
plt.xlabel(r'$z$', fontsize = 30)
plt.ylabel(r'$T$ [K]', fontsize = 30)
plt.xticks(fontsize=30)
plt.yticks(fontsize=30)
plt.tick_params(which='major', length=12, width=2, direction='in', top = True, bottom = True, left = True, right = True)
plt.tick_params(which='minor', length=5, width=2, direction='in', top = True, bottom = True, left = True, right = True)
plt.tick_params(axis="y", labelsize=30, pad = 10)
plt.tick_params(axis="x", labelsize=30, pad = 10)
plt.legend(fontsize=20, frameon = False)
plt.tight_layout()
plt.show()
[ ]:
plt.figure(figsize = (12, 6.75))
plt.semilogy(CoeffStructure.zintegral, CoeffStructure.SFRD_II_avg, color="#E64D4E", linewidth=3.0, label = 'Pop II')
plt.semilogy(CoeffStructure.zintegral, CoeffStructure.SFRD_III_avg, color="#0B92B1", linewidth=3.0, label = 'Pop III')
plt.semilogy(CoeffStructure.zintegral, CoeffStructure.SFRD_II_avg + CoeffStructure.SFRD_III_avg, color="#665191", linewidth=3.0, label = 'Pop II + III')
plt.xlim([10, 35])
plt.ylim(1e-7, 1e-1)
plt.xlabel(r'$z$', fontsize = 30)
plt.ylabel(r'$\dot{\overline{\rho}}_*(z)$', fontsize = 30)
plt.xticks(fontsize=30)
plt.yticks(fontsize=30)
plt.tick_params(which='major', length=12, width=2, direction='in', top = True, bottom = True, left = True, right = True)
plt.tick_params(which='minor', length=5, width=2, direction='in', top = True, bottom = True, left = True, right = True)
plt.tick_params(axis="y", labelsize=30, pad = 10)
plt.tick_params(axis="x", labelsize=30, pad = 10)
plt.legend(fontsize=20, frameon = False)
plt.tight_layout()
plt.show()
[ ]:
plt.figure(figsize = (12, 6.75))
plt.semilogy(CoeffStructure.zintegral, CoeffStructure.J_21_LW_II, color="#E64D4E", linewidth=3.0, label = 'Pop II')
plt.semilogy(CoeffStructure.zintegral, CoeffStructure.J_21_LW_III, color="#0B92B1", linewidth=3.0, label = 'Pop III')
plt.semilogy(CoeffStructure.zintegral, CoeffStructure.J_21_LW_II + CoeffStructure.J_21_LW_III, color="#665191", linewidth=3.0, label = 'Pop II + III')
plt.xlim([10, 35])
plt.ylim(1e-3, 1e2)
plt.xlabel(r'$z$', fontsize = 30)
plt.ylabel(r'$J_{21}$', fontsize = 30)
plt.xticks(fontsize=30)
plt.yticks(fontsize=30)
plt.tick_params(which='major', length=12, width=2, direction='in', top = True, bottom = True, left = True, right = True)
plt.tick_params(which='minor', length=5, width=2, direction='in', top = True, bottom = True, left = True, right = True)
plt.tick_params(axis="y", labelsize=30, pad = 10)
plt.tick_params(axis="x", labelsize=30, pad = 10)
plt.legend(fontsize=20, frameon = False)
plt.tight_layout()
plt.show()
[ ]:
plt.figure(figsize = (12, 6.75))
plt.semilogy(CoeffStructure.zintegral, CoeffStructure.xa_avg, color="#665191", linewidth=3.0, label = 'Pop II + III')
plt.xlim([10, 35])
plt.ylim(1e-3, 1e2)
plt.xlabel(r'$z$', fontsize = 30)
plt.ylabel(r'$x_\alpha$', fontsize = 30)
plt.xticks(fontsize=30)
plt.yticks(fontsize=30)
plt.tick_params(which='major', length=12, width=2, direction='in', top = True, bottom = True, left = True, right = True)
plt.tick_params(which='minor', length=5, width=2, direction='in', top = True, bottom = True, left = True, right = True)
plt.tick_params(axis="y", labelsize=30, pad = 10)
plt.tick_params(axis="x", labelsize=30, pad = 10)
plt.legend(fontsize=20, frameon = False)
plt.tight_layout()
plt.show()
[ ]:
plt.figure(figsize = (12, 6.75))
plt.semilogy(CoeffStructure.zintegral,CoeffStructure.xe_avg, color="#E64D4E", linewidth=3.0)
plt.xlim([10, 35])
plt.ylim([1e-4, 1e-3])
plt.xlabel(r'$z$', fontsize = 30)
plt.ylabel(r'$x_e$', fontsize = 30)
plt.xticks(fontsize=30)
plt.yticks(fontsize=30)
plt.tick_params(which='major', length=12, width=2, direction='in', top = True, bottom = True, left = True, right = True)
plt.tick_params(which='minor', length=5, width=2, direction='in', top = True, bottom = True, left = True, right = True)
plt.tick_params(axis="y", labelsize=30, pad = 10)
plt.tick_params(axis="x", labelsize=30, pad = 10)
# plt.legend(fontsize=20, frameon = False)
plt.tight_layout()
plt.show()
Step 4: Fluctuations
The last structure to run is Power_Spectra, which uses our lognormal and log-chisquare prescription to compute the statistics of IGM fluctuations in the form of two-point functions.
Note that the new version of Zeus includes a modified prescription for computing power spectra. The form for the 21-cm brightness temperature used in the old version of Zeus is
\begin{equation} T_{21}=T_0(z)\left(1+\delta-\delta_v\right) x_{\mathrm{HI}}\left(\frac{x_\alpha}{1+x_\alpha}\right)\left(1-\frac{T_{\mathrm{CMB}}}{T_c}\right) \end{equation}
where the leftmost \(\delta\) term is the total matter overdensity. Because \(T_{21}\) tracks neutral hydrogen fluctuations, we make the substitution \(\delta \rightarrow \delta_b\). Thus, the large-scale structure (LSS) terms that fold in to computing \(\Delta^2_{21}(k,z)\) now depend on the baryon power spectrum \(P_\mathrm{b}(k,z)\) and the baryon-dark matter cross power spectrum \(P_\mathrm{b,cdm}(k,z)\).
This new feature can be toggled with UserParams.USE_BARYON_FLAG = 1 for the new \(\delta_b\) prescription or UserParams.USE_BARYON_FLAG = 0 for the old \(\delta\) form.
[ ]:
RSDMODE = 1
UserParams.USE_BARYON_FLAG = 1
PS21 = zeus21.Power_Spectra(UserParams, CosmoParams, AstroParams, CoeffStructure, RSD_MODE = RSDMODE)
klist = PS21.klist_PS
Plotting Results
[ ]:
from scipy.interpolate import RegularGridInterpolator
interp = RegularGridInterpolator((zlist, klist), PS21.Deltasq_T21, method = 'cubic', bounds_error=False, fill_value=0)
interpLin = RegularGridInterpolator((zlist, klist), PS21.Deltasq_T21_lin, method = 'cubic', bounds_error=False, fill_value=0)
[ ]:
#choose a z to plot
kchoose=0.3
zlistHighRes = np.geomspace(zlist[0], zlist[-1], 1000)
powerSpectrum = interp((zlistHighRes, kchoose))
powerSpectrumLin = interpLin((zlistHighRes, kchoose))
plt.figure(figsize = (12, 6.75))
plt.semilogy(zlistHighRes, powerSpectrum, color="#E64D4E", linewidth=3.0, label = 'Nonlinear')
plt.semilogy(zlistHighRes, powerSpectrumLin, color="#0B92B1", linewidth=3.0, linestyle = 'dashed', label = 'Linear')
plt.xlim([10, 25])
plt.ylim([1,200])
plt.xlabel(r'$z$', fontsize = 30)
plt.ylabel(r'$\Delta^2_{21}\,\rm[mK^2]$', fontsize = 30)
plt.xticks(fontsize=30)
plt.yticks(fontsize=30)
plt.tick_params(which='major', length=12, width=2, direction='in', top = True, bottom = True, left = True, right = True)
plt.tick_params(which='minor', length=5, width=2, direction='in', top = True, bottom = True, left = True, right = True)
plt.tick_params(axis="y", labelsize=30, pad = 10)
plt.tick_params(axis="x", labelsize=30, pad = 10)
plt.legend(fontsize=20, frameon = False)
plt.tight_layout()
plt.show()
[ ]:
#choose a z to plot
zchoose=16
klistHighRes = np.geomspace(1e-4, 1e1, 1000)
powerSpectrum = interp((zchoose, klistHighRes))
powerSpectrumLin = interpLin((zchoose, klistHighRes))
plt.figure(figsize = (12, 6.75))
plt.loglog(klistHighRes, powerSpectrum, color="#E64D4E", linewidth=3.0, label = 'Nonlinear')
plt.loglog(klistHighRes, powerSpectrumLin, color="#0B92B1", linewidth=3.0, linestyle = 'dashed', label = 'Linear')
plt.xlim([1e-3,1e0])
plt.ylim([1e-1,1e3])
plt.xlabel(r'$k\,\rm [Mpc^{-1}]$', fontsize = 30)
plt.ylabel(r'$\Delta^2_{21}\,\rm[mK^2]$', fontsize = 30)
plt.xticks(fontsize=30)
plt.yticks(fontsize=30)
plt.tick_params(which='major', length=12, width=2, direction='in', top = True, bottom = True, left = True, right = True)
plt.tick_params(which='minor', length=5, width=2, direction='in', top = True, bottom = True, left = True, right = True)
plt.tick_params(axis="y", labelsize=30, pad = 10)
plt.tick_params(axis="x", labelsize=30, pad = 10)
plt.legend(fontsize=20, frameon = False)
plt.tight_layout()
plt.show()
For questions, don’t hesitate to reach out to hcruz2@jhu.edu!
[ ]: