Preparations for all

Boosted Maxwellian

frMB[γ_]:=((Exp[-γ/Θ])/(Θ*BesselK[2,1/Θ]))
fBMB[γ_,βp_,Θ_,l_]:=((Exp[-γ*γp/Θ])/(γp*Θ*BesselK[2,1/Θ]))*Sqrt[(π/2)*(Θ/(pp*p))]*BesselI[1/2+l,((pp*p)/Θ)]/.pp->Sqrt[γp^2-1]/.p->Sqrt[γ^2-1]/.γp->(1/(Sqrt[1-βp^2]))

ExpfMBboost[l_,kβ_]:=Normal[Simplify[Series[fBMB[γ,βp,Θ,l]/frMB[γ],{βp,0,kβ}]/.γ->Sqrt[1+p^2],Θ>0&&p>0&&βp>0]/.(Sqrt[1+(p^2)])->γ]

ExpfMBboost[1,3]

<p^k>

Mompk[k_]:=((2*(2*Θ)^(k/2)*BesselK[(k+4)/2,1/Θ])/(Sqrt[Pi]*BesselK[2,1/Θ]))*Gamma[(k+3)/2]

Mompk[4]

ReplpM[k_]:=Reverse[Table[p^m->Simplify[Normal[Series[Mompk[m],{Θ,0,k}]]],{m,1,2*k}]]

ReplpM[4]

OO^k

OOf[k_]:=(-1)^k*Sum[StirlingS2[k,m]*DD[m],{m,1,k}]

OOf[4]

ReplOOf[k_]:=Reverse[Table[OO^m->OOf[m],{m,1,2*k}]]

ReplOOf[4]

Df[k_]:=x^k*D[(1/(Exp[x]-1)),{x,k}]

Dx^k

GenerateDxkList[kmax_]:=Table[(ExpandAll[(OO^2-3*OO)^k])==Dx[k]/.OO^n_->Ox[n],{k,1,kmax}]/.OO->Ox[1]
GenerateOOkList[kmax_]:=Table[Ox[k],{k,1,kmax}]//ExpandAll

ReplacementsOODx[kmax_]:=Module[{k,LL=GenerateDxkList[kmax],RR={},OOfun},

OOfun[n_]:=OO^n;
Dxfun[n_]:=Dx^n;

For[k=kmax,k>=1,k--,
RR=Append[RR,Solve[LL[[k]],Ox[2*k]][[1]][[1]]];
];

Return[RR/.Ox->OOfun/.Dx->Dxfun];
];

ReplacementsOODx[5]

Aberration kernel computations and Doppler operator

lmaxv=5; (* set value for kernel setup *)

+/-2 spin-harmonic generator

p2Ylm[l_,m_,θv_,ϕv_]:=Sqrt[(((l-2)!)/((l+2)!))]*(D[SphericalHarmonicY[l,m,θ,ϕ],{θ,2}]-Cot[θ]*D[SphericalHarmonicY[l,m,θ,ϕ],{θ,1}]+((2*I)/(Sin[θ]))*D[D[SphericalHarmonicY[l,m,θ,ϕ],{θ,1}]-Cot[θ]*SphericalHarmonicY[l,m,θ,ϕ],{ϕ,1}]-((D[SphericalHarmonicY[l,m,θ,ϕ],{ϕ,2}])/(Sin[θ]^2)))/.θ->θv/.ϕ->ϕv

m2Ylm[l_,m_,θv_,ϕv_]:=Sqrt[(((l-2)!)/((l+2)!))]*(D[SphericalHarmonicY[l,m,θ,ϕ],{θ,2}]-Cot[θ]*D[SphericalHarmonicY[l,m,θ,ϕ],{θ,1}]-((2*I)/(Sin[θ]))*D[D[SphericalHarmonicY[l,m,θ,ϕ],{θ,1}]-Cot[θ]*SphericalHarmonicY[l,m,θ,ϕ],{ϕ,1}]-((D[SphericalHarmonicY[l,m,θ,ϕ],{ϕ,2}])/(Sin[θ]^2)))/.θ->θv/.ϕ->ϕv

SphericalHarmonicY[2,1,θ,ϕ]
p2Ylm[2,1,θ,ϕ]//FullSimplify
m2Ylm[2,1,θ,ϕ]//FullSimplify

+/-s Y_lm

Y2mmf[m_,μp_]:=(-1)^m(2m-1)!!*Sqrt[(((2m+1))/(4*Pi*(2*m)!))*(((m-2)!)/((m+2)!))]*(((m*(m-1) (1-μp) (((1-((μp)^2)))^(m/2)))/(1+μp)))
Ym2mmf[m_,μp_]:=(-1)^m(2m-1)!!*Sqrt[(((2m+1))/(4*Pi*(2*m)!))*(((m-2)!)/((m+2)!))]*(((m*(m-1) (1+μp) (((1-((μp)^2)))^(m/2)))/(1-μp)))

(*
p2Ylm[3,3,ArcCos[μp],0]-Y2mmf[3,μp]//FullSimplify
m2Ylm[3,3,ArcCos[μp],0]-Ym2mmf[3,μp]//FullSimplify
tt=Integrate[((Y2mmf[m,μp]*Y2mmf[m,μ])/((γ*(1+β*μp))^d))/.μ->(((μp+β))/(1+β*μp))/.β->p/γ,{μp,-1,1},Assumptions->γ>p>0&&γ>1&&m∈Integers&&m>=2]
FullSimplify[tt/.γ->Sqrt[1+p^2]]
*)

specific kernel functions

K00[d_]:=(((((γ+p))^(1-d))-(((γ-p))^(1-d)))/(2 (1-d) p));

K20f[d_]:=-(((Sqrt[5]) (-(((-p+γ))^(1-d))+(((p+γ))^(1-d))))/(4 (1-d) p))-(((Sqrt[15]) (-(((Sqrt[3]) γ (((γ (-(((-p+γ))^(1-d))+(((p+γ))^(1-d))))/(2 (1-d) p))-((-(((-p+γ))^(2-d))+(((p+γ))^(2-d)))/(2 (2-d) p))))/p)+(((Sqrt[3]) (((γ (-(((-p+γ))^(2-d))+(((p+γ))^(2-d))))/(2 (2-d) p))-((-(((-p+γ))^(3-d))+(((p+γ))^(3-d)))/(2 (3-d) p))))/p)))/(2 p));
K22s[d_]:=(15/(8*(p^5)))*(-(( ((pp^(3-d))-(pm^((3-d)))))/((d-1)(d-2)(d-3)))+((2 ((pp^(1-d))-(pm^((1-d)))))/(d(d-1)(d-2)))+((((pp^(1+d))-(pm^((1+d)))))/((d+1)d(d-1))))/.pp->γ+p/.pm->γ-p

initial condition Kmm^m for general d and s=0 [Eq. (39) in Dai & Chluba 2014 + generalization]

Kmmmd1[m_]:=((Hypergeometric2F1[(m/2)+(1/2),(m/2)+1,m+(3/2),((p/γ))^2])/(γ^(m+1))) 
Kmmm[m_,d_]:=((Hypergeometric2F1[((d+m)/2),(1/2) (1+d+m),(3/2)+m,((p^2)/(γ^2))])/(γ^(d+m)))

derivation

(* Integrate[2*Pi*((2m-1)!!)^2*(((2m+1))/(4*Pi*(2*m)!))*(((1-μp^2)^(m/2)*(1-μ^2)^(m/2))/((γ*(1+β*μp))^d))/.μ->(((μp+β))/(1+β*μp))/.β->p/γ,{μp,-1,1},Assumptions->γ>p>0&&γ>1&&m∈Integers&&m>=0]*)

some checks

(* Series[Kmmm[1,d]-Kllpmf[1,1,1,d,SphericalHarmonicY]/.γ->Sqrt[1+p^2],{p,0,10}] *)

Series[Kmmm[0,d]-(((((γ+p))^(1-d))-(((γ-p))^(1-d)))/(2 (1-d) p))/.γ->Sqrt[1+p^2],{p,0,10}]

examples

FullSimplify[FunctionExpand[Kmmm[0,d]/.γ->Sqrt[1+p^2],p>0],p>0]/.(Sqrt[1+(p^2)])->γ

FullSimplify[FunctionExpand[Kmmm[1,d]/.γ->Sqrt[1+p^2],p>0],p>0]/.(Sqrt[1+(p^2)])->γ

FullSimplify[FunctionExpand[Kmmm[2,d]/.γ->Sqrt[1+p^2],p>0],p>0]/.(Sqrt[1+(p^2)])->γ

brute force initial conditions sK22^m (remember sKllp^m --> -s Ylm)

(* 
Kp2T={Kllpmf[2,2,0,d,m2Ylm],Kllpmf[2,2,1,d,m2Ylm],Kllpmf[2,2,2,d,m2Ylm]};
FullSimplify[FullSimplify[FunctionExpand[Kp2T/.γ->Sqrt[1+p^2]],p>0]/.(Sqrt[1+(p^2)])->γ]
*)

