Linear perturbations slowly rotating, viscous stars

Quit[]

SetDirectory@NotebookDirectory[]
<<EDCRGTCcode.m

SSimplify[expr_]:=Series[expr,{ϵ,0,1}]//Normal//Simplify;
SSolve[eq_,var_]:= {var->(SSimplify[var /. First@Solve[eq,var]])};
LO0[expr_] :=(SeriesCoefficient[ SeriesCoefficient[expr, {δ,0,1}],{ϵ,0,0}]//Normal)/. ϵ->1 /.δ->1;
LO1[expr_] :=(SeriesCoefficient[ SeriesCoefficient[expr, {δ,0,1}],{ϵ,0,1}]//Normal)/. ϵ->1 /.δ->1;;

1. Background

1.1. Construct background

coord={t,r,θ,ϕ};
met = DiagonalMatrix[{-Exp[ν[r]],Exp[λ[r]],(r^2),(r^2)((Sin[θ])^2)}]+((O[ϵ])^2);
met[[1,4]]=met[[4,1]]=-ϵ (r^2)ω[r]((Sin[θ])^2) + ((O[ϵ])^2);
RGtensors[met,coord,{0,0,0}];

1.2. Define velocity field, stress energy tensor and particle current

PD[a_]:= covD[a].uU;

uU={Exp[-ν[r]/2],0,0,Exp[-ν[r ]/2] ϵ Ω0};
ud=gdd.uU;
ud.uU//SSimplify
Jd = n[r]ud;
projdd=gdd+Outer[Times, ud,ud];
projUd=gUU.projdd;
projUU=projUd.gUU;
ϑ=Contract[covD[ud].gUU,{1,2}];
cdud=covD[ud];
sheardd =Table[Sum[ (1/2)(projUd[[a,μ]]projUd[[b,ν]]+projUd[[a,ν]]projUd[[b,μ]]-(2/3)projdd[[μ,ν]]projUU[[a,b]])cdud[[a,b]],{a,1,4},{b,1,4}],{μ,1,4},{ν,1,4}]//SSimplify;
ℰ=ρ[r] + τe[r] (PD[ρ[r]]+(ρ[r]+p[r])ϑ);
𝒫 = p[r]-ζ[r] ϑ+τp[r](PD[ρ[r]]+(ρ[r]+p[r])ϑ);
Qd=τq[r](ρ[r]+p[r]) PD[ud]+βe[r]covD[ρ[r]].projUd+βn[r]covD[n[r]].projUd /. βn[r]->0 /. βe[r]->τq[r] ((cs[r])^2)/.βn[r]->0;
Tdd = ℰ Outer[Times, ud, ud] + 𝒫 projdd + Outer[Times, ud, Qd]+Outer[Times, Qd,ud]-2 η[r] sheardd;

Notice that the expansion and shear both vanish...

ϑ
sheardd//Normal

1.3. TOV + Hartle equation

Construct relevant equations

Contract[covD[Jd ].gUU,{1,2}]
cset=Contract[covD[Tdd].gUU,{2,3}]//SSimplify;
eeqs=Rdd-(1/2)R gdd - 8π Tdd //SSimplify;

Eliminate sound speed

csRule={ρ'[r]->((p'[r])/((cs[r])^2))};

Usual TOV equations

