aainstitutetext: School of Mathematics and Statistics, University College
Dublin, Belfield, Dublin 4, Irelandbbinstitutetext: Centre for Research and
Development in Mathematics and Applications (CIDMA),
Department of Mathematics, University of Aveiro, 3810-193 Aveiro, Portugal

#  The pseudospectra of black holes in AdS

Bradley Cownden<EMAIL_ADDRESS>a Christiana Pantelidou
<EMAIL_ADDRESS>b and Miguel Zilhão<EMAIL_ADDRESS>
###### Abstract

We study the stability of quasinormal modes (QNMs) in electrically charged
black brane spacetimes that asymptote to AdS by means of the pseudospectrum.
Methodologically, we adopt ingoing Eddington–Finkelstein coordinates to cast
QNMs in terms of a generalised eigenvalue problem involving a non-selfadjoint
operator; this simplifies the computation significantly in comparison with
previous results in the literature. Our analysis reveals spectral instability
for (neutral) scalar as well as gravitoelectric perturbations. This is
significant as it contains information about the dual strongly coupled system
via the gauge/gravity correspondence. Particular attention is given on the
hydrodynamic modes, where we see the spectral lines crossing to the upper half
plane indicating transient instability. We also investigate the asymptotic
structure of pseudospectral contour levels and we find remarkable universality
across all sectors, persistent in the extremal limit.

##  1 Introduction

For self-adjoint operators present in non-dissipative systems, the spectral
theorem guarantees the stability of the spectrum under perturbations, meaning
that small-scale perturbations of the operator will lead to a small movement
of the operator’s eigenvalues on the complex plane. In contrast, in the case
of dissipative systems such as black holes, the spectral theorem is not
applicable leading to lack of completeness in the set of eigenfunction and
their orthogonality, potentially giving rise to instabilities. Such
instabilities will manifest themselves as strong sensitivity of the
eigenvalues (i.e. large eigenvalue migrations) to small-scale perturbations of
the operator.

Mathematical tools for detecting such instabilities without the need to
compute eigenvalues have been extensively used in the literature, including
fields like quantum mechanics and hydrodynamics. One such tool is the
pseudospectrum trefethen2005spectra . A particularly interesting application
of the pseudospectrum comes from the study of the stability of fluid flows
treften90hydro . It is well known that fluid flows that are laminar (stable)
at low speeds become unstable and then turbulent at higher speeds. This
phenomenon was traditionally investigated by linearising the Navier-Stokes
equation and testing for unstable eigenvalues of the linearised problem, but
the results agreed poorly in many cases with experiments. Considerations of
the pseudospectrum of the linearised problem resolved this tension
treften90hydro , by suggesting that small perturbations of the smooth flow may
be amplified substantially by linear mechanisms.

Motivated by the seminal work of Nollert and Price Nollert:1998ys ;
Nollert:1996rf , the pseudospectrum was recently introduced in the context of
black hole physics in asymptotically flat Jaramillo:2020tuu ;
Jaramillo:2021tmt ; Destounis:2021lum , in asymptotically anti-de Sitter (AdS)
Arean:2023ejh  and in asymptotically de Sitter (dS) geometries Sarkar:2023rhp
aiming to address black hole quasinormal mode (QNM) instabilities. In the case
of asymptotically flat spacetimes, the study of the pseudospectrum was
motivated by the black hole spectroscopy program. In the near future, high-
precision gravitational wave observations from binary black hole mergers will
allow for a detailed inspection of the spectra of astrophysical black holes
during the ringdown phase Nollert:1999ji ;
lisaconsortiumwaveformworkinggroup2023waveform . Since these black holes are
not isolated objects, it is necessary to understand if the QNM spectrum itself
is affected by the astrophysical environment, quantum corrections, or other
generic modifications. The pseudospectrum analysis for the Schwarzschild and
Reissner–Nordström (RN) black holes was carried out in Jaramillo:2020tuu  and
Destounis:2021lum  respectively, and found spectral instability of scalar and
polar gravitoelectric quasinormal modes in subextremal and extremal black
holes. The response of the Schwarzschild black hole spectrum to specific
potential perturbations was studied in and it was shown that even the
fundamental mode can be unstable under a certain class of perturbations
Cheung:2021bol ; Jaramillo:2020tuu . This raised concerns about the robustness
of the spectroscopy program, which let to a detailed study of the physical
behaviour of time-domain signals Berti:2022xfj . The pseudospectrum was also
studied for exotic horizonless objects Boyanov:2022ark .

In the case of asymptotically AdS black branes in 5 dimensions, the
pseudospectrum for neutral scalar and transverse electromagnetic perturbations
of the planar Schwarzschild black hole was studied in Arean:2023ejh . In this
case the aim was to shed a new light on the behaviour of strongly coupled
quantum many-body systems modelled by the gauge/gravity duality. Compared to
the asymptotically flat case, dissipation in AdS is reduced given that the AdS
boundary works as a box not allowing energy to escape to infinity. However,
due to the presence of the event horizon, the system is still dissipative and
the spectral operator is still non-self adjoint. As such, similarly to the
asymptotically flat case, the pseudospectrum revealed spectral instability.

The computation we perform in this paper largely follows the computations
already done in the literature. An important difference is the coordinate
basis used. In asymptotically flat spacetimes the pseudospectrum computation
is carried out in hyperboloidal coordinates Jaramillo:2020tuu ;
Destounis:2021lum ; Boyanov:2022ark . Other than compactifying the domain,
this set of coordinates has the added benefit of making manifest ingoing
boundary conditions at the horizon and outgoing boundary conditions at future
null infinity. In asymptotically AdS, Arean:2023ejh  used a set of “regular
compact coordinates” that asymptote to ingoing Eddington–Finkelstein
coordinates close to the horizon. This made manifest the ingoing boundary
conditions at the horizon, while at the boundary of AdS they imposed that the
source for the scalar field vanishes. In both of these cases, the equations of
motion for the perturbations were second order in time and the spectrum was
obtained by solving an eigenvalue problem on a spectral grid. Note that in
Arean:2023ejh  the authors also considered formulating the problem in purely
ingoing Eddington–Finkelstein coordinates, but ultimately decided to proceed
using the above coordinate system.

In this work, we expand the previous work done for asymptotically AdS black
holes by (i) considering perturbations of a massless scalar as well as
electromagnetic and gravitational perturbations around the Reissner–Nordström
(RN) black hole and (ii) formulating the pseudospectrum computation as a
generalised eigenvalue problem in ingoing Eddington–Finkelstein coordinates.
The QNM spectrum for black holes/branes in AdS has been extensively studied in
the literature Policastro:2002se ; Kovtun:2005ev , including more mathematical
work Warnick:2013hba . Thanks to the coordinate system used, ingoing boundary
conditions at the horizon become manifest.111In addition to ingoing boundary
conditions at the horizon, we also need to ensure the absence of sources at
infinity. Moreover, the equations of motion for the perturbation are now first
order in time and in fact the time derivative appears as a mixed $(v,z)$
derivative, where $v$ in the ingoing Eddington–Finkelstein coordinate and $z$
is a compactified radial coordinate. Solving for the spectrum can then be
recast as a generalised eigenvalue problem instead of a standard one
trefethen2005spectra . Computations of the pseudospectrum in the case of
generalised eigenvalue problems are less developed in the literature. In
particular, there are different definitions of the pseudospectrum depending on
which operators are perturbed and one needs to resort to physical arguments in
order to select the definition most suitable to the problem at hand. In our
case, we perturb only the physically meaningful part of the operator.

Our results for the scalar field are presented in section 4.1 as function of
the momentum $k$ and charge $Q$ (or equivalently, the temperature $T$) and
naturally extend the results of Arean:2023ejh . In particular, we see that all
modes in the spectrum are unstable for sufficiently large perturbations.

The pseudospectrum for metric and gauge field perturbations is computed in
section 4.2. Note that in the case of the AdS-RN black hole the gravitational
and electromagnetic perturbations do not decouple from each other. However, it
is well understood that they can be split into three sectors: the scalar,
vector and tensor222The tensor channel is present only in bulk dimensions
$D\geq 5$. sectors, depending on how the perturbations transform under
rotations on the plane transverse to the momentum. Each of the scalar and
vector sectors contain two channels: in the scalar sector these are the sound
and charge diffusion channels; in the vector sector these are the shear and
transverse gauge channels. The tensor sector, on the other hand, only includes
one channel and in fact behaves as a neutral massless scalar. Thus, its
pseudospectrum is captured by the computation of section 4.1. Our
pseudospectral analysis demonstrates generic instability in all channels.
Particular attention is given to the channels that contain hydrodynamic modes,
namely modes that reside at the origin of the complex frequency plane at zero
momentum; these are the shear, sound and charge diffusion. We see that
pseudospectral contour lines around the hydrodynamic modes cross to the upper
half plan, which translates to an unstable perturbed spectrum and hints
towards some connection with transient black hole instabilities. We also study
the pseudospectrum in the extremal limit and compare and contrast with similar
studies in different asymptotics Destounis:2021lum ; Yang:2022wlm .

We have also investigated the asymptotic behaviour of pseudospectral contour
levels at large real values of the spectral frequency, finding universality
across different sectors. Similar universality was seen in 4-dimensional
asymptotically flat spacetimes with either Schwarzschild Jaramillo:2020tuu ;
Jaramillo:2021tmt  or RN black holes Destounis:2021lum .

The rest of the paper is organised as follows. In section 2 we give the formal
definition for the spectrum and pseudospectrum for a standard eigenvalue
problem as well as a generalised eigenvalue problem – as explained earlier,
when formulating the calculation of the spectrum in ingoing
Eddington–Finkelstein coordinates one needs to solve a generalised eigenvalue
problem. An important measure when constructing the pseudospectrum is the
choice of norm, and in section 2 we also motivate our choice of the _energy
norm_ as the fundamental scale of the system. In section 3 we describe details
of our numerical implementation, including the discretisation method used and
the matrix representation of the energy norm in terms of the Gram matrix. Our
results for the pseudospectrum of the 5-dimensional planar AdS-RN black brane
are then presented in section 4. In particular, in section 4.1 we discuss the
pseudospectrum of a neutral, massless scalar field as well as spectral
stability under specific potential perturbations, while in 4.2 we consider
gravitoelectric perturbation in all three channels. The universality of the
asymptotic behaviour of the pseudospectral contour lines is discussed in 4.3.
Finally, we conclude with an outlook and future directions in section 5. There
are also two appendices. In appendix A we study the scalar spectrum of global
AdS4, where we find spectral stability associated with the lack of energy
dissipation, while in appendix B we discuss the numerical convergence of our
QNM computation.

##  2 Spectra and pseudospectra

The notion of _spectrum_ of an operator is a familiar one in many areas of
Physics, Mathematics and Engineering. It can be understood as a generalisation
of the set of eigenvalues of a matrix. Recall that for a matrix $A$, its
eigenvalues are the values of $\lambda$ that satisfy the equation

$\det(A-\lambda\mathbb{1})=0$

where $\mathbb{1}$ is the identity matrix. If $\lambda$ is an eigenvalue of
$A$, the operator $A-\lambda\mathbb{1}$ is not one-to-one and therefore its
inverse $(A-\lambda\mathbb{1})^{-1}$ is not defined. We then define the
_spectrum of A_ as

$\sigma(A)=\\{\lambda\in\mathbb{C}:\lambda\mathbb{1}-A\text{ is not
invertible}\\}.$ (1)

It is also useful to define the _resolvent_
$R_{A}(\lambda)=(\lambda\mathbb{1}-A)^{-1}$ of the operator $A$ and note that
if the spectrum $\sigma(A)$ is empty $R_{A}(\lambda)$ is defined everywhere in
the complex plane.

The spectrum is an intrinsic property of an operator, and it is an incredibly
useful concept. Nevertheless, there are classes of problems where the
associated methods may fail to give a useful answer. In particular, it is only
for _normal_ operators $N$ satisfying $NN^{\dagger}=N^{\dagger}N$, where †
denotes the adjoint, that one can show that the spectrum of the operator is
_stable_. For non-normal operators, in general, small perturbations to the
operator can drastically alter its spectrum.

