Hartle-Thorne formalism for slowly-rotating geometries and perturbations
Reference: Richard Brito, Vitor Cardoso and Paolo Pani, “Black-hole superradiance: from fundamental physics to astrophysics” (2015)

This notebook comprises of 3 parts:
1a) We first derive the equations governing a spinning, self-gravitating perfect fluid in GR to second order in the spin
1b) We compute the Klein-Gordon equation for a probe scalar field propagating on this fixed geometry
2) We separate the angular and radial parts of the KG equation to second order in the spin, in order to obtain a master radial equation describing scalar perturbations on this geometry
3a) The background equations are integrated
3b) The final radial scalar equation is integrated to find the proper modes of oscillations.

Quit[];

Slow-rotation [Hartle-Thorne formalism]:

Background equations

ClearAll[GM0,x,g,ϕ,ϕ0,Λ,Ψ,α,γ,Y,Σ,Δ,NN,ω,Δ,Σ,gff,gtt,ϕ,M,F,V,f0,f1,f2,f3,f4,Ω,A,ϕ,B,ζ,ϕ,Φ1,ω,X,r,u]
$Assumptions={r≥0,θ≥0,θ≤π,J≥0,NN≥0,rh≥0,X≥0};
x={t,r,θ,φ};

ω[r_]:=Ω-ω1[r];

P2[θ_]:=((3(1-X[θ])-1)/2);
K0[r_]:=0;(*gauge choice*)


gdd0={{-Exp[ν[r]](1+((ϵa)^2)2(h0[r]+h2[r] P2[θ]))+((ϵa)^2) ((ω[r])^2)(r^2)X[θ],0,0,-ϵa ω[r] X[θ] (r^2)},{0,(((1-2((M[r])/r)))^(-1))(1+2((ϵa)^2)(((m0[r]+m2[r] P2[θ]))/(r-2 M[r]))),0,0},{0,0,(r^2)(1+2((ϵa)^2)(v2[r]-h2[r])P2[θ]),0},{-ϵa ω[r] X[θ] (r^2),0,0,(r^2)(1+2((ϵa)^2)(v2[r]-h2[r])P2[θ])X[θ]}};

ORDERa=2;
ORDERp=1;

(*$HistoryLength=0;*)

Rule 0th order:

