Fluctuations in the Cosmic Microwave Background

Turn off some messages

Off[Power::infy]

Off[∞::indet]

Set cosmological parameters (Planck 2018 arXiv:1807.06209)

Clear[
(* redshift 1+(Subscript[z, eq]) of matter-radiation equality *)
z1eq,
(* baryon to matter fraction fb ≡ ((Subscript[ρ, b])/((Subscript[ρ, c])+(Subscript[ρ, b]))) *)
fb,
(* neutrino to radiation fraction fν ≡ ((Subscript[ρ, ν])/((Subscript[ρ, γ])+(Subscript[ρ, ν]))) *)
fν,
(* matter (non-baryonic cold dark matter + baryons) *)
Ωm0,Ωm,
(* non-baryonic cold dark matter *)
Ωc0,Ωc,
(* baryons *)
Ωb0,Ωb,
(* radiation (photons + neutrinos) *)
Ωr0,Ωr,
(* photons *)
Ωγ0,Ωγ,
(* neutrinos *)
Ων0,Ων,
(* curvature *)
Ωm0,Ωk,
(* cosmological constant *)
ΩΛ0,ΩΛ
];

(* redshift 1+(Subscript[z, eq]) of matter-radiation equality *)
z1eq=3400;

Clear[aUnit];
(* cosmic scale factor aUnit at which a = 1 and H = 1, relative to present-day cosmic scale factor (Subscript[a, 0]) *)
(* aUnit: 1 = present, (1/z1eq) = matter-radiation equality *)
aUnit=(1/z1eq)

Clear[HUnit];
(* value of H at aUnit, that is, where a = 1 *)
HUnit=1