To study such problems, Trefethen and collaborators reddy1990lax ;
trefethen2005spectra  introduced the concept of _pseudospectrum_. To motivate
its definition let us start by noting that, from definition (1), whenever
$\lambda$ is an eigenvalue of $A$, the operator $\lambda\mathbb{1}-A$ is
singular. Equivalently, the resolvent
$R_{A}(\lambda)=(\lambda\mathbb{1}-A)^{-1}$ is not defined. Let us then define
the _$\epsilon$ -pseudospectrum_ as the set of points in the complex plane
such that trefethen2005spectra

$\sigma^{\epsilon}(A)=\\{\lambda\in\mathbb{C}:\|R_{A}(\lambda)\|=\|(\lambda\mathbb{1}-A)^{-1}\|>1/\epsilon\\}$
(2)

for some choice of norm $\|\cdot\|$. Note that $\|R_{A}(\lambda)\|=\infty$
whenever $\lambda$ is an eigenvalue of $A$ and thus, with this convention,
$\sigma^{0}(A)=\sigma(A)$.

It can also be shown reddy1990lax  that the definition (2) is equivalent to

$\sigma^{\epsilon}(A)=\\{\lambda\in\mathbb{C},\,\exists\,\delta A\text{ with
}\|\delta A\|<\epsilon:\lambda\in\sigma(A+\delta A)\\},$ (3)

which means that the points in $\sigma^{\epsilon}(A)$ are the eigenvalues of
an operator $A+\delta A$, that is, of some $\epsilon$-perturbation of the
operator $A$. Note however that, while the spectrum is an intrinsic property
of the operator, the pseudospectrum is not, as it depends also on the choice
of the operator norm. Indeed, the choice of norm is crucial since it
quantifies how large are the perturbations $\delta A$, and thus a small
perturbation in a given norm could actually correspond to a large perturbation
under another norm – see also the relevant discussion in Jaramillo:2020tuu ;
Gasperin:2021kfv .

A related measure of the stability of a spectrum is the _condition number_
$\kappa$ of each eigenvalue. Let us consider a matrix $A\in\mathbb{C}^{N\times
N}$ with a set of distinct eigenvalues $\\{\lambda_{1},\ldots,\lambda_{N}\\}$,
which implies the existence of a set of right eigenvectors ($v_{j}$) and left
eigenvectors ($u_{j}$) that satisfy

$\displaystyle Av_{j}=\lambda_{j}v_{j}\qquad\text{and}\qquad
u^{*}_{j}A=\lambda_{j}u^{*}_{j}\,.$ (4)

Under the perturbation $A\to A(t)=A+t\,\delta A$ for333The definition of the
matrix and vector norms will become central to our discussion in later
sections. For now we leave the definition of $||\cdot||$ general. $\|\delta
A\|=1$ and $t\in\mathbb{C},|t|<1$, the perturbed eigenpair
$(\lambda_{1}(t),v_{1}(t))$ satisfy

$\displaystyle A(t)v_{1}(t)=\lambda_{1}(t)v_{1}(t)\,.$ (5)

Expanding $\lambda_{1}(t)$ and $v_{1}(t)$ as Taylor series in $t$
trefethen2005spectra ,

$\displaystyle\lambda_{1}(t)$
$\displaystyle=\lambda_{1}+\alpha_{1}t+\alpha_{2}t^{2}+\ldots$ (6)
$\displaystyle v_{1}(t)$
$\displaystyle=v_{1}+\sum^{N}_{j=1}\left(t\beta^{(1)}_{j}+t^{2}\beta^{(2)}_{j}+\ldots\right)v_{j}$
(7)

we can take the inner product of (5) with $u^{*}_{1}$ and see that

$\displaystyle\lambda_{1}(t)\langle
u^{*}_{1},v_{1}(t)\rangle=\lambda_{1}(t)\langle u^{*}_{1},v_{1}\rangle$
$\displaystyle=\langle u^{*}_{1},A(t)v_{1}(t)\rangle$ $\displaystyle=\langle
u^{*}_{1},Av_{1}(t)\rangle+t\langle u^{*}_{1},\delta Av_{1}(t)\rangle$
$\displaystyle=\lambda_{1}\langle u^{*}_{1},v_{1}\rangle+t\langle
u^{*}_{1},\delta Av_{1}(t)\rangle\,.$

Thus, the perturbation of $A$ has moved the eigenvalue $\lambda_{1}$ by

$\displaystyle|\lambda_{1}-\lambda_{1}(t)|$ $\displaystyle=|t|\frac{|\langle
u^{*}_{1},\delta Av_{1}(t)\rangle|}{|\langle
u^{*}_{1},v_{1}\rangle|}\leq|t|\frac{\|u_{1}\|\,\|v_{1}(t)\|}{|\langle
u^{*}_{1},v_{1}\rangle|}$
$\displaystyle=|t|\frac{\|u_{1}^{*}\|\,\|v_{1}\|}{|\langle
u^{*}_{1},v_{1}\rangle|}+\mathcal{O}(t^{2})\,,$ (8)

where we have used the Cauchy-Schwarz inequality and the fact that $\|\delta
A\|=1$. Defining the condition number of the $i$-th eigenvalue by

$\displaystyle\kappa_{i}=\frac{\|u_{j}\|\,\|v_{j}\|}{|\langle
u^{*}_{j},v_{j}\rangle|}$ (9)

we see that the magnitude of $\kappa_{i}$ gives us an indication of the
stability of $\lambda_{i}$ as it determines whether the shift in the
eigenvalue remains of the order of the size of the perturbation $|t|$. Indeed,
for normal operators with $A^{\dagger}=A$, $u_{j}$ is colinear with $v_{j}$
and so $\kappa_{j}=1$. However, non-normal operators, i.e. unstable operators,
may have non-colinear ${(u_{j},v_{j})}$, resulting in arbitrarily small values
for $|\langle u^{*}_{j},v_{j}\rangle|$ and arbitrarily large values of
$\kappa_{j}$.

###  2.1 Pseudospectra for matrix pencils

The definitions (2) and (3) stem from an eigenvalue problem of the form
$Av=\lambda v$. As we will see, the computation of the spectrum (and,
correspondingly, also of the pseudospectrum) for black holes in AdS becomes
more natural in ingoing Eddington–Finkelstein coordinates. This coordinate
system naturally leads to a _generalised_ eigenvalue problem, in the form

$Av=\lambda Bv.$ (10)

The matrix $A-\lambda B$ is also known as a _matrix pencil $(A,B)$_.

When considering the pseudospectrum in the form of a generalised eigenvalue
problem, as we will do in this work, the available techniques are not as well-
developed as for the standard eigenvalue case. Herein we will follow the
approach of reference van1997pseudospectra , which we now outline.

Take the generalised eigenvalue problem, equation (10), with $B$ potentially
singular. The (generalised) eigenvalues of the matrix pencil $(A,B)$ are the
set of complex numbers $\lambda$ for which $\det(A-\lambda B)=0$. This set of
eigenvalues, as usual, is called the _spectrum_ and is denoted by
$\sigma(A,B)$. Analogously, the _$\epsilon$ -pseudospectrum of the matrix
pencil $(A,B)$_ is defined as riedel1994generalized ; van1997pseudospectra

$\sigma^{\epsilon}(A,B)=\\{\lambda\in\mathbb{C}:\|(\lambda
B-A)^{-1}\|>1/\epsilon\\}$ (11)

for some choice of norm, which naturally reduces to (2) when $B=\mathbb{1}$.
Furthermore, it is shown in van1997pseudospectra  that definition (11) is
equivalent to

$\sigma^{\epsilon}(A,B)=\\{\lambda\in\mathbb{C},\,\exists\,\delta A\text{ with
}\|\delta A\|<\epsilon:\lambda\in\sigma(A+\delta A,B)\\}.$ (12)

In other words, $\sigma^{\epsilon}(A,B)$ is the actual spectrum for some
$\epsilon$-perturbation of $A$, and we recover the (usual) definition (3) for
the pseudospectra of an operator $A$ when $B=\mathbb{1}$. This definition does
not consider perturbations to the operator $B$ and is therefore not as general
as others considered in the literature trefethen2005spectra . However, as we
will argue later on, for our present purposes this is the definition of
interest since all relevant physical perturbations will be fully contained
within the operator $A$.

The notation of the condition number also extends to matrix pencils. While
defining a condition number in the case of a singular $B$ matrix is more
subtle (see Stewart1975GershgorinTF ; James1993  for discussions on this
topic), when $B$ is non-singular the definition of a condition number follows
from applying the mechanism outlined in section 2 to the matrix product
$B^{-1}A$.

###  2.2 Energy norm

Following the discussion of the previous section, we need to define an
appropriate norm $\|\cdot\|$ to quantify the size of perturbations. A natural
definition is to use the _energy norm_ associated with the energy of each mode
present in the spectrum Jaramillo:2020tuu .

An intuitive way to obtain such a norm is as follows. We start by considering
the equation governing the spectrum of interest, which in a wide range of
cases (and in particular for all cases we will consider) reduces to a complex
scalar field that obeys the wave equation in 1+1 Minkowski spacetime with a
potential $V$.444In the context of black hole physics, it can be shown that a
wide class of perturbations around a black hole background take this form. The
coordinate $y$ is identified with the tortoise coordinate $y=r_{*}$ and the
potential depends, among other things, on the concrete black hole solution
being perturbed and the specific type of perturbation Berti:2009kk .  In
Cartesian coordinates $(t,y)$ this reads

$\left(\square_{{}_{2}}-V\right)\phi(t,y)=0,$ (13)

where $\square_{{}_{2}}=-\partial^{2}_{t}+\partial^{2}_{y}$. Multiplying the
latter by $\partial_{t}\,\bar{\phi}$, where ${}~{}\bar{~}{}{}~{}$ denotes
complex conjugation, and summing with its complex conjugate, we obtain
bartnik1996wave

$\partial_{t}\left(\partial_{t}\bar{\phi}\,\partial_{t}\phi+\partial_{y}\bar{\phi}\,\partial_{y}\phi+V\,\bar{\phi}\phi\right)=\partial_{y}\left(\partial_{t}\bar{\phi}\,\partial_{y}\phi+\partial_{t}\phi\,\partial_{y}\bar{\phi}\right).$

Integrating both sides in $y$, we find that solutions of equation (13) obey

$\frac{d}{dt}E(t)=\frac{1}{2}\partial_{t}\bar{\phi}\,\partial_{y}\phi\Big{|}_{{y=a}}^{{y=b}}+\frac{1}{2}\partial_{t}\phi\,\partial_{y}\bar{\phi}\Big{|}_{{y=a}}^{{y=b}}$
(14)

where

$E(t)=\frac{1}{2}\int_{a}^{b}\left(\partial_{t}\bar{\phi}\,\partial_{t}\phi+\partial_{y}\bar{\phi}\,\partial_{y}\phi+V\,\bar{\phi}\phi\right)dy$
(15)

and $a,b$ are the boundaries of our domain. $E(t)$ can be interpreted as the
“total energy” contained in a spatial slice $t=\text{constant}$, and we see
from equation (14) that the energy is conserved if there are no modes entering
or leaving the physical domain at $y=a,b$.

Another way of obtaining the energy norm that does not rely on a coordinate
system was introduced in Jaramillo:2020tuu . There one starts by defining the
usual stress-energy tensor of the scalar field $\phi$

$T_{ab}=\frac{1}{2}\nabla_{a}\bar{\phi}\nabla_{b}\phi+\frac{1}{2}\nabla_{a}\phi\nabla_{b}\bar{\phi}-\frac{\eta_{ab}}{2}\left(\nabla_{c}\phi\nabla^{c}\bar{\phi}+V\phi\bar{\phi}\right),$
(16)

where $\eta_{ab}$ is the 2-dimensional Minkowski metric (in arbitrary
coordinates). We can then define the energy on a hypersurface $\Sigma$ through

$E=\int_{\Sigma}T_{ab}\,\xi^{a}\kappa^{b}\,d\Sigma\,.$ (17)

$\xi$ is the timelike Killing vector associated with stationarity and $\kappa$
is the normal vector to the hypersurface $\Sigma$. In $(t,y)$ coordinates and
when the hypersurface $\Sigma$ is chosen to correspond to a
$t=\text{constant}$ slice, $\xi=\partial_{t}$ and it is easy to show that (17)
reduces to expression (15).

It is illuminating to also consider (17) on a null hypersurface.555 In the
context of black holes, a natural choice for a null surface is a
$v=\text{constant}$ slice, where $v$ corresponds to the ingoing
Eddington–Finkelstein coordinate defined as $v=t+y$, where $y=r_{*}$ is the
tortoise coordinate. In particular, let us consider the coordinate system
$(v,Y)$ where