(***)
(*BISERIES[X_]:=(Normal[Series[X/.{ϵa->eps ϵa,ϵp->eps ϵp},{eps,0,1}]]/.eps->1);*)
rule01orders={};
BISERIES[Z_]:=Z/.Union[{((X'[θ])^2)->4 X[θ](1-X[θ]),X''[θ]->2 (1-2X[θ])}];

Metric and co.

gdd=Simplify[BISERIES[Normal[Series[gdd0,{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]]];
guu=Simplify[BISERIES[Normal[Series[Inverse[gdd],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]]];
DET=Simplify[BISERIES[Normal[Series[Det[gdd],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]]];
dim=Length[gdd];
Do[g[i,j]=Simplify[gdd[[i]][[j]]];,{i,1,dim},{j,1,dim}];
MatrixForm[gdd]
MatrixForm[guu]

Basic definitions (it can take some minutes!)

Do[Christ[i,j,k]=Simplify[BISERIES[Normal[Series[Sum[(1/2)guu[[k]][[l]](D[gdd[[l]][[i]],x[[j]]]+D[g[l,j],x[[i]]]-D[g[j,i],x[[l]]]),{l,1,dim}],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]]];,{i,1,dim},{j,1,dim},{k,1,dim}];

Other things (only needed for grav pert)

Timing[Do[Ruddd[a,i,j,k]=Simplify[BISERIES[Normal[Series[D[Christ[i,k,a],x[[j]]]-D[Christ[i,j,a],x[[k]]]+Sum[Christ[i,k,β]Christ[β,j,a]-Christ[i,j,β]Christ[β,k,a],{β,1,dim}],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]]];,{a,1,dim},{i,1,dim},{j,1,dim},{k,1,dim}]]

Timing[Do[Ricci[i,k]=Simplify[BISERIES[Normal[Series[Sum[Ruddd[l,i,l,k],{l,1,dim}],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]]];,{i,1,dim},{k,1,dim}];]
Do[Ruu[i,k]=Simplify[BISERIES[Normal[Series[Sum[Ricci[l,m]guu[[i]][[l]]guu[[k]][[m]],{l,1,dim},{m,1,dim}],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]]];,{i,1,dim},{k,1,dim}];
SCAL=Simplify[BISERIES[Normal[Series[Sum[guu[[i]][[j]]Ricci[i,j],{i,1,dim},{j,1,dim}],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]]];
Do[Eins[i,j]=Simplify[BISERIES[Normal[Series[Ricci[i,j]-(1/2)SCAL g[i,j],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]]];,{i,1,dim},{j,1,dim}];

Klein-Gordon equation

ϕ=ϵp Y[θ,φ]ϕp[r,t];
eqS=Simplify[BISERIES[Normal[Series[(1/(Sqrt[-DET]))Sum[D[(Sqrt[-DET])D[ϕ,x[[i]]]guu[[i]][[j]],x[[j]]],{i,1,dim},{j,1,dim}],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]]];

eqS

SetDirectory[NotebookDirectory[]];
Put[eqS,"KG.m"];

Einstein equations with perfect fluid:

u={(1/(Sqrt[-(g[1,1]+2 ϵa Ω g[1,4]+((ϵa)^2)(Ω^2)g[4,4])])),0,0,((ϵa Ω)/(Sqrt[-(g[1,1]+2 ϵa Ω g[1,4]+((ϵa)^2)(Ω^2)g[4,4])]))};
ΔP=((ϵa)^2)(ρ[r]+P[r])(p0[r]+p2[r] P2[θ]);
Δρ=((ϵa)^2)(ρ[r]+P[r])((ρ'[r])/(P'[r]))(p0[r]+p2[r] P2[θ]);
Do[Tuufluid[i,j]=Simplify[Normal[Series[(P[r]+ρ[r]+ΔP+Δρ)u[[i]]u[[j]]+guu[[i]][[j]] (P[r]+ΔP),{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]];,{i,1,dim},{j,1,dim}];

(*Here, P and ρ are assumed in the Jordan frame ρJ=ρE (A^4)*)
Do[Tfluid[i,j]=Simplify[Normal[Series[(Sum[Tuufluid[a,b]g[a,i]g[b,j],{a,1,dim},{b,1,dim}]),{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]];,{i,1,dim},{j,1,dim}];

MatrixForm[Table[Tfluid[i,j],{i,1,dim},{j,1,dim}]]

Reducing to this particular case

f0[x_]:=(1/(16 π));(*κ=8πG*)
γ[x_]:=0;
Α[x_]:=1;
V[x_]:=0;

Some definitions

Clear[ϕ,Φ1];
ϕ=0;
(*Do[A[i]=Simplify[Normal[Series[Ad[[i]],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]],{i,1,dim}];
Do[FF[i,j]=Simplify[Normal[Series[D[A[j],x[[i]]]-D[A[i],x[[j]]],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]],{i,1,dim},{j,1,dim}];*)
(****)
Do[T[i,j]=Simplify[Normal[Series[(((Α[ϕ])^4)Tfluid[i,j]+γ[ϕ](D[ϕ,x[[i]]]D[ϕ,x[[j]]]+D[ϕ,x[[j]]]D[ϕ,x[[i]]]-g[i,j]Sum[guu[[a]][[b]]D[ϕ,x[[a]]]D[ϕ,x[[b]]],{a,1,dim},{b,1,dim}])-g[i,j] V[ϕ]+2((D[D[f0[ϕ],x[[i]]],x[[j]]]-Sum[Christ[i,j,l]D[f0[ϕ],x[[l]]],{l,1,dim}])-g[i,j]Sum[guu[[a]][[b]] (D[D[f0[ϕ],x[[a]]],x[[b]]]-Sum[Christ[a,b,l]D[f0[ϕ],x[[l]]],{l,1,dim}]),{a,1,dim},{b,1,dim}])),{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]],{i,1,dim},{j,1,dim}]

Einstein equations:

(*Do[EQEINST[i,j]=Simplify[Ricci[i,j]-(1/2)g[i,j] SCAL+(((f0[ϕ]))^(-1))(Η[i,j]+Ι[i,j]+J[i,j]+Κ[i,j]-(1/2) T[i,j])];,{i,1,dim},{j,1,dim}];*)
Do[EQEINST[i,j]=Simplify[BISERIES[Normal[Series[Eins[i,j]+(((f0[ϕ]))^(-1))(-(1/2) T[i,j]),{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]]];,{i,1,dim},{j,1,dim}];

Do[DivE[i]=Simplify[BISERIES[Normal[Series[Sum[guu[[a]][[j]]guu[[b]][[i]](D[EQEINST[a,b],x[[j]]]-Sum[Christ[a,j,l]EQEINST[l,b]+Christ[b,j,l]EQEINST[l,a],{l,1,dim}]),{a,1,dim},{b,1,dim},{j,1,dim}],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]]];,{i,1,dim}]

eq11=EQEINST[1,1];
eq12=EQEINST[1,2];
eq13=EQEINST[1,3];
eq14=EQEINST[1,4];
eq22=EQEINST[2,2];
eq23=EQEINST[2,3];
eq24=EQEINST[2,4];
eq33=EQEINST[3,3];
eq34=EQEINST[3,4];
eq44=EQEINST[4,4];
div1=DivE[1];
div2=DivE[2];
div3=DivE[3];
div4=DivE[4];

(*FullSimplify[DivE[2]]*)

SetDirectory[NotebookDirectory[]];
Put[{eq11,eq12,eq13,eq14,eq22,eq23,eq24,eq33,eq34,eq44,div1,div2,div3,div4},"eqs_HT_2nd_order.m"];

Field equations

Quit[];

Clear[r,X];
SetDirectory[NotebookDirectory[]];
EQS=Get["eqs_HT_2nd_order.m"];
{eq11,eq12,eq13,eq14,eq22,eq23,eq24,eq33,eq34,eq44,div1,div2,div3,div4}=EQS;

X[θ_]:=((Sin[θ])^2);
SS[YY_,ORD_]:=Normal[Series[YY,{ϵa,0,ORD}]];

EQs={eq11,eq12,eq13,eq14,eq22,eq23,eq24,eq33,eq34,eq44,div1,div2,div3,div4};
(*j[r_]:=Exp[-ν[r]/2](Sqrt[1-2((M[r])/r)]);*)
(*h0[r_]:=0;
K0[r_]:=0;
m0[r_]:=0;
(*****)
rulerho={ρ[r]->ρ[R],P[r]->P[R](*,ρ'[r]->ρ'[R],P'[r]->P'[R],ρ''[r]->ρ''[R],P''[r]->P''[R]*)};
(*ruleNOzero={h0[r]->0,K0[r]->0,m0[r]->0};*)
ruleNOzero={};
ruleNOzero=Union[ruleNOzero,D[ruleNOzero,r],D[ruleNOzero,{r,2}]];
*)

syst0=Solve[{SS[eq11,1]==0,SS[eq22,1]==0,SS[div2,1]==0},{M'[r],ν'[r],P'[r]}][[1]]//FullSimplify
syst0=Union[syst0,D[syst0,r],D[syst0,{r,2}]];
syst1=Solve[(SS[eq14,1]//.syst0)==0,ω1''[r]][[1]]//FullSimplify
syst1=Union[syst1,D[syst1,r],D[syst1,{r,2}]];
syst01=Union[syst0,syst1];

systl2=Solve[{(SS[div3,2]//.syst01)==0,(SS[(eq44/((Sin[θ])^2))-eq33,2]//.syst01)==0,(SS[eq23,2]//.syst01)==0},{p2[r],m2[r],v2'[r]}][[1]]//FullSimplify
systl2=Union[systl2,D[systl2,r],D[systl2,{r,2}]];

systl0=Solve[{(SS[div2,2]//.Union[syst01,systl2])==0},{p0'[r]}][[1]]//FullSimplify
systl0=Union[systl0,D[systl0,r],D[systl0,{r,2}]];

AAA=Simplify[SS[eq22,2]//.Union[syst01,systl2,systl0]];
C2=Coefficient[TrigExpand[AAA],Cos[θ],2]//Simplify;
S2=Coefficient[TrigExpand[AAA],Sin[θ],2]//Simplify;
FullSimplify[AAA-((Cos[θ])^2)C2-((Sin[θ])^2)S2]

h2p=Solve[S2==0,h2'[r]][[1]]//FullSimplify
RULE=Union[syst01,systl2,systl0,h2p,D[h2p,r],D[h2p,{r,2}]];
h0p=FullSimplify[Solve[{(AAA//.RULE)==0,(eq11//.RULE)==0},{m0'[r],h0'[r]}][[1]]]
RULE=Union[RULE,h0p,D[h0p,r],D[h0p,{r,2}]];

Put[RULE,"background_HT.m"];

Check equations:

Simplify[{eq11,eq12,eq13,eq14,eq22,eq23,eq24,eq33,eq34,eq44,div1,div2,div3,div4}//.RULE]

E01=(M'[r])-(M'[r]//.RULE);
E02=(ν'[r])-(ν'[r]//.RULE);
E03=(P'[r])-(P'[r]//.RULE);
E11=(ω1''[r])-(ω1''[r]//.RULE);
E21=(p0'[r])-(p0'[r]//.RULE);
E22=(m0'[r])-(m0'[r]//.RULE);
E23=(h0'[r])-(h0'[r]//.RULE);
E24=(p2[r])-(p2[r]//.RULE);
E25=(m2[r])-(m2[r]//.RULE);
E26=(v2'[r])-(v2'[r]//.RULE);
E27=(h2'[r])-(h2'[r]//.RULE);

ruleh={v2->v2h,h2->h2h};
c1h2=Coefficient[E26,h2[r]];
c1v2=Coefficient[E26,v2[r]];
S1=FullSimplify[E26-h2[r]c1h2-v2[r] c1v2-v2'[r]];
E26h=FullSimplify[E26-S1/.ruleh]
(***)
c2h2=Coefficient[E27,h2[r]];
c2v2=Coefficient[E27,v2[r]];
S2=FullSimplify[E27-h2[r]c2h2-v2[r] c2v2-h2'[r]];
E27h=FullSimplify[(E27-S2/.ruleh)]

Put[{E01,E02,E03,E11,E21,E22,E23,E24,E25,E26,E27,E26h,E27h},"field_eqs_GR.m"];

CHECKS

CHECK original Hartle-Thorne equations: OK

Zeroth and first order:

eq3a=P'[r]-(-(ρ[r]+P[r])((M[r]+4 π (r^3)P[r])/(r(r-2M[r]))));
eq3b=M'[r]-4 π (r^2)ρ[r];
eq3d=ν'[r]-(((-2)/(ρ[r]+P[r]))P'[r]);
(*****)
j[r_]:=Exp[((-ν[r])/2)](Sqrt[1-2((M[r])/r)]);
eq9=(1/(r^4))D[(r^4)j[r] ω1'[r],r]+(4/r)j'[r] ω1[r];
Simplify[{eq3a,eq3b,eq3d,eq9}//.RULE]

Second order:

eq15a=m0'[r]-(4 π (r^2)((ρ'[r])/(P'[r]))(ρ[r]+P[r])p0[r]+(1/12)((j[r])^2)(r^4)((ω1'[r])^2)-((r^3)/3)D[((j[r])^2),r]((ω1[r])^2));
eq15b=p0'[r]-(-m0[r](((1+8π (r^2)P[r]))/(((r-2M[r]))^2))-((4π(ρ[r]+P[r])(r^2))/(r-2M[r]))p0[r]+(1/12)(((r^4)((j[r])^2))/(r-2M[r]))((ω1'[r])^2)+(1/3)D[(((r^3)((j[r])^2)((ω1[r])^2))/(r-2M[r])),r]);
eq21a=v2'[r]-(-ν'[r]h2[r]+((1/r)+((ν'[r])/2))(((-(r^3))/3)D[((j[r])^2),r]((ω1[r])^2)+((((j[r])^2)(r^4))/6)((ω1'[r])^2)));
eq21b=h2'[r]-((-ν'[r]+(r/(r-2M[r]))((ν'[r])^(-1))(8 π(ρ[r]+P[r])-((4M[r])/(r^3))))h2[r]-((4 v2[r])/(r(r-2M[r])))((ν'[r])^(-1))+(1/6)(((ν'[r])/2)r-(((ν'[r])^(-1))/(r-2M[r])))(r^3)((j[r])^2)((ω1'[r])^2)-(1/3)(((ν'[r])/2)r+(((ν'[r])^(-1))/(r-2M[r])))(r^2)D[((j[r])^2),r]((ω1[r])^2));
eq23a=m2[r]-((r-2M[r])(-h2[r]-((r^3)/3)D[((j[r])^2),r]((ω1[r])^2)+(r^4)(((j[r])^2)/6)((ω1'[r])^2)));
eq23b=p2[r]-(-h2[r]-((r^2)/3)Exp[-ν[r]]((ω1[r])^2));
Simplify[{eq15a,eq15b,eq21a,eq21b,eq23a,eq23b}//.RULE]

h0'[r]//.RULE

CHECK KERR (note that we are not using Boyer-Lindquist coordinates!)

Quit[];

SetDirectory[NotebookDirectory[]];
EQSGR=Get["field_eqs_GR.m"];
{E01,E02,E03,E11,E21,E22,E23,E24,E25,E26,E27,E26h,E27h}=EQSGR;

ruleKERR={ρ->Function[r,0],P->Function[r,0],M->Function[r,MASS],ν->Function[r,Log[1-2MASS/r]],ω1->Function[r,Ω-2(J/(r^3))],m0->Function[{r},-((J^2)/(r^3))],h0->Function[{r},((J^2)/((r^3) (-2 MASS+r)))],m2->Function[{r},((-12 (J^2) MASS+6 (J^2) r+2 MASS (r^4) h2[r]-(r^5) h2[r])/(r^4))],h2->Function[r,(((J^2) (MASS+r))/(MASS (r^4)))],v2->Function[r,-((J^2)/(r^4))]};Simplify[EQSGR//.ruleKERR]

ERGO=-(E^(ν[r]))+((ϵa)^2) (-2 (E^(ν[r])) h0[r]+(E^(ν[r])) h2[r] (-2+3 ((Sin[θ])^2))+(r^2) ((Sin[θ])^2) ((ω1[r])^2))//FullSimplify
(*ERGO=(E^(ν[r]))(1+2((ϵa)^2)h0[r])-((ϵa)^2)((r^2)  ((ω1[r])^2)+Exp[ν[r]]h2[r])((Sin[θ])^2);*)

Deriving the scalar field equation in final form:

Quit[];

SetDirectory[NotebookDirectory[]];
eqS=Get["KG.m"];
RULE=Get["background_HT.m"];

rule0=RULE;
ϕp[r_,t_]:=((ψ[l][r])/(Sqrt[(r^2)+((ϵa)^2)f[r]])) Exp[-I ω t];
SJ[X_]:=Normal[Series[X,{ϵa,0,2}]];

Angular part:

$Assumptions={r≥0,θ≥0,θ≤π};
Sθ[θ_,φ_]:=((-1)/(Sin[θ]))D[Y[θ,φ],φ];
Sφ[θ_,φ_]:=Sin[θ]D[Y[θ,φ],θ];
X[θ_]:=((Sin[θ])^2);
ruleANG={Derivative[2,0][Y][θ,φ]->-l(l+1)Y[θ,φ]-(1/((Sin[θ])^2))Derivative[0,2][Y][θ,φ]-Cot[θ]Derivative[1,0][Y][θ,φ],Derivative[3,0][Y][θ,φ]->l (1+l) Cot[θ] Y[θ,φ]-(1+l+(l^2)) Derivative[1,0][Y][θ,φ]+((Csc[θ])^2) (3 Cot[θ] Derivative[0,2][Y][θ,φ]+2 Derivative[1,0][Y][θ,φ]-Derivative[1,2][Y][θ,φ]),Derivative[4,0][Y][θ,φ]->l (1+l) (1+l+(l^2)-3 ((Csc[θ])^2)) Y[θ,φ]+Csc[θ] (Csc[θ] (7+l+(l^2)-11 ((Csc[θ])^2)) Derivative[0,2][Y][θ,φ]+((1+2 l (1+l)) Cos[θ]-6 Cot[θ] Csc[θ]) Derivative[1,0][Y][θ,φ]+Csc[θ] (5 Cot[θ] Derivative[1,2][Y][θ,φ]-Derivative[2,2][Y][θ,φ])),Derivative[2,1][Y][θ,φ]->-l (1+l) Derivative[0,1][Y][θ,φ]-((Csc[θ])^2) Derivative[0,3][Y][θ,φ]-Cot[θ] Derivative[1,1][Y][θ,φ],Derivative[2,2][Y][θ,φ]->-l (1+l) Derivative[0,2][Y][θ,φ]-((Csc[θ])^2) Derivative[0,4][Y][θ,φ]-Cot[θ] Derivative[1,2][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][θ,φ]))};

ruleNOφ={Derivative[0,1][Y][θ,φ]->I m Y[θ,φ],Derivative[0,2][Y][θ,φ]->-(m^2) Y[θ,φ],Derivative[1,1][Y][θ,φ]->I m Derivative[1,0][Y][θ,φ],Derivative[1,2][Y][θ,φ]->-(m^2)Derivative[1,0][Y][θ,φ],Derivative[0,3][Y][θ,φ]->-I (m^3) Y[θ,φ]};
ruleT=Union[ruleANG,ruleNOφ,{(Sqrt[((Sin[θ])^2)])->Sin[θ]}];

Decomposition:

EQ=Simplify[SJ[(((((E^(-ν[r]))/(2 (r^4)))(E^(-I t ω)) ϵp ))^(-1))eqS//.ruleT]]

EQ=TrigExpand[EQ];
cY=Expand[Coefficient[EQ,Y[θ,φ]],Trig->True];
cYc2θ=Simplify[Coefficient[cY,Cos[θ],2]];
cYs2θ=Simplify[Coefficient[cY,Sin[θ],2]];
cYcθ=Simplify[Coefficient[cY,Cos[θ],1]];
cY0=Simplify[cY-(cYcθ Cos[θ]+cYc2θ ((Cos[θ])^2)+cYs2θ ((Sin[θ])^2))];
Simplify[cY-((cY0+cYs2θ)+cYcθ Cos[θ]+(cYc2θ-cYs2θ) ((Cos[θ])^2))]
(****)
(***)
cYθ=Simplify[Coefficient[EQ,Derivative[1,0][Y][θ,φ]]]//TrigExpand;
cYθsθ=Simplify[Coefficient[cYθ,Sin[θ]]];
cYθsθcθ=Simplify[Coefficient[cYθsθ,Cos[θ]]];
cYθsθ0=Simplify[cYθsθ-cYθsθcθ Cos[θ]];
(****)
Print["Checks:"];
(*CHECK=Simplify[EQ-((cY0-cYc2θ) Y[θ,φ]+cYc Cos[θ]Y[θ,φ]+2cYc2θ ((Cos[θ])^2)Y[θ,φ]+ (cYθsθ0 Sin[θ]+cYθsθcθ Sin[θ]Cos[θ]) Derivative[1,0][Y][θ,φ])]*)
CHECK=Simplify[EQ-((cY0+cYs2θ) Y[θ,φ]+cYcθ Cos[θ]Y[θ,φ]+(cYc2θ-cYs2θ) ((Cos[θ])^2)Y[θ,φ]+ (cYθsθ0 Sin[θ]+cYθsθcθ Sin[θ]Cos[θ]) Derivative[1,0][Y][θ,φ])]
vect={cY0,cYs2θ,cYcθ,cYc2θ,cYθsθ0,cYθsθcθ};
D[vect,θ]==Table[0,{i,1,Length[vect]}]
D[vect,φ]==Table[0,{i,1,Length[vect]}]

A0=Simplify[cY0+cYs2θ];
At0=cYcθ;
B0=cYθsθ0;
D0=2cYc2θ;
Bt0=cYθsθcθ;

{At0,B0,D0}

Simplify[EQ-((A0+At0 Cos[θ]+D0 ((Cos[θ])^2))Y[θ,φ]+(B0+Bt0 Cos[θ])Sin[θ]Derivative[1,0][Y][θ,φ])]

Massaging

(*ruleTORTOISE={u1[l_]'[r]->((u1[l]'[r])/(f[r])),u1[l_]''[r]->(1/((f[r])^2))(u1[l]''[r]-u1[l]'[r]D[f[r],r])};*)
(******)
rulel={(*u1->u1l,u2->u2l,u3->u3l,u4->u4l,*)l->l};
rulelp1={(*u1->u1lp1,u2->u2lp1,u3->u3lp1,u4->u4lp1,*)l->l+1};
rulelm1={(*u1->u1lm1,u2->u2lm1,u3->u3lm1,u4->u4lm1,*)l->l-1};
rulelp2={(*u1->u1lp2,u2->u2lp2,u3->u3lp2,u4->u4lp2,*)l->l+2};
rulelm2={(*u1->u1lm2,u2->u2lm2,u3->u3lm2,u4->u4lm2,*)l->l-2};

eq0l=(A0+ ((Q[l+1])^2)(D0+l Bt0)+((Q[l])^2)(D0-(l+1)Bt0))/.rulel;
eq0lp1=(Q[l+1](At0+l B0))/.rulelm1;
eq0lm1=(Q[l](At0-(l+1)B0))/.rulelp1;
eq0lp2=(Q[l+1]Q[l+2](D0+l Bt0))/.rulelm2;
eq0lm2=(Q[l]Q[l-1](D0-(l+1)Bt0))/.rulelp2;

EQfD=Simplify[eq0l+eq0lm1+eq0lp1+eq0lp2+eq0lm2]

Solve l+2, l-2:

EQfDlp2=Solve[(EQfD/.l->l+2)==0,(ψ[l+2])''[r]][[1]];
EQfDlm2=Solve[(EQfD/.l->l-2)==0,(ψ[l-2])''[r]][[1]];

EQfDF=Simplify[Normal[Simplify[Series[(EQfD/.EQfDlp2/.EQfDlm2),{ϵa,0,2}]]]];

RULE=Union[RULE,{f'[r]->(1/(6 r (((r-2 M[r]))^2) (M[r]+4 π (r^3) P[r])))(E^(-ν[r])) (12 (E^(ν[r])) f[r] (((r-2 M[r]))^2) (M[r]+4 π (r^3) P[r])+(r^2) (-12 (E^(ν[r])) h2[r] (r-2 M[r]) (((M[r]+4 π (r^3) P[r]))^2) (-1+3 ((Q[l])^2)+3 ((Q[1+l])^2))+r (-20 (r^2) ((M[r])^4) (-1+3 ((Q[l])^2)+3 ((Q[1+l])^2)) (((ω1)'[r])^2)-4 (r^3) ((M[r])^3) (-8 π r ((ω1[r])^2) ρ'[r]+24 π r ((Q[l])^2) ((ω1[r])^2) ρ'[r]+24 π r ((Q[1+l])^2) ((ω1[r])^2) ρ'[r]+9 (((ω1)'[r])^2)-24 ((Q[l])^2) (((ω1)'[r])^2)-24 ((Q[1+l])^2) (((ω1)'[r])^2)+4 π (-1+3 ((Q[l])^2)+3 ((Q[1+l])^2)) ρ[r] ω1[r] (9 ω1[r]+8 r (ω1)'[r])+4 π P[r] (-1+3 ((Q[l])^2)+3 ((Q[1+l])^2)) (9 ((ω1[r])^2)+8 r ω1[r] (ω1)'[r]+4 (r^2) (((ω1)'[r])^2)))+(r^2) ((M[r])^2) (-48 (E^(ν[r])) π p0[r] (ρ[r]-2 r ρ'[r])+64 (π^2) (r^4) ((P[r])^2) (-1+3 ((Q[l])^2)+3 ((Q[1+l])^2)) (-14 ((ω1[r])^2)-8 r ω1[r] (ω1)'[r]+(r^2) (((ω1)'[r])^2))+(r^2) (32 π r (-1+3 ((Q[l])^2)+3 ((Q[1+l])^2)) ((ω1[r])^2) ρ'[r]-3 (-7+17 ((Q[l])^2)+17 ((Q[1+l])^2)) (((ω1)'[r])^2)+8 π ρ[r] ω1[r] ((-9+39 ((Q[l])^2)+39 ((Q[1+l])^2)) ω1[r]+16 r (-1+3 ((Q[l])^2)+3 ((Q[1+l])^2)) (ω1)'[r]))-8 π P[r] (6 (E^(ν[r])) p0[r]+(r^2) (((ω1[r])^2) (9-112 π (r^2) ρ[r]-16 π (r^3) ρ'[r]+((Q[l])^2) (-39+336 π (r^2) ρ[r]+48 π (r^3) ρ'[r])+((Q[1+l])^2) (-39+336 π (r^2) ρ[r]+48 π (r^3) ρ'[r]))+16 r (-1+3 ((Q[l])^2)+3 ((Q[1+l])^2)) (-1+4 π (r^2) ρ[r]) ω1[r] (ω1)'[r]-2 (r^2) (-8+21 ((Q[l])^2)+21 ((Q[1+l])^2)) (((ω1)'[r])^2))))+4 π (r^3) (12 (E^(ν[r])) m0[r] P[r] (1+4 π (r^2) P[r]-4 π (r^2) ρ[r])+(r^2) (160 (π^2) (r^5) ((P[r])^3) (-1+3 ((Q[l])^2)+3 ((Q[1+l])^2)) ((ω1[r])^2)+6 (E^(ν[r])) p0[r] ρ'[r]+4 π r ((P[r])^2) (6 (E^(ν[r])) p0[r]+(r^2) (3 ((Q[l])^2) ((-4+40 π (r^2) ρ[r]) ((ω1[r])^2)-8 r ω1[r] (ω1)'[r]+(r^2) (((ω1)'[r])^2))+3 ((Q[1+l])^2) ((-4+40 π (r^2) ρ[r]) ((ω1[r])^2)-8 r ω1[r] (ω1)'[r]+(r^2) (((ω1)'[r])^2))-r (40 π r ρ[r] ((ω1[r])^2)+(ω1)'[r] (-8 ω1[r]+r (ω1)'[r]))))-r P[r] (-24 (E^(ν[r])) π p0[r] ρ[r]+(r^2) (4 (-2 π r ((ω1[r])^2) ρ'[r]-8 π r ρ[r] ω1[r] (ω1)'[r]+(((ω1)'[r])^2))+3 ((Q[l])^2) (8 π r ((ω1[r])^2) ρ'[r]-3 (((ω1)'[r])^2)+16 π ρ[r] ω1[r] (ω1[r]+2 r (ω1)'[r]))+3 ((Q[1+l])^2) (8 π r ((ω1[r])^2) ρ'[r]-3 (((ω1)'[r])^2)+16 π ρ[r] ω1[r] (ω1[r]+2 r (ω1)'[r]))))))+M[r] (12 (E^(ν[r])) m0[r] (1+4 π (r^2) P[r]-4 π (r^2) ρ[r])-(r^3) (-24 (E^(ν[r])) π p0[r] ρ[r]+1280 (π^3) (r^6) ((P[r])^3) (-1+3 ((Q[l])^2)+3 ((Q[1+l])^2)) ((ω1[r])^2)+48 π (r^2) ((Q[l])^2) ρ[r] ((ω1[r])^2)+48 π (r^2) ((Q[1+l])^2) ρ[r] ((ω1[r])^2)+96 (E^(ν[r])) π r p0[r] ρ'[r]-8 π (r^3) ((ω1[r])^2) ρ'[r]+24 π (r^3) ((Q[l])^2) ((ω1[r])^2) ρ'[r]+24 π (r^3) ((Q[1+l])^2) ((ω1[r])^2) ρ'[r]-32 π (r^3) ρ[r] ω1[r] (ω1)'[r]+96 π (r^3) ((Q[l])^2) ρ[r] ω1[r] (ω1)'[r]+96 π (r^3) ((Q[1+l])^2) ρ[r] ω1[r] (ω1)'[r]+4 (r^2) (((ω1)'[r])^2)-9 (r^2) ((Q[l])^2) (((ω1)'[r])^2)-9 (r^2) ((Q[1+l])^2) (((ω1)'[r])^2)+64 (π^2) (r^2) ((P[r])^2) (3 (E^(ν[r])) p0[r]+(r^2) ((7-20 π (r^2) ρ[r]+3 ((Q[l])^2) (-9+20 π (r^2) ρ[r])+3 ((Q[1+l])^2) (-9+20 π (r^2) ρ[r])) ((ω1[r])^2)-8 r (-1+3 ((Q[l])^2)+3 ((Q[1+l])^2)) ω1[r] (ω1)'[r]+(r^2) (-1+3 ((Q[l])^2)+3 ((Q[1+l])^2)) (((ω1)'[r])^2)))+8 π P[r] (3 (E^(ν[r])) p0[r] (-1+8 π (r^2) ρ[r])-2 (r^2) (3 ((Q[l])^2) (((ω1[r])^2) (-1+36 π (r^2) ρ[r]+8 π (r^3) ρ'[r])+2 r (-1+16 π (r^2) ρ[r]) ω1[r] (ω1)'[r]-4 (r^2) (((ω1)'[r])^2))+3 ((Q[1+l])^2) (((ω1[r])^2) (-1+36 π (r^2) ρ[r]+8 π (r^3) ρ'[r])+2 r (-1+16 π (r^2) ρ[r]) ω1[r] (ω1)'[r]-4 (r^2) (((ω1)'[r])^2))+r (-8 π (r^2) ((ω1[r])^2) ρ'[r]+2 ω1[r] (ω1)'[r]+5 r (((ω1)'[r])^2)-4 π r ρ[r] ω1[r] (7 ω1[r]+8 r (ω1)'[r])))))))))}];

g[r_]:=Exp[((ν[r])/2)](Sqrt[1-2M[r]/r]);
(*f[r_]:=1;
P[r_]:=0;
ρ[r_]:=0;
M[r_]:=MM;
ω1[r_]:=2J/(r^3);*)
Simplify[Series[(((g[r])^2)/(g'[r]g[r]))((Coefficient[EQfDF,ψ[l]'[r]])/(Coefficient[EQfDF,ψ[l]''[r]])),{ϵa,0,2}]//.RULE]
Clear[f,P,ρ,M,ω1]

Vs=Simplify[Series[(((g[r])^2)/1)((Coefficient[EQfDF,ψ[l][r]])/(Coefficient[EQfDF,ψ[l]''[r]])),{ϵa,0,2}]//.RULE];
Vs=Simplify[Vs//.Union[RULE,D[RULE,r]]];
Clp2=Simplify[Series[(((g[r])^2)/1)((Coefficient[EQfDF,ψ[l+2][r]])/(Coefficient[EQfDF,ψ[l]''[r]])),{ϵa,0,2}]//.Union[RULE,D[RULE,r]]]
Clm2=Simplify[Series[(((g[r])^2)/1)((Coefficient[EQfDF,ψ[l-2][r]])/(Coefficient[EQfDF,ψ[l]''[r]])),{ϵa,0,2}]//.Union[RULE,D[RULE,r]]]
Clpp2=Simplify[Series[(((g[r])^2)/1)((Coefficient[EQfDF,ψ[l+2]'[r]])/(Coefficient[EQfDF,ψ[l]''[r]])),{ϵa,0,2}]//.Union[RULE,D[RULE,r]]]
Clpm2=Simplify[Series[(((g[r])^2)/1)((Coefficient[EQfDF,ψ[l-2]'[r]])/(Coefficient[EQfDF,ψ[l]''[r]])),{ϵa,0,2}]//.Union[RULE,D[RULE,r]]]

EQF=((g[r])^2)ψ[l]''[r]+g'[r]g[r]ψ[l]'[r]+Vs ψ[l][r]+(Clp2 ψ[l+2][r]+Clm2 ψ[l-2][r]+Clpp2 ψ[l+2]'[r]+Clpm2 ψ[l-2]'[r])//.Union[RULE,D[RULE,r]];
EQFnocoup=((g[r])^2)ψ[l]''[r]+g'[r]g[r]ψ[l]'[r]+Vs ψ[l][r]+0(Clp2 ψ[l+2][r]+Clm2 ψ[l-2][r]+Clpp2 ψ[l+2]'[r]+Clpm2 ψ[l-2]'[r])//.Union[RULE,D[RULE,r]]/.ψ[l]->ψ//Normal;

Simplify[Series[(EQF/((g[r])^2))-(EQfDF/(Coefficient[EQfDF,ψ[l]''[r]]))//.Union[RULE,D[RULE,r]],{ϵa,0,2}]]

Put[{EQfDF,Vs,Clp2,Clm2,Clpp2,Clpm2,EQF,EQFnocoup},"eq_scalar_final.m"];

Final potential

Quit[];

SetDirectory[NotebookDirectory[]];

VV=Get["eq_scalar_final.m"][[2]];
V0=((SeriesCoefficient[VV,0])/1)//FullSimplify
V1=((SeriesCoefficient[VV,1])/1)//Simplify
V2=((SeriesCoefficient[VV,2])/1);

Second-order correction for a constant-density star:

Note (Subscript[V, 2]) is defined slightly differently than in the review (there’s a term (e^(ν[r]))in front)

Simplify[V2//.{ρ->Function[r,ρ0],M[r]->((4π)/3) ρ0 (r^3),l->1,m->1,Q[L_]->(Sqrt[(((l^2)-(m^2))/(4(l^2)-1))])}]

In vacuum:

Simplify[V2//.{P->Function[r,0],ρ->Function[r,0],M[r]->mass,l->1,m->1,Q[L_]->(Sqrt[(((l^2)-(m^2))/(4(l^2)-1))]),ν->Function[r,Log[1-2mass/r]]}]

Numerical Integration:

Quit[];

SetDirectory[NotebookDirectory[]];
EQSGR=Get["field_eqs_GR.m"];
{E01,E02,E03,E11,E21,E22,E23,E24,E25,E26,E27,E26h,E27h}=EQSGR;
(*l0=1;*)
EQscalar=Get["eq_scalar_final.m"][[8]];

(*EQscalar=Get["eq_scalar_final.m"][[8]]//.{l->l0,m->1,Q[L_]->(Sqrt[(((L^2)-(m^2))/(4(L^2)-1))])}//Normal;
Put[EQscalar,"eq_scalar_final_l1m1.m"];*)

UNITS

kg=(1/( 1.98892 (10^30)));(*Solar mass units*)
second=(1/(4.92549232189886339689643862 (10^(-6))));(*Solar mass units*)
meter=(1/(2.99792458 (10^8)))second;
cm=(10^(-2))meter;
dyne=(((10^(-3))kg)/(second^2))cm;
cx= (((10^(-3))kg)/(cm^3));
cy=(dyne/(cm^2));
mb=1.66 (10^(-27)) kg;

cz=(cm^(-3));
gcm3=1.6188*(10^(-18));(*g/(cm^3) in Msun=1 units*)
km=0.67722;(*km in Msun=1 units*)

EOS [tabulated]

$Assumptions={η[r]≥0,η0≥0,r≥0};
SetDirectory["/home/paolo/SVN/ILoveQST/EOS/"];
fileEOS={"eosAPR","eosFPS","eosAPRb","eosA","eospoly","eosMS1"};

fileEOS=fileEOS[[1]];
INTORD=3;
(***)
EOSv=Import[fileEOS,"Table"];


(**************third point just to be safe *********************)
EOS=Join[Table[{EOSv[[i]][[2]] cy,EOSv[[i]][[1]]cx},{i,3,Length[EOSv]}]];
EOSinv=Join[Table[{EOSv[[i]][[1]]cx,EOSv[[i]][[2]] cy},{i,3,Length[EOSv]}]];
BN=Table[{EOSv[[i]][[2]] cy,EOSv[[i]][[4]] cz},{i,3,Length[EOSv]}];
(***********************************)
ρP[x_]:=Interpolation[EOS,x,InterpolationOrder->INTORD];
Pρ[x_]:=Interpolation[EOSinv,x,InterpolationOrder->INTORD];
barnum[x_]:=Interpolation[BN,x,InterpolationOrder->INTORD];
(***********************************)
ρPT=Table[{EOS[[i]][[1]],ρP[EOS[[i]][[1]]]},{i,1,Length[EOS]}];
prho=ListLogLogPlot[{EOS},Joined->True,Frame->True,PlotRange->All,Mesh->Automatic,FrameLabel->{"P","ρ(P)"}];
pbn=ListLogLogPlot[BN,Joined->True,Frame->True,PlotRange->All,Mesh->Automatic,FrameLabel->{"P","(Subscript[n, B])(P)"}];
GraphicsGrid[{{prho,pbn}},ImageSize->700]
x0=EOSv[[2]][[2]]cy;
(*rhoP[x_]:=Piecewise[{{Interpolation[EOS,x,InterpolationOrder->5],x>x0},{0,x<x0}}];*)

EOSin={ρ[r]->ρP[P[r]]};
(********************)
fileEOS
EOS[[Length[EOS]]]
ρP[(10^(-5))]
barnum[1.1 (10^(-2))]

Polytropic EOS: P(ρ)=k (ρ^(1+1/n));

rulePOLY={k->(0.0016191) (0.00370515)^(-(1+n)/n),n->0.5};
rulePOLY={k->45.86269284000001`,n->1};(*this is to agree with Jan*)
rulePOLY={k->100.,n->1};(*this is to agree with Jan*)
ρP[x_]:=((((x/k)))^(n/(1+n)))//.rulePOLY;
Pρ[x_]:=k (x^(1+1/n))//.rulePOLY;
EOSin={ρ[r]->ρP[P[r]]};

Constant density

ρP[x_]:=1;
EOSin={ρ[r]->ρP[P[r]]};

Series Origin

eqh0=h0'[r]-(((E^(-ν[r])) (12 (E^(ν[r])) (m0[r] (1+8 π (r^2) P[r])+4 π (r^2) (r-2 M[r]) p0[r] (P[r]+ρ[r]))-(r^3) (((r-2 M[r]))^2) (((ω1)'[r])^2)))/(12 (((r-2 M[r]))^2)));

EQs={E01,E02,E03,E11,E21,E22,E26,E27,E26h,E27h,eqh0};
(**************)
ORD=3;
RULEc={};
(**************)
M[r_]:=Sum[Mc[i](r^i),{i,3,ORD+5}];
ν[r_]:=Sum[νc[i](r^i),{i,0,ORD+5}];
P[r_]:=Sum[Pc[i](r^i),{i,0,ORD+5}];
ρ[r_]:=rhoP[P[r]];
ω1[r_]:=Sum[ω1c[i](r^i),{i,0,ORD+5}];
m0[r_]:=Sum[m0c[i](r^i),{i,5,ORD+5}];
p0[r_]:=Sum[p0c[i](r^i),{i,2,ORD+5}];
v2[r_]:=Sum[v2c[i](r^i),{i,4,ORD+5}];
h2[r_]:=Sum[h2c[i](r^i),{i,2,ORD+5}];
v2h[r_]:=Sum[v2hc[i](r^i),{i,4,ORD+5}];
h2h[r_]:=Sum[h2hc[i](r^i),{i,2,ORD+5}];
(***)
h0[r_]:=Sum[h0c[i](r^i),{i,3,ORD+5}];

ss0=Series[EQs//.RULEc,{r,0,ORD+3}];

EQc0=Union[Table[SeriesCoefficient[ss0[[1]],i]==0,{i,2,ORD+1}],Table[SeriesCoefficient[ss0[[2]],i]==0,{i,0,ORD}],Table[SeriesCoefficient[ss0[[3]],i]==0,{i,0,ORD}],Table[SeriesCoefficient[ss0[[4]],i]==0,{i,-1,ORD}],Table[SeriesCoefficient[ss0[[5]],i]==0,{i,1,ORD-2}],Table[SeriesCoefficient[ss0[[6]],i]==0,{i,4,ORD+1}],Table[SeriesCoefficient[ss0[[7]],i]==0,{i,3,ORD+1}],Table[SeriesCoefficient[ss0[[8]],i]==0,{i,1,ORD-1}],Table[SeriesCoefficient[ss0[[9]],i]==0,{i,3,ORD+1}],Table[SeriesCoefficient[ss0[[10]],i]==0,{i,1,ORD-1}],Table[SeriesCoefficient[ss0[[11]],i]==0,{i,3,ORD+1}]];

(***)
yc0=Union[Table[Mc[i],{i,3,ORD+2}],Table[νc[i],{i,1,ORD+1}],Table[Pc[i],{i,1,ORD+1}],Table[ω1c[i],{i,1,ORD+2}],Table[m0c[i],{i,5,ORD+2}],Table[p0c[i],{i,2,ORD-1}],Table[v2c[i],{i,4,ORD+2}],Table[h2c[i],{i,3,ORD}],Table[v2hc[i],{i,4,ORD+2}],Table[h2hc[i],{i,3,ORD}],Table[h0c[i],{i,4,ORD+2}]];
Length[EQc0]
Length[yc0]
systc0=Solve[EQc0,yc0][[1]]//Simplify;
Length[systc0];
(******************)
(*Series[{EQs[[9]]}//.systc0,{r,0,4}]//Simplify
Series[{EQs[[10]]}//.systc0,{r,0,4}]//Simplify*)
ruleSc={ψc[1]->0,ψc[2]->(((E^(-νc[0])) ψc[0] (4 (E^(νc[0])) l (3+2 l) π rhoP[Pc[0]]-3 ((ω^2)+4 (E^(νc[0])) l π Pc[0]+2 m ϵa ω (-Ω+ω1c[0]))))/(6 (3+2 l)))};

ψ[r_]:=(r^(l+1))Sum[ψc[i](r^i),{i,0,ORD+5}];
ss1=Series[Series[(EQscalar/(r^(-2+l)))//.RULEc//.systc0//.ruleSc,{r,0,2}]//Normal//Expand,{r,0,0}]

(*Series[M[r]//.systc0,{r,0,5}]//Simplify;
Series[ν[r]//.systc0,{r,0,4}]//Simplify;
Series[P[r]//.systc0,{r,0,4}]//Simplify;
Series[ω1[r]//.systc0,{r,0,4}]//Simplify;
Series[p0[r]//.systc0,{r,0,4}]//Simplify;
Series[m0[r]//.systc0,{r,0,5}]//Simplify;
Series[v2[r]//.systc0,{r,0,4}]//Simplify;
Series[h2[r]//.systc0,{r,0,2}]//Simplify;
Series[v2h[r]//.systc0,{r,0,4}]//Simplify;
Series[h2h[r]//.systc0,{r,0,2}]//Simplify;*)

Clear[M,ν,P,ω1,m0,p0,v2,h2,ρ,v2h,h2h,h0,ψ];

SetDirectory[NotebookDirectory[]];
Put[{systc0,ORD,ruleSc},"series_center_GR.m"];

INTEGRATION

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Needs["GUIKit`"];

Clear[A];
SetDirectory[NotebookDirectory[]];
INT[Pc0_,nu0_,ωc0_,Rs0_,EPS_,AG_,PG_]:=(
ruleP={AccuracyGoal->AG,PrecisionGoal->PG,MaxSteps->(10^5),MaxStepSize->Automatic};
rc=EPS;
param={Pc[0]->Pc0,νc[0]->nu0(*-0.5542070640705589`*),ω1c[0]->ωc0(*0.6950795652873236`*),h2c[2]->1,h2hc[2]->1,rhoP[x_]->ρP[x],rhoP'[x_]->ρP'[x]};

systc0=Get["series_center_GR.m"][[1]];
ORD=Get["series_center_GR.m"][[2]];
ruleSc=Get["series_center_GR.m"][[3]];
Mg=Sum[Mc[i](r^i),{i,3,ORD+2}];
νg=Sum[νc[i](r^i),{i,0,ORD+1}];
Pg=Sum[Pc[i](r^i),{i,0,ORD+1}];
ω1g=Sum[ω1c[i](r^i),{i,0,ORD+2}];
p0g=Sum[p0c[i](r^i),{i,2,ORD-1}];
m0g=Sum[m0c[i](r^i),{i,5,ORD+2}];
v2g=Sum[v2c[i](r^i),{i,4,ORD+2}];
h2g=Sum[h2c[i](r^i),{i,3,ORD}];
v2hg=Sum[v2hc[i](r^i),{i,4,ORD+2}];
h2hg=Sum[h2hc[i](r^i),{i,3,ORD}];
h0g=Sum[h0c[i](r^i),{i,4,ORD+1}];


eqh0=h0'[r]-(1/(12 (((r-2 M[r]))^2)))(E^(-ν[r])) (12 (E^(ν[r])) (m0[r] (1+8 π (r^2) P[r])+4 π (r^2) (r-2 M[r]) p0[r] (P[r]+ρ[r]))-(r^3) (((r-2 M[r]))^2) ((ω1'[r])^2));

BCs={M[rc]==Mg/.r->rc,ν[rc]==νg/.r->rc,P[rc]==Pg/.r->rc,ω1[rc]==ω1g/.r->rc,ω1'[rc]==D[ω1g,r]/.r->rc,m0[rc]==m0g/.r->rc,p0[rc]==p0g/.r->rc,v2[rc]==v2g/.r->rc,h2[rc]==h2g/.r->rc,v2h[rc]==v2hg/.r->rc,h2h[rc]==h2hg/.r->rc,h0[rc]==h0g/.r->rc}//.Union[systc0,ruleSc,param];
EQSGR=Get["field_eqs_GR.m"];
{E01,E02,E03,E11,E21,E22,E23,E24,E25,E26,E27,E26h,E27h}=EQSGR;
EQs={E01==0,E02==0,E03==0,E11==0,E21==0,E22==0,E26==0,E27==0,E26h==0,E27h==0,eqh0==0}//.Union[param,EOSin,D[EOSin,r]];
variables={M,ν,P,ω1,m0,p0,v2,h2,v2h,h2h,h0};
sol=NDSolve[Union[EQs,BCs],variables,{r,rc,Rs0},ruleP,Method->{"EventLocator","Event"->P[r]==(10^(-9))Pc0(*Pρ[ 5 (10^9)gcm3]*)}(*,Method->"ExplicitRungeKutta"*)];
Mn=M[r]/.sol[[1]];
νn=ν[r]/.sol[[1]];
Pn=P[r]/.sol[[1]];
ω1n=ω1[r]/.sol[[1]];
m0n=m0[r]/.sol[[1]];
p0n=p0[r]/.sol[[1]];
h0n=h0[r]/.sol[[1]];
v2n=v2[r]/.sol[[1]];
h2n=h2[r]/.sol[[1]];
v2hn=v2h[r]/.sol[[1]];
h2hn=h2h[r]/.sol[[1]];
(****)
ρn=ρP[Pn];

(***)


(********************)
R=InterpolatingFunctionDomain[First[M/.sol]][[1,-1]];
If[Abs[(R/Rs0)-1]<(10^(-8)),Print["WARNING!!!! NO RADIUS!!!!!"]];
(***********************************)
MASS=Mn/.r->R;
Ωs=(Sqrt[MASS/(R^3)]);
(*omren=(1/(ω1n+(R/3)D[ω1n,r]))/.r->R;(*this is 1/Ω*)*)
nurenorm=Log[1-2(MASS/R)]-νn/.r->R;
νN=νn+nurenorm;
(******)
ω1n=ω1n Ωs;
m0n=m0n ((Ωs)^2);
p0n=p0n ((Ωs)^2);
h0n=h0n ((Ωs)^2);
v2n=v2n ((Ωs)^2);
h2n=h2n ((Ωs)^2);
v2hn=v2hn ((Ωs)^2);
h2hn=h2hn ((Ωs)^2);
(******)
omrescal=(((ω1n+(R/3)D[ω1n,r]))/(Ωs))/.r->R;
g00rescal=Exp[-νn](1-2(MASS/R))/.r->R;
rescal=g00rescal (omrescal^2);
h2T=Ap h2hn+ (h2n/(rescal^1));
v2T=Ap v2hn+ (v2n/(rescal^1));
m0n=(m0n/rescal);
p0n=(p0n/rescal);
h0n=(h0n/rescal);
ω1n=((ω1n)/omrescal);

(*Print[rescal];*)
(*A=mb NIntegrate[barnum[Pn/.r->x](1-2((Mn/.r->x)/x))^(-1/2)4 π (x^2),{x,rc,R},MaxRecursion->15];*)

(*EB=A-MASS;*)

(*****)
J=((R^4)/6)D[ω1n,r]/.r->R;
Omega=ω1n+((2J)/(R^3))/.r->R;(*by construction this is (Subscript[Ω, s])*)
Inertia=J/Omega;
Rg=(Sqrt[(J/(Omega MASS))]);
(*****)
δM=m0n+((J^2)/(R^3))/.r->R;
Mt=MASS+δM;
h0n=h0n-(h0n/.r->R)-(((m0n/.r->R))/(R-2MASS));

(******)
δR=-p0n(((ρn+Pn))/(D[Pn,r]))/.r->R;
ζ2=-(-h2T-((r^2)/3)((Exp[-νn]((ω1n)^2))/g00rescal))(((ρn+Pn))/(D[Pn,r]));(*only in GR*)

systAK=Solve[{h2T==((J^2)/(MASS (R^3)))(1+(MASS/R))+KK ((3/2)((z^2)-1)Log[((z+1)/(z-1))]-((3(z^3)-5z)/((z^2)-1)))//.{r->R,z->(R/MASS)-1},v2T==((-(J^2))/(R^4))+ ((2 KK MASS)/(Sqrt[R(R-2MASS)]))((((z^2)-1))^(1/2))(((3(z^2)-2)/((z^2)-1))-(3/2)z Log[((z+1)/(z-1))])}/.{r->R,z->(R/MASS)-1},{Ap,KK}];
(*KK=AA/.[[1]];*)
QQ=(8/5)KK (MASS^3)+((J^2)/MASS)/.systAK[[1]];
ecc=(Sqrt[-3(v2T-h2T+ζ2/r)])/.r->R/.systAK[[1]];

(***)
Ibar=Inertia/(MASS^3);
χ=J/(MASS^2);
Qbar=(QQ/((MASS^3)(χ^2)));
Ibar2=Inertia/(Mt^3);Qbar2=(QQ/((Mt^3)(((J/(Mt^2)))^2)));
ruleEXT=Get["ruleEXT.m"];
paramR=Union[{mass->MASS,ω1R-> ω1n,ω1pR->D[ ω1n,r],h0R->((-m0R)/(R-2MASS)),h2R->h2T,m0R->m0n}/.r->R,systAK[[1]]];


{MASS,R,J,Ibar,Qbar}

);
Print["Done"];

Pcc=0.4;
(*Pcc=2.62;
Pcc=7.5;*)
(*Pcc=12.5;*)
(*Pcc=30;*)
(*Pcc=Pρ[3 (10^12)kg/(cm^3)]*)
Timing[INT[Pcc,0,1,25,(10^(-4)),13,13];];
Timing[INT[Pcc,nurenorm,(1/omrescal),25,(10^(-4)),13,13];];
(*{((ρn)/((10^14)gcm3 ))/.r->rc,(R/km),MASS,EB,Ωs km,Rg/R,δR/R,δM/MASS,δEB/EB,(Q/(MASS (R^2))),ecc}*)
{(R/km),MASS, ((2π)/(Ωs MASS)),(MASS/R),R/MASS,rescal,omrescal}

ff=Piecewise[{{Mn,r<R},{M[r]/.ruleEXT//.paramR,r>R}}];
ff=Piecewise[{{Exp[νn],r<R},{Exp[ν[r]]/.ruleEXT//.paramR,r>R}}];

ff=Piecewise[{{h0n,r<R},{h0[r]/.ruleEXT//.paramR,r>R}}];
ff=Piecewise[{{h2T,r<R},{h2[r]/.ruleEXT//.paramR,r>R}}];

ff=Piecewise[{{m0n,r<R},{m0[r]/.ruleEXT//.paramR,r>R}}];
ff=Piecewise[{{v2T,r<R},{v2[r]/.ruleEXT//.paramR,r>R}}];
ff=Piecewise[{{p0n,r<R},{p0[r]/.ruleEXT//.paramR,r>R}}];
ff=Piecewise[{{Omega-ω1n,r<R},{Omega-ω1[r]/.ruleEXT//.paramR,r>R}}];

Plot[{ff},{r,rc,5R},PlotRange->All]

ERGOREGION: plots

ϵ0=0.2;
θ0=π/2;
consist=0;
ERGO=Piecewise[{{Exp[νN](1+2consist (ϵ^2)h0n)-(ϵ^2)((r^2)(((Omega-ω1n))^2)((Sin[θ])^2)+consist Exp[νN]h2T(3((Sin[θ])^2)-2))/.systAK[[1]],r≤R},{Exp[ν[r]](1+2consist (ϵ^2)h0[r])-(ϵ^2)((r^2)(((Omega-ω1[r]))^2)((Sin[θ])^2)+consist Exp[ν[r]]h2[r](3((Sin[θ])^2)-2))//.Union[ruleEXT,paramR],r>R}}];
Plot[{ERGO/.ϵ->ϵ0/.θ->θ0},{r,rc,2R},Frame->True]

epsofr[RR_]:=ϵ/.FindRoot[ERGO/.θ->θ0/.r->RR R,{ϵ,0.2,0,100}];
tt=Table[{N[xx],epsofr[xx]},{xx,rc/R,2,((1.5-rc/R)/400)}];

tt0b=tt;

ListPlot[{tt0,tt1,tt2,tt3,tt4},Joined->True,Frame->True,PlotRange->{{0,1.5},{0,1}}]

Scalar perturbations

Here we integrate the KG equation on the spinning geometry to second-order in the spin. The radial and angular parts of the KG equation are separated in the section above.

OOO=2;
SS[X_]:=Normal[Series[X,{ϵa,0,OOO}]];
ΔV1=SS[Get["deltaV.m"][[1]]];
eqscal1=SS[EQscalar];
ruleSc1={ψc[1]->SS[ψc[1]/.ruleSc],ψc[2]->SS[ψc[2]/.ruleSc]};
(***)

eqscal-(eqscal/.{ω->-ω,m->-m})

(*eqscal=Get["eq_scalar_final_l1m1.m"];*)
count=0;
SCAL[w0_,l0_,m0_,ϵ0_,INF_,ORDER_]:=(
eqscal=eqscal1;
ΔV=ΔV1;
(*****)
rinf=INF R;
param2={ω->w0,l->l0,m->m0,Ω->Ωs,ϵa->ϵ0,ρ[r]->ρP[Pn],ρ'[r]->D[ρP[Pn],r],ρ''[r]->D[ρP[Pn],{r,2}],ψc[0]->1,Q[L_]->(Sqrt[(((L^2)-(m^2))/(4(L^2)-1))])};

ψg=(r^(l0+1))Sum[ψc[i](r^i),{i,0,2}];
BCs={ψ[rc]==ψg/.r->rc,ψ'[rc]==D[ψg,r]/.r->rc}//.Union[systc0,ruleSc1,param,param2];
EQs={eqscal==0}//.param//.sol[[1]]//.param2;

variables={ψ};
solS=NDSolve[Union[EQs,BCs],variables,{r,rc,R},ruleP];
ψn=ψ[r]/.solS[[1]];
ψpn=D[ψn,r];
(**** EXTERNAL ***)
BCs={ψ[R]==ψn,ψ'[R]==ψpn-((ΔV)/(Exp[ν[r]](1-2M[r]/r))) ψn}//.Union[sol[[1]],param,param2,paramR,{r->R}];(*this is to take into account a possible jump in the scalar potential*)
EQs={eqscal==0}//.ruleEXT//.Union[param,param2,paramR];
variables={ψ};
solSext=NDSolve[Union[EQs,BCs],variables,{r,R,rinf},ruleP];
ψext=ψ[r]/.solSext[[1]];
ψT=Piecewise[{{ψn,rc<r<R},{ψext,r>R}}];


If[Re[w0]>0,k=w0,k=-w0];
Ain=(((E^(I kk (r+2 mass Log[-1+(r/(2 mass))]))) (kk r YY-2 I mass YYp+I r YYp))/(2 kk r))//.Union[paramR,param,param2,{r->rinf,YY->ψext,YYp->D[ψext,r],kk->k}];
Aout=(((E^(-I kk (r+2 mass Log[-1+(r/(2 mass))]))) (kk r YY+I (2 mass-r) YYp))/(2 kk r))//.Union[paramR,param,param2,{r->rinf,YY->ψext,YYp->D[ψext,r],kk->k}];
count++;

Ain
);

SCAL[0.004950093785661908`-7.49622928240463`*^-12,1,-1,0.32,30,OOO]
SCAL[-0.004950093785661908`+7.49622928240463`*^-12,1,1,0.32,30,OOO]

Vs=Get["deltaV.m"][[2]];
Plot[-Vs/.ω->0//.Union[sol[[1]],param,param2],{r,rc,R}]

count=0;
guess=(0.0049500937860579046` -7.503765744984121`*^-12 I)/MASS;
guess=0.0344794470326424` -5.324700487091181`*^-12 I;
guess=0.3320552882816647` +1.6957281819376626`*^-9 I;
guess=1.3510630890716415` +6.747998072336702`*^-6 I;(*spin=1, l=m=1, R/M~2.6*)
fun[x_?NumberQ]:=SCAL[x,1,1,1,200,OOO]
Monitor[root=FindRoot[fun[x],{x,guess}],count]
Plot[ψT//Abs,{r,rc,5R},PlotRange->All]