νsol=First@Solve[(eeqs[[2,2]]//Normal)==0, ν'[r]]//Simplify
λsol=First@Solve[(eeqs[[1,1]]//Normal)==0, λ'[r]]//Simplify
psol=First@Solve[(cset[[2]]/.νsol/.csRule//Normal//Simplify)==0,p'[r]]//Simplify

Check -only thing left is equation for ω(r)

eeqs //.csRule//.D[νsol,r]//.νsol//.λsol//.psol//Simplify;
%//MatrixForm//Normal

SeriesCoefficient[eeqs[[1,4]],{ϵ,0,1}]/.νsol/.D[νsol,r]/.psol //. λsol//Simplify;
ωsol=Collect[First@Solve[%==0, ω''[r]]//FullSimplify,{ω'[r],ω[r]}]

1.4. Collect equations -write in terms of m

λ2m= FuncRepRules[λ[r],-Log[1-2m[r]/r]];
(λ'[r]/.FuncRepRules[λ[r],-Log[1-2m[r]/r]])==(λ'[r]/.λsol)/.FuncRepRules[λ[r],-Log[1-2m[r]/r]]//Simplify;
mSOL=First@Solve[%, m'[r]]//Simplify;
νSOL=νsol/.λ2m//Simplify;
pSOL=psol/.λ2m//FullSimplify;
ωSOL=ωsol/.λ2m//Simplify;

Check against Hartle

jHartle[r_] := Exp[-(ν[r]+λ[r])/2];
0==(1/(r^4))D[(r^4)jHartle[r] D[Ω0-ω[r],r],r]+(4/r)jHartle'[r] (Ω0-ω[r])/.λsol//.νsol//.ωsol//Simplify

2. Linear Perturbations

2.0. Utils

vars = {h0[r],h0'[r],h0''[r],h1[r],h1'[r],h1''[r],Z[r],Z'[r],Z''[r],ψ[r],ψ'[r],ψ''[r]};
ν2Rule = D[νsol,r] //.csRule//.νsol//.psol//.λsol;
λ2Rule = D[λsol,r] //.csRule//.νsol//.psol//.λsol;
p2Rule = D[psol,r] //.csRule//.νsol//.psol//.λsol;
cs2Rule=D[csRule,r]/.p2Rule/.psol//Simplify;
bgRules=Flatten@Join[{ν2Rule,λ2Rule,p2Rule,cs2Rule,csRule,νsol,λsol,psol,ωsol}];
ORG[expr_,vrs_:vars] := Collect[expr,vrs,Simplify[#/.bgRules]&];

ZeroCoefficientTable[exprList_,vrs_:vars]:=Module[{bigTab,vrsString},
bigTab=Table[Table[If[TrueQ[!PossibleZeroQ[Coefficient[expr,vrs[[i]]]]],"●"," "],{i,Length[vrs]}],{expr,exprList}];
vrsString=Table[ToString[vr],{vr,vrs}];
bigTab=Prepend[bigTab,vrsString];
Grid[bigTab,Frame->False,Dividers->{False,Table[True,{i,1,Length[exprList]+2}]},Alignment->Center]
]

TorStar[expr_,var_] := expr /. D[var,r] -> Exp[(λ[r]-ν[r])/2]dz[var] /. D[var,{r,2}] -> Exp[λ[r]-ν[r]]*(dz[dz[var]]-dz[var]*(((E^((1/2) (-λ[r]+ν[r]))) (-1+(E^(λ[r]))+4 (E^(λ[r])) π (r^2) p[r]-4 (E^(λ[r])) π (r^2) ρ[r]))/r));
PF ={ η[r]->0,η'[r]->0 , τq[r]->0, τq'[r]->0,τe[r]->0,ζ[r]->0};
NR= {ϵ->0};
dz[0]:=0;

EXTERIOR = { η[r]->0,η'[r]->0 , τq[r]->0, τq'[r]->0,τe[r]->0,ζ[r]->0, ρ[r]->0, p[r]->0, ν[r]-> Log[1-((2 MS)/r)], λ[r]->-Log[1-((2 MS)/r)], ω[r]-> 2(JS/(r^3)), ω'[r]->- 6(JS/(r^4))};
ψToΨ = FuncRepRules[ψ[r], Ψ[r](1-((𝓂 ϵ ω[r] )/σ))];
h1ToMaster = FuncRepRules[h1[r], r Exp[((- ν[r])/2)+ ((λ[r])/2)] ψ[r]];

2.1. Vector Spherical Harmonics

Xd={Xθ[θ,ϕ], Xϕ[θ,ϕ]};
XU = {Xθ[θ,ϕ], (1/((Sin[θ])^2))Xϕ[θ,ϕ]};
XToYW = Flatten@{FuncRepRules[Xθ[θ,ϕ], ((-1)/(Sin[θ]))D[Y[θ,ϕ],ϕ]],FuncRepRules[Xϕ[θ,ϕ], Sin[θ]D[Y[θ,ϕ],θ]], D[Y[θ,ϕ],{θ,2}]->(1/2)(W[θ,ϕ]-ℓ(ℓ+1) Y[θ,ϕ])};
XToY=Flatten@{FuncRepRules[Xθ[θ,ϕ], ((-1)/(Sin[θ]))D[Y[θ,ϕ],ϕ]],FuncRepRules[Xϕ[θ,ϕ], Sin[θ]D[Y[θ,ϕ],θ]]};
WToY = {W[θ,ϕ]-> 2 D[Y[θ,ϕ],{θ,2}]+ℓ(ℓ+1) Y[θ,ϕ]};
YToW={D[Y[θ,ϕ],{θ,2}]->(1/2)(W[θ,ϕ]-ℓ(ℓ+1) Y[θ,ϕ])};
YEQ={
D[Y[θ,ϕ],{θ,2}]-> -Cot[θ]D[Y[θ,ϕ],θ]-((D[Y[θ,ϕ],{ϕ,2}])/((Sin[θ])^2))-ℓ(ℓ+1)Y[θ,ϕ],
D[Y[θ,ϕ],{θ,3}]->ℓ (1+ℓ) Cot[θ] Y[θ,ϕ]+3 Cot[θ] ((Csc[θ])^2) Derivative[0,2][Y][θ,ϕ]-ℓ Derivative[1,0][Y][θ,ϕ]-(ℓ^2) Derivative[1,0][Y][θ,ϕ]+((Cot[θ])^2) Derivative[1,0][Y][θ,ϕ]+((Csc[θ])^2) Derivative[1,0][Y][θ,ϕ]-((Csc[θ])^2) Derivative[1,2][Y][θ,ϕ],
D[Y[θ,ϕ],θ,θ,ϕ]-> -Cot[θ]D[Y[θ,ϕ],θ,ϕ]-((D[Y[θ,ϕ],{ϕ,3}])/((Sin[θ])^2))-ℓ(ℓ+1)D[Y[θ,ϕ],ϕ],
Derivative[0,2][Y][θ,ϕ]->-(1/2) Sin[θ] (Sin[θ] (W[θ,ϕ]+ℓ (1+ℓ)  Y[θ,ϕ])+2 Cos[θ] Derivative[1,0][Y][θ,ϕ])
};
ruleANG2={Derivative[1,1][Y][θ,ϕ]->(1/2) XX[θ,ϕ]+Cot[θ] Derivative[0,1][Y][θ,ϕ],Derivative[0,2][Y][θ,ϕ]->-(1/2) Sin[θ] (Sin[θ] (W[θ,ϕ]+l (1+l) Y[θ,ϕ])+2 Cos[θ] Derivative[1,0][Y][θ,ϕ]),Derivative[2,0][Y][θ,ϕ]->(1/2) (W[θ,ϕ]-l (1+l) Y[θ,ϕ]),Derivative[3,0][Y][θ,ϕ]->(1/2) (-3 Cot[θ] W[θ,ϕ]-l (1+l) Cot[θ] Y[θ,ϕ]-2 (-2+l+(l^2)+((Csc[θ])^2)) Derivative[1,0][Y][θ,ϕ]-2 ((Csc[θ])^2) Derivative[1,2][Y][θ,ϕ]),Derivative[4,0][Y][θ,ϕ]->(1/2) (-(7+2 l (1+l)-11 ((Csc[θ])^2)) W[θ,ϕ]+2 ((Csc[θ])^4) Derivative[0,4][Y][θ,ϕ]+Cot[θ] (5 l (1+l) Cot[θ] Y[θ,ϕ]+2 (-6+5 ((Csc[θ])^2)) Derivative[1,0][Y][θ,ϕ]+12 ((Csc[θ])^2) Derivative[1,2][Y][θ,ϕ])),Derivative[2,1][Y][θ,ϕ]->-(1/2) Cot[θ] XX[θ,ϕ]-(l+(l^2)+((Cot[θ])^2)) Derivative[0,1][Y][θ,ϕ]-((Csc[θ])^2) Derivative[0,3][Y][θ,ϕ],Derivative[2,2][Y][θ,ϕ]->-((Csc[θ])^2) Derivative[0,4][Y][θ,ϕ]+(1/2) l (1+l) Sin[θ] (Sin[θ] (W[θ,ϕ]+l (1+l) Y[θ,ϕ])+2 Cos[θ] Derivative[1,0][Y][θ,ϕ])-Cot[θ] Derivative[1,2][Y][θ,ϕ],Derivative[3,1][Y][θ,ϕ]->(1/4) (-2 (1+l+(l^2)-2 ((Csc[θ])^2)) XX[θ,ϕ]+((Csc[θ])^3) ((7 Cos[θ]+Cos[3 θ]) Derivative[0,1][Y][θ,ϕ]+12 Cos[θ] Derivative[0,3][Y][θ,ϕ]-4 Sin[θ] Derivative[1,3][Y][θ,ϕ]))}//.l->ℓ;
XX[θ,ϕ] = 2 (D[Y[θ,ϕ],θ,ϕ]-Cot[θ]D[Y[θ,ϕ],ϕ]);
YToSphHarm[ll_,mm_] := FuncRepRules[Y[θ,ϕ], SphericalHarmonicY[ll,mm,θ,ϕ]];
$Assumptions={θ∈Reals, ϕ∈Reals, ℓ∈Integers, ℓ!=0, ℓ>=2};
YToSphHarmC[ll_,mm_] := FuncRepRules[Y[θ,ϕ],Simplify[Conjugate[ SphericalHarmonicY[ll,mm,θ,ϕ]]]];

ProjectXθ[expr_]:=Module[
{exprclean0,exprclean,expanded,α,β,η,χ,expr2,expr3,expr4},
exprclean0=SeriesCoefficient[expr,{δ,0,1}];
exprclean=SeriesCoefficient[exprclean0,{ϵ,0,0}]+SeriesCoefficient[exprclean0,{ϵ,0,1}]ϵ+((O[ϵ])^2);
expanded= exprclean//.XToYW //ExpandAll;
α = Coefficient[expanded, Cos[θ]D[Y[θ,ϕ],θ]]//Simplify;
expr2 = (expanded-α Cos[θ]D[Y[θ,ϕ],θ] //Simplify)//Expand;
β=Select[(Coefficient[-Sin[θ]expr2, D[Y[θ,ϕ],ϕ]]//Expand)+θ, FreeQ[#,θ]&]//Simplify;
expr3 = (expr2 +  β ((D[Y[θ,ϕ],ϕ])/(Sin[θ]))//Simplify)//Expand;
η = Coefficient[expr3, Sin[θ] Y[θ,ϕ]]//Simplify;
expr4 =  (expr3 -η Y[θ,ϕ]Sin[θ]//Simplify)//Expand;
χ= Coefficient[expr4, Sin[θ] W[θ,ϕ]] //Simplify;
Return[{α,β,η,χ}]
];
XEquation[{α_,β_,η_,χ_}] := ℓ(ℓ+1)β +I 𝓂 ((ℓ-1)(ℓ+2)χ+α+η);
ProjectII[expr_]:=XEquation[ProjectXθ[expr]];
ProjectXθϕ[expr_] := Module[
{exprclean0,exprclean,expanded,tt,gg,expr2},
exprclean0=SeriesCoefficient[(expr/(Sin[θ])),{δ,0,1}];
exprclean=SeriesCoefficient[exprclean0,{ϵ,0,0}]+SeriesCoefficient[exprclean0,{ϵ,0,1}]ϵ+((O[ϵ])^2);
expanded= exprclean//.XToYW //ExpandAll;
tt =Select[( Coefficient[expanded, W[θ,ϕ]]//Expand)+θ, FreeQ[#,θ]&]//Simplify;
expr2 = (expanded-tt W[θ,ϕ]//Simplify)//Expand;
gg =Select[( Coefficient[expr2,D[ Y[θ,ϕ],ϕ]]//Expand)+θ+1, FreeQ[#,θ]&]-1//Simplify;
Return[{tt,gg}]
];
X2Equation[{tt_,gg_}] := 2(ℓ(ℓ+1)tt +I 𝓂 gg); (* Factor of 2 to match the xAct`Harmonics` convention (Brizuela+). Without the factor of 2, it matches the Gerlach-Sengupta *)
ProjectIII[expr_]:=X2Equation[ProjectXθϕ[expr]];

2.2. Metric and fluid perturbation

hdd = 0*IdentityMatrix[4];
hdd[[1,3]] = hdd[[3,1]] = h0[r]Xd[[1]]Exp[-I σ t];
hdd[[1,4]] = hdd[[4,1]] = h0[r]Xd[[2]]Exp[-I σ t];
hdd[[2,3]] = hdd[[3,2]] = h1[r]Xd[[1]]Exp[-I σ t];
hdd[[2,4]] = hdd[[4,2]] = h1[r]Xd[[2]]Exp[-I σ t];
hUd = gUU.hdd;
hUU=hUd.gUU;
δuU = Exp[-ν[r]/2]Exp[-I σ t]{ϵ Exp[-ν[r]](Ω0 h0[r]+Z[r](Ω0-ω[r])) Xd[[2]],0, ((Z[r])/(r^2))Xd[[1]], Z[r] ((Xd[[2]])/((r^2)((Sin[θ])^2)))};
δud=gdd.δuU//SSimplify;

uU.(hdd.uU)+2(gdd.uU).δuU//SSimplify//Normal

2.3. Brute force it

Pmet = met + δ hdd + ((O[δ])^2)+((O[ϵ])^2);
RGtensors[Pmet, coord,{0,0,0}];

2.4. Einstein and SET tensors

Gdd1=Rdd-(1/2)R gdd;

δuU =Exp[-ν[r]/2]Exp[-I σ t]{ϵ Exp[-ν[r]](Ω0 h0[r]+Z[r](Ω0-ω[r])) Xd[[2]],0, ((Z[r])/(r^2))Xd[[1]], Z[r] ((Xd[[2]])/((r^2)((Sin[θ])^2)))};
PuU = Table[SeriesData[δ,0,{uU[[a]], δuU[[a]]},0,2,1],{a,1,4}];
Pud = gdd.PuU;
PD[a_]:=covD[a].PuU;
projdd=gdd+Outer[Times,Pud,Pud];
projUd=gUU.projdd;
projUU=projUd.gUU;
ϑ1=Contract[covD[Pud].gUU,{1,2}];
cdud=covD[Pud];
sheardd1 =Table[Sum[ (1/2)(projUd[[a,μ]]projUd[[b,ζ]]+projUd[[a,ζ]]projUd[[b,μ]]-(2/3)projdd[[μ,ζ]]projUU[[a,b]])cdud[[a,b]],{a,1,4},{b,1,4}],{μ,1,4},{ζ,1,4}]//SSimplify;
ℰ1=ρ[r] + τe[r] (PD[ρ[r]]+(ρ[r]+p[r])ϑ1);
𝒫1 = p[r]-ζ[r] ϑ1+τp[r](PD[ρ[r]]+(ρ[r]+p[r])ϑ1);
Qd1=τq[r](ρ[r]+p[r]) PD[Pud]+βe[r]covD[ρ[r]].projUd+βn[r]covD[n[r]].projUd /. βn[r]->0 /. βe[r]->τq[r] ((cs[r])^2)/.βn[r]->0;
Tdd1 = ℰ1 Outer[Times, Pud, Pud] + 𝒫1 projdd + Outer[Times, Pud, Qd1]+Outer[Times, Qd1,Pud]-2 η[r] sheardd1;

2.4.5. Check Israel condition

Tdd1[[2,3]]//Simplify;
Series[%,{ϵ,0,1}]//Normal;
%/.δ->1;
%/.p[r]->0 /. ρ[r]->0 /. τq[r]->0 /. η[r]->0

Tdd1[[2,4]]//Simplify;
Series[%,{ϵ,0,1}]//Normal;
%/.δ->1;
%/.p[r]->0 /. ρ[r]->0 /. τq[r]->0 /. η[r]->0

Tdd1[[2,3]]//Simplify;
Series[%,{ϵ,0,1}]//Normal;
%/.δ->1/. τq[r]->0

Tdd1[[2,4]]//Simplify;
Series[%,{ϵ,0,1}]//Normal;
%/.δ->1 /. τq[r]->0

2.5. Einstein eqs & conservation of SET (takes a couple minutes to evaluate)

Useful check: when applying bgRules, the Einstein equations in the background must be satisfied

T0 = (Tdd1//Normal//.δ->0);
Gdd1-8π T0//Normal//.δ->0;
%//.bgRules//.XToY//.δ->0;
%//.λ2m//Simplify//MatrixForm

CSET = Contract[covD[Tdd1].gUU,{2,3}];

e1=ProjectII[Gdd1[[1,3]]-8π Tdd1[[1,3]]//.XToY//.YEQ]/.δ->1;
e2=ProjectII[Gdd1[[2,3]]-8π Tdd1[[2,3]]//.XToY//.YEQ]/.δ->1;
e3=ProjectIII[Gdd1[[3,4]]-8 π Tdd1[[3,4]]//.XToY//.YEQ]/.δ->1;

Be patient with this one!!!

e4 = ProjectII[CSET[[3]]//.XToY//.YEQ]/.δ->1;

3. Checks

3.0. Check equations against non-rotating limit

Match linearized Einstein Equations

OLDTONOW={
x[]->r,
Λ[r]->λ[r],
Φ[r]->ν[r],
𝓅[]->p[r],
ℯ[]->ρ[r],
(ℒ^2)->ℓ(ℓ+1),
η[]->η[r],
kn[]->Exp[- λ[r]/2] h1[r],
kl[]->Exp[-ν[r]/2]h0[r],
β[]->Exp[-ν[r]/2]Z[r],
ψ[]->ψ[r],
τQ[]->τq[r]
};
OLDTONOW2={
PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][a_]]:>D[a,{r,2}],
PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][a_]]:>-(σ^2)a,
PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][a_]]:>-I σ D[a,r],
PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][a_]:> -I σ a,
PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][a_]:>  D[a,r]};

tmp1=ORG[((2 e1/.ϵ->0)/(ℓ(ℓ+1) Exp[-I t σ]))];
(kl[]*(5+2*E^Λ[x[]]*(-3+2*ℒ^2+8*Pi*x[]^2*ℯ[])+E^(2*Λ[x[]])*(1+8*Pi*x[]^2*𝓅[])^2))/(4*E^Λ[x[]]*x[]^2)+16*Pi*(ℯ[]+𝓅[])*β[]+(16*Pi*(ℯ[]+𝓅[])*τQ[]*PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][kl[]])/E^(Φ[x[]]/2)-(E^(-1/2*Λ[x[]]-Φ[x[]]/2)*(-5+E^Λ[x[]]*(1+8*Pi*x[]^2*𝓅[]))*PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][kn[]])/(2*x[])+(16*Pi*(ℯ[]+𝓅[])*τQ[]*PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][β[]])/E^(Φ[x[]]/2)+E^(-1/2*Λ[x[]]-Φ[x[]]/2)*PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][kn[]]]+((-1+E^(-Λ[x[]])+4*Pi*x[]^2*(ℯ[]-𝓅[]))*PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][kl[]])/x[]-PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][kl[]]]/E^Λ[x[]]//.OLDTONOW //.OLDTONOW2;
tmp2=ORG[Exp[ν[r]/2]*%];
ORG[tmp1-tmp2]==0

tmp1=ORG@ORG[((2 e2/.ϵ->0)/(ℓ(ℓ+1) Exp[-I t σ]))];
((-2+ℒ^2)*kn[])/x[]^2+(8*Pi*(-5+E^Λ[x[]]*(1+8*Pi*x[]^2*𝓅[]))*β[]*η[])/(E^(Λ[x[]]/2)*x[])-(E^(-1/2*Λ[x[]]-Φ[x[]]/2)*(-5+E^Λ[x[]]*(1+8*Pi*x[]^2*𝓅[]))*PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][kl[]])/(2*x[])+(16*Pi*η[]*PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][kn[]])/E^(Φ[x[]]/2)+PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][kn[]]]/E^Φ[x[]]-E^(-1/2*Λ[x[]]-Φ[x[]]/2)*PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][kl[]]]+(16*Pi*η[]*PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][β[]])/E^(Λ[x[]]/2)//.OLDTONOW //.OLDTONOW2;
tmp2=ORG[Exp[λ[r]/2]*%];
ORG[tmp1-tmp2]==0

ORG[(( e3/.ϵ->0)/(ℓ(ℓ+1) Exp[-I t σ]))];
tmp1=ORG[((%)/(Coefficient[%, h1'[r]]))];
E^(Φ[x[]]/2)*kn[]*(-1+E^Λ[x[]]*(1+8*Pi*x[]^2*𝓅[]))+32*E^(Λ[x[]]/2+Φ[x[]]/2)*Pi*x[]*β[]*η[]-2*E^(Λ[x[]]/2)*x[]*PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][kl[]]+2*E^(Φ[x[]]/2)*x[]*PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][kn[]] //.OLDTONOW//.OLDTONOW2;
tmp2=ORG[((%)/(Coefficient[%, h1'[r]]))];
ORG[tmp1-tmp2]==0

Match conservation of SET

ORG[(( e4/.ϵ->0)/(ℓ(ℓ+1) Exp[-I t σ]))];
tmp1=ORG[((%)/(Coefficient[%, Z''[r]]))];
-(((ℯ[]+𝓅[])*PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][kl[]])/E^(Φ[x[]]/2))-((ℯ[]+𝓅[])*PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][β[]])/E^(Φ[x[]]/2)-((ℯ[]+𝓅[])*τQ[]*PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][kl[]]])/E^Φ[x[]]-((ℯ[]+𝓅[])*τQ[]*PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][β[]]])/E^Φ[x[]]+(E^(-1/2*Λ[x[]]-Φ[x[]]/2)*PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][kn[]]*(2*η[]+x[]*PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][η[]]))/x[]+(PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][β[]]*(-((3+E^Λ[x[]]*(-3+8*Pi*x[]^2*(ℯ[]-2*𝓅[])))*η[])+2*x[]*PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][η[]]))/(2*E^Λ[x[]]*x[])+(β[]*(-2*ℒ^2*η[]+8*Pi*x[]^2*(3*ℯ[]+𝓅[])*η[]+(x[]-(5*x[])/E^Λ[x[]]+8*Pi*x[]^3*𝓅[])*PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][η[]]))/(2*x[]^2)+E^(-1/2*Λ[x[]]-Φ[x[]]/2)*η[]*PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][PDxPertGSGM`Private`SS[{0,-xPertGSGM`Private`SS}][kn[]]]+(η[]*PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][PDxPertGSGM`Private`SS[{1,-xPertGSGM`Private`SS}][β[]]])/E^Λ[x[]]//.OLDTONOW //.OLDTONOW2;
tmp2=ORG[((%)/(Coefficient[%, Z''[r]]))];
ORG[tmp1-tmp2]==0