$\displaystyle v=t+y\,,\qquad Y=y$ (18)

and the metric element reduces to

$ds^{2}=-dv^{2}+2dv\,dY\,.$

Let us also choose the hypersurface $\Sigma_{v}$ to correspond to
$v=\text{constant}$. In this case, $\xi=\partial_{t}=\partial_{v}$ and
$\kappa$ is the normal vector to $\Sigma_{v}$. Since
$\kappa_{a}=-\partial_{a}v$, we have $\kappa=-\partial_{Y}$. We then have
$T_{ab}\,\xi^{a}\kappa^{b}=-T_{v\small{Y}}=\frac{1}{2}\partial_{Y}\phi\partial_{Y}\bar{\phi}+\frac{1}{2}V\phi\bar{\phi}$.
Finally, since $d\Sigma_{v}=-\kappa_{v}dY=dY$, we obtain (19).

$E[\phi]=\frac{1}{2}\int_{\Sigma_{v}}\left(\partial_{Y}\bar{\phi}\,\partial_{Y}\phi+V\,\bar{\phi}\phi\right)dY.$
(19)

This expression motivates the definition of the scalar product

$\langle\phi_{1},\phi_{2}\rangle_{{}_{E}}=\int_{a}^{b}\left(\partial_{Y}\bar{\phi}_{1}\,\partial_{Y}\phi_{2}+V\bar{\phi}_{1}\,\phi_{2}\right)dY,$
(20)

where $a,b$ are now the boundary points in the $Y$ coordinate.

##  3 Matrix formulation

In this section we give details about the discretisation approach we follow in
order to construct the spectrum and the pseudospectrum of our systems.

###  3.1 Choosing the grid

Based on the geometry of the problem, different choices of abscissa are
available Grandclement:2007sb . The most common choices are the Gauss-
Chebyshev points, given by ${x^{N}_{i}=\cos(\pi(2i-1)/2N)}$, ${i=1,\ldots,N}$,
and the Chebyshev-Lobatto points, ${x^{N}_{i}=\cos(\pi i/N),\,i=0,\ldots,N}$.
Note that the Chebyshev-Lobatto abscissa includes the endpoints $x=\pm 1$
while the Gauss-Chebyshev do not; this is relevant whenever the functions
being evaluated are singular at the ends of the grid. In this work we choose
to use the Chebyshev-Lobatto points (denoted $\bar{x}_{i}$) in part because we
have specific boundary conditions that must be enforced and also because of a
convenient overlap between a grid of $N$ points, $\bar{x}^{N}_{i}$, and a grid
of $\bar{x}^{2N}_{i}$ that allows for easy interpolation. This will be
important during the numerical evaluation of the energy inner product, as we
will explain below.

First, for completeness, we note that there are known expressions for the
discretised first and second derivative operators in terms of simple functions
of the gridpoints $\bar{x}_{i}$ when using Chebyshev collocation methods. In
particular, the first and second derivative matrices for the Chebyshev-Lobatto
abscissa are

$\displaystyle\mathbb{D}^{(1)}$
$\displaystyle=\begin{dcases}-\frac{2N^{2}+1}{6}&\text{if }i=j=0\\\
\frac{2N^{2}+1}{6}&\text{if }i=j=N\\\
-\frac{\bar{x}_{j}}{2(1-\bar{x}_{j}^{2})}&\text{if }0<i=j<N\\\
\frac{\kappa_{i}}{\kappa_{j}}\frac{(-1)^{i-j}}{(\bar{x}_{i}-\bar{x}_{j})}&\text{if
}i\neq j\end{dcases}$ (21)

and

$\displaystyle\mathbb{D}^{(2)}$
$\displaystyle=\begin{dcases}\frac{N^{4}-1}{15}&\text{if }i=j=0,N\\\
-\frac{1}{(1-\bar{x}_{i}^{2})^{2}}-\frac{N^{2}-1}{3(1-\bar{x}_{i}^{2})}&\text{if
}0<i=j<N\\\
\frac{2}{3}\frac{(-1)^{j}}{\kappa_{j}}\left[\frac{(2N^{2}+1)(1-\bar{x}_{j})-6}{(1-\bar{x}_{j})^{2}}\right]&\text{if
}i=0,~{}0<j\leq N\\\
\frac{2}{3}\frac{(-1)^{N+j}}{\kappa_{j}}\left[\frac{(2N^{2}+1)(1+\bar{x}_{j})-6}{(1+\bar{x}_{j})^{2}}\right]&\text{if
}i=N,~{}0\leq j<N\\\
\frac{(-1)^{i-j}}{\kappa_{j}}\left[\frac{\bar{x}_{i}^{2}+\bar{x}_{i}\bar{x}_{j}-2}{(1-\bar{x}_{i}^{2})(\bar{x}_{j}-\bar{x}_{i})^{2}}\right]&\text{if
}i\neq j,~{}0<i<N\end{dcases}$ (22)

where

$\kappa_{i}=\begin{dcases}2,&i=0,~{}N\\\ 1,&0<i<N\end{dcases}$

###  3.2 Chebyshev integration

Numerical evaluation of the energy inner product defined in (20) involves
evaluating integrals of the form

$\displaystyle\mathcal{I}(f)=\int^{1}_{-1}f(x)dx\,.$ (23)

To do so, we first expand the function in a basis of Chebyshev
polynomials666Exact equivalence is achieved in the limit $N\to\infty$, but at
finite $N$ the relation remains approximate.,

$\displaystyle f(x)\simeq
f_{N}(x)=\frac{c_{0}}{2}+\sum^{N}_{i=1}c_{i}T_{i}(x)\,,$ (24)

and then, using properties of the $T_{i}(x)$, find that

$\displaystyle\int^{1}_{-1}f(x)dx\approx c_{0}-\sum_{k=1}^{\lfloor
N/2\rfloor}\frac{c_{2k}}{4k^{2}-1}\,.$ (25)

More practical, however, is the writing of the coefficients $c_{i}$ in terms
of the function values at the collocation points $\bar{x}$:

$\displaystyle
c_{i}=\frac{(2-\delta^{i}_{N})}{2N}\left(f(\bar{x}_{0})+(-1)^{i}f(\bar{x}_{N})+2\sum_{j=1}^{N-1}f(\bar{x}_{j})T_{i}(\bar{x}_{j})\right)\,.$
(26)

Combining these expressions, the integral $\mathcal{I}(f)$ can be written as

$\displaystyle\mathcal{I}(f)\approx\mathcal{I}_{N}(f)=\int^{1}_{-1}f(x)dx=f^{T}(\bar{x})\cdot
C(\bar{x})$ (27)

where $f^{T}(\bar{x})$ is a row vector of the function evaluated on the
$\bar{x}$ abscissa and the vector $C$ has entries

$\displaystyle C_{i}=\frac{2}{\kappa_{i}N}\left[1-\sum_{k=1}^{\lfloor
N/2\rfloor}\frac{T_{2k}(\bar{x}_{i})(2-\delta^{N}_{2k})}{4k^{2}-1}\right]\,.$
(28)

When the integrand is a product of functions, each of which is described in
terms of an interpolation over the basis of Chebyshev functions, the integral

$\displaystyle\mathcal{I}(f(x)g(x))=\int^{1}_{-1}f(x)g(x)d\mu$ (29)

with respect to the weight $d\mu=\mu(x)dx$ is a simple extension of the
integral of a single function,

$\displaystyle\mathcal{I}(f(x)g(x)\mu(x))\simeq f^{T}(\bar{x})\cdot
C_{\mu}(\bar{x})\cdot g(\bar{x})$ (30)

where $C_{\mu}(\bar{x})$ is now a diagonal ${(N+1)\times(N+1)}$ matrix with
components

$\displaystyle(C_{\mu})_{ii}=\frac{2\mu(\bar{x}_{i})}{\kappa_{i}N}\left[1-\sum_{k=1}^{\lfloor
N/2\rfloor}\frac{T_{2k}(\bar{x}_{i})(2-\delta^{N}_{2k})}{4k^{2}-1}\right]\,.$
(31)

Thus, an inner product of the form in (20) can be approximated by the discrete
expression

$\langle\phi_{1},\phi_{2}\rangle_{{}_{G}}=(\phi_{1}^{*})^{i}\,G_{ij}\,(\phi_{2})^{j}=\phi_{1}^{*}\cdot
G\cdot\phi_{2},$ (32)

where ${\phi}^{*}=\bar{\phi}^{t}$. The matrix $G$ is known as the _Gram
matrix_ , corresponding to the inner product of interest.

Finally, we address an important technical note regarding integral quadrature
and grid interpolation: integral quadrature over two functions – each of which
are approximated by polynomials of degree $N+1$ – requires a grid of $2N+2$
points to be spectrally accurate boyd2013chebyshev . Therefore, when
evaluating an inner product of the form (32), where $\phi_{1}$ and $\phi_{2}$
are each approximated by $N+1$ Chebyshev polynomials, the Gram matrix $G$
containing the quadrature weights must be of size $(2N+2)\times(2N+2)$.
Practically speaking, we achieve this by evaluating (31) at double the
spectral resolution and then interpolating down to $(N+1)\times(N+1)$ by the
application of an interpolation matrix $\mathfrak{I}$ via
$\mathfrak{I}^{t}\cdot G\cdot\mathfrak{I}$. The form of $\mathfrak{I}$ can be
determined by evaluating the expression for the spectral coefficients in (26)
at the double-resolution collocation points, thereby giving the elements of
the (Chebyshev-Lobatto) interpolation matrix

$\displaystyle\mathfrak{I}_{\bar{i}i}=\frac{1}{\kappa_{i}N}\left(1+\sum_{k=1}^{N}(2-\delta_{k,N})T_{k}(\bar{x}_{\bar{i}})T_{k}(\bar{x}_{i})\right)\,,$
(33)

where $\bar{i}\in[0,2N+1]$.

###  3.3 Reduction of energy norm to the $\ell^{2}$ norm

Once discretised and interpolated by the methods described above, we arrive at
the Gram matrix $G$ for the energy inner product described by (32). Since $G$
will be Hermitian and positive definite, it admits a Cholesky decomposition of
the form $G=F^{*}F$. We can then follow trefethen2005spectra ; Arean:2023ejh
to express the (discretised) energy norm
$\|\cdot\|^{2}_{{}_{G}}=\langle\cdot,\cdot\rangle_{{}_{G}}$ (which we will
refer to as the $G$ norm) as an $\ell^{2}$ norm777For completeness, note that
the $\ell^{2}$ matrix norm (induced by the vector $\ell^{2}$ norm) is given by
the largest singular value of the matrix, i.e. the spectral norm.:

$\displaystyle\|M\|^{2}_{{}_{G}}$
$\displaystyle=\underset{x\in\mathbb{C}^{n}}{\max}\frac{\|Mx\|^{2}_{{}_{G}}}{\|x\|^{2}_{{}_{G}}}=\underset{x\in\mathbb{C}^{n}}{\max}\frac{\langle
Mx,Mx\rangle_{{}_{G}}}{\langle
x,x\rangle_{{}_{G}}}=\underset{x\in\mathbb{C}^{n}}{\max}\frac{x^{*}M^{*}F^{*}FMx}{x^{*}F^{*}Fx}$
$\displaystyle=\underset{y\in\mathbb{C}^{n}}{\max}\frac{y^{*}{F^{-1}}^{*}M^{*}F^{*}FMF^{-1}y}{y^{*}y}=\underset{y\in\mathbb{C}^{n}}{\max}\frac{\|FMF^{-1}y\|^{2}_{{}_{2}}}{\|y\|^{2}_{{}_{2}}}=\|FMF^{-1}\|^{2}_{{}_{2}}\,.$
(34)

When applied to the norm of the resolvent $R_{A,B}(\lambda)\equiv\left(\lambda
B-A\right)^{-1}$ used in equation (11), we can straightforwardly see that

$\left\|\left(\lambda
B-A\right)^{-1}\right\|_{{}_{G}}=\left\|\left(\lambda\tilde{B}-\tilde{A}\right)^{-1}\right\|_{{}_{2}}$
(35)

