Notation

M0,qPP: mass and charge of point particle
gamma= dTdtau of point particle
M,Q: mass and charge of black hole

Geodesics :

XX={t[tau],r[tau],theta,phi};

(*ND is the number of spacetime dimensions in the problem. Has to be in accord with metric; //g is covariant metric*)
ND=4;

(*Equation for Spherical Harmonics, we set m=0 since one can always rotate axis*)

m=0;
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'[θ],
Y''[theta]->-Cos[theta]/Sin[theta]* Y'[theta]-l*(l+1)*Y[theta]+m^2/Sin[theta]^2*Y[theta],Derivative[3][Y][theta]->Cot[theta] (l+(l^2)-3 (m^2) ((Csc[theta])^2)) Y[theta]+(-l (1+l)+((Cot[theta])^2)+(1+(m^2)) ((Csc[theta])^2)) Y'[theta]};

fe=1-2*M/r[tau];

g0=DiagonalMatrix[{-fe,fe^(-1),r[tau]^2,r[tau]^2*Sin[theta]^2}];

AMax0={Q/r[tau]*0,0,0,0};

A=Inverse[g0];

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

(*COVARIANT MAXWELL TENSOR:*)FMax0=Simplify[Table[D[AMax0[[nu]],XX[[mu]]]-D[AMax0[[mu]],XX[[nu]]],{mu,1,ND},{nu,1,ND}]];

(*(Subscript[∇, c])(X^( a))=(∂_c)(X^( a))+(Subscript[(Γ^( a)), bc])(X^( b)); the EOM is (4) in http://arxiv.org/pdf/1103.1807.pdf*)Geo=Table[Sum[D[XX[[a]],tau]*D[XX[[beta]],XX[[a]]],{a,1,ND}]-Sum[D[XX[[a]],tau]*Γ[[beta,b,a]]*D[XX[[b]],tau],{a,1,ND},{b,1,ND}]-Sum[qPP/M0*FMax0[[deta,gamma]]*A[[deta,beta]]*D[XX[[gamma]],tau],{gamma,1,ND},{deta,1,ND}],{beta,1,ND}];Simplify[Geo]

Solve[{Geo[[1]]==0,Geo[[2]]==0},{t'[tau],r'[tau]}]

Equations (note: 1/Sqrt(2pi) for each FT)

(*ND is the number of spacetime dimensions in the problem. Has to be in accord with metric; //g is covariant metric*)
ND=4;

fe=1-2*M/r;

g0=DiagonalMatrix[{-fe,fe^(-1),r^2,r^2*Sin[theta]^2}];

(*this relation comes down from Einstein equations, I am using it already:*)
H2[r_]:=H0[r];

g=g0+eps*{{fe*H0[r]*Y[theta],H1[r]*Y[theta],0,0},{H1[r]*Y[theta],1/fe*H2[r]*Y[theta],0,0},{0,0,r^2*K[r]*Y[theta],0},{0,0,0,r^2*K[r]*Sin[theta]^2*Y[theta]}}*Exp[-I*w*t];

(*A é a métrica contravariante*)A=Inverse[g];
X={t,r,theta,phi};

(*Γ are the metric connections (Christoffel symbols) ,the first index is up, the next two are down*)Γ=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,ND}],{c,1,ND},{a,1,ND},{b,1,ND}];

(*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,ND}]-Sum[Γ[[e,b,c]]Γ[[f,e,d]],{e,1,ND}]),{f,1,ND}],{a,1,ND},{b,1,ND},{c,1,ND},{d,1,ND}];

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

(*Ricci Scalar:*)Riccie0=Sum[A[[a,b]]Riccit[[a,b]],{a,1,ND},{b,1,ND}];Riccie=Normal[Series[Riccie0,{eps,0,1}]];

MatrixForm[g]

