Superradiance and superradiant instabilities in spinning conducting stars - v1.0 [13 April 2017]
“Superradiance in rotating stars and pulsar-timing constraints on dark photons”, Vitor Cardoso, Paolo Pani, Tien-Tien Yu (submitted to PRD, 2017)

This notebook derives and solve the linear perturbation equations for a Maxwell or Proca field on the background of a slowly-spinning conducting star in General Relativity. (or a review on the decoupling procedure in the slow-rot approximation, see https://arxiv.org/abs/1305.6759)

Quit[];

Deriving the Field equations

ClearAll[GM0,x,g,ϕ,ϕ0,Λ,Ψ,α,γ,Y,Σ,Δ,NN,ω,Δ,Σ,gff,gtt,ϕ,M,F,V,f0,f1,f2,f3,f4,Ω,A,ϕ,B,ζ,ϕ,Φ1]
$Assumptions={r≥0,θ≥0,θ≤π,J≥0,NN≥0,rh≥0,X≥0};
x={t,r,θ,φ};
gdd0={{-F[r],0,0,-X[θ]ϵa*w[r](r^2)},{0,((B[r])^(-1)),0,0},{0,0,(r^2),0},{-X[θ]ϵa*w[r](r^2),0,0,(r^2)X[θ]}};
gddPOLAR=0ϵp {{H0[r,t] Y[θ,φ],H1[r,t] Y[θ,φ],0,0},{H1[r,t] Y[θ,φ],H2[r,t] Y[θ,φ],0,0},{0,0,(r^2)KK[r,t]Y[θ,φ],0},{0,0,0,(r^2)KK[r,t]Y[θ,φ]X[θ]}};
gddAXIAL=0ϵp {{0,0,h0[r,t]Sθ[θ,φ],h0[r,t]Sφ[θ,φ]},{0,0,h1[r,t]Sθ[θ,φ],h1[r,t]Sφ[θ,φ]},{h0[r,t]Sθ[θ,φ],h1[r,t]Sθ[θ,φ],0,0},{h0[r,t]Sφ[θ,φ],h1[r,t]Sφ[θ,φ],0,0}};
A1a={0,0,((aa[r,t])/(Sin[θ]))D[Y[θ,φ],φ],-aa[r,t]Sin[θ]D[Y[θ,φ],θ]};
A1b={ff[r,t]Y[θ,φ],hh[r,t]Y[θ,φ],kk[r,t]D[Y[θ,φ],θ],kk[r,t]D[Y[θ,φ],φ]};
(****)
Ad=0{A0[r,t],0,0,0}+ϵp (A1a+A1b);

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

Basic definitions

Metric and co.

gdd=Simplify[BISERIES[Normal[Series[gdd0+gddPOLAR+gddAXIAL,{ϵ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]

Curvature tensors

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}];

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}]]

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}];

