#RECFAST + REIONIZATION tanh phenomenological model #original RECFAST Fortran version: see Seager, S., Sasselov, D. & Scott, D., 2000, ApJS, 128, 407 (astro-ph/9912182). #adapted for w0wa cosmology #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 ===================== #=====================USER DATA SECTION============================ Omega_c = 0.2622 Omega_b = 0.04886 Omega_k = 0. w0 = -1. wa = 0. Hubble_reduced = 0.6748 #h T_CMB = 2.726 #CMB temperature Yp = 0.25 #Helium mass abundance Neff = 3.046 #effective Neutrino species tau = 0.058 #optical depth; set to zero if do not want reionization #=================================================================== ##DO NOT CHANGE ANYTHING BELOW UNLESS YOU WANT TO MODIFY RECFAST PHYSICS=== #===================== 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_r = (blackbody_a*(1.e0+fnu)*T_CMB**4)/(3.e0*(HubbleConstant*lightspeed)**2/(8.e0*np.pi*NewtonG)) Omega_L = 1.-Omega_c-Omega_b-Omega_k-Omega_r z_eq = Omega_m/Omega_r - 1.e0 #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 #============================================================================================== ## rho_DE/rho_{DE0} as a function of z def rhoDE_ratio(z): return np.exp(3.*((1.+w0+wa)*np.log(1.+z)-wa*z/(1+z))) #Hubble parameter H = dot a / a as a function of z def Hofz(z): return HubbleConstant*np.sqrt(Omega_L*rhoDE_ratio(z) + Omega_m*(1+z)**3 + Omega_k * (1+z)**2 + Omega_r*(1+z)**4) def dHdzofz(z): return (Hofz(z+0.05)-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 = Hofz(z) dHdz = 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: tauHe_s = A2P_s*CK_He*3 *n_He*(1 -x_He)/Hz pHe_s = (1 - np.exp(-tauHe_s))/tauHe_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): tauHe_t = A2P_t*n_He*(1 -x_He)*3 tauHe_t = tauHe_t /(8 *np.pi*Hz*L_He_2Pt**(3 )) pHe_t = (1 - np.exp(-tauHe_t))/tauHe_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/Hofz(z) def optical_depth(zre, delta_zre): tau, err = quad(optical_depth_int, 0., zre + 5.*delta_zre, args = (zre, delta_zre) ) return tau def zre_from_tau(tau, delta_zre): zre_min = 0. zre_max = 25. tau_min = optical_depth(zre_min, delta_zre) tau_max = optical_depth(zre_max, delta_zre) if(tau_min > tau): return zre_min while(tau_max < tau and zre_max < 100.): zre_max += 10. tau_max = optical_depth(zre_max, delta_zre) while zre_max - zre_min > 0.01: zre_mid = (zre_max + zre_min)/2. tau_mid = optical_depth(zre_mid, delta_zre) if(tau_mid > tau): zre_max = zre_mid tau_max = tau_mid else: zre_min = zre_mid tau_min = tau_mid if tau_max - tau_min > 0.0001: return (zre_max*(tau-tau_min) + zre_min*(tau_max - tau))/(tau_max - tau_min) else: return (zre_max + zre_min)/2. if tau > 0.: #add reionization tanh model delta_zre = 0.5 zre = zre_from_tau(tau, 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") #plot z_sample = 10.**np.linspace(0., 4., 500) xe_sample = np.array([ xeofz(z) for z in z_sample]) plt.xscale('log') plt.xlabel('$z$') plt.ylabel('$X_e$') plt.plot(z_sample, xe_sample) plt.show()