3.1. Check: no viscosity & no rotation

EqsPF=ORG[{e1,e2,e3,e4}//.PF//.NR];
ZeroCoefficientTable[EqsPF]

h0PF = First@Solve[(EqsPF[[3]]==0)//FullSimplify, h0[r]]//Simplify;
ORG[EqsPF[[2]]//.FuncRepRules[h0[r], h0[r]/.h0PF]];
%//.h1ToMaster//.bgRules;
TorStar[%,ψ[r]]//FullSimplify;
MASTER$PF=ORG[((%)/(Coefficient[%,dz[dz[ψ[r]]]]))/.λ2m,{ψ[r],dz[ψ[r]],dz[dz[ψ[r]]]}]

It matches!

PreviousEq = dz[dz[ψ[r]]] + ((σ^2) -((Exp[ν[r]])/(r^2))(ℓ(ℓ+1)-6((m[r])/r)+4π (r^2)(-p[r]+ρ[r])))ψ[r];
MASTER$PF== PreviousEq//Simplify

3.2. Check: vacuum exterior

Eqs = ORG[{e1,e2,e3,e4}//.bgRules//.EXTERIOR];
ZeroCoefficientTable[Eqs]

CoefficientList[Eqs[[2]], h0'[r]];
h0psol = {h0'[r]-> -%[[1]]/%[[2]]}//Simplify;
Eqs[[1]]//.D[h0psol,r]//.h0psol;
ORG[%];
CoefficientList[%, h0[r]]//Simplify;
h0sol=FuncRepRules[h0[r],-%[[1]]/%[[2]]//Simplify];
RW1=ORG@Simplify@Normal@Series[Eqs[[2]] /. h0sol,{ϵ,0,1}];
RW1/.h1ToMaster//.bgRules//.EXTERIOR//Simplify;
RW2=ORG[((%)/(Coefficient[%, ψ''[r]]))]/.ψToΨ//.FuncRepRules[ ω[r], 2(JS/(r^3))];
MASTER=ORG[Series[((%)/(Coefficient[%, Ψ''[r]])),{ϵ,0,1}]//Normal//Simplify];

Check against Pani

Pani = (1/((F[r])^2))(((F[r])^2)φ''[r]+F'[r]F[r]φ'[r]+((σ^2)-((2 ϵa m ζ[r] σ)/(r^2))-F[r](((l(l+1))/(r^2))-((6MS)/(r^3))+((24 JS m (-7 MS+3 r) ϵa)/(l (1+l) (r^6) σ))))φ[r] ) /. FuncRepRules[φ[r], Ψ[r]] /. FuncRepRules[F[r], 1-2MS/r]/.FuncRepRules[ζ[r], 2JS/r]/.m->𝓂/.l->ℓ/.ϵa->ϵ;
MASTER==Pani//Simplify

3.3. Perfect fluid (rotating)

Eqs = ORG[{e1,e2,e3,e4}//.bgRules//.PF];
ZeroCoefficientTable[Eqs]
ZeroCoefficientTable[Eqs//.NR]

Write in terms of abstract coefficients (careful, explicit out that certain terms are 𝒪(ϵ)!)

e11 = a1[r]h0[r]+a2[r]h0'[r]+a3[r]h0''[r]+a4[r]h1[r]+a5[r]h1'[r]+a6[r]Z[r];
e22 = b1[r]h0[r]+b2[r]h0'[r]+b3[r]h1[r];
e33 = c1[r]h0[r]+ϵ c21[r]h0'[r]+c3[r]h1[r]+c4[r]h1'[r]+ϵ c51[r]Z[r];
e44=d1[r]h0[r]+d2[r]Z[r];

Map coefficients!

Clear[a1,a2,a3,a4,a5,a6,b1,b2,b3,c1,c2,c3,c4,c5,d1,d2];

ACOEFSSERIES0={
a1[r]-> a10[r]+ ϵ a11[r],
a2[r]-> a20[r]+ϵ a21[r],
a3[r] -> a30[r]+ϵ a31[r],
a4[r]->a40[r]+ϵ a41[r],
a5[r]->a50[r]+ϵ a51[r],
a6[r]->a60[r]+ϵ a61[r],
b1[r]->b10[r]+ϵ b11[r],
b2[r]->b20[r]+ϵ b21[r],
b3[r]->b30[r]+ϵ b31[r],
c1[r]->c10[r]+ϵ c11[r],
c3[r]->c30[r]+ϵ c31[r],
c4[r]->c40[r]+ϵ c41[r],
d1[r]->d10[r]+ϵ d11[r],
d2[r]->d20[r]+ϵ d21[r]
};
ACOEFSERIES=Flatten@Join[{ACOEFSSERIES0, D[ACOEFSSERIES0,r], D[ACOEFSSERIES0,{r,2}]}];

abstract$coefs = {
a10[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h0[r]],{ϵ,0,0}],
a11[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h0[r]],{ϵ,0,1}],
a20[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h0'[r]],{ϵ,0,0}],
a21[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h0'[r]],{ϵ,0,1}],
a30[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h0''[r]],{ϵ,0,0}],
a31[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h0''[r]],{ϵ,0,1}],
a40[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h1[r]],{ϵ,0,0}],
a41[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h1[r]],{ϵ,0,1}],
a50[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h1'[r]],{ϵ,0,0}],
a51[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h1'[r]],{ϵ,0,1}],
a60[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], Z[r]],{ϵ,0,0}],
a61[r] -> SeriesCoefficient[Coefficient[Eqs[[1]],Z[r]],{ϵ,0,1}],
b10[r] -> SeriesCoefficient[Coefficient[Eqs[[2]], h0[r]],{ϵ,0,0}],
b11[r] -> SeriesCoefficient[Coefficient[Eqs[[2]], h0[r]],{ϵ,0,1}],
b20[r] -> SeriesCoefficient[Coefficient[Eqs[[2]], h0'[r]],{ϵ,0,0}],
b21[r] -> SeriesCoefficient[Coefficient[Eqs[[2]], h0'[r]],{ϵ,0,1}],
b30[r] -> SeriesCoefficient[Coefficient[Eqs[[2]], h1[r]],{ϵ,0,0}],
b31[r] -> SeriesCoefficient[Coefficient[Eqs[[2]], h1[r]],{ϵ,0,1}],
c10[r] -> SeriesCoefficient[Coefficient[Eqs[[3]], h0[r]],{ϵ,0,0}],
c11[r] -> SeriesCoefficient[Coefficient[Eqs[[3]], h0[r]],{ϵ,0,1}],
c21[r] -> SeriesCoefficient[Coefficient[Eqs[[3]], h0'[r]],{ϵ,0,1}],
c30[r] -> SeriesCoefficient[Coefficient[Eqs[[3]], h1[r]],{ϵ,0,0}],
c31[r] -> SeriesCoefficient[Coefficient[Eqs[[3]], h1[r]],{ϵ,0,1}],
c40[r] -> SeriesCoefficient[Coefficient[Eqs[[3]], h1'[r]],{ϵ,0,0}],
c41[r] -> SeriesCoefficient[Coefficient[Eqs[[3]], h1'[r]],{ϵ,0,1}],
c51[r] -> SeriesCoefficient[Coefficient[Eqs[[3]], Z[r]],{ϵ,0,1}],
d10[r] -> SeriesCoefficient[Coefficient[Eqs[[4]], h0[r]],{ϵ,0,0}],
d11[r] -> SeriesCoefficient[Coefficient[Eqs[[4]], h0[r]],{ϵ,0,1}],
d20[r] -> SeriesCoefficient[Coefficient[Eqs[[4]],Z[r]],{ϵ,0,0}],
d21[r] -> SeriesCoefficient[Coefficient[Eqs[[4]], Z[r]],{ϵ,0,1}]
};
dabstract$coefs = D[abstract$coefs,r]//.bgRules//Simplify;
d2abstract$coefs=D[dabstract$coefs,r]//.bgRules//Simplify;

Eliminate Z & h0

Zsol=SSolve[e44==0, Z[r]];
h0psol=SSolve[e22==0, h0'[r]];
e11  /.h0psol /. D[h0psol,r]/.h0psol/.Zsol//Simplify;
h0sol=SSolve[%==0, h0[r]]//Simplify;
h0psol=D[h0sol,r]//Simplify;

Put together in equation for h1

e22 /.Zsol//.h0sol/.h0psol //Simplify;
RW0=ORG[%];

RW1=RW0/.h1ToMaster//.bgRules;
RW2=ORG[(RW1/(Coefficient[RW1, ψ''[r]]))]/.ψToΨ//.bgRules;

RW2/.ACOEFSERIES;
tmp=Series[%,{ϵ,0,1}]//Normal;
RW3 =Series[ (tmp/(Coefficient[tmp, Ψ''[r]])),{ϵ,0,1}]//Normal;

Evaluate on the actual coefficients

MASTER=ORG[RW3 //. abstract$coefs //. dabstract$coefs //. d2abstract$coefs];

POT=(-Exp[(ν[r]-λ[r])]Coefficient[MASTER, Ψ[r]]+(σ^2)-2 𝓂 ϵ ω[r]σ  )//FullSimplify
POT0 = POT/.NR//Simplify

POT-POT0//Simplify;
r ℓ(ℓ+1)((σ )/(2 𝓂 ϵ Exp[ν[r]]))%//Simplify

4. Decouple Equations

4.1.  Get Equations

Eqs = ORG[{e1,e2,e3,e4}//.bgRules//Simplify];
ZeroCoefficientTable[Eqs]

4.2. Write in terms of abstract coefficients

e11 = a1[r]h0[r]+a2[r]h0'[r]+a3[r]h0''[r]+a4[r]h1[r]+a5[r]h1'[r]+a6[r]Z[r];
e22 = b1[r]h0[r]+b2[r]h0'[r]+b3[r]h1[r]+b4[r] Z[r]+b5[r]Z'[r];
e33 = c1[r]h0[r]+ϵ c21[r]h0'[r]+c3[r]h1[r]+c4[r]h1'[r]+ϵ c51[r]Z[r];
e44=d1[r]h0[r]+d2[r]h1[r]+d3[r]h1'[r]+d4[r]Z[r]+d5[r]Z'[r]+d6[r]Z''[r];

ACOEFSSERIES0={
a1[r]-> a10[r]+ ϵ a11[r],
a2[r]-> a20[r]+ϵ a21[r],
a3[r] -> a30[r]+ϵ a31[r],
a4[r]->a40[r]+ϵ a41[r],
a5[r]->a50[r]+ϵ a51[r],
a6[r]->a60[r]+ϵ a61[r],
b1[r]->b10[r]+ϵ b11[r],
b2[r]->b20[r]+ϵ b21[r],
b3[r]->b30[r]+ϵ b31[r],
b4[r]->b40[r]+ϵ b41[r],
b5[r]->b50[r]+ϵ b51[r],
c1[r]->c10[r]+ϵ c11[r],
c3[r]->c30[r]+ϵ c31[r],
c4[r]->c40[r]+ϵ c41[r],
d1[r]->d10[r]+ϵ d11[r],
d2[r]->d20[r]+ϵ d21[r],
d3[r]->d30[r]+ϵ d31[r],
d4[r]->d40[r]+ϵ d41[r],
d5[r]->d50[r]+ϵ d51[r],
d6[r]->d60[r]+ϵ d61[r]
};
ACOEFSERIES=Flatten@Join[{ACOEFSSERIES0, D[ACOEFSSERIES0,r], D[ACOEFSSERIES0,{r,2}]}];

abstract$coefs = {
a10[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h0[r]],{ϵ,0,0}],
a11[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h0[r]],{ϵ,0,1}],
a20[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h0'[r]],{ϵ,0,0}],
a21[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h0'[r]],{ϵ,0,1}],
a30[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h0''[r]],{ϵ,0,0}],
a31[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h0''[r]],{ϵ,0,1}],
a40[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h1[r]],{ϵ,0,0}],
a41[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h1[r]],{ϵ,0,1}],
a50[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h1'[r]],{ϵ,0,0}],
a51[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], h1'[r]],{ϵ,0,1}],
a60[r] -> SeriesCoefficient[Coefficient[Eqs[[1]], Z[r]],{ϵ,0,0}],
a61[r] -> SeriesCoefficient[Coefficient[Eqs[[1]],Z[r]],{ϵ,0,1}],
b10[r] -> SeriesCoefficient[Coefficient[Eqs[[2]], h0[r]],{ϵ,0,0}],
b11[r] -> SeriesCoefficient[Coefficient[Eqs[[2]], h0[r]],{ϵ,0,1}],
b20[r] -> SeriesCoefficient[Coefficient[Eqs[[2]], h0'[r]],{ϵ,0,0}],
b21[r] -> SeriesCoefficient[Coefficient[Eqs[[2]], h0'[r]],{ϵ,0,1}],
b30[r] -> SeriesCoefficient[Coefficient[Eqs[[2]], h1[r]],{ϵ,0,0}],
b31[r] -> SeriesCoefficient[Coefficient[Eqs[[2]], h1[r]],{ϵ,0,1}],
b40[r] -> SeriesCoefficient[Coefficient[Eqs[[2]], Z[r]],{ϵ,0,0}],
b41[r] -> SeriesCoefficient[Coefficient[Eqs[[2]], Z[r]],{ϵ,0,1}],
b50[r] -> SeriesCoefficient[Coefficient[Eqs[[2]], Z'[r]],{ϵ,0,0}],
b51[r] -> SeriesCoefficient[Coefficient[Eqs[[2]], Z'[r]],{ϵ,0,1}],
c10[r] -> SeriesCoefficient[Coefficient[Eqs[[3]], h0[r]],{ϵ,0,0}],
c11[r] -> SeriesCoefficient[Coefficient[Eqs[[3]], h0[r]],{ϵ,0,1}],
c21[r] -> SeriesCoefficient[Coefficient[Eqs[[3]], h0'[r]],{ϵ,0,1}],
c30[r] -> SeriesCoefficient[Coefficient[Eqs[[3]], h1[r]],{ϵ,0,0}],
c31[r] -> SeriesCoefficient[Coefficient[Eqs[[3]], h1[r]],{ϵ,0,1}],
c40[r] -> SeriesCoefficient[Coefficient[Eqs[[3]], h1'[r]],{ϵ,0,0}],
c41[r] -> SeriesCoefficient[Coefficient[Eqs[[3]], h1'[r]],{ϵ,0,1}],
c51[r] -> SeriesCoefficient[Coefficient[Eqs[[3]], Z[r]],{ϵ,0,1}],
d10[r] -> SeriesCoefficient[Coefficient[Eqs[[4]], h0[r]],{ϵ,0,0}],
d11[r] -> SeriesCoefficient[Coefficient[Eqs[[4]], h0[r]],{ϵ,0,1}],
d20[r] -> SeriesCoefficient[Coefficient[Eqs[[4]],h1[r]],{ϵ,0,0}],
d21[r] -> SeriesCoefficient[Coefficient[Eqs[[4]], h1[r]],{ϵ,0,1}],
d30[r] -> SeriesCoefficient[Coefficient[Eqs[[4]],h1'[r]],{ϵ,0,0}],
d31[r] -> SeriesCoefficient[Coefficient[Eqs[[4]], h1'[r]],{ϵ,0,1}],
d40[r] -> SeriesCoefficient[Coefficient[Eqs[[4]],Z[r]],{ϵ,0,0}],
d41[r] -> SeriesCoefficient[Coefficient[Eqs[[4]], Z[r]],{ϵ,0,1}],
d50[r] -> SeriesCoefficient[Coefficient[Eqs[[4]],Z'[r]],{ϵ,0,0}],
d51[r] -> SeriesCoefficient[Coefficient[Eqs[[4]], Z'[r]],{ϵ,0,1}],
d60[r] -> SeriesCoefficient[Coefficient[Eqs[[4]],Z''[r]],{ϵ,0,0}],
d61[r] -> SeriesCoefficient[Coefficient[Eqs[[4]], Z''[r]],{ϵ,0,1}]
};
dabstract$coefs = D[abstract$coefs,r]//.bgRules//Simplify;
d2abstract$coefs=D[dabstract$coefs,r]//.bgRules//Simplify;

4.3. Eliminate h0

Zppsol=SSolve[e44==0, Z''[r]];
h0psol=SSolve[e22==0, h0'[r]];
e11  /.h0psol /. D[h0psol,r]/.h0psol/.Zppsol//Simplify;
h0sol=SSolve[%==0, h0[r]]//Simplify;
h0psol=D[h0sol,r]/.Zppsol//Simplify;

4.4. Equation for Z

e44 //.h0sol/.h0psol //Simplify;
FLUID0=ORG[%];
FLUID1=(FLUID0/(Coefficient[FLUID0, Z''[r]]));

FLUID2=  FLUID1 /. h1ToMaster /. ψToΨ //. bgRules;

tmp=FLUID2/.ACOEFSERIES;
Series[tmp,{ϵ,0,1}]//Normal;
FLUID3=ORG[((%)/(Coefficient[%, Z''[r]]))];

FLUID = ORG[FLUID3 /.abstract$coefs/.dabstract$coefs /. d2abstract$coefs];

Put[FLUID, "fluideq.mx"];

4.5. Equation for Ψ

e22 //.h0sol/.h0psol /.Zppsol/.h0sol//Simplify;
RW0=ORG[%];

RW1=RW0/.h1ToMaster//.bgRules;
RW2=ORG[(RW1/(Coefficient[RW1, ψ''[r]]))]/.ψToΨ//.bgRules;

RW2/.ACOEFSERIES;
tmp=Series[%,{ϵ,0,1}]//Normal;
RW3 =Series[ (tmp/(Coefficient[tmp, Ψ''[r]])),{ϵ,0,1}]//Normal;

This evaluation takes long, be patient...

MASTER=ORG[RW3 //. abstract$coefs //. dabstract$coefs //. d2abstract$coefs];

Put[MASTER, "gweq.mx"];

Coefficient[MASTER, Ψ[r]];
Coefficient[%//Expand, (σ^2)]//Simplify

4.6. Non-trivial check!

η[r] * FLUID//Expand;
%//.PF//Simplify;
CoefficientList[%, Z[r]];
ZToΨPF = {Z[r]->Series[((-%[[1]])/(%[[2]])),{ϵ,0,1}]//Normal//Simplify};

MASTER//.PF//.ZToΨPF/.D[ZToΨPF,r]//.bgRules;
ORG[%,{Ψ[r],Ψ'[r],Ψ''[r]}];
%/Coefficient[%, Ψ''[r]]//SSimplify;
ORG[%,{Ψ[r],Ψ'[r],Ψ''[r]}]

4.7. Transform to Tortoise coordinate

tmp=TorStar[TorStar[FLUID, Z[r]],Ψ[r]];
FLUIDT= ORG[(tmp/(Coefficient[tmp,dz[dz[Z[r]]]])),{Z[r],dz[Z[r]],dz[dz[Z[r]]],Ψ[r],dz[Ψ[r]],dz[dz[Ψ[r]]]}];

tmp=TorStar[TorStar[MASTER, Z[r]],Ψ[r]];
MASTERT= ORG[(tmp/(Coefficient[tmp,dz[dz[Ψ[r]]]])),{Z[r],dz[Z[r]],dz[dz[Z[r]]],Ψ[r],dz[Ψ[r]],dz[dz[Ψ[r]]]}];

4.8. Propagation Speeds

Coefficient[FLUIDT, Z[r]];
Coefficient[%//Expand, (σ^2)]//Simplify;
cZ2=Series[1/%,{ϵ,0,1}]//Normal;
((SeriesCoefficient[%,{ϵ,0,1}])/((η[r])/((p[r]+ρ[r]) τq[r])))//FullSimplify;
Collect[%,{Ω0-ω[r]}]
%+((16 I (E^(-((ν[r])/2))) π (r^2) 𝓂 (p[r]+ρ[r])( τq[r]- ((η[r])/(p[r]+ρ[r]))) (Ω0-ω[r]))/(ℓ (1+ℓ)))//FullSimplify;
% * (((16 I (E^(-((ν[r])/2))) π (r^2) (𝓂/(ℓ (1+ℓ)))))^(-1))

Coefficient[MASTERT,Ψ[r]]//Simplify;
cψ2=1/Coefficient[%//Expand, (σ^2)]//Simplify;
Series[%,{ϵ,0,1}]//Normal

4.9. To LaTeX

Unprotect[TeXForm];
TeXForm[τq[r]]:="\\tau_\\Q";
Protected[TeXForm];

ToTeX[expr_] := TeXForm[% /. λ[r]->λ /. ν[r]->ν /. τq[r]->(Subscript[τ, qq])/. τq'[r]-> (Subscript[τ, qq])' /. η'[r]-> η' /. η[r]-> η/. 𝓂->m /. ϵ->1 /. cs[r]-> (Subscript[c, s]) /. ω[r]->ϖ[r] /. σ->ω /. Ω0->Ω /. ρ[r]->ρ /. p[r]->p];

Coefficient[MASTERT, dz[Ψ[r]]]//Simplify;
FullSimplify[((%)/((1/(ℓ (1+ℓ) σ ((cs[r])^2)))8 (E^(-((λ[r])/2))) π r 𝓂 ϵ ))];
ToTeX[%]