The QNMs of Kerr black holes

Computing the QN frequencies

ar = 0.200;br = Sqrt[1 - 4*ar^2];rminus = (1 - br)/2;rplus = (1 + br)/2;s = -2;m = 2;l = 2;k1 = 1/2*Abs[m - s];k2 = 1/2*Abs[m + s];Alm = l*(l + 1) - s*(s + 1)  ;

winit = 2*(0.37 - 0.08 I);

\(ie = 0; 
    While[ie < 7, NITMAX = 130; c0 = 1 - s - I*w - 2*I/br*\((w/2 - ar*m)\); 
      c1 = \(-4\) + 2*I*w*\((2 + br)\) + 4*I/br*\((w/2 - ar*m)\); 
      c2 = s + 3 - 3*I*w - 2*I/br*\((w/2 - ar*m)\); 
      c3 = w^2*\((4 + 2*br - ar^2)\) - 2*ar*m*w - s - 1 + \((2 + br)\)*I*w - 
          Alm + \((4*w + 2*I)\)/br*\((w/2 - ar*m)\); 
      c4 = s + 1 - 
          2*w^2 - \((2*s + 3)\)*I*
            w - \((4*w + 2*I)\)/br*\((w/2 - ar*m)\); γ = 
        Function[n, 
          n\^2 + \((c2 - 3)\)*n + c4 - c2 + 2]; 
β = 
        Function[
          n, \(-2\)*n^2 + \((c1 + 2)\)*n + c3]; 
α = 
        Function[n, n\^2 + \((c0 + 1)\)*n + c0]; 
Leaver31[
          w_] := Module[{Rn}, 
For[{n = NITMAX; 
              Rn = \(-1.0\);}, 
n > 
              0, 
{
Rn = γ[n]\/\(β[n] - α[n]*Rn\); \(n--\);}
]; Rn]; 
      Leaver33[w_] := β[0]\/α[0] - Leaver31[w]; 
      wang = \(\(FindRoot[
              Leaver33[w] == 0, {w, winit}]\)[\([1]\)]\)[\([2]\)]; 
      witer[ie] = wang; γang = 
        Function[n, 2*ar*wang*\((n + k1 + k2 + s)\)]; 
βang = 
        Function[n, 
          n*\((n - 1)\) + 
            2*n*\((k1 + k2 + 1 - 2*ar*wang)\) - \((2*ar*
                  wang*\((2*k1 + s + 1)\) - \((k1 + k2)\)*\((k1 + k2 + 
                      1)\))\) - \((ar^2*wang^2 + s*\((s + 1)\) + 
                Sep)\)]; 
αang = 
        Function[n, \(-2\)*\((n + 1)\)*\((n + 2*k1 + 1)\)]; 
      Leaver31Ang[Sep_] := 
        Module[{Rn}, 
For[{n = NITMAX; 
              Rn = \(-1.0\);}, 
n > 
              0, 
{
Rn = γang[n]\/\(βang[n] - αang[n]*Rn\); \(n--\);}
]; Rn]; 
      Leaver33Ang[Sep_] := βang[0]\/αang[0] - Leaver31Ang[Sep]; 
      Alm = \(\(FindRoot[
              Leaver33Ang[Sep] == 0, {Sep, Alm}]\)[\([1]\)]\)[\([2]\)]; 
      ie = ie + 1];\)

Table[witer[n], {n, 0, 6}]

Angular wavefunctions

Sep = Alm;Nmax = 150;w = wang;sima = (w*rplus - ar*m)/br;ae[0] = 1;ae[1] = -βang[0]/αang[0];i = 1;While[i < Nmax, 
        ae[i + 1] = (-βang[i]*ae[i] - γang[i]*
                  ae[i - 1])/(αang[i]); i = i + 1];Sang = 
        Exp[ar*w*
                u]*(1 + u)^(1/2*Abs[m - s])*(1 - u)^(1/2*
                  Abs[m + s])*
            Sum[ae[i]*(1 + u)^(i), {i, 0, Nmax - 1}] /. u -> u;Normalization=1/Sqrt[NIntegrate[Sang*Conjugate[Sang],{u, -1, 1},MaxRecursion->12,PrecisionGoal -> 16,AccuracyGoal -> 16]];SangNormal = Normalization*Sang;