Do[Rdddd[i,j,k,l]=Simplify[BISERIES[Normal[Series[Sum[g[i,a]Ruddd[a,j,k,l],{a,1,dim}],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]]];,{i,1,dim},{j,1,dim},{k,1,dim},{l,1,dim}];

Do[ETA[i,j,k,l]=Simplify[BISERIES[Normal[Series[AA/(√(-DET))LeviCivitaTensor[dim][[i]][[j]][[k]][[l]],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]]];,{i,1,dim},{j,1,dim},{k,1,dim},{l,1,dim}];

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])]))};
Do[Tuufluid[i,j]=Simplify[Normal[Series[(P[r]+ρ[r])u[[i]]u[[j]]+guu[[i]][[j]] P[r],{ϵ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}]]

Series[u,{ϵa,0,1}]

Proca equation

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[FFuu[i,j]=Simplify[Normal[Series[Sum[FF[p,q]guu[[p]][[i]]guu[[q]][[j]],{q,1,dim},{p,1,dim}],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]],{i,1,dim},{j,1,dim}];
Do[jj[i]=Simplify[Normal[Series[Sum[(σ[r]) g[i,j]guu[[j]][[p]]guu[[b]][[q]]FF[p,q]u[[l]]g[b,l],{b,1,dim},{l,1,dim},{j,1,dim},{p,1,dim},{q,1,dim}]+ϵp ρρ[r,t]Y[θ,φ]Sum[u[[l]]g[i,l],{l,1,dim}],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]],{i,1,dim}];

Do[EQPROCA[i]=Simplify[Normal[Series[BISERIES[Sum[guu[[j]][[l]](D[FF[j,i],x[[l]]]-Sum[Christ[l,j,p]FF[p,i]+Christ[l,i,p]FF[j,p],{p,1,dim}]),{l,1,dim},{j,1,dim}]-(μ^2)A[i]+4π ee jj[i]],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]],{i,1,dim}];

(****)
EQP1=EQPROCA[1];
EQP2=EQPROCA[2];
EQP3=EQPROCA[3];
EQP4=EQPROCA[4];

DivF=Normal[Series[Sum[guu[[a]][[i]](D[EQPROCA[i],x[[a]]]-Sum[Christ[a,i,l]EQPROCA[l],{l,1,dim}]),{a,1,dim},{i,1,dim}],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]];
(*Simplify[DivF]*)

Put[{EQP1,EQP2,EQP3,EQP4,DivF,divJpartial},"Proca.m"];

Tmunu of fluid and of Maxwell/Proca field

Clear[ϕ,Φ1];
ϕ=0;
Α[X_]:=1;
V[X_]:=0;

(****)
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}]

Table[T[i,j],{i,dim},{j,dim}]

F2=Simplify[Normal[Series[Sum[guu[[i]][[p]]guu[[j]][[q]]FF[p,q]FF[i,j],{i,1,dim},{j,1,dim},{p,1,dim},{q,1,dim}],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]];
A2=Sum[guu[[l]][[m]]A[l]A[m],{l,1,dim},{m,1,dim}];
Do[TEM[i,j]=Simplify[Normal[Series[Sum[guu[[p]][[q]]FF[i,p]FF[q,j],{p,1,dim},{q,1,dim}]+α((g[i,j])/4)F2+β g[i,j]((μ^2)/2)A2+(μ^2)A[i]A[j],{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]],{i,1,dim},{j,1,dim}];

Table[TEM[i,j],{i,dim},{j,dim}]

Einstein’s eqs and background solution

Do[EQEINST[i,j]=Simplify[BISERIES[Normal[Series[Ricci[i,j]-(1/2)g[i,j] SCAL+κ(-(1/2) T[i,j]+TEM[i,j]),{ϵa,0,ORDERa},{ϵp,0,ORDERp}]]]],{i,1,dim},{j,1,dim}];

Do[DivE[i]=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}]

(*FullSimplify[DivE[2]]*)

SOL={M'[r]->4π (r^2)  ρ[r],F'[r]->((2 F[r] (M[r]+4 π (r^3) P[r]))/(r (r-2 M[r]))),P'[r]->-(((M[r]+4 π (r^3) P[r]) (P[r]+ρ[r]))/(r (r-2 M[r]))),wb''[r]->((4 π r (P[r]+ρ[r]) (4 wb[r]+r wb'[r]))/(r-2 M[r]))-(4/r)wb'[r]};
Simplify[Join[Table[EQEINST[i,j],{i,1,dim},{j,1,dim}],Table[DivE[i],{i,1,dim}]]//.Union[SOL,D[SOL,r],{w->Function[r,Ω-wb[r]],B->Function[r,1-2((M[r])/r)],κ->16π}]]

eqTOV1=M'[r]-(M'[r]//.SOL);
eqTOV2=F'[r]-(F'[r]//.SOL);
eqTOV3=P'[r]-(P'[r]//.SOL);
eqTOV4=wb''[r]-(wb''[r]//.SOL);

SOLTOV={ρ->Function[r,ρ0],M->Function[r,((4π)/3) ρ0 (r^3)],F->Function[r,(((R-2 MASS) ((((Sqrt[3-8 π (r^2) ρ0])-3 (Sqrt[3-8 π (R^2) ρ0])))^2))/(4 R (3-8 π (R^2) ρ0)))],P->Function[r,((ρ0 ((Sqrt[3-8 π (R^2) ρ0])-(Sqrt[3-8 π (r^2) ρ0])))/((Sqrt[3-8 π (r^2) ρ0])-3 (Sqrt[3-8 π (R^2) ρ0])))],ρ0->(MASS/(4(π/3)(R^3))),MASS->c R};Simplify[{eqTOV1,eqTOV2,eqTOV3,eqTOV4}//.SOLTOV]

To Kojima Form

Quit[];

ff[r_,t_]:=((u1[l][r])/1) Exp[-I ω t];
hh[r_,t_]:=((u2[l][r])/1)Exp[-I ω t];
kk[r_,t_]:=((u3[l][r])/1) Exp[-I ω t];
aa[r_,t_]:=((u4[l][r])/1)Exp[-I ω t];
ρρ[r_,t_]:=((ρEM[r])/1)Exp[-I ω t];
{EQP1,EQP2,EQP3,EQP4,DivF,divJpartial}=Get["Proca.m"];
SJ[X_]:=Normal[Series[X,{ϵa,0,1}]];

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[θ]}];

divJpartial:

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

cY=Simplify[Coefficient[EQ,Y[θ,φ]]];
cYc2θ=Simplify[Coefficient[cY,Cos[2θ]]];
cYc=Simplify[Coefficient[cY,Cos[θ]]];
cY0=Simplify[cY-cYc Cos[θ]-cYc2θ Cos[2θ]];
(***)
cYθ=Simplify[Coefficient[EQ,Derivative[1,0][Y][θ,φ]]];
cYθsθ=FullSimplify[Coefficient[cYθ,Sin[θ]]];
cYθsθcθ=Simplify[Coefficient[cYθsθ,Cos[θ]]];
cYθsθ0=FullSimplify[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][θ,φ])]
vect={cY0,cYc,cYc2θ,cYθsθ0,cYθsθcθ};
D[vect,θ]==Table[0,{i,1,Length[vect]}]/.ϵa->0//Simplify
D[vect,φ]==Table[0,{i,1,Length[vect]}]

AL=Simplify[cY0-cYc2θ];
AtL=cYc;
BL=cYθsθ0;
DL=2cYc2θ;
BtL=cYθsθcθ;

Simplify[EQ-((AL+AtL Cos[θ]+DL ((Cos[θ])^2))Y[θ,φ]+(BL+BtL Cos[θ])Sin[θ]Derivative[1,0][Y][θ,φ])]

Lorenz:

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

cY=Simplify[Coefficient[EQ,Y[θ,φ]]];
cYc2θ=Simplify[Coefficient[cY,Cos[2θ]]];
cYc=Simplify[Coefficient[cY,Cos[θ]]];
cY0=Simplify[cY-cYc Cos[θ]-cYc2θ Cos[2θ]];
(***)
cYθ=Simplify[Coefficient[EQ,Derivative[1,0][Y][θ,φ]]];
cYθsθ=FullSimplify[Coefficient[cYθ,Sin[θ]]];
cYθsθcθ=Simplify[Coefficient[cYθsθ,Cos[θ]]];
cYθsθ0=FullSimplify[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][θ,φ])]
vect={cY0,cYc,cYc2θ,cYθsθ0,cYθsθcθ};
D[vect,θ]==Table[0,{i,1,Length[vect]}]
D[vect,φ]==Table[0,{i,1,Length[vect]}]