with $\tilde{B}=FBF^{-1}$, $\tilde{A}=FAF^{-1}$, and thus the
$\epsilon$-pseudospectrum (11) in the $G$ norm $\sigma^{\epsilon}_{{}_{G}}$
can be expressed as the “usual” $\epsilon$-pseudospectrum in the $\ell_{2}$
norm $\sigma^{\epsilon}_{{}_{2}}$,

$\sigma^{\epsilon}_{{}_{G}}(A,B)=\sigma^{\epsilon}_{{}_{2}}(\tilde{A},\tilde{B}).$
(36)

The advantage of this is that it can be shown (see, e.g., trefethen2005spectra
; van1997pseudospectra ) that the $\epsilon$-pseudospectrum can then be
computed from

$\sigma^{\epsilon}_{{}_{2}}(\tilde{A},\tilde{B})=\\{\lambda\in\mathbb{C}:s^{\min}(\lambda\tilde{B}-\tilde{A})<\epsilon\\}$
(37)

where $s^{\min}$ denotes the smallest singular value888Numerically, this is
achieved by using the randomized subspace method outlined in halko2011finding
to compute only the smallest value rather than the entire set.. This is the
definition that we will be using for practical computations of the
pseudospectrum, since it can be readily obtained from the singular value
decomposition of the matrix $M\equiv\lambda FBF^{-1}-FAF^{-1}$.

##  4 Scalar, electromagnetic and gravitational perturbations for planar AdS-
RN

Let us now consider our physical setup, consisting of Einstein-Maxwell theory
in $D=5$ spacetime dimensions with action given by

$\displaystyle S=\frac{1}{16\pi G_{5}}\int
dx^{5}\sqrt{-g}\left(R-2\Lambda-\frac{1}{4}F_{\mu\nu}\,F^{\mu\nu}\right)$ (38)

where $F\equiv dA$ and $\Lambda=-6$ is the cosmological constant. The
variation of the action (38) gives rise to the following field equations

$\displaystyle
R_{\mu\nu}-\frac{R}{2}g_{\mu\nu}-6g_{\mu\nu}-\frac{1}{2}\left(F_{\mu\rho}F_{\nu}{}^{\rho}-\frac{1}{4}\,g_{\mu\nu}F_{\lambda\rho}F^{\lambda\rho}\right)=0\,,$
(39) $\displaystyle\nabla_{\mu}F^{\mu\nu}=0\,.$

The above equations admit what is known as the AdS Reissner–Nordström (RN)
black-brane solution. Concretely, the metric and gauge field, in Poincaré
coordinates, take the form

$\displaystyle ds^{2}$
$\displaystyle=-r^{2}f(r)dt^{2}+\frac{dr^{2}}{r^{2}f(r)}+r^{2}\left(dx_{1}^{2}+dx_{2}^{2}+dx_{3}^{2}\right)$
(40) $\displaystyle A$
$\displaystyle=\mu\left(1-\frac{r_{\mathrm{h}}^{2}}{r^{2}}\right)dt$
$\displaystyle f$
$\displaystyle=1-\frac{r_{\mathrm{h}}^{4}}{r^{4}}(1+Q^{2})+\frac{r_{\mathrm{h}}^{6}}{r^{6}}Q^{2}\,,$

where $r_{\mathrm{h}}$ is the black hole horizon, $Q$ is the electric charge
carried by the black hole and $\mu=\sqrt{3}\,r_{\mathrm{h}}\,Q$ is the
chemical potential. Without loss of generality, we will henceforth fix units
where the horizon radius is at $r_{\mathrm{h}}=1$ – the explicit dependence
can be reinstated through

$r\to\frac{r}{r_{\mathrm{h}}},\qquad t\to t\,r_{\mathrm{h}},\qquad x_{i}\to
x_{i}\,r_{\mathrm{h}},\qquad\mu\to\frac{\mu}{r_{\mathrm{h}}}\,.$ (41)

We can now perform the following coordinate transformation

$\displaystyle t=v-r_{*}\,,\qquad Y=r_{*},$ (42)

where $dY=dr_{*}=\tfrac{dr}{r^{2}f(r)}$. For convenience, we also compactify
the radial coordinate $r$ by introducing

$z=\frac{1}{r}$

so that our physical domain spans $z=0$ (the AdS boundary, $r\to\infty$) to
$z=1$ (the black hole horizon, $r=1$). To fix notation, we will refer to this
coordinate system $(v,z)$ as the ingoing Eddington–Finkelstein coordinates.

Overall, we obtain

$\displaystyle ds^{2}$
$\displaystyle=\frac{1}{z^{2}}\left(-f(z)dv^{2}-2dv\,dz+dx_{1}^{2}+dx_{2}^{2}+dx_{3}^{2}\right),$
(43) $\displaystyle A$
$\displaystyle=\mu(1-z^{2})\,\left(dv+\frac{dz}{f(z)}\right)\,,$

where

$f(z)=1-z^{4}\left(1+Q^{2}\right)+z^{6}Q^{2}.$

The (dimensionless) temperature of the black hole is given by

$\displaystyle\frac{T}{\mu}=\frac{1}{2\pi\,\mu}\left[2-Q^{2}\right]\,.$ (44)

When $Q=0$ the system is neutral and the bulk solution (43) reduces to the
AdS-Schwarzschild black brane, while $Q=Q_{\mathrm{max}}=\sqrt{2}$ corresponds
to the extremal limit $T/\mu\to 0$.

###  4.1 Massless neutral scalar perturbations

We now consider a neutral, massless scalar field $\Phi$, in the background of
the 5-dimensional planar AdS-RN solution described above. Such a perturbation
obeys the Klein-Gordon equation

$\square\Phi=0\,,$ (45)

which can then be tackled in different coordinate systems.

##### Poincaré coordinates:

It is not difficult to show that, in Poincaré coordinates $(t,r)$, the above
equation takes the form (13). To see that, we decompose the scalar as
Horowitz:1999jd

$\Phi=\frac{1}{r^{\frac{3}{2}}}\phi(t,r)\,e^{i\vec{k}\cdot\vec{x}}\,$

and we substitute in equation (45) to find

$\displaystyle\left(-\partial^{2}_{t}+\partial^{2}_{r_{*}}-V(r)\right)\phi(t,r)=0,$
(46) $\displaystyle V(r)\equiv
f(r)\left(k^{2}+\frac{3}{4}f(r)+\frac{3}{2}r\frac{df(r)}{dr}\right)$

and $\partial_{r_{*}}=r^{2}f(r)\partial_{r}$. These coordinates are not
horizon-penetrating, however, and are therefore unsuitable to perform the
pseudospectrum analysis. This fact has been noted also in Arean:2023ejh .

##### Ingoing Eddington–Finkelstein coordinates:

We consider instead ingoing Eddington–Finkelstein coordinates, which are
regular at the black hole horizon and also allow one to reach the AdS
boundary. With these coordinates, one can naturally impose the required
boundary conditions at the two ends of the physical domain following surfaces
of constant $v$. Note that this is not possible in asymptotically flat
spacetimes, where – due to the different asymptotics – surfaces of constant
$v$ lead one to past null infinity $\mathscr{I}^{-}$ instead. Indeed, that is
the major motivation for using the hyperboloidal coordinates originally
developed in Jaramillo:2020tuu  and used thereafter in the literature.999We
thank Rodrigo Panosso Macedo for clarification of these points. See e.g.
reference PanossoMacedo:2023qzp  for an overview of hyperboloidal coordinates
and their applications. In reference Arean:2023ejh , the authors worked with
“compactified regular coordinates” – which is more similar to the
hyperboloidal coordinates used in asymptotically flat cases. We find this to
be unnecessary, and prefer to keep the elegance and simplicity of the ingoing
Eddington–Finkelstein system, at the expense of dealing with a generalised
eigenvalue problem (as opposed to a standard eigenvalue problem). As discussed
in section 2.1, this poses no problems in practice.

With the ansatz

$\Phi=z^{\frac{3}{2}}\phi(v,z)\,e^{i\vec{k}\cdot\vec{x}}\,,$

the Klein-Gordon equation (45) becomes

$\displaystyle\left[-2\,\partial_{vz}+\partial_{z}\left(f(z)\partial_{z}\right)-\hat{V}(z)\right]\phi(v,z)=0,$
(47)
$\displaystyle\hat{V}(z)\equiv\frac{V(z)}{f(z)}=k^{2}+\frac{15}{4}\frac{f(z)}{z^{2}}-\frac{3}{2}\frac{f^{\prime}(z)}{z}=k^{2}+\frac{15}{4z^{2}}+\frac{9}{4}\left(1+Q^{2}\right)z^{2}-\frac{21}{4}Q^{2}z^{4}$

where the prime ′ denotes derivative with respect to $z$.

Since solutions to equation (47) behave as $\phi\sim z^{5/2}$ close to the AdS
boundary $z=0$, we proceed to redefine

$\phi(v,z)=z^{\frac{5}{2}}\psi(v,z).$ (48)

With this redefinition, the required boundary conditions reduce to demanding
regularity of the fields at the two boundaries. The equation of motion now
takes the form

$\displaystyle\left[-2\,\partial_{vz}-\frac{5}{z}\partial_{v}+\partial_{z}\left(f(z)\partial_{z}\right)+\frac{5f(z)}{z}\partial_{z}-\hat{U}(z)\right]\psi(v,z)=0,$
(49)
$\displaystyle\hat{U}(z)\equiv\frac{U(z)}{f(z)}=k^{2}-4\frac{f^{\prime}(z)}{z}=k^{2}+8z^{2}(26(1+Q^{2})-3Q^{2}z^{2})\,.$

The latter can be rewritten in the form of a generalised eigenvalue problem

$L_{2}\,\partial_{v}\psi=iL_{1}\psi$ (50)

where

$\displaystyle iL_{1}$
$\displaystyle=\left[\partial_{z}\left(zf(z)\partial_{z}\right)-\left(k^{2}z-4f^{\prime}(z)\right)\mathbb{1}\right]+4f(z)\partial_{z},$
(51) $\displaystyle L_{2}$ $\displaystyle=2z\partial_{z}+5\mathbb{1}.$ (52)

Note that the $iL_{1}$ operator contains the full physics of the problem (thus
justifying our definition (12) for the pseudospectrum) and that the term in
square brackets has the structure of a Sturm–Liouville operator.

Equation (19) now takes the form

$\displaystyle E[\psi]=\frac{1}{2}\int_{0}^{1}dz\bigg{[}$ $\displaystyle
z^{5}f(z)\partial_{z}\bar{\psi}\,\partial_{z}\psi+\frac{5}{2}z^{4}f(z)\left(\bar{\psi}\,\partial_{z}\psi+\partial_{z}\bar{\psi}\,\psi\right)$
$\displaystyle{}+z^{3}\left(10f(z)-\frac{3}{2}zf^{\prime}(z)+k^{2}z^{2}\right)\bar{\psi}\,\psi\bigg{]}.$
(53)

We note that every term in $L_{1}$, $L_{2}$ and $E[\psi]$ is explicitly
regular for $z\in[0,1]$. This contrasts with the approach of Arean:2023ejh ,
where divergent terms at $z=0$ need special care. We can now define the scalar
product

$\displaystyle\langle\psi_{1},\psi_{2}\rangle_{{}_{E}}=\frac{1}{2}\int_{0}^{1}dz\bigg{[}$
$\displaystyle
z^{5}f(z)\,\partial_{z}\bar{\psi}_{1}\,\partial_{z}\psi_{2}+\frac{5}{2}z^{4}f(z)\left(\bar{\psi}_{1}\,\partial_{z}\psi_{2}+\partial_{z}\bar{\psi}_{1}\,\psi_{2}\right)$
$\displaystyle{}+z^{3}\left(10f(z)-\frac{3}{2}zf^{\prime}(z)+k^{2}z^{2}\right)\bar{\psi}_{1}\,\psi_{2}\bigg{]}$
(54)

as well as the energy norm
$\|\cdot\|^{2}_{{}_{E}}=\langle\cdot,\cdot\rangle_{{}_{E}}$. Following section
3.2, we can also build its discretised counterpart

$\langle\psi_{1},\psi_{2}\rangle_{{}_{G}}=\psi_{1}^{*}\cdot
G^{E}\cdot\psi_{2},$ (55)

with the Gram matrix

$G^{E}=\frac{1}{2}\left(C_{V}+C_{1}\cdot\mathbb{D}+\mathbb{D}^{t}\cdot
C_{1}+\mathbb{D}^{t}\cdot C_{2}\cdot\mathbb{D}\right)$ (56)

