CMB tempurature spectral density

There are two ways we will compute CMB spectra in this class: pypico and Class. Class directly computes the spectral densities directly (and therefore takes some time to compute). Pypico, in contrast, fits a functional interpolation scheme to precomputed spectra for fast computation (useful when running MCMC on the LCDM parameters).

This notebook will show how to use pypico and use it to explore how the tempurature spectral density changes with respect to LCDM parameters.

Installing and wrapping pypico

The first thing you need to do to get pypico running is to follow the instructions for installation in python (here is a link).

Part of this installation requires you to download the data interpolation files for pypico. You can get the files by following the link here. Once this data file is on your computer change the following variable to point to your downloaded file.


In [1]:
picodata_path = "/Users/ethananderes/Dropbox/Courses/STA250CMB/data/pypico/pico3_tailmonty_v34.dat"


Out[1]:
"/Users/ethananderes/Dropbox/Courses/STA250CMB/data/pypico/pico3_tailmonty_v34.dat"

Now we can load the data file into pypico using PyCall


In [2]:
using PyCall
@pyimport pypico
const picoload = pypico.load_pico(picodata_path)

# --------  wrap pico
function pico(;
        omega_b     = 0.0224567, 
        omega_cdm   = 0.118489, 
        tau_reio    = 0.128312, 
        theta_s     = 0.0104098, 
        logA_s_1010 = 3.29056, 
        n_s         =  0.968602 
    )
    plout::Dict{ASCIIString, Array{Float64,1}} = picoload[:get](;
        :re_optical_depth => tau_reio,
        symbol("scalar_amp(1)") =>  1e-10*exp(logA_s_1010),
        :theta => theta_s,
        :ombh2 => omega_b,
        :omch2 => omega_cdm,
        symbol("scalar_spectral_index(1)") => n_s,
        :massive_neutrinos => 3.04,
        :helium_fraction => 0.25,
        :omnuh2 => 0.0,
        symbol("scalar_nrun(1)") => 0.0,
        :force     => true
    )
    return plout 
end


Out[2]:
pico (generic function with 1 method)

In [3]:
cls = pico()


