Computing scalar QNM’s of RNdS

Setup

Scalar QNM equation in final rescaled version, using coordinate x = ((r-rh)/(rc-rh)), and depending on parameters rh, rc, Λ, l and ω

eqϕ=(-(1/(3 (-2+rc+rh)))(rc-rh) (-3 l (-2+rc+rh) x-3 (l^2) (-2+rc+rh) x+2 ((rc^3) (x^3)+(rc^2) (-1+rh+2 rh (x^2)-2 (x^3)-rh (x^3))-rc rh (2+4 (x^2)-4 (x^3)+rh (-2+(x^3)))+(rh^2) (-3+4 (x^2)-2 (x^3)+rh (1-2 (x^2)+(x^3)))) Λ)-(I (2 (rh^4) (-1+x)+(rc^5) (x^3)+2 rc (rh^3) (-2+rh+x-rh x)-(rc^4) x (-2+(2+rh) (x^2))+(rc^2) (rh^2) (-2 (3-4 x+(x^3))+rh (4-4 x+(x^3)))-(rc^3) rh (2 x-4 (x^3)+rh (-2+2 x+(x^3)))) ω)/((rc^3)+rc (-2+rh) rh-(rh^2)+(rc^2) (-3+2 rh))-(3 (rc^2) (-2+rc+rh) x (-2 rc (rh^3) (2+rh (-1+x))+2 (rh^4) (-1+x)+(rc^4) (2-rh (-2+x)-2 x) x+(rc^5) (x^2)+(rc^3) rh (4 (-2+x) x+rh (2+2 x-(x^2)))+(rc^2) (rh^2) (-2 (3-2 x+(x^2))+rh (4-2 x+(x^2)))) (ω^2))/((rc-rh) ((((rc^3)+rc (-2+rh) rh-(rh^2)+(rc^2) (-3+2 rh)))^2) Λ)) ϕ[x]+(2 x ((rc^7) (x^3) Λ+(rc^6) ((-1+rh+2 rh (x^2)-5 (x^3)) Λ-3 I (x^3) ω)+(rc^5) ((3+6 (x^3)+(rh^2) (3+2 (x^2)-3 (x^3))+2 rh (-3-5 (x^2)+3 (x^3))) Λ-3 I x (-2+x+3 rh x-4 (x^2)) ω)+(rh^4) ((rh^2) (1-2 (x^2)+(x^3)) Λ+6 I (((-1+x))^2) ω-rh ((3-4 (x^2)+2 (x^3)) Λ+3 I (((-1+x))^2) ω))-rc (rh^3) ((rh^3) (1-2 (x^2)+(x^3)) Λ+12 I (-1+x) ω-(rh^2) (2 (3-3 (x^2)+(x^3)) Λ+3 I (((-1+x))^2) ω)+rh ((5-2 (x^3)) Λ+3 I (5-8 x+3 (x^2)) ω))+(rc^4) ((rh^3) (2-4 (x^2)) Λ-6 I x (2-x+2 (x^2)) ω+rh ((5+12 (x^2)-14 (x^3)) Λ+3 I (13-4 x) (x^2) ω)+(rh^2) ((-9+2 (x^2)+7 (x^3)) Λ+3 I (1-2 x-4 (x^2)+2 (x^3)) ω))+(rc^2) (rh^2) ((rh^3) (-3+2 (x^2)) Λ-6 I (-3+8 x-6 (x^2)+2 (x^3)) ω-3 (rh^2) ((-3+(x^3)) Λ+I (-3+6 x-4 (x^2)+(x^3)) ω)+3 I rh (2 I Λ+(-9+18 x-12 (x^2)+4 (x^3)) ω))+(rc^3) rh ((rh^3) (-2-4 (x^2)+3 (x^3)) Λ+12 I x (1-4 x+2 (x^2)) ω+(rh^2) (-8 (x^3) Λ+2 (x^2) (8 Λ+3 I ω)+9 I ω-18 I x ω)+rh (2 (3-8 (x^2)+4 (x^3)) Λ-3 I (5-10 x-4 (x^2)+4 (x^3)) ω))) ϕ'[x])/(3 ((rc^4)-(-2+rh) (rh^2)+(rc^3) (-5+3 rh)+rc rh (4-5 rh+(rh^2))+3 (rc^2) (2-3 rh+(rh^2))))+(1/6) (rc-rh) (-1+x) (x^2) (rc+3 rh+((Sqrt[(-2+rc+rh) ((rc^3)-(rc^2) (-6+rh)-rc (-4+rh) rh+(rh^2) (6+rh))])/(-2+rc+rh))+2 rc x-2 rh x) (rh+(1/2) (rc+rh-((Sqrt[(-2+rc+rh) ((rc^3)-(rc^2) (-6+rh)-rc (-4+rh) rh+(rh^2) (6+rh))])/(-2+rc+rh)))+rc x-rh x) Λ ϕ''[x];

Get (rh,rc) from (Λ,Q)

There are only 2 independent parameters, which we want to take as Λ and Q (in units where M=1, equivalently: ((ΛM)^2) and Q/M)

Λsol=((3 (-2+rc+rh))/((rc+rh) ((rc^2)+(rh^2))));
Qsqsol=((rc rh (-(rc^2) (-2+rh)-rc (-2+rh) rh+2 (rh^2)))/((rc+rh) ((rc^2)+(rh^2))));

ClearAll[ΛQtorhrc]
ΛQtorhrc[Λval_,Qval_,prec_ : 40]:=Block[{Λ=SetPrecision[Λval,prec+10],Qsq=SetPrecision[(Qval^2),prec+10]},
({rh,rc}/.NSolve[{Λ==Λsol,Qsq==Qsqsol},{rh,rc},WorkingPrecision->prec])//
Select[#,#[[1]]>0&&#[[2]]>0&&#[[2]]>#[[1]]&]&//
Flatten//Union//
DeleteDuplicatesBy[SetPrecision[#,10]&]//
Take[#,-2]&//(Thread[{rh,rc}->#])&
]~Join~{Λ->SetPrecision[Λval,prec]}

Numerics

Load package (for download and usage instructions, go to https://github.com/APJansen/QNMspectral )

Needs["QNMspectral`"]

Code to compute modes

RNdSmode[lval_,Λ_,Q_,acc_ : {{40},{60}}]:=RNdSmode[lval,Λ,Q]=Block[{$QNMQuiet=True},
GetAccurateModes[eqϕ/.ΛQtorhrc[Λ,Q]/.l->lval,Sequence@@acc]
]

Some Examples

Example 1: To be compared with figure 1.

modesTest=RNdSmode[1,(6/100),(2/10)];

modesTest//ShowModes

Note: the unstable mode does not satisfy the boundary conditions: (they are scaled to go to 0 as x->0)

modesTestEfs=GetAccurateModes[eqϕ/.ΛQtorhrc[(6/100),(2/10)]/.l->1,{40},{60},Eigenfunctions->"Later"];

modesTestEfs[[1;;3]]//PlotEigenfunctions

Of the correct modes, the dominant is a photon sphere mode, and then a purely imaginary dS mode, as we expect from fig. 1

Example 2: to be compared with figure 3

Some necessary parameters

This is the maximal Q^2 as a function of Λ

QsqMax=((Re[((((((-2))^(2/3))-(((-2))^(1/3)) (((I (Sqrt[2-9 Λ])+3 (Sqrt[Λ])))^(2/3))+6 (((3 (Λ^2)+(Sqrt[(Λ^3) (-2+9 Λ)])))^(1/3))) ((1+I (Sqrt[3])) (((3 (Λ^2)+(Sqrt[(Λ^3) (-2+9 Λ)])))^(2/3))+Λ Root[16+((#1)^3)&,2]))/(Λ (((3 (Λ^2)+(Sqrt[(Λ^3) (-2+9 Λ)])))^(2/3))))])/(8 (2^(2/3))));

Qmax=Sqrt[QsqMax]/.Λ->6/100//N[#,50]&

This is the surface gravity of the inner horizon:

κi=((Sqrt[(-2+rc+rh) ((rc^3)-(rc^2) (-6+rh)-rc (-4+rh) rh+(rh^2) (6+rh))]) (2 rc-(rc^2)+6 rh-4 rc rh-3 (rh^2)+(Sqrt[(-2+rc+rh) ((rc^3)-(rc^2) (-6+rh)-rc (-4+rh) rh+(rh^2) (6+rh))])) (6 rc-3 (rc^2)+2 rh-4 rc rh-(rh^2)+(Sqrt[(-2+rc+rh) ((rc^3)-(rc^2) (-6+rh)-rc (-4+rh) rh+(rh^2) (6+rh))])) Λ)/(6 (-2+rc+rh) (((-2 rc+(rc^2)-2 rh+2 rc rh+(rh^2)-(Sqrt[(-2+rc+rh) ((rc^3)-(rc^2) (-6+rh)-rc (-4+rh) rh+(rh^2) (6+rh))])))^2));

κiHere=κi/.ΛQtorhrc[(6/100),(992/1000)Qmax];

From the bottom left of Fig 3. we see that at Q/Qmax ~ 0.9920:

The dominant mode is a photon sphere mode with large l (computed at l=10), at  Im(ω)/(SubMinus[κ])~ -0.5

modesExample2=RNdSmode[10,(6/100),(992/1000)Qmax];

modesExample2/κiHere//MakeTable

Note again the unstable mode doesn’t satisfy the boundary conditions.

There is a dS mode with l=1 at Im(ω)/(SubMinus[κ])~-1

modesExample2dS=RNdSmode[1,(6/100),(992/1000)Qmax];

modesExample2dS/κiHere//MakeTable

Note indeed the purely imaginary mode at ~ -1, and also that the photon sphere mode with l=1 is quite close to that with l=10 and already dominates over the dS mode.

There is a near-extremal mode with l=0 in between:

modesExample2NE=RNdSmode[0,(6/100),(992/1000)Qmax];

modesExample2NE/κiHere//MakeTable

This is not so successful, the imaginary mode at -0.8261 might be what we’re looking for.
Note though that it’s less accurate than the dominant modes of the other families, so let’s be sure that this one is correct by improving the accuracy:

modesExample2NEv2=RNdSmode[0,(6/100),(992/1000)Qmax,{{120},{160}}];//AbsoluteTiming

modesExample2NEv2/κiHere//MakeTable[#,NModes->10]&

It is now more accurate but agrees with the previous result, so this is correct.