General set-up

Remove[Evaluate[Context[]<>"*"]];

$Assumptions={k>0,t>0,(Subscript[R, ν])>0,(Subscript[R, ν])<1,t0>0,t1>0,t>0,α>0,rt>0};

mkpure[ex_]:=ex/.(f_[x_]->expr_)->(f->Function[{x},expr])

Bardeen equation (from combining all the Einstein equations)

bd1:=ϕ''[t]+4H[t]ϕ'[t]+(k^2)ϕ[t]/3==(((H[t])^2)/(k^2))(((-(k^2))/3)Π[t]-2((H[t])^2)Π[t]+H[t]Π'[t])

Subsitutions for psi and phi (from the Einstein equation), also substitution for the comoving curvature

psisub:=ψ[t]->ϕ[t]-(((H[t])^2)/(k^2))Π[t]

phisub:=ϕ[t]->ψ[t]+(((H[t])^2)/(k^2))Π[t]

zsub:=ζ[t]->ϕ[t]+(1/2)(ψ[t]+((ϕ'[t])/(H[t])))

Pre-decoupling solutions

Radiation domination subs

subs:={H[t]->1/t,Π[t]->(Subscript[R, γ])(Subscript[π, B])}

Get Bardeen equation for pre-decoupling

bd2=Simplify[bd1/.mkpure[subs]]

Solve for phi,

phisol=Flatten[DSolve[bd2,ϕ[t],t]]

Expand to get most significant contributions

ps2=ϕ[t]->Normal[Simplify[Series[ϕ[t]/.phisol,{t,0,0}]]]

Expand[ps2/.C[2]->C[2]/((k^(3/2))(Sqrt[2/π])(3^(3/4)))]

Solve for initial conditions, enforce comoving density continuous (equivalent to ϕ continuous), and comoving curvature continuous (equivalent to ϕ'+H ψ continuous).

ps3=Simplify[ps2/.Solve[{ϕ[t0]==0,ϕ'[t0]+H[t0]ψ[t0]==0}//.mkpure[Flatten[{ps2,subs,psisub}]],{C[1],C[2]}]]

Superhorizon scales get lowest order terms in k

ps4=ϕ[t]->Normal[Simplify[Series[ϕ[t]/.ps3,{k,0,0}]]]

Expand[ps4]

Solve for comoving curvature

zsol=FullSimplify[zsub//.mkpure[Flatten[{subs,ps4,psisub}]]]

Solve for synchronous gauge curvature

etaeq:=η[t]+(1/2)((η'[t])/(H[t]))==ζ[t]

etaeq2=Simplify[etaeq//.Flatten[{subs,zsol}]]

etasol=Flatten[Simplify[DSolve[{etaeq2},η[t],t]]]

et2=etasol[[1]]/.Flatten[Solve[η[t0]==0//.mkpure[etasol]]]

Check the relative divergences of the matter perturbations. This will be important for the deriving the post decoupling solutions

vsub:=V[t]->(1/2)k t(t ϕ'[t]+ψ[t])

dsub:=Δ[t]->-2((k^2)(t^2)ϕ[t]/3+t ϕ'[t]+ψ[t])

dsol=Δ[t]->Simplify[Normal[Series[Simplify[Δ[t]/.dsub/.psisub/.mkpure[subs]/.mkpure[ps4]],{k,0,0}]]]

vsol=V[t]->Normal[Series[Simplify[V[t]/.vsub/.psisub/.mkpure[subs]/.mkpure[ps4]],{k,0,0}]]

Post-decoupling

Neutrino Boltzmann equations

nueq1:=(Subscript[π, ν])'[t]+3(k/5)((3/7)(Subscript[f, 3])[t]-(2/3)(Subscript[f, 1])[t])==0

nueq2:=(Subscript[f, 1])'[t]+(k/5)((10/3)(Subscript[π, ν])[t]-5(Subscript[f, 0])[t])==4k ψ[t]

nueq3:=(Subscript[f, 0])'[t]+(k/3)(Subscript[f, 1])[t]==4ϕ'[t]

nueq4:=(Subscript[f, l])'[t]+k((l+1)/(2l+3))(Subscript[f, l+1])[t]-k (l/(2l-1)) (Subscript[f, l-1])[t]==0

Pattern for substituting our higher l derivatives

nupat1:={Derivative[n_][(Subscript[f, l_/;l>1])][t]->-k((l+1)/(2l+3))Derivative[n-1][(Subscript[f, l+1])][t]+k (l/(2l-1)) Derivative[n-1][(Subscript[f, l-1])][t],(Subscript[f, 2])[t]->5(Subscript[π, ν])[t]/3}

A few handy functions

appbs[ls_List,f_]:=Map[appbs[#,f]&,ls];
appbs[a_==b_,f_]:=f[a]==f[b];
appbs[a_->b_,f_]:=f[a]->f[b];
appbs[expr_,f_]:=f[expr];

diffexp[eq_,t_]:=appbs[eq,D[#,t]&]

mulexp[eq_,m_]:=appbs[eq,m*#&]

addexp[eq_,m_]:=appbs[eq,m+#&]

mksubj[eq_,term_]:=Flatten[Solve[eq,term]]

Post decoupling substitutions

subs2:={H[t]->1/t,Π[t]->(Subscript[R, ν]) (Subscript[π, ν])[t]+(Subscript[R, γ])(Subscript[π, B])}

Define the variable w (total anisotropic stress that we will solve for in general)

pvsub:=(Subscript[π, ν])[t]->w[t]/(Subscript[R, ν])-(Subscript[R, γ])(Subscript[π, B])/(Subscript[R, ν])

Combine neutrino equations and solve for ϕ

nueqd3=FullSimplify[diffexp[nueq1,t]//.mksubj[nueq2,(Subscript[f, 1])'[t]]//.nupat1]

Get substitution for ψ, beware of (Subscript[f, 0])term (as diverges like (k ((t))^(-2)))

pnusub=Simplify[mksubj[nueqd3,ψ[t]]]

Substitution for ϕ. Similar caveat for (Subscript[f, 0]) as with ψ

pnusub2=ϕ[t]->Simplify[ϕ[t]/.phisub//.Flatten[{subs2,pnusub}]]

Substitution for ϕ'. Remove (Subscript[f, 0])' term using Bolztmann equation. Remove (Subscript[f, l])' terms by further substitution, the rearrange to get out ϕ'.

pnusub3=mksubj[ϕ'[t]==diffexp[(ϕ[t]//.pnusub2),t]//.nupat1//.mksubj[nueq3,(Subscript[f, 0])'[t]],ϕ'[t]]

Substitution for ϕ''. Again substitute out derivatives, similar to ϕ'.

pnusub4=Simplify[diffexp[pnusub3,t]//.mksubj[nueq3,(Subscript[f, 0])'[t]]//.mksubj[nueq2,(Subscript[f, 1])'[t]]//.psisub//.mkpure[subs2]//.pnusub2//.nupat1]

Substitute everything into the Bardeen equation, giving equation in π. Then substitute to get equation in w, the total anisotropic stress.

neq1=FullSimplify[mulexp[Simplify[bd1//.Flatten[{mkpure[subs2],pnusub2,pnusub3,pnusub4}]],1]/.mkpure[pvsub]]

Check location k dependence of terms. (Subscript[f, 1])appears at (k^3), fine as it should be ~(((kt))^(-1)). (Subscript[f, 0])at order (k^4).

Simplify[Series[neq1,{k,0,4}]]

Take superhorizon limit. Should be alright as diverging terms are at high enough order.

neq2=Simplify[Normal[Series[neq1,{k,0,0}]]]

Solve for w.

ws1=FullSimplify[DSolve[neq2,w[t],t]]

Solution agrees with Kojima, other than terms C[3], and C[4]. Comes from the fact our equation is fourth order.

ns1=pvsub/.ws1

This seems to be correct (with C[3], C[4] zero), beating down the equations to actually solve for the constants takes a while, that's all below.

Fixing constants

To solve for all the constants, essentially we're making use of the fact that Δ, V, π and ϕ are all continuous through decoupling (though not necessarily smooth).

Rewrite in a slightly nicer form (redefining constants).

nt1=Simplify[PowerExpand[Simplify[(Subscript[R, ν])(Subscript[π, ν])[t]/.ns1/.(Subscript[R, ν])->(5/32)(1+(α^2))/.t->t1 x/.{C[2]->(Subscript[π, B])(Subscript[R, γ])(t1^(1/2-I α/2))C[2],C[1]->(Subscript[π, B])(Subscript[R, γ])(t1^(1/2+I α/2))C[1],C[3]->(Subscript[π, B])(Subscript[R, γ])C[3]/t1,C[4]->(Subscript[π, B])(Subscript[R, γ])C[4]/t1^2},{α>0,t1>0}]]]

ns4=(Subscript[π, ν])[t]->(nt1/(Subscript[R, ν])/.x->t/t1)

First make sure neutrino stress is zero at decoupling (time t1).

cn0=Simplify[((Subscript[π, ν])[t1]==0)/.mkpure[ns4]]

First ensure l=2 Boltzmann eq is satisfied, by matching velocity before and after decoupling

cn1=Simplify[Normal[Series[nueq1/.t->t1/.(Subscript[f, 1])[t1]->4V[t1]/.mkpure[vsol]/.mkpure[ns4],{k,0,0}]]]

Second ensure derivative combination of l=1 and l=2 Boltzmann eqs are satisfied, by matching delta before and after decoupling

cn3=FullSimplify[Normal[Simplify[Series[nueqd3/.psisub//.mkpure[subs2]/.t->t1//.mkpure[ps4]//.mkpure[ns4]//.(Subscript[f, 0])[t1]->Δ[t1]//.mkpure[dsol]/.(Subscript[R, γ])->1-(Subscript[R, ν]),{k,0,0}]]],t1>0]

Ensure that ζ is continuous at leading order.

nj1=Normal[Simplify[Series[ζ[t]/.mkpure[zsub]//.mkpure[subs]//.Flatten[{pnusub,pnusub2,pnusub3}],{k,0,0}]/.t->t1//.mkpure[ns4]]]

cn5=Simplify[Normal[Series[Simplify[mulexp[nj1==0,k^2]/.(Subscript[f, 0])[t1]->Δ[t1]//.mkpure[dsol]],{k,0,0}]]]

alsubf:=α->(Sqrt[(32/5)(Subscript[R, ν])-1])

alsubi:=(Subscript[R, ν])->(5/32)(1+(α^2))

consts=Flatten[FullSimplify[Solve[Simplify[{cn0,cn1,cn3,cn5}/.(Subscript[R, γ])->1-(Subscript[R, ν])],{C[1],C[3],C[4],C[2]}]//.alsubf/.alsubi,{t1>0,α>0}]]

ComplexExpand[consts]

FullSimplify[ComplexExpand[FullSimplify[(Subscript[π, ν])[t]/.ns4//.consts//.t->t1 Exp[x],x>0]]]

pisub1=(Subscript[π, ν])[t]->FullSimplify[ComplexExpand[FullSimplify[(Subscript[π, ν])[t]/.ns4//.consts]]]

Expand[pisub1]

Potential growth

Expand expression for the derivative of the comoving curvature. To get this at zeroth order in k, we are going to need to tackle the O((k^(-2))) terms below. Can eliminate with Bardeen equation.

zfq1=Series[Simplify[ζ'[t]//.mkpure[zsub]//.mkpure[psisub]//.mkpure[subs2]//.Flatten[{pnusub,pnusub2,pnusub3,pnusub4}]],{k,0,4}]

Get substitution for w

ws=w[t]->Simplify[((Subscript[R, ν])(Subscript[π, ν])[t]+(Subscript[R, γ])(Subscript[π, B]))//.pisub1]

Using the substituted Bardeen equation (neq1), form an expression for the fourth derivative of w.

w4sub=mksubj[Normal[Simplify[Series[neq1,{k,0,6}]]],w''''[t]]

Form differential equation for ζ from zfq1 and substitution for w4. This allows us to relate the O((k^(-2))) terms with the O(1) terms below.

zmeq1=ζ'[t]==Series[Simplify[Simplify[zfq1//.mkpure[pvsub]]//.mksubj[nueq1,(Subscript[f, 1])[t]]//.mkpure[pvsub]//.w4sub//.mksubj[nueq1,(Subscript[f, 1])[t]]],{k,0,6}]

As (Subscript[f, 0]) proportional to (k^(-2))need to be careful and substitute it out. Result has divergent terms at high enough order to ignore.

zmeq2=appbs[mulexp[diffexp[mulexp[Normal[zmeq1],1/t],t]//.mksubj[nueq3,(Subscript[f, 0])'[t]]//.pnusub3//.nupat1//.mkpure[pvsub],(t^2)],Simplify]

Substitute in solution for w.

zmeq3=Simplify[Normal[Series[zmeq2,{k,0,0}]]//.mkpure[ws]]

Solve equation for ζ.

zmsol1=DSolve[zmeq3,ζ[t],t]

Eliminate growing mode.

zmsol2=zmsol1/.C[1]->0

Fix, such that ζ is continuous at decoupling, and matches onto pre-decoupling solution.

zmsol3=Flatten[Simplify[zmsol2/.Simplify[Solve[(ζ[t1]/.Flatten[mkpure[zmsol2]])==(ζ[t1]/.mkpure[zsol]),C[2]]/.alsubf]]]

Find the asymptotic limit as t tends to infinity.

Expand[Limit[Simplify[ζ[t]/.zmsol3],t->Infinity]]

Tensor evolution

Pre-decoupling

tensoreq:=(Subscript[h, T])''[t]+2H[t](Subscript[h, T])'[t]+(k^2)(Subscript[h, T])[t]==3((H[t])^2)(Subscript[w, s])[t](Subscript[Π, T])[t]

subst:={H[t]->1/t,(Subscript[w, s])[t]->1/3,(Subscript[Π, T])[t]->(Subscript[R, γ])(Subscript[π, B])}

teq1=tensoreq//.subst

tsol1=Simplify[Flatten[DSolve[{teq1,(Subscript[h, T])[t0]==0,(Subscript[h, T])'[t0]==0},(Subscript[h, T])[t],t]]]

tsol2=(Subscript[h, T])[t]->Simplify[Normal[Simplify[Series[(Subscript[h, T])[t]/.tsol1,{k,0,0}]]]]

Post-decoupling

Neutrino equations and subsitutions

nuteq1:=(Subscript[π, ν])'[t]+2(k/(7(Sqrt[5])))(Subscript[f, 3])[t]==-(8/5)(Subscript[h, T])'[t]

nutpat1:={Derivative[n_][(Subscript[f, l_/;l>1])][t]->-k((Sqrt[(((l+1))^2)-4])/(2l+3))Derivative[n-1][(Subscript[f, l+1])][t]+k ((Sqrt[(l^2)-4])/(2l-1)) Derivative[n-1][(Subscript[f, l-1])][t],(Subscript[f, 2])[t]->5(Subscript[π, ν])[t]/2}

subst2:={H[t]->1/t,(Subscript[w, s])[t]->1/3,(Subscript[Π, T])[t]->(Subscript[R, γ])(Subscript[π, B])+(Subscript[R, ν])(Subscript[π, ν])[t]}

Substitutions for derivatives of (Subscript[h, T]).

hnsub1=mksubj[nuteq1,(Subscript[h, T])'[t]]

hnsub2=Simplify[diffexp[hnsub1,t]//.nupat1]

hnsub3=Simplify[diffexp[hnsub2,t]//.nupat1//.mksubj[nuteq1,(Subscript[π, ν])'[t]]//.hnsub1]

Form tensor equation

tpeq1=FullSimplify[diffexp[tensoreq,t]//.mkpure[subst2]//.Flatten[{hnsub1,hnsub2,hnsub3}]]

Take superhorizon limit

tpeq2=Simplify[Normal[Series[tpeq1,{k,0,0}]]]

Solve for the neutrino growth.

tpsol1=DSolve[{tpeq2,(Subscript[π, ν])[t1]==0},(Subscript[π, ν])[t],t]

Potential Growth

npteq=Normal[Series[nuteq1,{k,0,0}]]

Simplify[Flatten[DSolve[{npteq,(Subscript[h, T])[t1]==((Subscript[h, T])[t1]/.mkpure[tsol2])},(Subscript[h, T])[t],t]/.mkpure[tpsol1]]/.C[3]->0/.alsubi]

Expand[Simplify[Limit[(Subscript[h, T])[t]/.Out[178],t->Infinity]/.alsubf/.t1->rt t0]/.rt->t1/t0]