Analytical template for gravitational-wave echoes from spinning exotic compact objects. 
Obtained by solving Teukolsky’s equations in the low-frequency limit.
Reference: Elisa Maggio, Adriano Testa, Swetha Bhagwat, Paolo Pani, “Analytical model for gravitational-wave echoes from spinning remnants”, arXiv:1907.03091

Quit[];

BH reflection coefficient at low frequency [Eq. (A12)]

k=ω-((m*a)/(2 M rp));
γ=rm/rp;
ξ=1+(Sqrt[1-(χ^2)]);
ζ=((I (2 ω M-m(Sqrt[γ])) (γ+1) ξ)/(γ-1));
RBHlf=(-8 M k  (E^((ζ(γ-1))/(γ+1)))((2 M k-I (γ-1))/((((γ-1))^2) ))(((((-M (γ-1) ξ)/L)))^(ζ(γ-1)))   (((16 (k^2) (M^2))/(((γ-1))^2))+1) Gamma[-2+ζ] Gamma[-1-ζ] (1800 I Gamma[-2-ζ]+(((ω M (γ-1) ξ))^5) Gamma[3-ζ]))/(Gamma[-2-ζ] Gamma[3-ζ] (1800 I Gamma[-2+ζ]+(((ω M (γ-1) ξ))^5) Gamma[3+ζ]));

Interpolating function from low to high frequencies [Eq. (23)]

kH=(1/2)(rp-rm)/((rp^2)+(a^2));

RBH=RBHlf*((Exp[((-2 π ωR)/kH)]+1)/(Exp[((2 π(Abs[ω]-ωR))/kH)]+1));

BH ringdown at infinity [Eq. (26)]

ωQNM=ωR+I*ωI;
α1=I*Exp[-I(ϕ+t0*ωQNM)];
α1p=α1/.{ϕ->ϕp};
α1c=α1/.{ϕ->ϕc};
α2p=-Conjugate[α1p];
α2c=-Conjugate[α1c];

ZBHp=((Exp[I ω t0])/(2 (Sqrt[2 π])))(((α1p*Ap+α1c*Ac)/(ω-ωQNM))+((α2p*Ap-α2c*Ac)/(ω+Conjugate[ωQNM])));

Transfer function [Eq. (30)]

x=1/(√(-(a^2)+(M^2)))(√(-(a^2)+(M^2)) r+M (M+√(-(a^2)+(M^2))) Log[((-M-√(-(a^2)+(M^2))+r)/L)]+M (-M+√(-(a^2)+(M^2))) Log[((-M+√(-(a^2)+(M^2))+r)/L)]);
x0=x//.{r->r0,r0-> rp(1+ϵ)};

KZBHm=((RBH*R*Exp[-2 I k x0])/(1-RBH*R*Exp[-2 I k x0]))ZBHp;

Final template [Eq. (14)]

Zp=ZBHp+KZBHm;

Example: evaluation for some parameters

param={L->LL*M,rp->M(1+(Sqrt[1-(χ^2)])),rm->M(1-(Sqrt[1-(χ^2)])),a->χ*M,ϵ-> (Sqrt[1-(χ^2)])((δ^2)/(4(rp^2)))};
sost={LL->1,M->1,m->2,χ->0,ωR->0.37367,ωI-> -0.08896,t0->0,Ap->1,Ac->0,ϕp->0,ϕc->0,R->1,δ->(10^(-5))M};

Zp//.param//.sost/.ω->0.1