(*This notebook separates and decouples scalar, vector and gravitational perturbations in the Schwarzschild geometry.*)

Scalar and Vector

Metric and coordinates

g=DiagonalMatrix[{-(1-2*M/r),(1-2*M/r)^(-1),r^2,(r^2)*(Sin[θ])^2}];

(*Equation for Spherical Harmonics*)par={Y''[θ]->-Cos[θ]/Sin[θ]* Y'[θ]-l*(l+1)*Y[θ]+m^2/Sin[θ]^2*Y[θ],Derivative[3][Y][θ]->Cot[θ] (l+(l^2)-3 (m^2) ((Csc[θ])^2)) Y[θ]+(-l (1+l)+((Cot[θ])^2)+(1+(m^2)) ((Csc[θ])^2)) Y'[θ]};


(*A is contravariant metric*)A=Inverse[g];

(*Γ are metric connections, first index is contravariant*)

X={t,r,θ,ϕ};

Γ=Table[Sum[1/2A[[c,d]]*(D[g[[b,d]],X[[a]]]+D[g[[a,d]],X[[b]]]-D[g[[a,b]],X[[d]]]),{d,1,4}],{c,1,4},{a,1,4},{b,1,4}];

(*Riemann tensor with all indices down: Rie[[a,b,c,d]]*)Rie=Table[Sum[g[[f,a]](D[Γ[[f,b,d]],X[[c]]]-D[Γ[[f,b,c]],X[[d]]]+Sum[Γ[[e,b,d]]Γ[[f,e,c]],{e,1,4}]-Sum[Γ[[e,b,c]]Γ[[f,e,d]],{e,1,4}]),{f,1,4}],{a,1,4},{b,1,4},{c,1,4},{d,1,4}];

(*Ricci tensor covariant (indices down)*)Riccit=Table[Sum[D[Γ[[a,b,d]],X[[a]]]-D[Γ[[a,b,a]],X[[d]]]+Sum[Γ[[e,b,d]]Γ[[a,e,a]],{e,1,4}]-Sum[Γ[[e,b,a]]Γ[[a,e,d]],{e,1,4}],{a,1,4}],{b,1,4},{d,1,4}];

(*Ricci scalar*)Riccie=Sum[A[[a,b]]Riccit[[a,b]],{a,1,4},{b,1,4}];

Simplify[Riccie]

(*Covariant Einstein tensor*)G=Table[Riccit[[a,b]]-1/2*g[[a,b]]*Riccie,{a,1,4},{b,1,4}];

Klein-Gordon Field

(*Expansion of field in spherical harmonics*)Psi=Etax[r]/r*Exp[-I*w*t+I*m*ϕ]*Y[θ];

KG=Sum[D[Sqrt[-Det[g]]*A[[a,b]]*D[Psi,X[[a]]],X[[b]]],{a,1,4},{b,1,4}];

FullSimplify[KG/.par]

Maxwell set II

SPH = Y[θ]*Exp[I*m*ϕ];
AMaxaxial={0,0,alm[t,r]/Sin[θ]*D[Y[θ]*Exp[I*m*ϕ],ϕ],-alm[t,r]*Sin[θ]*D[Y[θ]*Exp[I*m*ϕ],θ]};

AMaxpolar = {fe[t, r]*SPH, he[t, r]*SPH, klm[t, r]*D[SPH, θ], klm[t, r]*D[SPH, ϕ]};

AMax=AMaxaxial+AMaxpolar;


(*COVARIANT MAXWELL TENSOR:*)FMax = Table[D[AMax[[mu]], X[[nu]]] - D[AMax[[nu]], X[[mu]]], {mu, 1, 4}, {nu, 1, 4}];

FMaxcont = Table[Sum[A[[a, b]]*A[[c, d]]*FMax[[b, d]], {b, 1, 4}, {d, 1, 4}], {a, 1, 4}, {c, 1, 4}];

(*MAXWELL equations*)MEq = Table[Sum[D[FMaxcont[[a, b]], X[[a]]], {a, 1, 4}] + Sum[Γ[[a, d, a]]*FMaxcont[[d, b]] + Γ[[b, d, a]]*FMaxcont[[a, d]], {a, 1, 4}, {d, 1, 4}], {b, 1, 4}];

MEq2 = FullSimplify[MEq /. par]

(*Clearly, the angular dependence of the axial is different and must be set to zero separately*)Collect[MEq2[[4]],Y'[θ]]

EqAxial=(((E^(I m ϕ)) Csc[θ] Y'[θ] ((2 M-r) (l (1+l) alm[t,r]-2 M Derivative[0,1][alm][t,r]+(2 M-r) r Derivative[0,2][alm][t,r])-(r^3) Derivative[2,0][alm][t,r]))/((r^4) (-2 M+r)));

alm2=Solve[EqAxial==0,Derivative[2,0][alm][t,r]]

(*In fact, alm disappears now from the equations*)MEq2polar=Simplify[MEq2/.alm2[[1]]]

(*let's find the coefficients for which k drops out*)MEqa = Simplify[D[AA[r]*MEq2polar[[1]], r] + BB[r] D[MEq2polar[[2]], t]]

Collect[MEqa,Derivative[1,0][klm][t,r]]

DSolve[ (((l (1+l) r AA[r])/(2 M-r))-((l (1+l) (r^2) AA[r])/(((-2 M+r))^2))-((l (1+l) (r^2) AA'[r])/(2 M-r))) ==0,AA,r]

MEqb = Simplify[D[r (-2 M + r)*MEq2[[1]], r] + BB[r] D[MEq2[[2]], t]]

Collect[MEqb,Derivative[1,1][klm][t,r]]

Solve[(((E^(I m ϕ)) (l (r^3)+(l^2) (r^3)+2 l M BB[r]+2 (l^2) M BB[r]-l r BB[r]-(l^2) r BB[r]) Y[θ] Derivative[1,1][klm][t,r])/(r^3))==0,BB[r]]

MEqc=FullSimplify[D[r (-2 M+r)*MEq2[[1]],r]-((r^3)/(2 M-r))D[MEq2[[2]],t]]

Derivative[1,0][he][t,r]=l*(l+1)*PSI[t,r]/r^2+Derivative[0,1][fe][t,r];
Derivative[1,1][he][t,r]=D[l*(l+1)*PSI[t,r]/r^2+Derivative[0,1][fe][t,r],r];
Derivative[1,2][he][t,r]=D[l*(l+1)*PSI[t,r]/r^2+Derivative[0,1][fe][t,r],{r,2}];
Derivative[3,0][he][t,r]=D[l*(l+1)*PSI[t,r]/r^2+Derivative[0,1][fe][t,r],{t,2}];

EqPolar=FullSimplify[MEqc]

Bianchi=Simplify[Table[D[FMax[[mu,nu]],X[[lambda]]]+D[FMax[[lambda,mu]],X[[nu]]]+D[FMax[[nu,lambda]],X[[mu]]],{mu,1,ND},{nu,1,ND},{lambda,1,ND}]]

Gravitational perturbations

Axial

m=0;

Z=Y[θ]*Exp[-I*w*t+I*m*ϕ];

h1[r_]:=Qaxial[r]*r/(1-2*M/r);

g0=DiagonalMatrix[{-(1-2*M/r),(1-2*M/r)^(-1),r^2,(r^2)*(Sin[θ])^2}];

g=g0+eps*{{0,0,-h0[r]/Sin[θ]D[Z,ϕ],h0[r]*Sin[θ]D[Z,θ]},{0,0,-h1[r]/Sin[θ]D[Z,ϕ],h1[r]*Sin[θ]D[Z,θ]},{-h0[r]/Sin[θ]D[Z,ϕ],-h1[r]/Sin[θ]D[Z,ϕ],0,0},{h0[r]*Sin[θ]D[Z,θ],h1[r]*Sin[θ]D[Z,θ],0,0}};

(*Equation for Spherical Harmonics*)par={Y''[θ]->-Cos[θ]/Sin[θ]* Y'[θ]-l*(l+1)*Y[θ]+m^2/Sin[θ]^2*Y[θ],Derivative[3][Y][θ]->Cot[θ] (l+(l^2)-3 (m^2) ((Csc[θ])^2)) Y[θ]+(-l (1+l)+((Cot[θ])^2)+(1+(m^2)) ((Csc[θ])^2)) Y'[θ]};


(*A is contravariant metric*)A=Inverse[g];

(*Γ are metric connections, first index is contravariant*)

X={t,r,θ,ϕ};

Γ=Table[Sum[1/2A[[c,d]]*(D[g[[b,d]],X[[a]]]+D[g[[a,d]],X[[b]]]-D[g[[a,b]],X[[d]]]),{d,1,4}],{c,1,4},{a,1,4},{b,1,4}];
Γ=Normal[Series[Γ,{eps,0,1}]];
(*Riemann tensor with all indices down: Rie[[a,b,c,d]]*)Rie=Table[Sum[g[[f,a]](D[Γ[[f,b,d]],X[[c]]]-D[Γ[[f,b,c]],X[[d]]]+Sum[Γ[[e,b,d]]Γ[[f,e,c]],{e,1,4}]-Sum[Γ[[e,b,c]]Γ[[f,e,d]],{e,1,4}]),{f,1,4}],{a,1,4},{b,1,4},{c,1,4},{d,1,4}];

(*Ricci tensor covariant (indices down)*)Riccit=Table[Sum[D[Γ[[a,b,d]],X[[a]]]-D[Γ[[a,b,a]],X[[d]]]+Sum[Γ[[e,b,d]]Γ[[a,e,a]],{e,1,4}]-Sum[Γ[[e,b,a]]Γ[[a,e,d]],{e,1,4}],{a,1,4}],{b,1,4},{d,1,4}];
Riccit=Normal[Series[Riccit,{eps,0,1}]];

(*Ricci scalar*)Riccie=Sum[A[[a,b]]Riccit[[a,b]],{a,1,4},{b,1,4}];

MatrixForm[g-g0]

(*Covariant Einstein tensor*)G=Table[Riccit[[a,b]]-1/2*g[[a,b]]*Riccie,{a,1,4},{b,1,4}];

G2=Simplify[Normal[Series[G,{eps,0,1}]/.par]]

h0sol=Solve[G2[[3,4]]==0,h0[r]][[1]]

h0p=Simplify[D[h0sol,r]]

FullSimplify[G2/.h0sol/.h0p]

Polar

Done in "Schwarzschild_infall_PP_ICTP_Cardoso.nb"with a source term

References

[1]V. Cardoso, L. Gualtieri, C. Herdeiro, U. Sperhake, "Exploring New Physics Frontiers Through Numerical Relativity", Living Rev.Relativity 18 (2015) 1 
[2](The classical original work)T. Regge and J. A. Wheeler, "Stability of a Schwarzschild singularity", Phys.Rev.108 (1957) 1063-1069

[3] F. Zerilli,  "Gravitational field of a particle falling in a schwarzschild geometry analyzed in tensor harmonics", Phys.Rev.D2 (1970) 2141-2160

[6]https://centra.tecnico.ulisboa.pt/network/grit/files/