and $\left(C_{a}\right)_{ij}=\mu_{a}(z_{i})\,W_{i}\,\delta_{ij}$ (no summing
on $i$). $W_{i}$ are the quadrature weights introduced in section 3.2,
$\delta_{ij}$ is the Kronecker delta, and

$\displaystyle\mu_{V}(z)=z^{3}\left(10f(z)-\frac{3}{2}zf^{\prime}(z)+k^{2}z^{2}\right),\qquad\mu_{1}(z)=\frac{5}{2}z^{4}f(z),\qquad\mu_{2}(z)=z^{5}f(z).$

####  4.1.1 Results

With the operators $L_{1}$, $L_{2}$ and the Gram matrix $G^{E}$ just defined,
we are now ready to compute the $\epsilon$-pseudospectrum $\sigma^{\epsilon}$
of this system following the approach outlined in section 3, in particular
equation (37).

(a) $k=0$, $Q=0$.

(b) $k=10$, $Q=0$.

(c) $k=10$, $Q=0.1$.

(d) $k=10$, $Q=0.99Q_{\mathrm{max}}$.

Figure 1: The $\epsilon$-pseudospectrum $\sigma^{\epsilon}$ for various planar
black hole solutions coupled to a massless scalar field. Recall that the
extremal solution is obtained for $Q_{\mathrm{max}}=\sqrt{2}$. Above each
plot, the condition numbers $\kappa$ of the eigenvalues are plotted with
increasing overtone number. The red crosses mark the spectrum points.

Our results are displayed in figure 1 for a choice of $Q,k$, while a zoomed-in
and a zoomed-out plot of the pseudospectrum for the uncharged $k=0$ case are
shown in figure 2. The crosses correspond to the locations of the QNMs, where
$\epsilon=0$. In figure 1 we also display the condition numbers101010Except
for figure 1(d) where $\kappa_{0},\kappa_{1}$ are the condition numbers for
the fundamental mode, $\kappa_{2}$ is for the eigenvalue with $\Re(\omega)=0$,
and $\kappa_{i},~{}i\geq 3$ are for the overtones. for each case, as defined
in section 2. Recall that the condition number gives an indication of the
relative size of the change in an eigenvalue with respect to the magnitude of
a perturbation: quantitatively, for a completely stable spectrum $\kappa=1$
for all modes; for non-conservative systems typically the condition number
increases with increasing overtone, while for severely unstable spectra
$\kappa\gg 1$ for all modes. Thus, the condition number can be used along with
the pseudospectrum to understand the degree of instability of a spectrum.

Figure 2: Zoomed-in and zoomed-out view of the $\epsilon$-pseudospectrum
$\sigma^{\epsilon}$ for the $k=0$ case of the $Q=0$ black hole solution from
equation (43). The contour lines correspond to level sets of the quantity
$\epsilon$, as defined in equation (36). _Left:_ a wide view of the
pseudospectrum demonstrates the flat behaviour of contour lines away from the
quasinormal modes. _Right:_ a zoomed view of the pseudospectrum near the
fundamental mode $\omega_{0}=3.119-2.747i$ shows that small perturbations
around this mode remain close to the $\omega_{0}$ value.

As can be observed from these figures, although the pseudospectra still form
circular sets if we zoom arbitrarily close to the spectrum, its large scale
global structure present open sets. These extend into large regions of the
complex plane even for small $\epsilon$, indicating spectral instability. This
property implies that small-scale perturbations can lead to perturbed spectra
which migrate into regions that are much further away than $\epsilon$.
Interpreted through the gauge/gravity correspondence, this indicates spectral
instability of the dual strongly coupled system.

We also observe that modes that reside further from the real axis are more
unstable; for the dual theory, this implies that short-lived excitation are
more unstable. Given this observation, we find that as the momentum $k$
increases and the modes move closer to the real axis, their stability
increases. In a similar spirit, as the charge $Q$ increases, the modes move
marginally further away from the real axis and thus their stability marginally
decreases.

In comparing our results for $Q=0$ with those of Arean:2023ejh , we see
qualitative agreement. There, the computation was performed in a different
coordinate system, namely using “regular compact coordinates”, while herein we
use ingoing Eddington–Finkelstein coordinates. Nevertheless, the computed
pseudospectrum has the same qualitative features, indicating that the
pseudospectrum is “spacetime slicing”-independent. Note, however, even modulo
the overall scale which depends on the normalisation, we do not expect to see
a precise quantitative agreement since we are using a different definition for
the pseudospectrum.

Finally, in the near-extremal limit we observed the existence of a series of
equally-spaced quasinormal modes with zero real part, i.e. non-oscillatory
modes. In figure 1(d) we can see the first of these modes. Unlike the spectrum
of the vector and scalar channel of gravitoelectric perturbations (see for
example figure 4(d)), where non-oscillatory modes collect from $-i\infty$ to
form a branch cut, these modes do not accumulate as we approach extremality.
Curiously, they remain stable against both infrared and ultraviolet
perturbations, as we will show in the next section.

####  4.1.2 Testing the stability of scalar QNMs

While the pseudospectrum allows us to visualize the stability of the QNM
spectrum to general perturbations of order $\epsilon$, we may also be
interested in how _specific_ perturbations to the background affect the
spectrum Cheung:2021bol . These perturbations occur in the potential and can
be of any form in general, although perturbations informed by physical
processes are particularly important. For example, changes to the potential
due to deviations from pure Schwarzschild geometry in the near-horizon region
are expected to manifest as changes in the first few overtones
Konoplya:2022pbc . With this in mind, we consider changes in the potential in
(49) of the form $\hat{U}\to\hat{U}+\delta\hat{U}$ with
$\delta\hat{U}(z)=\epsilon\cos(\pi\omega(z-1))$. In terms of the generalised
eigenvalue problem

$\left(L_{1}-\lambda L_{2}\right)\psi=0,$

this corresponds to a perturbation of the form $L_{1}\to L_{1}+\delta\hat{U}$.

It is important to note that adding a perturbation of size
$\|\delta\hat{U}\|_{{}_{E}}\sim\epsilon$ relies on the definition of a norm
just as the pseudospectrum does. Thus, to ensure that the size of
$\delta\hat{U}$ remains of order $\epsilon$ relative to $L_{1}$ we must
convert the energy norm to the 2-norm via the procedure outlined in subsection
3.3, normalise $\|\delta\hat{U}\|_{{}_{2}}$, and finally add the overall
perturbation to $L_{1}$.

In figure 3, we show the effects of adding a sinusoidal perturbation of
magnitude $\epsilon=10^{-8}$ to the potential. In each plot, the unperturbed
spectrum is denoted by red dots; the spectra resulting from ${L_{1}\to
L_{1}+\delta\hat{U}}$ are overlaid in circles, crosses, and triangles. Contour
lines of the $\epsilon$-pseudospectra are included for reference, with the
magnitudes indicated in the plot. In all cases, higher overtones are more
prone to instability than the fundamental modes or lower overtones. As the
perturbation frequency as increased, lower overtones become shifted. Note that
the shifted spectrum follows the contours of the $\epsilon$-psuedospectrum.

(a) $k=0,Q=0$.

(b) $k=10,Q=0$.

(c) $k=10,Q=0.1$.

(d) $k=10,Q=0.99\,Q_{\mathrm{max}}$.

Figure 3: The perturbed spectrum of QNMs for the massless scalar field for
various AdS-RN configurations (c.f. spectra in figure 1). Recall that the
extremal solution is obtained for $Q_{\mathrm{max}}=\sqrt{2}$. The
perturbation in the potential is taken to be of the form
$\delta\hat{U}=\epsilon\,\cos(\pi\omega(z-1))$ with $\epsilon=10^{-8}$. Red
dots mark the spectrum points.

###  4.2 Metric and gauge field perturbations

We now move on to study gravitational and electromagnetic fluctuations around
the background configuration (43). Gravitational QNMs for the AdS-
Schwarzschild and AdS-RN black branes were discussed previously in the
literature Policastro:2002se ; Kovtun:2005ev . In particular, for non-zero
charge the electromagnetic and gravitational perturbations couple to each
other. However, it is well known that these split into 3 sectors, depending on
how they transform under the rotation group transverse to the momentum of the
fluctuations: the tensor (for $D>4$), the vector and the scalar sectors.

In what follows we will use the formalism of master equations, which was first
derived for AdS-RN in Kodama:2003kk  and later generalised in Jansen:2019wag .
The idea is to construct gauge invariant combinations of the perturbations,
and then combine them into master scalars that satisfy Klein-Gordon equations
with certain potentials. The tensor sector contains a single master scalar
$\Phi_{2}$. The vector sector contains two, whose equations can be decoupled
by taking certain linear combinations that we denote $\Phi_{1\pm}$: these
correspond to the shear channel (positive sign) and the transverse gauge field
(negative sign). Finally, the scalar sector contains two more master scalars
that can be decoupled, which we denote $\Phi_{0\pm}$: these correspond to the
sound channel (positive sign) and the diffusion channel (negative sign). Out
of these 5 types of fluctuations only the shear, sound and diffusion contain
hydrodynamic modes, meaning modes that obey dispersion relations such that the
frequency approaches zero as the momentum is decreased.

Each of these master scalars satisfy a master equation of the form

$\displaystyle\square\Phi_{i}-Z_{i}(z)\Phi_{i}$ $\displaystyle=0\,,$ (57)

where $\square$ is the Laplacian on the spacetime (43). Without loss of
generality, we decompose the master scalars into plane waves with the momentum
pointing in the $x_{1}$ direction

$\Phi_{i}(v,z,x_{1})=e^{i\,kx_{1}}\Phi_{i}(v,z)\,,$ (58)

where $k$ is the momentum in the $x_{1}$ direction.

The potentials for the tensor, vector and scalar sectors can be written
respectively as

$\displaystyle Z_{2}(z)$ $\displaystyle=0$ (59) $\displaystyle Z_{1\pm}(z)$
$\displaystyle=\frac{1}{2}(Z^{(1)}_{11}+Z^{(1)}_{22})\pm\frac{1}{2}\sqrt{(Z^{(1)}_{11}-Z^{(1)}_{22})^{2}+4(Z^{(1)}_{12})^{2}}$
$\displaystyle Z_{0\pm}(z)$
$\displaystyle=\frac{1}{2}(Z^{(0)}_{11}+Z^{(0)}_{22})\pm\frac{1}{2}\sqrt{(Z^{(0)}_{11}-Z^{(0)}_{22})^{2}+4(Z^{(0)}_{12})^{2}}$

where

$\displaystyle Z^{(1)}_{11}=zf^{\prime}-3f+z^{4}a^{\prime 2}$ (60)
$\displaystyle Z^{(1)}_{22}=3zf^{\prime}-3f$ $\displaystyle
Z^{(1)}_{12}=|k|z^{3}a^{\prime}$

and

$\displaystyle Z^{(0)}_{11}=\frac{a^{\prime
2}z^{2}}{\mathcal{D}^{2}}\left(18fa^{\prime
2}z^{2}-9\mathcal{F}\left(2f+zf^{\prime}\right)+36k^{2}f+4k^{4}z^{2}\right)-2(2f-zf^{\prime})$
(61) $\displaystyle
Z^{(0)}_{22}=\frac{2}{\mathcal{D}^{2}}\left(12k^{2}z^{2}fa^{\prime
2}+12k^{2}\mathcal{F}\left(nf-
zf^{\prime}\right)-8k^{4}\left(f-zf^{\prime}\right)\right)$ $\displaystyle
Z^{(0)}_{12}=\frac{|k|za^{\prime}}{\mathcal{D}^{2}}\frac{2}{\sqrt{3}}\left(18z^{2}fa^{\prime
2}+9\mathcal{F}(z)(2f-zf^{\prime})+12k^{2}f+4k^{4}z^{2}\right)$
$\displaystyle\mathcal{F}=f^{\prime}/z\,,\quad\mathcal{D}=2k^{2}-3\mathcal{F}\,.$

Note that the tensor sector corresponds to a massless scalar with no
potential, which is precisely what we studied in the previous subsection.
Hence, in what follows we will focus only on the vector and scalar sectors. It
is also important to emphasise that the two channels within the vector and the
scalar sectors decouple from each other precisely because the following
combination is constant

