**This notebook makes use of the MaTeX package for displaying plot labels.
Run the 1st cell below to install MaTeX from Wolfram’s resource cloud; if MaTeX already installed, run the 2nd cell. Run the 3rd cell to apply Ghostscript processing when running MaTeX.

Load MaTeX package

(* Install MaTeX if not installed *)
ResourceFunction["MaTeXInstall"][]

(* If MaTeX installed, load package *)
<<MaTeX`

(* Configure MaTeX using Ghostscript -- must set path of gswin64c.exe *)
ghostscriptPath="C:\\Program Files\\gs\\gs9.50\\bin\\gswin64c.exe"
ConfigureMaTeX["Ghostscript"->ghostscriptPath]

Functions

Tools

(*Function to join piecewise interpolating functions*)
JoinInterpolatingFunction[intervals_List,flist_List]:=Module[{getGrid},getGrid[f_InterpolatingFunction,min_?NumericQ,max_?NumericQ]:={{min,f[min]}}~Join~(Transpose@{f["Grid"]//Flatten,f["ValuesOnGrid"]}//Select[#,(min<#[[1]]<max)&]&)~Join~{{max,f[max]}};
Interpolation[Table[getGrid[flist[[i]],intervals[[i]],intervals[[i+1]]],{i,Length@flist}]//Flatten[#,1]&//DeleteDuplicates[#,(#1[[1]]==#2[[1]])&]&,InterpolationOrder->2]]

Naive vacuum state functions

(*Function to solve ODE*)
CfuncNaive[kval_,χcval_,prec_]:=Block[{χ0=-EllipticK[1/2]/2,χ1=EllipticK[1/2]/2,k,χc,sol,solc,α,deC1,deC2,deC1c,deC2c,initialC1,initialC2,f1,f2,c1χc,c2χc,f1c,f2c,prec2,intervals,C1,C2},
k=SetPrecision[kval,Infinity];
χc=SetPrecision[χcval,Infinity];
α[χ_]:=JacobiCN[2χ,1/2]/Sqrt[1-Sqrt[2]JacobiSN[2χ,1/2]JacobiDN[2χ,1/2]];
deC1[kk_]=c1'[χ]==-I/kk*((χ1-χ)^(-4)-α[χ]^2(χ1-χ)^(-2))*(c1[χ]*((χ1-χ)^2+1/kk^2)+c2[χ]*Exp[2*I*kk*(χ-χ1)]*(χ1-χ-I/kk)^2);
deC2[kk_]=c2'[χ]==I/kk*((χ1-χ)^(-4)-α[χ]^2(χ1-χ)^(-2))*(c1[χ]*Exp[-2*I*kk*(χ-χ1)](χ1-χ+I/kk)^2+c2[χ]*((χ1-χ)^2+1/kk^2));
deC1c[kk_]=c1'[χ]==-I/kk*(1/5-(χ1-χ)^4/75+2/4875*(χ1-χ)^8)*(c1[χ]*((χ1-χ)^2+1/kk^2)+c2[χ]*Exp[2*I*kk*(χ-χ1)]*(χ1-χ-I/kk)^2);
deC2c[kk_]=c2'[χ]==I/kk*(1/5-(χ1-χ)^4/75+2/4875*(χ1-χ)^8)*(c1[χ]*Exp[-2*I*kk*(χ-χ1)](χ1-χ+I/kk)^2+c2[χ]*((χ1-χ)^2+1/kk^2));
initialC1[kk_]=c1[χ0]==Exp[I*kk*(χ0-χ1)]*(1-I/kk/(χ1-χ0)-1/2/kk^2/(χ1-χ0)^2);
initialC2[kk_]=c2[χ0]==-Exp[-I*kk*(χ0-χ1)]/2/kk^2/(χ1-χ0)^2;
sol=NDSolve[
{deC1[k],deC2[k],initialC1[k],initialC2[k]},
{c1,c2},
{χ,χ0,χc},
AccuracyGoal->10,
WorkingPrecision->prec];
f1=c1/.sol[[1,1]];
f2=c2/.sol[[1,2]];
c1χc=SetPrecision[f1[χc],Infinity];
c2χc=SetPrecision[f2[χc],Infinity];
solc=NDSolve[
{deC1c[k],deC2c[k],c1[χc]==c1χc,c2[χc]==c2χc},
{c1,c2},
{χ,χc,χ1},
AccuracyGoal->10,
WorkingPrecision->prec];
f1c=c1/.solc[[1,1]];
f2c=c2/.solc[[1,2]];
intervals={χ0,χc,χ1};
C1=JoinInterpolatingFunction[intervals,{f1,f1c}];
C2=JoinInterpolatingFunction[intervals,{f2,f2c}];
{C1,C2}]
(*Function to calculate final C2 for a given k*)
CfinalNaive[kval_,χcval_,prec_]:=Block[{sol,χ1=EllipticK[1/2]/2},
sol=CfuncNaive[kval,χcval,prec];
{sol[[1]][χ1],sol[[2]][χ1]}]
(*Function to calculate power spectrum*)
PphiNaive[k_,prec_]:=Block[{χc=0.9,cf},
cf=CfinalNaive[k,χc,prec];
Abs[cf[[1]]-cf[[2]]]^2]
(*Function to build interpolting function for power spectrum*)
PinterpNaive[k_List]:=Block[{ptable,prec=20},
ptable=ParallelTable[{k[[i]],PphiNaive[k[[i]],prec]},{i,1,Length@k}];
Interpolation[ptable,InterpolationOrder->1]]

Adiabatic (WKB) state functions

(*Function to solve ODE*)
CfuncWKB[kval_,χcval_,xm_,W_,dW_,order_]:=Block[{prec=13,χm,χ0=-EllipticK[1/2]/2,χ1=EllipticK[1/2]/2,k,χc,psim,dpsim,v,dv,sol,solc,α,deC1,deC2,deC1c,deC2c,initialC1,initialC2,f1,f2,c1χc,c2χc,f1c,f2c,prec2,intervals,C1,C2},
k=SetPrecision[kval,Infinity];
χm=SetPrecision[xm,Infinity];
χc=SetPrecision[χcval,Infinity];
deC1[kk_]=c1'[χ]==-I/kk*((χ1-χ)^(-4)-α[χ]^2(χ1-χ)^(-2))*(c1[χ]*((χ1-χ)^2+1/kk^2)+c2[χ]*Exp[2*I*kk*(χ-χ1)]*(χ1-χ-I/kk)^2);
deC2[kk_]=c2'[χ]==I/kk*((χ1-χ)^(-4)-α[χ]^2(χ1-χ)^(-2))*(c1[χ]*Exp[-2*I*kk*(χ-χ1)](χ1-χ+I/kk)^2+c2[χ]*((χ1-χ)^2+1/kk^2));
deC1c[kk_]=c1'[χ]==-I/kk*(1/5-(χ1-χ)^4/75+2/4875*(χ1-χ)^8)*(c1[χ]*((χ1-χ)^2+1/kk^2)+c2[χ]*Exp[2*I*kk*(χ-χ1)]*(χ1-χ-I/kk)^2);
deC2c[kk_]=c2'[χ]==I/kk*(1/5-(χ1-χ)^4/75+2/4875*(χ1-χ)^8)*(c1[χ]*Exp[-2*I*kk*(χ-χ1)](χ1-χ+I/kk)^2+c2[χ]*((χ1-χ)^2+1/kk^2));
α[χ_]:=JacobiCN[2χ,1/2]/Sqrt[1-Sqrt[2]JacobiSN[2χ,1/2]JacobiDN[2χ,1/2]];
psim[χval_]:=(1/Sqrt[2W])/.{χ->χval,e->1};
dpsim[χval_]:=(1/Sqrt[2](-I*Sqrt[W]-dW/2/(W)^(3/2)))/.{χ->χval,e->1};
initialC1[kk_]=c1[χm]==-I*Exp[I*kk*(χm-χ1)]/Sqrt[2kk](psim[χm](I*kk-1/(χm-χ1)-I/kk/(χm-χ1)^2)-dpsim[χm](1+I/kk/(χm-χ1)));
initialC2[kk_]=c2[χm]==I*Exp[-I*kk*(χm-χ1)]/Sqrt[2kk](psim[χm](-I*kk-1/(χm-χ1)+I/kk/(χm-χ1)^2)-dpsim[χm](1-I/kk/(χm-χ1)));
sol=NDSolve[
{deC1[k],deC2[k],initialC1[k],initialC2[k]},
{c1,c2},
{χ,χm,χc},
AccuracyGoal->10,
WorkingPrecision->prec];
f1=c1/.sol[[1,1]];
f2=c2/.sol[[1,2]];
c1χc=SetPrecision[f1[χc],Infinity];
c2χc=SetPrecision[f2[χc],Infinity];
solc=NDSolve[
{deC1c[k],deC2c[k],c1[χc]==c1χc,c2[χc]==c2χc},
{c1,c2},
{χ,χc,χ1},
AccuracyGoal->10,
WorkingPrecision->prec];
f1c=c1/.solc[[1,1]];
f2c=c2/.solc[[1,2]];
intervals={χm,χc,χ1};
C1=JoinInterpolatingFunction[intervals,{f1,f1c}];
C2=JoinInterpolatingFunction[intervals,{f2,f2c}];
{C1,C2}]
(*Function to calculate final C2 for a given k*)
CfinalWKB[kval_,χcval_,xm_,W_,dW_,order_]:=Block[{sol,χ1=EllipticK[1/2]/2},
sol=CfuncWKB[kval,χcval,xm,W,dW,order];
{sol[[1]][χ1],sol[[2]][χ1]}]
(*Function to calculate power spectrum for a given k*)
PphiWKB[k_,xm_,W_,dW_,order_]:=Block[{χc=0.9,cf},
cf=CfinalWKB[k,χc,xm,W,dW,order];
Abs[cf[[1]]-cf[[2]]]^2]
(*Function to build interpolating functino for power spectrum*)
PinterpWKB[k_List,xm_,order_]:=Block[{ptable,xval,W,α,Worder,dWorder},
xval=SetPrecision[xm,Infinity];
α[χ_]:=JacobiCN[2χ,1/2]/Sqrt[1-Sqrt[2]JacobiSN[2χ,1/2]JacobiDN[2χ,1/2]];
W[kk_,0]=kk;
W[kk_,n_]:=Series[Sqrt[kk^2-2*e^2*α[χ]^2-e^2*D[W[kk,n-2],χ,χ]/W[kk,n-2]/2+3/4*e^2*D[W[kk,n-2],χ]^2/W[kk,n-2]^2],{e,0,n}]//Normal;
Worder=W[kk,order]/.χ->xval;
dWorder=(D[W[kk,order],χ])/.χ->xval;
ptable=ParallelTable[{k[[i]],PphiWKB[k[[i]],xm,(Worder/.kk->k[[i]]),(dWorder/.kk->k[[i]]),order]},{i,1,Length@k}];
Interpolation[ptable,InterpolationOrder->1]]

Fourth order adiabatic (WKB) state functions

(*Function to solve ODE*)
Cfunc4[kval_,xm_,χcval_,prec_]:=Block[{χ0=-EllipticK[1/2]/2,χ1=EllipticK[1/2]/2,k,χc,χm,sol,solc,α,deC1,deC2,deC1c,deC2c,psim,dpsim,initialC1,initialC2,f1,f2,c1χc,c2χc,f1c,f2c,prec2,intervals,C1,C2},
k=SetPrecision[kval,Infinity];
χc=SetPrecision[χcval,Infinity];
χm=SetPrecision[xm,Infinity];
α[χ_]:=JacobiCN[2χ,1/2]/Sqrt[1-Sqrt[2]JacobiSN[2χ,1/2]JacobiDN[2χ,1/2]];
deC1[kk_]=c1'[χ]==-I/kk*((χ1-χ)^(-4)-α[χ]^2(χ1-χ)^(-2))*(c1[χ]*((χ1-χ)^2+1/kk^2)+c2[χ]*Exp[2*I*kk*(χ-χ1)]*(χ1-χ-I/kk)^2);
deC2[kk_]=c2'[χ]==I/kk*((χ1-χ)^(-4)-α[χ]^2(χ1-χ)^(-2))*(c1[χ]*Exp[-2*I*kk*(χ-χ1)](χ1-χ+I/kk)^2+c2[χ]*((χ1-χ)^2+1/kk^2));
deC1c[kk_]=c1'[χ]==-I/kk*(1/5-(χ1-χ)^4/75+2/4875*(χ1-χ)^8)*(c1[χ]*((χ1-χ)^2+1/kk^2)+c2[χ]*Exp[2*I*kk*(χ-χ1)]*(χ1-χ-I/kk)^2);
deC2c[kk_]=c2'[χ]==I/kk*(1/5-(χ1-χ)^4/75+2/4875*(χ1-χ)^8)*(c1[χ]*Exp[-2*I*kk*(χ-χ1)](χ1-χ+I/kk)^2+c2[χ]*((χ1-χ)^2+1/kk^2));
psim[x_]:=SetPrecision[Exp[-I*NIntegrate[k-α[y]^2/k+(2α[y]^4+1)/2/k^3,{y,χ0,x}]]/Sqrt[2k-2α[x]^2/k+(1+2α[x]^4)/k^3],Infinity];
dpsim[x_]:=SetPrecision[-I*Exp[-I*NIntegrate[k-α[y]^2/k+(2α[y]^4+1)/2/k^3,{y,χ0,x}]]*Sqrt[k/2-α[x]^2/2/k+(1+2α[x]^4)/4/k^3],Infinity];
sol=NDSolve[
{deC1[k],deC2[k],
c1[χm]==-I*Exp[I*k*(χm-χ1)]/Sqrt[2k](psim[χm](I*k-1/(χm-χ1)-I/k/(χm-χ1)^2)-dpsim[χm](1+I/k/(χm-χ1))),
c2[χm]==I*Exp[-I*k*(χm-χ1)]/Sqrt[2k](psim[χm](-I*k-1/(χm-χ1)+I/k/(χm-χ1)^2)-dpsim[χm](1-I/k/(χm-χ1)))},
{c1,c2},
{χ,χm,χc},
AccuracyGoal->10,
WorkingPrecision->prec];
f1=c1/.sol[[1,1]];
f2=c2/.sol[[1,2]];
c1χc=SetPrecision[f1[χc],Infinity];
c2χc=SetPrecision[f2[χc],Infinity];
solc=NDSolve[
{deC1c[k],deC2c[k],c1[χc]==c1χc,c2[χc]==c2χc},
{c1,c2},
{χ,χc,χ1},
AccuracyGoal->10,
WorkingPrecision->prec];
f1c=c1/.solc[[1,1]];
f2c=c2/.solc[[1,2]];
intervals={χm,χc,χ1};
C1=JoinInterpolatingFunction[intervals,{f1,f1c}];
C2=JoinInterpolatingFunction[intervals,{f2,f2c}];
{C1,C2}]
(*Function to calculate final C2 for a given k*)
Cfinal4[kval_,xm_,χcval_,prec_]:=Block[{sol,χ1=EllipticK[1/2]/2},
sol=Cfunc4[kval,xm,χcval,prec];
{sol[[1]][χ1],sol[[2]][χ1]}]
(*Function to calculate power spectrum*)
Pphi4[k_,xm_,prec_]:=Block[{χc=0.9,cf},
cf=Cfinal4[k,xm,χc,prec];
Abs[cf[[1]]-cf[[2]]]^2]
(*Function to build interpolting function for power spectrum*)
Pinterp4[k_List,xm_]:=Block[{ptable,prec=13},
ptable=ParallelTable[{k[[i]],Pphi4[k[[i]],xm,prec]},{i,1,Length@k}];
Interpolation[ptable,InterpolationOrder->1]]

Sudden approximation functions

(* Function to get c1 and c2 with sudden approximation matching at time (Subscript[χ, m]) *)
CSudden[k_,χm_]:=Block[{xm,x1=EllipticK[1/2]/2,x0=-EllipticK[1/2]/2,s1,s2,c1,c2},
xm=SetPrecision[χm,Infinity];
s1[x_]:=c1*Exp[-I*k*(x-x1)](1-I/k/(x-x1))+c2*Exp[I*k(x-x1)](1+I/k/(x-x1));
s2[x_]:=Exp[-I*k*(x-x0)];
{c1,c2}/.Solve[{s1[xm]==s2[xm],s1'[xm]==s2'[xm]},{c1,c2},WorkingPrecision->100][[1]]]
(* Function to get complete power spectrum function using sudden approximation matching at time (Subscript[χ, m]) *)
PSudden[k_,xm_]:=Abs[CSudden[k,xm][[1]]-CSudden[k,xm][[2]]]^2

Angular power spectrum functions

Using power spectrum as a list of values

(*Function to get (Subscript[C, l]) from a given power spectrum interpolating function*)
Cl[l_,pinterp_,s_,prec_]:=Block[{func,kmin,kmax},
{kmin,kmax}=pinterp["Domain"][[1]];
func[kk_]=Piecewise[
{{pinterp[kmin]/kmin*kk,kk<kmin},
{pinterp[kk],kk≥kmin&&kk≤kmax},
{Interpolation[SetPrecision[{{kmax,1},{2*kmax,1}},prec],InterpolationOrder->1][kk],kk>kmax&&kk≤2*kmax}},
1];
2*(l^2+l)*NIntegrate[
func[s*k]*SphericalBesselJ[l,k]^2/k,{k,0,Infinity},Method->"LevinRule"]]
Cl[l_,pinterp_,klower_,s_,prec_]:=Block[{func,kmax},
kmax=pinterp["Domain"][[1,2]];
func[kk_]=Piecewise[
{{pinterp[kk],kk≤kmax},
{Interpolation[SetPrecision[{{kmax,1},{2*kmax,1}},prec],InterpolationOrder->1][kk],kk>kmax&&kk≤2*kmax}},
1];
2*(l^2+l)*NIntegrate[
func[s*k]*SphericalBesselJ[l,k]^2/k,{k,klower,Infinity},Method->"LevinRule"]]
(*Function to get Cl values in a list*)
Cllist[l_List,pinterp_,s_,prec_]:=ParallelTable[{ll,Cl[ll,pinterp,s,prec]},{ll,l}]
Cllist[l_List,pinterp_,klower_,s_,prec_]:=ParallelTable[{ll,Cl[ll,pinterp,klower,s,prec]},{ll,l}]

Using power spectrum as a pure function

(* Function to get (Subscript[C, l]) from a given power spectrum function with a complete domain (0,∞) (a.k.a. using a sudden approximation) *)
Clf[l_,pinterp_,s_,prec_]:=2*(l^2+l)*NIntegrate[
pinterp[s*k]*SphericalBesselJ[l,k]^2/k,{k,0,Infinity},Method->"LevinRule"]
(* Function to get (Subscript[C, l]) values in a list from a complete power spectrum function (a.k.a. using sudden approximation) *)
Clflist[l_List,pinterp_,s_,prec_]:=ParallelTable[{ll,Clf[ll,pinterp,s,prec]},{ll,l}]

Workspace

Initialize starting times

x1=x/.FindRoot[JacobiCN[2x,1/2]/Sqrt[1-Sqrt[2]JacobiSN[2x,1/2]JacobiDN[2x,1/2]]==10^(-1),{x,-1},WorkingPrecision->100];
k1min=Sqrt[2]*10^(-1);
x2=x/.FindRoot[JacobiCN[2x,1/2]/Sqrt[1-Sqrt[2]JacobiSN[2x,1/2]JacobiDN[2x,1/2]]==10^(0),{x,0},WorkingPrecision->100];
k2min=Sqrt[2]*10^(0);

Build power spectrum functions

Sudden approximation (-S suffixed variables)

xmS=0;
pS[k_]=PSudden[k,xmS];
svalS=0.3;
llistS=Table[l,{l,2,20,0.25}];
ClS=Clflist[llistS,pS,svalS,13];
ClmarkersS=Delete[ClS,Table[{(Flatten@Table[4l+{{2},{3},{4}},{l,0,18}])[[i]]},{i,1,3*18}]];

Fourth order adiabatic state (-4 suffixed variables)

npoints4=128;
xm4=-EllipticK[1/2]/2;
klist4=Table[Tan[n*Pi/2/npoints4],{n,1,npoints4-1}];
p4=Pinterp4[klist4,x1];
p4late=Pinterp4[klist4,x2];

sval4=0.3;
llist4=Table[l,{l,2,20,0.25}];
Cl4=Cllist[llist4,p4,k1min,sval4,13];
Clmarkers4=Delete[Cl4,Table[{(Flatten@Table[4l+{{2},{3},{4}},{l,0,18}])[[i]]},{i,1,3*18}]];

Plot[{1,p4[k],pS[k]},{k,0.5,13},PlotRange->All,PlotStyle->{Black,RGBColor[112/255,173/255,17/255],RGBColor[68/255,114/255,196/255]},
PlotLegends->{Style["Asymptotic",Black, 18],Style["Numerical Solution",Black,18],Style["Sudden Approximation",Black,18]},ImageSize->Large,FrameStyle->Directive[Black,12],
FrameLabel->{MaTeX["\\text{\\textbf{Scaled} }\\boldsymbol{k}",FontSize->18]},Frame->True]

ListPlot[{Table[1,{l,llistS}],Cl4,Clmarkers4,ClS,ClmarkersS},PlotRange->{{1,Last@llistS},{0,1.1}},Joined->{True,True,False,True,False},PlotStyle->{Black,RGBColor[112/255,173/255,17/255],Black,RGBColor[68/255,114/255,196/255],Black},ImageSize->Large,
AxesStyle->Directive[Black, 12],FrameStyle->Directive[Black,12],
FrameLabel->{MaTeX["\\boldsymbol{\\ell}",FontSize->18]},Frame->True]

Zeroth order adiabatic state (-0 suffixed variables)

npoints0=128;
xm0=-EllipticK[1/2]/2*.5;
klist0=Table[Tan[n*Pi/2/npoints0],{n,1,npoints0-1}];
p0=PinterpWKB[klist0,x1,0];
p0late=PinterpWKB[klist0,x2,0];

sval0=0.3;
llist0=Table[l,{l,2,20,0.25}];
Cl0=Cllist[llist0,p0,k1min,sval0,13];
Clmarkers0=Delete[Cl0,Table[{(Flatten@Table[4l+{{2},{3},{4}},{l,0,18}])[[i]]},{i,1,3*18}]];

Second order adiabatic state (-2 suffixed variables)

npoints2=128;
xm2=-EllipticK[1/2]/2*.5;
klist2=Table[Tan[n*Pi/2/npoints2],{n,1,npoints2-1}];
p2=PinterpWKB[klist2,x1,2];p2late=PinterpWKB[Select[klist2,#>1&],x2,2];

sval2=0.3;
llist2=Table[l,{l,2,20,0.25}];
Cl2=Cllist[llist2,p2,k1min,sval2,13];
Clmarkers2=Delete[Cl2,Table[{(Flatten@Table[4l+{{2},{3},{4}},{l,0,18}])[[i]]},{i,1,3*18}]];

Plot power spectrum

Plot Style Directives

styleBD=Directive[Black,Dotted];
styleS=Directive[Black,Dashing[Large]];
style4=Directive[Black,Normal];
style0=Directive[Black,DotDashed];
style2=Directive[Black,Dashing[{0.02,0.01}]];

Power spectrum (Subscript[P, ϕ]) for Sudden Approximation

xlabelPowerSpectrum=MaTeX["\\kappa",FontSize->20];
ylabelPowerSpectrum=MaTeX["\\frac{4\\pi^{2}}{H_{\\Lambda}^{2}}P_{\\phi}",FontSize->16];
Plot[{1,pS[k]},{k,0,10},
PlotRange->{{0,10},{0,1.2}},
PlotStyle->{Directive[Black,Dotted],Directive[Black,Normal]},Frame->True,
ImageSize->600,
FrameLabel->{xlabelPowerSpectrum,ylabelPowerSpectrum}]

Power spectrum (Subscript[P, ϕ]) comparison at α=0.1

xlabelPowerSpectrum=MaTeX["\\kappa",FontSize->20];
ylabelPowerSpectrum=MaTeX["\\frac{4\\pi^{2}}{H_{\\Lambda}^{2}}P_{\\phi}",FontSize->16];
legendPowerSpectrum=LineLegend[{"Sudden approximation","0th order","4th order"},LabelStyle->{Black,14},LegendLayout->{"Column",1}];
legendPositionPowerSpectrum={0.72,0.38};
Show[Plot[1,{k,0,10},PlotStyle->Directive[Black,Dotted]],
Plot[{pS[k],p0late[k],p4late[k]},{k,1.5,10},
PlotRange->All,
PlotStyle->{Directive[Black,Dashing[Large]],Directive[Black,DotDashed],Directive[Black,Normal]},
PlotLegends->Placed[legendPowerSpectrum,legendPositionPowerSpectrum]],
ImageSize->600,
Frame->True,
FrameStyle->Directive[Black,14],
PlotRange->{{1,10},{0.5,1.2}},
PlotRangePadding->0,
FrameLabel->{xlabelPowerSpectrum,ylabelPowerSpectrum}]

Power spectrum (Subscript[P, ϕ]) comparison at α=1

xlabelPowerSpectrum=MaTeX["\\kappa",FontSize->20];
ylabelPowerSpectrum=MaTeX["\\frac{4\\pi^{2}}{H_{\\Lambda}^{2}}P_{\\phi}",FontSize->16];
legendPowerSpectrum=LineLegend[{"Sudden approximation","0th order","4th order"},LabelStyle->{Black,14},LegendLayout->{"Column",1}];
legendPositionPowerSpectrum={0.72,0.38};
Show[Plot[1,{k,0,10},PlotStyle->Directive[Black,Dotted]],
Plot[{pS[k],p0late[k],p4late[k]},{k,1.5,10},
PlotRange->All,
PlotStyle->{Directive[Black,Dashing[Large]],Directive[Black,DotDashed],Directive[Black,Normal]},
PlotLegends->Placed[legendPowerSpectrum,legendPositionPowerSpectrum]],
ImageSize->600,
Frame->True,
FrameStyle->Directive[Black,14],
PlotRange->{{0,10},{0.5,1.2}},
PlotRangePadding->0,
FrameLabel->{xlabelPowerSpectrum,ylabelPowerSpectrum}]

Power spectrum (Subscript[P, ϕ]) for adiabatic states

xlabelAdPowerSpectrum=MaTeX["\\kappa",FontSize->20];
ylabelAdPowerSpectrum=MaTeX["\\frac{4\\pi^{2}}{H_{\\Lambda}^{2}}P_{\\phi}",FontSize->16];
legendAdPowerSpectrum=LineLegend[{"0th order","2nd order","4th order"},LabelStyle->{Black,14},LegendLayout->{"Column",1}];
legendAdPositionPowerSpectrum={0.72,0.4};
Show[Plot[1,{k,0,10},PlotStyle->styleBD],
Plot[{p0[k],p2[k],p4[k]},{k,0.2,10},
PlotRange->All,
PlotStyle->{style0,style2,style4},
PlotLegends->Placed[legendAdPowerSpectrum,legendAdPositionPowerSpectrum]],
PlotRange->{{0,10},{0,1.15}},
PlotRangePadding->0,
ImageSize->600,
Frame->True,
FrameStyle->Directive[Black,14],
FrameLabel->{xlabelAdPowerSpectrum,ylabelAdPowerSpectrum},
Epilog->Inset[
Plot[{1,p0[k],p2[k],p4[k]},{k,2,10},
PlotRange->{{1.8,10},Full},
ImageSize->300,
Frame->True,
FrameStyle->Directive[Black,14],
PlotStyle->{styleBD,style0,Directive[Black,Dashing[{0.04,0.02}]],style4}
],Scaled[{0.6,0.42}]]
]

xlabelAdPowerSpectrum=MaTeX["\\kappa",FontSize->20];
ylabelAdPowerSpectrum=MaTeX["\\frac{4\\pi^{2}}{H_{\\Lambda}^{2}}P_{\\phi}",FontSize->16];
legendAdPowerSpectrum=LineLegend[{"0th order","2nd order","4th order"},LabelStyle->{Black,14},LegendLayout->{"Column",1}];
legendAdPositionPowerSpectrum={0.72,0.32};
Show[Plot[1,{k,0,10},PlotStyle->styleBD],
Plot[{p0late[k],p2late[k],p4late[k]},{k,1.5,10},
PlotRange->All,
PlotStyle->{style0,style2,style4},
PlotLegends->Placed[legendAdPowerSpectrum,legendAdPositionPowerSpectrum]],
PlotRange->{{0,10},{0,1.15}},
PlotRangePadding->0,
ImageSize->600,
Frame->True,
FrameStyle->Directive[Black,14],
FrameLabel->{xlabelAdPowerSpectrum,ylabelAdPowerSpectrum},
Epilog->Inset[
Plot[{1,p0late[k],p2late[k],p4late[k]},{k,2,10},
PlotRange->{{1.8,10},Full},
ImageSize->300,
Frame->True,
FrameStyle->Directive[Black,14],
PlotStyle->{styleBD,style0,Directive[Black,Dashing[{0.04,0.02}]],style4}
],Scaled[{0.6,0.42}]]
]

Angular spectrum (Subscript[C, ℓ])

xlabelAngularSpectrum=MaTeX["\\boldsymbol{\\ell}",FontSize->18];
ylabelAngularSpectrum=MaTeX["\\frac{18\\pi}{H_{\\Lambda}^{2}}(\\ell^{2}+\\ell)C_{\\ell}",FontSize->16];
legendAngularSpectrum=LineLegend[{None,"Sudden approximation","0th order","4th order"},LabelStyle->{Black,14},LegendLayout->{"Column",1}];
legendPositionAngularSpectrum={0.72,0.38};
BDpoints={{0,1},{30,1}};
ListPlot[Join[{BDpoints},{ClS},{Cl0},{Cl4},{ClmarkersS},{Clmarkers0},{Clmarkers4}],
Joined->{True,True,True,True,False,False,False},
PlotRange->{{1.8,20.2},{0,1.1}},
ImageSize->600,
Frame->True,
FrameStyle->Directive[Black,14],
FrameLabel->{xlabelAngularSpectrum,ylabelAngularSpectrum},
FrameTicks->{{Automatic,Automatic},{Table[l,{l,2,20}],Table[{l,""},{l,2,20}]}},
PlotStyle->{Directive[Black,Dashing[Tiny]],Directive[Black,Dashing[Large]],Directive[Black,DotDashed],Directive[Black,Normal]},
PlotMarkers->Table[Graphics[],{i,1,4}]~Join~Table[{Graphics[Disk[],PlotRangePadding->1],0.045},{i,1,3}],
PlotLegends->Placed[legendAngularSpectrum,legendPositionAngularSpectrum]]

Angular spectrum (Subscript[C, ℓ]) at different s values

llist=Table[l,{l,2,20,0.25}];
Cl4Table=Table[Cllist[llist,p4,k1min,s,13],{s,{0.5,0.3,0.2,0.1,0.05}}];

Clmarkers4s=Table[Delete[Cl4Table[[c]],Table[{(Flatten@Table[4l+{{2},{3},{4}},{l,0,18}])[[i]]},{i,1,3*18}]],{c,1,Length@Cl4Table}];

xlabelAngularSpectrum4=MaTeX["\\boldsymbol{\\ell}",FontSize->18];
ylabelAngularSpectrum4=MaTeX["\\frac{18\\pi}{H_{\\Lambda}^{2}}(\\ell^{2}+\\ell)C_{\\ell}",FontSize->16];

ListPlot[Join[{{{0,1},{30,1}}},Cl4Table,Clmarkers4s],PlotStyle->Join[{Directive[Black,Dashing[Tiny]]},Table[Directive[Black,Normal],{i,1,Length@Cl4Table}],Table[Directive[Black,16],{i,1,Length@Clmarkers4s}]],
Joined->Join[Table[True,{t,1,1+Length@Cl4Table}],Table[False,{t,1,Length@Clmarkers4s}]],
PlotMarkers->Join[Table[Graphics[],{i,1,1+Length@Cl4Table}],Table[{Graphics[Disk[],PlotRangePadding->1],0.045},{i,1,Length@Clmarkers4s}]],
PlotRange->{{1.8,20.2},{0,1.1}},
ImageSize->600,
Frame->True,
FrameStyle->Directive[Black,16],
FrameTicks->{{Automatic,None},{Table[l,{l,2,20}],None}},
BaseStyle->{FontFamily->"Times"},
FrameLabel->{xlabelAngularSpectrum4,ylabelAngularSpectrum4}]