AL=Simplify[cY0-cYc2θ];
AtL=cYc;
BL=cYθsθ0;
DL=2cYc2θ;
BtL=cYθsθcθ;

Simplify[EQ-((AL+AtL Cos[θ]+DL ((Cos[θ])^2))Y[θ,φ]+(BL+BtL Cos[θ])Sin[θ]Derivative[1,0][Y][θ,φ])]

t-component:

EQ=Simplify[SJ[(((1))^(-1))EQP1//.ruleT]]

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,Y[θ,φ]]];
cYc2θ=Simplify[Coefficient[cY,Cos[2θ]]];
cYc=Simplify[Coefficient[cY,Cos[θ]]];
cY0=Simplify[cY-cYc Cos[θ]-cYc2θ Cos[2θ]];*)
(***)
cYθ=Simplify[Coefficient[EQ,Derivative[1,0][Y][θ,φ]]];
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θ;

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

r-component:

EQ=Simplify[SJ[(((1))^(-1))EQP2//.ruleT]]

cY=Simplify[Coefficient[EQ,Y[θ,φ]]];
cYc2θ=Simplify[Coefficient[cY,Cos[2θ]]];
cYc=Simplify[Coefficient[cY,Cos[θ]]];
cY0=FullSimplify[cY-cYc Cos[θ]-cYc2θ Cos[2θ]];
(***)
cYθ=Simplify[Coefficient[EQ,Derivative[1,0][Y][θ,φ]]];
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][θ,φ])]
vect={cY0,cYc,cYc2θ,cYθsθ0,cYθsθcθ};
D[vect,θ]==Table[0,{i,1,Length[vect]}]
D[vect,φ]==Table[0,{i,1,Length[vect]}]