(* 2K22^0, 2K22^1, 2K22^2 *)
Kp2T={((15 (-8 (p-γ) (((p+γ))^d) (3+(3+(-2+d) d) (p^2)-3 (-1+d) p γ)-8 (((-p+γ))^d) (p+γ) (3+(3+(-2+d) d) (p^2)+3 (-1+d) p γ)))/(16 (-3+d) (-2+d) (-1+d) d (1+d) (p^5))),((5 (((1+(p^2)))^(d/2)) (((p+γ))^(1-2 d)) (((1+(p^2)-p γ))^(-d)) (6 γ+3 ((-1+d) p-2 γ) (((p-γ))^2) (((p+γ))^(2 d))+(-1+d) p (9+2 (3+(-2+d) d) (p^2)+6 (-1+d) p γ)))/(2 (-3+d) (-2+d) (-1+d) d (1+d) (p^5))),((((1+(p^2)))^((1/2) (-1-d))) (15 (((1+(p^2)))^((1+d)/2)) (((p+γ))^d)-5 (((1+(p^2)-p γ))^d) (2 d p (1+(p^2)) (3+2 (-1+(d^2)) (p^2) (1+(-2+d) (p^2)))+(3+2 d (p^2) (3 d+d (-1+(d^2)) (p^2)+2 (-2+d) (-1+d) (1+d) (p^4))) γ)))/(4 (-3+d) (-2+d) (-1+d) d (1+d) (p^5) (p+γ))};

initial condition sK22^m for general d (remember sKllp^m --> -s Ylm)

Intξ[k_]:=(((pp^k)-(pm^k))/k)

Km0R=(15/(16*p^5))*(Intξ[1+d]-4*γ*Intξ[d]+2(1+2γ^2)*Intξ[1-d]-4*γ*Intξ[2-d]+Intξ[3-d]);
Km1R=-(5/(8(p^5)))*(pm*Intξ[1+d]-(3+pm^2)*Intξ[d]+6γ*Intξ[1-d]-(3+pp^2)*Intξ[2-d]+pp*Intξ[3-d]);
Km2R=(5/(32(p^5)))*(pm^2*Intξ[1+d]-4*pm*Intξ[d]+6*Intξ[1-d]-4*pp*Intξ[2-d]+pp^2*Intξ[3-d]);

Kp2Tm={Km0R,Km1R,Km2R}/.pm->γ-p/.pp->γ+p;

Kernels22m[s_,m_]:=Module[{},
If[Abs[s]!=2||Abs[m]>2,Return[0],];
If[s==2,If[m>=0,Return[Kp2Tm[[m+1]]];,Return[Kp2Tm[[-m+1]]/.p->-p];];];
If[s==-2,If[m>=0,Return[Kp2Tm[[m+1]]/.p->-p];,Return[Kp2Tm[[-m+1]]];];];
]

derivation

Km0I=FullSimplify[2*Pi*dKllpmf[2,2,0,d,m2Ylm]/.(γ^2)->1+p^2,γ>1&&γ>p>0&&-1<μp<1]
Km1I=FullSimplify[2*Pi*dKllpmf[2,2,1,d,m2Ylm]/.(γ^2)->1+p^2,γ>1&&γ>p>0&&-1<μp<1]
Km2I=FullSimplify[2*Pi*dKllpmf[2,2,2,d,m2Ylm]/.(γ^2)->1+p^2,γ>1&&γ>p>0&&-1<μp<1]

Km0ξ=Collect[Collect[Collect[FullSimplify[Km0I/.μp->((-γ+ξ)/p)]/((15/16)(1/((p^4)*(ξ^(2+d)))))/.γ->Sqrt[1+p^2],ξ,FullSimplify],ξ,FullSimplify]/.(Sqrt[1+(p^2)])->γ/.(p^2)->γ^2-1,ξ,FullSimplify]
Km1ξ=Collect[Collect[FullSimplify[Km1I/.μp->((-γ+ξ)/p)]/(-(5/8)(1/((p^4)*(ξ^(2+d)))))/.γ->Sqrt[1+p^2],ξ,FullSimplify],ξ,FullSimplify]/.(Sqrt[1+(p^2)])->γ/.(p+γ)->pp/.(-p+γ)->pm/.(-4+2 p (-p+γ))->-2*(2-p*pm)
Km2ξ=Collect[Collect[FullSimplify[Km2I/.μp->((-γ+ξ)/p)]/((5/32)(1/((p^4)*(ξ^(2+d)))))/.γ->Sqrt[1+p^2],ξ,FullSimplify],ξ,FullSimplify]/.(Sqrt[1+(p^2)])->γ/.(p+γ)->pp/.(-p+γ)->pm/.(p-γ)->-pm/.(4 p-4 γ)->-4*pm/.(1-2 p pm)->pm^2

integral, recursion & identities

FullSimplify[(k+1)Intξ[k+1]-(2*γ*k*Intξ[k]-(k-1)*Intξ[k-1])/.pp->γ+p/.pm->γ-p/.γ->Sqrt[1+p^2]]

FullSimplify[(1-2 p pm)-pm^2/.pm->γ-p/.γ->Sqrt[1+p^2],p>0]
FullSimplify[(1+2 p pp)-pp^2/.pp->γ+p/.γ->Sqrt[1+p^2],p>0]

checks

Series[Kp2T[[1]]-Kp2Tm[[1]]/.γ->Sqrt[1+p^2],{p,0,10}]
Series[Kp2T[[2]]-Kp2Tm[[2]]/.γ->Sqrt[1+p^2],{p,0,10}]
Series[Kp2T[[3]]-Kp2Tm[[3]]/.γ->Sqrt[1+p^2],{p,0,10}]

(* 
mv=-2;
Series[Kernels22m[2,mv]-Kllpmf[2,2,mv,d,m2Ylm]/.γ->Sqrt[1+p^2],{p,0,10}]//FullSimplify
Series[Kernels22m[-2,mv]-Kllpmf[2,2,mv,d,p2Ylm]/.γ->Sqrt[1+p^2],{p,0,10}]//FullSimplify 
*)

aberration kernel elements for s=0 using recursions (up to l=3 for m=0) [brute force method JC]

Kraisef[d_,l_,Klm1f_,Klm2f_]:=-Sqrt[((4l^2-1)/(l^2))]*((γ*Klm1f[d]-Klm1f[d-1])/p)-((l-1)/l)*Sqrt[((2l+1)/(2l-3))]*Klm2f[d]
Kswitchf[d_,Klf_]:=Klf[2-d]/.p->-p

K00[d_]:=(((((γ+p))^(1-d))-(((γ-p))^(1-d)))/(2 (1-d) p));

K10[d_]:=Kraisef[d,1,K00,K00]
K01[d_]:=Kswitchf[d,K10]

K11[d_]:=Kraisef[d,1,K01,K01]
K21[d_]:=Kraisef[d,2,K11,K01]
K12[d_]:=Kswitchf[d,K21]
K22[d_]:=Kraisef[d,2,K12,K02]
K32[d_]:=Kraisef[d,3,K22,K12]

K20[d_]:=Kraisef[d,2,K10,K00]
K02[d_]:=Kswitchf[d,K20]

K30[d_]:=Kraisef[d,3,K20,K10]
K03[d_]:=Kswitchf[d,K30]

T00[d_]:=K00[d-1]/γ/.p->-p
T10[d_]:=K10[d-1]/γ/.p->-p
T20[d_]:=K20[d-1]/γ/.p->-p
T30[d_]:=K30[d-1]/γ/.p->-p

T02[d_]:=K02[d-1]/γ/.p->-p
T12[d_]:=K12[d-1]/γ/.p->-p
T22[d_]:=K22[d-1]/γ/.p->-p
T32[d_]:=K32[d-1]/γ/.p->-p

some derivations

pp==γ+p
pm==γ-p

K00[d]/.p->(pp-pm)/2/.γ->(pp+pm)/2//Simplify
K10[d]/.p->(pp-pm)/2/.γ->(pp+pm)/2//Simplify
K20[d]/.p->(pp-pm)/2/.γ->(pp+pm)/2//Simplify
K30[d]/.p->(pp-pm)/2/.γ->(pp+pm)/2//Simplify

Series[K00[d]/.γ->Sqrt[1+p^2],{p,0,2}]
Series[K10[d]/.γ->Sqrt[1+p^2],{p,0,2}]
Series[K20[d]/.γ->Sqrt[1+p^2],{p,0,2}]
Series[K30[d]/.γ->Sqrt[1+p^2],{p,0,2}]

compute kernels explicitly for given Ylm function pointer (this is for checking mainly as quite slow)

dKllpmf[l_,lp_,m_,d_,Ylm_]:=Simplify[((ComplexExpand[Conjugate[Ylm[l,m,θp,ϕ]]]*Ylm[lp,m,θ,ϕ])/((γ*(1+β*μp))^d))/.θp->ArcCos[μp]/.θ->ArcCos[μ]/.μ->(((μp+β))/(1+β*μp))/.β->p/γ]

Kllpmf[l_,lp_,m_,d_,Ylm_]:=Integrate[2*Pi*dKllpmf[l,lp,m,d,Ylm],{μp,-1,1},Assumptions->γ>p>0&&γ>1];

Kllpf[l_,lp_,d_,Ylm_]:=Table[Kllpmf[l,lp,m,d,Ylm],{m,-Min[l,lp],Min[l,lp]}];

