Description

The code calculates \(S_\kappa(z)\) for nonzero curvation, nonzero matter with \(\Lambda=0\). Also, \(\Omega_r\) is negligible. Also, there is a code to calculate sound horizon distance \(s(z_\star)\), and the first acoustic peak multipole number of CMB anisotropy \(l_{acoustic}\).

Formula

\[\begin{equation} \begin{split} S_{\kappa}(z) &= R\sin{\frac{r(z)}{R}} \quad (\kappa=+1) \\ &= r(z) \quad (\kappa=0) \\ &= R\sinh{\frac{r(z)}{R}} \quad (\kappa=-1) \end{split} \end{equation}\] \[\begin{equation} 1 = \Omega_m(z) + \Omega_\kappa(z) \end{equation}\] \[\begin{equation} H(z) = H_0 \{ \Omega_{m,0}(1+z)^3 + \Omega_{\kappa,0}(1+z)^2 \} ^\frac{1}{2} \end{equation}\] \[\begin{equation} R = \sqrt{|\frac{\kappa}{\Omega_{\kappa,0}}|}\frac{c}{H_0} \end{equation}\] \[\begin{equation} r(z_\star) = \int_{0}^{z_\star}\frac{cdz}{H(z)} \end{equation}\] \[\begin{equation} s(z_\star) = \int_{z_\star}^{\infty}\frac{c_sdz}{H(z)} \end{equation}\] \[\begin{equation} l_{acoustic} = 4\frac{S_\kappa(z_\star)}{s(z_\star)} \end{equation}\]

Code

#---------------------------------------
# CalculatSkz calculates comoving distance, radius of curvature, and angular factor of FRW metric.
# Inputs:
# omegam0 = density parameter of matter at reference time.
# hubbleh = small h Hubble's parameter.
# redshift = redshift Z.
# Outputs:
# x[1] = comoving distance (in Mpc)
# x[2] = radius of curvature (in Mpc)
# x[3] = angular factor of FRW metric (in Mpc)
#---------------------------------------
CalculatSkz <- function(omegam0,hubbleh,redshift){
     
     #---------------------------------------
     # Find kappa
     if (1-omegam0 == 0){
          kappa = 0
     }
     if (1-omegam0 > 0){
          kappa = -1
     }
     if (1-omegam0 < 0){
          kappa = 1
     }
     #---------------------------------------
     
     #---------------------------------------
     # Find r(z)
     integrand = function(z){(omegam0*(1+z)^3 + (1-omegam0)*(1+z)^2)^(-0.5)}
     v = integrate(integrand, lower = 0, upper = redshift)
     rz = v$value * 3e8/(100e3 * hubbleh)
     #---------------------------------------
     
     #---------------------------------------
     # If kappa != 0, find radius
     # else set radius = 0
     if (kappa != 0){
          radius = sqrt(1/abs(1-omegam0)) * (3e8/(100e3 * hubbleh))
     }
     if (kappa == 0){
          radius = 0
     }
     #---------------------------------------
     
     #---------------------------------------
     # Find Sk(z)
     if (kappa == 0){
          Skz = rz
     }
     if (kappa == 1){
          Skz = radius * sin(rz/radius)
     }
     if (kappa == -1){
          Skz = radius * sinh(rz/radius)
     }
     #---------------------------------------
         
     return(c(rz,radius,Skz))
      
}

#---------------------------------------
# CalculateSoundHorizon calculates sound horizon distance by assuming constant sound speed.
# Inputs:
# soundspeed = constant sound speed (in m/s)
# omegam0 = density parameter of matter at reference time.
# hubbleh = small h Hubble's parameter.
# redshift = redshift Z.
# Outputs:
# x = sound horizon distance (in Mpc)
#---------------------------------------
CalculatSoundHorizon <- function(soundspeed,omegam0,hubbleh,redshift){

     integrand = function(z){(omegam0*(1+z)^3 + (1-omegam0)*(1+z)^2)^(-0.5)}
     v = integrate(integrand, lower = redshift, upper = Inf)
     x = v$value * soundspeed/(100e3 * hubbleh)
     return(c(x))
      
}

Test case

This test case is from part e of the problem.

omegam0 = 1
hubbleh = 1
redshift = 1000
soundspeed = 3e8/sqrt(3)
x = CalculatSkz(omegam0,hubbleh,redshift)
cat('Result from CalculateSkz: ',x)
## Result from CalculateSkz:  5810.358 0 5810.358
s = CalculatSoundHorizon(soundspeed,omegam0,hubbleh,redshift)
cat('Sound horizon (in Mpc) = ',s)
## Sound horizon (in Mpc) =  109.4897
cat('l acoustic = ',4*x[3]/s)
## l acoustic =  212.2704

part g

Part g of the problem asks to find \(l_{acoustic}\) for a \(\Lambda=0\) universe with \(\Omega_{m,0}=0.2,0.4,1\), using the code from part f. And use \(h=0.7\).

omegam0 = c(0.2,0.4,1)
hubbleh = 0.7
redshift = 1000
soundspeed = 3e8/sqrt(3)
result = c()
for (i in omegam0){
     x = CalculatSkz(i,hubbleh,redshift)
     s = CalculatSoundHorizon(soundspeed,i,hubbleh,redshift)
     result = c(result,4*x[3]/s)
}
print(result)
## [1] 431.8690 319.9684 212.2704

As we can see, comparing between the test case where \(h=1\) and this case where \(h=0.7\) at \(\Omega_{m,0}=1\), the result is insensitive to \(h\). Comparing among this case, the result is sensitive to \(\Omega_{m,0}\). We also visualize the result in the graph below.

ngrid = 100
omegam0 = c(1:ngrid)/ngrid
hubbleh = 0.7
redshift = 1000
soundspeed = 3e8/sqrt(3)
result = c()
for (i in omegam0){
     x = CalculatSkz(i,hubbleh,redshift)
     s = CalculatSoundHorizon(soundspeed,i,hubbleh,redshift)
     result = c(result,4*x[3]/s)
}
plot(omegam0,result,xlab = 'Omega_m,0',ylab = 'l_acoustic',type = 'l')