A1=Simplify[cY0-cYc2θ];
At1=cYc;
B1=cYθsθ0;
D1=2cYc2θ;
Bt1=cYθsθcθ;

Simplify[EQ-((A1+At1 Cos[θ]+D1 ((Cos[θ])^2))Y[θ,φ]+(B1+Bt1 Cos[θ])Sin[θ]Derivative[1,0][Y][θ,φ])]

θ component:

EQQ=Series[(Normal[SJ[EQP3]])//.ruleT,{ϵp,0,1}];
CHECK3=Simplify[SeriesCoefficient[EQQ,0]]

EQ=Simplify[(((1))^(-1))SeriesCoefficient[EQQ,1]];

cY=Expand[SJ[Coefficient[EQ,Y[θ,φ]]],Trig->True];
cYsθcθ=Coefficient[cY,Sin[θ]Cos[θ]];
cYcscθ=Coefficient[cY,Csc[θ]]//Simplify;
cYcθcotθ=Coefficient[cY,Cos[θ]Cot[θ]]//Simplify;
cYc2θcotθ=Coefficient[cY,((Cos[θ])^2)Cot[θ]]//Simplify;
cYsθ=Coefficient[cY-Sin[θ]Cos[θ]cYsθcθ,Sin[θ]]//Simplify;
cYcotθ=Coefficient[cY-cYc2θcotθ ((Cos[θ])^2)Cot[θ]-cYcθcotθ Cos[θ]Cot[θ],Cot[θ]]//Simplify;
(*CHECK*)
check=Simplify[cY-(cYsθ Sin[θ]+cYcscθ Csc[θ]+cYsθcθ Cos[θ]Sin[θ]+cYc2θcotθ ((Cos[θ])^2)Cot[θ]+cYcθcotθ Cos[θ]Cot[θ]+cYcotθ Cot[θ])]
vect={cYsθ,cYcscθ,cYsθcθ,cYc2θcotθ,cYcθcotθ,cYcotθ};
D[vect,θ]==Table[0,{i,1,Length[vect]}]
D[vect,φ]==Table[0,{i,1,Length[vect]}]

cYθ=Expand[SJ[Coefficient[EQ,Derivative[1,0][Y][θ,φ]]],Trig->True];
cYθc2θ=Coefficient[cYθ,Cos[θ],2]//Simplify;
cYθs2θ=Coefficient[cYθ,Sin[θ],2]//Simplify;
cYθ0=Simplify[cYθ-(cYθs2θ ((Sin[θ])^2)+cYθc2θ ((Cos[θ])^2))];
(*CHECK*)
check=Simplify[cYθ-((cYθ0+cYθs2θ)+(cYθc2θ-cYθs2θ)((Cos[θ])^2))]
vect={cYθ0,cYθs2θ,cYθc2θ};
D[vect,θ]==Table[0,{i,1,Length[vect]}]
D[vect,φ]==Table[0,{i,1,Length[vect]}]

Simplify[EQ-(Y[θ,φ](cYsθ Sin[θ]+cYcscθ Csc[θ]+cYsθcθ Cos[θ]Sin[θ]+cYc2θcotθ ((Cos[θ])^2)Cot[θ]+cYcθcotθ Cos[θ]Cot[θ]+cYcotθ Cot[θ])+Derivative[1,0][Y][θ,φ]((cYθ0+cYθs2θ)+(cYθc2θ-cYθs2θ)((Cos[θ])^2)))]

α0=Simplify[cYθ0+cYθc2θ];
ρ0=Simplify[-cYθc2θ+cYθs2θ];
η0=Simplify[cYsθ-cYcθcotθ];
σ0=Simplify[cYsθcθ+cYcotθ];
β0=Simplify[(I/m)(cYcscθ+cYcθcotθ)];
(*τ0=Simplify[(I/m)cYcθcotθ];*)

(****)
Simplify[EQ-((α0+ρ0 ((Sin[θ])^2))Derivative[1,0][Y][θ,φ]+((η0+σ0 Cos[θ])Sin[θ]+((-I m β0)/(Sin[θ])))Y[θ,φ])]

φ component:

EQQ=Series[(Normal[SJ[EQP4]])//.ruleT,{ϵp,0,1}];
CHECK3=Simplify[SeriesCoefficient[EQQ,0]]

EQ=Simplify[(((1))^(-1))(1/(Sin[θ]))SeriesCoefficient[EQQ,1]]

cY=Expand[SJ[Coefficient[EQ,Y[θ,φ]]],Trig->True];
cYsθcθ=Coefficient[cY,Sin[θ]Cos[θ]];
cYcscθ=Coefficient[cY,Csc[θ]]//Simplify;
cYcθcotθ=Coefficient[cY,Cos[θ]Cot[θ]]//Simplify;
cYc2θcotθ=Coefficient[cY,((Cos[θ])^2)Cot[θ]]//Simplify;
cYsθ=Coefficient[cY-Sin[θ]Cos[θ]cYsθcθ,Sin[θ]]//Simplify;
cYcotθ=Coefficient[cY-cYc2θcotθ ((Cos[θ])^2)Cot[θ]-cYcθcotθ Cos[θ]Cot[θ],Cot[θ]]//Simplify;
(*CHECK*)
check=Simplify[cY-(cYsθ Sin[θ]+cYcscθ Csc[θ]+cYsθcθ Cos[θ]Sin[θ]+cYc2θcotθ ((Cos[θ])^2)Cot[θ]+cYcθcotθ Cos[θ]Cot[θ]+cYcotθ Cot[θ])]
vect={cYsθ,cYcscθ,cYsθcθ,cYc2θcotθ,cYcθcotθ,cYcotθ};
D[vect,θ]==Table[0,{i,1,Length[vect]}]
D[vect,φ]==Table[0,{i,1,Length[vect]}]

cYθ=Expand[SJ[Coefficient[EQ,Derivative[1,0][Y][θ,φ]]],Trig->True];
cYθc2θ=Coefficient[cYθ,Cos[θ],2]//Simplify;
cYθs2θ=Coefficient[cYθ,Sin[θ],2]//Simplify;
cYθ0=Simplify[cYθ-(cYθs2θ ((Sin[θ])^2)+cYθc2θ ((Cos[θ])^2))];
(*CHECK*)
check=Simplify[cYθ-((cYθ0+cYθs2θ)+(cYθc2θ-cYθs2θ)((Cos[θ])^2))]
vect={cYθ0,cYθs2θ,cYθc2θ};
D[vect,θ]==Table[0,{i,1,Length[vect]}]
D[vect,φ]==Table[0,{i,1,Length[vect]}]

Simplify[EQ-(Y[θ,φ](cYsθ Sin[θ]+cYcscθ Csc[θ]+cYsθcθ Cos[θ]Sin[θ]+cYc2θcotθ ((Cos[θ])^2)Cot[θ]+cYcθcotθ Cos[θ]Cot[θ]+cYcotθ Cot[θ])+Derivative[1,0][Y][θ,φ]((cYθ0+cYθs2θ)+(cYθc2θ-cYθs2θ)((Cos[θ])^2)))]

α1=Simplify[cYθ0+cYθc2θ];
γ0=Simplify[-cYθc2θ+cYθs2θ];
ζ0=Simplify[cYsθ-cYcθcotθ];
λ0=Simplify[cYsθcθ+cYcotθ];
β1=Simplify[(I/m)(cYcscθ+cYcθcotθ)];
(*τ1=Simplify[(I/m)cYcθcotθ];*)
FullSimplify[(β0)-(α1)]
FullSimplify[α0+β1]

(****)
Simplify[EQ-((β0+γ0 ((Sin[θ])^2))Derivative[1,0][Y][θ,φ]+((ζ0+λ0 Cos[θ])Sin[θ]+((I m α0)/(Sin[θ])))Y[θ,φ])]

Store coefficients

Put[{AL,AtL,BL,BtL,DL,A0,At0,B0,Bt0,D0,A1,At1,B1,Bt1,D1,α0,β0,η0,ζ0,γ0,σ0,ρ0,λ0},"coefficients.m"];

Proca perturbation equations

Quit[];

SetDirectory[NotebookDirectory[]];
{AL,AtL,BL,BtL,DL,A0,At0,B0,Bt0,D0,A1,At1,B1,Bt1,D1,α0,β0,η0,ζ0,γ0,σ0,ρ0,λ0}=Get["coefficients.m"];
SJ[X_]:=Normal[Series[X,{ϵa,0,1}]];
(***)

(*CHECK of the parity*)
Clear[u1,u2,u3,u4]
u1[l][r_]:=polar u1v[r];
u2[l][r_]:=polar u2v[r];
u3[l][r_]:=polar u3v[r];
u4[l][r_]:=axial u4v[r];
D[{A0,A1,AL,α0,ζ0,ρ0,σ0,Bt0,Bt1,BtL,D0,D1,DL},axial](*POLAR FUNCTIONS*)
D[{At0,At1,AtL,B0,B1,BL,β0,η0,λ0,γ0},polar](*AXIAL FUNCTIONS*)
Clear[u1,u2,u3,u4]

Tortoise and rules:

ruleTORTOISE={};
(******)
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};

t-component (tortoise)

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;

EQf0=Simplify[eq0l+eq0lm1+eq0lp1+eq0lp2+eq0lm2];

EQtorto0=SJ[EQf0/.ruleTORTOISE]//Simplify;

r-component (tortoise)

eq1l=(A1+ ((Q[l+1])^2)(D1+l Bt1)+((Q[l])^2)(D1-(l+1)Bt1))/.rulel;
eq1lp1=(Q[l+1](At1+l B1))/.rulelm1;
eq1lm1=(Q[l](At1-(l+1)B1))/.rulelp1;
eq1lp2=(Q[l+1]Q[l+2](D1+l Bt1))/.rulelm2;
eq1lm2=(Q[l]Q[l-1](D1-(l+1)Bt1))/.rulelp2;

EQf1=Simplify[eq1l+eq1lm1+eq1lp1+eq1lp2+eq1lm2];

EQtorto1=SJ[EQf1/.ruleTORTOISE]//Simplify;

Lorenz-component (tortoise)

eqLl=(AL+ ((Q[l+1])^2)(DL+l BtL)+((Q[l])^2)(DL-(l+1)BtL))/.rulel;
eqLlp1=(Q[l+1](AtL+l BL))/.rulelm1;
eqLlm1=(Q[l](AtL-(l+1)BL))/.rulelp1;
eqLlp2=(Q[l+1]Q[l+2](DL+l BtL))/.rulelm2;
eqLlm2=(Q[l]Q[l-1](DL-(l+1)BtL))/.rulelp2;

EQfL=Simplify[eqLl+eqLlm1+eqLlp1+eqLlp2+eqLlm2];

EQtortoL=SJ[EQfL/.ruleTORTOISE]//Simplify;

θ-component (tortoise)

eq3l=(l(l+1)α0-I m ζ0 +((Q[l+1])^2)l(l ρ0+σ0)+((Q[l])^2)(l+1)((l+1)ρ0-σ0))/.rulel;
eq3lm1=Q[l](-(l+1)(η0/.rulelm1)-I m((l-1)(γ0/.rulelm1)+(λ0/.rulelm1)));
eq3lp1=Q[l+1](l(η0/.rulelp1)+I m((l+2)(γ0/.rulelp1)-(λ0/.rulelp1)));
eq3lm2=Q[l-1]Q[l](l+1)(-(l-2)(ρ0/.rulelm2)-(σ0/.rulelm2));
eq3lp2=Q[l+2]Q[l+1]l(-(l+3)(ρ0/.rulelp2)+(σ0/.rulelp2));

EQf3=Simplify[eq3l+eq3lm1+eq3lp1+eq3lp2+eq3lm2];

EQtorto3=SJ[EQf3/.ruleTORTOISE]//Simplify;

φ-component (tortoise)

eq4l=(l(l+1)β0+I m η0 +((Q[l+1])^2)l(l γ0+λ0)+((Q[l])^2)(l+1)((l+1)γ0-λ0))/.rulel;
eq4lm1=Q[l](-(l+1)(ζ0/.rulelm1)+I m((l-1)(ρ0/.rulelm1)+(σ0/.rulelm1)));
eq4lp1=Q[l+1](l(ζ0/.rulelp1)-I m((l+2)(ρ0/.rulelp1)-(σ0/.rulelp1)));
eq4lm2=Q[l-1]Q[l](l+1)(-(l-2)(γ0/.rulelm2)-(λ0/.rulelm2));
eq4lp2=Q[l+2]Q[l+1]l(-(l+3)(γ0/.rulelp2)+(λ0/.rulelp2));

EQf4=Simplify[eq4l+eq4lm1+eq4lp1+eq4lp2+eq4lm2];

EQtorto4=SJ[EQf4/.ruleTORTOISE]//Simplify;

Store eqs

Put[{EQf0,EQf1,EQf3,EQf4,EQfL,EQtorto0,EQtorto1,EQtorto3,EQtorto4,EQtortoL},"final_eqs_axialpolar.m"];

Massaging Proca

Quit[];

{EQf0,EQf1,EQf3,EQf4,EQfL,EQtorto0,EQtorto1,EQtorto3,EQtorto4,EQtortoL}=Get["final_eqs_axialpolar.m"];

SOLbkg={M'[r]->4π (r^2)  ρ[r],F'[r]->((2 F[r] (M[r]+4 π (r^3) P[r]))/(r (r-2 M[r]))),P'[r]->-(((M[r]+4 π (r^3) P[r]) (P[r]+ρ[r]))/(r (r-2 M[r]))),wb''[r]->((4 π r (P[r]+ρ[r]) (4 wb[r]+r wb'[r]))/(r-2 M[r]))-(4/r)wb'[r]};
SOLbkg=Union[SOLbkg,D[SOLbkg,r],{w->Function[r,Ω-wb[r]],B->Function[r,1-2((M[r])/r)],κ->16π}];

Axial

RULEA={u1[n_]->Function[r,0],u2[n_]->Function[r,0],u3[n_]->Function[r,0]};
EQa=(EQf4/(-(((E^(-I t ω)) l (1+l))/(F[r]))))//.RULEA//Simplify

Coefficient[(EQa/(B[r]F[r])),u4[l]''[r]]
Coefficient[(EQa/((B'[r]F[r]+B[r]F'[r])/2)),u4[l]'[r]]//Simplify
V=Coefficient[(EQa/1),u4[l][r]]//Simplify

(*V=((4 I ee π (r^2) σ (ω-m ϵa Ω) (Sqrt[F[r]])-(l+(l^2)+(r^2) (μ^2)) F[r]+(r^2) ω (ω-2 m ϵa w[r]))/(r^2));*)

Simplify[(((ω^2)-2m ϵa ω w[r]-V)/(F[r]))-(((l(l+1))/(r^2))+(μ^2))]