$\displaystyle\frac{Z^{(i)}_{11}-Z^{(i)}_{2}}{Z^{(i)}_{12}}=-c_{i}\frac{1+Q^{2}}{|k|Q},\qquad
i=0,1,$ (62)

with $c_{0}=2,c_{1}=\tfrac{4}{\sqrt{3}}$.

Following the same methodology as for the scalar field perturbations, we
consider the following rescalings. First, we rescale the master scalars as

$\displaystyle\Phi_{i}(v,z)=z^{\frac{3}{2}}\phi_{i}(v,z)$ (63)

to bring the equations into the form (47)

$\displaystyle-2f\partial_{vz}\phi_{i}+f\partial_{z}(f\partial_{z}\phi_{i})-V_{i}(z)\phi_{i}=0$
(64) $\displaystyle V_{i}(z)\equiv
f\hat{V}_{i}(z)=\frac{f}{z^{2}}\,(Z_{i}(z)+k^{2}z^{2})+\frac{3}{2}\,\frac{f}{z^{2}}\left(-zf^{\prime}+\frac{5}{2}f\right).$

However, as before, for the purposes of the numerical implementation we also
consider the redefinition

$\displaystyle\Phi_{i}(v,z)=z^{\Delta_{i}}\psi_{i}(v,z)\quad\text{ where
}\quad\Delta_{2}=4,\quad\Delta_{1\pm}=3,\quad\Delta_{0\pm}=2$ (65)

so that, given the absence of a source, the leading UV behaviour is
$\psi=\mathrm{VEV}+\mathcal{O}(z)$. In this case the equations for the master
scalars become

$\displaystyle-2f\partial_{vz}\psi_{i}-(2\Delta_{i}-3)\frac{f}{z}\partial_{v}\psi+f\partial_{z}(f\partial_{z}\psi_{i})+(2\Delta_{i}-3)\frac{f^{2}}{z}\partial_{z}\psi-
U_{i}(z)\psi_{i}=0$ (66) $\displaystyle U_{i}(z)\equiv
f\hat{U_{i}}(z)=\frac{f}{z^{2}}(Z_{i}(z)+k^{2}z^{2})+\Delta_{i}\,\frac{f}{z^{2}}\left(-zf^{\prime}+(4-\Delta_{i})f\right).$

The next step is to bring the above equation in the form

$L_{2}\,\partial_{v}\psi_{i}=iL_{1}\psi_{i}$ (67)

where

$\displaystyle iL_{1}$
$\displaystyle=zf\partial^{2}_{z}+zf^{\prime}\partial_{z}+(2\Delta_{i}-3)f\partial_{z}-z\hat{U_{i}}(z)$
(68) $\displaystyle L_{2}$
$\displaystyle=2z\partial_{z}+(2\Delta_{i}-3)\mathbb{1}\,.$ (69)

We again note that the tensor channel corresponds to a massless scalar and one
can easily check that the above expressions agree with those in the previous
subsection.

For the energy norm, our starting point is equation (19). After performing the
redefinition $\phi_{i}=z^{\Delta_{i}-3/2}\,\psi_{i}$, we obtain the following
expression in terms of $\psi_{i}$,

$E[\psi_{i}]=\int
dz\frac{z^{2\Delta_{i}-3}}{4}\left(\left(\hat{V_{i}}+(2\Delta_{i}-3)^{2}\frac{f}{4z^{2}}\right)\psi_{i}^{2}+\frac{f}{z}(2\Delta_{i}-3)\psi_{i}\partial_{z}\psi_{i}+f\partial_{z}\psi_{i}^{2}\right)$
(70)

We can now, in a matter entirely analogous to that of the previous subsection,
extract our Gram matrix $G^{E}$ and compute the pseudospectrum.

####  4.2.1 Results

(a) $k=0$, $Q=0$.

(b) $k=0$, $Q=0.1$.

(c) $k=1$, $Q=0.1$.

(d) $k=1$, $Q=1=0.7\,Q_{\mathrm{max}}$.

Figure 4: The $\epsilon$-pseudospectrum for various values of momentum $k$ and
charge $Q$ in the vector channel, shear sector. Recall that the extremal
solution is obtained for $Q_{\mathrm{max}}=\sqrt{2}$. Above each plot, the
condition numbers $\kappa$ of the eigenvalues are plotted with increasing
overtone number. The red crosses mark the spectrum points. Figure 5: Zoomed-in
plot of the $\epsilon$-pseudospectrum in the shear sector for $Q=0$ and $k=0$.
We see the spectral lines crossing in the upper half plane indicating
transient instability. The red crosses mark the spectrum points.

(a) $k=0$, $Q=0$.

(b) $k=1$, $Q=0$.

(c) $k=1$, $Q=0.1$.

(d) $k=1$, $Q=1=0.7\,Q_{\mathrm{max}}$.

Figure 6: The $\epsilon$-pseudospectrum for various values of momentum $k$ and
charge $Q$ in the vector channel, transverse gauge sector. Recall that the
extremal solution is obtained for $Q_{\mathrm{max}}=\sqrt{2}$. Above each
plot, the condition numbers $\kappa$ of the eigenvalues are plotted with
increasing overtone number. The red crosses mark the spectrum points.

(a) $k=1$, $Q=0$.

(b) $k=1.5$, $Q=0$.

(c) $k=1$, $Q=0.8$.

(d) $k=1$, $Q=1=0.7\,Q_{\mathrm{max}}$.

Figure 7: The $\epsilon$-pseudospectrum for various values of momentum $k$ and
charge $Q$ in the scalar channel, charge diffusion sector. Recall that the
extremal solution is obtained for $Q_{\mathrm{max}}=\sqrt{2}$. Above each
plot, the condition numbers $\kappa$ of the eigenvalues are plotted with
increasing overtone number. The red crosses mark the spectrum points.

(a) $k=1$, $Q=0$.

(b) $k=5$, $Q=0$.

(c) $k=1$, $Q=0.5$.

(d) $k=1$, $Q=1=0.7\,Q_{\mathrm{max}}$.

Figure 8: The $\epsilon$-pseudospectrum for various values of momentum $k$ and
charge $Q$ in the scalar channel, sound sector. Recall that the extremal
solution is obtained for $Q_{\mathrm{max}}=\sqrt{2}$. Above each plot, the
condition numbers $\kappa$ of the eigenvalues are plotted with increasing
overtone number. The red crosses mark the spectrum points.

The pseudospectrum for the shear, transverse gauge, charge diffusion and sound
channels are shown respectively in figures 4, 5, 6, 7 and 8 for different
values of the momentum $k$ and the charge $Q$. The red crosses correspond to
the locations of the QNMs for the respective sector. Just like for the scalar
field discussed in section 4.1, the pseudospectra form circular sets if we
zoom arbitrarily close to the spectrum, but their large scale global structure
present open sets signaling spectral instability in all channels. We also
observe that instability increases as the imaginary part of the QNM frequency
increases, with implications for the various modes as the momentum $k$ and the
charge $Q$ are tuned. From the dual field theory side, this can be interpreted
as shorter-lived excitations being more unstable.

With the exception of the transverse gauge, all other sectors contain
hydrodynamic modes, namely modes that reside in the origin of the complex
frequency plane for zero momentum, i.e. $\omega\to 0$ as $k\to 0$. A zoomed-in
version of the pseudospectrum for the shear mode for $Q=0,k=0$ is shown in
figure 5 – similar behaviour appears for hydrodynamic sound and charge
diffusion modes. We observe that the pseudospectral contour levels cross to
the upper half complex plane, which translates to an unstable perturbed
spectrum. A similar behaviour was observed in 4-dimensional asymptotically
flat RN black hole in the extremal limit Destounis:2021lum , where a zero-
frequency QNM exists associated to the Aretakis instability Aretakis:2011hc .
In general, the existence of a zero-modes is of key importance to transient
instabilities that are typically resolved with pseudospectra: see for example
the discussion of the transition to turbulence in hydrodynamics treften90hydro
.

Particular attention was also given to the extremal limit $T/\mu\to 0$ (or
equivalently $Q\to Q_{\mathrm{max}}$, where $Q_{\mathrm{max}}=\sqrt{2}$). In
terms of the spectrum for the scalar and vector channels, we see an
accumulation of poles along the negative imaginary axis that moved upwards
from $-i\infty$; the consensus is that these modes will give rise to a branch
cut at $T=0$ Edalati:2010hk ; Edalati:2010pn . We see that the pseudospectrum
is deformed in that region to reflect the existence of these new QNMs, while
its qualitative behaviour away from the negative imaginary axis remains mostly
unchanged; note, in particular, that the asymptotic behaviour of the contour
lines at large real part of the frequency stays the same. No other significant
features appear.

The pseudospectrum in the extremal limit of asymptotically flat RN black hole
was studied in Destounis:2021lum . A significant difference there is that the
spectrum contains a branch cut at all values of the temperature and at
extremality a zero-frequency mode, associated to the Aretakis instability
Aretakis:2011hc , appears. It was observed that the pseudospectral contour
levels cross to the upper half plane around the origin of the complex
frequency plane, where the Aretakis mode lives. The extremal limit of Kerr was
considered in Yang:2022wlm  for energetically infinitesimal perturbations
finding instability of the zero-damping modes. Interestingly, Yang:2022wlm
found stability of zero-damping modes in near-extremal RN-de Sitter black
holes under energetically infinitesimal perturbations.

####  4.2.2 Testing the stability of gravitoelectric QNMs

(a) Sound sector: $k=1$, $Q=0$.

(b) Sound sector: $k=1$, $Q=0.5$.

(c) Charge diffusion sector: $k=1$, $Q=0$.

(d) Charge diffusion sector: $k=1$, $Q=0.8$.

Figure 9: The perturbed spectrum of QNMs in the scalar sector, for both the
sound and charge diffusion channels for various AdS-RN configurations. The
perturbation in the potential is taken to be of the form
${\delta\hat{U}=\epsilon\,\cos(\pi\omega(z-1))}$ with $\epsilon=10^{-8}$. Red
dots mark the unperturbed spectrum points.

Just like the scalar case discussed in the previous subsection, we consider
the response of the spectrum under specific perturbations in the potential. In
figure 9, just like for the scalar field case, we show the effects of adding a
sinusoidal perturbation of magnitude $\epsilon=10^{-8}$ to the potential in
the sound and charge diffusion sectors (scalar channel). We see that the
hydrodynamic mode remains stable within this class of perturbations as the
perturbation frequency increases up to and including $\omega=100$ for all
values of $k,Q$ that we have considered. Note that in all cases the
hydrodynamic modes have smaller condition numbers compared to the non-
hydrodynamic modes, hence their increased stability is to be expected. The
non-hydrodynamic overtones are progressively more sensitive as the frequency
and the overtone number increase. Note that, once again, that the shifted QNMs
follow the contours of the $\epsilon$-pseudospectrum.

###  4.3 Universality of the pseudospectrum

In this subsection we comment on the behaviour of pseudospectral contour
levels at large real values of the spectral frequency. In these regions, the
pseudospectral levels are known to belong to one of four “universal” classes
Zworski

$\Im(\omega)\sim F(\Re(\omega))\,,\qquad\Re(\omega)>>1$ (71)

where

$F(x)=\begin{cases}(i)\quad e^{ax}\,,\;a>0\\\ (ii)\quad C\\\ (iii)\quad
B\,\ln(x)\\\
(iv)\quad\gamma\,x^{b}\,,\;\gamma>0\,,\;b\in\mathbb{R}\end{cases}$

where $a,b,\gamma,B,C$ are constant and are controlled by the qualitative
properties of the underlying system.

Figure 10 displays the asymptotic structure of the pseudospectrum, for the
scalar and vector channels for $Q=0.5,k=1$. We see that the pseudospectral
levels coincide with each other at large $\Re(\omega)$, demonstrating
asymptotic universality. This observation is in agreement with the analysis
done in 4-dimensional asymptotically flat spacetime for the Schwarzschild
Jaramillo:2020tuu ; Jaramillo:2021tmt  and RN black holes Destounis:2021lum ,
where universality was also observed.

Figure 10: Overlay of asymptotic pseudospectral contour levels for the
$Q=0.5,k=1$ AdS-RN black hole. The levels correspond to
$\sigma^{\epsilon}=\\{100,1,0.005,0.0005,0.00001\\}$ level curves. _Left:_ For
the sound (red) and charge diffusion (black) sectors. _Right:_ For the shear
(blue) and transverse gauge (green) sectors.