checks (--> passed)

(* FullSimplify[Series[Kllpf[0,0,d,SphericalHarmonicY]-K00[d]/.γ->Sqrt[1+p^2],{p,0,8}],d>0&&p>0]
FullSimplify[Series[Kllpf[1,0,d,SphericalHarmonicY]-K10[d]/.γ->Sqrt[1+p^2],{p,0,8}],d>0&&p>0]
FullSimplify[Series[Kllpf[0,1,d,SphericalHarmonicY]-K01[d]/.γ->Sqrt[1+p^2],{p,0,8}],d>0&&p>0]
FullSimplify[Series[Kllpf[2,0,d,SphericalHarmonicY]-K20[d]/.γ->Sqrt[1+p^2],{p,0,8}],d>0&&p>0]
FullSimplify[Series[Kllpf[0,2,d,SphericalHarmonicY]-K02[d]/.γ->Sqrt[1+p^2],{p,0,8}],d>0&&p>0] *)

(* FullSimplify[Series[Kllpmf[2,2,0,d,p2Ylm]-K22s[d]/.γ->Sqrt[1+p^2],{p,0,8}],d>0&&p>0]
FullSimplify[Series[Kllpmf[2,2,0,d,m2Ylm]-K22s[d]/.γ->Sqrt[1+p^2],{p,0,8}],d>0&&p>0]
FullSimplify[Series[Kllpmf[3,2,0,d,p2Ylm]-Ks2f[d,3,p]/.γ->Sqrt[1+p^2],{p,0,8}],d>0&&p>0]
FullSimplify[Series[Kllpmf[3,2,0,d,m2Ylm]-Ks2f[d,3,p]/.γ->Sqrt[1+p^2],{p,0,8}],d>0&&p>0]
FullSimplify[Series[Kllpmf[4,2,0,d,p2Ylm]-Ks2f[d,4,p]/.γ->Sqrt[1+p^2],{p,0,8}],d>0&&p>0]
FullSimplify[Series[Kllpmf[4,2,0,d,m2Ylm]-Ks2f[d,4,p]/.γ->Sqrt[1+p^2],{p,0,8}],d>0&&p>0] *)

generate general Doppler operator (s=0) [courtesy JC, Alex Hoey & Jacob Long]

recursion module (general m)

Cfunc[l_,m_,s_]:=If[l==0&&m==s==0,1,If[l>=m&&l>=s,Sqrt[(l^2-m^2)*(l^2-s^2)/(4*l^2-1)]/l,0],0];

KllT[m_,lmax_]:=Module[{l,lp,rT=ConstantArray[0,{lmax+1,lmax+1}],GetK,SetK,Kmmmin,Kp1, Ks},

(* make sure the variable d is not set to some value *)
Clear[d];

(* define functions to access right array elements *)
GetK[l_,lp_]:=If[l>=m&&lp>=m,rT[[l+1,lp+1]],0];
SetK[l_,lp_,val_]:=If[l>=m&&lp>=m,rT[[l+1,lp+1]]=val;,rT[[l+1,lp+1]]=0];

(* initial condition K[m,m] *)
Kmmmin=FullSimplify[FunctionExpand[Kmmm[m,d]/.γ->Sqrt[1+p^2],p>0],p>0]/.(Sqrt[1+(p^2)])->γ;
SetK[m,m,Kmmmin];

(* K[l,lp] will be generated starting by raising K[m,m] --> K[m+1,m] --> ... *) 
For[lp=m,lp<=lmax,lp++,

For[l=Max[lp,m+1],l<=lmax,l++,

If[l>0,
Kp1=-(1/(Cfunc[l,m,0]))*((γ*GetK[l-1,lp]-(GetK[l-1,lp]/.d->d-1))/p)-((Cfunc[l-1,m,0])/(Cfunc[l,m,0]))*GetK[l-2,lp];
SetK[l,lp,Kp1];                                          (* save value *)
];

If[l!=lp,SetK[lp,l,Kp1/.d->2-d/.p->-p]; (* directly set switched kernel value *)];

];
];

Return[rT];
];

setup kernels and access elements all m

KernelTm = Table[KllT[m,lmaxv],{m,0,lmaxv}];

KernelElementmpos[dv_,pv_,l_,lp_,m_]:=KernelTm[[m+1]][[l+1]][[lp+1]]/.γ->Sqrt[1+p^2]/.d->dv/.p->pv (* m>0 *)

KernelElement[dv_,pv_,l_,lp_,m_]:=If[m>=0,KernelElementmpos[dv,pv,l,lp,m],KernelElementmpos[dv,pv,l,lp,-m]];

check explicitly

