#Obtain Matter Power Spectrum by solving Boltzmann code for the simplest massless neutrinos + LCDM model #by Zhiqi Huang, huangzhq25@mail.sysu.edu.cn import numpy as np from scipy.integrate import quad, odeint from scipy.interpolate import interp1d import matplotlib.pyplot as plt #******************************* cosmological parameters ********************************** #*******************************Modify this if you want **************************************** Omega_c = 0.26 Omega_b = 0.048 Hubble_reduced = 0.6748 #h T_CMB = 2.726 #CMB temperature Yp = 0.25 #Helium mass abundance Neff = 3.046 #effective Neutrino species opt_re = 0.058 #optical depth; set to zero if do not want reionization A_s = 2.1e-9 #primordial power amplitude n_s = 0.965 #primordial power index khMpc_pivot = 0.05/Hubble_reduced #pivot k in h Mpc^{-1} #============================================================================================ #*********************************RECFAST**************************************************************************** #original RECFAST Fortran version: see Seager, S., Sasselov, D. & Scott, D., 2000, ApJS, 128, 407 (astro-ph/9912182). #*********************************use SI units**************************************************************************** #===================== constants for unit conversion etc. =========================================================== Hunit =1.e5/3.0856775807e22 #100km/s/Mpc in SI unit lightspeed = 2.99792458e8 #speed of light k_B = 1.380658e-23 #Boltzmann constant h_P = 6.6260755e-34 #Planck constant m_e = 9.1093897e-31 #electron mass m_H = 1.673575e-27 #H mass He_H_massratio = 3.9715 #He/H mass ratio sigma_T = 6.6524616e-29 #Thomson cross-section blackbody_a = 7.565914e-16 #radiation constant u = aT^4 NewtonG = 6.6742e-11 #Newton gravitational constant A2s1s_H = 8.2245809 #H 2s-> 1s two photon rate A2s1s_He = 51.3 #He 2s-> 1s two photon rate L_H_ion = 1.096787737e7 #level for H ionization in m^-1 L_H_alpha = 8.225916453e6 #level for H Ly alpha in m^-1, averaged over 2 levels L_He1_ion = 1.98310772e7 #level for HeI ionization, from Drake (1993) L_He2_ion = 4.389088863e7 #level for HeII ionization, from JPhysChemRefData (1987) L_He_2s = 1.66277434e7 #level for HeI 2s, from Drake (1993) L_He_2p = 1.71134891e7 #level for He 2p (21P1-11S0) in m^-1, from Drake (1993) A2P_s = 1.798287e9 #Einstein A coefficient for He 21P1-11S0 A2P_t = 177.58 #Einstein A coefficient for He 23P1-11S0 L_He_2Pt = 1.690871466e7 #level for 23P012-11S0 in m^-1 L_He_2st = 1.5985597526e7 #level for 23S1-11S0 in m^-1 L_He2St_ion = 3.8454693845e6 #level for 23S1-continuum in m^-1 sigma_He_2Ps = 1.436289e-22 #H ionization x-section at HeI 21P1-11S0 freq. in m^2 sigma_He_2Pt = 1.484872e-22 #H ionization x-section at HeI 23P1-11S0 freq. in m^2 #H fudging fitting AGauss1 = -0.14 #Amplitude of 1st Gaussian AGauss2 = 0.079 #Amplitude of 2nd Gaussian zGauss1 = 7.28 #ln(1+z) of 1st Gaussian zGauss2 = 6.73 #ln(1+z) of 2nd Gaussian wGauss1 = 0.18 #Width of 1st Gaussian wGauss2 = 0.33 #Width of 2nd Gaussian #========================= derived parameters ============================================= #derived cosmological parameters Omega_m = Omega_c+Omega_b HubbleConstant = Hubble_reduced * Hunit mu_H = 1.e0/(1.e0-Yp) #Mass per H atom mu_T = He_H_massratio/(He_H_massratio-(He_H_massratio-1.e0)*Yp) #Mass per atom fHe = Yp/(He_H_massratio*(1.e0-Yp)) #n_He_tot / n_H_tot n_H_now = 3.e0*HubbleConstant**2*Omega_b/(8.e0*np.pi*NewtonG*mu_H*m_H) #H number density fnu = (7.e0/8.e0*Neff)*(4.e0/11.e0)**(4.e0/3.e0) Omega_g = (blackbody_a*T_CMB**4)/(3.e0*(HubbleConstant*lightspeed)**2/(8.e0*np.pi*NewtonG)) Omega_nu = fnu*Omega_g Omega_r = Omega_g + Omega_nu Omega_L = 1.-Omega_c-Omega_b-Omega_r ##assuming flat universe z_eq = Omega_m/Omega_r - 1.e0 a_eq = Omega_r/Omega_m #derived atomic parameters Lalpha = 1.e0/L_H_alpha Lalpha_He = 1.e0/L_He_2p DeltaB = h_P*lightspeed*(L_H_ion-L_H_alpha) CDB = DeltaB/k_B DeltaB_He = h_P*lightspeed*(L_He1_ion-L_He_2s) #2s, not 2p CDB_He = DeltaB_He/k_B CB1 = h_P*lightspeed*L_H_ion/k_B CB1_He1 = h_P*lightspeed*L_He1_ion/k_B #ionization for HeI CB1_He2 = h_P*lightspeed*L_He2_ion/k_B #ionization for HeII CR = 2.e0*np.pi*(m_e/h_P)*(k_B/h_P) CK = Lalpha**3/(8.e0*np.pi) CK_He = Lalpha_He**3/(8.e0*np.pi) CL = lightspeed*h_P/(k_B*Lalpha) CL_He = lightspeed*h_P/(k_B/L_He_2s) #comes from det.bal. of 2s-1s CT = (8.e0/3.e0)*(sigma_T/(m_e*lightspeed))*blackbody_a Bfact = h_P*lightspeed*(L_He_2p-L_He_2s)/k_B H_frac = 1.e-3 fu=1.125 #H fudge b_He = 0.86 #Helium fudge #=========================================================================================================================== #Hubble parameter H = dot a / a as a function of z, SI unit def recfast_Hofz(z): return HubbleConstant*np.sqrt(Omega_L + (Omega_m + Omega_r*(1+z))*(1.+z)**3) #dHdzofz, SI unit def recfast_dHdzofz(z): return (recfast_Hofz(z+0.05)-recfast_Hofz(z-0.05))/0.1 #I use a numeric one so you don't need to rewrite #main RECFAST ODE def ion_eq(y, z): f = np.empty(3) a_PPB = 4.309 b_PPB = -0.6166 c_PPB = 0.6703 d_PPB = 0.5300 a_VF = 10 **(-16.744) b_VF = 0.711 T_0 = 10 **(0.477121) #3K T_1 = 10 **(5.114) a_trip = 10 **(-16.306) b_trip = 0.761 x_H = y[0] x_He = y[1] x = x_H + fHe * x_He Tmat = y[2] n = n_H_now * (1 +z)**3 n_He = fHe * n Trad = T_CMB * (1 +z) Hz = recfast_Hofz(z) dHdz = recfast_dHdzofz(z) Rdown=1.e-19*a_PPB*(Tmat/1.e4)**b_PPB /(1 +c_PPB*(Tmat/1.e4)**d_PPB) Rup = Rdown * (CR*Tmat)**(1.5)*np.exp(-CDB/Tmat) sq_0 = np.sqrt(Tmat/T_0) sq_1 = np.sqrt(Tmat/T_1) Rdown_He = a_VF/(sq_0*(1 +sq_0)**(1 -b_VF)) Rdown_He = Rdown_He/(1 +sq_1)**(1 +b_VF) Rup_He = Rdown_He*(CR*Tmat)**(1.5)*np.exp(-CDB_He/Tmat) Rup_He = 4 *Rup_He #statistical weights factor for HeI if((Bfact/Tmat) > 680.): He_Boltz = np.exp(680.) else: He_Boltz = np.exp(Bfact/Tmat) K = CK/Hz*(1.0+AGauss1*np.exp(-((np.log(1.0+z)-zGauss1)/wGauss1)**2 )+AGauss2*np.exp(-((np.log(1.0+z)-zGauss2)/wGauss2)**2 )) Rdown_trip = a_trip/(sq_0*(1 +sq_0)**(1.0-b_trip)) Rdown_trip = Rdown_trip/((1 +sq_1)**(1 +b_trip)) Rup_trip = Rdown_trip*np.exp(-h_P*lightspeed*L_He2St_ion/(k_B*Tmat)) Rup_trip = Rup_trip*((CR*Tmat)**(1.5))*(4 /3 ) if ((x_He < 5.e-9) or (x_He > 0.980)): Heflag = 0 else: Heflag = 6 if (Heflag==0): K_He = CK_He/Hz else: optHe_s = A2P_s*CK_He*3 *n_He*(1 -x_He)/Hz pHe_s = (1 - np.exp(-optHe_s))/optHe_s K_He = 1 /(A2P_s*pHe_s*3 *n_He*(1 -x_He)) if ( Heflag >= 5 and (x_H < 0.9999999)): Doppler = 2 *k_B*Tmat/(m_H*He_H_massratio*lightspeed**2) Doppler = lightspeed*L_He_2p*np.sqrt(Doppler) gamma_2Ps = 3 *A2P_s*fHe*(1 -x_He)*lightspeed**2/(np.sqrt(np.pi)*sigma_He_2Ps*8 *np.pi*Doppler*(1 -x_H))/((lightspeed*L_He_2p)**2 ) pb = 0.36 #value from KIV (2007) qb = b_He AHcon = A2P_s/(1 +pb*(gamma_2Ps**qb)) K_He=1 /((A2P_s*pHe_s+AHcon)*3 *n_He*(1 -x_He)) if (Heflag >= 3): optHe_t = A2P_t*n_He*(1 -x_He)*3 optHe_t = optHe_t /(8 *np.pi*Hz*L_He_2Pt**(3 )) pHe_t = (1 - np.exp(-optHe_t))/optHe_t CL_PSt = h_P*lightspeed*(L_He_2Pt - L_He_2st)/k_B if (x_H > 0.99999): CfHe_t = A2P_t*pHe_t*np.exp(-CL_PSt/Tmat) CfHe_t = CfHe_t/(Rup_trip+CfHe_t) else: Doppler = 2 *k_B*Tmat/(m_H*He_H_massratio*lightspeed**2) Doppler = lightspeed*L_He_2Pt*np.sqrt(Doppler) gamma_2Pt = 3 *A2P_t*fHe*(1 -x_He)*lightspeed**2/(np.sqrt(np.pi)*sigma_He_2Pt*8 *np.pi*Doppler*(1 -x_H))/((lightspeed*L_He_2Pt)**2 ) pb = 0.66 qb = 0.9 AHcon = A2P_t/(1 +pb*gamma_2Pt**qb)/3 CfHe_t = (A2P_t*pHe_t+AHcon)*np.exp(-CL_PSt/Tmat) CfHe_t = CfHe_t/(Rup_trip+CfHe_t) #"C" factor for triplets timeTh=(1. /(CT*Trad**4))*(1 +x+fHe)/x #Thomson time timeH=2. /(3. *HubbleConstant*(1 +z)**1.5) #Hubble time if (x_H > 0.99): f[0] = 0. elif (x_H > 0.985): f[0] = (x*x_H*n*Rdown - Rup*(1 -x_H)*np.exp(-CL/Tmat))/(Hz*(1 +z)) factor=(1 + K*A2s1s_H*n*(1 -x_H))/(Hz*(1 +z)*(1 +K*A2s1s_H*n*(1 -x)+K*Rup*n*(1 -x))) else: #use full rate for H f[0] = ((x*x_H*n*Rdown - Rup*(1 -x_H)*np.exp(-CL/Tmat))*(1 + K*A2s1s_H*n*(1 -x_H)))/(Hz*(1 +z)*(1 /fu+K*A2s1s_H*n*(1 -x_H)/fu+K*Rup*n*(1 -x_H))) if (x_He < 1.e-15): f[1]=0. else: f[1] = ((x*x_He*n*Rdown_He - Rup_He*(1 -x_He)*np.exp(-CL_He/Tmat))*(1 + K_He*A2s1s_He*n_He*(1 -x_He)*He_Boltz))/(Hz*(1 +z)* (1 + K_He*(A2s1s_He+Rup_He)*n_He*(1 -x_He)*He_Boltz)) if (Heflag >= 3): f[1] = f[1]+ (x*x_He*n*Rdown_trip- (1 -x_He)*3 *Rup_trip*np.exp(-h_P*lightspeed*L_He_2st/(k_B*Tmat)))*CfHe_t/(Hz*(1 +z)) if (timeTh < H_frac*timeH): epsilon = Hz*(1 +x+fHe)/(CT*Trad**3*x) f[2] = T_CMB+ epsilon*((1 +fHe)/(1 +fHe+x))*((f[0]+fHe*f[1])/x)- epsilon* dHdz/Hz + 3.0*epsilon/(1 +z) else: f[2]= CT * (Trad**4) * x / (1 +x+fHe)* (Tmat-Trad) / (Hz*(1 +z)) + 2 *Tmat/(1 +z) return f #======================== main routine to compute RECFAST X_e(z)======================================================== Nz=1999 #number of z samples zinitial = 1.e4 zfinal = 0.1 z_arr = np.exp(np.linspace(np.log(zinitial), np.log(zfinal), Nz)) xe_arr = np.empty(Nz) z = zinitial y = np.empty(3) y[2] = T_CMB*(1. +z) #Initial rad. & mat. temperature Tmat = y[2] x_H0 = 1. x_He0 = 1. x0 = 1. +2. *fHe y[0] = x_H0 y[1] = x_He0 for i in range(Nz): zstart = z z = z_arr[i] if z > 8000: x_H0 = 1. x_He0 = 1. x0 = 1 +2 *fHe y[0] = x_H0 y[1] = x_He0 y[2] = T_CMB*(1 +z) elif z > 5000.: x_H0 = 1. x_He0 = 1. rhs = np.exp( 1.5 * np.log(CR*T_CMB/(1 +z)) - CB1_He2/(T_CMB*(1 +z)) ) / n_H_now rhs = rhs*1 #ratio of g's is 1 for He++ <-> He+ x0 = 0.5 * ( np.sqrt( (rhs-1 -fHe)**2 + 4 *(1 +2 *fHe)*rhs) - (rhs-1 -fHe) ) y[0] = x_H0 y[1] = x_He0 y[2] = T_CMB*(1 +z) elif z > 3500.: x_H0 = 1. x_He0 = 1. x0 = x_H0 + fHe*x_He0 y[0] = x_H0 y[1] = x_He0 y[2] = T_CMB*(1 +z) elif y[1] > 0.99: x_H0 = 1. rhs = np.exp( 1.5 * np.log(CR*T_CMB/(1 +z)) - CB1_He1/(T_CMB*(1 +z)) ) / n_H_now rhs = rhs*4 #ratio of g's is 4 for He+ <-> He0 x_He0 = 0.5 * ( np.sqrt( (rhs-1 )**2 + 4 *(1 +fHe)*rhs ) - (rhs-1 )) x0 = x_He0 x_He0 = (x0 - 1 )/fHe y[0] = x_H0 y[1] = x_He0 y[2] = T_CMB*(1 +z) elif y[0] > 0.99: rhs = np.exp( 1.5 * np.log(CR*T_CMB/(1 +z)) - CB1/(T_CMB*(1 +z)) ) / n_H_now x_H0 = 0.5 * (np.sqrt( rhs**2+4 *rhs ) - rhs ) newy = odeint(ion_eq, y, [zstart, z]) y = newy[1, :] y[0] = x_H0 x0 = y[0] + fHe*y[1] else: newy = odeint(ion_eq, y, [zstart, z]) y = newy[1, :] x0 = y[0] + fHe*y[1] Trad = T_CMB * (1 +z) Tmat = y[2] x_H = y[0] x_He = y[1] x = x0 xe_arr[i] = x #==============================================reionization tanh model================================================== def reion_xe(z, zre, delta_zre): return (1.+fHe)/2. * (1.+np.tanh(((1.+zre)**1.5 - (1.+z)**1.5)/(1.5*np.sqrt(1.+zre)*delta_zre))) def optical_depth_int(z, zre, delta_zre): return n_H_now*reion_xe(z, zre, delta_zre)*lightspeed*sigma_T*(1.+z)**2/recfast_Hofz(z) def optical_depth(zre, delta_zre): opt, err = quad(optical_depth_int, 0., zre + 5.*delta_zre, args = (zre, delta_zre) ) return opt def zre_from_opt(opt, delta_zre): zre_min = 0. zre_max = 25. opt_min = optical_depth(zre_min, delta_zre) opt_max = optical_depth(zre_max, delta_zre) if(opt_min > opt): return zre_min while(opt_max < opt and zre_max < 100.): zre_max += 10. opt_max = optical_depth(zre_max, delta_zre) while zre_max - zre_min > 0.01: zre_mid = (zre_max + zre_min)/2. opt_mid = optical_depth(zre_mid, delta_zre) if(opt_mid > opt): zre_max = zre_mid opt_max = opt_mid else: zre_min = zre_mid opt_min = opt_mid if opt_max - opt_min > 0.0001: return (zre_max*(opt-opt_min) + zre_min*(opt_max - opt))/(opt_max - opt_min) else: return (zre_max + zre_min)/2. if opt_re > 0.: #add reionization tanh model delta_zre = 0.5 zre = zre_from_opt(opt_re, delta_zre) for i in range(Nz): xe_arr[i] += reion_xe(z_arr[i], zre, delta_zre) #interpolate to get function X_e(z) xeofz = interp1d(z_arr, xe_arr, bounds_error = False, fill_value = "extrapolate") ##This is the main output def dopt_dlna_of_z(z): return xeofz(z)*n_H_now*(1.+z)**3*sigma_T*lightspeed/recfast_Hofz(z) #*************************************** END OF RECFAST************************************************* #**************************************BACKGROUND MODULE************************************************** #============== time unit: H_0^{-1}; length unit: c/H_0; density unit: (3H_0^2)/(8\pi G) ==================== #************************************************************************************************************ def H_of_a(a): return np.sqrt((Omega_m+Omega_r/a)/a**3+Omega_L) def dtaudsqrta(a): return 2./np.sqrt(Omega_m+Omega_L*a**3) def tau_of_a(a): #a very inaccurate routine to compute tau(a); this is only used for hierarchy truncation so should be ok acut = 0.2 if(a <= acut): return 2./np.sqrt(Omega_m)*(np.sqrt(a+a_eq)-np.sqrt(a_eq)) else: return 2./np.sqrt(Omega_m)*(np.sqrt(acut+a_eq)-np.sqrt(a_eq)) + (dtaudsqrta(acut) + 4*dtaudsqrta(acut*0.75+a*0.25)+2.*dtaudsqrta((acut+a)/2.) + 4*dtaudsqrta(acut*0.25+a*0.75) + dtaudsqrta(a) )*(np.sqrt(a)-np.sqrt(acut))/12. def tauc_of_a(a): #~collision conformal time return 1./(dopt_dlna_of_z(1./a-1.)*a*H_of_a(a)) def background_of_a(a): background = {} z = 1./a-1. rho_b = Omega_b/a**3 rho_c = Omega_c/a**3 rho_L = Omega_L rho_g = Omega_g/a**4 rho_nu = Omega_nu/a**4 rho_tot = rho_b+rho_c+rho_L+rho_g+rho_nu p_tot = -rho_L+(rho_g+rho_nu)/3. background["Omega_b"] = rho_b/rho_tot background["Omega_c"] = rho_c/rho_tot background["Omega_g"] = rho_g/rho_tot background["Omega_nu"] = rho_nu/rho_tot background["aH"] = a*np.sqrt(rho_tot) background["tau"] = tau_of_a(a) background["dotHbyHsq"] = -1.5 * ( 1. + p_tot/rho_tot) background["dopt_dlna"] = dopt_dlna_of_z(z) background["dtauc_dlna"] = (tauc_of_a(a*1.001)-tauc_of_a(a/1.001))/0.002 background["R"] = 0.75*rho_b/rho_g return background #******************************************************************************************************** #**************************************** PERTURBATION MODULE ****************************************** #******************************************************************************************************** #indice of variables lmax = 15 index_Psi = 0 index_PsiN = 1 index_delta_c = 2 index_v_c = 3 index_delta_b = 4 index_v_b = 5 index_F_nu = range(6, 7+lmax) index_F_g = range(index_F_nu[lmax]+1,index_F_nu[lmax]+lmax+2) index_E = range(index_F_g[lmax]-1,index_F_g[lmax]+lmax) num_vars = index_E[lmax]+1 def get_Phi(y, lna, k): a = np.exp(lna) bg = background_of_a(a) kbyaH = k/bg["aH"] aniso = 0.6/kbyaH**2*(bg["Omega_g"]*y[index_F_g[2]] + bg["Omega_nu"]*y[index_F_nu[2]]) return y[index_Psi] - aniso def perturb_ode(y, lna, k): #tight coupling approximation tc_cut = 50. #for k> k_suppression turn off late-time radiation and neutrino hierarchies (physically radiation and neutrinos are smooth on small scales) #this is not only for speeding up, but also for stability for large k k_suppression = 500. a = np.exp(lna) bg = background_of_a(a) aH=bg["aH"] kbyaH = k/aH doptdlna = bg["dopt_dlna"] suppr = (k/k_suppression/(1.+doptdlna))**6 #this is only significant for k>~ k_suppr and doptdlna<~1, set this to zero if you want to compare the exact scheme and the approx scheme R = bg["R"] eps = - bg["dotHbyHsq"] omg = bg["Omega_g"] omnu = bg["Omega_nu"] omb = bg["Omega_b"] omc = bg["Omega_c"] aHtau = bg["aH"]*bg["tau"] ktauc = kbyaH/doptdlna yp = np.zeros(num_vars) yp[index_Psi] = y[index_PsiN] yp[index_delta_c] = -y[index_v_c]*kbyaH + 3. * y[index_PsiN] yp[index_delta_b] = -y[index_v_b]*kbyaH + 3. * y[index_PsiN] yp[index_F_nu[0]] = -y[index_F_nu[1]]*kbyaH/3. + 4.*y[index_PsiN] yp[index_F_g[0]] = -y[index_F_g[1]]*kbyaH/3. + 4.*y[index_PsiN] aniso = 0.6/kbyaH**2*(omg*y[index_F_g[2]] + omnu*y[index_F_nu[2]]) Phi = y[index_Psi] - aniso yp[index_v_c] = -y[index_v_c] + kbyaH * Phi yp[index_F_nu[1]] = (y[index_F_nu[0]] + 4. * Phi - 0.4 * y[index_F_nu[2]]) * kbyaH for l in range(2, lmax): yp[index_F_nu[l]] = kbyaH*(l/(2*l-1.)*y[index_F_nu[l-1]] - (l+1.)/(2.*l+3.)*y[index_F_nu[l+1]]) - suppr*y[index_F_nu[l]] yp[index_F_nu[lmax]] = (lmax+0.5)/(lmax-0.5)*kbyaH*y[index_F_nu[lmax-1]] - ((lmax+1.)/aHtau + suppr)*y[index_F_nu[lmax]] if(doptdlna > tc_cut and ktauc < 0.1): #tight coupling ktaucN = k*bg["dtauc_dlna"] Uterm = - y[index_v_b] + (- y[index_F_g[0]]/4.0 + y[index_F_g[2]]/10.0)*kbyaH vterm = kbyaH * (1+R)/R + ktaucN slip = ktauc * (Uterm - ktauc*(Uterm*(-kbyaH/R + (1-eps)*ktaucN)/vterm)/vterm)/vterm # v_b - v_g accurate to (k tau_c)^2 yp[index_v_b] = -y[index_v_b] + kbyaH * Phi - slip*doptdlna/R yp[index_F_g[1]] = (y[index_F_g[0]] + 4. * Phi - 0.4 * y[index_F_g[2]]) * kbyaH + 4.*slip*doptdlna yp[index_F_g[2]] = (8./9.)*(ktauc* yp[index_F_g[1]] + y[index_F_g[1]]*ktaucN) yp[index_E[2]] = -np.sqrt(6.)/4. * yp[index_F_g[2]] else: slip = y[index_v_b] - y[index_F_g[1]]/4. yp[index_v_b] = -y[index_v_b] + kbyaH * Phi - slip*doptdlna/R yp[index_F_g[1]] = (y[index_F_g[0]] + 4. * Phi - 0.4 * y[index_F_g[2]]) * kbyaH + 4.*slip*doptdlna yp[index_F_g[2]] = kbyaH*(2./3.*y[index_F_g[1]] - 3./7.*y[index_F_g[3]])- doptdlna*(0.9* y[index_F_g[2]]+np.sqrt(6.)/10.*y[index_E[2]]) - suppr*y[index_F_g[2]] for l in range(3, lmax): yp[index_F_g[l]] = kbyaH*(l/(2*l-1.)*y[index_F_g[l-1]] - (l+1.)/(2.*l+3.)*y[index_F_g[l+1]])-(doptdlna+suppr)*y[index_F_g[l]] yp[index_F_g[lmax]] = (lmax+0.5)/(lmax-0.5)*kbyaH*y[index_F_g[lmax-1]] - ((lmax+1.)/aHtau+doptdlna+suppr)*y[index_F_g[lmax]] yp[index_E[2]] = -kbyaH* np.sqrt(5.)/7.*y[index_E[3]] - doptdlna*(0.4*y[index_E[2]] + np.sqrt(6.)/10.* y[index_F_g[2]]) - suppr*y[index_E[2]] for l in range(3, lmax): yp[index_E[l]] = kbyaH * (np.sqrt(l**2-4.)/(2*l-1.)*y[index_E[l-1]] - np.sqrt((l+1.)**2-4.)/(2*l+3.)*y[index_E[l+1]])- (doptdlna+suppr)*y[index_E[l]] yp[index_E[lmax]] = (lmax-2./lmax+0.5)/(lmax-2./lmax-0.5)*kbyaH*y[index_E[lmax-1]]-(doptdlna+(lmax-2./lmax+1.)/aHtau+suppr)*y[index_E[lmax]] G00 = -kbyaH**2/3.*y[index_Psi]-(y[index_PsiN]+Phi) T00 = 0.5*(omb*y[index_delta_b]+omc*y[index_delta_c]+omg*y[index_F_g[0]]+omnu*y[index_F_nu[0]]) yp[index_PsiN] = - 2.*y[index_Psi]-2.*(1.-eps)*Phi+0.6/kbyaH**2*(omg*yp[index_F_g[2]] + omnu*yp[index_F_nu[2]]) - (5.-eps)*y[index_PsiN] - (y[index_delta_b]*omb+y[index_delta_c]*omc)/2. - kbyaH**2/3.*(2.*y[index_Psi]-Phi) + (T00-G00)*(1.-omg-omnu) return yp def get_perturbs(k, n, a_ini, a_end): lna_list = np.linspace(np.log(a_ini), np.log(a_end), n) #-------------set adiabatic initial conditions for perturbations ------- bg_ini = background_of_a(a_ini) kbyaH_ini = k/bg_ini["aH"] Rnu = Omega_nu / Omega_r Phi_ini = -1./(1.5+0.4*Rnu) tau_ini = bg_ini["tau"] y = np.zeros(num_vars) y[index_Psi] = (1.+0.4*Rnu)*Phi_ini y[index_delta_b] = -1.5*Phi_ini y[index_delta_c] = y[index_delta_b] y[index_F_g[0]] = -2.*Phi_ini y[index_F_nu[0]] = y[index_F_g[0]] y[index_v_b] = Phi_ini/2.*(k*tau_ini) y[index_v_c] = y[index_v_b] y[index_F_nu[1]] = 4.* y[index_v_b] y[index_F_g[1]] = y[index_F_nu[1]] y[index_F_g[2]] = (8./9.)*y[index_F_g[1]]*(kbyaH_ini/bg_ini["dopt_dlna"]) y[index_F_nu[2]] = (2./3.)*Phi_ini*(k*tau_ini)**2 #yp = perturb_ode(y, lna_list[0], k) #print(yp[index_PsiN], yp[index_F_nu[0]], yp[index_F_nu[1]], y[index_F_nu[1]], yp[index_F_nu[2]], y[index_F_nu[2]]) #exit() pert = odeint(perturb_ode, y, lna_list, args = ( k, ) ) Phi_list = np.empty(n) for i in range(n): Phi_list[i] = get_Phi(pert[i, :], lna_list[i], k) return lna_list, pert, Phi_list #********************************END of PERTURBATION MODULE************************************************************************ #===============================primordial power ==================== def Primordial_D2_of_khMpc(khMpc): return A_s*(khMpc/khMpc_pivot)**(n_s-1.) #********************************Example: modify this section if you want ********************************************************** z_list = np.array( [0. ] ) #the redshifts where you want to compute the matter power spectrum nz = len(z_list) khMpc_list = np.array( [ 1.e-4, 3.e-4, 1.e-3, 2.e-3, 4.e-3, 6.e-3, 8.e-3, 9.e-3, 0.01, 0.011, 0.012, 0.014, 0.017, 0.02, 0.023, 0.026, 0.03, 0.035, 0.04, 0.045, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.18, 0.2, 0.22, 0.24, 0.26, 0.28, 0.3, 0.35, 0.4, 0.5, 0.6, 0.8, 1. ]) #wavenumbers where you want to compute the matter power spectrum; nk = len(khMpc_list) H0_in_hMpc = 1.e5/lightspeed k_list = khMpc_list / H0_in_hMpc Pm_norm = 8.*np.pi**2/9./Omega_m**2/H0_in_hMpc**4 Pm = np.empty((nz, nk)) #allocate space to save the matter power a_ini = 1.e-8 #initial a a_end = 1. n_steps = 600 #do not set this too small, otherwise the code may miss the switch between tight coupling and exact schemes for ik in range(nk): k = k_list[ik] print("solving for k/h = ", khMpc_list[ik]) lna_list, pert, Phi_list = get_perturbs(k, n_steps, a_ini, a_end) Phi_of_lna = interp1d(lna_list, Phi_list) for iz in range(nz): Pm[iz, ik] = Phi_of_lna(-np.log(1.+z_list[iz]))**2 * Primordial_D2_of_khMpc(khMpc_list[ik])*khMpc_list[ik]/(1.+z_list[iz])**2*Pm_norm for iz in range(nz): print("Matter Power at z = ", z_list[iz], ":") print("---------------------------------------") for ik in range(nk): print(khMpc_list[ik], Pm[iz, ik]) print("---------------------------------------") plt.xlabel(r'$k [h\,\mathrm{Mpc}^{-1}]$') plt.ylabel(r'$P_m [h^{-3}\mathrm{Mpc}^3]$') plt.xscale('log') plt.yscale('log') plt.plot(khMpc_list, Pm[0, :], r'--') #only plot for z[0] plt.show()