If[ValueQ[aUnit],
(* (Subscript[Ω, m]) in matter today (non-baryonic cold dark matter + baryons) *)
Ωm0=0.30;
(* (Subscript[Ω, b]) in baryons *)
Ωb0=.049;
(* baryon to matter fraction fb ≡ ((Subscript[ρ, b])/((Subscript[ρ, c])+(Subscript[ρ, b]))) *)
fb=((Ωb0)/(Ωm0));
(* neutrino to radiation fraction fν ≡ ((Subscript[ρ, ν])/((Subscript[ρ, γ])+(Subscript[ρ, ν]))) *)
fν=N[((6 (7/8)((((4/11)))^(4/3)))/(2+6 (7/8)((((4/11)))^(4/3))))];
(* (Subscript[Ω, r]) in radiation today (photons + neutrinos relativistic at matter-radiation equality) *)
Ωr0=((Ωm0)/z1eq);
(* (Subscript[Ω, k]) in curvature today *)
Ωk0=0;
(* (Subscript[Ω, Λ]) in cosmological constant today *)
ΩΛ0=1-Ωr0-Ωm0-Ωk0;
(* (Subscript[Ω, c]) in cold dark matter *)
Ωc0=(1-fb)Ωm0;
(* (Subscript[Ω, γ]) in photons *)
Ωγ0=(1-fν)Ωr0;
(* (Subscript[Ω, ν]) in neutrinos *)
Ων0=fν Ωr0;
(* scale Ω's to their values at aUnit *)
Ωr=Ωr0 (aUnit^(-4));
Ωm=Ωm0 (aUnit^(-3));
Ωk=Ωk0 (aUnit^(-2));
ΩΛ=ΩΛ0;
Ωtot=Ωm+Ωr+Ωk+ΩΛ;
Ωm/=Ωtot;
Ωr/=Ωtot;
Ωk/=Ωtot;
ΩΛ/=Ωtot;
Ωtot=1;
(* (Subscript[Ω, c]) in cold dark matter *)
Ωc=(1-fb)Ωm;
(* (Subscript[Ω, c]) in baryons *)
Ωb=fb Ωm;
(* (Subscript[Ω, γ]) in photons *)
Ωγ=(1-fν)Ωr;
(* (Subscript[Ω, ν]) in neutrinos *)
Ων=fν Ωγ;
]

fν

(* values of Ω today *)
Simplify[{Ωγ0,Ων0,Ωm0,Ωc0,Ωb0,ΩΛ0,Ωk0}]

Set inflationary initial conditions

Clear[InitialConditionOption,Unspecified,Adiabatic,Isocurvature];
(* Options for InitialConditionOption are: Unspecified, Adiabatic, Isocurvature *)
InitialConditionOption=Adiabatic;

Clear[ζc,ζb,ζγ,ζν];
Switch[InitialConditionOption
,Unspecified,
,Adiabatic,
ζc=ζb=ζγ=ζν=1;
,Isocurvature,
ζc=ζb=1;
ζγ=ζν=0;
]

Choose integration variable x

Clear[x,a,lna,η,lnη];
(* Options for integration variable x are:
a = cosmic scale factor; lna = Log[cosmic scale factor]; η = conformal time; lnη = Log[Conformal time] *) 
x=lna

(* define conformal time derivative Dη ≡ (d/(dη)) *)
Clear[Dη];
Dη[expr__]:=Switch[x
,a,a aH[a](∂_a)(expr)
,lna,aH[(E^lna)](∂_lna)(expr)
,η,(∂_η)(expr)
,lnη,(E^(-lnη))(∂_(lnη))(expr)
]

(* initial value of η; note (OverDot[a]) ≡ (da/(dη)) = (a^2)H is constant in initial radiation-dominated regime *)
Clear[ηinit];
ηinit[xinit_]:=Switch[x
,a,(1/(Limit[a aH[a],a->0]))xinit
,lna,(1/(Limit[a aH[a],a->0]))(E^xinit)
,η,xinit
,lnη,(E^xinit)
]

(* initial value of a *)
Clear[ainit];
ainit[xinit_]:=Switch[x
,a,xinit
,lna,(E^xinit)
,η,Limit[a aH[a],a->0]xinit
,lnη,Limit[a aH[a],a->0](E^xinit)
]

Clear dependent variables

(* Dependent variables *)
Clear[
(* conformal time η (use ηx to avoid conflict with integration variable η) *)
ηx,
(* cosmic scale factor a (use ax to avoid conflict with integration variable a) *)
ax,
(* non-baryonic cold dark matter, subscripted c *)
(* (Subscript[δ, c]) - 3Φ *)
δcm3Φ,
(* (Subscript[v, c]) *)
vc,
(* baryons, subscriped b *)
(* (Subscript[δ, b]) - 3Φ *)
δbm3Φ,
(* (Subscript[v, b]) *)
vb,
(* photons, l'th harmonic designated (Subscript[Θ, l]) *)
(* (Subscript[Θ, 0]) - Φ *)
Θ0mΦ,
(* (Subscript[Θ, l]) *)
Θ1,Θ2,Θ,
(* neutrinos, l'th harmonic designated (Subscript[𝒩, l]) *)
(* (Subscript[𝒩, 0]) - Φ *)
𝒩0mΦ,
(* (Subscript[𝒩, l]) *)
𝒩1,𝒩2,𝒩,
(* scalar gravitational potentials Ψ and Φ *)
Ψ,Φ
]

Hubble parameter

(* Hubble parameter H = (da/(a dt)) = (da/((a^2)dη)) = ((OverDot[a])/(a^2)), scaled to H = 1 at a = 1 *)
Clear[H];
H[a_]:=HUnit(Sqrt[Ωr (a^(-4))+Ωm (a^(-3))+Ωk (a^(-2))+ΩΛ])

(* aH = a H = (da/dt) = (da/(a dη)) = ((OverDot[a])/a) *)
Clear[aH];
aH[a_]:=a H[a];

Einstein tensor

(* define F ≡ ((OverDot[a])/a)Ψ + (OverDot[Φ]) *)
Clear[F];
F[k_][x_]:=aH[ax[x]]Ψ[k][x]+Dη[Φ[k][x]];

(* Einstein tensor *)
Clear[EinsteinEnergy,EinsteinScalarMomentum,EinsteinQuadrupolePressure];
(* Einstein energy equation *)
EinsteinEnergy[k_]:=Simplify[-3aH[ax[x]]F[k][x]-(k^2)Φ[k][x],ax[x]>0];
(* Einstein momentum equation *)
EinsteinScalarMomentum[k_]:=Simplify[-k F[k][x],ax[x]>0];
(* Einstein quadrupole pressure equation *)
EinsteinQuadrupolePressure[k_]:=Simplify[(k^2)(Ψ[k][x]-Φ[k][x]),ax[x]>0];

Differential equations to solve (simple approximation)

(* Dependent variables in differential equations *)
Clear[SimpleVars];
SimpleVars[klist_]:=Flatten[{ηx,ax,Table[{δcm3Φ[k],vc[k],Θ0mΦ[k],Θ1[k],Ψ[k],Φ[k]},{k,klist}]}]

(* Einstein equations *)
Clear[SimpleEnergy,SimpleScalarMomentum,SimpleQuadrupolePressure];
(* energy equation *)
SimpleEnergy[k_]:=Simplify[(3/2)((ax[x])^2)(HUnit^2)(4Ωr ((ax[x])^(-4))(Θ0mΦ[k][x]+Φ[k][x])+Ωm ((ax[x])^(-3))(δcm3Φ[k][x]+3Φ[k][x])),ax[x]>0];
(* Einstein momentum equation *)
SimpleScalarMomentum[k_]:=Simplify[(3/2)((ax[x])^2)(HUnit^2)(4Ωr ((ax[x])^(-4))Θ1[k][x]+Ωm ((ax[x])^(-3))vc[k][x]),ax[x]>0];
(* Einstein quadrupole pressure equation *)
SimpleQuadrupolePressure[k_]:=0;

(* introduce artificial damping coefficient to damp radiation when k (Subscript[η, s]) is large *)
DampingRate[k_][x_]:=(10^(-4))k ((ηx[x])/(Sqrt[3]));

(* Define differential equations to solve *)
Clear[SimpleDiffEqs];
SimpleDiffEqs[klist_]:=Simplify[Flatten[{
(* ηx = η[x] is conformal time *)
Dη[ηx[x]]==1,
(* ax = a[x] is cosmic scale factor *)
Dη[ax[x]]==ax[x]aH[ax[x]],
Table[{
(* non-baryonic cold dark matter *)
Dη[δcm3Φ[k][x]]-k vc[k][x]==0,
Dη[vc[k][x]]+aH[ax[x]] vc[k][x]==-k Ψ[k][x],
(* photons *)
Dη[Θ0mΦ[k][x]]-k Θ1[k][x]==0,
Dη[Θ1[k][x]]+((2k)/3)DampingRate[k][x]Θ1[k][x]+(k/3)Θ0mΦ[k][x]==-(k/3)(Ψ[k][x]+Φ[k][x]),
(* Einstein equations *)
EinsteinEnergy[k]==SimpleEnergy[k],
EinsteinQuadrupolePressure[k]==SimpleQuadrupolePressure[k]
},{k,klist}]
}],ax[x]>0]

(* Initial conditions of dependent variables *)
Clear[SimpleInitConds];
SimpleInitConds[xinit_,klist_]:=Flatten[{
ηx[xinit]==ηinit[xinit],
ax[xinit]==ainit[xinit],
Table[{
δcm3Φ[k][xinit]==3ζc,
vc[k][xinit]==0,
Θ0mΦ[k][xinit]==ζγ,
Θ1[k][xinit]==0,
Ψ[k][xinit]==-(2/3)ζγ,
Φ[k][xinit]==-(2/3)ζγ
},{k,klist}]
}]

(* check initial conditions *)
Simplify[Simplify[Simplify[((ax[x])^2)(EinsteinEnergy[k]-SimpleEnergy[k])/.{Ωm->0,Ωc->0,Ωb->0,Ωk->0,ΩΛ->0}/.{ax[x]->a},a>0]/.a->0]/.{Θ0mΦ[k][x_]->ζγ,Ψ[k][x_]->-(2/3)ζγ,Φ[k][x_]->-(2/3)ζγ}]

Solve differential equations

(* macro to solve differential equations *)
Clear[SimpleSoln];
SimpleSoln[klist_,xinit_,xfin_]:=SimpleSoln[klist,xinit,xfin]=NDSolve[Flatten[{SimpleDiffEqs[klist],SimpleInitConds[xinit,klist]}],SimpleVars[klist],{x,xinit,xfin},MaxSteps->(10^5)][[1]]

(* macro to prepare variables to plot *)
Clear[PlotSimpleVars];
PlotSimpleVars[exprlist_,klist_,xinit_,xfin_]:=(Table[expr,{expr,exprlist},{k,klist}]/.SimpleSoln[klist,xinit,xfin])

Clear[amin,amax];
(* amax = present-day a *)
amax=(1/aUnit);
amin=(10^(-9))amax;
Clear[xmin,xmax];
xmin=Switch[x
,a,amin
,lna,Log[amin]
];
xmax=Switch[x
,a,amax
,lna,Log[amax]
];

(* (c/H0) in units (h^(-1))Mpc *)
Clear[cH0];
cH0=2997.92458;

(PlotSimpleVars[{1/H[ax[x]]},Table[(10^logk),{logk,{0}}],xmin,xmax]/.lna->Log[amax])[[1,1]]

(* distance (c/(a0 H0)) in program units *)
cH0prog=(PlotSimpleVars[{(1/(aH[ax[x]]))},Table[(10^logk),{logk,{0}}],xmin,xmax]/.lna->Log[amax])[[1]]

(* (h^(-1))Mpc in program units *)
Clear[hMpc];
hMpc=(cH0/cH0prog)

(* horizon distance ((a H η)/c) at amax in units of Hubble distance (c/H) *)
(PlotSimpleVars[{aH[ax[x]]ηx[x]},Table[(10^logk),{logk,{0}}],xmin,xmax]/.lna->Log[amax])[[1]]

(* η versus a *)
LogLogPlot[Evaluate[PlotSimpleVars[{ηx[x]},Table[(10^logk),{logk,{1}}],xmin,xmax]/.lna->Log[a]],{a,(10^(-6))amax,amax},PlotRange->All]

(* -Φ versus a *)
LogLogPlot[Evaluate[PlotSimpleVars[{-Φ[k][x]},Table[(10^logk),{logk,-2,-2}],xmin,xmax]/.lna->Log[a]],{a,(10^(-6))amax,amax}]

(* (Subscript[δ, c]) - 3Φ versus a *)
LogLogPlot[Evaluate[PlotSimpleVars[{δcm3Φ[k][x]},Table[(10^logk),{logk,-1,2}],xmin,xmax]/.lna->Log[a]],{a,(10^(-6))amax,amax}]

(* (Subscript[Θ, 0]) - Φ and -2Φ versus a *)
LogLinearPlot[Evaluate[PlotSimpleVars[{Θ0mΦ[k][x],-2Φ[k][x]},Table[(10^logk),{logk,-1,2}],xmin,xmax]/.lna->Log[a]],{a,(10^(-6))amax,amax},PlotRange->All]

(* -Φ versus k at fixed a *)
ListLogPlot[Evaluate[PlotSimpleVars[{{Log10[k],-Φ[k][x]}},Table[(10^logk),{logk,-2,2,.1}],xmin,xmax]/.lna->Log[a]/.a->amax],Joined->True]