(* lv=2; lpv=4; mv=-1;
a=KernelElement[d,p,lv,lpv,mv]//Simplify
b=Kllpmf[lv,lpv,mv,d,SphericalHarmonicY]
Series[a-b/.γ->Sqrt[1+p^2],{p,0,10}]//Simplify *)

generate general Doppler operator (s=+/-2 and -2<= m <= 2)

recursion module (general m)

Cfunc[l_,m_,s_]:=If[l==0&&m==s==0,1,If[l>=m&&l>=s,Sqrt[(l^2-m^2)*(l^2-s^2)/(4*l^2-1)]/l,0],0];

KllTs[s_,m_,lmax_]:=Module[{l,lp,rT=ConstantArray[0,{lmax+1,lmax+1}],GetK,SetK,Kmmmin,Kp1, Ks},

(* make sure the variable d is not set to some value *)
Clear[d];

(* define functions to access right array elements *)
GetK[l_,lp_]:=If[l>=m&&lp>=m,rT[[l+1,lp+1]],0];
SetK[l_,lp_,val_]:=If[l>=m&&lp>=m,rT[[l+1,lp+1]]=val;,rT[[l+1,lp+1]]=0];

(* initial condition K[m,m] *)
Kmmmin=Kernels22m[s,m]/.(Sqrt[1+(p^2)])->γ;
SetK[2,2,Kmmmin];

(* K[l,lp] will be generated starting by raising K[m,m] --> K[m+1,m] --> ... *) 
For[lp=2,lp<=lmax,lp++,

For[l=Max[lp,2+1],l<=lmax,l++,

Kp1=-(1/(Cfunc[l,m,s]))*((γ*GetK[l-1,lp]-(GetK[l-1,lp]/.d->d-1))/p)-((s*m/(l*(l-1)))/(Cfunc[l,m,s]))*GetK[l-1,lp]-((Cfunc[l-1,m,s])/(Cfunc[l,m,s]))*GetK[l-2,lp];
SetK[l,lp,Kp1];                                          (* save value *)

If[l!=lp,SetK[lp,l,Kp1/.d->2-d/.p->-p]; (* directly set switched kernel value *)];

];
];

Return[rT];
];

setup kernels and access elements all m

KernelT2m = Table[KllTs[2,m,lmaxv],{m,0,2}];

KernelElementp2mpos[dv_,pv_,l_,lp_,m_]:=KernelT2m[[m+1]][[l+1]][[lp+1]]/.γ->Sqrt[1+p^2]/.d->dv/.p->pv (* m>0 *)

KernelElementp2m[dv_,pv_,l_,lp_,m_]:=If[m>=0,KernelElementp2mpos[dv,pv,l,lp,m],(-1)^(l+lp)*KernelElementp2mpos[dv,-pv,l,lp,-m]];

KernelElement2[dv_,pv_,l_,lp_,s_,m_]:=If[s==2,KernelElementp2m[dv,pv,l,lp,m],(-1)^(l+lp)*KernelElementp2m[dv,-pv,l,lp,m]]

check explicitly

(* 
mv=-1;
Series[KernelElement2[d,p,2,2,2,mv]-Kllpmf[2,2,mv,d,m2Ylm]/.γ->Sqrt[1+p^2],{p,0,10}]//FullSimplify
Series[KernelElement2[d,p,3,2,2,mv]-Kllpmf[3,2,mv,d,m2Ylm]/.γ->Sqrt[1+p^2],{p,0,10}]//FullSimplify
Series[KernelElement2[d,p,2,3,2,mv]-Kllpmf[2,3,mv,d,m2Ylm]/.γ->Sqrt[1+p^2],{p,0,10}]//FullSimplify
Series[KernelElement2[d,p,4,2,2,mv]-Kllpmf[4,2,mv,d,m2Ylm]/.γ->Sqrt[1+p^2],{p,0,10}]//FullSimplify
Series[KernelElement2[d,p,2,4,2,mv]-Kllpmf[2,4,mv,d,m2Ylm]/.γ->Sqrt[1+p^2],{p,0,10}]//FullSimplify

Series[KernelElement2[d,p,2,2,-2,mv]-Kllpmf[2,2,mv,d,p2Ylm]/.γ->Sqrt[1+p^2],{p,0,10}]//FullSimplify
Series[KernelElement2[d,p,3,2,-2,mv]-Kllpmf[3,2,mv,d,p2Ylm]/.γ->Sqrt[1+p^2],{p,0,10}]//FullSimplify
Series[KernelElement2[d,p,2,3,-2,mv]-Kllpmf[2,3,mv,d,p2Ylm]/.γ->Sqrt[1+p^2],{p,0,10}]//FullSimplify
Series[KernelElement2[d,p,4,2,-2,mv]-Kllpmf[4,2,mv,d,p2Ylm]/.γ->Sqrt[1+p^2],{p,0,10}]//FullSimplify
Series[KernelElement2[d,p,2,4,-2,mv]-Kllpmf[2,4,mv,d,p2Ylm]/.γ->Sqrt[1+p^2],{p,0,10}]//FullSimplify
*)

Doppler operators

KernelElement[dv_,pv_,l_,lp_,s_,m_]:=If[s==0,KernelElement[dv,pv,l,lp,m],KernelElement2[dv,pv,l,lp,s,m]]

Doppler[dv_,pv_,l_,li_,lp_]:=((KernelElement[dv+OO,-pv,l,li,0,0]*KernelElement[OO,pv,li,lp,0,0])/(Sqrt[1+pv^2]))
Doppler[dv_,pv_,l_,li_,lp_,m_]:=((KernelElement[dv+OO,-pv,l,li,0,m]*KernelElement[OO,pv,li,lp,0,m])/(Sqrt[1+pv^2]))
DopplerAv[dv_,pv_,l_,li_,lp_]:=Sum[Doppler[dv,pv,l,li,l,m],{m,-Min[l,li,lp],Min[l,li,lp]}]

Doppler[dv_,pv_,l_,li_,lp_,si_,sj_,m_]:=((KernelElement[dv+OO,-pv,l,li,si,m]*KernelElement[OO,pv,li,lp,sj,m])/(Sqrt[1+pv^2]))
DopplerAv[dv_,pv_,l_,li_,lp_,si_,sj_]:=Sum[Doppler[dv,pv,l,li,l,si,sj,m],{m,-Min[l,li,lp],Min[l,li,lp]}]

SZ Operator approach to the thermal & kinematic SZ effect 
Chluba & Rosenberg https://arxiv.org/abs/2508.20659

SZ operators single momenta

single momentum operator (l=0)

D00=Doppler[-1,p,0,0,0];
D02=Doppler[-1,p,0,2,0];
Stot0=D00+(D02/10)-1;

ExpS0[k_]:=Collect[FullSimplify[Normal[Series[Stot0/.β->p/γ/.γ->Sqrt[1+p^2],{p,0,2*k}]]],{p,OO}]

ExpS0[4]

Collect[ExpS0[4]/.ReplOOf[4],{p,DD[_]},FullSimplify]

reorganizing terms I

(((7 DD[2])/5)+((14 DD[3])/25)+((7 DD[4])/150))/((7 )/5)//Expand

(((88 DD[3])/75)+((44 DD[4])/75)+((44 DD[5])/525)+((11 DD[6])/3150))/((88 )/75)//Expand

(((256 DD[4])/315)+((1024 DD[5])/2205)+((64 DD[6])/735)+((128 DD[7])/19845)+((16 DD[8])/99225))/((256/315))//Expand

reorganizing terms II

(-OO+((OO^2)/3))/(1/3)//Expand
(((14 OO)/25)+((7 (OO^2))/30)-((7 (OO^3))/25)+((7 (OO^4))/150))/(1/15)//Expand
(-((44 OO)/105)-((473 (OO^2))/1575)+((209 (OO^3))/1050)+((143 (OO^4))/3150)-((11 (OO^5))/350)+((11 (OO^6))/3150))/(-1/105)//Expand

(((256 OO)/735)+((3392 (OO^2))/11025)-((4736 (OO^3))/33075)-((7856 (OO^4))/99225)+((64 (OO^5))/2205)+((352 (OO^6))/99225)-((64 (OO^7))/33075)+((16 (OO^8))/99225))/(1/945)//Expand

(((4 DD[1])/3)+((DD[2])/3))/(1/3)//Expand

(((7 DD[2])/5)+((14 DD[3])/25)+((7 DD[4])/150))/((1 )/15)//Expand

(((88 DD[3])/75)+((44 DD[4])/75)+((44 DD[5])/525)+((11 DD[6])/3150))/(1/105)//Expand

(((256 DD[4])/315)+((1024 DD[5])/2205)+((64 DD[6])/735)+((128 DD[7])/19845)+((16 DD[8])/99225))/(1/945)//Expand

single momentum operator (l=1)

D10=Doppler[-1,p,1,0,0];
D12=Doppler[-1,p,1,2,0];
Stot1=Sqrt[3]*(D10+(D12/10))+β;

ExpS1[k_]:=Collect[FullSimplify[Normal[Series[Stot1/.β->p/γ/.γ->Sqrt[1+p^2],{p,0,2*k}]]],{p,OO}]

ExpS1[4]

Collect[ExpS1[4]/.ReplOOf[4],{p,DD[_]},FullSimplify]

single momentum operator (l=2)

D20=Doppler[-1,p,2,0,0];
D22=Doppler[-1,p,2,2,0];
Stot2=Sqrt[5]*(D20+(D22/10));

ExpS2[k_]:=Collect[FullSimplify[Normal[Series[Stot2/.β->p/γ/.γ->Sqrt[1+p^2],{p,0,2*k}]]],{p,OO}]

ExpS2[4]

Collect[ExpS2[4]/.ReplOOf[4],{p,DD[_]},FullSimplify]

single momentum operator (l=3)

D30=Doppler[-1,p,3,0,0];
D32=Doppler[-1,p,3,2,0];
Stot3=Sqrt[7]*(D30+(D32/10));

ExpS3[k_]:=Collect[FullSimplify[Normal[Series[Stot3/.β->p/γ/.γ->Sqrt[1+p^2],{p,0,2*k}]]],{p,OO}]

ExpS3[4]

Collect[ExpS3[4]/.ReplOOf[4],{p,DD[_]},FullSimplify]

averaged thermal SZ operator

final result (compare to Table  B1 of first SZpack paper)

OΘ=20;
JJ=ExpandAll[ExpandAll[Simplify[ExpS0[OΘ]]]/.ReplpM[OΘ]/.ReplOOf[OΘ]];
Collect[JJ,{Θ,DD[_]}]

SZ operators averaged operators with kSZ

OΘ=4;
Oβ=3;

ExpfS[kp_,kβ_,S_,l_]:=Collect[Normal[Series[ExpfMBboost[l,kβ]*S/.β->p/γ/.γ->Sqrt[1+p^2]/.p->η*Sqrt[Θ],{Θ,0,kp}]]/.η->p/Sqrt[Θ],{βp,p,OO},Simplify]

final result (l=0)

Irappr=Normal[Series[Collect[ExpfS[OΘ,Oβ,Stot0,0],p]/.ReplpM[OΘ+Oβ],{Θ,0,OΘ}]];
JJ=ExpandAll[ExpandAll[FullSimplify[Irappr]]/.ReplOOf[OΘ+Oβ]];
Collect[JJ,{βp,Θ,DD[_]},FullSimplify]

final result (l=1)

Irappr=Normal[Series[Collect[ExpfS[OΘ,Oβ,Stot1,1],p]/.ReplpM[OΘ+Oβ],{Θ,0,OΘ}]];
JJ=ExpandAll[ExpandAll[FullSimplify[Irappr]]/.ReplOOf[OΘ+Oβ]];
Collect[JJ,{βp,Θ,DD[_]},FullSimplify]

check terms

FullSimplify[(Y+(-10 DD[1]-((47 DD[2])/5)-((7 DD[3])/5))/.DD->Df)-(((-138)/5)*G-(42/5)*Y+(7/5)*x*G*(6*G+x))/.Y->G*(x*Coth[x/2]-4)/.G->x*Exp[x]/(Exp[x]-1)^2]

final result (l=2)

Irappr=Normal[Series[Collect[ExpfS[OΘ,Oβ,Stot2,2],p]/.ReplpM[OΘ+Oβ],{Θ,0,OΘ}]];
JJ=ExpandAll[ExpandAll[FullSimplify[Irappr]]/.ReplOOf[OΘ+Oβ]];
Collect[JJ,{βp,Θ,DD[_]},FullSimplify]

check terms

FullSimplify[(2/3)*G+(((2 DD[1])/3)+((11 DD[2])/30)/.DD->Df)-((11Df[2])/30)/.Y->G*(x*Coth[x/2]-4)/.G->x*Exp[x]/(Exp[x]-1)^2]

FullSimplify[((4 DD[1]+12 DD[2]+6 DD[3]+((19 DD[4])/30))/.DD->Df)-(44*G+12*Y+6*Df[3]+(19/30)*Df[4])/.Y->G*(x*Coth[x/2]-4)/.G->x*Exp[x]/(Exp[x]-1)^2]

final result (l=3)

Irappr=Normal[Series[Collect[ExpfS[OΘ,Oβ,Stot3,3],p]/.ReplpM[OΘ+Oβ],{Θ,0,OΘ}]];
JJ=ExpandAll[ExpandAll[FullSimplify[Irappr]]/.ReplOOf[OΘ+Oβ]];
Collect[JJ,{βp,Θ,DD[_]},FullSimplify]

Comparison to Nozawa et al 1998

Ci

Clear[x]

RepXS={XX->x*Coth[x/2],SS->x/Sinh[x/2]};

C1=10-(47/5)*XX+(7/5)*XX^2+(7/10)*SS^2/.RepXS
C2=25-(1117/10)*XX+(847/10)*XX^2-(183/10)*XX^3+(11/10)*XX^4+SS^2*((847/20)-(183/5)*XX+(121/20)*XX^2)+(11/10)*SS^4/.RepXS

DDf[k_]:=x^k*D[(1/(Exp[x]-1)),{x,k}]
Gx=((x*Exp[x])/((Exp[x]-1)^2));
C1CR=(-10 DD[1]-((47 DD[2])/5)-((7 DD[3])/5))/Gx/.DD->DDf
C2CR=(-25 DD[1]-((1117 DD[2])/10)-((847 DD[3])/10)-((183 DD[4])/10)-((11 DD[5])/10))/Gx/.DD->DDf

FullSimplify[C1-C1CR]
FullSimplify[C2-C2CR]

Di

Clear[x]

RepXS={XX->x*Coth[x/2],SS->x/Sinh[x/2]};

D0=-(2/3)+(11/30)*XX/.RepXS
D1=-4+12*XX-6*XX^2+(19/30)*XX^3+SS^2*(-3+(19/15)*XX)/.RepXS

DDf[k_]:=x^k*D[(1/(Exp[x]-1)),{x,k}]
Gx=((x*Exp[x])/((Exp[x]-1)^2));
D0CR=(((2 DD[1])/3)+((11 DD[2])/30))/Gx/.DD->DDf
D1CR=(4 DD[1]+12 DD[2]+6 DD[3]+((19 DD[4])/30))/Gx/.DD->DDf

FullSimplify[D0-D0CR]
FullSimplify[D1-D1CR]

Polarized SZ effect 
Rosenberg &  Chluba (Sect. 3.1)

DOs2[l_,p_]:=Doppler[-1,p,l,2,0,2,0,0]

l=2 only kinematic (chi factor not included)

Sqrt[((4*Pi)/5)]*m2Ylm[2,0,θ,ϕ]//FullSimplify

ty=FullSimplify[Series[DOs2[2,p]/.p->β*γ/.γ->1/Sqrt[1-β^2],{β,0,6}]]

Collect[ty/(((OO (1+OO))/(3*(Sqrt[5])))),β,ExpandAll]

(-(3/7)-((3 OO)/7)+((OO^2)/7))/(-(1/7))//ExpandAll

(-(11/49)-((2 OO)/21)+((50 (OO^2))/441)-((8 (OO^3))/147)+((4 (OO^4))/441))/(-1)//ExpandAll

Collect[ExpandAll[ty/(((β^2)/(3 (Sqrt[5])))OO*(1+OO))]//.ReplD,β,Simplify]

ExpandAll[ty]/.ReplOOf[porder]//FullSimplify

tyu=Collect[Sqrt[3/8]*(-Sqrt[6]/10*Sqrt[5])(((DD[2] (β^2))/(3 (Sqrt[5])))+(((7 DD[2]+8 DD[3]+DD[4]) (β^4))/(21 (Sqrt[5])))+(((441 DD[2]+1008 DD[3]+486 DD[4]+80 DD[5]+4 DD[6]) (β^6))/(1323 (Sqrt[5])))+((O[β])^7)),β,ExpandAll]

Df[2]/Df[1]//FullSimplify

Collect[tyu/.DD->Df,β,Simplify]

Collect[tyu/(-(1/20) (β^2) ),β,ExpandAll]

l=3 only kinematic (chi factor not included)

Sqrt[((4*Pi)/7)]*m2Ylm[3,0,θ,ϕ]//FullSimplify

ty=FullSimplify[Series[DOs2[3,p]/.p->β*γ/.γ->1/Sqrt[1-β^2],{β,0,6}]]

Collect[ExpandAll[ty/(((OO (1+OO)^2 (β^3))/(3 (Sqrt[35]))))]//.ReplD,β,ExpandAll]

tyD=ExpandAll[ty]/.ReplOOf[8]//FullSimplify

tyu=Collect[Sqrt[15/8]*(-Sqrt[6]/10*Sqrt[7])*tyD/((β^3)/20),β,ExpandAll]

Collect[tyu/.DD->Df,β,Simplify]

l=4 only kinematic (chi factor not included)

Y402=Sqrt[4*Pi/9]*m2Ylm[4,0,θ,ϕ]//FullSimplify

Y402-(Sqrt[(5/18)])*(1+(7/2)*P2)*(1-μ^2)/.P2->LegendreP[2,μ]/.μ->Cos[θ]//FullSimplify

ty=FullSimplify[Series[DOs2[4,p]/.p->β*γ/.γ->1/Sqrt[1-β^2],{β,0,6}]]

Collect[ty/(((OO (((1+OO))^2) (2+OO) (β^4))/(21 (Sqrt[15])))),β,ExpandAll]

tyD=ExpandAll[ty]/.ReplOOf[8]//FullSimplify

tyu=Collect[Sqrt[5/18]*(-Sqrt[6]/10*Sqrt[9])*tyD/(-(β^4)/210),β,ExpandAll]

Collect[tyu/.DD->Df,β,Simplify]

thermal averaged pSZ

definitions of functions

RpMoments=ReplpM[20];
pOrder=Length[RpMoments];

FGen[l_,kβ_]:=ExpfMBboost[l,kβ]

SOf[l_]:=-((Sqrt[6])/10)*Sqrt[2l+1]*DOs2[l,p];

SpSZ[lmax_,kβ_]:=Sum[SOf[l]*FGen[l,kβ]*PG[l]/.γ->Sqrt[1+p^2],{l,2,lmax}]

results

yyu=Normal[Series[SpSZ[lmaxv,2],{βp,0,2},{p,0,pOrder}]]//Simplify;

rrt=Collect[Series[ExpandAll[yyu]/.RpMoments,{Θ,0,10}],{βp,Θ},Expand];

tyD=Collect[ExpandAll[rrt]/.ReplOOf[40],{βp,Θ},ExpandAll];

Collect[(tyD/.PG[2]->(1/2) (Sqrt[(3/2)]) ((βperp^2)/(βp^2)))/(-(1/20) ((βperp)^2)),Θ,Expand]

yyu=Normal[Series[DOs2[2,p],{p,0,6}]]

Collect[FullSimplify[yyu/((((OO+(OO^2)) (p^2))/(3 (Sqrt[5]))))],p,ExpandAll]

yui=Collect[ExpandAll[FullSimplify[yyu]]/.ReplOOf[10],p,ExpandAll]

Collect[yui/(((p^2)/(3 (Sqrt[5])))),p,ExpandAll]

Collect[(((βp)^2) (-((DD[2] PG[2])/(5 (Sqrt[6])))-((Θ (6 DD[2]+8 DD[3]+DD[4]) PG[2])/(5 (Sqrt[6])))-(((Θ^2) (210 DD[2]+1176 DD[3]+867 DD[4]+160 DD[5]+8 DD[6]) PG[2])/(70 (Sqrt[6])))))/(((βp)^2) (-((PG[2])/(5 (Sqrt[6])))))//ExpandAll,Θ,ExpandAll]

Comparison to Itoh et al

first order term

F1JC=(6 DD[2]+8 DD[3]+DD[4])/(((x*Exp[x])/((Exp[x]-1)^2)))/.DD->Df//FullSimplify

F1I=(3/10)*X-(2/5)*(X^2+((S^2)/2))+(1/20)*X*(X^2+2*S^2)/.X->x*Coth[x/2]/.S->x/Sinh[x/2]

(F1JC/20)-F1I//FullSimplify

second order term

F1JC=(15 DD[2]+84 DD[3]+((867 DD[4])/14)+(4/7) (20 DD[5]+DD[6]))/(((x*Exp[x])/((Exp[x]-1)^2)))/.DD->Df//FullSimplify

F1I=(3/4)*X-(21/5)*(X^2+((S^2)/2))+(867/280)*X*(X^2+2*S^2)-(4/7)*(X^4+(11/2)X^2*S^2+S^4)+(1/35)*X*(X^4+13X^2*S^2+(17/2)*S^4)/.X->x*Coth[x/2]/.S->x/Sinh[x/2]

(F1JC/20)-F1I//FullSimplify

plots

functions

YYDef=DD[2]+Θ (6 DD[2]+8 DD[3]+DD[4])+(Θ^2) (15 DD[2]+84 DD[3]+((867 DD[4])/14)+((80 DD[5])/7)+((4 DD[6])/7))+(Θ^3) (((45 DD[2])/4)+399 DD[3]+((48873 DD[4])/56)+((3480 DD[5])/7)+((724 DD[6])/7)+((60 DD[7])/7)+((5 DD[8])/21))+(Θ^4) (-((45 DD[2])/4)+945 DD[3]+((352215 DD[4])/56)+((59100 DD[5])/7)+((12290 DD[6])/3)+((18610 DD[7])/21)+((3865 DD[8])/42)+((40 DD[9])/9)+((5 DD[10])/63))+(Θ^5) (((405 DD[2])/64)+((11655 DD[3])/16)+((23582385 DD[4])/896)+((572925 DD[5])/7)+((1124155 DD[6])/14)+((483297 DD[7])/14)+((1261835 DD[8])/168)+880 DD[9]+((390 DD[10])/7)+((16 DD[11])/9)+((DD[12])/45))+(Θ^6) (((135 DD[2])/16)-((2835 DD[3])/4)+((13458555 DD[4])/224)+((3521475 DD[5])/7)+((1953955 DD[6])/2)+((10687647 DD[7])/14)+((548726887 DD[8])/1848)+((703048 DD[9])/11)+((617725 DD[10])/77)+((59000 DD[11])/99)+((25327 DD[12])/990)+((32 DD[13])/55)+((8 DD[14])/1485))+(Θ^7) (-((22275 DD[2])/512)+((44415 DD[3])/128)+((338307705 DD[4])/7168)+((219802725 DD[5])/112)+((256479605 DD[6])/32)+((2500520425 DD[7])/224)+((71938284587 DD[8])/9856)+((85884082 DD[9])/33)+((507442445 DD[10])/924)+((2366318 DD[11])/33)+((1557899 DD[12])/264)+((16736 DD[13])/55)+((4688 DD[14])/495)+((16 DD[15])/99)+((4 DD[16])/3465))+(Θ^8) (((2025 DD[2])/16)+((2835 DD[3])/4)-((10088955 DD[4])/224)+((122491575 DD[5])/28)+((361098135 DD[6])/8)+((6455349795 DD[7])/56)+((306913646233 DD[8])/2464)+((773681434 DD[9])/11)+((92913071645 DD[10])/4004)+((186233854 DD[11])/39)+((198699571 DD[12])/312)+((40118448 DD[13])/715)+((21062984 DD[14])/6435)+((53336 DD[15])/429)+((44134 DD[16])/15015)+((32 DD[17])/819)+((2 DD[18])/9009))+(Θ^9) (-((5716575 DD[2])/16384)-((13310325 DD[3])/4096)+((4510219725 DD[4])/229376)+((3114514125 DD[5])/896)+((43157406525 DD[6])/256)+((1540227048945 DD[7])/1792)+((122726719610755 DD[8])/78848)+((120688199435 DD[9])/88)+((43790377396475 DD[10])/64064)+((5048356431605 DD[11])/24024)+((16159991401009 DD[12])/384384)+((514951236 DD[13])/91)+((671936774 DD[14])/1287)+((14252930 DD[15])/429)+((8699725 DD[16])/6006)+((34624 DD[17])/819)+((3812 DD[18])/4851)+((16 DD[19])/1911)+((2 DD[20])/51597))+(Θ^10) (((273375 DD[2])/256)+((609525 DD[3])/64)+((190890675 DD[4])/3584)-((91560375 DD[5])/28)+((2958020325 DD[6])/8)+((256517254125 DD[7])/56)+((35887348355535 DD[8])/2464)+((443570077785 DD[9])/22)+((240041457856425 DD[10])/16016)+((13411957255245 DD[11])/2002)+((61816577837225 DD[12])/32032)+((34098012604 DD[13])/91)+((50537987718 DD[14])/1001)+((14429201650 DD[15])/3003)+((1955682835 DD[16])/6006)+((4284256 DD[17])/273)+((33510770 DD[18])/63063)+((70520 DD[19])/5733)+((3175 DD[20])/17199)+((32 DD[21])/19845)+((8 DD[22])/1289925))/.DD->Df;

YYf[xv_,Θv_]:=YYDef/.x->xv/.Θ->Θv;

YYf0[x_,Θ_]:=YYf[x,Θv]/.{((Θv)^10)->0,((Θv)^9)->0,((Θv)^8)->0,((Θv)^7)->0,((Θv)^6)->0,((Θv)^5)->0,((Θv)^4)->0,((Θv)^3)->0,((Θv)^2)->0,Θv->0}
YYf2[x_,Θ_]:=YYf[x,Θv]/.{((Θv)^10)->0,((Θv)^9)->0,((Θv)^8)->0,((Θv)^7)->0,((Θv)^6)->0,((Θv)^5)->0,((Θv)^4)->0,((Θv)^3)->0}/.Θv->Θ

Gx[x_]:=x*Exp[x]/(Exp[x]-1)^2

plots

LogLinearPlot[{x^3*YYf[x,0.03],x^3*YYf0[x,0.03],x^3*YYf2[x,0.03]},{x,0.1,20},WorkingPrecision->100]

LogLinearPlot[{x^3*YYf[x,0.001],x^3*YYf[x,0.01],x^3*YYf[x,0.02],x^3*YYf[x,0.03]},{x,0.1,20},WorkingPrecision->100]

LogLinearPlot[{YYf[x,0.03]/Gx[x],YYf2[x,0.03]/Gx[x],YYf0[x,0.03]/Gx[x]},{x,0.1,20},WorkingPrecision->100]

scattering of CMB anisotropies with polarization 
Rosenberg &  Chluba (Sect. 3.2)

porder=4;
ReplD=ReplacementsOODx[porder];

scaling of p terms

Scalepterms[DD_,porder_]:=Module[{lmax,pArray,pScale,CList=CoefficientList[ExpandAll[DD],p]/.{}->ConstantArray[0,porder+1]},

lmax=Length[DD]-1;
pArray=Table[p^k,{k,0,porder}];
pScale=ConstantArray[1,porder+6];
pScale[[3]]=3; pScale[[5]]=15; pScale[[7]]=105; pScale[[9]]=945;

Return[Table[Sum[ExpandAll[CList[[l+1]][[k+1]]*pScale[[k+1]]]*((pArray[[k+1]])/(pScale[[k+1]])),{k,0,porder}],{l,0,lmax}]];
];

generate Doppler operator Dl0l

Dl0l=Table[Collect[Series[((DopplerAv[-1,p,l,0,l])/(2*l+1))/.γ->Sqrt[1+p^2],{p,0,porder}]//FullSimplify,p,Expand],{l,0,lmaxv}];
Scalepterms[Dl0l,porder]

generate Doppler operator Dl2l

Dl2l=Table[Collect[Series[((DopplerAv[-1,p,l,2,l])/(10*(2*l+1)))/.γ->Sqrt[1+p^2],{p,0,porder}]//Simplify,p,Expand],{l,0,lmaxv}];
Scalepterms[Dl2l,porder]

generate Doppler operator Dl2ls20

Dl2ls20=Table[Collect[Series[((DopplerAv[-1,p,l,2,l,2,0])/(10*(2*l+1)))/.γ->Sqrt[1+p^2],{p,0,porder}]//Simplify,p,Expand],{l,0,lmaxv}];
Scalepterms[Dl2ls20,porder]

compact expression

yyy=Collect[ExpandAll[Dl2ls20]//.ReplD/.ReplpM[5],{Θ,Dx},FullSimplify]

CoefficientList[yyy[[3]]*10*G[x]/.Table[(Dx^n)->(Y[n+1]-3Y[n])/G[x],{n,5,2,-1}]/.Dx->Y[1]/G[x],Θ]//ExpandAll

CoefficientList[yyy[[4]]*G[x]/(-((1 )/(14 (Sqrt[5]))))/.Table[(Dx^n)->(Y[n+1]-3Y[n])/G[x],{n,5,2,-1}]/.Dx->Y[1]/G[x],Θ]//ExpandAll

generate Doppler operator Dl2ls02

Dl2ls02=Table[Collect[Series[((DopplerAv[-1,p,l,2,l,0,2])/(10*(2*l+1)))/.γ->Sqrt[1+p^2],{p,0,porder}]//Simplify,p,Expand],{l,0,lmaxv}];
Scalepterms[Dl2ls02,porder]

Dl2ls20-Dl2ls02

generate Doppler operator Dl2ls22

Dl2ls22=Table[Collect[Series[((DopplerAv[-1,p,l,2,l,2,2])/(10*(2*l+1)))/.γ->Sqrt[1+p^2],{p,0,porder}]//Simplify,p,Expand],{l,0,lmaxv}];
Scalepterms[Dl2ls22,porder]

generate Doppler operator Dl2ls-2-2

Dl2lsm2m2=Table[Collect[Series[((DopplerAv[-1,p,l,2,l,-2,-2])/(10*(2*l+1)))/.γ->Sqrt[1+p^2],{p,0,porder}]//Simplify,p,Expand],{l,0,lmaxv}];
Scalepterms[Dl2lsm2m2,porder]

Dl2ls22-Dl2lsm2m2

(-(6/7)-((11 OO)/14)+((13 (OO^2))/84)+((OO^3)/14)-((OO^4)/84))*-1

generate Doppler operator Dl2ls-2+2

Dl2lsm2p2=Table[Collect[Series[((DopplerAv[-1,p,l,2,l,-2,2])/(10*(2*l+1)))/.γ->Sqrt[1+p^2],{p,0,porder}]//Simplify,p,Expand],{l,0,lmaxv}];
Scalepterms[Dl2lsm2p2,porder]

Scalepterms[Dl2lsm2p2-Dl2lsm2m2,porder]

generate Doppler operator Dl2ls+2-2

Dl2lsp2m2=Table[Collect[Series[((DopplerAv[-1,p,l,2,l,2,-2])/(10*(2*l+1)))/.γ->Sqrt[1+p^2],{p,0,porder}]//Simplify,p,Expand],{l,0,lmaxv}];
Scalepterms[Dl2lsp2m2,porder]

Dl2lsm2p2-Dl2lsp2m2

Evolution equation

simplified using Dx

intensity only

Dl0ld=Collect[Dl0l//.ReplD,p,FullSimplify]

Dl2ld=Collect[Dl2l//.ReplD,p,FullSimplify]

Dlld=Collect[Dl0l+Dl2l//.ReplD,p,FullSimplify]

Dllr={1+((Dx (p^2))/3)-(7/10) (4-Dx) Dx ((p^4)/15)+((11 (10-Dx) (4-Dx) Dx )/30)((p^6)/105),-(2/5)(1+Dx) ((p^2)/3)+ (2/5)(2+((7/2)-Dx) Dx) ((p^4)/15)-(((-10+Dx) (-4+Dx) (2+17 Dx))/70)((p^6)/105),(1/10)-(1/10) (6-Dx) ((p^2)/3)+((144/35)-Dx +((Dx^2)/7)) ((p^4)/15)+(((-10+Dx) (720+Dx (-146+23 Dx)) (p^6))/22050),(3/70)(4-Dx) ((p^2)/3)-(3/70) (12-Dx) (4-Dx) ((p^4)/15)-(((-4+Dx) (3600+Dx (-482+23 Dx)) (p^6))/66150),(((10-Dx) (4-Dx))/105)((p^4)/15)-(((20-Dx) (10-Dx) (4-Dx) (p^6))/11025),(((18-Dx) (10-Dx) (4-Dx) )/693)((p^6)/105)};

Dlld-Dllr//Simplify

polarization

Dl2ls20d=Scalepterms[ExpandAll[Dl2ls20]//.ReplD,porder]

630/15

Collect[Dl2ls20d,p,FullSimplify]

Dl2ls20r={0, 0,(1/10)-(1/10) (4-Dx) ((p^2)/3)+(2/35) (10-Dx) (4-Dx) ((p^4)/15)+(((-18+Dx) (-10+Dx) (-4+Dx) (p^6))/4410),(1/(14(Sqrt[5])))(4-Dx) ((p^2)/3)-(((10-Dx) (4-Dx))/(14 (Sqrt[5]))) ((p^4)/15)-(((-18+Dx) (-10+Dx) (-4+Dx) (p^6))/(2646 (Sqrt[5]))),(((10-Dx) (4-Dx))/(42 (Sqrt[15])))((p^4)/15)+(((-18+Dx) (-10+Dx) (-4+Dx) (p^6))/(4410 (Sqrt[15]))) ,-(((-18+Dx) (-10+Dx) (-4+Dx) (p^6))/(20790 (Sqrt[35])))};

Dl2ls20d-Dl2ls20r//Simplify

Dl2ls22d=Collect[Dl2ls22//.ReplD,p,FullSimplify]

Dl2ls22r={0, 0,(1/10)-(1/30)(10-Dx) ((p^2)/3)+(((18-Dx) (10-Dx))/105) ((p^4)/15)+(((-28+Dx) (-18+Dx) (-10+Dx) (p^6))/44100),(1/42) (4-Dx) ((p^2)/3)-(((18-Dx) (4-Dx))/84)((p^4)/15)-(((-28+Dx) (-18+Dx) (-4+Dx) (p^6))/26460),(((10-Dx) (4-Dx))/252)((p^4)/15)+(((-28+Dx) (-10+Dx) (-4+Dx) (p^6))/44100),-(((-18+Dx) (-10+Dx) (-4+Dx) (p^6))/207900)};

Dl2ls22d-Dl2ls22r//Simplify

DDl2lsm2p2d=Collect[Dl2lsm2p2-Dl2ls22//.ReplD,p,FullSimplify]

DDl2lsm2p2r={0, 0,(2/15) (2+Dx) ((p^2)/3)-(2/15) (5-Dx) (2+Dx) ((p^4)/15)+(((-10+Dx) (2+Dx) (-96+17 Dx) (p^6))/22050),(1/42) (4-Dx) (2+Dx) ((p^4)/15)-(((-12+Dx) (-4+Dx) (2+Dx) (p^6))/4410),(((-10+Dx) (-4+Dx) (2+Dx) (p^6))/33075),0};

DDl2lsm2p2d-DDl2lsm2p2r//Simplify

630/15

DDl2lsm2p2r={0, 0,(1/3)((4/5)+((2Dx)/5)) ((p^2)/3)-(((4 )/3)+((2 )/5)(Dx-(DDx/3))) ((p^4)/15),(2/21)(2+((1 )/2)(Dx-(DDx/2))) ((p^4)/15),0,0};

DDl2lsm2p2d-DDl2lsm2p2r//Simplify

matrices

LDlld=CoefficientList[Dlld/.p^4->15*Θ^2/.p^2->3*Θ*(1+5/2*Θ),Θ]/.{}->ConstantArray[0,3]//ExpandAll

LDl2ls20d=CoefficientList[Dl2ls20d/.p^4->15*Θ^2/.p^2->3*Θ*(1+5/2*Θ),Θ]/.{}->ConstantArray[0,3]//ExpandAll

LDl2ls22d=CoefficientList[Dl2ls22d/.p^4->15*Θ^2/.p^2->3*Θ*(1+5/2*Θ),Θ]/.{}->ConstantArray[0,3]//ExpandAll

LDl2lsm2p2d=CoefficientList[(DDl2lsm2p2d+Dl2ls22d)/.p^4->15*Θ^2/.p^2->3*Θ*(1+5/2*Θ),Θ]/.{}->ConstantArray[0,3]//ExpandAll

lowest order terms (Θ=0)

CreateA[l_,el_]:={{LDlld[[l+1]][[el+1]],-((Sqrt[6])/2)*LDl2ls20d[[l+1]][[el+1]],-((Sqrt[6])/2)*LDl2ls20d[[l+1]][[el+1]]},{-Sqrt[6]*LDl2ls20d[[l+1]][[el+1]],3*LDl2ls22d[[l+1]][[el+1]],3*LDl2lsm2p2d[[l+1]][[el+1]]},{-Sqrt[6]*LDl2ls20d[[l+1]][[el+1]],3*LDl2lsm2p2d[[l+1]][[el+1]],3*LDl2ls22d[[l+1]][[el+1]]}}

CreateA[0,0]//MatrixForm
CreateA[1,0]//MatrixForm
CreateA[2,0]/(1/10)//MatrixForm
CreateA[3,0]//MatrixForm
CreateA[4,0]//MatrixForm

lowest order terms (Θ)

MatrixForm[A0=CreateA[0,1]/.Dx->0//FullSimplify]
MatrixForm[A1=CreateA[1,1]/-(2/5)/.Dx->0//FullSimplify]
MatrixForm[A2=CreateA[2,1]/-(3/5)/.Dx->0//FullSimplify]
MatrixForm[A3=CreateA[3,1]/(6/35)/.Dx->0//FullSimplify]
MatrixForm[A4=CreateA[4,1]/.Dx->0//FullSimplify]

CreateA[0,1]-A0/.DDx->0//FullSimplify//MatrixForm
(CreateA[1,1]/-(2/5)-A1)/.DDx->0//FullSimplify//MatrixForm
(CreateA[2,1]/-(3/5)-A2)*(-(3/5))/(1/10)/.DDx->0//FullSimplify//MatrixForm
(CreateA[3,1]/(6/35)-A3)*(6/35)/-(3/70)/.DDx->0//FullSimplify//MatrixForm
(CreateA[4,1]-A4)/.DDx->0//FullSimplify//MatrixForm

ExpandAll[OO*(OO^2-3*OO)]//.ReplOOf[4]//ExpandAll

MatrixForm[(CreateA[2,0]+Θ*CreateA[2,1]-DiagonalMatrix[{1,1,1}]).{n2m,0,0}-{n2m (-(9/10)-(3/5)Θ+(Dx/10) Θ),n2m (-((Sqrt[6])/10)+(3/5)*(Sqrt[(8/3)])Θ-(((Sqrt[6])Dx)/10)Θ),n2m (-((Sqrt[(3/2)])/5)-(Sqrt[6]) (-(2/5)+(Dx/10)) Θ)}//Simplify]

MatrixForm[(CreateA[2,0]+Θ*CreateA[2,1]-DiagonalMatrix[{1,1,1}]).{0,np2m,nm2m}-{ (-((Sqrt[6])/10)+(3/5)*(Sqrt[(8/3)])Θ-(((Sqrt[6])Dx)/10)Θ)*((np2m+nm2m)/2),(-(7/5)-2*Θ+(Dx/5)Θ)*((np2m+nm2m)/2)+(1+(4/5)*Θ+((2Dx)/5)Θ)*nm2m,(-(7/5)-2*Θ+(Dx/5)Θ)*((np2m+nm2m)/2)+(1+(4/5)*Θ+((2Dx)/5)Θ)*np2m}]//Simplify

FullSimplify[((-(7/5)-2*Θ+(Dx/5)Θ)*((np2m+nm2m)/2)+(1+(4/5)*Θ+((2Dx)/5)Θ)*nm2m+(-(7/5)-2*Θ+(Dx/5)Θ)*((np2m+nm2m)/2)+(1+(4/5)*Θ+((2Dx)/5)Θ)*np2m)/2]/.(nm2m+np2m)->2//Expand

FullSimplify[((-(7/5)-2*Θ+(Dx/5)Θ)*((np2m+nm2m)/2)+(1+(4/5)*Θ+((2Dx)/5)Θ)*nm2m-((-(7/5)-2*Θ+(Dx/5)Θ)*((np2m+nm2m)/2)+(1+(4/5)*Θ+((2Dx)/5)Θ)*np2m))/2]/.(nm2m-np2m)->-2//Expand

(-(2/5)-((6 Θ)/5)+((3 Dx Θ)/5))/-(2/5)//ExpandAll

terms in evolution equation T=0

AM=CreateA[2,0]-DiagonalMatrix[{1,1,1}];
MatrixForm[AM]

AMn=AM.{nI,nQ+I nU,nQ-I nU}//FullSimplify

{{1,0,0},{0,1/2,1/2},{0,1/2/I,-1/2/I}}.AMn//FullSimplify

terms in evolution equation T>0

AM=CreateA[2,0]+Θ*CreateA[2,1]-DiagonalMatrix[{1,1,1}];
MatrixForm[AM]

AMn=AM.{nI,nQ+I nU,nQ-I nU}//FullSimplify

Sn={{1,0,0},{0,1/2,1/2},{0,1/2/I,-1/2/I}}.AMn//FullSimplify

MatrixForm[Collect[Sn,{nI,nQ,nU,Θ},FullSimplify]]

Challinor 2000

first order terms

(-(2/5)(1-(Dx/4)) ((p^2)/3)+((4 )/5)*((20/7)-(Dx-(DDx/14))) ((p^4)/15))*10/.(p^4)->0/.(p^2)->3//ExpandAll

FullSimplify[(ExpandAll[OO*(OO^2-3*OO)]//.ReplOOf[4]/.DD->Df)/Gx-4-(-6*F+(1/2)*(2*F^2+G^2))/.F->x*Coth[x/2]/.G->x/Sinh[x/2]]

second order term

(-(2/5)(1-(Dx/4)) ((p^2)/3)+((4 )/5)*((20/7)-(Dx-(DDx/14))) ((p^4)/15))*10/.(p^4)->15/.(p^2)->3*5/2//ExpandAll

Y1f=ExpandAll[OO*(OO^2-3*OO)]//.ReplOOf[4]//Expand
Y2f=ExpandAll[OO^2*(OO^2-3*OO)]//.ReplOOf[4]//Expand
Y3f=ExpandAll[OO^3*(OO^2-3*OO)]//.ReplOOf[4]//Expand

FullSimplify[(((90/7)*Gx-(11/2)*Y1f-(12/7)*Y2f+(4/7)*Y3f)/Gx/.DD->Df)-(-15*F+(69/4)*(2*F^2+G^2)-(64/7)*F*(F^2+2*G^2)+(4/7)*(F^4+(1/2)*G^2*(11*F^2+2*G^2)))/.F->x*Coth[x/2]/.G->x/Sinh[x/2]]

Paper plots 
Rosenberg &  Chluba Figures 1 & 2

velocity pSZ only

B1=DD[2]/.DD->Df;
B2=DD[2]+(4/3)*DD[3]+(5/21)*DD[4]/.DD->Df;
B3=-(DD[2]+DD[3])/.DD->Df;
B4=2*DD[3]+DD[4]/.DD->Df;

LogLinearPlot[{x^3*B1,x^3*B2/2,x^3*B3/4,x^3*B4/20},{x,0.1,50},WorkingPrecision->20,PlotRange->All]

xmin=0.1;
xmax=50;
np=500;
DDEx=Table[N[{x,x^3*B1,x^3*B2/2,x^3*B3/4,x^3*B4/20},10]/.x->Exp[lx],{lx,Log[xmin],Log[xmax],Log[xmax/xmin]/np}];
DDEx=Append[DDEx,{}];

Export["/Users/user/Desktop/pSZ-Paper/Fig1.dat",DDEx];

v^2 term with temperature corrections

define function

Dn[xv_,Θv_,order_]:=(DD[2]+Θ (6 DD[2]+8 DD[3]+DD[4])+(Θ^2) (15 DD[2]+84 DD[3]+((867 DD[4])/14)+((80 DD[5])/7)+((4 DD[6])/7))+(Θ^3) (((45 DD[2])/4)+399 DD[3]+((48873 DD[4])/56)+((3480 DD[5])/7)+((724 DD[6])/7)+((60 DD[7])/7)+((5 DD[8])/21))+(Θ^4) (-((45 DD[2])/4)+945 DD[3]+((352215 DD[4])/56)+((59100 DD[5])/7)+((12290 DD[6])/3)+((18610 DD[7])/21)+((3865 DD[8])/42)+((40 DD[9])/9)+((5 DD[10])/63))+(Θ^5) (((405 DD[2])/64)+((11655 DD[3])/16)+((23582385 DD[4])/896)+((572925 DD[5])/7)+((1124155 DD[6])/14)+((483297 DD[7])/14)+((1261835 DD[8])/168)+880 DD[9]+((390 DD[10])/7)+((16 DD[11])/9)+((DD[12])/45))+(Θ^6) (((135 DD[2])/16)-((2835 DD[3])/4)+((13458555 DD[4])/224)+((3521475 DD[5])/7)+((1953955 DD[6])/2)+((10687647 DD[7])/14)+((548726887 DD[8])/1848)+((703048 DD[9])/11)+((617725 DD[10])/77)+((59000 DD[11])/99)+((25327 DD[12])/990)+((32 DD[13])/55)+((8 DD[14])/1485))+(Θ^7) (-((22275 DD[2])/512)+((44415 DD[3])/128)+((338307705 DD[4])/7168)+((219802725 DD[5])/112)+((256479605 DD[6])/32)+((2500520425 DD[7])/224)+((71938284587 DD[8])/9856)+((85884082 DD[9])/33)+((507442445 DD[10])/924)+((2366318 DD[11])/33)+((1557899 DD[12])/264)+((16736 DD[13])/55)+((4688 DD[14])/495)+((16 DD[15])/99)+((4 DD[16])/3465))+(Θ^8) (((2025 DD[2])/16)+((2835 DD[3])/4)-((10088955 DD[4])/224)+((122491575 DD[5])/28)+((361098135 DD[6])/8)+((6455349795 DD[7])/56)+((306913646233 DD[8])/2464)+((773681434 DD[9])/11)+((92913071645 DD[10])/4004)+((186233854 DD[11])/39)+((198699571 DD[12])/312)+((40118448 DD[13])/715)+((21062984 DD[14])/6435)+((53336 DD[15])/429)+((44134 DD[16])/15015)+((32 DD[17])/819)+((2 DD[18])/9009))+(Θ^9) (-((5716575 DD[2])/16384)-((13310325 DD[3])/4096)+((4510219725 DD[4])/229376)+((3114514125 DD[5])/896)+((43157406525 DD[6])/256)+((1540227048945 DD[7])/1792)+((122726719610755 DD[8])/78848)+((120688199435 DD[9])/88)+((43790377396475 DD[10])/64064)+((5048356431605 DD[11])/24024)+((16159991401009 DD[12])/384384)+((514951236 DD[13])/91)+((671936774 DD[14])/1287)+((14252930 DD[15])/429)+((8699725 DD[16])/6006)+((34624 DD[17])/819)+((3812 DD[18])/4851)+((16 DD[19])/1911)+((2 DD[20])/51597))+(Θ^10) (((273375 DD[2])/256)+((609525 DD[3])/64)+((190890675 DD[4])/3584)-((91560375 DD[5])/28)+((2958020325 DD[6])/8)+((256517254125 DD[7])/56)+((35887348355535 DD[8])/2464)+((443570077785 DD[9])/22)+((240041457856425 DD[10])/16016)+((13411957255245 DD[11])/2002)+((61816577837225 DD[12])/32032)+((34098012604 DD[13])/91)+((50537987718 DD[14])/1001)+((14429201650 DD[15])/3003)+((1955682835 DD[16])/6006)+((4284256 DD[17])/273)+((33510770 DD[18])/63063)+((70520 DD[19])/5733)+((3175 DD[20])/17199)+((32 DD[21])/19845)+((8 DD[22])/1289925))/.DD->Df/.Table[Θ^k->0,{k,10,order+1,-1}])/.Θ->Θv/.x->xv//Simplify

plots

Θv=2/100;
Dn0=Dn[x,Θv,0];
Dn2=Dn[x,Θv,2];
Dn10=Dn[x,Θv,10];

LogLinearPlot[{x^3*Dn0,x^3*Dn2,x^3*Dn10},{x,0.1,50},WorkingPrecision->20,PlotRange->All]

xmin=1/10;
xmax=50;
np=500;

Θv=2/100;
Dn0=Dn[x,Θv,0];
Dn2=Dn[x,Θv,2];
Dn10=Dn[x,Θv,10];

DDEx=Table[N[{x,x^3*Dn0,x^3*Dn2,x^3*Dn10}/.x->Exp[lx],10],{lx,Log[xmin],Log[xmax],Log[xmax/xmin]/np}];
DDEx=Append[DDEx,{}];

Export["/Users/user/Desktop/pSZ-Paper/Fig2.dat",DDEx];