Given the universality observed, it is natural to try to fit the
pseudospectrum contour line in an attempt to extract this asymptotic
behaviour. Typical behaviours in our setting belong to case (iii) or (iv)
Zworski , with logarithmic asymptotics (case (iii)) appearing when considering
potentials and/or boundary conditions allowing low regularity and polynomial
asymptotics (case (iv)) in settings with enhanced regularity.

Analysis of the asymptotic structure of contour levels in 4-dimensional
asymptotically flat spacetime for both Schwarzschild Jaramillo:2020tuu  and RN
Destounis:2021lum  revealed a logarithmic function, namely case (iii). A
preliminary analysis of the contour lines in our case seems to hint towards a
polynomial behaviour (case (iv))

$\Im(\omega)\sim C_{1}+C_{2}\Re(\omega)^{\alpha}\,$ (72)

where $C_{1},C_{2},\alpha$ are constants, with $\alpha>0$. Note that this
behaviour seems to persist even in the extremal limit, and also using a
different coordinate system.111111We have checked this using the “compactified
regular coordinates” of Arean:2023ejh . More analysis is required to reach a
stronger conclusion.

##  5 Conclusions and outlook

In this work we have studied the spectral stability of AdS RN black branes as
a function of momenta and charge (or equivalently, $T/\mu$) by means of the
pseudospectrum trefethen2005spectra . We considered perturbations
corresponding to a neutral, massless scalar field as well as gravitoelectric
perturbations in all sectors Kovtun:2005ev ; Jansen:2019wag . This
calculation, other than shedding light on the properties of AdS black branes,
is also relevant for strongly coupled systems via the gauge/gravity
correspondence. Within this context, QNMs are interpreted as the poles of the
retarded Green’s function of the dual strongly coupled theory Kovtun:2005ev ;
Horowitz:1999jd , and they have been pivotal in identifying phase transitions
holographically. The spectral instability found herein implies that the
spectrum of the dual strongly coupled field theory is also unstable under
perturbations and that one should focus more on correlation functions in
position-space rather than momentum-space.

For all the perturbations considered we find spectral instability, including
the channels that contain hydrodynamic modes. This agrees and naturally
extends the results of Arean:2023ejh . For the hydrodynamic modes in
particular, we find that the pseudospectral contour lines cross to the upper
half plane (when $k$ is small enough) indicating transient instabilities. Such
instabilities are in fact known to arise in hydrodynamics and it has been
argued that they explain the early onset of turbulence treften90hydro .

We have also investigated the asymptotic behaviour of pseudospectral contours
and we find universality across different perturbation sectors. This extends
the universality observed for asymptotically flat black holes
Jaramillo:2020tuu ; Destounis:2021lum . Some preliminary analysis indicated
that our curves may be fitted better by a polynomial at large real frequency,
rather than a logarithmic function as in asymptotic behaviour; further
analysis is required in order to reach a solid conclusion. The difference
might be due to the enhanced regularity of the computation in AdS Zworski .

On a technical note, we have carried out our computation making use of ingoing
Eddington-Finkelstein coordinates. The equation of motion for the
perturbations become first order in time and the time derivative appears in a
mixed derivative term, which naturally leads to a formulation of the
pseudospectrum in terms of a generalised eigenvalue problem
trefethen2005spectra . This choice of coordinates together with a suitable
rescaling of the field not only substantially simplify the problem and the
imposition of boundary conditions, but also significantly reduced the
condition number of the associated matrices (defined as the ratio of the
largest and smallest eigenvalue) making the numerical implementation less
demanding. Comparison of our results with those of Arean:2023ejh  supports the
statement that the pseudospectrum is gauge-independent.

It would be interesting to investigate the implications of spectral
instability in the context of the phenomenon of pole-skipping
Grozdanov:2017ajz . Holographically, pole-skipping has been studied for the
AdS-Schwarzschild black brane Grozdanov:2017ajz , AdS-RN Jansen:2020hfd  as
well as for a neutral black hole that explicitly breaks translation invariance
Blake:2019otz . This phenomenon is now understood to be generic, manifesting
itself in all hydrodynamic and non-hydrodynamic channels Grozdanov:2019uhi ;
Blake:2019otz , and it has been associated with absence of a unique ingoing
solution close to the black brane horizon Blake:2019otz . Pole-skipping is
particularly interesting in the case of the hydrodynamic sound mode, where it
was shown that when the sound mode is driven to instability by a choice of a
specific value of imaginary momentum $k=ik_{*}$, then the retarded two-point
function exhibits an exponential growth related to chaos; the frequency and
momentum are given by the holographic Lyapunov exponent and the butterfly
velocity Grozdanov:2017ajz . If spectral instability persists for imaginary
$k$, we expect the physical relevance of this phenomenon to be substantially
reduced and one would need to understand the implications of this in the
context of quantum chaos. In any case, it would be interesting to study the
contour lines of the pseudospectrum slightly before and slightly after the
pole-skipping point. A technical complication here is that the momentum $k$ is
taken to be imaginary, but one needs to ensure that the energy norm remains
real.

In addition, it would also be interesting to investigate the impact of these
results in the context of the radius of convergence of the hydrodynamic
expansion, within the framework of holography. Strictly speaking, in order to
compute the radius of convergence one needs to perform a perturbative
expansion of the dispersion relations in small momenta to a very high order
Withers:2018srf ; Jansen:2020hfd , extract the coefficients and then perform a
fit. In Jansen:2020hfd  it was found that the radius of convergence for a
charged holographic fluid extracted this way agrees quantitatively with the
locations of the branch points obtained by considering the analyticity
properties of the dispersion relations of the hydrodynamic modes on the
complex frequency and momentum plane. Understanding how the two computations
are modified given the instability of the spectrum and how they compare will
be valuable.

In this work we have also probed the extremal limit of the RN black brane, by
considering $T/\mu\to 0$ Edalati:2010pn ; Edalati:2010hk . We found that the
pseudospectrum was sufficiently deformed along the negative imaginary axis to
accommodate for the additional modes that appear there. It would be
interesting to repeat the pseudospectrum analysis starting directly at zero
temperature.

##  Appendix A Massless scalar perturbations in global AdS4

To contrast with the cases presented in the main text, it is illuminating to
analyse a system where there is energy conservation. One such example is that
of pure AdS with reflecting boundary conditions. Since there is no black hole
horizon, and the (timelike) AdS boundary acts like a box, no energy is lost
through the boundaries of the domain.

Let us then consider the familiar case of global AdS4, with metric

$\displaystyle
ds^{2}=-f(r)dt^{2}+\frac{dr^{2}}{f(r)}+r^{2}(d\theta^{2}+\sin^{2}\theta\,d\varphi^{2}),$
(73)

where $f(r)=r^{2}+1$. As before, the equation to solve is

$\square\Phi=0.$ (74)

In this case, the ingoing Eddington-Finkelstein coordinates used in the main
text are not particularly helpful, so we will proceed with the global $(t,r)$
coordinates instead. With the ansatz

$\Phi=\frac{1}{r}\phi(t,r)Y^{m}_{\ell}(\theta,\varphi)$

and coordinate transformation

$\frac{dr}{dx}=f(r)\equiv r^{2}+1,$ (75)

the Klein-Gordon equation reduces to

$\displaystyle\left(-\partial^{2}_{t}+\partial^{2}_{x}-V(r)\right)\phi(t,r)=0,$
(76) $\displaystyle V(r)\equiv f(r)\left(2+\frac{\ell(\ell+1)}{r^{2}}\right).$

We can explicitly solve equation (75) to find

$r=\tan(x),$ (77)

so that $x=0$ at $r=0$ and $x=\frac{\pi}{2}$ at $r=\infty$.

Due to the lack of energy dissipation, this spacetime features _normal_ modes,
and in fact the spectrum for scalar perturbations is given by Berti:2009kk

$\omega=2n+1+\ell,\qquad n=0,1,2,\ldots$ (78)

To write equation (76) in a first-order (in time) form, we introduce
$\Pi\equiv\partial_{t}\phi$. It is also convenient to have the compact
$x$-domain spanning the interval $x\in[0,1]$, so we further do
$x\to\frac{\pi}{2}x$. Equation (76) then becomes

$\partial_{t}\begin{pmatrix}\Pi\\\
\phi\end{pmatrix}=\begin{pmatrix}0&\frac{4}{\pi^{2}}\partial^{2}_{x}-V\mathbb{1}\\\
\mathbb{1}&0\end{pmatrix}\begin{pmatrix}\Pi\\\ \phi\end{pmatrix}.$ (79)

Defining $u=(\Pi,\phi)^{T}$, we can write this as

$\partial_{t}u=i\hat{L}\,u,\qquad
i\hat{L}\equiv\begin{pmatrix}0&\frac{4}{\pi^{2}}\partial^{2}_{x}-V\mathbb{1}\\\
\mathbb{1}&0\end{pmatrix}$ (80)

where

$V\equiv(r^{2}+1)\left(2+\frac{\ell(\ell+1)}{r^{2}}\right),\quad
r=\tan\left(\frac{\pi}{2}x\right)$ (81)

$V$ becomes singular at $x=0$ and $x=1$, but note that $V\phi$ is regular (and
indeed zero) at these points since

$\phi\underset{r\sim\infty}{\sim}\frac{1}{r^{2}},\qquad\phi\underset{r\sim
0}{\sim}r^{\ell+1}.$ (82)

As in section 2.2, we can define an energy norm given by

$E=\frac{\pi}{4}\int_{0}^{1}\left(\partial_{t}\bar{\phi}\,\partial_{t}\phi+\frac{4}{\pi^{2}}\partial_{x}\bar{\phi}\,\partial_{x}\phi+V\bar{\phi}\,\phi\right)dx$
(83)

which motivates the definition of the scalar product

$\langle
u_{1},u_{2}\rangle=\frac{\pi}{4}\int_{0}^{1}\left(\bar{\Pi}_{1}\,\Pi_{2}+\frac{4}{\pi^{2}}\partial_{x}\bar{\phi}_{1}\,\partial_{x}\phi_{2}+V\bar{\phi}_{1}\,\phi_{2}\right)dx\,.$
(84)

With this definition and our choice of boundary conditions $u(x=0)=0=u(x=1)$,
we can check that $\langle
u_{1},\hat{L}u_{2}\rangle=\langle\hat{L}u_{1},u_{2}\rangle$, and so
$\hat{L}^{\dagger}=\hat{L}$ as expected.

We now follow the approach outlined in section 3, where we define the Gram
matrix $G^{E}$ for the discretised scalar product. In particular, we define

$\displaystyle\langle u_{1},u_{2}\rangle_{{}_{G}}$
$\displaystyle=u_{1}^{*}\cdot G^{E}\cdot u_{2}=(\bar{\Pi}_{1}\
\bar{\phi}_{1})\left(\begin{array}[]{c|c}G^{E}_{1}&0\\\ \hline\cr
0&G^{E}_{2}\end{array}\right)\begin{pmatrix}\Pi_{2}\\\ \phi_{2}\end{pmatrix}$
(87)

where

$\displaystyle G^{E}_{1}=\frac{\pi}{4}C_{1},\qquad
G^{E}_{2}=\frac{\pi}{4}\left(C_{V}+\mathbb{D}^{t}\cdot
C_{2}\cdot\mathbb{D}\right)$
$\displaystyle(C_{1})_{ij}=W_{i}\,\delta_{ij},\qquad(C_{2})_{ij}=\frac{4}{\pi^{2}}W_{i}\,\delta_{ij},\qquad(C_{V})_{ij}=V(x_{i})\,W_{i}\,\delta_{ij}$

(no sum in $i$) and $W_{i}$ are the quadrature weights introduced in section
3.2. Finally, note that the end points where the potential $V$ diverges
contribute to zero to the integration due to the boundary conditions used.
Therefore the corresponding rows and columns can be removed from the $G^{E}$
matrix. We can now proceed as outlined in section 3.3 to compute the
pseudospectrum, using equation (37).

Figure 11: The $\epsilon$-pseudospectrum for the $\ell=2$ case of global AdS4.
The contour lines correspond to level sets of $\epsilon$, as defined in
equation (36), and the red crosses mark the spectrum points from equation
(78). _Left:_ a wide view of the pseudospectrum. _Right:_ a zoomed view of the
pseudospectrum near the mode $\omega_{0}=3$.

In figure 11 we observe an example of this computation. In sharp contrast with
the case of figure 2, the contour levels are now perfectly spherical,
indicating stability of the spectrum.

##  Appendix B Convergence

To measure the convergence of the QNM spectra, the drift ratio $\delta_{r}$
was adopted as defined in BOYD199611 . This measures the weighted distance
between eigenvalues at spectral order $N_{1}$ and $N_{2}$ with $N_{1}<N_{2}$
via

$\displaystyle\delta^{(j)}_{r}=\frac{\min(|\lambda_{j}|,\sigma_{j})}{|\lambda_{j}-\bar{\lambda}_{2}|}$
(88)

where $\bar{\lambda}_{2}$ is the nearest eigenvalue in the $N_{2}$ set to the
current eigenvalue $\lambda_{j}$ of the $N_{1}$ set based on
$\min|\lambda_{j}-\lambda_{2}|$, and the weights $\sigma_{j}$ are given by

$\displaystyle\sigma_{j}=\begin{dcases}|\lambda_{j}-\lambda_{1}|&\text{if}\;j=0\\\
|\lambda_{N-1}-\lambda_{j}|&\text{if}\;j=N\\\
\frac{1}{2}\left(|\lambda_{j-1}-\lambda_{j}|+|\lambda_{j+1}-\lambda_{j}|\right)&\text{if}\;0<j<N\end{dcases}$
(89)

In figure 12, we demonstrate this measure of convergence when applied to the
scalar field spectrum shown in figures 1(a) and 1(d). We plot the logarithm of
the inverse of the drift ratio against the mode number for several choices of
spectral degree, $N$. Each set of points is compared with the spectrum of the
next-largest spectral degree, i.e. when ${N_{1}=30}$, ${N_{2}=40}$, etc.
Convergent eigenvalues have large inverse drift ratio compared to the
background. In the uncharged case, eigenvalues appear as complex conjugate
pairs while in the near-extremal case there is a single non-oscillating mode
before the appearance of conjugate pairs. The benefit of the drift ratio
measure is that it is robust against spurious modes arising from numerical
errors, as these will never remain close to meaningful eigenvalues as the
spectral degree changes.

(a) Convergence of the spectrum for the ${k=0}$, ${Q=0}$ planar black hole in
AdS5.

(b) Convergence of the spectrum for the extremal $k=0$,
$Q=0.99\,Q_{\mathrm{max}}$ planar black hole in AdS5.

Figure 12: Convergence of QNM spectra for different planar black holes in AdS5
in terms of the drift ratio, $\delta_{r}$. In each case, 96 bits of precision
were used. Here $Q_{\mathrm{max}}=\sqrt{2}$.

###### Acknowledgements.

It is a pleasure to thank Rodrigo Panosso Macedo and Benjamin Withers for
discussions. C.P. thanks the Princeton Center for Theoretical Science for
hospitality during the programme “Workshop on Nonlinear Aspects of General
Relativity”. C.P. and B.C. acknowledge support from a Royal Society – Science
Foundation Ireland University Research Fellowship via grant URF/R1/211027.
M.Z. acknowledges financial support by the Center for Research and Development
in Mathematics and Applications (CIDMA) through the Portuguese Foundation for
Science and Technology (FCT – Fundação para a Ciência e a Tecnologia) –
references UIDB/04106/2020 and UIDP/04106/2020 – as well as FCT projects
2022.00721.CEECIND, CERN/FIS-PAR/0027/2019, PTDC/FIS-AST/3041/2020, CERN/FIS-
PAR/0024/2021 and 2022.04560.PTDC. This work has further been supported by the
European Horizon Europe staff exchange (SE) programme HORIZON-MSCA-2021-SE-01
Grant No. NewFunFiCO-101086251.

## References

  * (1) L. Trefethen and M. Embree, _Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators_ , Princeton University Press (2005).
  * (2) R. Trefethen, Trefethen and Driscoll, _Hydrodynamic stability without eigenvalues_ , _Sciene_ 261 (1993) 578.
  * (3) H.-P. Nollert and R.H. Price, _Quantifying excitations of quasinormal mode systems_ , _J. Math. Phys._ 40 (1999) 980 [gr-qc/9810074].
  * (4) H.-P. Nollert, _About the significance of quasinormal modes of black holes_ , _Phys. Rev. D_ 53 (1996) 4397 [gr-qc/9602032].
  * (5) J.L. Jaramillo, R. Panosso Macedo and L. Al Sheikh, _Pseudospectrum and Black Hole Quasinormal Mode Instability_ , _Phys. Rev. X_ 11 (2021) 031003 [2004.06434].
  * (6) J.L. Jaramillo, R. Panosso Macedo and L.A. Sheikh, _Gravitational Wave Signatures of Black Hole Quasinormal Mode Instability_ , _Phys. Rev. Lett._ 128 (2022) 211102 [2105.03451].
  * (7) K. Destounis, R.P. Macedo, E. Berti, V. Cardoso and J.L. Jaramillo, _Pseudospectrum of Reissner-Nordström black holes: Quasinormal mode instability and universality_ , _Phys. Rev. D_ 104 (2021) 084091 [2107.09673].
  * (8) D. Arean, D. Garcia-Fariña and K. Landsteiner, _Pseudospectra of Holographic Quasinormal Modes_ , 2307.08751.
  * (9) S. Sarkar, M. Rahman and S. Chakraborty, _Perturbing the perturbed: Stability of quasinormal modes in presence of a positive cosmological constant_ , _Phys. Rev. D_ 108 (2023) 104002 [2304.06829].
  * (10) H.-P. Nollert, _TOPICAL REVIEW: Quasinormal modes: the characteristic ‘sound’ of black holes and neutron stars_ , _Class. Quant. Grav._ 16 (1999) R159.
  * (11) LISA Consortium Waveform Working Group, _Waveform modelling for the laser interferometer space antenna_ , 2311.01300.
  * (12) M.H.-Y. Cheung, K. Destounis, R.P. Macedo, E. Berti and V. Cardoso, _Destabilizing the Fundamental Mode of Black Holes: The Elephant and the Flea_ , _Phys. Rev. Lett._ 128 (2022) 111103 [2111.05415].
  * (13) E. Berti, V. Cardoso, M.H.-Y. Cheung, F. Di Filippo, F. Duque, P. Martens et al., _Stability of the fundamental quasinormal mode in time-domain observations against small perturbations_ , _Phys. Rev. D_ 106 (2022) 084011 [2205.08547].
  * (14) V. Boyanov, K. Destounis, R. Panosso Macedo, V. Cardoso and J.L. Jaramillo, _Pseudospectrum of horizonless compact objects: A bootstrap instability mechanism_ , _Phys. Rev. D_ 107 (2023) 064012 [2209.12950].
  * (15) G. Policastro, D.T. Son and A.O. Starinets, _From AdS / CFT correspondence to hydrodynamics_ , _JHEP_ 09 (2002) 043 [hep-th/0205052].
  * (16) P.K. Kovtun and A.O. Starinets, _Quasinormal modes and holography_ , _Phys. Rev. D_ 72 (2005) 086009 [hep-th/0506184].
  * (17) C.M. Warnick, _On quasinormal modes of asymptotically anti-de Sitter black holes_ , _Commun. Math. Phys._ 333 (2015) 959 [1306.5760].
  * (18) H. Yang and J. Zhang, _Spectral stability of near-extremal spacetimes_ , _Phys. Rev. D_ 107 (2023) 064045 [2210.01724].
  * (19) S.C. Reddy and L.N. Trefethen, _Lax-stability of fully discrete spectral methods via stability regions and pseudo-eigenvalues_ , _Computer methods in applied mechanics and engineering_ 80 (1990) 147.
  * (20) E. Gasperin and J.L. Jaramillo, _Energy scales and black hole pseudospectra: the structural role of the scalar product_ , _Class. Quant. Grav._ 39 (2022) 115010 [2107.12865].
  * (21) J.L. Van Dorsselaer, _Pseudospectra for matrix pencils and stability of equilibria_ , _BIT Numerical Mathematics_ 37 (1997) 833.
  * (22) K.S. Riedel, _Generalized epsilon-pseudospectra_ , _SIAM journal on numerical analysis_ 31 (1994) 1219.
  * (23) B.G.W. Stewart, _Gershgorin theory for the generalized eigenvalue problem_ , _Mathematics of Computation_ 29 (1975) 600.
  * (24) J. Demmel and B. Kågström, _The generalized schur decomposition of an arbitrary pencil a– $\lambda$b—robust software with error bounds and applications. part i: Theory and algorithms_, _ACM Trans. Math. Softw._ 19 (1993) 160–174.
  * (25) E. Berti, V. Cardoso and A.O. Starinets, _Quasinormal modes of black holes and black branes_ , _Class. Quant. Grav._ 26 (2009) 163001 [0905.2975].
  * (26) R. Bartnik, _The wave equation_ , in _Instructional Workshop on Analysis and Geometry, Part 2_ , vol. 34, pp. 185–195, Australian National University, Mathematical Sciences Institute, 1996.
  * (27) P. Grandclement and J. Novak, _Spectral methods for numerical relativity_ , _Living Rev. Rel._ 12 (2009) 1 [0706.2286].
  * (28) J. Boyd, _Chebyshev and Fourier Spectral Methods: Second Revised Edition_ , Dover Books on Mathematics, Dover Publications (2013).
  * (29) N. Halko, P.-G. Martinsson and J.A. Tropp, _Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions_ , _SIAM review_ 53 (2011) 217.
  * (30) G.T. Horowitz and V.E. Hubeny, _Quasinormal modes of AdS black holes and the approach to thermal equilibrium_ , _Phys. Rev. D_ 62 (2000) 024027 [hep-th/9909056].
  * (31) R. Panosso Macedo, _Hyperboloidal approach for static spherically symmetric spacetimes: a didactical introduction and applications in black-hole physics_ , 2307.15735.
  * (32) R.A. Konoplya and A. Zhidenko, _First few overtones probe the event horizon geometry_ , 2209.00679.
  * (33) H. Kodama and A. Ishibashi, _Master equations for perturbations of generalized static black holes with charge in higher dimensions_ , _Prog. Theor. Phys._ 111 (2004) 29 [hep-th/0308128].
  * (34) A. Jansen, A. Rostworowski and M. Rutkowski, _Master equations and stability of Einstein-Maxwell-scalar black holes_ , _JHEP_ 12 (2019) 036 [1909.04049].
  * (35) S. Aretakis, _Stability and Instability of Extreme Reissner-Nordström Black Hole Spacetimes for Linear Scalar Perturbations II_ , _Annales Henri Poincare_ 12 (2011) 1491 [1110.2009].
  * (36) M. Edalati, J.I. Jottar and R.G. Leigh, _Shear Modes, Criticality and Extremal Black Holes_ , _JHEP_ 04 (2010) 075 [1001.0779].
  * (37) M. Edalati, J.I. Jottar and R.G. Leigh, _Holography and the sound of criticality_ , _JHEP_ 10 (2010) 058 [1005.4075].
  * (38) M. Zworski, _Mathematical study of scattering resonances_ , _Bull. Math. Sci._ 7 (2017) 1–85.
  * (39) S. Grozdanov, K. Schalm and V. Scopelliti, _Black hole scrambling from hydrodynamics_ , _Phys. Rev. Lett._ 120 (2018) 231601 [1710.00921].
  * (40) A. Jansen and C. Pantelidou, _Quasinormal modes in charged fluids at complex momentum_ , _JHEP_ 10 (2020) 121 [2007.14418].
  * (41) M. Blake, R.A. Davison and D. Vegh, _Horizon constraints on holographic Green’s functions_ , _JHEP_ 01 (2020) 077 [1904.12883].
  * (42) S. Grozdanov, P.K. Kovtun, A.O. Starinets and P. Tadić, _The complex life of hydrodynamic modes_ , _JHEP_ 11 (2019) 097 [1904.12862].
  * (43) B. Withers, _Short-lived modes from hydrodynamic dispersion relations_ , _JHEP_ 06 (2018) 059 [1803.08058].
  * (44) J.P. Boyd, _Traps and snares in eigenvalue calculations with application to pseudospectral computations of ocean tides in a basin bounded by meridians_ , _Journal of Computational Physics_ 126 (1996) 11.