pypico.datafiles.3e624ff23f85df8ba6355a4e18135245:326: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.
pypico.datafiles.3e624ff23f85df8ba6355a4e18135245:116: RuntimeWarning: divide by zero encountered in divide
Out[3]:
Dict{ASCIIString,Array{Float64,1}} with 15 entries:
  "scalar_TT" => [1295.07,1270.87,1213.14,1138.87,1065.0,1005.1,959.377,925.765…
  "lensed_TE" => [-0.0,1.43604,4.36515,5.82392,6.34472,6.22473,5.73343,5.05908,…
  "pk"        => [1.98537e7,1.98535e7,1.98532e7,1.9853e7,1.98527e7,1.98524e7,1.…
  "lensed_BB" => [-2.46179e-8,6.71008e-7,2.02271e-6,4.02825e-6,6.68535e-6,9.992…
  "cl_EE"     => [0.0,0.00162267,0.0768062,0.14593,0.192574,0.203976,0.183408,0…
  "lensed_EE" => [0.0,0.00162267,0.0768062,0.14593,0.192574,0.203976,0.183408,0…
  "scalar_EE" => [0.0,0.00162095,0.0768044,0.145925,0.192565,0.203965,0.183397,…
  "cl_pp"     => [-3.16208e5,4.69586e5,1.17908e6,1.81575e6,2.38305e6,2.88509e6,…
  "cl_TE"     => [-0.0,1.43604,4.36515,5.82392,6.34472,6.22473,5.73343,5.05908,…
  "k"         => [9.99e-5,0.000104633,0.00010959,0.000114782,0.00012022,0.00012…
  "cl_TT"     => [1295.05,1270.87,1213.15,1138.87,1065.0,1005.12,959.408,925.79…
  "cl_pT"     => [3550.74,22176.5,34006.0,40755.1,44139.8,45691.3,46169.2,45985…
  "scalar_TE" => [-0.0,1.43528,4.36516,5.82393,6.34474,6.22482,5.73346,5.05917,…
  "lensed_TT" => [1295.05,1270.87,1213.15,1138.87,1065.0,1005.12,959.408,925.79…
  "cl_BB"     => [-2.46179e-8,6.71008e-7,2.02271e-6,4.02825e-6,6.68535e-6,9.992…

The output of pypico is a dictionary with a bunch of spectral densities. Here are the corresponding ell values for the TT spectral density.


In [5]:
ells =  collect(0:length(cls["cl_TT"])-1)


Out[5]:
5626-element Array{Int64,1}:
    0
    1
    2
    3
    4
    5
    6
    7
    8
    9
   10
   11
   12
    ⋮
 5614
 5615
 5616
 5617
 5618
 5619
 5620
 5621
 5622
 5623
 5624
 5625

Note: the CMB spectral densities already have a factor ell(ell+1)/(2pi).

Therefore to get the raw tempurature spectral density you need to divide it out.


In [6]:
ellsc  = 2π ./ ells ./ (ells + 1) |> x -> (x[1] = 0; x);
rawClTT = cls["cl_TT"] .* ellsc


Out[6]:
5626-element Array{Float64,1}:
    0.0       
 3992.57      
 1270.41      
  596.308     
  334.58      
  210.511     
  143.527     
  103.874     
   78.6844    
   61.806     
   49.9961    
   41.4209    
   35.0209    
    ⋮         
    5.98833e-8
    5.97982e-8
    5.97133e-8
    5.96285e-8
    5.95439e-8
    5.94595e-8
    5.93753e-8
    5.92912e-8
    5.92073e-8
    5.91236e-8
    5.90401e-8
    5.89567e-8

Common scales for plotting C_l^TT


In [11]:
using PyPlot
plot(ells, rawClTT, linewidth = 4) #<---plotting the raw cls looks like a spike


Out[11]:
1-element Array{Any,1}:
 PyObject <matplotlib.lines.Line2D object at 0x31fd3ecd0>

In [14]:
semilogy(ells, rawClTT, linewidth = 2)


Out[14]:
1-element Array{Any,1}:
 PyObject <matplotlib.lines.Line2D object at 0x31faa2d50>

In [15]:
plot(ells, ells.*(ells+1).* rawClTT ./ (2π), linewidth = 2)


Out[15]:
1-element Array{Any,1}:
 PyObject <matplotlib.lines.Line2D object at 0x32d3d9b50>

In [16]:
semilogx(ells, ells.*(ells+1).* rawClTT ./ (2π), linewidth = 2)


Out[16]:
1-element Array{Any,1}:
 PyObject <matplotlib.lines.Line2D object at 0x32d768e50>

Response to varying the LCDM parameters


In [33]:
# Omega matter: omega_b ≈ 0.0223
μ = 0.0223
σ = 0.00014
l_max = 3000
for θ in linspace(μ-10σ, μ+10σ, 4)
    cl  =  pico(omega_b = θ)["cl_TT"]
    clμ =  pico(omega_b = μ)["cl_TT"]
    subplot(2,1,1)
    plot(ells[1:l_max], cl[1:l_max], label = "theta = ")
    subplot(2,1,2)
    plot(ells[1:l_max], (cl-clμ)[1:l_max])
    
end
subplot(2,1,1)
title("omega_b cl response")
legend()
subplot(2,1,2)
title("delta cl = cl - clmu")


Out[33]:
PyObject <matplotlib.text.Text object at 0x32f67fd90>

In [34]:
# Omega dark matter: omega_cdm ≈ 0.118489
μ = 0.118489
σ = 0.00014
l_max = 3000
for θ in linspace(μ-20σ, μ+20σ, 4)
    cl  =  pico(omega_cdm = θ)["cl_TT"]
    clμ =  pico(omega_cdm = μ)["cl_TT"]
    subplot(2,1,1)
    plot(ells[1:l_max], cl[1:l_max], label = "theta = ")
    subplot(2,1,2)
    plot(ells[1:l_max], (cl-clμ)[1:l_max])
    
end
subplot(2,1,1)
title("omega_cdm cl response")
legend()
subplot(2,1,2)
title("delta cl = cl - clmu")


Out[34]:
PyObject <matplotlib.text.Text object at 0x33123fa50>

In [35]:
# tau_reio ≈ 0.128312
μ = 0.118489
σ = 0.00014
l_max = 3000
for θ in linspace(μ-20σ, μ+20σ, 4)
    cl  =  pico(tau_reio = θ)["cl_TT"]
    clμ =  pico(tau_reio = μ)["cl_TT"]
    subplot(2,1,1)
    plot(ells[1:l_max], cl[1:l_max], label = "theta = ")
    subplot(2,1,2)
    plot(ells[1:l_max], (cl-clμ)[1:l_max])
    
end
subplot(2,1,1)
title("tau_reio cl response")
legend()
subplot(2,1,2)
title("delta cl = cl - clmu")


Out[35]:
PyObject <matplotlib.text.Text object at 0x3317997d0>

In [39]:
#  theta_s ≈ 0.0104098, 
μ = 0.0104098
σ = 0.0000014
l_max = 3000
for θ in linspace(μ-15σ, μ+15σ, 4)
    cl  =  pico(theta_s = θ)["cl_TT"]
    clμ =  pico(theta_s = μ)["cl_TT"]
    subplot(2,1,1)
    plot(ells[1:l_max], cl[1:l_max], label = "theta = ")
    subplot(2,1,2)
    plot(ells[1:l_max], (cl-clμ)[1:l_max])
    
end
subplot(2,1,1)
title("theta_s cl response")
legend()
subplot(2,1,2)
title("delta cl = cl - clmu")


Out[39]:
PyObject <matplotlib.text.Text object at 0x332838cd0>

In [46]:
#  logA_s_1010 ≈ 3.29056
μ = 3.29056
σ = 0.0014
l_max = 3000
for θ in linspace(μ-10σ, μ+10σ, 4)
    cl  =  pico(logA_s_1010 = θ)["cl_TT"]
    clμ =  pico(logA_s_1010 = μ)["cl_TT"]
    subplot(2,1,1)
    plot(ells[1:l_max], cl[1:l_max], label = "theta = ")
    subplot(2,1,2)
    plot(ells[1:l_max], (cl-clμ)[1:l_max])
    
end
subplot(2,1,1)
title("logA_s_1010 cl response")
legend()
subplot(2,1,2)
title("delta cl = cl - clmu")


Out[46]:
PyObject <matplotlib.text.Text object at 0x3341cfc10>

In [45]:
# n_s ≈  0.968602 
μ = 0.968602 
σ = 0.0062
l_max = 3000
for θ in linspace(μ-2σ, μ+2σ, 4)
    cl  =  pico(n_s = θ)["cl_TT"]
    clμ =  pico(n_s = μ)["cl_TT"]
    subplot(2,1,1)
    plot(ells[1:l_max], cl[1:l_max], label = "theta = ")
    subplot(2,1,2)
    plot(ells[1:l_max], (cl-clμ)[1:l_max])
    
end
subplot(2,1,1)
title("n_s cl response")
legend()
subplot(2,1,2)
title("delta cl = cl - clmu")


Out[45]:
PyObject <matplotlib.text.Text object at 0x334077590>

In [8]:
using PyPlot
plot((1:5626).^2.*cls["scalar_TT"])  # unlensed
plot((1:5626).^2.*cls["cl_TT"])      # lensed.


Out[8]:
1-element Array{Any,1}:
 PyObject <matplotlib.lines.Line2D object at 0x326453c10>

In [ ]: