dRGT linear field equations

This notebook calculates the field equations for the linearisation of the dRGT massive gravity theory (arXiv:1208.1515). We do so through the Mathematica package xAct (http://www.xact.es/), which must be installed to run the notebook. These calculations were performed for the paper “Gravitational waves in massive gravity theories: waveforms, fluxes and constraints from extreme-mass-ratio mergers”, by Vítor Cardoso, Gonçalo Castro and Andrea Maselli, still to appear.

We would like to mention the notebook RW_Zerilli_Tutorial.nb fromSayan Chakrabarti, Jan Steinhoff and Jorge Rocha, which was of great help to introduce us to the problem of the perturbation of a Schwarzschild black hole.

Sections 1 and 2 must always be run. Afterwards, one can compute the equations for any of the multipolar orders in sections 3 and 4, afterwards obtaining its source terms in section 5, either for a radial infall or a circular orbit.

We will eventually update the notebook so that each section can be run independently, through the use of some .m files generated in a first run of the notebook.

Section 1 - Basic Definitions

<<xAct`xCoba`;
<<xAct`xPert`;

DefManifold[M,4,IndexRange[a,l]];
(*DefMetric[-1,metric[-a,-b],CD,{";","▽"},PrintAs->"g"];*)

DefChart[swc,M,{0,1,2,3},{t[],r[],θ[],ϕ[]}];

DefConstantSymbol[ϵ];

linearize[expr_]:=Collect[Series[expr,{ϵ,0,1}]//Normal,ϵ,Together];

DefConstantSymbol/@{ω,ll,m};
(*m=0;*)
DefScalarFunction[PL];
PL''[θ_]:=-(ll-m^2 Csc[θ]^2) PL[θ]-Cot[θ]PL'[θ];
YLM:=E^(I m ϕ[]-I ω t[])PL[θ[]];
XLM:=2D[D[YLM,θ[]]-Cot[θ[]]YLM,ϕ[]]
WLM:=D[YLM,{θ[],2}]-Cot[θ[]]D[YLM,θ[]]-Csc[θ[]]^2D[YLM,{ϕ[],2}]

DefScalarFunction/@{F};
metricBG=DiagonalMatrix[{-F[r[]],1/F[r[]],r[]^2,r[]^2Sin[θ[]]^2}];
metricS=CTensor[DiagonalMatrix[{-F[r[]],1/F[r[]],r[]^2,r[]^2Sin[θ[]]^2}],{-swc,-swc}];
DefConstantSymbol[Mass];
backg={F->Function[r,1-2Mass/r]};
DefScalarFunction/@{h0,h1,h2,H0,H1,H2,η1,η0,G,K};
DefConstantSymbol[pol];
(*this is the definition component by component*)
alm0={{YLM,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}};
alm1=I/Sqrt[2]{{0,YLM,0,0},{YLM,0,0,0},{0,0,0,0},{0,0,0,0}};
alm={{0,0,0,0},{0,YLM,0,0},{0,0,0,0},{0,0,0,0}};
blm0=I r[]1/Sqrt[2ll]{{0,0,D[YLM,θ[]],D[YLM,ϕ[]]},{0,0,0,0},{D[YLM,θ[]],0,0,0},{D[YLM,ϕ[]],0,0,0}};
blm=r[]1/Sqrt[2ll]{{0,0,0,0},{0,0,D[YLM,θ[]],D[YLM,ϕ[]]},{0,D[YLM,θ[]],0,0},{0,D[YLM,ϕ[]],0,0}};
glm=r[]^2/Sqrt[2]{{0,0,0,0},{0,0,0,0},{0,0,YLM,0},{0,0,0,Sin[θ[]]^2YLM}};
flm=r[]^2/Sqrt[2ll(ll-2)]{{0,0,0,0},{0,0,0,0},{0,0,WLM,XLM},{0,0,XLM,-Sin[θ[]]^2WLM}};
dlm=I r[]^2 /Sqrt[2ll(ll-2)]{{0,0,0,0},{0,0,0,0},{0,0,-Csc[θ[]]XLM,Sin[θ[]]WLM},{0,0,Sin[θ[]]WLM,Sin[θ[]]XLM}};
clm=I r[]1/Sqrt[2ll]{{0,0,0,0},{0,0,Csc[θ[]]D[YLM,ϕ[]],-Sin[θ[]]D[YLM,θ[]]},{0,Csc[θ[]]D[YLM,ϕ[]],0,0},{0,-Sin[θ[]]D[YLM,θ[]],0,0}};
clm0=r[]1/Sqrt[2ll]{{0,0,Csc[θ[]]D[YLM,ϕ[]],-Sin[θ[]]D[YLM,θ[]]},{0,0,0,0},{Csc[θ[]]D[YLM,ϕ[]],0,0,0},{-Sin[θ[]]D[YLM,θ[]],0,0,0}};
(*this is the definition in a single formula; the polar is altered so I can try something*)
metricPertAxialF=ϵ(-Sqrt[2ll]/r[]h0[r[]]clm0+I Sqrt[2ll]/r[]h1[r[]]clm+Sqrt[2ll(ll-2)]/(2r[])h2[r[]]dlm);
metricPertPolarF=ϵ pol(F[r[]]H0[r[]]alm0-I Sqrt[2]H1[r[]]alm1+H2[r[]]/F[r[]]alm-I/r[]Sqrt[2ll]η0[r[]]blm0+1/r[]Sqrt[2ll]η1[r[]]blm+(Sqrt[2]K[r[]]-ll/Sqrt[2]G[r[]])glm+Sqrt[ll(ll-2)/2]G[r[]]flm);

(*MetricInBasis[metric,-swc,metricBG];*)
SetCMetric[metricS,swc,SignatureOfMetric->{3,1,0}];

MetricCompute[metricS,swc,All];

cd=CovDOfMetric[metricS];

Section 2 - Field Equation

(*MetricCompute[metricFull,RW,All,CVSimplify->linearize]*)

(*GEinstein=Simplify[linearize[EinsteinCD[-b,-d]//ContractMetric//ToBasis[RW]//ComponentArray//ToValues//ToValues]/.backg/.r[]->r/.t[]->t/.θ[]->θ/.ϕ[]->ϕ]*)

DefScalarFunction/@{ALM1,ALM0,ALM,BLM,BLM0,FLM,GLM,QLM,QLM0,DLM};

H=CTensor[metricPertPolarF+metricPertAxialF,{-swc,-swc}];

metricTrace=H[a,-a];

SETensorAxial=ϵ*(QLM0[r[]]clm0+QLM[r[]]clm+DLM[r[]]dlm);
SETensorPolar=ϵ pol(ALM0[r[]]alm0+ALM[r[]]alm+ALM1[r[]]alm1+BLM[r[]]blm+BLM0[r[]]blm0+GLM[r[]]glm+FLM[r[]]flm);

(*SETensorAxial=ϵ *{{0,0,QLM0[r[]]Csc[θ[]]D[YLM,ϕ[]],-QLM0[r[]]Sin[θ[]]D[YLM,θ[]]},{0,0,QLM[r[]]Csc[θ[]]D[YLM,ϕ[]],-QLM[r[]]Sin[θ[]]D[YLM,θ[]]},{QLM0[r[]]Csc[θ[]]D[YLM,ϕ[]],QLM[r[]]Csc[θ[]]D[YLM,ϕ[]],-DLM[r[]]Csc[θ[]]XLM,DLM[r[]]Sin[θ[]]WLM},{-QLM0[r[]]Sin[θ[]]D[YLM,θ[]],-QLM[r[]]Sin[θ[]]D[YLM,θ[]],DLM[r[]]Sin[θ[]]WLM,DLM[r[]]Sin[θ[]]XLM}};
SETensorPolar=ϵ pol*{{ALM0[r[]]YLM,ALM1[r[]]YLM,BLM0[r[]]D[YLM,θ[]],BLM0[r[]]D[YLM,ϕ[]]},{ALM1[r[]]YLM,ALM[r[]]YLM,BLM[r[]]D[YLM,θ[]],BLM[r[]]D[YLM,ϕ[]]},{BLM0[r[]]D[YLM,θ[]],BLM[r[]]D[YLM,θ[]],r[]^2(GLM[r[]]YLM+FLM[r[]]WLM),r[]^2FLM[r[]]XLM},{BLM0[r[]]D[YLM,ϕ[]],BLM[r[]]D[YLM,ϕ[]],r[]^2GLM[r[]]XLM,r[]^2Sin[θ[]]^2(GLM[r[]]YLM-FLM[r[]]WLM)}};*)

SETensor=CTensor[SETensorPolar+SETensorAxial,{-swc,-swc}];

SETrace=SETensor[a,-a];

EFErhs=8*π(SETensor-1/3*SETrace*metricS)/.r[]->r/.t[]->t/.θ[]->θ/.ϕ[]->ϕ/.backg;

eqgaugelhs=(cd[a][H[-a,-b]]//ComponentArray)/.r[]->r/.t[]->t/.θ[]->θ/.ϕ[]->ϕ/.backg;

eqgaugelhs[[0]][{0,-swc}];

eqgaugerhs=(-16π/(3μ^2)cd[a][SETrace])//ComponentArray;

DefChart[RW,M,{0,1,2,3},{t[],r[],θ[],ϕ[]}];
(*MetricInBasis[metric,-RW,metricFull]*)
metricFull=CTensor[metricBG+metricPertPolarF+metricPertAxialF,{-RW,-RW}];
SetCMetric[metricFull,RW,SignatureOfMetric->{3,1,0},CVSimplify->linearize];

MetricCompute[metricFull,RW,"Einstein"[-1,-1],CVSimplify->linearize]//AbsoluteTiming

cdF=CovDOfMetric[metricFull];

GEinstein=Simplify[Einstein[cdF]];

EFElhs=Table[(GEinstein[{ii,-RW},{jj,-RW}]+μ^2/2*H[{ii,-swc},{jj,-swc}])/.r[]->r/.t[]->t/.θ[]->θ/.ϕ[]->ϕ/.backg,{ii,0,3},{jj,0,3}];

DefScalarFunction[rT,PrintAs->"(SubStar[r])"];
rTdef={rT'[r_]-> 1/(1-2Mass/r),rT''[r_]-> D[1/(1-2Mass/r),r],rT'''[r_]-> D[1/(1-2Mass/r),{r,2}],rT''''[r_]-> D[1/(1-2Mass/r),{r,3}]};

EFEs={EFElhs,EFErhs,eqgaugelhs[[0]],eqgaugerhs[[0]],metricTrace,SETrace};
Export["./EFE_dRGT.m",EFEs];

Section 3 - Polar Equation

3.0 - Initialization

Clear[m,ll];
Clear[ALM,ALM1,ALM0,BLM,BLM0,GLM,FLM,QLM,QLM0,DLM];
DefScalarFunction/@{ALM1,ALM0,ALM,BLM,BLM0,FLM,GLM,QLM,QLM0,DLM};

Clear[EFElhs,EFErhs,eqgaugelhs,eqgaugerhs,metricTrace,SETrace];
{EFElhs,EFErhs,gaugelhs,gaugerhs,metricTrace,SETrace}=Import["./EFE_dRGT.m"];

EFEpolarlhs=Coefficient[EFElhs,pol,1]//Together;
EFEpolarrhs=Coefficient[EFErhs,pol,1]//Together;

eqttpol=EFEpolarlhs[[1,1]]==Coefficient[EFErhs[{0,-swc},{0,-swc}],pol,1];
eqtrpol=EFEpolarlhs[[1,2]]==Coefficient[EFErhs[{0,-swc},{1,-swc}],pol,1];
eqtthpol=Simplify[EFEpolarlhs[[1,3]]==Coefficient[EFErhs[{0,-swc},{2,-swc}],pol,1]];
eqtphpol=EFEpolarlhs[[1,4]]==Coefficient[EFErhs[{0,-swc},{3,-swc}],pol,1];
eqrrpol=Simplify[EFEpolarlhs[[2,2]]==Coefficient[EFErhs[{1,-swc},{1,-swc}],pol,1]];
eqrthpol=EFEpolarlhs[[2,3]]==Coefficient[EFErhs[{1,-swc},{2,-swc}],pol,1];
eqrphpol=EFEpolarlhs[[2,4]]==Coefficient[EFErhs[{1,-swc},{3,-swc}],pol,1];
eqththpol=(EFEpolarlhs[[3,3]]==Coefficient[EFErhs[{2,-swc},{2,-swc}],pol,1])/.G[r]->0;
eqthphpol=EFEpolarlhs[[3,4]]==Coefficient[EFErhs[{2,-swc},{3,-swc}],pol,1];
eqphphpol=EFEpolarlhs[[4,4]]==Coefficient[EFErhs[{3,-swc},{3,-swc}],pol,1];

eqpoltrace=Coefficient[metricTrace==-16π/(3μ^2)*SETrace,pol,1]/.r[]->r/.t[]->t/.θ[]->θ/.ϕ[]->ϕ/.backg;
eqpolgauget=Coefficient[gaugelhs[{0,-swc}]==gaugerhs[{0,-swc}],pol,1]/.r[]->r/.t[]->t/.θ[]->θ/.ϕ[]->ϕ/.backg;
eqpolgauger=Coefficient[gaugelhs[{1,-swc}]==gaugerhs[{1,-swc}],pol,1]/.r[]->r/.t[]->t/.θ[]->θ/.ϕ[]->ϕ/.backg;
eqpolgaugeth=Coefficient[gaugelhs[{2,-swc}]==gaugerhs[{2,-swc}],pol,1]/.r[]->r/.t[]->t/.θ[]->θ/.ϕ[]->ϕ/.backg;
eqpolgaugeph=Coefficient[gaugelhs[{3,-swc}]==gaugerhs[{3,-swc}],pol,1]/.r[]->r/.t[]->t/.θ[]->θ/.ϕ[]->ϕ/.backg;

eqspolar={{eqttpol,eqtrpol,eqtthpol,eqtphpol,eqrrpol,eqrthpol,eqrphpol,eqththpol,eqthphpol,eqphphpol},{eqpoltrace,eqpolgauget,eqpolgauger,eqpolgaugeth,eqpolgaugeph}};
Export["./pol_dRGT.m",eqspolar];

Section 3.1 - l=0

Clear[eqttpol,eqtrpol,eqtthpol,eqtphpol,eqrrpol,eqrthpol,eqrphpol,eqththpol,eqthphpol,eqphphpol,eqpoltrace,eqpolgauget,eqpolgauger,eqpolgaugeth,eqpolgaugeph];
{{eqttpol,eqtrpol,eqtthpol,eqtphpol,eqrrpol,eqrthpol,eqrphpol,eqththpol,eqthphpol,eqphphpol},{eqpoltrace,eqpolgauget,eqpolgauger,eqpolgaugeth,eqpolgaugeph}}=Import["./pol_dRGT.m"];

ll=0;
m=0;

sH0=Simplify[Solve[eqpolgauget,H0[r]]//First//First];
repH0=H0->Function[r,H0[r]/.sH0//Evaluate];

sH2=Simplify[H2[r]/.(Solve[eqpoltrace/.repH0,H2[r]]//First//First)];
repH2=H2->Function[r,sH2//Evaluate];

rep02={repH2,repH0};

eqa=Collect[eqttpol/.rep02,{H1[r],K[r],H1'[r],K'[r],H1''[r],K''[r]},Simplify];
eqb=Collect[eqtrpol/.rep02,{H1[r],K[r],H1'[r],K'[r],H1''[r],K''[r]},Simplify];
eqc=Collect[eqrrpol/.rep02,{H1[r],K[r],H1'[r],K'[r],H1''[r],K''[r]},Simplify];
eqd=Collect[Coefficient[eqththpol/.rep02,PL[θ]],{H1[r],K[r],H1'[r],K'[r],H1''[r],K''[r]},Simplify];

repH1=Simplify[Solve[{eqa,eqb,(eqpolgauger/.rep02)},{H1[r],H1'[r],H1''[r]}]]//First;
eqfinal1=Collect[eqc/.repH1,{K[r],K'[r],K''[r]},Simplify];
Clear[g];
eqfinal2=eqfinal1/.K->Function[r,g[r]*KK[r]]/.KK->Function[r,KK[rT[r]]]/.rTdef;
repg=g->Function[r,g[r]/.((DSolve[Simplify[Coefficient[eqfinal2,KK'[rT[r]]]],g[r],r]/.C[1]->1)//First//First)//Evaluate];
eqfinal3=Collect[eqfinal2/.repg/.backg,{KK[rT[r]],KK'[rT[r]],KK''[rT[r]]},Simplify];
(*,ALM[r],ALM'[r],ALM0[r],ALM0'[r],ALM1[r],ALM1'[r],GLM[r],GLM'[r]*)
eqpoll0=Collect[eqfinal3/Coefficient[eqfinal3//First,KK''[rT[r]]],{KK[rT[r]],KK'[rT[r]],KK''[rT[r]]},Simplify];
Export["./pol_l0_dRGT.m",eqpoll0];

To obtain the formula for (dE/(dω)) in terms of the master functionϕ0:

K[r_]:=μ^2/r Exp[I Sqrt[ω^2-μ^2]r]ϕ0
a1=Simplify[Series[H0[r]/.repH0/.repH1/.{GLM[r]->0,ALM0[r]->0,ALM[r]->0,ALM'[r]->0,ALM0'[r]->0,ALM1[r]->0,GLM'[r]->0}/.Mass->0,{r,∞,1}]];
a2=Simplify[Series[H1[r]/.repH1/.{GLM[r]->0,ALM0[r]->0,ALM[r]->0,ALM'[r]->0,ALM0'[r]->0,ALM1[r]->0,GLM'[r]->0}/.Mass->0,{r,∞,1}]];
a3=Simplify[Series[H2[r]/.repH2/.repH1/.{GLM[r]->0,ALM0[r]->0,ALM[r]->0,ALM'[r]->0,ALM0'[r]->0,ALM1[r]->0,GLM'[r]->0}/.Mass->0,{r,∞,1}]];
Simplify[a1^2-2a2^2+a3^2+2K[r]^2]
Clear[K];

Section 3.2 - l=1

ll=2;
Clear[eqttpol,eqtrpol,eqtthpol,eqtphpol,eqrrpol,eqrthpol,eqrphpol,eqththpol,eqthphpol,eqphphpol,eqpoltrace,eqpolgauget,eqpolgauger,eqpolgaugeth,eqpolgaugeph];
{{eqttpol,eqtrpol,eqtthpol,eqtphpol,eqrrpol,eqrthpol,eqrphpol,eqththpol,eqthphpol,eqphphpol},{eqpoltrace,eqpolgauget,eqpolgauger,eqpolgaugeth,eqpolgaugeph}}=Import["./pol_dRGT.m"];

eqpolgauget

Clear[eqa,eqb,eqc,eqd,eq1,eq2,eq3,eq4,eq5,eq6,eq7];
Clear[m,ll];
(*ll=2;*)
listfl1={H0[r],H1[r],H2[r],η0[r],η1[r],K[r],H0'[r],H1'[r],H2'[r],η0'[r],η1'[r],K'[r],H0''[r],H1''[r],H2''[r],η0''[r],η1''[r],K''[r]};

repG={G[r]->0,G'[r]->0,G''[r]->0};
repFLM={FLM[r]->0,FLM'[r]->0,FLM''[r]->0};

soltrace=Simplify[Solve[eqpoltrace/.repG/.repFLM,H0[r]]]//First;
reptrace=H0->Function[r,(H0[r]/.soltrace)//Evaluate];

solgauge=Simplify[Solve[{eqpolgauger/.reptrace/.repG/.repFLM,eqpolgaugeth/.reptrace/.repG/.repFLM},{η0[r],H1[r]}]]//First;
repgauge={η0->Function[r,(η0[r]/.solgauge)//Evaluate],H1->Function[r,(H1[r]/.solgauge)//Evaluate]};

eq1=Collect[eqttpol/.reptrace/.repgauge/.repG/.repFLM,listfl1,FullSimplify];

eq2=Collect[eqtrpol/.reptrace/.repgauge/.repG/.repFLM,listfl1,FullSimplify];

eq3=Collect[eqrrpol/.reptrace/.repgauge/.repG/.repFLM,listfl1,FullSimplify];

eq4=Collect[eqrthpol/.reptrace/.repgauge/.repG/.repFLM,listfl1,FullSimplify];

repH2l1=Simplify[Solve[{eq2,eq4},{H2[r],H2'[r]}]]//First;

listl1rT={η1[rT[r]],η1'[rT[r]],η1''[rT[r]],K[rT[r]],K'[rT[r]],K''[rT[r]]};
listl1rT2={ψ1[rT[r]],ψ1'[rT[r]],ψ1''[rT[r]],ψ2[rT[r]],ψ2'[rT[r]],ψ2''[rT[r]]};
eqa=Collect[eq1/.repH2l1/.{η1->Function[r,r/(2Mass-r)ψ1[rT[r]]],K->Function[r,ψ2[rT[r]]/r]}/.rTdef/.backg,listl1rT2,Simplify];
eqb=Collect[eq3/.repH2l1/.{η1->Function[r,r/(2Mass-r)ψ1[rT[r]]],K->Function[r,ψ2[rT[r]]/r]}/.rTdef/.backg,listl1rT2,Simplify];
eqor1=Collect[eq1/.repH2l1/.rTdef/.backg,listfl1,Simplify];
eqor2=Collect[eq3/.repH2l1/.rTdef/.backg,listfl1,Simplify];

eqsor=Solve[{eqor1,eqor2},{K''[r],η1''[r]}]//First;
eql1η=Collect[η1''[r]/.eqsor,listfl1,Simplify];
eql1K=Collect[K''[r]/.eqsor,listfl1,Simplify];

Export["./eql1eta.m",eql1η];
Export["./eql1K.m",eql1K];

eqs=Solve[{eqa,eqb},{ψ1''[rT[r]],ψ2''[rT[r]]}]//First;
eql11=Collect[ψ1''[rT[r]]/.eqs,listl1rT2,Simplify];
eql12=Collect[ψ2''[rT[r]]/.eqs,listl1rT2,Simplify];
Export["./eql1_psi1.m",eql11];
Export["./eql1_psi2.m",eql12];

To obtain the formula for (dE/(dω)) in terms of the functions K and η1:

repsimp={Mass->0,ll->2,GLM[r]->0,ALM0[r]->0,GLM'[r]->0,ALM0'[r]->0,GLM''[r]->0,ALM0''[r]->0};
Clear[H0,H1,H2,K,η0,η1];
eqH0=H0[r]==(H0[r]/.reptrace/.repsimp);
eqH1=H1[r]==(H1[r]/.repgauge/.repsimp);
eqH2=H2[r]==(H2[r]/.repH2l1/.repsimp);
eqη0=η0[r]==(η0[r]/.repgauge/.repsimp);
K[r_]:=1/r*Kω*Exp[I*Sqrt[ω^2-μ^2]r]
H0[r_]:=1/r*H0ω*Exp[I*Sqrt[ω^2-μ^2]r]
H1[r_]:=1/r*H1ω*Exp[I*Sqrt[ω^2-μ^2]r]
H2[r_]:=1/r*H2ω*Exp[I*Sqrt[ω^2-μ^2]r]
η0[r_]:=η0ω*Exp[I*Sqrt[ω^2-μ^2]r]
η1[r_]:=η1ω*Exp[I*Sqrt[ω^2-μ^2]r]
repall=Solve[{eqH0,eqH1,eqH2,eqη0},{H0ω,H1ω,H2ω,η0ω}]//First;
Simplify[Series[H0ω/.repall,{r,∞,0}]]
Simplify[Series[H1ω/.repall,{r,∞,0}]]
Simplify[Series[H2ω/.repall,{r,∞,0}]]
Simplify[Series[η0ω/.repall,{r,∞,0}]]
Clear[H0,H1,H2,K,η0,η1];

Section 4 - Axial Equation

4.0 - Initialization

Clear[m,ll];
Clear[ALM,ALM1,ALM0,BLM,BLM0,GLM,FLM,QLM,QLM0,DLM];
DefScalarFunction/@{ALM1,ALM0,ALM,BLM,BLM0,FLM,GLM,QLM,QLM0,DLM};

Clear[EFElhs,EFErhs,eqgaugelhs,eqgaugerhs,metricTrace,SETrace];
{EFElhs,EFErhs,eqgaugelhs,eqgaugerhs,metricTrace,SETrace}=Import["./EFE_dRGT.m"];

EFEaxiallhs=Coefficient[EFElhs/.pol->0,ϵ,1]//Together;
EFEaxialrhs=Coefficient[EFErhs/.pol->0,ϵ,1]//Together;

eqttax=Simplify[EFEaxiallhs[[1,1]]==Coefficient[EFErhs[{0,-swc},{0,-swc}]/.pol->0,ϵ,1]];
eqtrax=Simplify[EFEaxiallhs[[1,2]]==Coefficient[EFErhs[{0,-swc},{1,-swc}]/.pol->0,ϵ,1]];
eqtthax=Simplify[EFEaxiallhs[[1,3]]==Coefficient[EFErhs[{0,-swc},{2,-swc}]/.pol->0,ϵ,1]];
eqtphax=Simplify[EFEaxiallhs[[1,4]]==Coefficient[EFErhs[{0,-swc},{3,-swc}]/.pol->0,ϵ,1]];
eqrrax=Simplify[EFEaxiallhs[[2,2]]==Coefficient[EFErhs[{1,-swc},{1,-swc}]/.pol->0,ϵ,1]];
eqrthax=Simplify[EFEaxiallhs[[2,3]]==Coefficient[EFErhs[{1,-swc},{2,-swc}]/.pol->0,ϵ,1]];
eqrphax=Simplify[EFEaxiallhs[[2,4]]==Coefficient[EFErhs[{1,-swc},{3,-swc}]/.pol->0,ϵ,1]];
eqththax=Simplify[EFEaxiallhs[[3,3]]==Coefficient[EFErhs[{2,-swc},{2,-swc}]/.pol->0,ϵ,1]];
eqthphax=Simplify[EFEaxiallhs[[3,4]]==Coefficient[EFErhs[{2,-swc},{3,-swc}]/.pol->0,ϵ,1]];
eqphphax=Simplify[EFEaxiallhs[[4,4]]==Coefficient[EFErhs[{3,-swc},{3,-swc}]/.pol->0,ϵ,1]];

eqaxtrace=Coefficient[Simplify[(metricTrace==-16π/(3μ^2)*SETrace)/.pol->0],ϵ,1]//Together;
eqaxgauget=Simplify[Coefficient[(gaugelhs[{0,-swc}]==gaugerhs[[0]][{0,-swc}])/.pol->0,ϵ,1]/.r[]->r/.t[]->t/.θ[]->θ/.ϕ[]->ϕ/.backg];
eqaxgauger=Simplify[Coefficient[(gaugelhs[{1,-swc}]==gaugerhs[{1,-swc}])/.pol->0,ϵ,1]/.r[]->r/.t[]->t/.θ[]->θ/.ϕ[]->ϕ/.backg];
eqaxgaugeth=Simplify[Coefficient[(gaugelhs[{2,-swc}]==gaugerhs[{2,-swc}])/.pol->0,ϵ,1]/.r[]->r/.t[]->t/.θ[]->θ/.ϕ[]->ϕ/.backg];
eqaxgaugeph=Simplify[Coefficient[(gaugelhs[{3,-swc}]==gaugerhs[{3,-swc}])/.pol->0,ϵ,1]/.r[]->r/.t[]->t/.θ[]->θ/.ϕ[]->ϕ/.backg];

eqsaxial={{eqttax,eqtrax,eqtthax,eqtphax,eqrrax,eqrthax,eqrphax,eqththax,eqthphax,eqphphax},{eqaxtrace,eqaxgauget,eqaxgauger,eqaxgaugeth,eqaxgaugeph}};
Export["./ax_dRGT.m",eqsaxial];

4.1 - l=1

Clear[eqttax,eqtrax,eqtthax,eqtphax,eqrrax,eqrthax,eqrphax,eqththax,eqthphax,eqphphax,eqaxtrace,eqaxgauget,eqaxgauger,eqaxgaugeth,eqaxgaugeph];
{{eqttax,eqtrax,eqtthax,eqtphax,eqrrax,eqrthax,eqrphax,eqththax,eqthphax,eqphphax},{eqaxtrace,eqaxgauget,eqaxgauger,eqaxgaugeth,eqaxgaugeph}}=Import["./ax_dRGT.m"];

ll=2;

sh0axl1=Simplify[Solve[eqaxgaugeph,h0[r]]//First//First];
reph0axl1=h0->Function[r,(h0[r]/.sh0axl1)//Evaluate];
eqh1axl1=Collect[eqrthax/.reph0axl1,{h1[r],h1'[r],h1''[r]},Simplify];
eqh1axl1rt=Collect[eqh1axl1/.h1->Function[r,1/(1-2Mass/r)Q[rT[r]]]/.rTdef,{Q[rT[r]],Q'[rT[r]],Q''[rT[r]]},Simplify];
eqaxl1=Collect[eqh1axl1rt/(Coefficient[eqh1axl1rt,Q''[rT[r]]]//First),{Q[rT[r]],Q'[rT[r]],Q''[rT[r]]},Simplify];

Export["./axial_l1_dRGT.m",eqaxl1];

4.2 - l≥2

Clear[h0,h1,h2,Q,Z]
Clear[m,ll];
Clear[eqttax,eqtrax,eqtthax,eqtphax,eqrrax,eqrthax,eqrphax,eqththax,eqthphax,eqphphax,eqaxtrace,eqaxgauget,eqaxgauger,eqaxgaugeth,eqaxgaugeph];
{{eqttax,eqtrax,eqtthax,eqtphax,eqrrax,eqrthax,eqrphax,eqththax,eqthphax,eqphphax},{eqaxtrace,eqaxgauget,eqaxgauger,eqaxgaugeth,eqaxgaugeph}}=Import["./ax_dRGT.m"];

sh0=Simplify[Solve[eqaxgaugeph,h0[r]]]//First//First;
reph0=h0->Function[r,(h0[r]/.sh0)//Evaluate];
{ddh1,ddh2}=Simplify[Solve[{eqrthax/.reph0,Coefficient[eqththax,Cot[θ]PL[θ]]/.reph0},{h1''[r],h2''[r]}]//First];
RWE1=Simplify[(h1''[r]-(h1''[r]/.ddh1)==0)];
RWE2=Simplify[(h2''[r]-(h2''[r]/.ddh2)==0)];
Collect[RWE1,{h1[r],h1'[r],h1''[r],h2[r],h2'[r],h2''[r]}];
Collect[RWE2,{h1[r],h1'[r],h1''[r],h2[r],h2'[r],h2''[r]}];
RWE1new=RWE1/.{h1->Function[r,f1[r]Q[r]+f2[r]Z[r]],h2->Function[r,f3[r]Q[r]+f4[r]Z[r]]};
RWE2new=RWE2/.{h1->Function[r,f1[r]Q[r]+f2[r]Z[r]],h2->Function[r,f3[r]Q[r]+f4[r]Z[r]]};
solf23=Simplify[Solve[{Coefficient[RWE1new,Z''[r]],Coefficient[RWE2new,Q''[r]]},{f2[r],f3[r]}]]//First;
solf2=Simplify[f2[r]/.solf23[[1]]];
solf3=Simplify[f3[r]/.solf23[[2]]];
RWE1aux1=RWE1new/.f2->Function[r,solf2//Evaluate]/.f3->Function[r,solf3//Evaluate];
RWE2aux1=RWE2new/.f2->Function[r,solf2//Evaluate]/.f3->Function[r,solf3//Evaluate];
RWE1aux2=Simplify[RWE1aux1/.{Q-> Function[r,Q[rT[r]]],Z-> Function[r,Z[rT[r]]]}/.rTdef];
RWE2aux2=Simplify[RWE2aux1/.{Q-> Function[r,Q[rT[r]]],Z-> Function[r,Z[rT[r]]]}/.rTdef];
repf1=Simplify[DSolve[Coefficient[RWE1aux2,Q'[rT[r]]],f1[r],r]/.C[1]->1]//First//First;
repf4=Simplify[DSolve[Coefficient[RWE2aux2,Z'[rT[r]]],f4[r],r]/.C[1]->1]//First//First;
solf1=f1[r]/.(repf1);
solf4=f4[r]/.(repf4);
RWE1new2=Simplify[RWE1aux2/.f1->Function[r,solf1//Evaluate]/.f4->Function[r,solf4//Evaluate]];
RWE2new2=Simplify[RWE2aux2/.f1->Function[r,solf1//Evaluate]]/.f4->Function[r,solf4//Evaluate];
RWE1f=Collect[(-RWE1new2*(2Mass-r)^2/r^4)/.ll->2(λ+1),{Q''[rT[r]],Q[rT[r]],Z[rT[r]],μ^2,ω^2},Simplify];
RWE2f=Collect[(RWE2new2*(2Mass-r)/r^3)/.ll->2(λ+1),{Q''[rT[r]],Q[rT[r]],Z[rT[r]],μ^2,ω^2},Simplify];

RWE={RWE1f,RWE2f};
Export["./axial_l2_dRGT.m",RWE];

To obtain the formula for (dE/(dω)) in terms of the master functions Q and Z:

Q[r_]:=Qω Exp[I Sqrt[ω^2-μ^2] r]
Z[r_]:=Zω Exp[I Sqrt[ω^2-μ^2]r]
h1[r_]:=Q[r]
h2[r_]:=Z[r]
Series[h0[r]/.reph0/.Mass->0,{r,∞,1}]
Clear[Q,Z,h1,h2];

Section 5 - Source terms

As a function of t:

$Assumptions={r[]∈Reals,θ[]∈Reals,ϕ[]∈Reals,ω∈Reals,t[]∈Reals,PL[_]∈Reals,m∈Reals,PL'[_]∈Reals,ll∈Reals,Θ[_]∈Reals,t∈Reals,r∈Reals,θ∈Reals,ϕ∈Reals,ωk∈Reals};

metricmink=DiagonalMatrix[{-1,1,r[]^2,r[]^2Sin[θ[]]^2}];
invmink=Inverse[metricmink];
vecschwarz=D[{t,R[t],Θ[t],Φ[t]},t];

normalm=1;
normalm0=1;
normalm1=-1;
normblm=1;
normblm0=-1;
normclm=1;
normclm0=-1;
normglm=1;(*VERIFICAR A NORMA DESTE*)
normdlm=-1;
normflm=1;

ALMt[r_]:=Refine[m0*γ*Integrate[Sum[Conjugate[alm[[i,j]]]invmink[[i,k]]metricBG[[k,l]]vecschwarz[[l]]invmink[[j,mm]]metricBG[[mm,n]]vecschwarz[[n]]/normalm,{i,1,4},{j,1,4},{k,1,4},{l,1,4},{mm,1,4},{n,1,4}]DiracDelta[θ[]-Θ[t]]DiracDelta[ϕ[]-Φ[t]],{θ[],-∞,∞},{ϕ[],-∞,∞}]*DiracDelta[r-R[t]]/r^2]
ALM1t[r_]:=Refine[m0*γ*Integrate[Sum[Conjugate[alm1[[i,j]]]invmink[[i,k]]metricBG[[k,l]]vecschwarz[[l]]invmink[[j,mm]]metricBG[[mm,n]]vecschwarz[[n]]/normalm1,{i,1,4},{j,1,4},{k,1,4},{l,1,4},{mm,1,4},{n,1,4}]DiracDelta[ϕ[]-Φ[t]]DiracDelta[θ[]-Θ[t]],{θ[],-∞,∞},{ϕ[],-∞,∞}]*DiracDelta[r-R[t]]/r^2]
ALM0t[r_]:=Refine[m0*γ*Integrate[Sum[Conjugate[alm0[[i,j]]]invmink[[i,k]]metricBG[[k,l]]vecschwarz[[l]]invmink[[j,mm]]metricBG[[mm,n]]vecschwarz[[n]],{i,1,4},{j,1,4},{k,1,4},{l,1,4},{mm,1,4},{n,1,4}]DiracDelta[ϕ[]-Φ[t]]DiracDelta[θ[]-Θ[t]],{θ[],-∞,∞},{ϕ[],-∞,∞}]*DiracDelta[r-R[t]]/r^2/normalm0]
BLM0t[r_]:=Refine[m0*γ*Integrate[Sum[Conjugate[blm0[[i,j]]]invmink[[i,k]]metricBG[[k,l]]vecschwarz[[l]]invmink[[j,mm]]metricBG[[mm,n]]vecschwarz[[n]]/normblm0,{i,1,4},{j,1,4},{k,1,4},{l,1,4},{mm,1,4},{n,1,4}]DiracDelta[ϕ[]-Φ[t]]DiracDelta[θ[]-Θ[t]],{θ[],-∞,∞},{ϕ[],-∞,∞}]*DiracDelta[r-R[t]]/r^2]
BLMt[r_]:=Refine[m0*γ*Integrate[Sum[Conjugate[blm[[i,j]]]invmink[[i,k]]metricBG[[k,l]]vecschwarz[[l]]invmink[[j,mm]]metricBG[[mm,n]]vecschwarz[[n]]/normblm,{i,1,4},{j,1,4},{k,1,4},{l,1,4},{mm,1,4},{n,1,4}]DiracDelta[ϕ[]-Φ[t]]DiracDelta[θ[]-Θ[t]],{θ[],-∞,∞},{ϕ[],-∞,∞}]*DiracDelta[r-R[t]]/r^2]
GLMt[r_]:=Refine[(m0*γ*Integrate[Sum[Conjugate[glm[[i,j]]]invmink[[i,k]]metricBG[[k,l]]vecschwarz[[l]]invmink[[j,mm]]metricBG[[mm,n]]vecschwarz[[n]]/normglm,{i,1,4},{j,1,4},{k,1,4},{l,1,4},{mm,1,4},{n,1,4}]DiracDelta[ϕ[]-Φ[t]]DiracDelta[θ[]-Θ[t]],{θ[],-∞,∞},{ϕ[],-∞,∞}]*DiracDelta[r-R[t]]/r^2)]
FLMt[r_]:=Refine[(m0*γ*Integrate[Sum[Conjugate[flm[[i,j]]]invmink[[i,k]]metricBG[[k,l]]vecschwarz[[l]]invmink[[j,mm]]metricBG[[mm,n]]vecschwarz[[n]]/normflm,{i,1,4},{j,1,4},{k,1,4},{l,1,4},{mm,1,4},{n,1,4}]DiracDelta[ϕ[]-Φ[t]]DiracDelta[θ[]-Θ[t]],{θ[],-∞,∞},{ϕ[],-∞,∞}]*DiracDelta[r-R[t]]/r^2)]
DLMt[r_]:=Refine[(m0*γ*Integrate[Sum[Conjugate[dlm[[i,j]]]invmink[[i,k]]metricBG[[k,l]]vecschwarz[[l]]invmink[[j,mm]]metricBG[[mm,n]]vecschwarz[[n]]/normdlm,{i,1,4},{j,1,4},{k,1,4},{l,1,4},{mm,1,4},{n,1,4}]DiracDelta[ϕ[]-Φ[t]]DiracDelta[θ[]-Θ[t]],{θ[],-∞,∞},{ϕ[],-∞,∞}]*DiracDelta[r-R[t]]/r^2)]
QLMt[r_]:=Refine[(m0*γ*Integrate[Sum[Conjugate[clm[[i,j]]]invmink[[i,k]]metricBG[[k,l]]vecschwarz[[l]]invmink[[j,mm]]metricBG[[mm,n]]vecschwarz[[n]]/normclm,{i,1,4},{j,1,4},{k,1,4},{l,1,4},{mm,1,4},{n,1,4}]DiracDelta[ϕ[]-Φ[t]]DiracDelta[θ[]-Θ[t]],{θ[],-∞,∞},{ϕ[],-∞,∞}]*DiracDelta[r-R[t]]/r^2)]
QLM0t[r_]:=Refine[(m0*γ*Integrate[Sum[Conjugate[clm0[[i,j]]]invmink[[i,k]]metricBG[[k,l]]vecschwarz[[l]]invmink[[j,mm]]metricBG[[mm,n]]vecschwarz[[n]]/normclm0,{i,1,4},{j,1,4},{k,1,4},{l,1,4},{mm,1,4},{n,1,4}]DiracDelta[ϕ[]-Φ[t]]DiracDelta[θ[]-Θ[t]],{θ[],-∞,∞},{ϕ[],-∞,∞}]*DiracDelta[r-R[t]]/r^2)]

For radial infall

Clear[R,Θ,Φ,γ,ll]
Clear[MyFourier,ALM1,ALM0,ALM,BLM,BLM0,GLM,FLM,DLM,QLM,QLM0];
R'[t_]:=-(1-2Mass/R[t])
Θ[t_]:=0;
Φ[t_]:=0;
T[r_]:=-rT[r];
γ=ee/(1-2Mass/r[]);
MyFourier[f_]:=Assuming[r∈Reals,Return[Simplify[(Integrate[f/.r[]->r/.t[]->t,{r,-∞,∞}]/(R'[t]Sqrt[2π]))/.R[t]->r/.t->T[r]]]]
ALM1[r_]=MyFourier[ALM1t[r]];
ALM[r_]=MyFourier[ALMt[r]];
ALM0[r_]=MyFourier[ALM0t[r]];
BLM[r_]=MyFourier[BLMt[r]];
BLM0[r_]=MyFourier[BLM0t[r]];
GLM[r_]=MyFourier[GLMt[r]];
FLM[r_]=MyFourier[FLMt[r]];
DLM[r_]=MyFourier[DLMt[r]];
QLM[r_]=MyFourier[QLMt[r]];
QLM0[r_]=MyFourier[QLM0t[r]];
sourcefunctionsradial={ALM0[r],ALM[r],ALM1[r],BLM0[r],BLM[r],GLM[r],FLM[r],QLM0[r],QLM[r],DLM[r]}
Export["./source_radial.m",sourcefunctionsradial];

For circular orbit

Clear[R,Θ,Φ,T,γ,ll];
Clear[MyFourier,ALM1,ALM0,ALM,BLM,BLM0,GLM,FLM,DLM,QLM,QLM0];
R[t_]:=Rc;
Θ[t_]:=π/2;
Φ[t_]:=Sqrt[Mass/r[]^3]*t;
(*ωk=Sqrt[Mass/Rc^3];*)
γ=1/Sqrt[1-3Mass/r[]];
MyFourier[f_]:=Assuming[t∈Reals,Return[Simplify[FourierTransform[(f/.r[]->r/.t[]->t)*E^(-I ω t),t,ω]]]]
ALM1[r_]=Simplify[MyFourier[ALM1t[r]]]/.DiracDelta[-m Sqrt[Mass/r^3]+ω]->δω;
ALM[r_]=Simplify[MyFourier[ALMt[r]]]/.DiracDelta[-m Sqrt[Mass/r^3]+ω]->δω;
ALM0[r_]=Simplify[MyFourier[ALM0t[r]]]/.DiracDelta[-m Sqrt[Mass/r^3]+ω]->δω;
BLM[r_]=Simplify[MyFourier[BLMt[r]]]/.DiracDelta[-m Sqrt[Mass/r^3]+ω]->δω;
BLM0[r_]=Simplify[MyFourier[BLM0t[r]]]/.DiracDelta[-m Sqrt[Mass/r^3]+ω]->δω;
GLM[r_]=Simplify[MyFourier[GLMt[r]]]/.DiracDelta[-m Sqrt[Mass/r^3]+ω]->δω;
FLM[r_]=Simplify[MyFourier[FLMt[r]]]/.DiracDelta[-m Sqrt[Mass/r^3]+ω]->δω;
DLM[r_]=Simplify[MyFourier[DLMt[r]]]/.DiracDelta[-m Sqrt[Mass/r^3]+ω]->δω;
QLM[r_]=Simplify[MyFourier[QLMt[r]]]/.DiracDelta[-m Sqrt[Mass/r^3]+ω]->δω;
QLM0[r_]=Simplify[MyFourier[QLM0t[r]]]/.DiracDelta[-m Sqrt[Mass/r^3]+ω]->δω;
sourcefunctionscircular={ALM0[r],ALM[r],ALM1[r],BLM0[r],BLM[r],GLM[r],FLM[r],QLM0[r],QLM[r],DLM[r]};
Export["./source_circular.m",sourcefunctionscircular];