cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C Integrator for Cosmic Recombination of Hydrogen and Helium, C developed by Douglas Scott (dscott@astro.ubc.ca) C based on calculations in the papers Seager, Sasselov & Scott C (ApJ, 523, L1, 1999; ApJS, 128, 407, 2000) C and "fudge" updates in Wong, Moss & Scott (2008). C C Permission to use, copy, modify and distribute without fee or royalty at C any tier, this software and its documentation, for any purpose and without C fee or royalty is hereby granted, provided that you agree to comply with C the following copyright notice and statements, including the disclaimer, C and that the same appear on ALL copies of the software and documentation, C including modifications that you make for internal use or for distribution: C C Copyright 1999-2010 by University of British Columbia. All rights reserved. C C THIS SOFTWARE IS PROVIDED "AS IS", AND U.B.C. MAKES NO C REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED. C BY WAY OF EXAMPLE, BUT NOT LIMITATION, C U.B.C. MAKES NO REPRESENTATIONS OR WARRANTIES OF C MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT C THE USE OF THE LICENSED SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE C ANY THIRD PARTY PATENTS, COPYRIGHTS, TRADEMARKS OR OTHER RIGHTS. C cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc CN Name: RECFAST CV Version: 1.5 C CP Purpose: Calculate ionised fraction as a function of redshift. CP Solves for H and He simultaneously, and includes CP H "fudge factor" for low z effect, as well as CP HeI fudge factor. C CD Description: Solves for ionisation history since recombination CD using the equations in Seager, Sasselov & Scott (ApJ, 1999). CD The Cosmological model can be flat or open. CD The matter temperature is also followed, with an update from CD Scott & Scott (2009). CD The values for \alpha_B for H are from Hummer (1994). CD The singlet HeI coefficient is a fit from the full code. CD Additional He "fudge factors" are as described in Wong, Moss CD and Scott (2008). CD Extra fitting function included (in optical depth) to account CD for extra H physics described in Rubino-Martin et al. (2010). CD Care is taken to use the most accurate constants. CD Note that some parameters are fixed (e.g. N_nu=3, nu's are CD massless, w=-1, etc.) - some users may want to explictly CD imput their own H(z) to account for extra physics. CD This is provided as a PROGRAM, which can be easily converted CD to a SUBROUTINE for use in CMB Boltzmann codes. C CA Arguments: CA Name, Description CA Double precision throughout CA CA z is redshift - W is sqrt(1+z), like conformal time CA x is total ionised fraction, relative to H CA x_H is ionized fraction of H - y(1) in R-K routine CA x_He is ionized fraction of He - y(2) in R-K routine CA (note that x_He=n_He+/n_He here and not n_He+/n_H) CA Tmat is matter temperature - y(3) in R-K routine CA f's are the derivatives of the Y's CA alphaB is case B recombination rate CA alpHe is the singlet only HeII recombination rate CA a_PPB is Pequignot, Petitjean & Boisson fitting parameter for Hydrogen CA b_PPB is Pequignot, Petitjean & Boisson fitting parameter for Hydrogen CA c_PPB is Pequignot, Petitjean & Boisson fitting parameter for Hydrogen CA d_PPB is Pequignot, Petitjean & Boisson fitting parameter for Hydrogen CA a_VF is Verner and Ferland type fitting parameter for Helium CA b_VF is Verner and Ferland type fitting parameter for Helium CA T_0 is Verner and Ferland type fitting parameter for Helium CA T_1 is Verner and Ferland type fitting parameter for Helium CA Tnow is the observed CMB temperature today CA OmegaT is the total Omega_0 CA OmegaL is the Omega_0 contribution from a Cosmological constant CA OmegaK is the Omega_0 contribution in curvature (1-O_T-O_L) CA OmegaB is Omega in baryons today CA OmegaC is the Omega_0 in (cold) dark matter: OmegaT=OmegaC+OmegaB CA Yp is the primordial helium abundace CA fHe is He/H number ratio = Yp/4(1-Yp) CA Trad and Tmat are radiation and matter temperatures CA epsilon is the approximate difference (=Trad-Tmat) at high z CA OmegaB is Omega in baryons today CA H is Hubble constant in units of 100 km/s/Mpc CA HOinp is input value of Hubble constant in units of 100 km/s/Mpc CA HO is Hubble constant in SI units CA bigH is 100 km/s/Mpc in SI units CA Hz is the value of H at the specific z (in ION) CA G is grvitational constant CA n is number density of hydrogen CA Nnow is number density today CA x0 is initial ionized fraction CA x_H0 is initial ionized fraction of Hydrogen CA x_He0 is initial ionized fraction of Helium CA rhs is dummy for calculating x0 CA zinitial and zfinal are starting and ending redshifts CA fnu is the contribution of neutrinos to the radn. energy density CA zeq is the redshift of matter-radiation equality CA zstart and zend are for each pass to the integrator CA w0 and w1 are conformal-time-like initial and final zi and zf's CA Lw0 and Lw1 are logs of w0 and w1 CA hw is the interval in W CA C,k_B,h_P: speed of light, Boltzmann's and Planck's constants CA m_e,m_H: electron mass and H atomic mass in SI CA not4: ratio of 4He atomic mass to 1H atomic mass CA sigma: Thomson cross-section CA a: radiation constant for u=aT^4 CA Pi: Pi CA Lambda: 2s-1s two photon rate for Hydrogen CA Lambda_He: 2s-1s two photon rate for Helium CA DeltaB: energy of first excited state from continuum = 3.4eV CA DeltaB_He: energy of first excited state from cont. for He = 3.4eV CA L_H_ion: level for H ionization in m^-1 CA L_H_alpha: level for H Ly alpha in m^-1 CA L_He1_ion: level for HeI ionization CA L_He2_ion: level for HeII ionization CA L_He_2s: level for HeI 2s CA L_He_2p: level for He 2p (21P1-11S0) in m^-1 CA Lalpha: Ly alpha wavelength in SI CA Lalpha_He: Helium I 2p-1s wavelength in SI CA mu_H,mu_T: mass per H atom and mass per particle CA H_frac: follow Tmat when t_Compton / t_Hubble > H_frac CA dHdz is the derivative of H at the specific z (in ION) CA CDB=DeltaB/k_B Constants derived from B1,B2,R CA CDB_He=DeltaB_He/k_B n=2-infinity for He in Kelvin CA CB1=CDB*4. Lalpha and sigma_Th, calculated CA CB1_He1: CB1 for HeI ionization potential CA CB1_He2: CB1 for HeII ionization potential CA CR=2*Pi*(m_e/h_P)*(k_B/h_P) once and passed in a common block CA CK=Lalpha**3/(8.*Pi) CA CK_He=Lalpha_He**3/(8.*Pi) CA CL=C*h_P/(k_B*Lalpha) CA CL_He=C*h_P/(k_B*Lalpha_He) CA CT=(8./3.)*(sigma/(m_e*C))*a CA Bfact=exp((E_2p-E_2s)/kT) Extra Boltzmann factor CA fu is a "fudge factor" for H, to approximate low z behaviour CA b_He is a "fudge factor" for HeI, to approximate higher z behaviour CA Heswitch is an integer for modifying HeI recombination CA Parameters and quantities to describe the extra triplet states CA and also the continuum opacity of H, with a fitting function CA suggested by KIV, astro-ph/0703438 CA a_trip: used to fit HeI triplet recombination rate CA b_trip: used to fit HeI triplet recombination rate CA L_He_2Pt: level for 23P012-11S0 in m^-1 CA L_He_2St: level for 23S1-11S0 in m^-1 CA L_He2St_ion: level for 23S1-continuum in m^-1 CA A2P_s: Einstein A coefficient for He 21P1-11S0 CA A2P_t: Einstein A coefficient for He 23P1-11S0 CA sigma_He_2Ps: H ionization x-section at HeI 21P1-11S0 freq. in m^2 CA sigma_He_2Pt: H ionization x-section at HeI 23P1-11S0 freq. in m^2 CA CL_PSt = h_P*C*(L_He_2Pt - L_He_2st)/k_B CA CfHe_t: triplet statistical correction CA Hswitch is an integer for modifying the H recombination CA AGauss1 is the amplitude of the 1st Gaussian for the H fudging CA AGauss2 is the amplitude of the 2nd Gaussian for the H fudging CA zGauss1 is the ln(1+z) central value of the 1st Gaussian CA zGauss2 is the ln(1+z) central value of the 2nd Gaussian CA wGauss1 is the width of the 1st Gaussian CA wGauss2 is the width of the 2nd Gaussian CA tol: tolerance for the integrator CA cw(24),w(3,9): work space for DVERK CA Ndim: number of d.e.'s to solve (integer) CA Nz: number of output redshitf (integer) CA I: loop index (integer) CA ind,nw: work-space for DVERK (integer) C CG Global data (common blocks) referenced: CG /zLIST/zinitial,zfinal,Nz CG /Cfund/C,k_B,h_P,m_e,m_H,not4,sigma,a,Pi CG /data/Lambda,H_frac,CB1,CDB,CR,CK,CL,CT, CG fHe,CB1_He1,CB1_He2,CDB_He,Lambda_He,Bfact,CK_He,CL_He CG /Cosmo/Tnow,HO,Nnow,z_eq,OmegaT,OmegaL,OmegaK CG /Hemod/b_He,A2P_s,A2P_t,sigma_He_2Ps,sigma_He_2Pt, CG L_He_2p,L_He_2Pt,L_He_2St,L_He2St_ion CG /Hmod/AGauss1,AGauss2,zGauss1,zGauss2,wGauss1,wGauss2 CG /Switch/Heswitch,Hswitch C CF File & device access: CF Unit /I,IO,O /Name (if known) C CM Modules called: CM DVERK (numerical integrator) CM GET_INIT (initial values for ionization fractions) CM ION (ionization and Temp derivatives) C CC Comments: CC none C CH History: CH CREATED (simplest version) 19th March 1989 CH RECREATED 11th January 1995 CH includes variable Cosmology CH uses DVERK integrator CH initial conditions are Saha CH TESTED a bunch, well, OK, not really CH MODIFIED January 1995 (include Hummer's 1994 alpha table) CH January 1995 (include new value for 2s-1s rate) CH January 1995 (expand comments) CH March 1995 (add Saha for Helium) CH August 1997 (add HeII alpha table) CH July 1998 (include OmegaT correction and H fudge factor) CH Nov 1998 (change Trad to Tmat in Rup) CH Jan 1999 (tidied up for public consumption) CH Sept 1999 (switch to formula for alpha's, fix glitch) CH Feb 2000 (fixed overflow problem in He_Boltz) CH Oct 2001 (fixed OmegaT in z_eq) CH June 2003 (fixed error in Rdown_He formula) CH June 2003 (fixed He recombination coefficients) CH June 2003 (comments to point out fixed N_nu etc.) CH Oct 2006 (included new value for G) CH Oct 2006 (improved m_He/m_H to be "not4") CH Oct 2006 (fixed error, x for x_H in part of f(1)) CH Jan 2008 (improved HeI recombination effects, CH including HeI rec. fudge factor) CH Feb 2008 (avoid calculating gamma_2Ps and CH gamma_2Pt when x_H close to 1.0) CH Aug 2008 (correction for x_H when Heflag=2 CH and Helfag>=5 to be smoother) CH Sept 2008 (added extra term to make transition CH smoother for Tmat evolution) CH Jan 2010 (added fitting function to modify K CH to match x_e(z) for new H physics) C- C =============================================================== PROGRAM recfast implicit none C --- Arguments real*8 Trad,Tmat real*8 OmegaT,OmegaB,H,HO,HOinp,bigH,G,OmegaL,OmegaK,OmegaC real*8 z,n,x,x0,rhs,x_H,x_He,x_H0,x_He0 real*8 Tnow,zinitial,zfinal,Nnow,z_eq,fnu real*8 zstart,zend,w0,w1,Lw0,Lw1,hw real*8 C,k_B,h_P,m_e,m_H,not4,sigma,a,Pi real*8 Lambda,DeltaB,DeltaB_He,Lalpha,mu_H,mu_T,H_frac real*8 Lambda_He,Lalpha_He,Bfact,CK_He,CL_He real*8 L_H_ion,L_H_alpha,L_He1_ion,L_He2_ion,L_He_2s,L_He_2p real*8 CB1,CDB,CR,CK,CL,CT,Yp,fHe,CB1_He1,CB1_He2,CDB_He,fu,b_He real*8 A2P_s,A2P_t,sigma_He_2Ps,sigma_He_2Pt real*8 L_He_2Pt,L_He_2St,L_He2St_ion real*8 AGauss1,AGauss2,zGauss1,zGauss2,wGauss1,wGauss2 real*8 tol real*8 cw(24),w(3,9) real*8 y(3) integer Ndim,Nz,I integer ind,nw integer Heswitch,Hswitch character*80 fileout C --- Parameter statements parameter(bigH=100.0D3/(1.0D6*3.0856775807D16)) !Ho in s-1 parameter(tol=1.D-5) !Tolerance for R-K external ION C --- Commons common/zLIST/zinitial,zfinal,Nz common/Cfund/C,k_B,h_P,m_e,m_H,not4,sigma,a,Pi common/Cdata/Lambda,H_frac,CB1,CDB,CR,CK,CL,CT, 1 fHe,CB1_He1,CB1_He2,CDB_He,Lambda_He,Bfact,CK_He,CL_He,fu common/Hemod/b_He,A2P_s,A2P_t,sigma_He_2Ps,sigma_He_2Pt, 1 L_He_2p,L_He_2Pt,L_He_2St,L_He2St_ion common/Hmod/AGauss1,AGauss2,zGauss1,zGauss2,wGauss1,wGauss2 common/Switch/Heswitch,Hswitch common/Cosmo/Tnow,HO,Nnow,z_eq,OmegaT,OmegaL,OmegaK C =============================================================== C --- Data data C,k_B,h_P /2.99792458D8,1.380658D-23,6.6260755D-34/ data m_e,m_H /9.1093897D-31,1.673575D-27/ !av. H atom c note: neglecting deuterium, making an O(e-5) effect data not4 /3.9715D0/ !mass He/H atom data sigma,a /6.6524616D-29,7.565914D-16/ data Pi /3.141592653589d0/ data G /6.6742D-11/ !new value C Fundamental constants in SI units C ("not4" pointed out by Gary Steigman) data Lambda /8.2245809d0/ data Lambda_He /51.3d0/ !new value from Dalgarno data L_H_ion /1.096787737D7/ !level for H ion. (in m^-1) data L_H_alpha /8.225916453D6/ !averaged over 2 levels data L_He1_ion /1.98310772D7/ !from Drake (1993) data L_He2_ion /4.389088863D7/ !from JPhysChemRefData (1987) data L_He_2s /1.66277434D7/ !from Drake (1993) data L_He_2p /1.71134891D7/ !from Drake (1993) C 2 photon rates and atomic levels in SI units data A2P_s /1.798287D9/ !Morton, Wu & Drake (2006) data A2P_t /177.58D0/ !Lach & Pachuski (2001) data L_He_2Pt /1.690871466D7/ !Drake & Morton (2007) data L_He_2St /1.5985597526D7/ !Drake & Morton (2007) data L_He2St_ion /3.8454693845D6/ !Drake & Morton (2007) data sigma_He_2Ps /1.436289D-22/ !Hummer & Storey (1998) data sigma_He_2Pt /1.484872D-22/ !Hummer & Storey (1998) C Atomic data for HeI data AGauss1 /-0.14D0/ !Amplitude of 1st Gaussian data AGauss2 /0.05D0/ !Amplitude of 2nd Gaussian data zGauss1 /7.28D0/ !ln(1+z) of 1st Gaussian data zGauss2 /6.75D0/ !ln(1+z) of 2nd Gaussian data wGauss1 /0.18D0/ !Width of 1st Gaussian data wGauss2 /0.33D0/ !Width of 2nd Gaussian C Gaussian fits for extra H physics (fit by Adam Moss) C dimensions for integrator Ndim = 3 write(*,*)'recfast version 1.5' write(*,*)'Using Hummer''s case B recombination rates for H' write(*,*)' with H fudge factor = 1.14 (or 1.105 plus high z fit),' write(*,*)' b_He fudge factor = 0.86,' write(*,*)' and a fit to tabulated HeII singlet recombination rates' write(*,*) c These are easy to inquire as input, but let's use simple values zinitial = 1.d4 z = zinitial zfinal=0.d0 c will output every 10 in z, but this is easily changed also write(*,*)'Enter output file name' read(*,'(a)')fileout write(*,*)'Enter Omega_B, Omega_DM, Omega_vac (e.g. 0.04 0.20 0.76)' read(*,*)OmegaB,OmegaC,OmegaL OmegaT=OmegaC+OmegaB !total dark matter + baryons OmegaK=1.d0-OmegaT-OmegaL !curvature write(*,'(1x,''Omega_K = '',f4.2)')OmegaK write(*,*) write(*,*)'Enter H_0 (in km/s/Mpc), T_0, Y_p (e.g. 70 2.725 0.25)' read(*,*)HOinp,Tnow,Yp c convert the Hubble constant units H = HOinp/100.d0 HO = H*bigH c sort out the helium abundance parameters mu_H = 1.d0/(1.d0-Yp) !Mass per H atom mu_T = not4/(not4-(not4-1.d0)*Yp) !Mass per atom fHe = Yp/(not4*(1.d0-Yp)) !n_He_tot / n_H_tot Nnow = 3.d0*HO*HO*OmegaB/(8.d0*Pi*G*mu_H*m_H) n = Nnow * (1.d0+z)**3 fnu = (21.d0/8.d0)*(4.d0/11.d0)**(4.d0/3.d0) c (this is explictly for 3 massless neutrinos - change if N_nu.ne.3) z_eq = (3.d0*(HO*C)**2/(8.d0*Pi*G*a*(1.d0+fnu)*Tnow**4))*OmegaT z_eq = z_eq - 1.d0 C Set up some constants so they don't have to be calculated later Lalpha = 1.d0/L_H_alpha Lalpha_He = 1.d0/L_He_2p DeltaB = h_P*C*(L_H_ion-L_H_alpha) CDB = DeltaB/k_B DeltaB_He = h_P*C*(L_He1_ion-L_He_2s) !2s, not 2p CDB_He = DeltaB_He/k_B CB1 = h_P*C*L_H_ion/k_B CB1_He1 = h_P*C*L_He1_ion/k_B !ionization for HeI CB1_He2 = h_P*C*L_He2_ion/k_B !ionization for HeII CR = 2.d0*Pi*(m_e/h_P)*(k_B/h_P) CK = Lalpha**3/(8.d0*Pi) CK_He = Lalpha_He**3/(8.d0*Pi) CL = C*h_P/(k_B*Lalpha) CL_He = C*h_P/(k_B/L_He_2s) !comes from det.bal. of 2s-1s CT = (8.d0/3.d0)*(sigma/(m_e*C))*a Bfact = h_P*C*(L_He_2p-L_He_2s)/k_B C Matter departs from radiation when t(Th) > H_frac * t(H) C choose some safely small number H_frac = 1.D-3 c Modification for H correction (Hswitch): write(*,*) 'Modification for H recombination:' write(*,*)'0) no change from old Recfast' write(*,*)'1) include correction' write(*,*)'Enter the choice of modification for H (0-1):' read(*,*)Hswitch C Fudge factor to approximate the low z out of equilibrium effect if (Hswitch .eq. 0) then fu=1.14d0 else fu=1.105d0 end if C Modification for HeI recombination (Heswitch): write(*,*)'Modification for HeI recombination:' write(*,*)'0) no change from old Recfast' write(*,*)'1) full expression for escape probability for singlet' write(*,*)' 1P-1S transition' write(*,*)'2) also including effect of contiuum opacity of H on HeI' write(*,*)' singlet (based in fitting formula suggested by' write(*,*)' Kholupenko, Ivanchik & Varshalovich, 2007)' write(*,*)'3) only including recombination through the triplets' write(*,*)'4) including 3 and the effect of the contiuum ' write(*,*)' (although this is probably negligible)' write(*,*)'5) including only 1, 2 and 3' write(*,*)'6) including all of 1 to 4' write(*,*)'Enter the choice of modification for HeI (0-6):' read(*,*)Heswitch c Set the He fudge factor cc if((Heswitch.eq. 2).or.(Heswitch.eq. 5).or.(Heswitch.eq.6))then cc write(*,*)'Enter the fudge factor b_He' cc read(*,*)b_He cc endif b_He = 0.86 c Set initial matter temperature y(3) = Tnow*(1.d0+z) !Initial rad. & mat. temperature Tmat = y(3) call get_init(z,x_H0,x_He0,x0) y(1) = x_H0 y(2) = x_He0 c OK that's the initial conditions, now start writing output file open(unit=7,status='new',form='formatted',file=fileout) write(7,'(1x,'' z '',1x,'' x_e '')') w0=1.d0/ dsqrt(1.d0 + zinitial) !like a conformal time w1=1.d0/ dsqrt(1.d0 + zfinal) Lw0 = dLog(w0) Lw1 = dLog(w1) Nz=1000 hW=(Lw1-Lw0)/dfloat(Nz) !interval in log of conf time c Set up work-space stuff for DVERK ind = 1 nw = 3 do i = 1,24 cw(i) = 0.d0 end do do i = 1,Nz C calculate the start and end redshift for the interval at each z C or just at each z zstart = zinitial + dfloat(i-1)*(zfinal-zinitial)/dfloat(Nz) zend = zinitial + dfloat(i)*(zfinal-zinitial)/dfloat(Nz) C Use Saha to get x_e, using the equation for x_e for ionized helium C and for neutral helium. C Everyb_trip ionized above z=8000. First ionization over by z=5000. C Assume He all singly ionized down to z=3500, then use He Saha until C He is 99% singly ionized, and *then* switch to joint H/He recombination. z = zend if (zend.gt.8000.d0) then x_H0 = 1.d0 x_He0 = 1.d0 x0 = 1.d0+2.d0*fHe y(1) = x_H0 y(2) = x_He0 y(3) = Tnow*(1.d0+z) else if(z.gt.5000.d0)then x_H0 = 1.d0 x_He0 = 1.d0 rhs = dexp( 1.5d0 * dLog(CR*Tnow/(1.d0+z)) 1 - CB1_He2/(Tnow*(1.d0+z)) ) / Nnow rhs = rhs*1.d0 !ratio of g's is 1 for He++ <-> He+ x0 = 0.5d0 * ( dsqrt( (rhs-1.d0-fHe)**2 1 + 4.d0*(1.d0+2.d0*fHe)*rhs) - (rhs-1.d0-fHe) ) y(1) = x_H0 y(2) = x_He0 y(3) = Tnow*(1.d0+z) else if(z.gt.3500.d0)then x_H0 = 1.d0 x_He0 = 1.d0 x0 = x_H0 + fHe*x_He0 y(1) = x_H0 y(2) = x_He0 y(3) = Tnow*(1.d0+z) else if(y(2).gt.0.99)then x_H0 = 1.d0 rhs = dexp( 1.5d0 * dLog(CR*Tnow/(1.d0+z)) 1 - CB1_He1/(Tnow*(1.d0+z)) ) / Nnow rhs = rhs*4.d0 !ratio of g's is 4 for He+ <-> He0 x_He0 = 0.5d0 * ( dsqrt( (rhs-1.d0)**2 + 4.d0*(1.d0+fHe)*rhs ) 1 - (rhs-1.d0)) x0 = x_He0 x_He0 = (x0 - 1.d0)/fHe y(1) = x_H0 y(2) = x_He0 y(3) = Tnow*(1.d0+z) else if (y(1).gt.0.99d0) then rhs = dexp( 1.5d0 * dLog(CR*Tnow/(1.d0+z)) 1 - CB1/(Tnow*(1.d0+z)) ) / Nnow x_H0 = 0.5d0 * (dsqrt( rhs**2+4.d0*rhs ) - rhs ) call DVERK(nw,ION,zstart,y,zend,tol,ind,cw,nw,w) y(1) = x_H0 x0 = y(1) + fHe*y(2) else call DVERK(nw,ION,zstart,y,zend,tol,ind,cw,nw,w) x0 = y(1) + fHe*y(2) end if Trad = Tnow * (1.d0+zend) Tmat = y(3) x_H = y(1) x_He = y(2) x = x0 write(7,'(1x,f8.2,2x,g15.8)') 1 zend,x end do stop end C =============================================================== subroutine GET_INIT(z,x_H0,x_He0,x0) C Set up the initial conditions so it will work for general, C but not pathological choices of zstart C Initial ionization fraction using Saha for relevant species implicit none real*8 OmegaT,HO,OmegaL,OmegaK real*8 z,x0,rhs,x_H0,x_He0 real*8 Tnow,Nnow,z_eq real*8 Lambda,H_frac real*8 Lambda_He,Bfact,CK_He,CL_He real*8 CB1,CDB,CR,CK,CL,CT,fHe,CB1_He1,CB1_He2,CDB_He,fu common/Cdata/Lambda,H_frac,CB1,CDB,CR,CK,CL,CT, 1 fHe,CB1_He1,CB1_He2,CDB_He,Lambda_He,Bfact,CK_He,CL_He,fu common/Cosmo/Tnow,HO,Nnow,z_eq,OmegaT,OmegaL,OmegaK C =============================================================== if(z.gt.8000.d0)then x_H0 = 1.d0 x_He0 = 1.d0 x0 = 1.d0+2.d0*fHe else if(z.gt.3500.d0)then x_H0 = 1.d0 x_He0 = 1.d0 rhs = dexp( 1.5d0 * dLog(CR*Tnow/(1.d0+z)) 1 - CB1_He2/(Tnow*(1.d0+z)) ) / Nnow rhs = rhs*1.d0 !ratio of g's is 1 for He++ <-> He+ x0 = 0.5d0 * ( dsqrt( (rhs-1.d0-fHe)**2 1 + 4.d0*(1.d0+2.d0*fHe)*rhs) - (rhs-1.d0-fHe) ) else if(z.gt.2000.d0)then x_H0 = 1.d0 rhs = dexp( 1.5d0 * dLog(CR*Tnow/(1.d0+z)) 1 - CB1_He1/(Tnow*(1.d0+z)) ) / Nnow rhs = rhs*4.d0 !ratio of g's is 4 for He+ <-> He0 x_He0 = 0.5d0 * ( dsqrt( (rhs-1.d0)**2 + 4.d0*(1.d0+fHe)*rhs ) 1 - (rhs-1.d0)) x0 = x_He0 x_He0 = (x0 - 1.d0)/fHe else rhs = dexp( 1.5d0 * dLog(CR*Tnow/(1.d0+z)) 1 - CB1/(Tnow*(1.d0+z)) ) / Nnow x_H0 = 0.5d0 * (dsqrt( rhs**2+4.d0*rhs ) - rhs ) x_He0 = 0.d0 x0 = x_H0 end if return end C =============================================================== subroutine ION(Ndim,z,Y,f) implicit none integer Ndim,Heflag,Heswitch,Hswitch real*8 z,x,n,n_He,Trad,Tmat,x_H,x_He real*8 y(Ndim),f(Ndim) real*8 C,k_B,h_P,m_e,m_H,not4,sigma,a,Pi real*8 Lambda,H_frac,Lambda_He real*8 Tnow,HO,Nnow,z_eq,Hz,OmegaT,OmegaL,OmegaK real*8 Rup,Rdown,K,K_He,Rup_He,Rdown_He,He_Boltz real*8 timeTh,timeH,factor real*8 CB1,CDB,CR,CK,CL,CT,fHe,CB1_He1,CB1_He2,CDB_He,fu,b_He real*8 Bfact,CK_He,CL_He real*8 a_VF,b_VF,T_0,T_1,sq_0,sq_1,a_PPB,b_PPB,c_PPB,d_PPB real*8 tauHe_s,pHe_s real*8 A2P_s,A2P_t,sigma_He_2Ps,sigma_He_2Pt real*8 Doppler,gamma_2Ps,pb,qb,AHcon real*8 L_He_2p,L_He_2Pt,L_He_2St,L_He2St_ion real*8 a_trip,b_trip,Rdown_trip,Rup_trip real*8 tauHe_t,pHe_t,CL_PSt,CfHe_t,gamma_2Pt real*8 AGauss1,AGauss2,zGauss1,zGauss2,wGauss1,wGauss2 real*8 dHdz,epsilon common/Cfund/C,k_B,h_P,m_e,m_H,not4,sigma,a,Pi common/Cdata/Lambda,H_frac,CB1,CDB,CR,CK,CL,CT, 1 fHe,CB1_He1,CB1_He2,CDB_He,Lambda_He,Bfact,CK_He,CL_He,fu common/Hemod/b_He,A2P_s,A2P_t,sigma_He_2Ps,sigma_He_2Pt, 1 L_He_2p,L_He_2Pt,L_He_2St,L_He2St_ion common/Hmod/AGauss1,AGauss2,zGauss1,zGauss2,wGauss1,wGauss2 common/Switch/Heswitch,Hswitch common/Cosmo/Tnow,HO,Nnow,z_eq,OmegaT,OmegaL,OmegaK C =============================================================== c the Pequignot, Petitjean & Boisson fitting parameters for Hydrogen a_PPB = 4.309d0 b_PPB = -0.6166d0 c_PPB = 0.6703d0 d_PPB = 0.5300d0 c the Verner and Ferland type fitting parameters for Helium c fixed to match those in the SSS papers, and now correct a_VF = 10.d0**(-16.744d0) b_VF = 0.711d0 T_0 = 10.d0**(0.477121d0) !3K T_1 = 10.d0**(5.114d0) c fitting parameters for HeI triplets c (matches Hummer's table with <1% error for 10^2.8 < T/K < 10^4) a_trip = 10.d0**(-16.306d0) b_trip = 0.761D0 x_H = y(1) x_He = y(2) x = x_H + fHe * x_He Tmat = y(3) n = Nnow * (1.d0+z)**3 n_He = fHe * Nnow * (1.d0+z)**3 Trad = Tnow * (1.d0+z) Hz = HO * dsqrt((1.d0+z)**4/(1.d0+z_eq)*OmegaT + OmegaT*(1.d0+z)**3 1 + OmegaK*(1.d0+z)**2 + OmegaL) c Also calculate derivative for use later dHdz = (HO**2/2.d0/Hz)*(4.d0*(1.d0+z)**3/(1.d0+z_eq)*OmegaT 1 + 3.d0*OmegaT*(1.d0+z)**2 + 2.d0*OmegaK*(1.d0+z) ) c Get the radiative rates using PPQ fit (identical to Hummer's table) Rdown=1.d-19*a_PPB*(Tmat/1.d4)**b_PPB 1 /(1.d0+c_PPB*(Tmat/1.d4)**d_PPB) Rup = Rdown * (CR*Tmat)**(1.5d0)*dexp(-CDB/Tmat) c calculate He using a fit to a Verner & Ferland type formula sq_0 = dsqrt(Tmat/T_0) sq_1 = dsqrt(Tmat/T_1) c typo here corrected by Wayne Hu and Savita Gahlaut Rdown_He = a_VF/(sq_0*(1.d0+sq_0)**(1.d0-b_VF)) Rdown_He = Rdown_He/(1.d0+sq_1)**(1.d0+b_VF) Rup_He = Rdown_He*(CR*Tmat)**(1.5d0)*dexp(-CDB_He/Tmat) Rup_He = 4.d0*Rup_He !statistical weights factor for HeI c Avoid overflow (pointed out by Jacques Roland) if((Bfact/Tmat).gt.680.d0)then He_Boltz = dexp(680.d0) else He_Boltz = dexp(Bfact/Tmat) end if c now deal with H and its fudges if (Hswitch.eq.0) then K = CK/Hz !Peebles coefficient K=lambda_a^3/8piH else c fit a double Gaussian correction function K = CK/Hz*(1.0d0 1 +AGauss1*dexp(-((log(1.0d0+z)-zGauss1)/wGauss1)**2.d0) 2 +AGauss2*dexp(-((log(1.0d0+z)-zGauss2)/wGauss2)**2.d0)) end if c add the HeI part, using same T_0 and T_1 values Rdown_trip = a_trip/(sq_0*(1.d0+sq_0)**(1.0-b_trip)) Rdown_trip = Rdown_trip/((1.d0+sq_1)**(1.d0+b_trip)) Rup_trip = Rdown_trip*dexp(-h_P*C*L_He2St_ion/(k_B*Tmat)) Rup_trip = Rup_trip*((CR*Tmat)**(1.5d0))*(4.d0/3.d0) c last factor here is the statistical weight c try to avoid "NaN" when x_He gets too small if ((x_He.lt.5.d-9) .or. (x_He.gt.0.980)) then Heflag = 0 else Heflag = Heswitch end if if (Heflag.eq.0)then !use Peebles coeff. for He K_He = CK_He/Hz else !for Heflag>0 !use Sobolev escape probability tauHe_s = A2P_s*CK_He*3.d0*n_He*(1.d0-x_He)/Hz pHe_s = (1.d0 - dexp(-tauHe_s))/tauHe_s K_He = 1.d0/(A2P_s*pHe_s*3.d0*n_He*(1.d0-x_He)) c smoother criterion here from Antony Lewis & Chad Fendt if (((Heflag.eq.2).or.(Heflag.ge.5)).and.(x_H.lt.0.9999999d0))then c use fitting formula for continuum opacity of H c first get the Doppler width parameter Doppler = 2.D0*k_B*Tmat/(m_H*not4*C*C) Doppler = C*L_He_2p*dsqrt(Doppler) gamma_2Ps = 3.d0*A2P_s*fHe*(1.d0-x_He)*C*C 1 /(dsqrt(Pi)*sigma_He_2Ps*8.d0*Pi*Doppler*(1.d0-x_H)) 2 /((C*L_He_2p)**2.d0) pb = 0.36d0 !value from KIV (2007) qb = b_He c calculate AHcon, the value of A*p_(con,H) for H continuum opacity AHcon = A2P_s/(1.d0+pb*(gamma_2Ps**qb)) K_He=1.d0/((A2P_s*pHe_s+AHcon)*3.d0*n_He*(1.d0-x_He)) end if if (Heflag.ge.3) then !include triplet effects tauHe_t = A2P_t*n_He*(1.d0-x_He)*3.d0 tauHe_t = tauHe_t /(8.d0*Pi*Hz*L_He_2Pt**(3.d0)) pHe_t = (1.d0 - dexp(-tauHe_t))/tauHe_t CL_PSt = h_P*C*(L_He_2Pt - L_He_2st)/k_B if ((Heflag.eq.3) .or. (Heflag.eq.5).or.(x_H.gt.0.99999d0)) then c no H cont. effect CfHe_t = A2P_t*pHe_t*dexp(-CL_PSt/Tmat) CfHe_t = CfHe_t/(Rup_trip+CfHe_t) !"C" factor for triplets else !include H cont. effect Doppler = 2.d0*k_B*Tmat/(m_H*not4*C*C) Doppler = C*L_He_2Pt*dsqrt(Doppler) gamma_2Pt = 3.d0*A2P_t*fHe*(1.d0-x_He)*C*C 1 /(dsqrt(Pi)*sigma_He_2Pt*8.d0*Pi*Doppler*(1.d0-x_H)) 2 /((C*L_He_2Pt)**2.d0) c use the fitting parameters from KIV (2007) in this case pb = 0.66d0 qb = 0.9d0 AHcon = A2P_t/(1.d0+pb*gamma_2Pt**qb)/3.d0 CfHe_t = (A2P_t*pHe_t+AHcon)*dexp(-CL_PSt/Tmat) CfHe_t = CfHe_t/(Rup_trip+CfHe_t) !"C" factor for triplets end if end if end if c Estimates of Thomson scattering time and Hubble time timeTh=(1.d0/(CT*Trad**4))*(1.d0+x+fHe)/x !Thomson time timeH=2.d0/(3.d0*HO*(1.d0+z)**1.5) !Hubble time c calculate the derivatives c turn on H only for x_H<0.99, and use Saha derivative for 0.98