(*geodesics for infalling particles see http://arxiv.org/pdf/1103.1807.pdf or http://arxiv.org/pdf/1011.5399v1.pdf*)

gamma=dTdtau[r];

dRdt=dRdtime[r];

Einstein

(*Y00 is the spherical harmonics evaluated at theta=phi=0, since the particle is falling along the z-axis*)

Y00=((Sqrt[(((2*l+1))/π)])/2);
Tmunuco=eps*{{M0/Sqrt[2*Pi]*gamma*fe^2/r^2*Y[theta],-M0/Sqrt[2*Pi]*gamma/r^2*dRdt*Y[theta],0,0},{-M0/Sqrt[2*Pi]*gamma/r^2*dRdt*Y[theta],M0/Sqrt[2*Pi]*gamma/fe^2*1/r^2*(dRdt)^2*Y[theta],0,0},{0,0,0,0},{0,0,0,0}}*Y00*Exp[I*w*T[r]]*Exp[-I*w*t]*1/(dRdt);(*note to self: the last 1/dRdt is to do the Fourier integral*)

(*Einstein's covariant tensor:*)G=Table[Riccit[[a,b]]-1/2*g[[a,b]]*Riccie,{a,1,ND},{b,1,ND}];

FieldSeries=Normal[Series[Einstot=G-8*Pi*(Tmunuco),{eps,0,1}]];

FullSimplify[FieldSeries[[1,1]]//.par]

Eq33=FullSimplify[FieldSeries[[3,3]]//.par];;

Eq23=FullSimplify[FieldSeries[[2,3]]//.par];;

Eq22=FullSimplify[FieldSeries[[2,2]]//.par];;

Eq13=FullSimplify[FieldSeries[[1,3]]//.par];

Eq12=FullSimplify[FieldSeries[[1,2]]//.par];

Eq11=FullSimplify[FieldSeries[[1,1]]//.par];

finding reduced system

solp=Simplify[Solve[{Eq12==0,Eq13==0,Eq22==0},{K'[r],H0'[r],H1'[r]}]][[1]]

H0sol=Simplify[Solve[FullSimplify[Eq23/.solp]==0,H0[r]]][[1]];H0sol1=Simplify[D[H0sol,r]/.solp];

(*These equations are the analogue of Zerilli 33 and 34 in the RN work*)
{Eqb,Eqc}=Simplify[{(-(1/(2 (r^2) (-2 M+r) (6 M+(-2+l+(l^2)) r)))(E^(-I t w)) eps Y[theta])^(-1)Eq12,((E^(-I t w)) eps Y'[theta])^(-1)Eq13}/.H0sol/.H0sol1]

(*we need an isolated "-K'" to identify alpha and beta*)
Collect[Eqb2=-Eqb*(2 I (2 M-r) (r^2) (6 M+(-2+l+(l^2)) r) w )^(-1),K'[r]]

xb=Simplify[Eqb2/.{H1[r]->H1[r]*0,K[r]->1,K'[r]->0,fLM'[r]->0,M0->0,fLM[r]->0}];

xb2=Series[xb,{w,0,4}]

alpha0=SeriesCoefficient[xb2,0]

alpha2=SeriesCoefficient[xb2,2]

xc=Simplify[Eqb2/.{H1[r]->1,K[r]->0,K'[r]->0,fLM'[r]->0,M0->0,fLM[r]->0}];

xx2=Series[xc,{w,0,4}]

beta0=SeriesCoefficient[xx2,-1]

beta2=SeriesCoefficient[xx2,1]

(*Here, purpose is to isolated -H1'/w*)
Collect[Eqc2=Eqc*r*2/(-w*(-2 M+r)),H1'[r]]

xd=Simplify[Eqc2/.{H1[r]->H1[r]*0,H1'[r]->0,K[r]->1,K'[r]->0,fLM'[r]->0,M0->0,fLM[r]->0}]

xy2=Series[xd,{w,0,4}]

gamma0=SeriesCoefficient[xy2,0];
gamma2=SeriesCoefficient[xy2,2]

xd2=Simplify[Eqc2/.{H1[r]->1,H1'[r]->0,K[r]->0,K'[r]->0,fLM'[r]->0,M0->0,fLM[r]->0}]

xy3=Series[xd2,{w,0,4}]

delta0=SeriesCoefficient[xy3,-1]
delta2=SeriesCoefficient[xy3,1]

Zerilli transformations

(*This checks that the transformations satisfy Zerilli's condition 41*)

alpha=alpha0+w^2*alpha2;
beta=beta0+w^2*beta2;
gamma=gamma0+w^2*gamma2;
delta=delta0+w^2*delta2;

Q=0;

h[r_]:=-I (1+(((Q^2)-M r)/((Q^2)-2 M r+(r^2)))+((8 (Q^2)-6 M r)/(-4 (Q^2)+r (6 M+(-2+l+(l^2)) r))));

ffe[r_]:=1-2*M/r+Q^2/r^2;

f[r_]=((16 (Q^4)-4 (Q^2) r (11 M+(-4+l+(l^2)) r)+(r^2) (24 (M^2)+6 (-2+l+(l^2)) M r+(-1+l) l (1+l) (2+l) (r^2)))/(2 (r^3) (-4 (Q^2)+r (6 M+(-2+l+(l^2)) r))));

k[r_]:=-((I r)/(ffe[r]));

V=(((Q^2)+r (-2 M+r)) (-32 (Q^6)+24 (Q^4) r (6 M+(-2+l+(l^2)) r)+8 (Q^2) (r^2) (-27 (M^2)-2 (-11+4 l (1+l)) M r+3 (-2+l+(l^2)) (r^2))+(r^3) (72 (M^3)+36 (-2+l+(l^2)) (M^2) r+6 (((-2+l+(l^2)))^2) M (r^2)+l (1+l) (((-2+l+(l^2)))^2) (r^3))))/((r^6) (((-4 (Q^2)+r (6 M+(-2+l+(l^2)) r)))^2));

A={{alpha,beta},{gamma,delta}};

F={{f[r],1},{h[r],k[r]}};

FullSimplify[ffe[r]*Inverse[F].(A.F-D[F,r])-{{0,1},{V-w^2,0}}]

Reduction to two second order ODEs

K[r_]:=f[r]*RLM[r]+LL[r];
H1[r_]:=w*h[r]*RLM[r]+w*k[r]*LL[r];

Eqc2=Simplify[Eqc]

Eqb2=Simplify[Eqb]

(* Let us use (36) and (37) in Zerilli *)

LLpsol=Simplify[Solve[Eqb2==0,LL'[r]]][[1]];

(*Eq 36*)Eqc3=Simplify[Eqc2/.LLpsol]

RLMpsol=Simplify[Solve[Eqc3==0,RLM'[r]]][[1]][[1]]

(*Eq37*)Eqb3=Simplify[Eqb2/.RLMpsol]

(*let us isolate "-(1-2*M/r)*R'"*)Collect[Eqc4=-(1-2*M/r)(-((I (2 M-r) (6 M+(-2+l+(l^2)) r)*w)/(4 r (-2 M+r))))^(-1)Eqc3,RLM'[r]]

S1hat=Eqc4/.{RLM'[r]->0,LL[r]->0}

(*let us isolate "-(1-2*M/r)*L'"*)Collect[Eqb4=-(1-2*M/r)(2 I (r^2) (12 (M^2)+2 (-5+l+(l^2)) M r-(-2+l+(l^2)) (r^2)) w )^(-1)Eqb3,LL'[r]]

S2hat=Simplify[Eqb4/.{LL'[r]->0,RLM[r]->0}]

Q=0;Vgrav=(((Q^2)+r (-2 M+r)) (-32 (Q^6)+24 (Q^4) r (6 M+(-2+l+(l^2)) r)+8 (Q^2) (r^2) (-27 (M^2)-2 (-11+4 l (1+l)) M r+3 (-2+l+(l^2)) (r^2))+(r^3) (72 (M^3)+36 (-2+l+(l^2)) (M^2) r+6 (((-2+l+(l^2)))^2) M (r^2)+l (1+l) (((-2+l+(l^2)))^2) (r^3))))/((r^6) (((-4 (Q^2)+r (6 M+(-2+l+(l^2)) r)))^2));

TotSource=FullSimplify[(1-2*M/r)*D[S1hat,r]+S2hat]

ff=1-2*M/r;
ff1=D[ff,r];

EQgrav0=ff^2*D[ygrav[r],{r,2}]+ff*ff1*D[ygrav[r],{r,1}]+(w^2-Vgrav)*ygrav[r]-TotSource

Final equations with motion plugged in

fe=1-2*M/r;

f=fe;f1=D[fe,r];

gamma=En/fe;

gamma1=Simplify[D[gamma,r]];

dRdt=-fe*Sqrt[En^2-fe]/En;

TotSource0=((4 (E^(I w T[r])) (Sqrt[2+4 l]) M0 (I r (((-2 M+r))^2) (6 M+(-2+l+(l^2)) r) dTdtau'[r]+dTdtau[r] (-2 I (2 M-r) (6 (M^2)+3 (-2+l+(l^2)) M r-(-2+l+(l^2)) (r^2))+r (6 M+(-2+l+(l^2)) r) w ((r^2) dRdtime[r]-(((-2 M+r))^2) T'[r]))))/((r^3) (((6 M+(-2+l+(l^2)) r))^2) w));

TotSource=Simplify[TotSource0/.{dTdtau[r]->gamma,dRdtime[r]->dRdt,dTdtau'[r]->gamma1,T'[r]->1/dRdt}]

(*When energy is very large, we have a photon coliding with a BH, the source term is*)Sourcelight=Normal[Series[TotSource,{En,Infinity,0}]]

(*Check that it gives the correctultra-relativistic limit, https://arxiv.org/pdf/gr-qc/0202019.pdf*)
sigma=Expand[(l-1)*(l+2)/2];
Simplify[Sourcelight+f*4*I*En*M0*(((E^(I w T[r]))  (Sqrt[2+4 l]) sigma)/((((3 M+sigma r))^2) w))]

References

[1] V. Cardoso and J. P. S. Lemos, "Gravitational radiation from collisions at the speed of light: A Massless particle falling into a Schwarzschild black hole". PPhys.Lett.B538 (2002) 1-5.

[2]E. Berti, V. Cardoso and C. M. Will, "On gravitational-wave spectroscopy of massive black holes with the space interferometer LISA", Phys. Rev. D73:064030 (2006). [arXiv: gr-qc/0512160].

[3]E. Berti, V. Cardoso, T. Hinderer, M. Lemos, F. Pretorius, U. Sperhake and N. Yunes, "Semianalytical estimates of scattering thresholds and gravitational radiation in ultrarelativistic black hole encounters", Phys.Rev.D81 (2010) 104048.

[4]V. Cardoso, L. Gualtieri, C. Herdeiro, U. Sperhake, "Exploring New Physics Frontiers Through Numerical Relativity", Living Rev.Relativity 18 (2015) 1 
[5](The classical original work)M. Davis, R. Ruffini, W. H. Press and R. Price, "Gravitational radiation from a particle falling radially into a schwarzschild black hole", Phys.Rev.Lett.27 (1971) 1466-1469

[6]V. Cardoso, C. Macedo,P. Pani and V. Ferrari, "Black holes and gravitational waves in models of minicharged dark matter", JCAP 1605 (2016) no.05,054

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