Quasinormal modes of Schwarzschild-(Subscript[AdS, 5]) black branes

Quasinormal modes are eigenfunctions of linearized perturbations around equilibrium or stationary solutions.  Their (generally complex) frequencies characterize the relaxation of perturbations away from a thermal equilibrium state.

The metric (ds^2)=2dt(dr - A(r) dt) + (r^2)((dx^2)+(dy^2)+(dz^2)), with A(r)=1/2 ((r^2)-m/(r^2)), describes a Schwarzschild-(Subscript[AdS, 5]) black brane which, via gauge/string duality, represents a thermal equilibrium state of a 3+1 dimensional quantum field theory at a temperature T=(m^(1/4))/π.  If an infinitesimal metric perturbation (Subscript[δg, xy])=(r^2)(e^(i(qz-ωt)))f(r) is added, then Einstein’s equations require that f(r) be a solution to the equation:

(lhs=(m-(r^4)) f''[r]+(m-5 (r^4)+2 I (r^3) ω) f'[r]/r+((q^2)+3 I r ω) f[r])==0;

The physical domain is r ∈ [(Subscript[r, h]),∞], where (Subscript[r, h])=(m^(1/4)) is the black brane horizon, and r = ∞ is the AdS boundary.  Physical solutions must be regular at the horizon and vanish at the boundary.  For a given value of the wavevector q, such solutions only exist for a discrete set of frequencies ω = (Subscript[ω, i]) known as quasinormal mode frequencies.

Your task: accurately compute the first few (smallest) quasinormal mode frequencies at q=0.

Warm-up :
Is the horizon a singular point of the equation?  Is the boundary?

The horizon is a singularity.  If one rewrites the equation with a unit coefficient of the highest (second) derivative term, then the coefficients of the 0 and 1st derivative terms blow up at r=(Subscript[r, h]):

lhs2=lhs/(m-(r^4))//Collect[#,{f''[r],f'[r],f[r]}]&

The boundary at r = ∞ is also a singularity.  To show this, make a change of variable which maps ∞ to a finite location and examine the result.  Starting with the large-r form of the equation, let z=1/r and again divide by the coefficient of the highest derivative:

lhs2/.q->0/.m->0/.(r^3)ω->0
%/.r->1/z/.f->(g[1/#]&)
%/(z^4)//Expand

The 0 and 1st derivative terms are singular at z = 0, showing the r = ∞ is a singular point of the original ODE.

What is the local behavior of the two linearly independent solutions near r = ∞?

Plug f(r) = (r^α)into the equation, and drop all but the leading power:

lhs/.f->(((#)^α)&)//Collect[#,{(r^α),r}]&

The coefficient of the leading (r^(α+2)) term is -(4α + (α^2)), which vanishes if α = 0 or α = -4.  So as r -> ∞ one solution behaves as (r^0) and approaches a constant, while the other solution vanishes as (r^(-4)).

What is the local behavior of the two linearly independent solutions near r=(Subscript[r, h])?

Plug f(r) = (((r-(Subscript[r, h])))^β)into the equation, and expand around r=(Subscript[r, h]), keeping only the leading term:

lhs/.f->((((#-(m^(1/4))))^β)&)//Series[#,{r,(m^(1/4)),1}]&//Simplify

The coefficient of the leading term,  -4 (m^(3/4)) (β^2)+2 i (Sqrt[m]) β ω, vanishes if β = 0 or β = i ω/2(m^(1/4)).  So as r -> (Subscript[r, h]) one solution approaches a constant, while the other solution, f ~ (((r-(Subscript[r, h])))^(iω/2(Subscript[r, h])))=(e^((iω/2(Subscript[r, h]))ln(r-(Subscript[r, h])))), wiggles infinitely many times and has a branch point singularity at (Subscript[r, h]).

Can you make a redefinition of f(r) which will map the semi-infinite domain into a finite interval (say [0,1]) and at the same time make the only needed boundary conditions regularity at either endpoint?

One choice: Let  f(r) = (r^(-4))g((Subscript[r, h])/r) and set r = (Subscript[r, h])/u.  Scaling out four powers of 1/r means that in the physical solutions of interest g(u) will approach a constant as u -> 0.

lhs3=(u^(-3))(m^(1/2))lhs/.f->(((#)^(-4))g[(m^(1/4))/#]&)/.r->(m^(1/4))/u//Collect[#,{g''[u],g'[u],g[u]},Simplify]&

An equally good alternative is to let f(r) = (r^(-3))h((Subscript[r, h])/r).  Scaling out three powers of u means that in the physical solutions of interest h(u) will vanish linearly as u -> 0.

lhs4=(u^(-1))(m^(1/4))lhs/.f->(((#)^(-3))h[(m^(1/4))/#]&)/.r->(m^(1/4))/u//Collect[#,{h''[u],h'[u],h[u]},Simplify]&

Based on answers to these warm-up questions, appropriately rewrite the equation in a more computationally amenable form, choose a spectral grid, and convert the equation into a finite dimensional generalized eigenvalue problem of the form α v = ω β v, where α and β are MxM matrices, v is the eigenvector, and ω is one of the desired eigenvalues.  (Mathematica, Matlab, and most any linear algebra library has routines to solve such generalized eigenvalue problems.)  Examine convergence with increasing truncation size M.

Load convenient spectral modules:

<<"~/ph506/SpectralCode.m"

Work with the above form lhs3 == 0.  Separate equation into pieces independent of, or proportional to ω, so equation becomes a generalized eigenvalue problem α g = λ β g, with λ = ω/(m^(1/4)):

α=lhs3/.q->0/.ω->0
β=-Coefficient[lhs3,ω/(m^(1/4))]

Use a Chebyshev basis with Lobatto (endpoints + extrema) spectral grid.  Since all basis functions are regular at both horizon and boundary, this vector space of functions automatically excludes the unphysical solutions (which are singular at horizon or boundary) --- no additional conditions are needed.

(Alternatively, if one works with a form of the equation containing inverse powers of u, then one must fix the spectral matrix to avoid nonsense at u=0.  One can (a) use the interior spectral grid to avoid having a grid point at u=0, or (b) use the rescaling f(r) = (r^(-3))h((Subscript[r, h])/r) so that h(u) vanishes linearly at the boundary, plus the exterior grid, and then just zero out the first row of the matrix α.)

Write a little module to do spectral approximation with M basis elements, and compute K smallest eigenvalues:

run[M_,K_]:=run[M,K]=Module[{},
{grid,d,dd,card}=SpectralChebyshevExterior[0.0,1.0,M];
{A,B}={α,β}/.{g''[u]->dd,g'[u]->d,g[u]->IdentityMatrix[M],u->grid};
Eigenvalues[{A,B},-K]//Reverse];

Table[Prepend[10k][run[10k,6]],{k,1,6}]//TableForm

Observe that eigenvalues come in pairs, with opposite real parts and identical imaginary parts.  (Why?)  So drop half of them to make a nicer table:

Table[Prepend[10k][Select[run[10k,6],Re[#]>0&]],{k,1,6}]//TableForm

First eigenvalue is stable by M=20, second one is stable at M=30, and the third one is roughly, but not fully stable.  Looks like precision loss is occurring, so try re-running using extended precision:

Clear[run];
run[M_,K_]:=run[M,K]=Module[{},
{grid,d,dd,card}=SpectralChebyshevExterior[0,1`30,M];
{A,B}={α,β}/.{g''[u]->dd,g'[u]->d,g[u]->IdentityMatrix[M],u->grid};
Eigenvalues[{A,B},-K]//Reverse];

Table[Prepend[10k][Select[run[10k,6],Re[#]>0&]],{k,1,6}]//TableForm

Now third eigenvalue is stable to 9 digits by M=40.