Plot[Re[SangNormal], {u, -0.99, 0.99}]

SangNormal /. {u -> Cos[1.3]}

Radial Wavefunctions

(*\(We\ now\ use\ Leaver' 
            s\ expansion\ to\ compute\ the\ radial\ wavefunctions\ themselves . \ The\ wavefunction\ is\ in - 
          going\ at\ the\ horizon, \ ψ == \(e\^\(iw\ \(r\_+\)\)\) \(\((\(r\_+\) - \(r\_-\))\)\^\(\(-1\) - s + 
                  iw + \(iσ\_+\)\)\) \((r - \(r\_+\))\)\^\(\(-s\) - \(iσ\_+\)\), \ r\  → \ \(r\_+\)\)  *)

Nmax = 1500; w = wang; sima = (w*rplus - ar*m)/br; ae[0] = 1; 
    ae[1] = -β[0]/α[0]; i = 1; 
    While[i < Nmax, 
      ae[i + 1] = (-β[i]*ae[i] - γ[i]*
                ae[i - 1])/(α[i]); i = i + 1];

Psi = 
        Exp[I*w*r]*(r - rminus)^(-1 - s + I*w + I*sima)*(r - 
                rplus)^(-s - I*sima)*
          Sum[ae[i]*((r - rplus)/(r - rminus))^(i), {i, 0, 
              Nmax - 1}];

R = 10; Psi /. r -> R

(*A sanity check is to compare with the direct numerical integration of Teukolsky' s wavefunction,  
        with the same boundary conditions at the horizon:*)

r1 = rplus + 10^(-6);r2 = 150;Delta = (r - rplus)*(r - rminus);Kr = w*(r^2 + ar^2) - ar*m;VgLe = 
        2*I*s*w*r - ar^2*w^2 - Alm + 
          1/Delta*((r^2 + ar^2)^2*w^2 - 2*ar*m*w*r + ar^2*m^2 + 
                I*s*(ar*m*(2*r - 1) - w*(r^2 - ar^2)));Vg = (Kr^2 - 2*I*s*Kr*(r - 1/2))/Delta + 4*I*s*w*r - 
          Alm + 2*ar*w*m - ar^2*w^2;sol=NDSolve[{Delta*D[y[r], {r, 2}] + (s + 1)*(2*r - 1)*
                      D[y[r], {r, 1}] + Vg*y[r] == 0,y[
                    r1] == (rplus - rminus)^(-1 - s + I*w + 
                          I*sima)*
                    Exp[I*w*r1]*(r1 - rplus)^(-s - I*sima),y'[r1]==(rplus - rminus)^(-1 - s + I*w + 
                          I*sima)*
                    Exp[I*w*r1]*(-s - 
                        I*sima)*(r1 - rplus)^(-s - I*sima - 
                          1)},y,{r, r1, r2},
            AccuracyGoal  ->  20, PrecisionGoal  ->  20,WorkingPrecision  ->  30,MaxSteps -> 340000];

y[R] /. sol

(*the two methods agree to a good precision*)

References

[1] E .  
          Leaver,  "An analytic representation for te quasi-normal modes \
of kerr black holes" .  Proc .  R .  Soc .  
        Lond . A402 :  285 (1985).

[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 and M .  
          Casals,  "Eigenvalues and eigenfunctions of spin-weighted \
spheroidal harmonics in four and higher dimensions",  
      Phys .  Rev .  D73 : 024013 (2006);  
      Erratum - 
        ibid . D73 : 
            109902 (2006) .  [arXiv :  gr - qc/0511111] .